1 #include "mrilib.h"
2
3 static void linear_filter_reflect( int ntap, float *wt, int npt, float *x ) ;
4 static void median5_filter_reflect( int npt, float *x ) ;
5
6 /*-------------------------------------------------------------------------*/
7
main(int argc,char * argv[])8 int main( int argc , char *argv[] )
9 {
10 THD_3dim_dataset *yset=NULL , *aset=NULL , *mset=NULL , *wset=NULL ;
11 MRI_IMAGE *fim=NULL, *qim,*tim, *pfim=NULL , *vim , *wim=NULL ;
12 float *flar , *qar,*tar, *par=NULL , *var , *war=NULL ;
13 MRI_IMARR *fimar=NULL ;
14 MRI_IMAGE *aim , *yim ; float *aar , *yar ;
15 int nt=0 , nxyz=0 , nvox=0 , nparam=0 , nqbase , polort=0 , ii,jj,kk,bb ;
16 byte *mask=NULL ; int nmask=0 , iarg ;
17 char *fname_out="-" ; /** equiv to stdout **/
18
19 float alpha=0.0f ;
20 int nfir =0 ; float firwt[5]={0.09f,0.25f,0.32f,0.25f,0.09f} ;
21 int nmed =0 ;
22 int nwt =0 ;
23
24 #define METHOD_C 3
25 #define METHOD_K 11
26 int method = METHOD_C ;
27
28 /**--- help the pitiful user? ---**/
29
30 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
31 printf(
32 "Usage: 3dInvFMRI [options]\n"
33 "Program to compute stimulus time series, given a 3D+time dataset\n"
34 "and an activation map (the inverse of the usual FMRI analysis problem).\n"
35 "-------------------------------------------------------------------\n"
36 "OPTIONS:\n"
37 "\n"
38 " -data yyy =\n"
39 " *OR* = Defines input 3D+time dataset [a non-optional option].\n"
40 " -input yyy =\n"
41 "\n"
42 " -map aaa = Defines activation map; 'aaa' should be a bucket dataset,\n"
43 " each sub-brick of which defines the beta weight map for\n"
44 " an unknown stimulus time series [also non-optional].\n"
45 "\n"
46 " -mapwt www = Defines a weighting factor to use for each element of\n"
47 " the map. The dataset 'www' can have either 1 sub-brick,\n"
48 " or the same number as in the -map dataset. In the\n"
49 " first case, in each voxel, each sub-brick of the map\n"
50 " gets the same weight in the least squares equations.\n"
51 " [default: all weights are 1]\n"
52 "\n"
53 " -mask mmm = Defines a mask dataset, to restrict input voxels from\n"
54 " -data and -map. [default: all voxels are used]\n"
55 "\n"
56 " -base fff = Each column of the 1D file 'fff' defines a baseline time\n"
57 " series; these columns should be the same length as\n"
58 " number of time points in 'yyy'. Multiple -base options\n"
59 " can be given.\n"
60 " -polort pp = Adds polynomials of order 'pp' to the baseline collection.\n"
61 " The default baseline model is '-polort 0' (constant).\n"
62 " To specify no baseline model at all, use '-polort -1'.\n"
63 "\n"
64 " -out vvv = Name of 1D output file will be 'vvv'.\n"
65 " [default = '-', which is stdout; probably not good]\n"
66 "\n"
67 " -method M = Determines the method to use. 'M' is a single letter:\n"
68 " -method C = least squares fit to data matrix Y [default]\n"
69 " -method K = least squares fit to activation matrix A\n"
70 "\n"
71 " -alpha aa = Set the 'alpha' factor to 'aa'; alpha is used to penalize\n"
72 " large values of the output vectors. Default is 0.\n"
73 " A large-ish value for alpha would be 0.1.\n"
74 "\n"
75 " -fir5 = Smooth the results with a 5 point lowpass FIR filter.\n"
76 " -median5 = Smooth the results with a 5 point median filter.\n"
77 " [default: no smoothing; only 1 of these can be used]\n"
78 "-------------------------------------------------------------------\n"
79 "METHODS:\n"
80 " Formulate the problem as\n"
81 " Y = V A' + F C' + errors\n"
82 " where Y = data matrix (N x M) [from -data]\n"
83 " V = stimulus (N x p) [to -out]\n"
84 " A = map matrix (M x p) [from -map]\n"
85 " F = baseline matrix (N x q) [from -base and -polort]\n"
86 " C = baseline weights (M x q) [not computed]\n"
87 " N = time series length = length of -data file\n"
88 " M = number of voxels in mask\n"
89 " p = number of stimulus time series to estimate\n"
90 " = number of parameters in -map file\n"
91 " q = number of baseline parameters\n"
92 " and ' = matrix transpose operator\n"
93 " Next, define matrix Z (Y detrended relative to columns of F) by\n"
94 " -1\n"
95 " Z = [I - F(F'F) F'] Y\n"
96 "-------------------------------------------------------------------\n"
97 " The method C solution is given by\n"
98 " -1\n"
99 " V0 = Z A [A'A]\n"
100 "\n"
101 " This solution minimizes the sum of squares over the N*M elements\n"
102 " of the matrix Y - V A' + F C' (N.B.: A' means A-transpose).\n"
103 "-------------------------------------------------------------------\n"
104 " The method K solution is given by\n"
105 " -1 -1\n"
106 " W = [Z Z'] Z A and then V = W [W'W]\n"
107 "\n"
108 " This solution minimizes the sum of squares of the difference between\n"
109 " the A(V) predicted from V and the input A, where A(V) is given by\n"
110 " -1\n"
111 " A(V) = Z' V [V'V] = Z'W\n"
112 "-------------------------------------------------------------------\n"
113 " Technically, the solution is unidentfiable up to an arbitrary\n"
114 " multiple of the columns of F (i.e., V = V0 + F G, where G is\n"
115 " an arbitrary q x p matrix); the solution above is the solution\n"
116 " that is orthogonal to the columns of F.\n"
117 "\n"
118 "-- RWCox - March 2006 - purely for experimental purposes!\n"
119 ) ;
120
121 printf("\n"
122 "===================== EXAMPLE USAGE =====================================\n"
123 "** Step 1: From a training dataset, generate activation map.\n"
124 " The input dataset has 4 runs, each 108 time points long. 3dDeconvolve\n"
125 " is used on the first 3 runs (time points 0..323) to generate the\n"
126 " activation map. There are two visual stimuli (Complex and Simple).\n"
127 "\n"
128 " 3dDeconvolve -x1D xout_short_two.1D -input rall_vr+orig'[0..323]' \\\n"
129 " -num_stimts 2 \\\n"
130 " -stim_file 1 hrf_complex.1D -stim_label 1 Complex \\\n"
131 " -stim_file 2 hrf_simple.1D -stim_label 2 Simple \\\n"
132 " -concat '1D:0,108,216' \\\n"
133 " -full_first -fout -tout \\\n"
134 " -bucket func_ht2_short_two -cbucket cbuc_ht2_short_two\n"
135 "\n"
136 " N.B.: You may want to de-spike, smooth, and register the 3D+time\n"
137 " dataset prior to the analysis (as usual). These steps are not\n"
138 " shown here -- I'm presuming you know how to use AFNI already.\n"
139 "\n"
140 "** Step 2: Create a mask of highly activated voxels.\n"
141 " The F statistic threshold is set to 30, corresponding to a voxel-wise\n"
142 " p = 1e-12 = very significant. The mask is also lightly clustered, and\n"
143 " restricted to brain voxels.\n"
144 "\n"
145 " 3dAutomask -prefix Amask rall_vr+orig\n"
146 " 3dcalc -a 'func_ht2_short+orig[0]' -b Amask+orig -datum byte \\\n"
147 " -nscale -expr 'step(a-30)*b' -prefix STmask300\n"
148 " 3dmerge -dxyz=1 -1clust 1.1 5 -prefix STmask300c STmask300+orig\n"
149 "\n"
150 "** Step 3: Run 3dInvFMRI to estimate the stimulus functions in run #4.\n"
151 " Run #4 is time points 324..431 of the 3D+time dataset (the -data\n"
152 " input below). The -map input is the beta weights extracted from\n"
153 " the -cbucket output of 3dDeconvolve.\n"
154 "\n"
155 " 3dInvFMRI -mask STmask300c+orig \\\n"
156 " -data rall_vr+orig'[324..431]' \\\n"
157 " -map cbuc_ht2_short_two+orig'[6..7]' \\\n"
158 " -polort 1 -alpha 0.01 -median5 -method K \\\n"
159 " -out ii300K_short_two.1D\n"
160 "\n"
161 " 3dInvFMRI -mask STmask300c+orig \\\n"
162 " -data rall_vr+orig'[324..431]' \\\n"
163 " -map cbuc_ht2_short_two+orig'[6..7]' \\\n"
164 " -polort 1 -alpha 0.01 -median5 -method C \\\n"
165 " -out ii300C_short_two.1D\n"
166 "\n"
167 "** Step 4: Plot the results, and get confused.\n"
168 "\n"
169 " 1dplot -ynames VV KK CC -xlabel Run#4 -ylabel ComplexStim \\\n"
170 " hrf_complex.1D'{324..432}' \\\n"
171 " ii300K_short_two.1D'[0]' \\\n"
172 " ii300C_short_two.1D'[0]'\n"
173 "\n"
174 " 1dplot -ynames VV KK CC -xlabel Run#4 -ylabel SimpleStim \\\n"
175 " hrf_simple.1D'{324..432}' \\\n"
176 " ii300K_short_two.1D'[1]' \\\n"
177 " ii300C_short_two.1D'[1]'\n"
178 "\n"
179 " N.B.: I've found that method K works better if MORE voxels are\n"
180 " included in the mask (lower threshold) and method C if\n"
181 " FEWER voxels are included. The above threshold gave 945\n"
182 " voxels being used to determine the 2 output time series.\n"
183 "=========================================================================\n"
184 ) ;
185
186 PRINT_COMPILE_DATE ; exit(0) ;
187 }
188
189 /**--- bureaucracy ---**/
190
191 mainENTRY("3dInvFMRI main"); machdep();
192 PRINT_VERSION("3dInvFMRI"); AUTHOR("Zhark");
193 AFNI_logger("3dInvFMRI",argc,argv) ;
194
195 /**--- scan command line ---**/
196
197 iarg = 1 ;
198 while( iarg < argc ){
199
200 if( strcmp(argv[iarg],"-method") == 0 ){
201 switch( argv[++iarg][0] ){
202 default:
203 WARNING_message("Ignoring illegal -method '%s'",argv[iarg]) ;
204 break ;
205 case 'C': method = METHOD_C ; break ;
206 case 'K': method = METHOD_K ; break ;
207 }
208 iarg++ ; continue ;
209 }
210
211 if( strcmp(argv[iarg],"-fir5") == 0 ){
212 if( nmed > 0 ) WARNING_message("Ignoring -fir5 in favor of -median5") ;
213 else nfir = 5 ;
214 iarg++ ; continue ;
215 }
216
217 if( strcmp(argv[iarg],"-median5") == 0 ){
218 if( nfir > 0 ) WARNING_message("Ignoring -median5 in favor of -fir5") ;
219 else nmed = 5 ;
220 iarg++ ; continue ;
221 }
222
223 if( strcmp(argv[iarg],"-alpha") == 0 ){
224 alpha = (float)strtod(argv[++iarg],NULL) ;
225 if( alpha <= 0.0f ){
226 alpha = 0.0f ; WARNING_message("-alpha '%s' ignored!",argv[iarg]) ;
227 }
228 iarg++ ; continue ;
229 }
230
231 if( strcmp(argv[iarg],"-data") == 0 || strcmp(argv[iarg],"-input") == 0 ){
232 if( yset != NULL ) ERROR_exit("Can't input 2 3D+time datasets") ;
233 yset = THD_open_dataset(argv[++iarg]) ;
234 CHECK_OPEN_ERROR(yset,argv[iarg]) ;
235 nt = DSET_NVALS(yset) ;
236 if( nt < 2 ) ERROR_exit("Only 1 sub-brick in dataset %s",argv[iarg]) ;
237 nxyz = DSET_NVOX(yset) ;
238 iarg++ ; continue ;
239 }
240
241 if( strcmp(argv[iarg],"-map") == 0 ){
242 if( aset != NULL ) ERROR_exit("Can't input 2 -map datasets") ;
243 aset = THD_open_dataset(argv[++iarg]) ;
244 CHECK_OPEN_ERROR(aset,argv[iarg]) ;
245 nparam = DSET_NVALS(aset) ;
246 iarg++ ; continue ;
247 }
248
249 if( strcmp(argv[iarg],"-mapwt") == 0 ){
250 if( wset != NULL ) ERROR_exit("Can't input 2 -mapwt datasets") ;
251 wset = THD_open_dataset(argv[++iarg]) ;
252 CHECK_OPEN_ERROR(wset,argv[iarg]) ;
253 iarg++ ; continue ;
254 }
255
256 if( strcmp(argv[iarg],"-mask") == 0 ){
257 if( mset != NULL ) ERROR_exit("Can't input 2 -mask datasets") ;
258 mset = THD_open_dataset(argv[++iarg]) ;
259 CHECK_OPEN_ERROR(mset,argv[iarg]) ;
260 iarg++ ; continue ;
261 }
262
263 if( strcmp(argv[iarg],"-polort") == 0 ){
264 char *cpt ;
265 polort = (int)strtod(argv[++iarg],&cpt) ;
266 if( *cpt != '\0' ) WARNING_message("Illegal non-numeric value after -polort") ;
267 iarg++ ; continue ;
268 }
269
270 if( strcmp(argv[iarg],"-out") == 0 ){
271 fname_out = strdup(argv[++iarg]) ;
272 if( !THD_filename_ok(fname_out) )
273 ERROR_exit("Bad -out filename '%s'",fname_out) ;
274 iarg++ ; continue ;
275 }
276
277 if( strcmp(argv[iarg],"-base") == 0 ){
278 if( fimar == NULL ) INIT_IMARR(fimar) ;
279 qim = mri_read_1D( argv[++iarg] ) ;
280 if( qim == NULL ) ERROR_exit("Can't read 1D file %s",argv[iarg]) ;
281 ADDTO_IMARR(fimar,qim) ;
282 iarg++ ; continue ;
283 }
284
285 ERROR_exit("Unrecognized option '%s'",argv[iarg]) ;
286 }
287
288 /**--- finish up processing options ---**/
289
290 if( yset == NULL ) ERROR_exit("No input 3D+time dataset?!") ;
291 if( aset == NULL ) ERROR_exit("No input FMRI -map dataset?!") ;
292
293 if( DSET_NVOX(aset) != nxyz )
294 ERROR_exit("Grid mismatch between -data and -map") ;
295
296 INFO_message("Loading dataset for Y") ;
297 DSET_load(yset); CHECK_LOAD_ERROR(yset) ;
298 INFO_message("Loading dataset for A") ;
299 DSET_load(aset); CHECK_LOAD_ERROR(aset) ;
300
301 if( wset != NULL ){
302 if( DSET_NVOX(wset) != nxyz )
303 ERROR_exit("Grid mismatch between -data and -mapwt") ;
304 nwt = DSET_NVALS(wset) ;
305 if( nwt > 1 && nwt != nparam )
306 ERROR_exit("Wrong number of values=%d in -mapwt; should be 1 or %d",
307 nwt , nparam ) ;
308 INFO_message("Loading dataset for mapwt") ;
309 DSET_load(wset); CHECK_LOAD_ERROR(wset) ;
310 }
311
312 if( mset != NULL ){
313 if( DSET_NVOX(mset) != nxyz )
314 ERROR_exit("Grid mismatch between -data and -mask") ;
315 INFO_message("Loading dataset for mask") ;
316 DSET_load(mset); CHECK_LOAD_ERROR(mset) ;
317 mask = THD_makemask( mset , 0 , 1.0f,-1.0f ); DSET_delete(mset);
318 nmask = THD_countmask( nxyz , mask ) ;
319 if( nmask < 3 ){
320 WARNING_message("Mask has %d voxels -- ignoring!",nmask) ;
321 free(mask) ; mask = NULL ; nmask = 0 ;
322 }
323 }
324
325 nvox = (nmask > 0) ? nmask : nxyz ;
326 INFO_message("N = time series length = %d",nt ) ;
327 INFO_message("M = number of voxels = %d",nvox ) ;
328 INFO_message("p = number of params = %d",nparam) ;
329
330 /**--- set up baseline funcs in one array ---*/
331
332 nqbase = (polort >= 0 ) ? polort+1 : 0 ;
333 if( fimar != NULL ){
334 for( kk=0 ; kk < IMARR_COUNT(fimar) ; kk++ ){
335 qim = IMARR_SUBIMAGE(fimar,kk) ;
336 if( qim != NULL && qim->nx != nt )
337 WARNING_message("-base #%d length=%d; data length=%d",kk+1,qim->nx,nt) ;
338 nqbase += qim->ny ;
339 }
340 }
341
342 INFO_message("q = number of baselines = %d",nqbase) ;
343
344 #undef F
345 #define F(i,j) flar[(i)+(j)*nt] /* nt X nqbase */
346 if( nqbase > 0 ){
347 fim = mri_new( nt , nqbase , MRI_float ) ; /* F matrix */
348 flar = MRI_FLOAT_PTR(fim) ;
349 bb = 0 ;
350 if( polort >= 0 ){ /** load polynomial baseline **/
351 double a = 2.0/(nt-1.0) ;
352 for( jj=0 ; jj <= polort ; jj++ ){
353 for( ii=0 ; ii < nt ; ii++ )
354 F(ii,jj) = (float)Plegendre( a*ii-1.0 , jj ) ;
355 }
356 bb = polort+1 ;
357 }
358 #undef Q
359 #define Q(i,j) qar[(i)+(j)*qim->nx] /* qim->nx X qim->ny */
360
361 if( fimar != NULL ){ /** load -base baseline columns **/
362 for( kk=0 ; kk < IMARR_COUNT(fimar) ; kk++ ){
363 qim = IMARR_SUBIMAGE(fimar,kk) ; qar = MRI_FLOAT_PTR(qim) ;
364 for( jj=0 ; jj < qim->ny ; jj++ ){
365 for( ii=0 ; ii < nt ; ii++ )
366 F(ii,bb+jj) = (ii < qim->nx) ? Q(ii,jj) : 0.0f ;
367 }
368 bb += qim->ny ;
369 }
370 DESTROY_IMARR(fimar) ; fimar=NULL ;
371 }
372
373 /* remove mean from each column after first? */
374
375 if( polort >= 0 && nqbase > 1 ){
376 float sum ;
377 for( jj=1 ; jj < nqbase ; jj++ ){
378 sum = 0.0f ;
379 for( ii=0 ; ii < nt ; ii++ ) sum += F(ii,jj) ;
380 sum /= nt ;
381 for( ii=0 ; ii < nt ; ii++ ) F(ii,jj) -= sum ;
382 }
383 }
384
385 /* compute pseudo-inverse of baseline matrix,
386 so we can project it out from the data time series */
387
388 /* -1 */
389 /* (F'F) F' matrix */
390
391 INFO_message("Computing pseudo-inverse of baseline matrix F") ;
392 pfim = mri_matrix_psinv(fim,NULL,0.0f) ; par = MRI_FLOAT_PTR(pfim) ;
393
394 #undef P
395 #define P(i,j) par[(i)+(j)*nqbase] /* nqbase X nt */
396
397 #if 0
398 qim = mri_matrix_transpose(pfim) ; /** save to disk? **/
399 mri_write_1D( "Fpsinv.1D" , qim ) ;
400 mri_free(qim) ;
401 #endif
402 }
403
404 /**--- set up map image into aim/aar = A matrix ---**/
405
406 #undef GOOD
407 #define GOOD(i) (mask==NULL || mask[i])
408
409 #undef A
410 #define A(i,j) aar[(i)+(j)*nvox] /* nvox X nparam */
411
412 INFO_message("Loading map matrix A") ;
413 aim = mri_new( nvox , nparam , MRI_float ); aar = MRI_FLOAT_PTR(aim);
414 for( jj=0 ; jj < nparam ; jj++ ){
415 for( ii=kk=0 ; ii < nxyz ; ii++ ){
416 if( GOOD(ii) ){ A(kk,jj) = THD_get_voxel(aset,ii,jj); kk++; }
417 }}
418 DSET_unload(aset) ;
419
420 /**--- set up map weight into wim/war ---**/
421
422 #undef WT
423 #define WT(i,j) war[(i)+(j)*nvox] /* nvox X nparam */
424
425 if( wset != NULL ){
426 int numneg=0 , numpos=0 ;
427 float fac ;
428
429 INFO_message("Loading map weight matrix") ;
430 wim = mri_new( nvox , nwt , MRI_float ) ; war = MRI_FLOAT_PTR(wim) ;
431 for( jj=0 ; jj < nwt ; jj++ ){
432 for( ii=kk=0 ; ii < nxyz ; ii++ ){
433 if( GOOD(ii) ){
434 WT(kk,jj) = THD_get_voxel(wset,ii,jj);
435 if( WT(kk,jj) > 0.0f ){ numpos++; WT(kk,jj) = sqrt(WT(kk,jj)); }
436 else if( WT(kk,jj) < 0.0f ){ numneg++; WT(kk,jj) = 0.0f; }
437 kk++;
438 }
439 }}
440 DSET_unload(wset) ;
441 if( numpos <= nparam )
442 WARNING_message("Only %d positive weights found in -wtmap!",numpos) ;
443 if( numneg > 0 )
444 WARNING_message("%d negative weights found in -wtmap!",numneg) ;
445
446 for( jj=0 ; jj < nwt ; jj++ ){
447 fac = 0.0f ;
448 for( kk=0 ; kk < nvox ; kk++ ) if( WT(kk,jj) > fac ) fac = WT(kk,jj) ;
449 if( fac > 0.0f ){
450 fac = 1.0f / fac ;
451 for( kk=0 ; kk < nvox ; kk++ ) WT(kk,jj) *= fac ;
452 }
453 }
454 }
455
456 /**--- set up data image into yim/yar = Y matrix ---**/
457
458 #undef Y
459 #define Y(i,j) yar[(i)+(j)*nt] /* nt X nvox */
460
461 INFO_message("Loading data matrix Y") ;
462 yim = mri_new( nt , nvox , MRI_float ); yar = MRI_FLOAT_PTR(yim);
463 for( ii=0 ; ii < nt ; ii++ ){
464 for( jj=kk=0 ; jj < nxyz ; jj++ ){
465 if( GOOD(jj) ){ Y(ii,kk) = THD_get_voxel(yset,jj,ii); kk++; }
466 }}
467 DSET_unload(yset) ;
468
469 /**--- project baseline out of data image = Z matrix ---**/
470
471 if( pfim != NULL ){
472 #undef T
473 #define T(i,j) tar[(i)+(j)*nt] /* nt X nvox */
474 INFO_message("Projecting baseline out of Y") ;
475 qim = mri_matrix_mult( pfim , yim ) ; /* nqbase X nvox */
476 tim = mri_matrix_mult( fim , qim ) ; /* nt X nvox */
477 tar = MRI_FLOAT_PTR(tim) ; /* Y projected onto baseline */
478 for( jj=0 ; jj < nvox ; jj++ )
479 for( ii=0 ; ii < nt ; ii++ ) Y(ii,jj) -= T(ii,jj) ;
480 mri_free(tim); mri_free(qim); mri_free(pfim); mri_free(fim);
481 }
482
483 /***** At this point:
484 matrix A is in aim,
485 matrix Z is in yim.
486 Solve for V into vim, using the chosen method *****/
487
488 switch( method ){
489 default: ERROR_exit("Illegal method code! WTF?") ; /* Huh? */
490
491 /*.....................................................................*/
492 case METHOD_C:
493 /**--- compute pseudo-inverse of A map ---**/
494
495 INFO_message("Method C: Computing pseudo-inverse of A") ;
496 if( wim != NULL ) WARNING_message("Ignoring -mapwt dataset") ;
497 pfim = mri_matrix_psinv(aim,NULL,alpha) ; /* nparam X nvox */
498 if( pfim == NULL ) ERROR_exit("mri_matrix_psinv() fails") ;
499 mri_free(aim) ;
500
501 /**--- and apply to data to get results ---*/
502
503 INFO_message("Computing result V") ;
504 vim = mri_matrix_multranB( yim , pfim ) ; /* nt x nparam */
505 mri_free(pfim) ; mri_free(yim) ;
506 break ;
507
508 /*.....................................................................*/
509 case METHOD_K:
510 /**--- compute pseudo-inverse of transposed Z ---*/
511
512 INFO_message("Method K: Computing pseudo-inverse of Z'") ;
513 if( nwt > 1 ){
514 WARNING_message("Ignoring -mapwt dataset: more than 1 sub-brick") ;
515 nwt = 0 ; mri_free(wim) ; wim = NULL ; war = NULL ;
516 }
517
518 if( nwt == 1 ){
519 float fac ;
520 for( kk=0 ; kk < nvox ; kk++ ){
521 fac = war[kk] ;
522 for( ii=0 ; ii < nt ; ii++ ) Y(ii,kk) *= fac ;
523 for( ii=0 ; ii < nparam ; ii++ ) A(kk,ii) *= fac ;
524 }
525 }
526
527 tim = mri_matrix_transpose(yim) ; mri_free(yim) ;
528 pfim = mri_matrix_psinv(tim,NULL,alpha) ; mri_free(tim) ;
529 if( pfim == NULL ) ERROR_exit("mri_matrix_psinv() fails") ;
530
531 INFO_message("Computing W") ;
532 tim = mri_matrix_mult( pfim , aim ) ;
533 mri_free(aim) ; mri_free(pfim) ;
534
535 INFO_message("Computing result V") ;
536 pfim = mri_matrix_psinv(tim,NULL,0.0f) ; mri_free(tim) ;
537 vim = mri_matrix_transpose(pfim) ; mri_free(pfim);
538 break ;
539
540 } /* end of switch on method */
541
542 if( wim != NULL ) mri_free(wim) ;
543
544 /**--- smooth? ---**/
545
546 if( nfir > 0 && vim->nx > nfir ){
547 INFO_message("FIR-5-ing result") ;
548 var = MRI_FLOAT_PTR(vim) ;
549 for( jj=0 ; jj < vim->ny ; jj++ )
550 linear_filter_reflect( nfir,firwt , vim->nx , var + (jj*vim->nx) ) ;
551 }
552
553 if( nmed > 0 && vim->nx > nmed ){
554 INFO_message("Median-5-ing result") ;
555 var = MRI_FLOAT_PTR(vim) ;
556 for( jj=0 ; jj < vim->ny ; jj++ )
557 median5_filter_reflect( vim->nx , var + (jj*vim->nx) ) ;
558 }
559
560 /**--- write results ---**/
561
562 INFO_message("Writing result to '%s'",fname_out) ;
563 mri_write_1D( fname_out , vim ) ;
564 exit(0) ;
565 }
566
567 /*-------------------------------------------------------------------------*/
568
linear_filter_reflect(int ntap,float * wt,int npt,float * x)569 static void linear_filter_reflect( int ntap, float *wt, int npt, float *x )
570 {
571 int ii , nt2=(ntap-1)/2 , jj , np1=npt-1 , np2=2*(npt-1) ;
572 register float sum ;
573 static int nfar=0 ;
574 static float *far=NULL ;
575
576 if( ntap >= npt || wt == NULL || x == NULL ) return ;
577
578 if( npt > nfar ){
579 if( far != NULL ) free(far) ;
580 far = (float *)malloc(sizeof(float)*npt) ; nfar = npt ;
581 }
582
583 #undef XX
584 #define XX(i) ( ((i)<0) ? far[-(i)] : ((i)>np1) ? far[np2-(i)] : far[(i)] )
585
586 memcpy( far , x , sizeof(float)*npt ) ;
587
588 for( ii=0 ; ii < npt ; ii++ ){
589 for( sum=0.0f,jj=0 ; jj < ntap ; jj++ ) sum += wt[jj] * XX(ii-nt2+jj) ;
590 x[ii] = sum ;
591 }
592
593 return ;
594 }
595
596 /*-------------------------------------------------------------------------*/
597
598 #undef SWAP
599 #define SWAP(x,y) (temp=x,x=y,y=temp)
600 #undef SORT2
601 #define SORT2(a,b) if(a>b) SWAP(a,b);
602
median_float5(float * p)603 static float median_float5(float *p)
604 {
605 register float temp ;
606 SORT2(p[0],p[1]) ; SORT2(p[3],p[4]) ; SORT2(p[0],p[3]) ;
607 SORT2(p[1],p[4]) ; SORT2(p[1],p[2]) ; SORT2(p[2],p[3]) ;
608 SORT2(p[1],p[2]) ; return(p[2]) ;
609 }
610
median5_filter_reflect(int npt,float * x)611 static void median5_filter_reflect( int npt, float *x )
612 {
613 int ii ;
614 float p[5] ;
615 static int nfar=0 ;
616 static float *far=NULL ;
617
618 if( x == NULL || npt < 5 ) return ;
619
620 if( npt > nfar ){
621 if( far != NULL ) free(far) ;
622 far = (float *)malloc(sizeof(float)*(npt+4)) ; nfar = npt ;
623 }
624
625 memcpy( far+2 , x , sizeof(float)*npt ) ;
626 far[0] = x[2]; far[1] = x[1]; far[npt+2] = x[npt-2]; far[npt+3] = x[npt-3];
627
628 for( ii=0 ; ii < npt ; ii++ ){
629 p[0] = far[ii]; p[1] = far[ii+1]; p[2] = far[ii+2];
630 p[3] = far[ii+3]; p[4] = far[ii+4];
631 x[ii] = median_float5(p) ;
632 }
633
634 return ;
635 }
636