1 #include "mrilib.h"
2
3 /** To do:
4 allow output for a range of temporal subsets
5 allow censoring
6 **/
7
8 #ifdef USE_OMP
9 #include <omp.h>
10 #endif
11
12 #include "parser.h"
13
14 /*----------------------------------------------------------------------------*/
etime(struct timeval * t,int Report)15 static float etime (struct timeval *t, int Report )
16 {/*etime*/
17 struct timeval tn;
18 float Time_Fact = 1000000.0;
19 float delta_t;
20
21 /* get time */
22 gettimeofday(&tn, NULL);
23
24 if (Report)
25 {
26 delta_t = ( ( (float)(tn.tv_sec - t->tv_sec)*Time_Fact) +
27 (float)(tn.tv_usec - t->tv_usec) ) /Time_Fact;
28 }
29 else
30 {
31 t->tv_sec = tn.tv_sec;
32 t->tv_usec = tn.tv_usec;
33 delta_t = 0.0;
34 }
35
36 return (delta_t);
37 }
38
39 /*----------------------------------------------------------------------------*/
40
vstep_print(void)41 static void vstep_print(void)
42 {
43 static char xx[10] = "0123456789" ; static int vn=0 ;
44 fprintf(stderr , "%c" , xx[vn%10] ) ;
45 if( vn%10 == 9) fprintf(stderr,".") ;
46 vn++ ;
47 }
48
49 /*----------------------------------------------------------------------------*/
50 /* Binary search for tt in a sorted integer array. */
51
mybsearch_int(int tt,int nar,int * ar)52 static int mybsearch_int( int tt , int nar , int *ar )
53 {
54 int targ , ii , jj , kk , nn ;
55
56 if( nar <= 0 || ar == NULL ) return -1 ; /* bad inputs */
57
58 targ = tt ; ii = 0 ; jj = nar-1 ; /* setup */
59
60 if( targ < ar[0] ) return -1 ; /* not found */
61 else if( targ == ar[0] ) return 0 ; /* at start! */
62
63 if( targ > ar[jj] ) return -1 ; /* not found */
64 else if( targ == ar[jj] ) return jj ; /* at end! */
65
66 /* at the start of this loop, we've already checked
67 indexes ii and jj, so check the middle of them (kk),
68 and if that doesn't work, make the middle the
69 new ii or the new jj -- so again we will have
70 checked both ii and jj when the loop iterates back. */
71
72 while( jj-ii > 1 ){
73 kk = (ii+jj) / 2 ; /* midpoint */
74 nn = ar[kk] - targ ; /* sign of difference */
75 if( nn == 0 ) return kk ; /* found it! */
76 if( nn < 0 ) ii = kk ; /* must be above kk */
77 else jj = kk ; /* must be below kk */
78 }
79
80 return -1 ;
81 }
82
usage_3dTcorrMap(int detail)83 void usage_3dTcorrMap(int detail) {
84 printf(
85 "Usage: 3dTcorrMap [options]\n"
86 "For each voxel time series, computes the correlation between it\n"
87 "and all other voxels, and combines this set of values into the\n"
88 "output dataset(s) in some way.\n"
89 "\n"
90 "Supposed to give a measure of how 'connected' each voxel is\n"
91 "to the rest of the brain. [[As if life were that simple.]]\n"
92 "\n"
93 "---------\n"
94 "WARNINGS:\n"
95 "---------\n"
96 "** This program takes a LONG time to run.\n"
97 "** This program will use a LOT of memory.\n"
98 "** Don't say I didn't warn you about these facts, and don't whine.\n"
99 "\n"
100 "--------------\n"
101 "Input Options:\n"
102 "--------------\n"
103 " -input dd = Read 3D+time dataset 'dd' (a mandatory option).\n"
104 " This provides the time series to be correlated\n"
105 " en masse.\n"
106 " ** This is a non-optional 'option': you MUST supply\n"
107 " and input dataset!\n"
108 "\n"
109 " -seed bb = Read 3D+time dataset 'bb'.\n"
110 " ** If you use this option, for each voxel in the\n"
111 " -seed dataset, its time series is correlated\n"
112 " with every voxel in the -input dataset, and\n"
113 " then that collection of correlations is processed\n"
114 " to produce the output for that voxel.\n"
115 " ** If you don't use -seed, then the -input dataset\n"
116 " is the -seed dataset [i.e., the normal usage].\n"
117 " ** The -seed and -input datasets must have the\n"
118 " same number of time points and the same number\n"
119 " of voxels!\n"
120 " ** Unlike the -input dataset, the -seed dataset is not\n"
121 " preprocessed (i.e., no detrending/bandpass or blur).\n"
122 " (The main purpose of this -seed option is to)\n"
123 " (allow you to preprocess the seed voxel time)\n"
124 " (series in some personalized and unique way.)\n"
125 "\n"
126 " -mask mmm = Read dataset 'mmm' as a voxel mask.\n"
127 " -automask = Create a mask from the input dataset.\n"
128 " ** -mask and -automask are mutually exclusive!\n"
129 " ** If you don't use one of these masking options, then\n"
130 " all voxels will be processed, and the program will\n"
131 " probably run for a VERY long time.\n"
132 " ** Voxels with constant time series will be automatically\n"
133 " excluded.\n"
134 "\n"
135 "----------------------------------\n"
136 "Time Series Preprocessing Options: (applied only to -input, not to -seed)\n"
137 "[[[[ In general, it would be better to pre-process with afni_proc.py ]]]]\n"
138 "----------------------------------\n"
139 "TEMPORAL FILTERING:\n"
140 "-------------------\n"
141 " -polort m = Remove polynomial trend of order 'm', for m=-1..19.\n"
142 " [default is m=1; removal is by least squares].\n"
143 " ** Using m=-1 means no detrending; this is only useful\n"
144 " for data/information that has been pre-processed\n"
145 " (e.g., using the 3dBandpass program).\n"
146 "\n"
147 " -bpass L H = Bandpass the data between frequencies L and H (in Hz).\n"
148 " ** If the input dataset does not have a time step defined,\n"
149 " then TR = 1 s will be assumed for this purpose.\n"
150 " **** -bpass and -polort are mutually exclusive!\n"
151 "\n"
152 " -ort ref = 1D file with other time series to be removed from -input\n"
153 " (via least squares regression) before correlation.\n"
154 " ** Each column in 'ref' will be regressed out of\n"
155 " each -input voxel time series.\n"
156 " ** -ort can be used with -polort and/or -bandpass.\n"
157 " ** You can use programs like 3dmaskave and 3dmaskSVD\n"
158 " to create reference files from regions of the\n"
159 " input dataset (e.g., white matter, CSF).\n"
160 #if 0
161 "\n"
162 " -PCort n mmm = From the -input dataset, extract the time series\n"
163 " from the mask 'mmm', then detrend them (as in\n"
164 " -polort/-bpass + -ort), then compute the first\n"
165 " 'n' principal components (SVD eigenvectors)\n"
166 " and regress these out of the -input data as well.\n"
167 " ** 'mmm' might be a white matter mask, for example.\n"
168 " ** Similar to using program '3dmaskSVD -sval n' and\n"
169 " using that output with '-ort' in this program.\n"
170 " **** This is the last preprocessing step before the\n"
171 " long long all-voxel-pairs correlation begins.\n"
172 #endif
173 "\n"
174 "SPATIAL FILTERING: (only for volumetric input datasets) \n"
175 "-----------------\n"
176 " -Gblur ff = Gaussian blur the -input dataset (inside the mask)\n"
177 " using a kernel width of 'ff' mm.\n"
178 " ** Uses the same approach as program 3dBlurInMask.\n"
179 "\n"
180 " -Mseed rr = When extracting the seed voxel time series from the\n"
181 " (preprocessed) -input dataset, average it over a radius\n"
182 " of 'rr' mm prior to doing the correlations with all\n"
183 " the voxel time series from the -input dataset.\n"
184 " ** This extra smoothing is said by some mystics to\n"
185 " improve and enhance the results. YMMV.\n"
186 " ** Only voxels inside the mask will be used.\n"
187 " ** A negative value for 'rr' means to treat the voxel\n"
188 " dimensions as all equal to 1.0 mm; thus, '-Mseed -1.0'\n"
189 " means to average a voxel with its 6 nearest\n"
190 " neighbors in the -input dataset 3D grid.\n"
191 " ** -Mseed and -seed are mutually exclusive!\n"
192 " (It makes NO sense to use both options.)\n"
193 #if 0
194 "\n"
195 " -Pseed rr = Similar to -Mseed, but instead of averaging, use the\n"
196 " first principal component (first SVD eigenvector),\n"
197 " instead of the mean vector.\n"
198 " **** -Mseed and -Pseed are mutually exclusive!\n"
199 #endif
200 "\n"
201 "---------------\n"
202 "Output Options: (at least one of these must be given!)\n"
203 "---------------\n"
204 " -Mean pp = Save average correlations into dataset prefix 'pp'\n"
205 " ** As pointed out to me by CC, '-Mean' is the same\n"
206 " as computing the correlation map with the 1D file\n"
207 " that is the mean of all the normalized time series\n"
208 " in the mask -- that is, a form of the global signal.\n"
209 " Such a calculation could be done much faster with\n"
210 " program 3dTcorr1D.\n"
211 " ** Nonlinear combinations of the correlations, as done by\n"
212 " the options below, can't be done in such a simple way.\n"
213 " -Zmean pp = Save tanh of mean arctanh(correlation) into 'pp'\n"
214 " -Qmean pp = Save RMS(correlation) into 'pp'\n"
215 " -Pmean pp = Save average of squared positive correlations into 'pp'\n"
216 " (negative correlations don't count in this calculation)\n"
217 " -Thresh tt pp\n"
218 " = Save the COUNT of how many voxels survived thresholding\n"
219 " at level abs(correlation) >= tt (for some tt > 0).\n"
220 "\n"
221 " -VarThresh t0 t1 dt pp\n"
222 " = Save the COUNT of how many voxels survive thresholding\n"
223 " at several levels abs(correlation) >= tt, for\n"
224 " tt = t0, t0+dt, ..., t1. This option produces\n"
225 " a multi-volume dataset, with prefix 'pp'.\n"
226 " -VarThreshN t0 t1 dt pp\n"
227 " = Like '-VarThresh', but the output counts are\n"
228 " 'Normalized' (divided) by the expected number\n"
229 " of such supra-threshold voxels that would occur\n"
230 " from white noise timeseries.\n"
231 " ** N.B.: You can't use '-VarThresh' and '-VarThreshN'\n"
232 " in the same run of the program!\n"
233 " -CorrMap pp\n"
234 " Output at each voxel the entire correlation map, into\n"
235 " a dataset with prefix 'pp'.\n"
236 " ** Essentially this does what 3dAutoTcorrelate would,\n"
237 " with some of the additional options offered here.\n"
238 " ** N.B.: Output dataset will be HUGE and BIG in most cases.\n"
239 " -CorrMask\n"
240 " By default, -CorrMap outputs a sub-brick for EACH\n"
241 " input dataset voxel, even those that are NOT in\n"
242 " the mask (such sub-bricks will be all zero).\n"
243 " If you want to eliminate these sub-bricks, use\n"
244 " this option.\n"
245 " ** N.B.: The label for the sub-brick that was seeded\n"
246 " from voxel (i,j,k) will be of the form\n"
247 " v032.021.003 (when i=32, j=21, k=3).\n"
248 "\n"
249 " --** The following 3 options let you create a customized **--\n"
250 " --** method of combining the correlations, if the above **--\n"
251 " --** techniques do not meet your needs. (Of course, you **--\n"
252 " --** could also use '-CorrMap' and then process the big **--\n"
253 " --** output dataset yourself later, in some clever way.) **--\n"
254 "\n"
255 " -Aexpr expr ppp\n"
256 " = For each correlation 'r', compute the calc-style\n"
257 " expression 'expr', and average these values to get\n"
258 " the output that goes into dataset 'ppp'.\n"
259 " -Cexpr expr ppp\n"
260 " = As in '-Aexpr', but only average together nonzero\n"
261 " values computed by 'expr'. Example:\n"
262 " -Cexpr 'step(r-0.3)*r' TCa03\n"
263 " would compute (for each voxel) the average of all\n"
264 " correlation coefficients larger than 0.3.\n"
265 " -Sexpr expr ppp\n"
266 " = As above, but the sum of the expressions is computed\n"
267 " rather than the average. Example:\n"
268 " -Sexpr 'step(r-0.3)' TCn03\n"
269 " would compute the number of voxels with correlation\n"
270 " coefficients larger than 0.3.\n"
271 " ** N.B.: At most one '-?expr' option can be used in\n"
272 " the same run of the program!\n"
273 " ** N.B.: Only the symbols 'r' and 'z' [=atanh(r)] have any\n"
274 " meaning in the expression; all other symbols will\n"
275 " be treated as zeroes.\n"
276 "\n"
277 " -Hist N ppp\n"
278 " = For each voxel, save a histogram of the correlation\n"
279 " coefficients into dataset ppp.\n"
280 " ** N values will be saved per voxel, with the i'th\n"
281 " sub-brick containing the count for the range\n"
282 " -1+i*D <= r < -1+(i+1)*D with D=2/N and i=0..N-1\n"
283 " ** N must be at least 20, and at most 1000.\n"
284 " * N=200 is good; then D=0.01, yielding a decent resolution.\n"
285 " ** The output dataset is short format; thus, the maximum\n"
286 " count in any bin will be 32767.\n"
287 " ** The output from this option will probably require further\n"
288 " processing before it can be useful -- but it is fun to\n"
289 " surf through these histograms in AFNI's graph viewer.\n"
290 "\n"
291 "----------------\n"
292 "Random Thoughts:\n"
293 "----------------\n"
294 "-- In all output calculations, the correlation of a voxel with itself\n"
295 " is ignored. If you don't understand why, step away from the keyboard.\n"
296 "-- This purely experimental program is somewhat time consuming.\n"
297 " (Of course, it's doing a LOT of calculations.)\n"
298 "-- For Kyle, AKA the new Pat (assuming such a thing were possible).\n"
299 "-- For Steve, AKA the new Kyle (which makes him the newest Pat).\n"
300 "-- RWCox - August 2008 et cetera.\n"
301 ) ;
302
303 PRINT_AFNI_OMP_USAGE("3dTcorrMap",NULL) ;
304 #ifndef USE_OMP
305 printf(" * You REALLY want to use an OpenMP version of this program,\n"
306 " and run on a fast multi-CPU computer system! Otherwise,\n"
307 " you will wait a LONG LONG time for the results.\n"
308 ) ;
309 #endif
310 PRINT_COMPILE_DATE ;
311 return;
312 }
313
314 /*----------------------------------------------------------------------------*/
315
main(int argc,char * argv[])316 int main( int argc , char *argv[] )
317 {
318 THD_3dim_dataset *xset=NULL ; int nx,ny,nz,nxy,nxyz ;
319 THD_3dim_dataset *sset=NULL ;
320 int nopt=1 , do_automask=0 ;
321 int nvox , nvals , ii,jj,kk , polort=1 , ntime ;
322 float *xsar ; float acc,cc,csum,Mcsum,Zcsum,Qcsum , Tcount ;
323 byte *mask=NULL ; int nxmask=0,nymask=0,nzmask=0 , nmask=0 , vstep=0 ;
324 int nref=0 , iv=0, N_iv=0, Tbone=0;
325 float **ref=NULL, t0=0.0, t1=0.0, ti=0.0;
326 MRI_IMAGE *ortim=NULL ; float *ortar=NULL ;
327 int *indx=NULL ; MRI_IMARR *timar=NULL ; MRI_IMARR *simar=NULL ;
328 char *Mprefix=NULL ; THD_3dim_dataset *Mset=NULL ; float *Mar=NULL ;
329 char *Zprefix=NULL ; THD_3dim_dataset *Zset=NULL ; float *Zar=NULL ;
330 char *Qprefix=NULL ; THD_3dim_dataset *Qset=NULL ; float *Qar=NULL ;
331 char *Tprefix=NULL ; THD_3dim_dataset *Tset=NULL ; float *Tar=NULL ;
332 char *Pprefix=NULL ; THD_3dim_dataset *Pset=NULL ; float *Par=NULL ;
333 char *COprefix=NULL; THD_3dim_dataset *COset=NULL ; short *COar=NULL ;
334 int COmask=0 ;
335 float Thresh=0.0f ;
336 char *Tvprefix=NULL ; THD_3dim_dataset *Tvset=NULL ;
337 float **Tvar=NULL ; float *Threshv=NULL, *Tvcount=NULL ;
338 char stmp[256];
339 int nout=0 ;
340 int isodd ; /* 29 Apr 2009: for unrolling innermost dot product */
341 struct timeval tt;
342 float dtt=0.0;
343 float *ccar=NULL ; /* 29 Apr 2009: for OpenMP usage */
344
345 float Pcsum ; int nPcsum ; /* 23 Jun 2009 */
346
347 char *expr_string=NULL , expr_type='\0' ; /* 23 Jun 2009 */
348 PARSER_code *expr_code=NULL ;
349 char *Eprefix=NULL ; THD_3dim_dataset *Eset=NULL ; float *Ear=NULL ;
350 double *atoz[26] ;
351 double *eear=NULL ; float Esum ; int nEsum ;
352 int expr_has_z=0 , expr_has_r=0 ;
353
354 char *HHprefix=NULL ; THD_3dim_dataset *HHset=NULL ;
355 int HHnum=0 ; short *HHist=NULL ; float HHdel=0.0f , HHdin=0.0f ;
356
357 float bpass_L=0.0f , bpass_H=0.0f , dtime ; int do_bpass=0 ;
358 double ctime ;
359 float gblur=0.0f ;
360 float Mseedr=0.0f , dx,dy,dz ; MCW_cluster *Mseed_nbhd=NULL ;
361 float *Mseedar=NULL ;
362
363 int PCortn=0 , PCnmask=0 ; byte *PCmask=NULL ; int PCnx=0, PCny=0, PCnz=0 ;
364 MRI_IMARR *PCimar=NULL ;
365
366 int need_acc = 0 ;
367
368 /*----*/
369
370 #if defined(USING_MCW_MALLOC) && !defined(USE_OMP)
371 enable_mcw_malloc() ;
372 #endif
373 if( argc > 2 ){
374 mainENTRY("3dTcorrMap main"); machdep(); PRINT_VERSION("3dTcorrMap");
375 AFNI_logger("3dTcorrMap",argc,argv);
376 THD_check_AFNI_version("3dTcorrMap") ;
377 AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */
378 }
379
380 /*-- option processing --*/
381
382 while( nopt < argc && argv[nopt][0] == '-' ){
383 if( strcasecmp(argv[nopt],"-h") == 0 ||
384 strcasecmp(argv[nopt],"-help") == 0 ) {
385 usage_3dTcorrMap(strlen(argv[1]) > 3 ? 2:1);
386 exit(0);
387 }
388 if( strcasecmp(argv[nopt],"-Mseed") == 0 ){
389 if( sset != NULL ) ERROR_exit("Can't use -Mseed with -seed!") ;
390 Mseedr = (float)strtod( argv[++nopt] , NULL ) ;
391 if( Mseedr == 0.0f ) ERROR_exit("Illegal value after -Mseed") ;
392 nopt++ ; continue ;
393 }
394
395 if( strcasecmp(argv[nopt],"-Gblur") == 0 ){
396 gblur = (float)strtod( argv[++nopt] , NULL ) ;
397 if( gblur < 0.0f ) ERROR_exit("Illegal value after -Gblur") ;
398 nopt++ ; continue ;
399 }
400
401 if( strcasecmp(argv[nopt],"-ort") == 0 ){
402 if( ortim != NULL ) ERROR_exit("Can't have 2 -ort options!") ;
403 ortim = mri_read_1D( argv[++nopt] ) ;
404 if( ortim == NULL ) ERROR_exit("Can't read file after -ort") ;
405 nopt++ ; continue ;
406 }
407
408 #if 0
409 if( strcasecmp(argv[nopt],"-PCort") == 0 ){
410 THD_3dim_dataset *mset ;
411 if( nopt+2 >= argc ) ERROR_exit("Need 2 arguments after '-PCort'") ;
412 if( PCortn > 0 ) ERROR_exit("Can't have 2 -PCort options") ;
413 PCortn = (int)strtod(argv[++nopt],NULL) ;
414 if( PCortn <= 0 ) ERROR_exit("Illegal 'n' after -PCort") ;
415
416 mset = THD_open_dataset( argv[++nopt] ); CHECK_OPEN_ERROR(mset,argv[nopt-1]);
417 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
418 PCnx = DSET_NX(mset); PCny = DSET_NY(mset); PCnz = DSET_NZ(mset);
419 PCmask = THD_makemask( mset , 0 , 0.5f, 0.0f ) ; DSET_delete(mset) ;
420 if( PCmask == NULL )
421 ERROR_exit("Can't make -PCort mask from dataset '%s'",argv[nopt]) ;
422 PCnmask = THD_countmask( PCnx*PCny*PCnz , PCmask ) ;
423 INFO_message("Number of voxels in -PCort mask = %d",PCnmask) ;
424 if( nmask < 9 ) ERROR_exit("-PCort mask is too small to process") ;
425 nopt++ ; continue ;
426 }
427 #endif
428
429 if( strcasecmp(argv[nopt],"-automask") == 0 ){
430 if( mask != NULL ) ERROR_exit("Can't use -automask and -mask!") ;
431 do_automask = 1 ; nopt++ ; continue ;
432 }
433
434 if( strcasecmp(argv[nopt],"-prefix") == 0 ){
435 WARNING_message("-prefix option is converted to -Mean in 3dTcorrMap") ;
436 Mprefix = argv[++nopt] ; nout++ ;
437 if( !THD_filename_ok(Mprefix) ) ERROR_exit("Illegal string after -prefix!\n") ;
438 nopt++ ; continue ;
439 }
440 if( strcasecmp(argv[nopt],"-Mean") == 0 || strcasecmp(argv[nopt],"-Mmean") == 0 ){
441 Mprefix = argv[++nopt] ; nout++ ;
442 if( !THD_filename_ok(Mprefix) ) ERROR_exit("Illegal prefix after %s!\n",argv[nopt-1]) ;
443 nopt++ ; continue ;
444 }
445 if( strcasecmp(argv[nopt],"-Zmean") == 0 ){
446 Zprefix = argv[++nopt] ; nout++ ;
447 if( !THD_filename_ok(Zprefix) ) ERROR_exit("Illegal prefix after -Zmean!\n") ;
448 nopt++ ; continue ;
449 }
450 if( strcasecmp(argv[nopt],"-Qmean") == 0 ){
451 Qprefix = argv[++nopt] ; nout++ ;
452 if( !THD_filename_ok(Qprefix) ) ERROR_exit("Illegal prefix after -Qmean!\n") ;
453 nopt++ ; continue ;
454 }
455 if( strcasecmp(argv[nopt],"-Pmean") == 0 ){
456 Pprefix = argv[++nopt] ; nout++ ;
457 if( !THD_filename_ok(Pprefix) ) ERROR_exit("Illegal prefix after -Pmean!\n") ;
458 nopt++ ; continue ;
459 }
460 if( strcasecmp(argv[nopt],"-Thresh") == 0 ){
461 Thresh = (float)strtod(argv[++nopt],NULL) ;
462 if( Thresh <= 0.0f || Thresh >= 0.99f )
463 ERROR_exit("Illegal -Thresh value %g",Thresh) ;
464 Tprefix = argv[++nopt] ; nout++ ;
465 if( !THD_filename_ok(Tprefix) )
466 ERROR_exit("Illegal prefix after -Thresh!\n") ;
467 nopt++ ; need_acc = 1 ; continue ;
468 }
469 if( strcasecmp(argv[nopt],"-VarThresh") == 0 ||
470 strcasecmp(argv[nopt],"-VarThreshN") == 0 ){
471
472 if (nopt+4 >= argc) {
473 ERROR_exit("Need three values and a prefix after -VarThresh*");
474 }
475
476 Tbone = !strcasecmp(argv[nopt],"-VarThreshN");
477
478 t0 = (float)strtod(argv[++nopt],NULL) ;
479 t1 = (float)strtod(argv[++nopt],NULL) ;
480 ti = (float)strtod(argv[++nopt],NULL) ;
481 N_iv = (int)((t1-t0)/ti)+1;
482 if( t0 <= 0.0f || t0 >= 0.99f )
483 ERROR_exit("Illegal 1st value of %g after -VarThresh* ",t0) ;
484 if( t1 <= 0.0f || t1 >= 0.99f )
485 ERROR_exit("Illegal 2nd value of %g after -VarThresh* ",t1) ;
486 if( ti <= 0.0f || ti >= 0.99f )
487 ERROR_exit("Illegal 3rd value of %g after -VarThresh* ",ti) ;
488 if (N_iv <= 0)
489 ERROR_exit("Bad combination of values after -VarThresh* ") ;
490
491 Threshv = (float *)calloc(N_iv+1, sizeof(float));
492 for (iv=0; iv < N_iv; ++iv)
493 Threshv[iv] = t0 + (float)iv*ti;
494
495 Tvprefix = argv[++nopt] ; nout++ ;
496 if( !THD_filename_ok(Tvprefix) )
497 ERROR_exit("Illegal prefix after %s",argv[nopt-4]) ;
498 if( Tvprefix[0] == '-' )
499 WARNING_message("%s prefix '%s' starts with '-' :: is this a mistake?",
500 argv[nopt-4] , Tvprefix ) ;
501
502 INFO_message("VarThresh mode with %d levels: %.3f .. %.3f\n",
503 N_iv , Threshv[0] , Threshv[N_iv-1] ) ;
504
505 nopt++ ; need_acc = 1 ; continue ;
506 }
507 if( strcasecmp(argv[nopt],"-CorrMask") == 0 ){ /* 25 Feb 2011 */
508 COmask++ ; nopt++ ; continue ;
509 }
510 if( strcasecmp(argv[nopt],"-CorrMap") == 0 ){
511
512 if (nopt+1 >= argc)
513 ERROR_exit("Need a prefix after -CorrMap");
514
515 COprefix = argv[++nopt] ; nout++ ;
516 if( !THD_filename_ok(COprefix) )
517 ERROR_exit("Illegal prefix after %s",argv[nopt-1]) ;
518 if( COprefix[0] == '-' )
519 WARNING_message(
520 "%s prefix '%s' starts with '-' :: is this a mistake?",
521 argv[nopt-1] , COprefix ) ;
522
523 nopt++ ; continue ;
524 }
525
526 if( strcasecmp(argv[nopt],"-Aexpr") == 0 ||
527 strcasecmp(argv[nopt],"-Cexpr") == 0 ||
528 strcasecmp(argv[nopt],"-Sexpr") == 0 ){
529
530 if( expr_type != '\0' ) ERROR_exit("Can't have 2 'expr' options!") ;
531 expr_type = argv[nopt][1] ; nopt++ ;
532 expr_string = strdup(argv[nopt]) ;
533 expr_code = PARSER_generate_code( expr_string ) ;
534 if( expr_code == NULL )
535 ERROR_exit("Illegal expression in option '%s'",argv[nopt-1]) ;
536 expr_has_r = PARSER_has_symbol("r",expr_code ) ;
537 expr_has_z = PARSER_has_symbol("z",expr_code ) ;
538 if( !expr_has_r && !expr_has_z )
539 ERROR_exit("Expression doesn't have 'r' and/or 'z' symbol!") ;
540 Eprefix = argv[++nopt] ; nout++ ;
541 if( !THD_filename_ok(Eprefix) )
542 ERROR_exit("Illegal prefix after '%s'",argv[nopt-2]) ;
543 nopt++ ; continue ;
544 }
545
546 if( strcasecmp(argv[nopt],"-Hist") == 0 ){
547 if( nopt+2 >= argc ) ERROR_exit("need 2 args after -Hist") ;
548 HHnum = (int)strtod(argv[++nopt],NULL) ;
549 if( HHnum < 20 ) ERROR_exit("-Hist count must be at least 20") ;
550 else if( HHnum > 1000 ) ERROR_exit("-Hist count must be at most 1000");
551 HHist = (short *)malloc(sizeof(short)*(HHnum+1)) ;
552 HHprefix = argv[++nopt] ; nout++ ;
553 HHdel = 2.0f / HHnum ; HHdin = HHnum / 2.0f ;
554 if( !THD_filename_ok(HHprefix) )
555 ERROR_exit("Illegal prefix after '%s'",argv[nopt-2]) ;
556 nopt++ ; continue ;
557 }
558
559 if( strcasecmp(argv[nopt],"-polort") == 0 ){
560 char *cpt ;
561 int val = (int)strtod(argv[++nopt],&cpt) ;
562 if( *cpt != '\0' || val < -1 || val > 19 )
563 ERROR_exit("Illegal value '%s' after -polort!",argv[nopt]) ;
564 polort = val ; nopt++ ; continue ;
565 }
566
567 if( strcasecmp(argv[nopt],"-bpass") == 0 ){
568 if( nopt+2 >= argc ) ERROR_exit("need 2 args after -bpass") ;
569 bpass_L = (float)strtod(argv[++nopt],NULL) ;
570 bpass_H = (float)strtod(argv[++nopt],NULL) ;
571 if( bpass_L < 0.0f || bpass_H <= bpass_L )
572 ERROR_exit("Illegal values after -bpass") ;
573 nopt++ ; continue ;
574 }
575
576 if( strcasecmp(argv[nopt],"-input") == 0 ){
577 if( xset != NULL ) ERROR_exit("Can't use -input twice!") ;
578 xset = THD_open_dataset(argv[++nopt]); CHECK_OPEN_ERROR(xset,argv[nopt]);
579 nopt++ ; continue ;
580 }
581
582 if( (strcasecmp(argv[nopt],"-base") == 0) ||
583 (strcasecmp(argv[nopt],"-seed") == 0) ){
584 /* Added -seed to fit with -help text. Don't know how -base
585 got in here. ZSS: Emergency patch. Jan 25 09 */
586 if( sset != NULL ) ERROR_exit("Can't use -seed twice!") ;
587 if( Mseedr != 0.0f ) ERROR_exit("Can't use -seed with -Mseed!") ;
588 sset = THD_open_dataset(argv[++nopt]); CHECK_OPEN_ERROR(sset,argv[nopt]);
589 nopt++ ; continue ;
590 }
591
592 if( strcasecmp(argv[nopt],"-mask") == 0 ){
593 THD_3dim_dataset *mset ;
594 if( do_automask ) ERROR_exit("Can't use -automask and -mask") ;
595 if( mask != NULL ) ERROR_exit("Can't use -mask twice") ;
596 if( ++nopt >= argc ) ERROR_exit("Need argument after '-mask'") ;
597 mset = THD_open_dataset( argv[nopt] ); CHECK_OPEN_ERROR(mset,argv[nopt]);
598 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
599 nxmask = DSET_NX(mset); nymask = DSET_NY(mset); nzmask = DSET_NZ(mset);
600 mask = THD_makemask( mset , 0 , 0.5f, 0.0f ) ; DSET_delete(mset) ;
601 if( mask == NULL )
602 ERROR_exit("Can't make mask from dataset '%s'",argv[nopt]) ;
603 nmask = THD_countmask( nxmask*nymask*nzmask , mask ) ;
604 INFO_message("Number of voxels in mask = %d",nmask) ;
605 if( nmask < 9 ) ERROR_exit("Mask is too small to process") ;
606 nopt++ ; continue ;
607 }
608
609 ERROR_message("Unknown option: %s",argv[nopt]) ;
610 suggest_best_prog_option(argv[0], argv[nopt]);
611 exit(1);
612 }
613 if( argc < 2 ){
614 ERROR_message("Too few options");
615 usage_3dTcorrMap(0);
616 exit(1) ;
617 }
618
619
620 /*-- open dataset, check for legality --*/
621
622 if( nout == 0 ) ERROR_exit("Don't you want any output datasets?") ;
623
624 if( nopt >= argc && xset == NULL)
625 ERROR_exit("No -input option and no dataset in final position?!") ;
626
627 if( xset == NULL ){
628 xset = THD_open_dataset(argv[nopt]); CHECK_OPEN_ERROR(xset,argv[nopt]);
629 }
630 ntime = DSET_NVALS(xset) ;
631 if( ntime < 9 )
632 ERROR_exit("Input dataset '%s' is too short: %d time points" ,
633 DSET_HEADNAME(xset) , ntime ) ;
634 if( ortim != NULL && ortim->nx < ntime )
635 ERROR_exit("-ort file is shorter than input dataset!") ;
636
637 DSET_UNMSEC(xset) ;
638 dtime = DSET_TR(xset) ;
639 if( dtime <= 0.0f ) dtime = 1.0f ;
640 do_bpass = (bpass_L < bpass_H) ;
641 if( do_bpass ){
642 kk = THD_bandpass_OK( ntime , dtime , bpass_L , bpass_H , 1 ) ;
643 if( kk <= 0 ) ERROR_exit("Can't continue since -bpass setup is illegal") ;
644 polort = -1 ;
645 }
646
647 if( sset != NULL ){
648 if( DSET_NVALS(sset) != ntime )
649 ERROR_exit("-seed dataset time series length %d doesn't match -input %d",
650 DSET_NVALS(sset) , ntime ) ;
651
652 if( DSET_NVOX(sset) != DSET_NVOX(xset) )
653 ERROR_exit("-seed dataset doesn't match -input dataset in space") ;
654 }
655
656 /*-- compute references, if any --*/
657
658 if( polort >= 0 ){
659 nref = polort + 1 ; ref = THD_build_polyref( nref , ntime ) ;
660 }
661 if( ortim != NULL ){
662 jj = ortim->ny ; ortar = MRI_FLOAT_PTR(ortim) ;
663 ref = (float **)realloc( ref , sizeof(float *)*(nref+jj) ) ;
664 for( ii=0 ; ii < jj ; ii++ )
665 ref[ii+nref] = ortar + (ii*ortim->nx) ;
666 nref += jj ;
667 }
668
669 /*-- compute mask array, if desired --*/
670
671 nx = DSET_NX(xset) ;
672 ny = DSET_NY(xset) ; nxy = nx*ny ;
673 nz = DSET_NZ(xset) ; nxyz = nvox = nxy*nz ;
674
675 if( do_automask ){
676 if (!DSET_IS_VOL(xset)) {
677 ERROR_exit("Can't use -automask with non-volumetric input datasets");
678 }
679 mask = THD_automask( xset ) ;
680 nmask = THD_countmask( nvox , mask ) ;
681 if( nmask > 9 )
682 INFO_message("%d voxels survive -automask",nmask) ;
683 else
684 ERROR_exit("only %d voxels survive -automask",nmask) ;
685 } else if( mask != NULL ){
686 if( nxmask != nx || nymask != ny || nzmask != nz )
687 ERROR_exit("-mask and -input datasets differ in voxel grids!") ;
688 } else {
689 nmask = nvox ;
690 mask = (byte *)malloc(sizeof(byte)*nmask) ;
691 memset( mask , 1 , sizeof(byte)*nmask ) ;
692 INFO_message("computing for all %d voxels!",nmask) ;
693 }
694
695 /*--- load input data ---*/
696
697 if( !DSET_LOADED(xset) ){
698 INFO_message("Loading input dataset") ;
699 DSET_load(xset) ; CHECK_LOAD_ERROR(xset) ;
700 }
701
702 /* check for constant voxels, and remove them from the mask */
703
704 xsar = (float *)malloc( sizeof(float)*ntime ) ;
705 for( ii=0 ; ii < nvox ; ii++ ){
706 if( mask[ii] == 0 ) continue ;
707 (void)THD_extract_array( ii , xset , 0 , xsar ) ;
708 for( kk=1 ; kk < ntime && xsar[kk] == xsar[0] ; kk++ ) ; /*nada*/
709 if( kk == ntime ) mask[ii] = 0 ; /* xsar is constant */
710 }
711 free(xsar) ;
712 ii = THD_countmask( nvox , mask ) ;
713 if( ii < nmask ){
714 if( ii > 9 ) ININFO_message("only %d voxels in dataset are non-constant" ,ii);
715 else ERROR_exit ("only %d voxels in dataset are non-constant!?",ii);
716 nmask = ii ;
717 }
718
719 /* create indx[jj] = voxel index in dataset whence
720 jj-th extracted time series comes from */
721
722 indx = (int *)malloc(sizeof(int)*nmask) ;
723 for( ii=jj=0 ; ii < nvox ; ii++ ) if( mask[ii] ) indx[jj++] = ii ;
724
725 /*-- check PCort mask (if any) for matchingness --*/
726
727 if( PCortn > 0 ){
728 if( PCnx != nx || PCny != ny || PCnz != nz )
729 ERROR_exit("-PCort mask doesn't match -input dataset voxel grid!") ;
730 }
731
732 /*-- create output datasets --*/
733
734 if( Mprefix != NULL ){
735 Mset = EDIT_empty_copy( xset ) ;
736 EDIT_dset_items( Mset ,
737 ADN_prefix , Mprefix ,
738 ADN_nvals , 1 ,
739 ADN_ntt , 0 ,
740 ADN_brick_fac , NULL ,
741 ADN_type , HEAD_FUNC_TYPE ,
742 ADN_func_type , FUNC_BUCK_TYPE ,
743 ADN_none ) ;
744 EDIT_substitute_brick( Mset , 0 , MRI_float , NULL ) ;
745 Mar = DSET_ARRAY(Mset,0) ; /* get array */
746 EDIT_BRICK_TO_NOSTAT(Mset,0) ;
747 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(Mset)) )
748 ERROR_exit("Output dataset %s already exists!",
749 DSET_HEADNAME(Mset)) ;
750 tross_Make_History( "3dTcorrMap" , argc,argv , Mset ) ;
751 }
752
753 if( COprefix != NULL ){
754 int iii, jjj, kkk, ijk , dig=0 , nv=(COmask) ? nmask : nxyz ;
755 float *lesfac=(float*)calloc(nv,sizeof(float));
756 for (ii=0; ii<nv; ++ii) lesfac[ii]=1/10000.0;
757 COset = EDIT_empty_copy( xset ) ;
758 EDIT_dset_items( COset ,
759 ADN_prefix , COprefix ,
760 ADN_nvals , nv ,
761 ADN_ntt , nv ,
762 ADN_ttdel , 1.0 ,
763 ADN_nsl , 0 ,
764 ADN_brick_fac , lesfac ,
765 ADN_type , HEAD_FUNC_TYPE,
766 ADN_func_type , FUNC_FIM_TYPE ,
767 ADN_none ) ;
768
769 INFO_message ("CorrMap: this might take a long while and a lot of memory,\n") ;
770 ININFO_message(" and will create big [%s bytes] output dataset ...",
771 approximate_number_string((double)DSET_TOTALBYTES(COset)) ) ;
772 free(lesfac); lesfac = NULL;
773 dig = (int)ceil(log((double)nv)/log(10.0));
774 for (ii=0; ii<nv; ++ii) {
775 EDIT_substitute_brick( COset , ii , MRI_short , NULL ) ;
776 EDIT_BRICK_TO_NOSTAT(COset,ii) ;
777 if (!DSET_IS_VOL(xset)) {
778 switch (dig) {
779 case 6: sprintf(stmp,"n%06d", ii); break;
780 case 5: sprintf(stmp,"n%05d", ii); break;
781 case 4: sprintf(stmp,"n%04d", ii); break;
782 case 3:
783 case 2:
784 case 1: sprintf(stmp,"n%03d", ii); break;
785 default: sprintf(stmp,"n%07d", ii); break;
786 }
787 } else {
788 ijk = (COmask) ? indx[ii] : ii ;
789 kkk = ijk / nxy ; jjj = ijk % nxy ; iii = jjj % nx ; jjj = jjj / nx ;
790 sprintf(stmp,"v%03d.%03d.%03d",iii,jjj,kkk) ;
791 }
792 EDIT_BRICK_LABEL(COset,ii,stmp) ;
793 }
794 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(COset)) )
795 ERROR_exit("Output dataset %s already exists!",
796 DSET_HEADNAME(COset)) ;
797 tross_Make_History( "3dTcorrMap" , argc,argv , COset ) ;
798 }
799
800 if( Pprefix != NULL ){
801 Pset = EDIT_empty_copy( xset ) ;
802 EDIT_dset_items( Pset ,
803 ADN_prefix , Pprefix ,
804 ADN_nvals , 1 ,
805 ADN_ntt , 0 ,
806 ADN_brick_fac , NULL ,
807 ADN_type , HEAD_FUNC_TYPE ,
808 ADN_func_type , FUNC_BUCK_TYPE ,
809 ADN_none ) ;
810 EDIT_substitute_brick( Pset , 0 , MRI_float , NULL ) ;
811 Par = DSET_ARRAY(Pset,0) ; /* get array */
812 EDIT_BRICK_TO_NOSTAT(Pset,0) ;
813 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(Pset)) )
814 ERROR_exit("Output dataset %s already exists!",
815 DSET_HEADNAME(Pset)) ;
816 tross_Make_History( "3dTcorrMap" , argc,argv , Pset ) ;
817 }
818
819 if( Eprefix != NULL ){
820 Eset = EDIT_empty_copy( xset ) ;
821 EDIT_dset_items( Eset ,
822 ADN_prefix , Eprefix ,
823 ADN_nvals , 1 ,
824 ADN_ntt , 0 ,
825 ADN_brick_fac , NULL ,
826 ADN_type , HEAD_FUNC_TYPE ,
827 ADN_func_type , FUNC_BUCK_TYPE ,
828 ADN_none ) ;
829 EDIT_substitute_brick( Eset , 0 , MRI_float , NULL ) ;
830 Ear = DSET_ARRAY(Eset,0) ; /* get array */
831 EDIT_BRICK_TO_NOSTAT(Eset,0) ;
832 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(Eset)) )
833 ERROR_exit("Output dataset %s already exists!",
834 DSET_HEADNAME(Eset)) ;
835 tross_Make_History( "3dTcorrMap" , argc,argv , Eset ) ;
836 }
837
838 if( Zprefix != NULL ){
839 Zset = EDIT_empty_copy( xset ) ;
840 EDIT_dset_items( Zset ,
841 ADN_prefix , Zprefix ,
842 ADN_nvals , 1 ,
843 ADN_ntt , 0 ,
844 ADN_brick_fac , NULL ,
845 ADN_type , HEAD_FUNC_TYPE ,
846 ADN_func_type , FUNC_BUCK_TYPE ,
847 ADN_none ) ;
848 EDIT_substitute_brick( Zset , 0 , MRI_float , NULL ) ;
849 Zar = DSET_ARRAY(Zset,0) ; /* get array */
850 EDIT_BRICK_TO_NOSTAT(Zset,0) ;
851 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(Zset)) )
852 ERROR_exit("Output dataset %s already exists!",
853 DSET_HEADNAME(Zset)) ;
854 tross_Make_History( "3dTcorrMap" , argc,argv , Zset ) ;
855 }
856
857 if( Qprefix != NULL ){
858 Qset = EDIT_empty_copy( xset ) ;
859 EDIT_dset_items( Qset ,
860 ADN_prefix , Qprefix ,
861 ADN_nvals , 1 ,
862 ADN_ntt , 0 ,
863 ADN_brick_fac , NULL ,
864 ADN_type , HEAD_FUNC_TYPE ,
865 ADN_func_type , FUNC_BUCK_TYPE ,
866 ADN_none ) ;
867 EDIT_substitute_brick( Qset , 0 , MRI_float , NULL ) ;
868 Qar = DSET_ARRAY(Qset,0) ; /* get array */
869 EDIT_BRICK_TO_NOSTAT(Qset,0) ;
870 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(Qset)) )
871 ERROR_exit("Output dataset %s already exists!",
872 DSET_HEADNAME(Qset)) ;
873 tross_Make_History( "3dTcorrMap" , argc,argv , Qset ) ;
874 }
875
876 if( Tprefix != NULL ){
877 Tset = EDIT_empty_copy( xset ) ;
878 EDIT_dset_items( Tset ,
879 ADN_prefix , Tprefix ,
880 ADN_nvals , 1 ,
881 ADN_ntt , 0 ,
882 ADN_brick_fac , NULL ,
883 ADN_type , HEAD_FUNC_TYPE ,
884 ADN_func_type , FUNC_BUCK_TYPE ,
885 ADN_none ) ;
886 EDIT_substitute_brick( Tset , 0 , MRI_float , NULL ) ;
887 Tar = DSET_ARRAY(Tset,0) ; /* get array */
888 EDIT_BRICK_TO_NOSTAT(Tset,0) ;
889 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(Tset)) )
890 ERROR_exit("Output dataset %s already exists!",
891 DSET_HEADNAME(Tset)) ;
892 tross_Make_History( "3dTcorrMap" , argc,argv , Tset ) ;
893 }
894
895 if( Tvprefix != NULL ){
896 Tvset = EDIT_empty_copy( xset ) ;
897 EDIT_dset_items( Tvset ,
898 ADN_prefix , Tvprefix ,
899 ADN_nvals , N_iv ,
900 ADN_ntt , 0 ,
901 ADN_brick_fac , NULL ,
902 ADN_type , HEAD_FUNC_TYPE ,
903 ADN_func_type , FUNC_BUCK_TYPE ,
904 ADN_none ) ;
905 Tvar = (float **)calloc(N_iv,sizeof(float*));
906 Tvcount = (float * )calloc(N_iv,sizeof(float ));
907 for( iv=0 ; iv < N_iv ; ++iv ){
908 EDIT_substitute_brick( Tvset , iv , MRI_float , NULL ) ;
909 Tvar[iv] = DSET_ARRAY(Tvset,iv) ; /* get array */
910 EDIT_BRICK_TO_NOSTAT(Tvset,iv) ;
911 if (Tbone) sprintf(stmp,"nTc%.3f", Threshv[iv]);
912 else sprintf(stmp,"Tc%.3f" , Threshv[iv]);
913 EDIT_BRICK_LABEL(Tvset, iv, stmp);
914 }
915 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(Tvset)) )
916 ERROR_exit("Output dataset %s already exists!",
917 DSET_HEADNAME(Tvset)) ;
918 tross_Make_History( "3dTcorrMap" , argc,argv , Tvset ) ;
919 }
920
921 if( HHprefix != NULL ){
922 HHset = EDIT_empty_copy( xset ) ;
923 EDIT_dset_items( HHset ,
924 ADN_prefix , HHprefix ,
925 ADN_nvals , HHnum ,
926 ADN_ntt , 0 ,
927 ADN_brick_fac , NULL ,
928 ADN_type , HEAD_FUNC_TYPE ,
929 ADN_func_type , FUNC_BUCK_TYPE ,
930 ADN_none ) ;
931 for( iv=0 ; iv < HHnum ; iv++ ){
932 EDIT_substitute_brick( HHset , iv , MRI_short , NULL ) ;
933 EDIT_BRICK_TO_NOSTAT(HHset,iv) ;
934 sprintf(stmp,"%+6.3f:%+6.3f",-1.0+iv*HHdel,-1.0+(iv+1)*HHdel) ;
935 EDIT_BRICK_LABEL(HHset,iv,stmp) ;
936 }
937 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(HHset)) )
938 ERROR_exit("Output dataset %s already exists!",
939 DSET_HEADNAME(HHset)) ;
940 tross_Make_History( "3dTcorrMap" , argc,argv , HHset ) ;
941 }
942
943 /* make neighborhood mask for -Mseed, if needed */
944
945 dx = fabsf(DSET_DX(xset)) ;
946 dy = fabsf(DSET_DY(xset)) ;
947 dz = fabsf(DSET_DZ(xset)) ;
948
949 if ( (Mseedr > 0.0f || Mseedr < 0.0f) &&
950 !DSET_IS_VOL(xset)) {
951 ERROR_exit("Can't use -Mseed with non-volumetric input datasets");
952 }
953 if( Mseedr > 0.0f){
954 Mseed_nbhd = MCW_spheremask( dx,dy,dz , Mseedr ) ;
955 if( Mseed_nbhd == NULL || Mseed_nbhd->num_pt < 2 )
956 ERROR_exit("Can't build -Mseed neighborhood mask!?") ;
957 ININFO_message("%d voxels in -Mseed neighborhood",Mseed_nbhd->num_pt) ;
958 } else if( Mseedr < 0.0f ){
959 Mseedr = -Mseedr ; if( Mseedr < 1.01f ) Mseedr = 1.01f ;
960 Mseed_nbhd = MCW_spheremask( 1.0f,1.0f,1.0f , Mseedr ) ;
961 if( Mseed_nbhd == NULL || Mseed_nbhd->num_pt < 2 )
962 ERROR_exit("Can't build -Mseed neighborhood mask!?") ;
963 ININFO_message("%d voxels in -Mseed neighborhood",Mseed_nbhd->num_pt) ;
964 }
965
966 if( Mseedr > 0.0f ) Mseedar = (float *)malloc(sizeof(float)*ntime) ;
967
968 /* Gblur inside the mask */
969
970 if( gblur > 0.0f ){
971 if (DSET_IS_VOL(xset)) {
972 float fx,fy,fz ; int nrep=0 ;
973 MRI_IMAGE *bim ; float *bar ;
974
975 if( nx == 1 ) dx = 0.0f ;
976 if( ny == 1 ) dy = 0.0f ;
977 if( nz == 1 ) dz = 0.0f ;
978
979 mri_blur3D_getfac( gblur , dx,dy,dz , &nrep , &fx,&fy,&fz ) ;
980 if( nrep > 0 && fx > 0.0f && fy > 0.0f && fz >= 0.0f ){
981 ININFO_message("-Gblur-ing parameters: "
982 "#iter=%d fx=%.5f fy=%.5f fz=%.5f" ,
983 nrep , fx , fy , fz ) ;
984 for( iv=0 ; iv < ntime ; iv++ ){
985 bim = THD_extract_float_brick( iv , xset ) ;
986 mri_blur3D_inmask( bim , mask , fx,fy,fz , nrep ) ;
987 EDIT_substitute_brick( xset , iv , MRI_float , MRI_FLOAT_PTR(bim) ) ;
988 EDIT_BRICK_FACTOR( xset,iv,0.0f ) ;
989 mri_clear_data_pointer(bim) ; mri_free(bim) ;
990 }
991 } else {
992 WARNING_message("Can't -Gblur input dataset for some reason") ;
993 }
994 } else {
995 if (!DSET_IS_VOL(xset)) {
996 ERROR_exit("Can't use -Gblur with non-volumetric "
997 "input datasets");
998 }
999 }
1000 }
1001
1002 /* remove trends (-polort and -ort) now (if not bandpassing later) */
1003
1004 if( nref > 0 && !do_bpass ){
1005 THD_3dim_dataset *yset ; byte *rmask=mask ;
1006 ININFO_message("Detrending input dataset with %d references",nref) ;
1007 if( PCmask != NULL ){
1008 rmask = (byte *)malloc(sizeof(byte)*nvox) ;
1009 for( ii=0 ; ii < nvox ; ii++ ) rmask[ii] = (mask[ii] || PCmask[ii]) ;
1010 }
1011 yset = THD_detrend_dataset( xset , nref,ref , 2,0 , rmask , NULL ) ;
1012 if( yset == NULL ) ERROR_exit("Detrending fails!?") ;
1013 DSET_delete(xset) ; xset = yset ;
1014 if( rmask != mask ) free(rmask) ;
1015 }
1016
1017 /* extract dataset time series into an array of
1018 time series vectors, for ease and speed of access */
1019
1020 ININFO_message("extracting mask-ed time series") ;
1021 timar = THD_extract_many_series( nmask , indx , xset ) ;
1022 if( timar == NULL ) ERROR_exit("-input extraction failed!?") ;
1023
1024 /* extract PCmask time series from xset */
1025 /* if -polort was invoked, then these will already be detrended above */
1026
1027 if( PCortn > 0 ){
1028 int *pcdx = (int *)malloc(sizeof(int)*PCnmask) ;
1029 for( ii=jj=0 ; ii < nvox ; ii++ ) if( PCmask[ii] ) pcdx[jj++] = ii ;
1030 PCimar = THD_extract_many_series( PCnmask , pcdx , xset ) ;
1031 if( PCimar == NULL ) /* should not happen */
1032 ERROR_exit("Can't get -PCort data from -input dataset!?") ;
1033 free(pcdx) ;
1034 }
1035
1036 /* done with -input dataset */
1037
1038 DSET_delete(xset) ;
1039
1040 /* if -seed was given, get that set of time series as well */
1041
1042 if( sset != NULL ){
1043 DSET_load(sset) ; CHECK_LOAD_ERROR(sset) ;
1044 simar = THD_extract_many_series( nmask , indx , sset ) ;
1045 if( simar == NULL ) ERROR_exit("-seed extraction failed!?") ;
1046 DSET_delete(sset) ;
1047 }
1048
1049 /* only need the mask later if doing spatial
1050 processing inside the voxel processing loop */
1051
1052 if( Mseedr == 0.0f ){ free(mask); mask = NULL; }
1053
1054 /* do bandpass now, if so desired */
1055
1056 if( do_bpass ){
1057 float **vec = (float **)malloc(sizeof(float *)*nmask) ;
1058 ININFO_message("Bandpassing input dataset") ;
1059 for( kk=0 ; kk < nmask ; kk++ )
1060 vec[kk] = MRI_FLOAT_PTR( IMARR_SUBIM(timar,kk) ) ;
1061
1062 (void)THD_bandpass_vectors( ntime , nmask , vec ,
1063 dtime , bpass_L , bpass_H ,
1064 2 , nref , ref ) ;
1065
1066 /* bandpass PCort vectors also */
1067
1068 if( PCortn > 0 ){
1069 vec = (float **)realloc(vec,sizeof(float *)*PCnmask) ;
1070 for( kk=0 ; kk < PCnmask ; kk++ )
1071 vec[kk] = MRI_FLOAT_PTR( IMARR_SUBIM(PCimar,kk) ) ;
1072
1073 (void)THD_bandpass_vectors( ntime , PCnmask , vec ,
1074 dtime , bpass_L , bpass_H ,
1075 2 , nref , ref ) ;
1076 }
1077
1078 free(vec) ;
1079 }
1080
1081 /* Now remove PCort trends from timar */
1082
1083 if( PCortn > 0 ){
1084 float *uvec = (float *)malloc(sizeof(float)*ntime*PCortn) ;
1085 float *svec = (float *)malloc(sizeof(float)*PCortn) ;
1086 int nvec ; register float sum ; float *usar ;
1087
1088 nvec = mri_principal_vectors( PCimar , PCortn , svec,uvec ) ;
1089 if( nvec <= 0 )
1090 ERROR_exit("Can't compute -PCort SVD!?") ;
1091 else if( nvec < PCortn )
1092 WARNING_message("-PCort: only have %d components, not %d",nvec,PCortn) ;
1093 DESTROY_IMARR(PCimar) ;
1094
1095 fprintf(stderr,"++ -PCort singular values:") ;
1096 for( kk=0 ; kk < nvec ; kk++ ) fprintf(stderr," %.4g",svec[kk]) ;
1097 fprintf(stderr,"\n") ; free(svec) ;
1098
1099 /* in the loop below, recall that the vectors [usar]
1100 are L2-normalized and orthogonal, so they can be
1101 projected out sequentially and by a simple dot product */
1102
1103 for( ii=0 ; ii < nmask ; ii++ ){
1104 xsar = MRI_FLOAT_PTR( IMARR_SUBIM(timar,kk) ) ;
1105 for( kk=0 ; kk < nvec ; kk++ ){
1106 sum = 0.0f ; usar = uvec + kk*ntime ;
1107 for( jj=0 ; jj < ntime ; jj++ ) sum += xsar[jj]*usar[jj] ;
1108 for( jj=0 ; jj < ntime ; jj++ ) xsar[jj] -= sum*usar[jj] ;
1109 }
1110 }
1111 free(uvec) ; free(PCmask) ;
1112 }
1113
1114 /* normalize -input vectors so sum of squares is 1,
1115 for speed in computing correlations (and we need speed) */
1116
1117 ININFO_message("normalizing extracted time series") ;
1118 for( ii=0 ; ii < nmask ; ii++ ){
1119 xsar = MRI_FLOAT_PTR( IMARR_SUBIM(timar,ii) ) ;
1120 csum = 0.0f ;
1121 for( jj=0 ; jj < ntime ; jj++ ) csum += xsar[jj]*xsar[jj] ;
1122 if( csum > 0.0f ){ /* should always be true */
1123 csum = 1.0f / sqrtf(csum) ;
1124 for( jj=0 ; jj < ntime ; jj++ ) xsar[jj] *= csum ;
1125 }
1126 }
1127 if( simar != NULL ){ /* ditto for -seed vectors */
1128 for( ii=0 ; ii < nmask ; ii++ ){
1129 xsar = MRI_FLOAT_PTR( IMARR_SUBIM(simar,ii) ) ;
1130 csum = 0.0f ;
1131 for( jj=0 ; jj < ntime ; jj++ ) csum += xsar[jj]*xsar[jj] ;
1132 if( csum > 0.0f ){ /* should always be true */
1133 csum = 1.0f / sqrtf(csum) ;
1134 for( jj=0 ; jj < ntime ; jj++ ) xsar[jj] *= csum ;
1135 }
1136 }
1137 }
1138
1139 /*-------------------------------------------------------------------*/
1140 /*--- loop over voxels, correlate (lots and lots of CPU time now) ---*/
1141
1142 vstep = (nmask > 999) ? nmask/50 : 0 ;
1143
1144 isodd = (ntime%2 == 1) ;
1145 ccar = (float *)malloc(sizeof(float)*nmask) ; /* 29 Apr 2009 */
1146
1147 for( ii=0 ; ii < 26 ; ii++ ) atoz[ii] = NULL ;
1148 if( expr_type != '\0' ){ /* 23 Jun 2009 */
1149 eear = (double *)malloc(sizeof(double)*nmask) ;
1150 if( expr_has_r )
1151 atoz[17] = (double *)malloc(sizeof(double)*nmask) ;
1152 if( expr_has_z )
1153 atoz[25] = (double *)malloc(sizeof(double)*nmask) ;
1154 }
1155
1156 /* initialize timer */
1157 etime(&tt,0);
1158
1159 /* 29 Apr 2009: print # of threads message to amuse the user */
1160
1161 ctime = 2.0 * (double)nmask * (double)nmask * (double)ntime
1162 + 100.0 * (double)nmask * (double)ntime ;
1163 if( HHnum > 0 )
1164 ctime += HHnum * (double)nmask * (double)nmask ;
1165 if( expr_type != '\0' )
1166 ctime += 111.0 * (double)nmask * (double)nmask ;
1167 if( Threshv )
1168 ctime += 2.0 * N_iv * (double)nmask * (double)nmask ;
1169
1170 #ifdef USE_OMP
1171 #pragma omp parallel
1172 {
1173 if( omp_get_thread_num() == 0 )
1174 INFO_message("Starting long long loop: OpenMP threads=%d %s Flops" ,
1175 omp_get_num_threads() , approximate_number_string(ctime) ) ;
1176 }
1177 #else
1178 INFO_message("Starting long long loop through all voxels: %s Flops" ,
1179 approximate_number_string(ctime) ) ;
1180 ININFO_message(" [You would be happier if you used an OpenMP-ized version]") ;
1181 ININFO_message(" [of this program, and had a multiple CPU computer system.]") ;
1182 #endif
1183
1184 for( ii=0 ; ii < nmask ; ii++ ){ /* outer loop over voxels: */
1185 /* seed time series to correlate with */
1186
1187 if( vstep && ii%vstep == vstep-1 ){ /* palliate the user's pain */
1188 if( ii < vstep ){
1189 dtt = etime(&tt,1) ;
1190 ININFO_message("Single loop (1/50) duration: %.1f secs = %.2f mins\n"
1191 " Remaining (49/50) time est.: %.0f secs = %.1f mins = %.2f hrs = %.3f days\n",
1192 dtt,dtt/60.0, dtt*49.0 , dtt/60.0*49.0 ,
1193 dtt/3600.0*49.0 , dtt/86400.0*49.0 ) ;
1194 fprintf(stderr,"++ Voxel loop: ") ;
1195 }
1196 vstep_print() ;
1197 }
1198
1199 /* get the ii-th time series */
1200
1201 if( simar != NULL ){ /* from -seed dataset */
1202 xsar = MRI_FLOAT_PTR( IMARR_SUBIM(simar,ii) ) ;
1203
1204 } else if( Mseedr == 0.0f ){ /* from -input dataset */
1205 xsar = MRI_FLOAT_PTR( IMARR_SUBIM(timar,ii) ) ;
1206
1207 } else { /* from -input dataset, then averaged */
1208 int mm , xc,yc,zc , xm,ym,zm , uu , nuu ;
1209 xsar = MRI_FLOAT_PTR( IMARR_SUBIM(timar,ii) ) ; /* central voxel */
1210 for( uu=0 ; uu < ntime ; uu++ ) Mseedar[uu] = xsar[uu] ;
1211 IJK_TO_THREE(indx[ii],xc,yc,zc,nx,nxy) ; /* 3D index of central voxel */
1212 for( nuu=mm=1 ; mm < Mseed_nbhd->num_pt ; mm++ ){
1213 xm = xc + Mseed_nbhd->i[mm] ; if( xm < 0 || xm >= nx ) continue ;
1214 ym = yc + Mseed_nbhd->j[mm] ; if( ym < 0 || ym >= ny ) continue ;
1215 zm = zc + Mseed_nbhd->k[mm] ; if( zm < 0 || zm >= nz ) continue ;
1216 jj = THREE_TO_IJK(xm,ym,zm,nx,nxy) ; /* 3D index of offset voxel */
1217 if( !mask[jj] ) continue ; /* not in mask ==> skip it */
1218 kk = mybsearch_int( jj , nmask , indx ) ; /* find in indx[] array */
1219 if( kk < 0 ) continue ; /* should not happen */
1220 xsar = MRI_FLOAT_PTR( IMARR_SUBIM(timar,kk) ) ; /* add add it in */
1221 for( uu=0 ; uu < ntime ; uu++ ) Mseedar[uu] += xsar[uu] ;
1222 nuu++ ; /* count of number of vectors added in */
1223 }
1224 if( nuu > 1 ) THD_normalize( ntime , Mseedar ) ; /* L2 normalize */
1225 xsar = Mseedar ; /* assign computed array to seed pointer */
1226 }
1227
1228 AFNI_OMP_START ;
1229 #pragma omp parallel
1230 { int vv,uu ; float *ysar ; float qcc ;
1231 #pragma omp for
1232 for( vv=0 ; vv < nmask ; vv++ ){ /*------- inner loop over voxels -------*/
1233
1234 if( vv==ii ){ ccar[vv] = 0.0f ; continue ; }
1235 ysar = MRI_FLOAT_PTR( IMARR_SUBIM(timar,vv) ) ;
1236
1237 /** dot products (unrolled by 2 on 29 Apr 2009) **/
1238
1239 if( isodd ){
1240 for( qcc=xsar[0]*ysar[0],uu=1 ; uu < ntime ; uu+=2 )
1241 qcc += xsar[uu]*ysar[uu] + xsar[uu+1]*ysar[uu+1] ;
1242 } else {
1243 for( qcc=0.0f,uu=0 ; uu < ntime ; uu+=2 )
1244 qcc += xsar[uu]*ysar[uu] + xsar[uu+1]*ysar[uu+1] ;
1245 }
1246 ccar[vv] = qcc ; /* save correlation in ccar for later (OpenMP mod) */
1247 } /* end of inner loop over voxels (vv) */
1248 } /* end OpenMP */
1249 AFNI_OMP_END ;
1250
1251 /** combine results in ccar to give output values **/
1252
1253 Tcount = Mcsum = Zcsum = Qcsum = 0.0f ;
1254 Pcsum = 0.0f ; nPcsum = 0 ;
1255 for( iv=0 ; iv < N_iv ; ++iv ) Tvcount[iv] = 0.0f ;
1256
1257 #if 0
1258 INFO_message("correlations for ii=%d indx=%d",ii,indx[ii]) ;
1259 for( jj=0 ; jj < nmask ; jj++ ) fprintf(stderr," %.3f",ccar[jj]) ;
1260 fprintf(stderr,"\n") ;
1261 #endif
1262
1263 if (COset) {
1264 int ijk = (COmask) ? ii : indx[ii] ;
1265 COar = DSET_ARRAY(COset,ijk);
1266 for( jj=0 ; jj < nmask ; jj++ ) {
1267 cc = ccar[jj] * 10000.0f; /* scale up because output is short */
1268 COar[indx[jj]] = cc < 0.0f ? (short)(cc-0.5f):(short)(cc+0.5f);
1269 }
1270 }
1271
1272 for( jj=0 ; jj < nmask ; jj++ ){
1273 if( jj == ii ) continue ; /* no self dealing [we're not in Congress] */
1274 cc = ccar[jj] ;
1275 Mcsum += cc ;
1276 Qcsum += cc*cc ;
1277 if( Par != NULL && cc > 0.0f ){ Pcsum += cc*cc ; nPcsum++ ; }
1278 if( Zar != NULL ) Zcsum += 0.5f * logf((1.0001f+cc)/(1.0001f-cc));
1279 if( need_acc ){
1280 acc = (cc < 0.0f) ? -cc : cc ;
1281 if( acc >= Thresh ) Tcount++ ;
1282 if( Threshv ){
1283 iv = N_iv - 1 ;
1284 while( iv > -1 ){
1285 if( acc >= Threshv[iv] ){
1286 do { Tvcount[iv--]++; } while (iv > -1);
1287 }
1288 --iv;
1289 }
1290 }
1291 }
1292 } /* end of combining */
1293
1294 if( Mar != NULL ) Mar[indx[ii]] = Mcsum / (nmask-1.0f) ;
1295 if( Zar != NULL ) Zar[indx[ii]] = tanh( Zcsum / (nmask-1.0f) ) ;
1296 if( Qar != NULL ) Qar[indx[ii]] = sqrt( Qcsum / (nmask-1.0f) ) ;
1297 if( Tar != NULL ) Tar[indx[ii]] = Tcount ;
1298 if( Tvar != NULL ){
1299 for(iv=0 ; iv < N_iv ; ++iv ) Tvar[iv][indx[ii]] = Tvcount[iv] ;
1300 }
1301
1302 if( Par != NULL ) /* 23 Jun 2009 */
1303 Par[indx[ii]] = Pcsum / MAX(1,nPcsum) ;
1304
1305 if( expr_type != '\0' ){ /* 23 Jun 2009 */
1306 if( expr_has_r ){
1307 for( jj=0 ; jj < nmask ; jj++ ) atoz[17][jj] = (double)ccar[jj] ;
1308 }
1309 if( expr_has_z ){
1310 for( jj=0 ; jj < nmask ; jj++ ) atoz[25][jj] = atanhf(ccar[jj]) ;
1311 }
1312 PARSER_evaluate_vector( expr_code , atoz , nmask, eear ) ;
1313 Esum = 0.0f ; nEsum = 0.0f ;
1314 for( jj=0 ; jj < nmask ; jj++ ){
1315 if( jj == ii ) continue ;
1316 Esum += eear[jj] ; if( eear[jj] != 0.0 ) nEsum++ ;
1317 }
1318 switch( expr_type ){
1319 case 'A': Ear[indx[ii]] = Esum / (nmask-1.0f) ; break ;
1320 case 'C': Ear[indx[ii]] = Esum / MAX(1,nEsum) ; break ;
1321 case 'S': Ear[indx[ii]] = Esum ; break ;
1322 }
1323 }
1324
1325 if( HHist != NULL ){
1326 memset( HHist , 0 , sizeof(short)*(HHnum+1) ) ;
1327 for( jj=0 ; jj < nmask ; jj++ ){
1328 if( jj == ii ) continue ;
1329 kk = (1.0f + ccar[jj])*HHdin ;
1330 if( HHist[kk] < MRI_maxshort ) HHist[kk]++ ;
1331 }
1332 THD_insert_series( indx[ii], HHset, HHnum, MRI_short, HHist, 1 ) ;
1333 }
1334
1335 } /* end of outer loop over voxels (ii) */
1336
1337 /*-------------------------------------------------------------------*/
1338
1339 if (Tbone) { /* scale by expected number of voxels by chance */
1340 float p[3], sc, pval;
1341 p[0] = (float)ntime;
1342 p[1] = 1.0;
1343 p[2] = (float)nref;
1344 for (iv=0; iv<N_iv; ++iv) {
1345 pval = THD_stat_to_pval (Threshv[iv], NI_STAT_CORREL, p);
1346 sc = (float) DSET_NVOX(Tvset) * pval;
1347 if (sc < 0.05) sc = 0.05;
1348 /*
1349 fprintf(stderr,"Scaling sb %d with threshold %.2f\n"
1350 "and p %.5f by %f\n",
1351 iv, Threshv[iv],
1352 pval, sc);
1353 */
1354 for( ii=0 ; ii < nmask ; ii++ ) Tvar[iv][indx[ii]] /= sc;
1355 }
1356 }
1357
1358 if( vstep ) fprintf(stderr,"!\n") ;
1359
1360 /*--- finito ---*/
1361
1362 free(indx) ; DESTROY_IMARR(timar) ; free(ccar) ;
1363 if( simar != NULL ) DESTROY_IMARR(simar) ;
1364 if( mask != NULL ) free(mask) ;
1365 if( Mseedar != NULL ) free(Mseedar) ;
1366
1367 if( atoz[17] ) free(atoz[17]) ;
1368 if( atoz[25] ) free(atoz[25]) ;
1369 if( eear ) free(eear) ;
1370 if( HHist ) free(HHist) ;
1371 if( Tvar ) free(Tvar) ;
1372 if( Threshv ) free(Threshv) ;
1373 if( Tvcount ) free(Tvcount) ;
1374
1375 if( Mset != NULL ){
1376 DSET_write(Mset) ; WROTE_DSET(Mset) ; DSET_delete(Mset) ;
1377 }
1378 if( Pset != NULL ){
1379 DSET_write(Pset) ; WROTE_DSET(Pset) ; DSET_delete(Pset) ;
1380 }
1381 if( Zset != NULL ){
1382 DSET_write(Zset) ; WROTE_DSET(Zset) ; DSET_delete(Zset) ;
1383 }
1384 if( Qset != NULL ){
1385 DSET_write(Qset) ; WROTE_DSET(Qset) ; DSET_delete(Qset) ;
1386 }
1387 if( Tset != NULL ){
1388 DSET_write(Tset) ; WROTE_DSET(Tset) ; DSET_delete(Tset) ;
1389 }
1390 if( Tvset != NULL ){
1391 DSET_write(Tvset) ; WROTE_DSET(Tvset) ; DSET_delete(Tvset) ;
1392 }
1393 if( Eset != NULL ){
1394 DSET_write(Eset) ; WROTE_DSET(Eset) ; DSET_delete(Eset) ;
1395 }
1396 if( HHset != NULL ){
1397 DSET_write(HHset) ; WROTE_DSET(HHset) ; DSET_delete(HHset) ;
1398 }
1399 if ( COset != NULL ) {
1400 THD_set_write_compression(COMPRESS_NONE) ;
1401 AFNI_setenv("AFNI_AUTOGZIP NO") ;
1402 if( DSET_TOTALBYTES(COset) > 1000000000 )
1403 ININFO_message("Writing big [%s bytes] -CorrMap output dataset",
1404 approximate_number_string((double)DSET_TOTALBYTES(COset)) ) ;
1405 DSET_write(COset) ; WROTE_DSET(COset) ; DSET_delete(COset) ;
1406 }
1407 INFO_message("total CPU time = %.2f s",COX_cpu_time()) ; exit(0) ;
1408 }
1409