Doxygen Source Code Documentation
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
|
|
|
|
|
|
|
|
Function Documentation
|
------------------------------------------------------------------ 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 } |
|
Definition at line 122 of file mri_lsqfit.c. 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 } |
|
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 } |
|
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 } |
|
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 } |
|
macros will be used in routines below, as well! * Definition at line 243 of file mri_lsqfit.c. 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 } |