1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 /*** 7D SAFE ***/
10 /*-----------------------------------------------------------------*/
11 /* Check ref vectors for all zero components [13 Mar 2017]
12 Return value
13 * 0 if all are OK.
14 * -1 if inputs are invalid
15 * > 0 is count of all zero ref vectors
16 *//*---------------------------------------------------------------*/
17
check_ref_vectors(int veclen,int nref,float * ref[],char * funcname)18 int check_ref_vectors( int veclen , int nref , float *ref[] , char *funcname )
19 {
20 int ii,jj , nbad=0 ; float *vv ;
21
22 if( veclen < 1 || nref < 1 || ref == NULL ){
23 if( funcname != NULL )
24 ERROR_message("%s :: bad inputs for data fitting",funcname) ;
25 return -1 ;
26 }
27
28 for( jj=0 ; jj < nref ; jj++ ){
29 vv = ref[jj] ;
30 for( ii=0 ; ii < veclen && vv[ii] == 0.0f ; ii++ ) ; /*nada*/
31 if( ii == veclen ) nbad++ ;
32 }
33
34 if( nbad > 0 && funcname != NULL )
35 ERROR_message("%s :: %d vector%s are all 0 for data fitting",
36 funcname , nbad , (nbad > 1) ? "s" : "\0" ) ;
37
38 return nbad ;
39 }
40
41 /*-----------------------------------------------------------------
42 Routine to compute the least squares fit of one image
43 (fitim) to a linear combination of other images (refim).
44 If not NULL, wtim is an image of weighting factors.
45 Returns a malloc-ed array of floats (the fit coefficients).
46 -------------------------------------------------------------------*/
47
mri_lsqfit(MRI_IMAGE * fitim,MRI_IMARR * refim,MRI_IMAGE * wtim)48 float * mri_lsqfit( MRI_IMAGE * fitim , MRI_IMARR * refim , MRI_IMAGE * wtim )
49 {
50 float *fit = NULL ; /* will be the output */
51 MRI_IMAGE *ffitim , *tim , *wim ; /* local versions of inputs */
52 MRI_IMARR *frefim ;
53
54 int ii , jj , npix,nref ;
55 float **refar , *fitar , *war ;
56
57 ENTRY("mri_lsqfit") ;
58
59 /****---- check inputs, convert to float type if needed ----****/
60
61 if( fitim == NULL ){
62 ERROR_message("mri_lsqfit has NULL fitim"); RETURN(NULL);
63 }
64
65 if( fitim->kind == MRI_float ) ffitim = fitim ;
66 else ffitim = mri_to_float( fitim ) ;
67 npix = ffitim->nvox ;
68 fitar = mri_data_pointer(ffitim) ;
69
70 if( wtim == NULL ){
71 wim = NULL ;
72 war = NULL ;
73 } else if( wtim->kind == MRI_float ){
74 wim = wtim ;
75 war = mri_data_pointer( wim ) ;
76 if( wim->nvox != npix ){
77 ERROR_message("mri_lsqfit: MISMATCH wtim->nvox=%d npix=%d",
78 wtim->nvox , npix );
79 RETURN(NULL);
80 }
81 } else {
82 wim = mri_to_float( wtim ) ;
83 war = mri_data_pointer( wim ) ;
84 if( wim->nvox != npix ){
85 ERROR_message("mri_lsqfit: MISMATCH wtim->nvox=%d npix=%d",
86 wtim->nvox , npix );
87 RETURN(NULL);
88 }
89 }
90
91 if( refim == NULL || refim->num < 1 ){
92 ERROR_message("mri_lsqfit: NULL refim"); RETURN(NULL);
93 }
94
95 nref = refim->num ;
96
97 INIT_IMARR(frefim) ;
98 refar = (float **) malloc( sizeof(float *) * nref ) ;
99 if( refar == NULL ){
100 fprintf(stderr,"** mri_lsqfit: malloc failure for refar"); RETURN(NULL);
101 }
102
103 for( ii=0 ; ii < nref ; ii++ ){
104 if( refim->imarr[ii] == NULL ){
105 ERROR_message("mri_lsqfit: NULL refim[%d]!",ii); RETURN(NULL);
106 }
107 if( refim->imarr[ii]->nvox != npix ){
108 ERROR_message("mri_lsqfit: MISMATCH refim[%d]!",ii); RETURN(NULL);
109 }
110 if( refim->imarr[ii]->kind == MRI_float ) tim = refim->imarr[ii] ;
111 else tim = mri_to_float(refim->imarr[ii]) ;
112 ADDTO_IMARR(frefim,tim) ;
113 refar[ii] = mri_data_pointer(tim) ;
114 }
115
116 /****---- get coefficients ----****/
117
118 fit = lsqfit( npix , fitar , war , nref , refar ) ;
119
120 /****---- clean up and exit ----****/
121
122 if( ffitim != fitim ) mri_free( ffitim ) ;
123 if( wim != NULL && wim != wtim ) mri_free( wim ) ;
124 for( ii=0 ; ii < nref ; ii++ ){
125 if( frefim->imarr[ii] != refim->imarr[ii] ) mri_free( frefim->imarr[ii] ) ;
126 }
127 FREE_IMARR(frefim) ;
128 free(refar) ;
129
130 RETURN(fit) ;
131 }
132
133 /*----------------------------------------------------------------
134 Routine to compute the weighted least square fit of
135 a data vector to a bunch of reference vectors:
136
137 veclen = length of vectors
138 data = array holding data vector (length veclen)
139 wt = array holding weight for each data point (length veclen)
140 [if NULL, then a vector of all 1's is used]
141 nref = number of reference vectors, each of length veclen;
142 must have 1 <= nref <= veclen
143 ref = array of pointers to reference vectors, so that
144 ref[i][k] is the k-th component of the i-th reference,
145 for i=0..nref-1, k=0..veclen-1
146
147 Return value is a pointer to a vector of length nref with the
148 weights of each reference; if NULL is returned, there was an
149 error. The array is allocated by malloc and so should be freed
150 when the caller is finished with it.
151
152 Input and output vectors are floats. Internal calculuations
153 are done with doubles.
154 ------------------------------------------------------------------*/
155
lsqfit(int veclen,float * data,float * wt,int nref,float * ref[])156 float * lsqfit( int veclen ,
157 float *data , float *wt , int nref , float *ref[] )
158 {
159 int ii , jj , kk ;
160 float *alpha = NULL ;
161 double *cc = NULL , *rr = NULL ;
162 double sum ;
163
164 /** macros will be used in routines below, as well! **/
165
166 #define DBLEVEC(ll) (double *) malloc( sizeof(double) * (ll) )
167 #define DISCARD(xx) if( xx != NULL ){free(xx); xx = NULL;}
168 #define CLEANUP {DISCARD(cc); DISCARD(rr);}
169
170 if( nref < 1 || veclen < nref || data == NULL || ref == NULL ) return NULL ;
171
172 /* 13 Mar 2017 */
173 { int nbad = check_ref_vectors( veclen , nref , ref , "lsqfit" ) ;
174 if( nbad ) return NULL ;
175 }
176
177 /*** form RHS vector into rr ***/
178
179 rr = DBLEVEC(nref) ;
180 cc = DBLEVEC(nref*nref) ;
181 if( rr == NULL || cc == NULL ){ CLEANUP ; return NULL ; }
182
183 if( wt != NULL ){
184 for( ii=0 ; ii < nref ; ii++ ){
185 sum = 0.0 ;
186 for( jj=0 ; jj < veclen ; jj++ ) sum += ref[ii][jj] * wt[jj] * data[jj] ;
187 rr[ii] = sum ;
188 }
189 } else {
190 for( ii=0 ; ii < nref ; ii++ ){
191 sum = 0.0 ;
192 for( jj=0 ; jj < veclen ; jj++ ) sum += ref[ii][jj] * data[jj] ;
193 rr[ii] = sum ;
194 }
195 }
196
197 /*** form coefficient matrix into cc */
198
199 #define CC(i,j) cc[(i)+(j)*nref]
200
201 if( wt != NULL ){
202 for( jj=0 ; jj < nref ; jj++ ){
203 for( ii=0 ; ii <= jj ; ii++ ){
204 sum = 0.0 ;
205 for( kk=0 ; kk < veclen ; kk++ ) sum += ref[ii][kk] * ref[jj][kk] * wt[kk] ;
206 CC(ii,jj) = CC(jj,ii) = sum ;
207 }
208 }
209 } else {
210 for( jj=0 ; jj < nref ; jj++ ){
211 for( ii=0 ; ii <= jj ; ii++ ){
212 sum = 0.0 ;
213 for( kk=0 ; kk < veclen ; kk++ ) sum += ref[ii][kk] * ref[jj][kk] ;
214 CC(ii,jj) = CC(jj,ii) = sum ;
215 }
216 }
217 }
218
219 /*** Choleski decompose cc ***/
220
221 for( ii=0 ; ii < nref ; ii++ ){
222 for( jj=0 ; jj < ii ; jj++ ){
223 sum = CC(ii,jj) ;
224 for( kk=0 ; kk < jj ; kk++ ) sum -= CC(ii,kk) * CC(jj,kk) ;
225 CC(ii,jj) = sum / CC(jj,jj) ;
226 }
227 sum = CC(ii,ii) ;
228 for( kk=0 ; kk < ii ; kk++ ) sum -= CC(ii,kk) * CC(ii,kk) ;
229 if( sum <= 0.0 ){ CLEANUP ; return NULL ; }
230 CC(ii,ii) = sqrt(sum) ;
231 }
232
233 /*** forward solve ***/
234
235 for( ii=0 ; ii < nref ; ii++ ){
236 sum = rr[ii] ;
237 for( jj=0 ; jj < ii ; jj++ ) sum -= CC(ii,jj) * rr[jj] ;
238 rr[ii] = sum / CC(ii,ii) ;
239 }
240
241 /*** backward solve ***/
242
243 for( ii=nref-1 ; ii >= 0 ; ii-- ){
244 sum = rr[ii] ;
245 for( jj=ii+1 ; jj < nref ; jj++ ) sum -= CC(jj,ii) * rr[jj] ;
246 rr[ii] = sum / CC(ii,ii) ;
247 }
248
249 /*** put result into alpha ***/
250
251 alpha = (float *) malloc( sizeof(float) * nref ) ;
252 for( ii=0 ; ii < nref ; ii++ ) alpha[ii] = rr[ii] ;
253
254 /*** cleanup and exit ***/
255
256 CLEANUP ; return alpha ;
257 }
258
259 /*----------------------------------------------------------------
260 Similar to above, but only produce the Choleski decomposition
261 and weight scaled ref vectors for later use in delayed_lsqfit.
262 This is to be used when fitting several data vectors to the
263 same references with the same weight factors.
264
265 veclen = length of vectors
266 wt = array holding weight for each data point (length veclen)
267 [if NULL, then a vector of all 1's is used]
268 nref = number of reference vectors, each of length veclen;
269 must have 1 <= nref <= veclen
270 ref = array of pointers to reference vectors, so that
271 ref[i][k] is the k-th component of the i-th reference,
272 for i=0..nref-1, k=0..veclen-1
273
274 Return value is a (double *) which points to malloc-ed memory
275 for the Choleksi decomposition. It should later be passed
276 to delayed_lsqfit. If wt != NULL, then ref[ii][jj] is
277 scaled by wt[jj] as well.
278
279 If NULL is returned, an error occured.
280 ------------------------------------------------------------------*/
281
startup_lsqfit(int veclen,float * wt,int nref,float * ref[])282 double * startup_lsqfit( int veclen ,
283 float *wt , int nref , float *ref[] )
284 {
285 int ii , jj , kk ;
286 double *cc = NULL ;
287 double sum ;
288
289 if( nref < 1 || veclen < nref || ref == NULL ){
290 ERROR_message("startup_lsqfit: nref=%d veclen=%d ref=%p",
291 nref , veclen , (void *)ref ) ;
292 return NULL ;
293 }
294
295 /* 13 Mar 2017 */
296 { int nbad = check_ref_vectors( veclen , nref , ref , "lsqfit" ) ;
297 if( nbad ) return NULL ;
298 }
299
300 /*** form coefficient matrix into cc */
301
302 cc = DBLEVEC(nref*nref) ;
303 if( cc == NULL ){
304 fprintf(stderr,"Can't malloc workspace in startup_lsqfit\n") ;
305 return NULL ;
306 }
307
308 if( wt != NULL ){
309 for( jj=0 ; jj < nref ; jj++ ){
310 for( ii=0 ; ii <= jj ; ii++ ){
311 sum = 0.0 ;
312 for( kk=0 ; kk < veclen ; kk++ ) sum += ref[ii][kk] * ref[jj][kk] * wt[kk] ;
313 CC(ii,jj) = CC(jj,ii) = sum ;
314 }
315 }
316 } else {
317 for( jj=0 ; jj < nref ; jj++ ){
318 for( ii=0 ; ii <= jj ; ii++ ){
319 sum = 0.0 ;
320 for( kk=0 ; kk < veclen ; kk++ ) sum += ref[ii][kk] * ref[jj][kk] ;
321 CC(ii,jj) = CC(jj,ii) = sum ;
322 }
323 }
324 }
325
326 /*** Choleski decompose cc ***/
327
328 for( ii=0 ; ii < nref ; ii++ ){
329 for( jj=0 ; jj < ii ; jj++ ){
330 sum = CC(ii,jj) ;
331 for( kk=0 ; kk < jj ; kk++ ) sum -= CC(ii,kk) * CC(jj,kk) ;
332 CC(ii,jj) = sum / CC(jj,jj) ;
333 }
334 sum = CC(ii,ii) ;
335 for( kk=0 ; kk < ii ; kk++ ) sum -= CC(ii,kk) * CC(ii,kk) ;
336 if( sum <= 0.0 ){
337 free(cc) ;
338 ERROR_message(
339 "Choleski factorization failure in startup_lsqfit [row=%d sum=%g]",ii,sum) ;
340 return NULL ;
341 }
342 CC(ii,ii) = sqrt(sum) ;
343 }
344
345 /*** scale ref by wt, if desired ***/
346
347 if( wt != NULL ){
348 for( ii=0 ; ii < nref ; ii++ )
349 for( jj=0 ; jj < veclen ; jj++ ) ref[ii][jj] *= wt[jj] ;
350 }
351
352 return cc ;
353 }
354
355 /**------------------------------------------------------------------
356 Given the data from startup_lsqfit, finish the job.
357 ------------------------------------------------------------------**/
358
delayed_lsqfit(int veclen,float * data,int nref,float * ref[],double * cc)359 float * delayed_lsqfit( int veclen ,
360 float * data , int nref , float *ref[] , double * cc )
361 {
362 int ii , jj ;
363 float *alpha = NULL ;
364 double *rr = NULL ;
365 register double sum ;
366
367 if( nref < 1 || veclen < nref ||
368 data == NULL || ref == NULL || cc == NULL ) return NULL ;
369
370 /* 13 Mar 2017 */
371 { int nbad = check_ref_vectors( veclen , nref , ref , "lsqfit" ) ;
372 if( nbad ) return NULL ;
373 }
374
375 /*** form RHS vector into rr ***/
376
377 rr = DBLEVEC(nref) ; if( rr == NULL ) return NULL ;
378
379 for( ii=0 ; ii < nref ; ii++ ){
380 sum = 0.0 ;
381 for( jj=0 ; jj < veclen ; jj++ ) sum += ref[ii][jj] * data[jj] ;
382 rr[ii] = sum ;
383 }
384
385 /*** forward solve ***/
386
387 for( ii=0 ; ii < nref ; ii++ ){
388 sum = rr[ii] ;
389 for( jj=0 ; jj < ii ; jj++ ) sum -= CC(ii,jj) * rr[jj] ;
390 rr[ii] = sum / CC(ii,ii) ;
391 }
392
393 /*** backward solve ***/
394
395 for( ii=nref-1 ; ii >= 0 ; ii-- ){
396 sum = rr[ii] ;
397 for( jj=ii+1 ; jj < nref ; jj++ ) sum -= CC(jj,ii) * rr[jj] ;
398 rr[ii] = sum / CC(ii,ii) ;
399 }
400
401 /*** put result into alpha ***/
402
403 alpha = (float *) malloc( sizeof(float) * nref ) ; if( alpha == NULL ) return NULL ;
404 for( ii=0 ; ii < nref ; ii++ ) alpha[ii] = rr[ii] ;
405
406 /*** cleanup and exit ***/
407
408 free(rr) ;
409 return alpha ;
410 }
411
412 /*-----------------------------------------------------------------
413 'startup' and 'delayed' versions for image fitting.
414 N.B.: unlike mri_lsqfit, all the images in refim and wtim must
415 be of MRI_float kind.
416 -------------------------------------------------------------------*/
417
mri_startup_lsqfit(MRI_IMARR * refim,MRI_IMAGE * wtim)418 double * mri_startup_lsqfit( MRI_IMARR * refim , MRI_IMAGE * wtim )
419 {
420 double *cc = NULL ; /* will be the output */
421 int ii , npix,nref ;
422 float * wtar , ** refar ;
423
424 ENTRY("mri_startup_lsqfit") ;
425
426 /****---- check inputs ----****/
427
428 if( wtim != NULL && wtim->kind != MRI_float ){
429 ERROR_message("mri_startup_lsqfit: non-float wtim") ; RETURN(NULL);
430 }
431 wtar = MRI_FLOAT_PTR(wtim) ;
432
433 if( refim == NULL || refim->num < 1 ){
434 ERROR_message("mri_startup_lsqfit: NULL refim") ; RETURN(NULL);
435 }
436
437 nref = refim->num ;
438 npix = refim->imarr[0]->nvox ;
439 refar = (float **) malloc( sizeof(float *) * nref ) ;
440 if( refar == NULL ){
441 fprintf(stderr,"mri_startup_lsqfit: malloc failure for refar") ;
442 RETURN(NULL);
443 }
444
445 for( ii=0 ; ii < nref ; ii++ ){
446 if( refim->imarr[ii] == NULL ){
447 ERROR_message("mri_startup_lsqfit: NULL refim[%d]",ii) ; RETURN(NULL);
448 }
449 if( refim->imarr[ii]->nvox != npix ){
450 ERROR_message("mri_startup_lsqfit: MISMATCH refim[%d]",ii) ; RETURN(NULL);
451 }
452 if( refim->imarr[ii]->kind != MRI_float ){
453 ERROR_message("mri_startup_lsqfit: non-float refim[%d]",ii) ; RETURN(NULL);
454 }
455 refar[ii] = MRI_FLOAT_PTR(refim->imarr[ii]) ;
456 }
457
458 /****---- get Choleski, send it out ----****/
459
460 cc = startup_lsqfit( npix , wtar , nref , refar ) ;
461 if( cc == NULL ){
462 ERROR_message("mri_startup_lsqfit: bad call to startup_lsqfit!\a\n") ;
463 RETURN(NULL);
464 }
465 free(refar) ;
466 RETURN(cc) ;
467 }
468
469 /*-----------------------------------------------------------------*/
470
mri_delayed_lsqfit(MRI_IMAGE * fitim,MRI_IMARR * refim,double * cc)471 float * mri_delayed_lsqfit( MRI_IMAGE * fitim , MRI_IMARR * refim , double * cc )
472 {
473 int ii , npix,nref ;
474 float * fit ;
475 static float ** refar = NULL ;
476 static int nrefar = -1 ;
477
478 ENTRY("mri_delayed_lsqfit") ;
479
480 nref = refim->num ;
481 npix = refim->imarr[0]->nvox ;
482
483 if( nrefar < nref ){
484 if( refar != NULL ) free(refar) ;
485 refar = (float **) malloc( sizeof(float *) * nref ) ;
486 nrefar = nref ;
487 }
488 if( refar == NULL ){
489 fprintf(stderr,"mri_delayed_lsqfit: malloc failure for refar");
490 RETURN(NULL);
491 }
492
493 for( ii=0 ; ii < nref ; ii++ )
494 refar[ii] = MRI_FLOAT_PTR(refim->imarr[ii]) ;
495
496 fit = delayed_lsqfit( npix , MRI_FLOAT_PTR(fitim) , nref , refar , cc ) ;
497 RETURN(fit) ;
498 }
499