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 }