00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 /*********************************************************************
00008    These are the template routines to form functional images
00009    recursively from various data types.
00010 
00011    To create actual routines for this purpose, you must compile
00012    the file with the preprocessor symbol DTYPE set to one of
00013    the following types:
00014 
00015       byte short int float
00016 
00017       cc -c -DDTYPE=short afni_pcor.c
00018       mv -f afni_pcor.o afni_pcor_short.o
00019 
00020    In the example above, the resulting routine will be named
00021    PCOR_update_short and will take as input a "short *"
00022    (plus the other stuff, which isn't DTYPE dependent).
00023 **********************************************************************/
00024 
00025 #ifndef DTYPE
00026 #error "Cannot compile, since DTYPE is undefined."
00027 #endif
00028 
00029 #undef MAIN
00030 #include "afni_pcor.h"
00031 
00032 /** macros for function names defined in this file **/
00033 
00034 #define PCOR_UPDATE TWO_TWO(PCOR_update_,DTYPE)
00035 
00036 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00037 
00038 /*** update all voxels with a new array of data
00039       inputs:  vdata = array[nvox] of new data for each voxel
00040                ref   = pointer to references structure to use
00041                vc    = pointer to correlation data structure
00042       output:  updated vc
00043 ***/
00044 
00045 void PCOR_UPDATE( DTYPE * vdata , PCOR_references * ref , PCOR_voxel_corr * vc )
00046 {
00047    int vox , jj ,       /* loop indices */
00048        nv = vc->nvox ,  /* number of voxels */
00049        nr = vc->nref ;  /* number of references */
00050 
00051    float *aaa = ref->alp ,
00052          *fff = ref->ff  ,
00053          *ggg = ref->gg  ;
00054 
00055    float zz , bq = ref->betasq ;
00056 
00057    /*** check inputs for OK-ness ***/
00058 
00059    if( vc->nref != ref->nref ){
00060       fprintf( stderr , "PCOR_UPDATE: reference size mismatch!\n" ) ;
00061       EXIT(1) ;
00062    }
00063 
00064 /**----------------------------------------------------------------------
00065     innermost loop expansion is for speedup if nref is small, if enabled
00066       UPZZ updates zz for each row element  > This pair is performed for
00067       UPCH updates the Cholesky row element > each non-diagonal element
00068       UPLL updates the Cholesky diagonal element
00069 -------------------------------------------------------------------------**/
00070 
00071 #ifdef EXPAND_UPDATE
00072 #  define UPZZ(j)  zz -= aaa[j] * VCH(vc,vox,j)
00073 #  define UPCH(j)  VCH(vc,vox,j) = fff[j] * VCH(vc,vox,j) + ggg[j] * zz
00074 #  define UPLL(j)  VCH(vc,vox,j) += bq * zz * zz
00075 
00076    switch( nr ){
00077    default:       /*** generic case: nested loops ***/
00078 #endif
00079 
00080    /*** for each voxel ***/
00081 
00082    for( vox=0 ; vox < nv ; vox++ ){
00083 
00084       /*** update last row of each Cholesky factor ***/
00085 
00086       zz = (float) vdata[vox] ;
00087       for( jj=0 ; jj < nr ; jj++ ){
00088          zz            -= aaa[jj] * VCH(vc,vox,jj) ;
00089          VCH(vc,vox,jj) = fff[jj] * VCH(vc,vox,jj) + ggg[jj] * zz ;
00090       }
00091       VCH(vc,vox,nr) += bq * zz * zz ; /* square of true Cholesky diagonal */
00092    }
00093 
00094 #ifdef EXPAND_UPDATE
00095    break ;
00096 
00097    /***------------------------------------------------------------
00098      Below here are the special cases for 1-9 reference waveforms:
00099         The above loop over jj is just unrolled manually
00100         (using the UP?? macros) in order to gain speed.
00101    ----------------------------------------------------------------***/
00102 
00103    case 1:
00104       for( vox=0 ; vox < nv ; vox++ ){
00105          zz = (float) vdata[vox] ;
00106          UPZZ(0) ; UPCH(0) ; UPLL(1) ;
00107       }
00108    break ;
00109 
00110    case 2:
00111       for( vox=0 ; vox < nv ; vox++ ){
00112          zz = (float) vdata[vox] ;
00113          UPZZ(0) ; UPCH(0) ;
00114          UPZZ(1) ; UPCH(1) ; UPLL(2) ;
00115       }
00116    break ;
00117 
00118    case 3:
00119       for( vox=0 ; vox < nv ; vox++ ){
00120          zz = (float) vdata[vox] ;
00121          UPZZ(0) ; UPCH(0) ;
00122          UPZZ(1) ; UPCH(1) ;
00123          UPZZ(2) ; UPCH(2) ; UPLL(3) ;
00124    }
00125    break ;
00126 
00127    case 4:
00128       for( vox=0 ; vox < nv ; vox++ ){
00129          zz = (float) vdata[vox] ;
00130          UPZZ(0) ; UPCH(0) ;
00131          UPZZ(1) ; UPCH(1) ;
00132          UPZZ(2) ; UPCH(2) ;
00133          UPZZ(3) ; UPCH(3) ; UPLL(4) ;
00134    }
00135    break ;
00136 
00137    case 5:
00138       for( vox=0 ; vox < nv ; vox++ ){
00139          zz = (float) vdata[vox] ;
00140          UPZZ(0) ; UPCH(0) ;
00141          UPZZ(1) ; UPCH(1) ;
00142          UPZZ(2) ; UPCH(2) ;
00143          UPZZ(3) ; UPCH(3) ;
00144          UPZZ(4) ; UPCH(4) ; UPLL(5) ;
00145    }
00146    break ;
00147 
00148    case 6:
00149       for( vox=0 ; vox < nv ; vox++ ){
00150          zz = (float) vdata[vox] ;
00151          UPZZ(0) ; UPCH(0) ;
00152          UPZZ(1) ; UPCH(1) ;
00153          UPZZ(2) ; UPCH(2) ;
00154          UPZZ(3) ; UPCH(3) ;
00155          UPZZ(4) ; UPCH(4) ;
00156          UPZZ(5) ; UPCH(5) ; UPLL(6) ;
00157    }
00158    break ;
00159 
00160    case 7:
00161       for( vox=0 ; vox < nv ; vox++ ){
00162          zz = (float) vdata[vox] ;
00163          UPZZ(0) ; UPCH(0) ;
00164          UPZZ(1) ; UPCH(1) ;
00165          UPZZ(2) ; UPCH(2) ;
00166          UPZZ(3) ; UPCH(3) ;
00167          UPZZ(4) ; UPCH(4) ;
00168          UPZZ(5) ; UPCH(5) ;
00169          UPZZ(6) ; UPCH(6) ; UPLL(7) ;
00170    }
00171    break ;
00172 
00173    case 8:
00174       for( vox=0 ; vox < nv ; vox++ ){
00175          zz = (float) vdata[vox] ;
00176          UPZZ(0) ; UPCH(0) ;
00177          UPZZ(1) ; UPCH(1) ;
00178          UPZZ(2) ; UPCH(2) ;
00179          UPZZ(3) ; UPCH(3) ;
00180          UPZZ(4) ; UPCH(4) ;
00181          UPZZ(5) ; UPCH(5) ;
00182          UPZZ(6) ; UPCH(6) ;
00183          UPZZ(7) ; UPCH(7) ; UPLL(8) ;
00184    }
00185    break ;
00186 
00187    case 9:
00188       for( vox=0 ; vox < nv ; vox++ ){
00189          zz = (float) vdata[vox] ;
00190          UPZZ(0) ; UPCH(0) ;
00191          UPZZ(1) ; UPCH(1) ;
00192          UPZZ(2) ; UPCH(2) ;
00193          UPZZ(3) ; UPCH(3) ;
00194          UPZZ(4) ; UPCH(4) ;
00195          UPZZ(5) ; UPCH(5) ;
00196          UPZZ(6) ; UPCH(6) ;
00197          UPZZ(7) ; UPCH(7) ;
00198          UPZZ(8) ; UPCH(8) ; UPLL(9) ;
00199    }
00200    break ;
00201 
00202    }
00203 #endif /* EXPAND_UPDATE */
00204 
00205    (vc->nupdate)++ ;  /* June 1995: another update completed! */
00206    return ;
00207 }