1 #include "mrilib.h"
2
3 #if 0
4 extern void mri_FWHM_1dif_moments( MRI_IMAGE *im , byte *mask ) ;
5 extern void THD_estimate_FWHM_moments_all( THD_3dim_dataset *dset,
6 byte *mask, int demed , int unif ) ;
7 #endif
8
9 #ifdef USE_OMP
10 # include <omp.h>
11 # include "mri_fwhm.c"
12 #else
13 extern void mri_fwhm_mom12_set_stdev_fac(double) ;
14 #endif
15
main(int argc,char * argv[])16 int main( int argc , char *argv[] )
17 {
18 THD_3dim_dataset *inset=NULL ; char *inset_prefix , *cpp ;
19 int iarg=1 , ii , nvals,nvox , ncon ;
20 MRI_IMAGE *outim ; float *outar ;
21 byte *mask=NULL ; int mask_nx=0,mask_ny=0,mask_nz=0 , automask=0 ;
22 char *outfile = NULL ;
23 double fx,fy,fz , cx,cy,cz , ccomb ; int nx,ny,nz , ncomb ;
24 int geom=1 , demed=0 , unif=0 , corder=0 , combine=0 ;
25 char *newprefix=NULL ;
26 int do_acf = -1 ; float acf_rad=0.0f ; int do_classic=0 ; int addcol5=0 ;
27 char *acf_fname="3dFWHMx.1D" ; MRI_IMAGE *acf_im=NULL ; float_quad acf_Epar ;
28 double ct ;
29
30 /*---- for the clueless who wish to become clueful ----*/
31
32 AFNI_SETUP_OMP(0) ;
33
34 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
35 printf(
36 "Usage: 3dFWHMx [options] dataset\n"
37 "\n"
38 "**** NOTICE ****\n"
39 " You should use the '-acf' option (which is what afni_proc.py uses now).\n"
40 " The 'Classic' method giving just a Gaussian FWHM can no longer be\n"
41 " considered reliable for FMRI statistical analyses!\n"
42 "****************\n"
43 "\n"
44 " >>>>> 20 July 2017: Results from the 'Classic' method are no longer output!\n"
45 " >>>>> If you want to see these values, you must give the\n"
46 " >>>>> command line option '-ShowMeClassicFWHM'.\n"
47 " >>>>> You no longer need to give the '-acf' option, as it\n"
48 " >>>>> is now the default method of calculation (and\n"
49 " >>>>> cannot be turned off). Note that if you need the\n"
50 " >>>>> FWHM estimate, the '-acf' method gives a value\n"
51 " >>>>> for that as its fourth output.\n"
52 " >>>>> Options and comments that only apply to the 'Classic' FWHM estimation\n"
53 " >>>>> method are now marked below with this '>>>>>' marker, to indicate that\n"
54 " >>>>> they are obsolete, archaic, and endangered (as well as fattening).\n"
55 "\n"
56 #if 1
57 ">>>>> Unlike the older 3dFWHM, this program computes FWHMs for all sub-bricks\n"
58 ">>>>> in the input dataset, each one separately. The output for each one is\n"
59 ">>>>> written to the file specified by '-out'. The mean (arithmetic or geometric)\n"
60 ">>>>> of all the FWHMs along each axis is written to stdout. (A non-positive\n"
61 ">>>>> output value indicates something bad happened; e.g., FWHM in z is meaningless\n"
62 ">>>>> for a 2D dataset; the estimation method computed incoherent intermediate results.)\n"
63 #endif
64 "\n"
65 "(Classic) METHOD: <<<<< NO LONGER OUTPUT -- SEE ABOVE >>>>>\n"
66 " - Calculate ratio of variance of first differences to data variance.\n"
67 " - Should be the same as 3dFWHM for a 1-brick dataset.\n"
68 " (But the output format is simpler to use in a script.)\n"
69 "\n"
70 "**----------------------------------------------------------------------------**\n"
71 "************* IMPORTANT NOTE [Dec 2015] ****************************************\n"
72 "**----------------------------------------------------------------------------**\n"
73 "A completely new method for estimating and using noise smoothness values is\n"
74 "now available in 3dFWHMx and 3dClustSim. This method is implemented in the\n"
75 "'-acf' options to both programs. 'ACF' stands for (spatial) AutoCorrelation\n"
76 "Function, and it is estimated by calculating moments of differences out to\n"
77 "a larger radius than before.\n"
78 "\n"
79 "Notably, real FMRI data does not actually have a Gaussian-shaped ACF, so the\n"
80 "estimated ACF is then fit (in 3dFWHMx) to a mixed model (Gaussian plus\n"
81 "mono-exponential) of the form\n"
82 " ACF(r) = a * exp(-r*r/(2*b*b)) + (1-a)*exp(-r/c)\n"
83 "where 'r' is the radius, and 'a', 'b', 'c' are the fitted parameters.\n"
84 "The apparent FWHM from this model is usually somewhat larger in real data\n"
85 "than the FWHM estimated from just the nearest-neighbor differences used\n"
86 "in the 'classic' analysis.\n"
87 "\n"
88 "The longer tails provided by the mono-exponential are also significant.\n"
89 "3dClustSim has also been modified to use the ACF model given above to generate\n"
90 "noise random fields.\n"
91 "\n"
92 "**----------------------------------------------------------------------------**\n"
93 "** The take-away (TL;DR or summary) message is that the 'classic' 3dFWHMx and **\n"
94 "** 3dClustSim analysis, using a pure Gaussian ACF, is not very correct for **\n"
95 "** FMRI data -- I cannot speak for PET or MEG data. **\n"
96 #if 0
97 "** **\n"
98 "** You should start using the '-acf' options in your own scripts. AFNI **\n"
99 "** scripts from afni_proc.py are moving away from the 'classic' method. At **\n"
100 "** some point in the future, '-acf' will become the default in both programs, **\n"
101 "** and you will have to actively specify '-classic' to use the older method. **\n"
102 #endif
103 "**----------------------------------------------------------------------------**\n"
104 "\n"
105 "OPTIONS:\n"
106 " -mask mmm = Use only voxels that are nonzero in dataset 'mmm'.\n"
107 " -automask = Compute a mask from THIS dataset, a la 3dAutomask.\n"
108 " [Default = use all voxels]\n"
109 "\n"
110 " -input ddd }=\n"
111 " *OR* }= Use dataset 'ddd' as the input.\n"
112 " -dset ddd }=\n"
113 "\n"
114 " -demed = If the input dataset has more than one sub-brick\n"
115 " (e.g., has a time axis), then subtract the median\n"
116 " of each voxel's time series before processing FWHM.\n"
117 " This will tend to remove intrinsic spatial structure\n"
118 " and leave behind the noise.\n"
119 " [Default = don't do this]\n"
120 " -unif = If the input dataset has more than one sub-brick,\n"
121 " then normalize each voxel's time series to have\n"
122 " the same MAD before processing FWHM. Implies -demed.\n"
123 " [Default = don't do this]\n"
124 " -detrend [q]= Instead of demed (0th order detrending), detrend to\n"
125 " order 'q'. If q is not given, the program picks q=NT/30.\n"
126 " -detrend disables -demed, and includes -unif.\n"
127 " **N.B.: I recommend this option IF you are running 3dFWHMx on\n"
128 " functional MRI time series that have NOT been processed\n"
129 " to remove any activation and/or physiological artifacts.\n"
130 " **** If you are running 3dFWHMx on the residual (errts) time\n"
131 " series from afni_proc.py, you don't need -detrend.\n"
132 " **N.B.: This is the same detrending as done in 3dDespike;\n"
133 " using 2*q+3 basis functions for q > 0.\n"
134 " ******* If you don't use '-detrend', the program checks\n"
135 " if a large number of voxels are have significant\n"
136 " nonzero means. If so, the program will print a warning\n"
137 " message suggesting the use of '-detrend', since inherent\n"
138 " spatial structure in the image will bias the estimation\n"
139 " of the FWHM of the image time series NOISE (which is usually\n"
140 " the point of using 3dFWHMx).\n"
141 " -detprefix d= Save the detrended file into a dataset with prefix 'd'.\n"
142 " Used mostly to figure out what the hell is going on,\n"
143 " when strange results transpire.\n"
144 "\n"
145 ">>>>>\n"
146 " -geom }= If the input dataset has more than one sub-brick,\n"
147 " *OR* }= compute the final estimate as the geometric mean\n"
148 " -arith }= or the arithmetic mean of the individual sub-brick\n"
149 " FWHM estimates. [Default = -geom, for no good reason]\n"
150 "\n"
151 ">>>>>\n"
152 " -combine = combine the final measurements along each axis into\n"
153 " one result\n"
154 "\n"
155 ">>>>>\n"
156 " -out ttt = Write output to file 'ttt' (3 columns of numbers).\n"
157 " If not given, the sub-brick outputs are not written.\n"
158 " Use '-out -' to write to stdout, if desired.\n"
159 "\n"
160 ">>>>>\n"
161 " -compat = Be compatible with the older 3dFWHM, where if a\n"
162 " voxel is in the mask, then its neighbors are used\n"
163 " for differencing, even if they are not themselves in\n"
164 " the mask. This was an error; now, neighbors must also\n"
165 " be in the mask to be used in the differencing.\n"
166 " Use '-compat' to use the older method.\n"
167 " ** NOT RECOMMENDED except for comparison purposes! **\n"
168 "\n"
169 " -ACF [anam] = ** new option Nov 2015 **\n"
170 " *or* The '-ACF' option computes the spatial autocorrelation\n"
171 " -acf [anam] of the data as a function of radius, then fits that\n"
172 " to a model of the form\n"
173 " ACF(r) = a * exp(-r*r/(2*b*b)) + (1-a)*exp(-r/c)\n"
174 " and outputs the 3 model parameters (a,b,c) to stdout.\n"
175 " * The model fit assumes spherical symmetry in the ACF.\n"
176 " * The results shown on stdout are in the format\n"
177 " >>>>> The first 2 lines below will only be output <<<<<\n"
178 " >>>>> if you use the option '-ShowMeClassicFWHM'. <<<<<\n"
179 " >>>>> Otherwise, the 'old-style' FWHM values will <<<<<\n"
180 " >>>>> show up as all zeros (0 0 0 0). <<<<<\n"
181 " # old-style FWHM parameters\n"
182 " 10.4069 10.3441 9.87341 10.2053\n"
183 " # ACF model parameters for a*exp(-r*r/(2*b*b))+(1-a)*exp(-r/c) plus effective FWHM\n"
184 " 0.578615 6.37267 14.402 16.1453\n"
185 " The lines that start with '#' are comments.\n"
186 " >>>>> The first numeric line contains the 'old style' FWHM estimates,\n"
187 " >>>>> FWHM_x FWHM_y FHWM_z FWHM_combined\n"
188 " The second numeric line contains the a,b,c parameters, plus the\n"
189 " combined estimated FWHM from those parameters. In this example,\n"
190 " the fit was about 58%% Gaussian shape, 42%% exponential shape,\n"
191 " and the effective FWHM from this fit was 16.14mm, versus 10.21mm\n"
192 " estimated in the 'old way'.\n"
193 " * If you use '-acf' instead of '-ACF', then the comment #lines\n"
194 " in the stdout information will be omitted. This might help\n"
195 " in parsing the output inside a script.\n"
196 " * The empirical ACF results are also written to the file\n"
197 #ifdef ADD_COL5
198 " 'anam' in 5 columns:\n"
199 " radius ACF(r) model(r) gaussian_NEWmodel(r) gaussian_OLDmodel(r)\n"
200 " where 'gaussian_NEWmodel' is the Gaussian with the FHWM estimated\n"
201 " from the ACF, and 'gaussian_OLDmodel' is with the FWHM estimated\n"
202 " via the 'classic' (Forman 1995) method.\n"
203 #else
204 " 'anam' in 4 columns:\n"
205 " radius ACF(r) model(r) gaussian_NEWmodel(r)(r)\n"
206 " where 'gaussian_NEWmodel' is the Gaussian with the FWHM estimated\n"
207 " from the ACF, NOT via the 'classic' (Forman 1995) method.\n"
208 #endif
209 " * If 'anam' is not given (that is, another option starting\n"
210 " with '-' immediately follows '-acf'), then '3dFWHMx.1D' will\n"
211 " be used for this filename. If 'anam' is set to 'NULL', then\n"
212 " the corresponding output files will not be saved.\n"
213 " * By default, the ACF is computed out to a radius based on\n"
214 " a multiple of the 'classic' FWHM estimate. If you want to\n"
215 " specify that radius (in mm), you can put that value after\n"
216 " the 'anam' parameter, as in '-acf something.1D 40.0'.\n"
217 " * In addition, a graph of these functions will be saved\n"
218 " into file 'anam'.png, for your pleasure and elucidation.\n"
219 " * Note that the ACF calculations are slower than the\n"
220 " 'classic' FWHM calculations.\n"
221 " To reduce this sloth, 3dFWHMx now uses OpenMP to speed things up.\n"
222 #ifndef USE_OMP
223 " (: Unfortunately, this version was NOT compiled to use OpenMP :)\n"
224 #endif
225 " * The ACF modeling is intended to enhance 3dClustSim, and\n"
226 " may or may not be useful for any other purpose!\n"
227 "\n"
228 ">>>>> SAMPLE USAGE: (tcsh)\n"
229 ">>>>> set zork = ( `3dFWHMx -automask -input junque+orig` )\n"
230 ">>>>> Captures the FWHM-x, FWHM-y, FWHM-z values into shell variable 'zork'.\n"
231 "\n"
232 "INPUT FILE RECOMMENDATIONS:\n"
233 "* For FMRI statistical purposes, you DO NOT want the FWHM or ACF to reflect\n"
234 " any spatial structure of the underlying anatomy. Rather, you want\n"
235 " the FWHM/ACF to reflect the spatial structure of the NOISE. This means\n"
236 " that the input dataset should not have anatomical (spatial) structure.\n"
237 "* One good form of input is the output of '3dDeconvolve -errts', which is\n"
238 " the dataset of residuals left over after the GLM fitted signal model is\n"
239 " subtracted out from each voxel's time series.\n"
240 "* If you don't want to go to that much trouble, use '-detrend' to approximately\n"
241 " subtract out the anatomical spatial structure, OR use the output of 3dDetrend\n"
242 " for the same purpose.\n"
243 "* If you do not use '-detrend', the program attempts to find non-zero spatial\n"
244 " structure in the input, and will print a warning message if it is detected.\n"
245 "\n"
246 " *** Do NOT use 3dFWHMx on the statistical results (e.g., '-bucket') from ***\n"
247 " *** 3dDeconvolve or 3dREMLfit!!! The function of 3dFWHMx is to estimate ***\n"
248 " *** the smoothness of the time series NOISE, not of the statistics. This ***\n"
249 " *** proscription is especially true if you plan to use 3dClustSim next!! ***\n"
250 "\n"
251 " *** ------------------- ***\n"
252 " *** NOTE FOR SPM USERS: ***\n"
253 " *** ------------------- ***\n"
254 " *** If you are using SPM for your analyses, and wish to use 3dFHWMX plus ***\n"
255 " *** 3dClustSim for cluster-level thresholds, you need to understand the ***\n"
256 " *** process that AFNI uses. Otherwise, you will likely make some simple ***\n"
257 " *** mistake (such as using 3dFWHMx on the statistical maps from SPM) ***\n"
258 " *** that will render your cluster-level thresholding completely wrong! ***\n"
259 #if 1
260 "\n"
261 ">>>>>\n"
262 "IF YOUR DATA HAS SMOOTH-ISH SPATIAL STRUCTURE YOU CAN'T GET RID OF:\n"
263 "For example, you only have 1 volume, say from PET imaging. In this case,\n"
264 "the standard estimate of the noise smoothness will be mixed in with the\n"
265 "structure of the background. An approximate way to avoid this problem\n"
266 "is provided with the semi-secret '-2difMAD' option, which uses a combination of\n"
267 "first-neighbor and second-neighbor differences to estimate the smoothness,\n"
268 "rather than just first-neighbor differences, and uses the MAD of the differences\n"
269 "rather than the standard deviation. (If you must know the details, read the\n"
270 "source code in mri_fwhm.c!) [For Jatin Vaidya, March 2010]\n"
271 #endif
272 #if 0
273 "\n"
274 "IF YOU WISH TO ALLOW FOR SPATIAL VARIABILITY IN NOISE SMOOTHNESS:\n"
275 "The semi-secret '-1difMOM' option uses moments of the first differences to\n"
276 "crudely estimate a single smoothness that is intended to represent the effect\n"
277 "of variable smoothness -- the idea being that larger smoothness has a bigger\n"
278 "effect on the 3dClustSim results than smaller smoothness does, so the results\n"
279 "from this option will tend to be larger than from the standard option.\n"
280 "** This option is intended for use with single subject data; it is probably\n"
281 " too conservative to use this option on group data when the smoothness\n"
282 " is adjusted upwards based on single subject noise statistics. [August 2015]\n"
283 "** The adjustment upwards involves shifting an intermediate estimate by\n"
284 " a fraction of its standard deviation. The default shift is 1.0 times the\n"
285 " standard deviation estimate. If you want to see the result without this\n"
286 " shift, use '-1difMOM 0.0' (smoothness values will be smaller).\n"
287 #endif
288 "\n"
289 "ALSO SEE:\n"
290 "* The older program 3dFWHM is now completely superseded by 3dFWHMx.\n"
291 "* The program 3dClustSim takes as input the ACF estimates and then\n"
292 " estimates the cluster sizes thresholds to help you get 'corrected'\n"
293 " (for multiple comparisons) p-values.\n"
294 #if 1
295 ">>>>>\n"
296 "* 3dLocalstat -stat FWHM will estimate the FWHM values at each voxel,\n"
297 " using the same first-difference algorithm as this program, but applied\n"
298 " only to a local neighborhood of each voxel in turn.\n"
299 #endif
300 "* 3dLocalACF will estimate the 3 ACF parameters in a local neighborhood\n"
301 " around each voxel.\n"
302 ">>>>>\n"
303 "* 3dBlurToFWHM will iteratively blur a dataset (inside a mask) to have\n"
304 " a given global FWHM. This program may or may not be useful :)\n"
305 "* 3dBlurInMask will blur a dataset inside a mask, but doesn't measure FWHM or ACF.\n"
306 "\n"
307 "-- Zhark, Ruler of the (Galactic) Cluster!\n"
308 ) ;
309 PRINT_AFNI_OMP_USAGE("3dFWHMx",NULL) ; PRINT_COMPILE_DATE ; exit(0) ;
310 }
311
312 /*---- official startup ---*/
313
314 PRINT_VERSION("3dFWHMx"); mainENTRY("3dFWHMx main"); machdep(); AUTHOR("The Bob");
315 AFNI_logger("3dFWHMx",argc,argv) ;
316
317 /*---- loop over options ----*/
318
319 while( iarg < argc && argv[iarg][0] == '-' ){
320
321 if( strcmp(argv[iarg],"-detrend") == 0 ){ /* 10 May 2007 */
322 corder = -1 ;
323 if( iarg < argc-1 && isdigit(argv[iarg+1][0]) ){
324 corder = (int)strtod(argv[++iarg],NULL) ;
325 if( corder == 0 ){
326 demed = 1 ; INFO_message("-detrend 0 replaced by -demed") ;
327 }
328 }
329 iarg++ ; continue ;
330 }
331
332 if( strcmp(argv[iarg],"-2dif") == 0 ){ /* 20 Nov 2006 */
333 mri_fwhm_setfester( mri_estimate_FWHM_12dif ) ; /* secret option */
334 iarg++ ; continue ;
335 }
336
337 if( strcmp(argv[iarg],"-2difMAD") == 0 ){ /* 24 Mar 2010 */
338 mri_fwhm_setfester( mri_estimate_FWHM_12dif_MAD ) ; /* secret option */
339 iarg++ ; continue ;
340 }
341
342 if( strcmp(argv[iarg],"-1difMOM") == 0 ){ /* 11 Aug 2015 */
343 mri_fwhm_setfester( mri_FWHM_1dif_mom12 ) ;
344 iarg++ ;
345 if( iarg < argc && argv[iarg][0] >= '0' && argv[iarg][0] <= '9' ){
346 double val = strtod(argv[iarg++],NULL) ;
347 if( val >= 0.0 && val <= 1.0 ) mri_fwhm_mom12_set_stdev_fac(val) ;
348 }
349 continue ;
350 }
351
352 if( strncmp(argv[iarg],"-geom",4) == 0 ){ /* 15 Nov 2006 */
353 geom = 1 ; iarg++ ; continue ;
354 }
355 if( strncmp(argv[iarg],"-arith",5) == 0 ){ /* 15 Nov 2006 */
356 geom = 0 ; iarg++ ; continue ;
357 }
358 if( strncmp(argv[iarg],"-demed",5) == 0 ){ /* 15 Nov 2006 */
359 demed = 1 ; iarg++ ; continue ;
360 }
361 if( strncmp(argv[iarg],"-unif",5) == 0 ){ /* 07 Dec 2006 */
362 unif = demed = 1 ; iarg++ ; continue ;
363 }
364 if( strncmp(argv[iarg],"-comb",4) == 0 ){ /* 24 Mar 2010 */
365 combine = 1 ; iarg++ ; continue ;
366 }
367
368 if( strncmp(argv[iarg],"-compat",6) == 0 ){ /* 09 Nov 2006 */
369 FHWM_1dif_dontcheckplus(1) ; iarg++ ; continue ;
370 }
371
372 if( strncasecmp(argv[iarg],"-ACF",4) == 0 ){ /* 09 Nov 2015 */
373 do_acf = 1 ; if( argv[iarg][1] == 'a' ) do_acf = -1 ;
374 iarg++ ;
375 if( iarg < argc && argv[iarg][0] != '-' ){
376 acf_fname = strdup(argv[iarg]) ;
377 if( !THD_filename_ok(acf_fname) )
378 ERROR_exit("filename after -ACF is invalid!") ;
379 iarg++ ;
380 if( iarg < argc && isdigit(argv[iarg][0]) ){ /* 07 Dec 2015 */
381 acf_rad = (float)strtod(argv[iarg],NULL) ;
382 iarg++ ;
383 }
384 }
385 continue ;
386 }
387
388 if( strcasecmp(argv[iarg],"-addcol5") == 0 ){
389 addcol5++ ; iarg++ ; continue ;
390 }
391
392 if( strcasecmp(argv[iarg],"-ShowMeClassicFWHM") == 0 ){ /* 20 Jul 2017 */
393 do_classic = 1 ; iarg++ ;
394 WARNING_message("Using the 'Classic' Gaussian FWHM is not recommended :(") ;
395 ININFO_message (" The '-acf' method gives a FWHM estimate which is more robust;") ;
396 ININFO_message (" however, assuming the spatial correlation of FMRI noise has") ;
397 ININFO_message (" a Gaussian shape is not a good model.") ;
398 continue ;
399 }
400
401 if( strncmp(argv[iarg],"-out",4) == 0 ){
402 if( ++iarg >= argc ) ERROR_exit("Need argument after '-out'") ;
403 outfile = argv[iarg] ;
404 if( strcasecmp(outfile,"NULL") == 0 ) outfile = NULL ;
405 else if( !THD_filename_ok(outfile) ) ERROR_exit("Illegal filename after '-out'") ;
406 iarg++ ; continue ;
407 }
408
409 if( strcmp(argv[iarg],"-input") == 0 || strcmp(argv[iarg],"-dset") == 0 ){
410 if( inset != NULL ) ERROR_exit("Can't have two -input options") ;
411 if( ++iarg >= argc ) ERROR_exit("Need argument after '-input'") ;
412 inset = THD_open_dataset( argv[iarg] ) ;
413 CHECK_OPEN_ERROR(inset,argv[iarg]) ;
414 iarg++ ; continue ;
415 }
416
417 if( strcmp(argv[iarg],"-mask") == 0 ){
418 THD_3dim_dataset *mset ; int mmm ;
419 if( ++iarg >= argc ) ERROR_exit("Need argument after '-mask'") ;
420 if( mask != NULL || automask ) ERROR_exit("Can't have two mask inputs") ;
421 mset = THD_open_dataset( argv[iarg] ) ;
422 CHECK_OPEN_ERROR(mset,argv[iarg]) ;
423 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
424 mask_nx = DSET_NX(mset); mask_ny = DSET_NY(mset); mask_nz = DSET_NZ(mset);
425 mask = THD_makemask( mset , 0 , 0.5f, 0.0f ) ; DSET_delete(mset) ;
426 if( mask == NULL ) ERROR_exit("Can't make mask from dataset '%s'",argv[iarg]) ;
427 mmm = THD_countmask( mask_nx*mask_ny*mask_nz , mask ) ;
428 INFO_message("Number of voxels in mask = %d",mmm) ;
429 if( mmm < 16 ) ERROR_exit("Mask is too small to process") ;
430 iarg++ ; continue ;
431 }
432
433 if( strcmp(argv[iarg],"-automask") == 0 ){
434 if( mask != NULL ) ERROR_exit("Can't have -automask and -mask") ;
435 automask = 1 ;
436 iarg++ ; continue ;
437 }
438
439 if( strcmp(argv[iarg],"-detprefix") == 0 ){
440 newprefix = argv[++iarg] ; iarg++ ; continue ;
441 }
442
443 ERROR_exit("Unknown option '%s'",argv[iarg]) ;
444
445 } /*--- end of loop over options ---*/
446
447 /*-------------- deal with input dataset --------------*/
448
449 if( inset == NULL ){
450 if( iarg >= argc ) ERROR_exit("No input dataset on command line?") ;
451 inset = THD_open_dataset( argv[iarg] ) ;
452 CHECK_OPEN_ERROR(inset,argv[iarg]) ;
453 }
454
455 inset_prefix = strdup( DSET_PREFIX(inset) ) ;
456 cpp = strstr(inset_prefix,".nii") ; if( cpp != NULL ) *cpp = '\0' ;
457 cpp = strstr(inset_prefix,"+orig") ; if( cpp != NULL ) *cpp = '\0' ;
458 cpp = strstr(inset_prefix,"+tlrc") ; if( cpp != NULL ) *cpp = '\0' ;
459 cpp = THD_trailname(inset_prefix,0) ; if( cpp != NULL ) inset_prefix = cpp ;
460
461 { THD_3dim_dataset *qset = THD_remove_allzero(inset) ; /* 25 Jul 2017 */
462 if( qset != NULL ){
463 DSET_delete(inset) ; inset = qset ;
464 }
465 }
466
467 if( (demed || unif || corder ) && DSET_NVALS(inset) < 4 ){
468 WARNING_message(
469 "-demed and/or -corder and/or -unif ignored: only %d input sub-bricks",
470 DSET_NVALS(inset) ) ;
471 demed = corder = unif = 0 ;
472 }
473
474 if( demed && corder ){
475 demed = 0 ; WARNING_message("-demed is overriden by -corder") ;
476 }
477
478 if( corder < 0 ){
479 corder = DSET_NVALS(inset) / 30 ;
480 if( corder == 0 ){
481 WARNING_message("Fewer than 30 time points ==> -corder converted to -unif") ;
482 unif = demed = 1 ;
483 }
484 } else if( corder > 0 && 2*corder+3 >= DSET_NVALS(inset) ){
485 ERROR_exit("-corder %d is too big for this dataset",corder) ;
486 }
487
488 DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;
489
490 nvals = DSET_NVALS(inset) ; nvox = DSET_NVOX(inset) ;
491
492 if( mask != NULL ){
493 if( mask_nx != DSET_NX(inset) ||
494 mask_ny != DSET_NY(inset) ||
495 mask_nz != DSET_NZ(inset) )
496 ERROR_exit("-mask dataset grid dimensions don't match input dataset") ;
497
498 } else if( automask ){
499 int mmm ;
500 mask = THD_automask( inset ) ;
501 if( mask == NULL )
502 ERROR_message("Can't create -automask from input dataset?") ;
503 mmm = THD_countmask( DSET_NVOX(inset) , mask ) ;
504 INFO_message("Number of voxels in automask = %d",mmm) ;
505 if( mmm < 16 ) ERROR_exit("Automask is too small to process") ;
506 } else {
507 mask = (byte *)malloc(sizeof(byte)*nvox) ;
508 memset(mask,1,sizeof(byte)*nvox) ;
509 }
510 if( nvals > 3 ){
511 for( ncon=ii=0 ; ii < nvox ; ii++ ){
512 if( mask[ii] && THD_voxel_is_constant(ii,inset) ){ mask[ii] = 0; ncon++; }
513 }
514 if( ncon > 0 )
515 WARNING_message("removed %d voxels from mask because they are constant in time",ncon) ;
516 }
517
518 /*-- if NOT detrending or de-median-ing, check if that's a good idea --*/
519
520 if( !(corder > 0 || demed) && nvals > 4 ){ /* 13 Aug 2010 */
521 MRI_IMARR *imar ;
522 imar = THD_medmad_bricks(inset) ;
523 if( imar != NULL ){
524 float *med , *mad ; int nchk=0,nbad=0 ;
525 med = MRI_FLOAT_PTR(IMARR_SUBIM(imar,0)) ;
526 mad = MRI_FLOAT_PTR(IMARR_SUBIM(imar,1)) ;
527 for( ii=0 ; ii < nvox ; ii++ ){
528 if( mask[ii] && mad[ii] > 0.0f ){
529 nchk++ ; if( fabsf(med[ii]) > 6.66f*mad[ii] ) nbad++ ;
530 }
531 }
532 DESTROY_IMARR(imar) ;
533 if( nbad > nchk/16 ){
534 WARNING_message("Suggestion: use the '-detrend' option:") ;
535 ININFO_message ("%d (out of %d) voxel time series have significant means",nbad,nchk) ;
536 }
537 }
538 }
539
540 /*-- if detrending, do that now --*/
541
542 ct = COX_cpu_time() ;
543 if( corder > 0 ){
544 int nref=2*corder+3 , jj,iv,kk ;
545 float **ref , tm,fac,fq ;
546 THD_3dim_dataset *newset ;
547
548 INFO_message("detrending start: %d baseline funcs, %d time points",nref,nvals) ;
549
550 ref = THD_build_trigref( corder , nvals ) ;
551 if( ref == NULL ) ERROR_exit("THD_build_trigref failed!") ;
552
553 newset = THD_detrend_dataset( inset , nref , ref , 2 , 1 , mask , NULL ) ;
554 if( newset == NULL ) ERROR_exit("detrending failed!") ;
555
556 for(jj=0;jj<nref;jj++) free(ref[jj]) ;
557 free(ref); DSET_delete(inset); inset=newset;
558 demed = 0 ; unif = 1 ;
559 ININFO_message("detrending done (%.2f CPU s thus far)",COX_cpu_time()-ct) ;
560
561 if( newprefix != NULL ){ /** for debugging **/
562 EDIT_dset_items(newset,ADN_prefix,newprefix,NULL) ;
563 (void)THD_deconflict_prefix(newset) ;
564 DSET_write(newset) ; WROTE_DSET(newset) ;
565 }
566 }
567
568 /*-- do the FWHM-izing work --*/
569
570 if( do_classic )
571 INFO_message("start Classic FWHM calculations") ;
572
573 outim = THD_estimate_FWHM_all( inset , mask , demed,unif ) ;
574
575 #if 0
576 if( AFNI_yesenv("MOMENTS") )
577 THD_estimate_FWHM_moments_all( inset , mask , demed,unif ) ;
578 #endif
579
580 if( !do_acf ) DSET_unload(inset) ;
581
582 if( outim == NULL ) ERROR_exit("Function THD_estimate_FWHM_all() fails?!") ;
583
584 if( outfile != NULL ) mri_write_ascii( outfile , outim ) ;
585
586 outar = MRI_FLOAT_PTR(outim) ;
587
588 nx = thd_floatscan( 3*nvals, outar ) ; /* 07 Dec 2006 */
589 if( nx > 0 ) WARNING_message("found %d non-finite FWHM array values!",nx);
590
591 nx = ny = nz = 0 ;
592 if( geom ){
593 cx = cy = cz = 0.0 ;
594 for( ii=0 ; ii < nvals ; ii++ ){
595 fx = outar[0+3*ii]; fy = outar[1+3*ii]; fz = outar[2+3*ii];
596 /* INFO_message("geom: fx=%g fy=%g fz=%g",fx,fy,fz) ; */
597 if( fx > 0.0 ){ cx += log(fx) ; nx++ ; }
598 if( fy > 0.0 ){ cy += log(fy) ; ny++ ; }
599 if( fz > 0.0 ){ cz += log(fz) ; nz++ ; }
600 }
601 cx = (nx == 0) ? 0.0 : exp(cx/nx) ;
602 cy = (ny == 0) ? 0.0 : exp(cy/ny) ;
603 cz = (nz == 0) ? 0.0 : exp(cz/nz) ;
604 ccomb = 1.0 ; ncomb = 0 ;
605 if( cx > 0.0 ){ ccomb *= cx ; ncomb++ ; }
606 if( cy > 0.0 ){ ccomb *= cy ; ncomb++ ; }
607 if( cz > 0.0 ){ ccomb *= cz ; ncomb++ ; }
608 if( ncomb == 2 ) ccomb = sqrt(ccomb) ;
609 else if( ncomb == 3 ) ccomb = cbrt(ccomb) ;
610 } else {
611 cx = cy = cz = 0.0 ;
612 for( ii=0 ; ii < nvals ; ii++ ){
613 fx = outar[0+3*ii]; fy = outar[1+3*ii]; fz = outar[2+3*ii];
614 if( fx > 0.0 ){ cx += fx ; nx++ ; }
615 if( fy > 0.0 ){ cy += fy ; ny++ ; }
616 if( fz > 0.0 ){ cz += fz ; nz++ ; }
617 }
618 cx = (nx == 0) ? 0.0 : cx/nx ;
619 cy = (ny == 0) ? 0.0 : cy/ny ;
620 cz = (nz == 0) ? 0.0 : cz/nz ;
621 /* fix arithmetic mean 19 Jul 2010 [rickr] */
622 ccomb = 0.0 ; ncomb = 0 ;
623 if( cx > 0.0 ){ ccomb += cx ; ncomb++ ; }
624 if( cy > 0.0 ){ ccomb += cy ; ncomb++ ; }
625 if( cz > 0.0 ){ ccomb += cz ; ncomb++ ; }
626 if( ncomb > 1 ) ccomb /= ncomb ;
627 }
628
629 if( do_classic )
630 ININFO_message("Classic FWHM done (%.2f CPU s thus far)",COX_cpu_time()-ct) ;
631
632 if( do_acf ){
633 MCW_cluster *acf ; int pp ; float dx,dy,dz,arr ;
634
635 if( acf_rad <= 0.0f ) acf_rad = 2.999f * ccomb ;
636 dx = fabsf(DSET_DX(inset)); dy = fabsf(DSET_DY(inset)); dz = fabsf(DSET_DZ(inset));
637 arr = 3.999f * cbrtf(dx*dy*dz) ; if( acf_rad < arr ) acf_rad = arr ;
638 INFO_message("start ACF calculations out to radius = %.2f mm",acf_rad) ;
639
640 acf = THD_estimate_ACF( inset , mask , demed,unif , acf_rad ) ;
641 if( acf == NULL ) ERROR_exit("Error calculating ACF :-(") ;
642 #if 0
643 printf("# ACF dx=%g dy=%g dz=%g\n", dx,dy,dz ) ;
644 printf("dx dy dz ACF\n") ;
645 printf("--- --- --- ------\n") ;
646 for( pp=0 ; pp < acf->num_pt ; pp++ )
647 printf("%3d %3d %3d %.5f\n",
648 acf->i[pp] , acf->j[pp] , acf->k[pp] , acf->mag[pp] ) ;
649 #endif
650
651 acf_Epar = ACF_cluster_to_modelE( acf, fabsf(DSET_DX(inset)) ,
652 fabsf(DSET_DY(inset)) ,
653 fabsf(DSET_DZ(inset)) ) ;
654
655 if( acf_fname != NULL && strcmp(acf_fname,"NULL") != 0 ) acf_im = ACF_get_1D() ;
656
657 ININFO_message("ACF done (%.2f CPU s thus far)",COX_cpu_time()-ct) ;
658
659 if( do_acf > 0 )
660 printf("# old-style FWHM parameters\n") ;
661 if( do_classic )
662 printf(" %g %g %g %g",cx,cy,cz,ccomb) ;
663 else
664 printf(" 0 0 0 0") ;
665 printf("\n") ;
666
667 if( do_acf > 0 )
668 printf("# ACF model parameters for a*exp(-r*r/(2*b*b))+(1-a)*exp(-r/c) plus effective FWHM\n") ;
669 printf(" %g %g %g %g\n",acf_Epar.a,acf_Epar.b,acf_Epar.c,acf_Epar.d) ;
670
671 if( acf_im != NULL ){
672 char cmd[4096] ;
673
674 if( addcol5 ){
675 MRI_IMAGE *qim,*pim ; float *rar, *qar, sig ; MRI_IMARR *imar ;
676 qim = mri_new( acf_im->nx , 1 , MRI_float ) ;
677 qar = MRI_FLOAT_PTR(qim) ; rar = MRI_FLOAT_PTR(acf_im) ;
678 sig = FWHM_TO_SIGMA(ccomb) ;
679 for( pp=0 ; pp < acf_im->nx ; pp++ )
680 qar[pp] = exp(-0.5*rar[pp]*rar[pp]/(sig*sig)) ;
681 INIT_IMARR(imar) ; ADDTO_IMARR(imar,acf_im) ; ADDTO_IMARR(imar,qim) ;
682 pim = mri_catvol_1D(imar,2) ; DESTROY_IMARR(imar) ; acf_im = pim ;
683 }
684
685 mri_write_1D( acf_fname , acf_im ) ;
686
687 if( addcol5 ){
688 INFO_message("ACF 1D file [radius ACF mixed_model gaussian_NEWmodel gaussian_OLDmodel] written to %s",acf_fname) ;
689 sprintf(cmd,
690 "1dplot -one -xlabel 'r (mm)'"
691 " -ylabel 'Autocorrelation \\small [FWHM=\\green %.2f \\blue %.2f\\black]'"
692 " -yaxis 0:1:10:2 -DAFNI_1DPLOT_BOXSIZE=0.004"
693 " -plabel '\\small\\noesc %s\\esc\\red %.2f*exp[-r^2/2*%.2f^2]+%.2f*exp[-r/%.2f]'"
694 " -box -png %s.png -x %s'[0]' %s'[1]' %s'[2]' %s'[3]' %s'[4]'" ,
695 acf_Epar.d , ccomb ,
696 inset_prefix ,
697 acf_Epar.a , acf_Epar.b , 1.0f-acf_Epar.a , acf_Epar.c ,
698 acf_fname, acf_fname, acf_fname, acf_fname, acf_fname, acf_fname ) ;
699 } else {
700 INFO_message("ACF 1D file [radius ACF mixed_model gaussian_NEWmodel] written to %s",acf_fname) ;
701 sprintf(cmd,
702 "1dplot -one -xlabel 'r (mm)'"
703 " -ylabel 'Autocorrelation \\small [FWHM=\\green %.2f\\black]'"
704 " -yaxis 0:1:10:2 -DAFNI_1DPLOT_BOXSIZE=0.004"
705 " -plabel '\\small\\noesc %s\\esc\\red %.2f*exp[-r^2/2*%.2f^2]+%.2f*exp[-r/%.2f]'"
706 " -box -png %s.png -x %s'[0]' %s'[1]' %s'[2]' %s'[3]'" ,
707 acf_Epar.d ,
708 inset_prefix ,
709 acf_Epar.a , acf_Epar.b , 1.0f-acf_Epar.a , acf_Epar.c ,
710 acf_fname, acf_fname, acf_fname, acf_fname, acf_fname ) ;
711 }
712
713 system(cmd) ;
714 ININFO_message("and 1dplot-ed to file %s.png",acf_fname) ;
715 }
716
717 } else if( do_classic ){ /* no ACF -- the OLD way */
718
719 printf(" %g %g %g",cx,cy,cz) ;
720 if( combine ) printf(" %g",ccomb) ;
721 printf("\n") ;
722 }
723 exit(0) ;
724 }
725