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 #include "mrilib.h"
00008 
00009 /*** 7D SAFE ***/
00010 
00011 /*-----------------------------------------------------------------
00012    Routine to compute the least squares fit of one image
00013    (fitim) to a linear combination of other images (refim).
00014    If not NULL, wtim is an image of weighting factors.
00015    Returns a malloc-ed array of floats (the fit coefficients).
00016 -------------------------------------------------------------------*/
00017 
00018 float * mri_lsqfit( MRI_IMAGE * fitim , MRI_IMARR * refim , MRI_IMAGE * wtim )
00019 {
00020    float *fit = NULL ;                 /* will be the output */
00021    MRI_IMAGE *ffitim , *tim , *wim ;   /* local versions of inputs */
00022    MRI_IMARR *frefim ;
00023 
00024    int ii , jj , npix,nref ;
00025    float **refar , *fitar , *war ;
00026 
00027 ENTRY("mri_lsqfit") ;
00028 
00029    /****---- check inputs, convert to float type if needed ----****/
00030 
00031    if( fitim == NULL ){
00032      fprintf(stderr,"mri_lsqfit: NULL fitim!\a\n"); RETURN(NULL);
00033    }
00034 
00035    if( fitim->kind == MRI_float ) ffitim = fitim ;
00036    else                           ffitim = mri_to_float( fitim ) ;
00037    npix  = ffitim->nvox ;
00038    fitar = mri_data_pointer(ffitim) ;
00039 
00040    if( wtim == NULL ){
00041       wim = NULL ;
00042       war = NULL ;
00043    } else if( wtim->kind == MRI_float ){
00044       wim = wtim ;
00045       war = mri_data_pointer( wim ) ;
00046       if( wim->nvox != npix ){
00047          fprintf(stderr,"mri_lsqfit: MISMATCH wtim\a\n"); RETURN(NULL);
00048       }
00049    } else {
00050       wim = mri_to_float( wtim ) ;
00051       war = mri_data_pointer( wim ) ;
00052       if( wim->nvox != npix ){
00053          fprintf(stderr,"mri_lsqfit: MISMATCH wtim\a\n"); RETURN(NULL);
00054       }
00055    }
00056 
00057    if( refim == NULL || refim->num < 1 ){
00058       fprintf(stderr,"mri_lsqfit: NULL refim!\a\n"); RETURN(NULL);
00059    }
00060 
00061    nref = refim->num ;
00062 
00063    INIT_IMARR(frefim) ;
00064    refar = (float **) malloc( sizeof(float *) * nref ) ;
00065    if( refar == NULL ){
00066       fprintf(stderr,"mri_lsqfit: malloc failure for refar!\a\n"); RETURN(NULL);
00067    }
00068 
00069    for( ii=0 ; ii < nref ; ii++ ){
00070       if( refim->imarr[ii] == NULL ){
00071          fprintf(stderr,"mri_lsqfit: NULL refim[%d]!\a\n",ii); RETURN(NULL);
00072       }
00073       if( refim->imarr[ii]->nvox != npix ){
00074          fprintf(stderr,"mri_lsqfit: MISMATCH refim[%d]!\a\n",ii); RETURN(NULL);
00075       }
00076       if( refim->imarr[ii]->kind == MRI_float ) tim = refim->imarr[ii] ;
00077       else                                      tim = mri_to_float(refim->imarr[ii]) ;
00078       ADDTO_IMARR(frefim,tim) ;
00079       refar[ii] = mri_data_pointer(tim) ;
00080    }
00081 
00082    /****---- get coefficients ----****/
00083 
00084    fit = lsqfit( npix , fitar , war , nref , refar ) ;
00085 
00086    /****---- clean up and exit ----****/
00087 
00088    if( ffitim != fitim ) mri_free( ffitim ) ;
00089    if( wim != NULL && wim != wtim ) mri_free( wim ) ;
00090    for( ii=0 ; ii < nref ; ii++ ){
00091       if( frefim->imarr[ii] != refim->imarr[ii] ) mri_free( frefim->imarr[ii] ) ;
00092    }
00093    FREE_IMARR(frefim) ;
00094    free(refar) ;
00095 
00096    RETURN(fit) ;
00097 }
00098 
00099 /*----------------------------------------------------------------
00100      Routine to compute the weighted least square fit of
00101      a data vector to a bunch of reference vectors:
00102 
00103  veclen = length of vectors
00104  data   = array holding data vector (length veclen)
00105  wt     = array holding weight for each data point (length veclen)
00106           [if NULL, then a vector of all 1's is used]
00107  nref   = number of reference vectors, each of length veclen;
00108           must have 1 <= nref <= veclen
00109  ref    = array of pointers to reference vectors, so that
00110           ref[i][k] is the k-th component of the i-th reference,
00111           for i=0..nref-1, k=0..veclen-1
00112 
00113  Return value is a pointer to a vector of length nref with the
00114  weights of each reference;  if NULL is returned, there was an
00115  error.  The array is allocated by malloc and so should be freed
00116  when the caller is finished with it.
00117 
00118  Input and output vectors are floats.  Internal calculuations
00119  are done with doubles.
00120 ------------------------------------------------------------------*/
00121 
00122 float * lsqfit( int veclen ,
00123                 float *data , float *wt , int nref , float *ref[] )
00124 {
00125    int    ii , jj , kk ;
00126    float  *alpha = NULL ;
00127    double *cc = NULL , *rr = NULL ;
00128    double sum ;
00129 
00130    /** macros will be used in routines below, as well! **/
00131 
00132 #define DBLEVEC(ll) (double *) malloc( sizeof(double) * (ll) )
00133 #define DISCARD(xx) if( xx != NULL ){free(xx); xx = NULL;}
00134 #define CLEANUP     {DISCARD(cc); DISCARD(rr);}
00135 
00136    if( nref < 1 || veclen < nref || data == NULL || ref == NULL ) return NULL ;
00137 
00138    /*** form RHS vector into rr ***/
00139 
00140    rr = DBLEVEC(nref) ;
00141    cc = DBLEVEC(nref*nref) ;
00142    if( rr == NULL || cc == NULL ){ CLEANUP ; return NULL ; }
00143 
00144    if( wt != NULL ){
00145       for( ii=0 ; ii < nref ; ii++ ){
00146          sum = 0.0 ;
00147          for( jj=0 ; jj < veclen ; jj++ ) sum += ref[ii][jj] * wt[jj] * data[jj] ;
00148          rr[ii] = sum ;
00149       }
00150    } else {
00151       for( ii=0 ; ii < nref ; ii++ ){
00152          sum = 0.0 ;
00153          for( jj=0 ; jj < veclen ; jj++ ) sum += ref[ii][jj] * data[jj] ;
00154          rr[ii] = sum ;
00155       }
00156    }
00157 
00158    /*** form coefficient matrix into cc */
00159 
00160 #define CC(i,j) cc[(i)+(j)*nref]
00161 
00162    if( wt != NULL ){
00163       for( jj=0 ; jj < nref ; jj++ ){
00164          for( ii=0 ; ii <= jj ; ii++ ){
00165             sum = 0.0 ;
00166             for( kk=0 ; kk < veclen ; kk++ ) sum += ref[ii][kk] * ref[jj][kk] * wt[kk] ;
00167             CC(ii,jj) = CC(jj,ii) = sum ;
00168          }
00169       }
00170    } else {
00171       for( jj=0 ; jj < nref ; jj++ ){
00172          for( ii=0 ; ii <= jj ; ii++ ){
00173             sum = 0.0 ;
00174             for( kk=0 ; kk < veclen ; kk++ ) sum += ref[ii][kk] * ref[jj][kk] ;
00175             CC(ii,jj) = CC(jj,ii) = sum ;
00176          }
00177       }
00178    }
00179 
00180    /*** Choleski decompose cc ***/
00181 
00182    for( ii=0 ; ii < nref ; ii++ ){
00183       for( jj=0 ; jj < ii ; jj++ ){
00184          sum = CC(ii,jj) ;
00185          for( kk=0 ; kk < jj ; kk++ ) sum -= CC(ii,kk) * CC(jj,kk) ;
00186          CC(ii,jj) = sum / CC(jj,jj) ;
00187       }
00188       sum = CC(ii,ii) ;
00189       for( kk=0 ; kk < ii ; kk++ ) sum -= CC(ii,kk) * CC(ii,kk) ;
00190       if( sum <= 0.0 ){ CLEANUP ; return NULL ; }
00191       CC(ii,ii) = sqrt(sum) ;
00192    }
00193 
00194    /*** forward solve ***/
00195 
00196    for( ii=0 ; ii < nref ; ii++ ){
00197       sum = rr[ii] ;
00198       for( jj=0 ; jj < ii ; jj++ ) sum -= CC(ii,jj) * rr[jj] ;
00199       rr[ii] = sum / CC(ii,ii) ;
00200    }
00201 
00202    /*** backward solve ***/
00203 
00204    for( ii=nref-1 ; ii >= 0 ; ii-- ){
00205       sum = rr[ii] ;
00206       for( jj=ii+1 ; jj < nref ; jj++ ) sum -= CC(jj,ii) * rr[jj] ;
00207       rr[ii] = sum / CC(ii,ii) ;
00208    }
00209 
00210    /*** put result into alpha ***/
00211 
00212    alpha = (float *) malloc( sizeof(float) * nref ) ;
00213    for( ii=0 ; ii < nref ; ii++ ) alpha[ii] = rr[ii] ;
00214 
00215    /*** cleanup and exit ***/
00216 
00217    CLEANUP ; return alpha ;
00218 }
00219 
00220 /*----------------------------------------------------------------
00221    Similar to above, but only produce the Choleski decomposition
00222    and weight scaled ref vectors for later use in delayed_lsqfit.
00223    This is to be used when fitting several data vectors to the
00224    same references with the same weight factors.
00225 
00226  veclen = length of vectors
00227  wt     = array holding weight for each data point (length veclen)
00228           [if NULL, then a vector of all 1's is used]
00229  nref   = number of reference vectors, each of length veclen;
00230           must have 1 <= nref <= veclen
00231  ref    = array of pointers to reference vectors, so that
00232           ref[i][k] is the k-th component of the i-th reference,
00233           for i=0..nref-1, k=0..veclen-1
00234 
00235    Return value is a (double *) which points to malloc-ed memory
00236    for the Choleksi decomposition.  It should later be passed
00237    to delayed_lsqfit.  If wt != NULL, then ref[ii][jj] is
00238    scaled by wt[jj] as well.
00239 
00240    If NULL is returned, an error occured.
00241 ------------------------------------------------------------------*/
00242 
00243 double * startup_lsqfit( int veclen ,
00244                          float *wt , int nref , float *ref[] )
00245 {
00246    int    ii , jj , kk ;
00247    double *cc = NULL ;
00248    double sum ;
00249 
00250    if( nref < 1 || veclen < nref || ref == NULL ){
00251       fprintf(stderr,"*** Illegal inputs to startup_lsqfit\n") ;
00252       return NULL ;
00253    }
00254 
00255    /*** form coefficient matrix into cc */
00256 
00257    cc = DBLEVEC(nref*nref) ;
00258    if( cc == NULL ){
00259       fprintf(stderr,"Can't malloc workspace in startup_lsqfit\n") ;
00260       return NULL ;
00261    }
00262 
00263    if( wt != NULL ){
00264       for( jj=0 ; jj < nref ; jj++ ){
00265          for( ii=0 ; ii <= jj ; ii++ ){
00266             sum = 0.0 ;
00267             for( kk=0 ; kk < veclen ; kk++ ) sum += ref[ii][kk] * ref[jj][kk] * wt[kk] ;
00268             CC(ii,jj) = CC(jj,ii) = sum ;
00269          }
00270       }
00271    } else {
00272       for( jj=0 ; jj < nref ; jj++ ){
00273          for( ii=0 ; ii <= jj ; ii++ ){
00274             sum = 0.0 ;
00275             for( kk=0 ; kk < veclen ; kk++ ) sum += ref[ii][kk] * ref[jj][kk] ;
00276             CC(ii,jj) = CC(jj,ii) = sum ;
00277          }
00278       }
00279    }
00280 
00281    /*** Choleski decompose cc ***/
00282 
00283    for( ii=0 ; ii < nref ; ii++ ){
00284       for( jj=0 ; jj < ii ; jj++ ){
00285          sum = CC(ii,jj) ;
00286          for( kk=0 ; kk < jj ; kk++ ) sum -= CC(ii,kk) * CC(jj,kk) ;
00287          CC(ii,jj) = sum / CC(jj,jj) ;
00288       }
00289       sum = CC(ii,ii) ;
00290       for( kk=0 ; kk < ii ; kk++ ) sum -= CC(ii,kk) * CC(ii,kk) ;
00291       if( sum <= 0.0 ){
00292          free(cc) ;
00293          fprintf(stderr,"Choleski factorization failure in startup_lsqfit\n") ;
00294          return NULL ;
00295       }
00296       CC(ii,ii) = sqrt(sum) ;
00297    }
00298 
00299    /*** scale ref by wt, if desired ***/
00300 
00301    if( wt != NULL ){
00302       for( ii=0 ; ii < nref ; ii++ )
00303          for( jj=0 ; jj < veclen ; jj++ ) ref[ii][jj] *= wt[jj] ;
00304    }
00305 
00306    return cc ;
00307 }
00308 
00309 /**------------------------------------------------------------------
00310   Given the data from startup_lsqfit, finish the job.
00311 ------------------------------------------------------------------**/
00312 
00313 float * delayed_lsqfit( int veclen ,
00314                         float * data , int nref , float *ref[] , double * cc )
00315 {
00316    int    ii , jj ;
00317    float  *alpha = NULL ;
00318    double *rr = NULL ;
00319    register double sum ;
00320 
00321    if( nref < 1 || veclen < nref ||
00322        data == NULL || ref == NULL || cc == NULL ) return NULL ;
00323 
00324    /*** form RHS vector into rr ***/
00325 
00326    rr = DBLEVEC(nref) ; if( rr == NULL ) return NULL ;
00327 
00328    for( ii=0 ; ii < nref ; ii++ ){
00329       sum = 0.0 ;
00330       for( jj=0 ; jj < veclen ; jj++ ) sum += ref[ii][jj] * data[jj] ;
00331       rr[ii] = sum ;
00332    }
00333 
00334    /*** forward solve ***/
00335 
00336    for( ii=0 ; ii < nref ; ii++ ){
00337       sum = rr[ii] ;
00338       for( jj=0 ; jj < ii ; jj++ ) sum -= CC(ii,jj) * rr[jj] ;
00339       rr[ii] = sum / CC(ii,ii) ;
00340    }
00341 
00342    /*** backward solve ***/
00343 
00344    for( ii=nref-1 ; ii >= 0 ; ii-- ){
00345       sum = rr[ii] ;
00346       for( jj=ii+1 ; jj < nref ; jj++ ) sum -= CC(jj,ii) * rr[jj] ;
00347       rr[ii] = sum / CC(ii,ii) ;
00348    }
00349 
00350    /*** put result into alpha ***/
00351 
00352    alpha = (float *) malloc( sizeof(float) * nref ) ; if( alpha == NULL ) return NULL ;
00353    for( ii=0 ; ii < nref ; ii++ ) alpha[ii] = rr[ii] ;
00354 
00355    /*** cleanup and exit ***/
00356 
00357    free(rr) ;
00358    return alpha ;
00359 }
00360 
00361 /*-----------------------------------------------------------------
00362    'startup' and 'delayed' versions for image fitting.
00363    N.B.: unlike mri_lsqfit, all the images in refim and wtim must
00364          be of MRI_float kind.
00365 -------------------------------------------------------------------*/
00366 
00367 double * mri_startup_lsqfit( MRI_IMARR * refim , MRI_IMAGE * wtim )
00368 {
00369    double *cc = NULL ;                 /* will be the output */
00370    int ii , npix,nref ;
00371    float * wtar , ** refar ;
00372 
00373 ENTRY("mri_startup_lsqfit") ;
00374 
00375    /****---- check inputs ----****/
00376 
00377    if( wtim != NULL && wtim->kind != MRI_float ){
00378       fprintf(stderr,"mri_startup_lsqfit: non-float wtim!\a\n") ; RETURN(NULL);
00379    }
00380    wtar = (wtim == NULL) ? (NULL) : (MRI_FLOAT_PTR(wtim)) ;
00381 
00382    if( refim == NULL || refim->num < 1 ){
00383       fprintf(stderr,"mri_startup_lsqfit: NULL refim!\a\n") ; RETURN(NULL);
00384    }
00385 
00386    nref  = refim->num ;
00387    npix  = refim->imarr[0]->nvox ;
00388    refar = (float **) malloc( sizeof(float *) * nref ) ;
00389    if( refar == NULL ){
00390       fprintf(stderr,"mri_startup_lsqfit: malloc failure for refar!\a\n") ; RETURN(NULL);
00391    }
00392 
00393    for( ii=0 ; ii < nref ; ii++ ){
00394       if( refim->imarr[ii] == NULL ){
00395          fprintf(stderr,"mri_startup_lsqfit: NULL refim[%d]!\a\n",ii) ; RETURN(NULL);
00396       }
00397       if( refim->imarr[ii]->nvox != npix ){
00398          fprintf(stderr,"mri_startup_lsqfit: MISMATCH refim[%d]!\a\n",ii) ; RETURN(NULL);
00399       }
00400       if( refim->imarr[ii]->kind != MRI_float ){
00401          fprintf(stderr,"mri_startup_lsqfit: non-float refim[%d]!\a\n",ii) ; RETURN(NULL);
00402       }
00403       refar[ii] = MRI_FLOAT_PTR(refim->imarr[ii]) ;
00404    }
00405 
00406    /****---- get Choleski, send it out ----****/
00407 
00408    cc = startup_lsqfit( npix , wtar , nref , refar ) ;
00409    if( cc == NULL ){
00410          fprintf(stderr,"mri_startup_lsqfit: bad call to startup_lsqfit!\a\n") ; RETURN(NULL);
00411    }
00412    free(refar) ;
00413    RETURN(cc) ;
00414 }
00415 
00416 /*-----------------------------------------------------------------*/
00417 
00418 float * mri_delayed_lsqfit( MRI_IMAGE * fitim , MRI_IMARR * refim , double * cc )
00419 {
00420    int ii , npix,nref ;
00421    float * fit ;
00422    static float ** refar = NULL ;
00423    static int     nrefar = -1 ;
00424 
00425 ENTRY("mri_delayed_lsqfit") ;
00426 
00427    nref = refim->num ;
00428    npix = refim->imarr[0]->nvox ;
00429 
00430    if( nrefar < nref ){
00431       if( refar != NULL ) free(refar) ;
00432       refar  = (float **) malloc( sizeof(float *) * nref ) ;
00433       nrefar = nref ;
00434    }
00435    if( refar == NULL ){
00436      fprintf(stderr,"mri_delayed_lsqfit: malloc failure for refar!\a\n"); RETURN(NULL);
00437    }
00438 
00439    for( ii=0 ; ii < nref ; ii++ )
00440       refar[ii] = MRI_FLOAT_PTR(refim->imarr[ii]) ;
00441 
00442    fit = delayed_lsqfit( npix , MRI_FLOAT_PTR(fitim) , nref , refar , cc ) ;
00443    RETURN(fit) ;
00444 }