00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 float * mri_lsqfit( MRI_IMAGE * fitim , MRI_IMARR * refim , MRI_IMAGE * wtim )
00019 {
00020 float *fit = NULL ;
00021 MRI_IMAGE *ffitim , *tim , *wim ;
00022 MRI_IMARR *frefim ;
00023
00024 int ii , jj , npix,nref ;
00025 float **refar , *fitar , *war ;
00026
00027 ENTRY("mri_lsqfit") ;
00028
00029
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
00083
00084 fit = lsqfit( npix , fitar , war , nref , refar ) ;
00085
00086
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
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
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
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
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
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
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
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
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
00211
00212 alpha = (float *) malloc( sizeof(float) * nref ) ;
00213 for( ii=0 ; ii < nref ; ii++ ) alpha[ii] = rr[ii] ;
00214
00215
00216
00217 CLEANUP ; return alpha ;
00218 }
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
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
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
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
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
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
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
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
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
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
00356
00357 free(rr) ;
00358 return alpha ;
00359 }
00360
00361
00362
00363
00364
00365
00366
00367 double * mri_startup_lsqfit( MRI_IMARR * refim , MRI_IMAGE * wtim )
00368 {
00369 double *cc = NULL ;
00370 int ii , npix,nref ;
00371 float * wtar , ** refar ;
00372
00373 ENTRY("mri_startup_lsqfit") ;
00374
00375
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
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 }