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 /*-------------------------------------------------------------------------*/
10
THD_is_zero(int npt,float * xx)11 int THD_is_zero( int npt , float *xx ) /* 20 May 2011 */
12 {
13 int ii ;
14
15 if( npt <= 0 || xx == NULL ) return 1 ;
16
17 for( ii=0 ; ii < npt ; ii++ ) if( xx[ii] != 0.0f ) return 0 ;
18 return 1 ;
19 }
20
21 /*-------------------------------------------------------------------------*/
22
THD_is_constant(int npt,float * xx)23 int THD_is_constant( int npt , float *xx ) /* 20 May 2011 */
24 {
25 int ii ; float val ;
26
27 if( npt <= 1 || xx == NULL ) return 1 ;
28
29 val = xx[0] ;
30 for( ii=1 ; ii < npt ; ii++ ) if( xx[ii] != val ) return 0 ;
31 return 1 ;
32 }
33
34 /*-------------------------------------------------------------------------*/
35
THD_const_detrend(int npt,float * xx,float * xx0)36 void THD_const_detrend( int npt, float *xx, float *xx0 ) /* 24 Aug 2001 */
37 {
38 int ii ; float xbar ;
39
40 if( npt < 2 || xx == NULL ) return ;
41
42 xbar = 0.0 ;
43 for( ii=0 ; ii < npt ; ii++ ) xbar += xx[ii] ;
44 xbar /= npt ;
45 for( ii=0 ; ii < npt ; ii++ ) xx[ii] -= xbar ;
46 if( xx0 != NULL ) *xx0 = xbar ;
47 return ;
48 }
49
50 /*-------------------------------------------------------------------------*/
51
get_linear_trend(int npt,float * xx,float * f0,float * f1)52 void get_linear_trend( int npt, float *xx, float *f0, float *f1 )
53 {
54 double t1,t3,t10 , x0,x1 ;
55 int ii ;
56
57 if( npt < 2 || xx == NULL || f0 == NULL || f1 == NULL ) return ;
58
59 x0 = xx[0] ; x1 = 0.0 ;
60 for( ii=1 ; ii < npt ; ii++ ){
61 x0 += xx[ii] ;
62 x1 += xx[ii] * ii ;
63 }
64
65 t1 = npt*x0; t3 = 1.0/npt; t10 = npt*npt;
66
67 *f0 = (float)(2.0/(npt+1.0)*t3*(2.0*t1-3.0*x1-x0));
68 *f1 = (float)(-6.0/(t10-1.0)*t3*(-x0-2.0*x1+t1));
69 return ;
70 }
71
72 /*-------------------------------------------------------------------------
73 Linear detrend a 1D float array in place.
74 If xx0 != NULL and xx1 != NULL, then the trend removed is
75 far[i] -= (*xx0) + (*xx1) * i, for i=0..npt-1
76 ---------------------------------------------------------------------------*/
77
THD_linear_detrend(int npt,float * far,float * xx0,float * xx1)78 void THD_linear_detrend( int npt, float *far, float *xx0, float *xx1 )
79 {
80 int ii ;
81 float f0=0.0, f1=0.0;
82
83 if( npt < 3 || far == NULL ) return ;
84
85 get_linear_trend( npt , far , &f0 , &f1 ) ;
86
87 far[0] -= f0 ;
88 for( ii=1 ; ii < npt ; ii++ ) far[ii] -= (f0 + f1*ii) ;
89
90 if( xx0 != NULL ) *xx0 = f0 ;
91 if( xx1 != NULL ) *xx1 = f1 ;
92
93 return ;
94 }
95
96 /*---------------------------------------------------------------------------*/
97
THD_linear_detrend_complex(int npt,complex * cx)98 void THD_linear_detrend_complex( int npt , complex *cx ) /* 05 Mar 2007 */
99 {
100 float *f ; int ii ;
101
102 if( npt < 3 || cx == NULL ) return ;
103
104 f = malloc(sizeof(float)*npt) ;
105 for( ii=0 ; ii < npt ; ii++ ) f[ii] = cx[ii].r ;
106 THD_linear_detrend( npt , f , NULL,NULL ) ;
107 for( ii=0 ; ii < npt ; ii++ ){
108 cx[ii].r = f[ii] ; f[ii] = cx[ii].i ;
109 }
110 THD_linear_detrend( npt , f , NULL,NULL ) ;
111 for( ii=0 ; ii < npt ; ii++ ) cx[ii].i = f[ii] ;
112 return ;
113 }
114
115 /*---------------------------------------------------------------------------
116 Given x[0..npt-1], return f0,f1,f2 as the least squares coefficients to
117 x[j] = f0 + f1*j + f2*j*j
118 -----------------------------------------------------------------------------*/
119
get_quadratic_trend(int npt,float * xx,float * f0,float * f1,float * f2)120 void get_quadratic_trend( int npt, float *xx, float *f0, float *f1, float *f2 )
121 {
122 double x0,x1,x2 ; double N=npt ;
123 int ii ;
124
125 if( npt < 3 || xx == NULL || f0 == NULL || f1 == NULL || f2 == NULL ) return;
126
127 x0 = xx[0] ; x1 = x2 = 0.0 ;
128 for( ii=1 ; ii < npt ; ii++ ){
129 x0 += xx[ii] ;
130 x1 += (xx[ii] * ii) ;
131 x2 += (xx[ii] * ii) * ii ;
132 }
133
134 *f0 = ( 3.0*(3.0*N*N-3.0*N+2.0) * x0
135 -18.0*(-1.0+2.0*N) * x1
136 +30.0 * x2 ) / (N*(N+2.0)*(N+1.0)) ;
137
138 *f1 = ( -18.0*(-1.0+2.0*N) * x0
139 +12.0*(-1.0+2.0*N)*(8.0*N-11.0) * x1 /((N-1.0)*(N-2.0))
140 -180.0 * x2 /(N-2.0) )
141 / (N*(N+2.0)*(N+1.0)) ;
142
143 *f2 = ( 30.0 * x0
144 -180.0 * x1 / (N-2.0)
145 +180.0 * x2 / ((N-1.0)*(N-2.0)) ) / (N*(N+2.0)*(N+1.0)) ;
146 return ;
147 }
148
149 /*-------------------------------------------------------------------------
150 Quadratic detrend a 1D float array in place.
151 If xx0 != NULL, xx1 != NULL, xx2 != NULL, then the trend removed is
152 far[i] -= (*xx0) + (*xx1) * i + (*xx2)*(i*i) , for i=0..npt-1
153 ---------------------------------------------------------------------------*/
154
THD_quadratic_detrend(int npt,float * far,float * xx0,float * xx1,float * xx2)155 void THD_quadratic_detrend( int npt, float *far,
156 float *xx0, float *xx1, float *xx2 )
157 {
158 int ii ;
159 float f0 , f1 , f2 ;
160
161 if( npt < 4 || far == NULL ) return ;
162
163 get_quadratic_trend( npt , far , &f0 , &f1 , &f2 ) ;
164
165 far[0] -= f0 ;
166 for( ii=1 ; ii < npt ; ii++ ) far[ii] -= ( (f2*ii + f1)*ii + f0 ) ;
167
168 if( xx0 != NULL ) *xx0 = f0 ;
169 if( xx1 != NULL ) *xx1 = f1 ;
170 if( xx2 != NULL ) *xx2 = f2 ;
171
172 return ;
173 }
174
175 /*------------------------------------------------------------------------------*/
176 /*! Cubic detrend a float array in place. */
177
THD_cubic_detrend(int npt,float * far)178 void THD_cubic_detrend( int npt , float *far ) /* 15 Nov 1999 */
179 {
180 int ii ;
181 float g0,g1,g2,g3 , f0,f1,f2,f3 , t1,t2,t5,t8 , t95,t56,t22,t25,txx ;
182
183 if( npt < 5 || far == NULL ) return ;
184
185 t8 = npt*npt ; t2 = npt-1.0 ; t5 = t2*(npt-2.0) ;
186 t95 = 0.05*t5*(npt-3.0) ;
187 t56 = 0.16666667*t5 ;
188 t22 = 0.5*t2 ;
189 t25 = 1.5*t2 ;
190 txx = 0.6*t8-1.5*npt+1.1 ;
191
192 g0=g1=g2=g3=0.0 ;
193 for( ii=0 ; ii < npt ; ii++ ){
194 t1 = ii*ii ;
195 f1 = ii - t22 ;
196 f2 = t1 - t2*ii + t56 ;
197 f3 = t1*(ii - t25) + txx*ii - t95 ;
198
199 g0 += far[ii] ;
200 g1 += far[ii] * f1 ;
201 g2 += far[ii] * f2 ;
202 g3 += far[ii] * f3 ;
203 }
204 g0 *= (1.0/npt) ;
205 g1 *= 12.0/(npt*(t8-1.0)) ;
206 g2 *= 180.0/(npt*(t8-1.0)*(t8-4.0)) ;
207 g3 *= 2800.0/(npt*(t8-1.0)*(t8-4.0)*(t8-9.0)) ;
208
209 for( ii=0 ; ii < npt ; ii++ ){
210 t1 = ii*ii ;
211 f1 = ii- t22 ;
212 f2 = t1 - t2*ii + t56 ;
213 f3 = t1*(ii - t25) + txx*ii - t95 ;
214
215 far[ii] -= ( g0 + g1*f1 + g2*f2 + g3*f3 ) ;
216 }
217
218 return ;
219 }
220
221 /*-------------------------------------------------------------------------
222 Make a vector have L2 norm 1
223 ---------------------------------------------------------------------------*/
224
THD_normalize(int npt,float * far)225 float THD_normalize( int npt , float *far )
226 {
227 int ii ;
228 float fac ;
229
230 if( npt <= 0 || far == NULL ) return 0.0f ;
231
232 fac = 0.0f ;
233 for( ii=0 ; ii < npt ; ii++ ) fac += far[ii]*far[ii] ;
234 if( fac <= 1.e-20f ) return 0.0f ;
235 fac = 1.0f / sqrtf(fac) ;
236 for( ii=0 ; ii < npt ; ii++ ) far[ii] *= fac ;
237 return fac ;
238 }
239
240 /*-------------------------------------------------------------------------
241 Make a vector have RMS value = 1
242 ---------------------------------------------------------------------------*/
243
THD_normRMS(int npt,float * far)244 void THD_normRMS( int npt , float *far )
245 {
246 int ii ;
247 float fac ;
248
249 if( npt <= 0 || far == NULL ) return ;
250
251 fac = 0.0f ;
252 for( ii=0 ; ii < npt ; ii++ ) fac += far[ii]*far[ii] ;
253 if( fac == 0.0f ) return ;
254 fac = 1.0f / sqrtf(fac/npt) ;
255 for( ii=0 ; ii < npt ; ii++ ) far[ii] *= fac ;
256 return ;
257 }
258
259 /*-----------------------------------------------------------------------*/
260 /* Make a vector have max abs 1 [26 Mar 2008]. */
261
THD_normmax(int npt,float * far)262 void THD_normmax( int npt , float *far )
263 {
264 int ii ;
265 float fac , val ;
266
267 if( npt <= 0 || far == NULL ) return ;
268 fac = 0.0f ;
269 for( ii=0 ; ii < npt ; ii++ ){ val = fabsf(far[ii]); fac = MAX(fac,val); }
270 if( fac == 0.0f ) return ;
271 fac = 1.0f / fac ;
272 for( ii=0 ; ii < npt ; ii++ ) far[ii] *= fac ;
273 return ;
274 }
275
276 /*-----------------------------------------------------------------------*/
277 /* Make a vector have max sum |far| = 1 [26 Mar 2008]. */
278
THD_normL1(int npt,float * far)279 void THD_normL1( int npt , float *far )
280 {
281 int ii ;
282 float fac , val ;
283
284 if( npt <= 0 || far == NULL ) return ;
285 fac = 0.0f ;
286 for( ii=0 ; ii < npt ; ii++ ){ val = fabsf(far[ii]); fac += val; }
287 if( fac == 0.0f ) return ;
288 fac = 1.0f / fac ;
289 for( ii=0 ; ii < npt ; ii++ ) far[ii] *= fac ;
290 return ;
291 }
292
293 /*-----------------------------------------------------------------------*/
294 /*! Detrend a vector with a given polort level, plus some others, using
295 least squares regression.
296 - npt = length of vector
297 - far = vector of data (will be modified)
298 - polort = polynomial order
299 - nort = number of extra orts in ort[] (can be 0)
300 - ort = array of extra time series to detrend:
301 ort[j][i] for j=0..nort-1, i=0..npt-1
302 - fit = array of length nref=polort+nort+1 to hold parameters of fit
303 (can be NULL)
304 -------------------------------------------------------------------------*/
305
THD_generic_detrend_LSQ(int npt,float * far,int polort,int nort,float ** ort,float * fit)306 void THD_generic_detrend_LSQ( int npt, float *far ,
307 int polort, int nort, float **ort , float *fit )
308 {
309 int ii,jj , nref ;
310 float **ref , *qfit , xmid , xfac , val ;
311
312 ENTRY("THD_generic_detrend_LSQ") ;
313
314 /* check inputs */
315
316 if( npt <= 1 || far == NULL ){
317 ERROR_message("THD_generic_detrend_LSQ :: bad inputs"); EXRETURN;
318 }
319 if( nort > 0 ){
320 if( ort == NULL ){
321 ERROR_message("THD_generic_detrend_LSQ :: ort == NULL"); EXRETURN;
322 }
323 for( jj=0 ; jj < nort ; jj++ ){
324 if( ort[jj] == NULL || THD_is_zero(npt,ort[jj]) ){
325 ERROR_message("THD_generic_detrend_LSQ :: ort[%d] is invalid",jj) ;
326 EXRETURN ;
327 }
328 }
329 }
330 if( polort < 0 ) polort = -1 ;
331 if( nort < 0 ) nort = 0 ;
332
333 nref = polort+1+nort ;
334 if( nref <= 0 || nref >= npt-1 ) EXRETURN ;
335
336 /* assemble all reference vectors */
337
338 ref = (float **) malloc( sizeof(float *) * nref ) ;
339 xmid = 0.5*(npt-1) ; xfac = 1.0 / xmid ;
340 for( jj=0 ; jj <= polort ; jj++ ){
341 ref[jj] = (float *) malloc( sizeof(float) * npt ) ;
342 for( ii=0 ; ii < npt ; ii++ )
343 ref[jj][ii] = (float)Plegendre(xfac*(ii-xmid),jj) ;
344 }
345 for( jj=0 ; jj < nort ; jj++ ) /* user supplied refs */
346 ref[polort+1+jj] = ort[jj] ;
347
348 qfit = lsqfit( npt , far , NULL , nref , ref ) ; /* actual fitting */
349
350 if( qfit != NULL ){ /* good */
351 for( ii=0 ; ii < npt ; ii++ ){
352 val = far[ii] ;
353 for( jj=0 ; jj < nref ; jj++ ) val -= qfit[jj] * ref[jj][ii] ;
354 far[ii] = val ;
355 }
356 if( fit != NULL ){ for( ii=0 ; ii < nref ; ii++ ) fit[ii] = qfit[ii] ; }
357 free(qfit) ;
358 } else {
359 ERROR_message("THD_generic_detrend_LSQ: fit fails - no detrending!") ;
360 if( fit != NULL ) memset(fit,0,sizeof(float)*nref) ;
361 }
362
363 for( jj=0 ; jj <= polort ; jj++ ) free(ref[jj]) ;
364 free(ref) ; EXRETURN ;
365 }
366
367 /*-----------------------------------------------------------------------*/
368 /*! Detrend a vector with a given polort level, plus some others, using
369 L1 regression.
370 - npt = length of vector
371 (if < 0, then abs(npt) is used, and the fit isn't subtracted)
372 - far = vector of data (will be modified)
373 - polort = polynomial order
374 - nort = number of extra orts in ort[] (can be 0)
375 - ort = array of extra time series to detrend:
376 ort[j][i] for j=0..nort-1, i=0..npt-1
377 - fit = array of length nref=polort+nort+1 to hold parameters of fit
378 (can be NULL)
379 -------------------------------------------------------------------------*/
380
THD_generic_detrend_L1(int npt,float * far,int polort,int nort,float ** ort,float * fit)381 void THD_generic_detrend_L1( int npt, float *far ,
382 int polort, int nort, float **ort , float *fit )
383 {
384 int ii,jj , nref , nosub=0 ;
385 float **ref , *qfit , xmid , xfac , val ;
386
387 ENTRY("THD_generic_detrend_L1") ;
388
389 /* check inputs */
390
391 if( npt < -1 ){ nosub = 1 ; npt = -npt ; } /* 02 Jan 2013 */
392
393 if( npt <= 1 || far == NULL ){
394 ERROR_message("THD_generic_detrend_L1 :: bad inputs"); EXRETURN;
395 }
396 if( nort > 0 ){
397 if( ort == NULL ){
398 ERROR_message("THD_generic_detrend_LSQ :: ort == NULL"); EXRETURN;
399 }
400 for( jj=0 ; jj < nort ; jj++ ){
401 if( ort[jj] == NULL || THD_is_zero(npt,ort[jj]) ){
402 ERROR_message("THD_generic_detrend_L1 :: ort[%d] is invalid",jj) ;
403 EXRETURN ;
404 }
405 }
406 }
407 if( polort < 0 ) polort = -1 ;
408 if( nort < 0 ) nort = 0 ;
409
410 nref = polort+1+nort ;
411 if( nref <= 0 || nref >= npt-1 ) EXRETURN ;
412
413 /* assemble all reference vectors */
414
415 ref = (float **) malloc( sizeof(float *) * nref ) ;
416 xmid = 0.5*(npt-1) ; xfac = 1.0 / xmid ;
417 for( jj=0 ; jj <= polort ; jj++ ){
418 ref[jj] = (float *) malloc( sizeof(float) * npt ) ;
419 for( ii=0 ; ii < npt ; ii++ )
420 ref[jj][ii] = (float)Plegendre(xfac*(ii-xmid),jj) ;
421 }
422 for( jj=0 ; jj < nort ; jj++ ) /* user supplied refs */
423 ref[polort+1+jj] = ort[jj] ;
424
425 qfit = (float *)malloc(sizeof(float)*nref) ;
426 val = cl1_solve( npt , nref , far , ref , qfit , 0 ) ;
427
428 if( val >= 0.0f ){ /* good */
429 if( !nosub ){
430 for( ii=0 ; ii < npt ; ii++ ){
431 val = far[ii] ;
432 for( jj=0 ; jj < nref ; jj++ ) val -= qfit[jj] * ref[jj][ii] ;
433 far[ii] = val ;
434 }
435 }
436 if( fit != NULL ){ for( ii=0 ; ii < nref ; ii++ ) fit[ii] = qfit[ii] ; }
437 } else {
438 ERROR_message("THD_generic_detrend_L1: fit fails - no detrending!") ;
439 if( fit != NULL ) memset(fit,0,sizeof(float)*nref) ;
440 }
441 free(qfit) ;
442
443 for( jj=0 ; jj <= polort ; jj++ ) free(ref[jj]) ;
444 free(ref) ; EXRETURN ;
445 }
446
447 /*-----------------------------------------------------------------------*/
448 /* Add the trend back into far; fit must not be NULL!
449 -------------------------------------------------------------------------*/
450
THD_generic_retrend(int npt,float * far,int polort,int nort,float ** ort,float * fit)451 void THD_generic_retrend( int npt , float *far ,
452 int polort, int nort, float **ort , float *fit )
453 {
454 int ii,jj , nref ;
455 float **ref , xmid , xfac , val ;
456
457 ENTRY("THD_generic_retrend") ;
458
459 /* check inputs */
460
461 if( npt <= 1 || far == NULL ){
462 ERROR_message("THD_generic_retrend :: bad inputs"); EXRETURN;
463 }
464 if( nort > 0 ){
465 if( ort == NULL ){
466 ERROR_message("THD_generic_detrend_LSQ :: ort == NULL"); EXRETURN;
467 }
468 for( jj=0 ; jj < nort ; jj++ ){
469 if( ort[jj] == NULL || THD_is_zero(npt,ort[jj]) ){
470 ERROR_message("THD_generic_retrend :: ort[%d] is invalid",jj) ;
471 EXRETURN ;
472 }
473 }
474 }
475 if( polort < 0 ) polort = -1 ;
476 if( nort < 0 ) nort = 0 ;
477
478 nref = polort+1+nort ;
479 if( nref == 0 || nref >= npt-1 ) EXRETURN ;
480
481 /* assemble all reference vectors */
482
483 ref = (float **) malloc( sizeof(float *) * nref ) ;
484 xmid = 0.5*(npt-1) ; xfac = 1.0 / xmid ;
485 for( jj=0 ; jj <= polort ; jj++ ){
486 ref[jj] = (float *) malloc( sizeof(float) * npt ) ;
487 switch( jj ){
488 case 0:
489 for( ii=0 ; ii < npt ; ii++ ) ref[jj][ii] = 1.0 ;
490 break ;
491
492 case 1:
493 for( ii=0 ; ii < npt ; ii++ ) ref[jj][ii] = xfac*(ii-xmid) ;
494 break ;
495
496 case 2:
497 for( ii=0 ; ii < npt ; ii++ ){
498 val = xfac*(ii-xmid) ; ref[jj][ii] = val*val ;
499 }
500 break ;
501
502 case 3:
503 for( ii=0 ; ii < npt ; ii++ ){
504 val = xfac*(ii-xmid) ; ref[jj][ii] = val*val*val ;
505 }
506 break ;
507
508 default:
509 for( ii=0 ; ii < npt ; ii++ ){
510 val = xfac*(ii-xmid) ; ref[jj][ii] = pow(val,(double)(jj)) ;
511 }
512 break ;
513 }
514 }
515 for( jj=0 ; jj < nort ; jj++ ) /* user supplied refs */
516 ref[polort+1+jj] = ort[jj] ;
517
518 for( ii=0 ; ii < npt ; ii++ ){
519 val = far[ii] ;
520 for( jj=0 ; jj < nref ; jj++ ) val += fit[jj] * ref[jj][ii] ;
521 far[ii] = val ;
522 }
523
524 for( jj=0 ; jj <= polort ; jj++ ) free(ref[jj]) ;
525 free(ref) ; EXRETURN ;
526 }
527
528 /*----------------------------------------------------------------------------*/
529 /*! Returns nref reference functions. */
530
THD_build_polyref(int nref,int nvals)531 float ** THD_build_polyref( int nref , int nvals ) /* 20 Sep 2007 */
532 {
533 int jj,iv,kk ;
534 float **ref ;
535 double fac , x ;
536
537 ENTRY("THD_build_polyref") ;
538
539 if( nref < 1 || nvals <= nref ){
540 ERROR_message("THD_build_polyref: nref=%d nvals=%d",nref,nvals) ;
541 RETURN(NULL) ;
542 }
543
544 fac = 2.0 / (nvals-1.0) ;
545 ref = (float **)malloc(sizeof(float *)*nref) ;
546 for( jj=0 ; jj < nref ; jj++ ){
547 ref[jj] = (float *) malloc( sizeof(float) * nvals ) ;
548 for( kk=0 ; kk < nvals ; kk++ ){
549 x = fac * kk - 1.0 ; ref[jj][kk] = (float)Plegendre( x , jj ) ;
550 }
551 }
552
553 RETURN(ref) ;
554 }
555
556 /*----------------------------------------------------------------------------*/
557 /*! Returns 2*corder+3 reference functions. */
558
THD_build_trigref(int corder,int nvals)559 float ** THD_build_trigref( int corder , int nvals )
560 {
561 int nref=2*corder+3 , jj,iv,kk ;
562 float **ref , tm,fac,fq ;
563
564 ENTRY("THD_build_trigref") ;
565
566 if( corder < 0 || nvals <= nref ){
567 ERROR_message("THD_build_trigref: corder=%d nvals=%d",corder,nvals) ;
568 RETURN(NULL) ;
569 }
570
571 ref = (float **)malloc(sizeof(float *)*nref) ;
572 for( jj=0 ; jj < nref ; jj++ )
573 ref[jj] = (float *) malloc( sizeof(float) * nvals ) ;
574
575 /* ref[0]: r(t) = 1 */
576 for( iv=0 ; iv < nvals ; iv++ ) ref[0][iv] = 1.0f ;
577
578 /* ref[1]: r(t) = t - tmid */
579 tm = 0.5f * (nvals-1.0f) ; fac = 2.0f / nvals ;
580 for( iv=0 ; iv < nvals ; iv++ ) ref[1][iv] = (iv-tm)*fac ;
581
582 /* ref[2]: r(t) = (t-tmid)**2 */
583 fac = fac*fac ;
584 for( iv=0 ; iv < nvals ; iv++ ) ref[2][iv] = (iv-tm)*(iv-tm)*fac ;
585
586 /* higher order ref[jj]: */
587 for( jj=3,kk=1 ; kk <= corder ; kk++ ){
588 fq = (2.0*PI*kk)/nvals ;
589
590 /* r(t) = sin(2*PI*k*t/N) */
591 for( iv=0 ; iv < nvals ; iv++ ) ref[jj][iv] = sin(fq*iv) ;
592 jj++ ;
593
594 /* r(t) = cos(2*PI*k*t/N) */
595 for( iv=0 ; iv < nvals ; iv++ ) ref[jj][iv] = cos(fq*iv) ;
596 jj++ ;
597 }
598
599 RETURN(ref) ;
600 }
601
602 /*----------------------------------------------------------------------------*/
603 #undef INMASK /* ZSS: Changed from */
604 #define INMASK(i) ((mask == NULL || mask[i])) /* mask != NULL */
605 /*----------------------------------------------------------------------------*/
606 /* Fits each voxel time series to a linear model:
607
608 nref-1
609 z[i][t] = SUM ref[q][t] * fit[q][i]
610 q=0
611
612 where
613 - i = voxel index 0..DSET_NVOX(dset)-1
614 - t = time index 0..DSET_NVALS(dset)-1
615 - q = reference index 0..nref-1
616 - Fitting is done in the L-p sense, where p=meth (1 or 2).
617 - The dataset itself is not modified.
618 _ Return value is the nref fit images, plus 1 extra image that is
619 the MAD of the residuals for meth=1 and the standard deviation for meth=2.
620 - In float format, of course.
621 - If NULL is returned, something bad happened. Be afraid.
622 - Also see function THD_medmad_bricks() for something simpler & similar.
623 ------------------------------------------------------------------------------*/
624
THD_time_fit_dataset(THD_3dim_dataset * dset,int nref,float ** ref,int meth,byte * mask)625 MRI_IMARR * THD_time_fit_dataset( THD_3dim_dataset *dset ,
626 int nref , float **ref , int meth , byte *mask )
627 {
628 int ii , nvox,nval , qq,tt ;
629 float *far , *fit , *var , val ;
630 MRI_IMARR *imar ; MRI_IMAGE *qim ; float **fitar ;
631
632 ENTRY("THD_time_fit_dataset") ;
633
634 if( !ISVALID_DSET(dset) ||
635 nref < 1 || nref >= DSET_NVALS(dset) || ref == NULL ) RETURN(NULL) ;
636 DSET_load(dset) ; if( !DSET_LOADED(dset) ) RETURN(NULL) ;
637
638 /* construct output images */
639
640 INIT_IMARR(imar) ;
641 fitar = (float **)malloc(sizeof(float *)*nref) ;
642 for( qq=0 ; qq < nref ; qq++ ){
643 qim = mri_new_conforming( DSET_BRICK(dset,0) , MRI_float ) ;
644 fitar[qq] = MRI_FLOAT_PTR(qim) ;
645 ADDTO_IMARR(imar,qim) ;
646 }
647 qim = mri_new_conforming( DSET_BRICK(dset,0) , MRI_float ) ;
648 var = MRI_FLOAT_PTR(qim) ; ADDTO_IMARR(imar,qim) ;
649
650 nvox = DSET_NVOX(dset) ; nval = DSET_NVALS(dset) ;
651 far = (float *)malloc(sizeof(float)*nval) ;
652 fit = (float *)malloc(sizeof(float)*nref) ;
653 for( ii=0 ; ii < nvox ; ii++ ){
654 if( !INMASK(ii) ) continue ;
655 qq = THD_extract_array( ii , dset , 0 , far ) ; /* get data */
656 if( qq == 0 ){
657 switch(meth){ /* get fit */
658 default:
659 case 2: THD_generic_detrend_LSQ( nval, far, -1, nref,ref, fit ); break;
660 case 1: THD_generic_detrend_L1 ( nval, far, -1, nref,ref, fit ); break;
661 }
662 for( qq=0 ; qq < nref ; qq++ ) fitar[qq][ii] = fit[qq] ; /* save fit */
663
664 /* at this point, far[] contains the residuals */
665
666 switch(meth){ /* get stdev or MAD */
667 default:
668 case 2:{
669 float mm,vv ;
670 for( mm=0.0,tt=0 ; tt < nval ; tt++ ) mm += far[tt] ;
671 mm /= nval ;
672 for( vv=0.0,tt=0 ; tt < nval ; tt++ ) vv += (far[tt]-mm)*(far[tt]-mm) ;
673 var[ii] = sqrtf( vv/(nval-1.0) ) ;
674 }
675 break ;
676
677 case 1:{
678 for( tt=0 ; tt < nval ; tt++ ) far[tt] = fabsf(far[tt]) ;
679 var[ii] = qmed_float( nval , far ) ;
680 }
681 break ;
682 }
683 }
684 }
685
686 free(fit); free(far); free(fitar); RETURN(imar);
687 }
688
689 /*----------------------------------------------------------------------------*/
690 /* Get the time series data array at voxel #ii, detrended by the fit
691 from THD_time_fit_dataset(). If scl != 0, also scale the result
692 by the reciprocal of the stdev or MAD (last sub-brick in imar).
693 ------------------------------------------------------------------------------*/
694
THD_extract_detrended_array(THD_3dim_dataset * dset,int nref,float ** ref,MRI_IMARR * imar,int ii,int scl,float * far)695 void THD_extract_detrended_array( THD_3dim_dataset *dset ,
696 int nref, float **ref, MRI_IMARR *imar,
697 int ii, int scl, float *far )
698 {
699 int tt , nval , qq ;
700 float val , **fitar , *var ;
701 MRI_IMAGE *qim ;
702
703 ENTRY("THD_extract_detrended_array") ;
704
705 if( !ISVALID_DSET(dset) ||
706 nref < 1 || ref == NULL ||
707 imar == NULL || IMARR_COUNT(imar) < nref+1 ||
708 ii < 0 || ii >= DSET_NVOX(dset) || far == NULL ) EXRETURN ;
709
710 qq = THD_extract_array( ii , dset , 0 , far ) ; /* get data */
711 if( qq < 0 ) EXRETURN ;
712 nval = DSET_NVALS(dset) ;
713
714 fitar = (float **)malloc(sizeof(float *)*nref) ;
715 for( qq=0 ; qq < nref ; qq++ ){
716 qim = IMARR_SUBIM(imar,qq) ; fitar[qq] = MRI_FLOAT_PTR(qim) ;
717 }
718 qim = IMARR_SUBIM(imar,nref) ; var = MRI_FLOAT_PTR(qim) ;
719
720 for( tt=0 ; tt < nval ; tt++ ){ /* get residuals */
721 val = far[tt] ;
722 for( qq=0 ; qq < nref ; qq++ ) val -= ref[qq][tt] * fitar[qq][ii] ;
723 far[tt] = val ;
724 }
725
726 if( scl && var[ii] > 0.0f ){
727 val = 1.0f / var[ii] ;
728 for( tt=0 ; tt < nval ; tt++ ) far[tt] *= val ;
729 }
730
731 /* ZSS: Need to free fitar */
732 free(fitar); fitar=NULL;
733 EXRETURN ;
734 }
735
736 /*----------------------------------------------------------------------------*/
737 /* Detrend a dataset:
738 - dset = input time series dataset (nvals long)
739 - nref = number of ref vectors
740 - ref = reference vectors [0..nref-1][0..nvals-1]
741 - meth = 1 or 2 (for L_p fitting, where p=meth)
742 - scl = if 1, scale data by reciprocal of MAD or stdev (meth=1 or 2)
743 - mask = if not NULL, byte mask; voxel outside the mask are set to 0
744 - imar = if not NULL, coefficient images from THD_time_fit_dataset()
745
746 Output is in float format, of course. If NULL is returned, something bad
747 transpired, and your mother will get an e-mail explaining how stupid
748 you were.
749 ------------------------------------------------------------------------------*/
750
THD_detrend_dataset(THD_3dim_dataset * dset,int nref,float ** ref,int meth,int scl,byte * mask,MRI_IMARR ** imar)751 THD_3dim_dataset * THD_detrend_dataset( THD_3dim_dataset *dset ,
752 int nref , float **ref ,
753 int meth , int scl ,
754 byte *mask , MRI_IMARR **imar )
755 {
756 MRI_IMARR *qmar ;
757 int ii,jj,kk , nvals,nvox , iv ;
758 float *var ;
759 THD_3dim_dataset *newset ;
760
761 ENTRY("THD_detrend_dataset") ;
762
763 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
764 nvals = DSET_NVALS(dset) ; nvox = DSET_NVOX(dset) ;
765
766 qmar = THD_time_fit_dataset( dset , nref,ref , meth , mask ) ;
767 if( qmar == NULL ) RETURN(NULL) ;
768
769 newset = EDIT_empty_copy(dset) ;
770 for( iv=0 ; iv < nvals ; iv++ ){
771 EDIT_substitute_brick( newset , iv , MRI_float , NULL ) ;
772 EDIT_BRICK_FACTOR( newset , iv , 0.0f ) ; /* 04 Jun 2007 */
773 }
774
775 var = (float *)malloc(sizeof(float)*nvals) ;
776 for( ii=0 ; ii < nvox ; ii++ ){
777 if( mask == NULL || mask[ii] )
778 THD_extract_detrended_array( dset , nref,ref , qmar , ii,scl , var ) ;
779 else
780 memset(var,0,sizeof(float)*nvals) ;
781 THD_insert_series( ii , newset , nvals , MRI_float , var , 0 ) ;
782 }
783
784 free(var) ;
785 if( imar != NULL ) *imar = qmar ;
786 else DESTROY_IMARR(qmar) ;
787
788 RETURN(newset) ;
789 }
790
791 /*----------------------------------------------------------------------------*/
792 /* Can be used as a sort-of inverse to THD_detrend_dataset() */
793
THD_retrend_dataset(THD_3dim_dataset * dset,int nref,float ** ref,int scl,byte * mask,MRI_IMARR * imar)794 int THD_retrend_dataset( THD_3dim_dataset *dset ,
795 int nref , float **ref ,
796 int scl , byte *mask , MRI_IMARR *imar )
797 {
798 int ii,qq,tt , nvals,nvox ;
799 MRI_IMAGE *qim ; float **fitar , *far , fac , *var , val ;
800
801 ENTRY("THD_retrend_dataset") ;
802
803 if( !ISVALID_DSET(dset) ||
804 nref < 1 || ref == NULL ||
805 imar == NULL || IMARR_COUNT(imar) <= nref ) RETURN(0) ;
806
807 nvals = DSET_NVALS(dset) ; nvox = DSET_NVOX(dset) ;
808
809 fitar = (float **)malloc(sizeof(float *)*nref) ;
810 for( qq=0 ; qq < nref ; qq++ ){
811 qim = IMARR_SUBIM(imar,qq) ; fitar[qq] = MRI_FLOAT_PTR(qim) ;
812 }
813 qim = IMARR_SUBIM(imar,nref) ; var = MRI_FLOAT_PTR(qim) ;
814
815 far = (float *)malloc(sizeof(float)*nvals) ;
816 for( ii=0 ; ii < nvox ; ii++ ){
817 if( mask != NULL && !mask[ii] ) continue ;
818 qq = THD_extract_array( ii , dset , 0 , far ) ; /* get data */
819 if( qq < 0 ) continue ;
820 fac = (scl) ? var[ii] : 1.0f ;
821 for( tt=0 ; tt < nvals ; tt++ ){ /* add fit back in */
822 val = far[tt] * fac ;
823 for( qq=0 ; qq < nref ; qq++ ) val += ref[qq][tt] * fitar[qq][ii] ;
824 far[tt] = val ;
825 }
826 THD_insert_series( ii , dset , nvals , MRI_float , far , 0 ) ;
827 }
828
829 free(far) ; free(fitar) ; RETURN(1) ;
830 }
831