Skip to content

AFNI/NIfTI Server

Sections
Personal tools
You are here: Home » AFNI » Documentation
 

Doxygen Source Code Documentation

Main Page   Alphabetical List   Data Structures   File List   Data Fields   Globals   Search  

mri_lsqfit.c File Reference

#include "mrilib.h"

Go to the source code of this file.


Defines

#define DBLEVEC(ll)   (double *) malloc( sizeof(double) * (ll) )
#define DISCARD(xx)   if( xx != NULL ){free(xx); xx = NULL;}
#define CLEANUP   {DISCARD(cc); DISCARD(rr);}
#define CC(i, j)   cc[(i)+(j)*nref]

Functions

float * mri_lsqfit (MRI_IMAGE *fitim, MRI_IMARR *refim, MRI_IMAGE *wtim)
float * lsqfit (int veclen, float *data, float *wt, int nref, float *ref[])
double * startup_lsqfit (int veclen, float *wt, int nref, float *ref[])
float * delayed_lsqfit (int veclen, float *data, int nref, float *ref[], double *cc)
double * mri_startup_lsqfit (MRI_IMARR *refim, MRI_IMAGE *wtim)
float * mri_delayed_lsqfit (MRI_IMAGE *fitim, MRI_IMARR *refim, double *cc)

Define Documentation

#define CC (  i,
j    )     cc[(i)+(j)*nref]
 

#define CLEANUP   {DISCARD(cc); DISCARD(rr);}
 

#define DBLEVEC (  ll    )     (double *) malloc( sizeof(double) * (ll) )
 

#define DISCARD (  xx    )     if( xx != NULL ){free(xx); xx = NULL;}
 

Function Documentation

float* delayed_lsqfit (  int    veclen,
float *    data,
int    nref,
float *    ref[],
double *    cc
) 
 

------------------------------------------------------------------ Given the data from startup_lsqfit, finish the job. ------------------------------------------------------------------*

Definition at line 313 of file mri_lsqfit.c.

References free, malloc, and ref.

Referenced by LSQ_worker(), main(), and mri_delayed_lsqfit(). 

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 }

float* lsqfit (  int    veclen,
float *    data,
float *    wt,
int    nref,
float *    ref[]
) 
 

Definition at line 122 of file mri_lsqfit.c.

References malloc, and ref.

Referenced by huber_func(), mri_lsqfit(), and THD_generic_detrend(). 

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 }

float* mri_delayed_lsqfit (  MRI_IMAGE *    fitim,
MRI_IMARR *    refim,
double *    cc
) 
 

Definition at line 418 of file mri_lsqfit.c.

References delayed_lsqfit(), ENTRY, fit, free, MRI_IMARR::imarr, malloc, MRI_FLOAT_PTR, MRI_IMARR::num, MRI_IMAGE::nvox, and RETURN.

Referenced by mri_2dalign_one(), mri_3dalign_one(), and mri_align_dfspace(). 

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 }

float* mri_lsqfit (  MRI_IMAGE *    fitim,
MRI_IMARR *    refim,
MRI_IMAGE *    wtim
) 
 

Definition at line 18 of file mri_lsqfit.c.

References ADDTO_IMARR, ENTRY, fit, free, FREE_IMARR, MRI_IMARR::imarr, INIT_IMARR, MRI_IMAGE::kind, lsqfit(), malloc, mri_data_pointer(), mri_free(), mri_to_float(), MRI_IMARR::num, MRI_IMAGE::nvox, and RETURN.

Referenced by main(), and mri_align_dfspace(). 

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 }

double* mri_startup_lsqfit (  MRI_IMARR *    refim,
MRI_IMAGE *    wtim
) 
 

Definition at line 367 of file mri_lsqfit.c.

References ENTRY, free, MRI_IMARR::imarr, MRI_IMAGE::kind, malloc, MRI_FLOAT_PTR, MRI_IMARR::num, MRI_IMAGE::nvox, RETURN, and startup_lsqfit().

Referenced by mri_2dalign_setup(), mri_3dalign_setup(), and mri_align_dfspace(). 

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 }

double* startup_lsqfit (  int    veclen,
float *    wt,
int    nref,
float *    ref[]
) 
 

macros will be used in routines below, as well! *

Definition at line 243 of file mri_lsqfit.c.

References free, and ref.

Referenced by LSQ_worker(), main(), and mri_startup_lsqfit(). 

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 }
 
Powered by Plone

This site conforms to the following standards: