1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*---------------------------------------------------------------------------*/
8 /*
9 This program performs cluster detection in 3D datasets.
10 */
11
12 /*---------------------------------------------------------------------------*/
13
14 #define PROGRAM_NAME "3dclust" /* name of this program */
15 #define PROGRAM_AUTHOR "RW Cox et alii" /* program author */
16 #define PROGRAM_DATE "12 Jul 2017" /* date of last program mod */
17
18 /*---------------------------------------------------------------------------*/
19
20 /* Modified 3/26/99 by BDW to enable -1erode and -1dilate options. */
21 /* Modified 1/19/99 by BDW to allow use of signed intensities in calculation
22 of cluster averages, etc. (-noabs option), as requested by H. Garavan */
23 /* Modified 1/24/97 by BDW to combine the previous modifications */
24 /* Modified 12/27/96 by BDW Corrections to the SEM calculations. */
25 /* Modified 4/19/96 by MSB to give cluster intensity and extent information */
26 /* Modified 11/1/96 by MSB to give cluster intensity standard error of the mean
27 (SEM) and average intensity and SEM for all voxels (globally)
28
29 Modified source files were
30 3dclust.c
31 Modified variables were:
32 sem (added) contains calculated SEM
33 sqsum (added) contains cumulative sum of squares
34 glmm (added) global mean
35 glsqsum (added) global sum of squares
36 Major modifications to code are noted with "MSB 11/1/96" and comments
37 Testing and Verification
38 Program was run on data file JKt3iravfm0.5cymc+orig
39 3dclust -1noneg -1thresh 0.50 15 100 gives a cluster of
40 138 15525.0 -32.8 46.2 -26.7 -16.9 39.4 10.0 338.7 12.4 984.0
41 Voxel values in this cluster were downloaded to Excel;
42 the avg, SE, and max were calculated and identical values were found.
43
44 information and summary information
45 Important User Notes:
46 - center of mass calculations makes use of the intensity at each point; this
47 should perhaps be selectable
48 - cluster calculations are all done on the absolute value of the intensity;
49 hence, positive and negative voxels can be grouped together into the same
50 cluster, which skews results; to prevent this, use the -1noneg option
51 - SEM values are not realistic for interpolated data sets (because comparisons
52 are not independent) a ROUGH correction is to multiply the interpolated SEM
53 of the interpolated data set by the square root of the number of interpolated
54 voxels per original voxel
55
56 */
57
58 /*-- 29 Nov 2001: RWCox adds the -prefix option --*/
59 /*-- 30 Apr 2002: RWCox adds the -mni option --*/
60 /*-- 21 Jul 2005: P Christidis modified -help menu --*/
61 /*-- 26 Aug 2011: RWCox adds the -savemask option --*/
62
63 /*---------------------------------------------------------------------------*/
64 #include <stdio.h>
65 #include <stdlib.h>
66
67 #include "mrilib.h"
68
69 static int do_NN = 0 ; /* 12 Jul 2017 */
70 static int NNvox = 0 ; /* 08 Sep 2017 */
71
72 static EDIT_options CL_edopt ;
73 static int CL_ivfim=-1 , CL_ivthr=-1 ;
74
75 static int CL_nopt ;
76
77 static int CL_summarize = 0 ;
78
79 static int CL_noabs = 0; /* BDW 19 Jan 1999 */
80
81 static int CL_verbose = 0 ; /* RWC 01 Nov 1999 */
82
83 static int CL_quiet = 0; /* MSB 02 Dec 1999 */
84
85 static char * CL_prefix = NULL ; /* 29 Nov 2001 -- RWCox */
86 static char * CL_savemask = NULL ; /* 26 Aug 2011 */
87 static int do_binary = 0 ; /* 26 May 2017 */
88
89 static int CL_do_mni = 0 ; /* 30 Apr 2002 -- RWCox */
90
91 static int CL_1Dform = 1 ; /* 02 Mar 2006 -- Zaid (it's hopeless)
92 Changed to '1' 23 Mar 2007 -- Said (still hopeless)*/
93
94 static int no_inmask = 1 ; /* 02 Aug 2011 */
95
96 /**-- RWCox: July 1997
97 Report directions based on AFNI_ORIENT environment --**/
98
99 static THD_coorder CL_cord ;
100
101 int compare_cluster( void * , void * ) ;
102 void CL_read_opts( int , char ** ) ;
103 #define CL_syntax(str) ERROR_exit("%s",(str))
104
105 void MCW_fc7( float qval , char * buf ) ;
106
107 /*---------------------------------------------------------------------------*/
108
main(int argc,char * argv[])109 int main( int argc , char * argv[] )
110 {
111 float rmm , vmul ;
112 int iarg ;
113 THD_3dim_dataset *dset=NULL ;
114 void * vfim ;
115 int nx,ny,nz , nxy,nxyz , ivfim ,
116 iclu , ptmin , ipt , ii,jj,kk , ndet , nopt ;
117 float dx,dy,dz , xx,yy,zz,mm ,ms, fimfac,
118 xxmax,yymax,zzmax, mmmax , msmax ,
119 RLmax, RLmin, APmax, APmin, ISmax, ISmin;
120 double xxsum,yysum,zzsum,mmsum , volsum , mssum ;
121 double mean, sem, sqsum, glmmsum, glsqsum, glmssum,
122 glmean, glxxsum, glyysum, glzzsum;
123 MCW_cluster_array * clar , * clbig ;
124 MCW_cluster * cl ;
125 THD_fvec3 fv ;
126 int nvox_total ;
127 float vol_total ;
128 char buf1[16],buf2[16],buf3[16] ;
129 float dxf,dyf,dzf ; /* 24 Jan 2001: for -dxyz=1 option */
130 int do_mni ; /* 30 Apr 2002 */
131 char c1d[2] = {""}, c1dn[2] = {""};
132 byte *mask=NULL ; int nmask=0 ; /* 02 Aug 2011 */
133
134 mainENTRY("3dclust"); machdep();
135
136 if( argc < 4 || strncmp(argv[1],"-help",4) == 0 ){
137 printf ("\n\n");
138 printf ("Program: %s \n", PROGRAM_NAME);
139 printf ("Author: %s \n", PROGRAM_AUTHOR);
140 printf ("Date: %s \n", PROGRAM_DATE);
141 printf ("\n");
142 printf(
143 "3dclust - performs simple-minded cluster detection in 3D datasets \n"
144 "\n"
145 " *** PLEASE NOTE THAT THE NEWER PROGRAM 3dClusterize ***\n"
146 " *** IS BETTER AND YOU SHOULD USE THAT FROM NOW ON!! ***\n"
147 " \n"
148 " This program can be used to find clusters of 'active' voxels and \n"
149 " print out a report about them. \n"
150 " * 'Active' refers to nonzero voxels that survive the threshold \n"
151 " that you (the user) have specified \n"
152 " * Clusters are defined by a connectivity radius parameter 'rmm' \n"
153 " *OR*\n"
154 " Clusters are defined by how close neighboring voxels must\n"
155 " be in the 3D grid:\n"
156 " first nearest neighbors (-NN1)\n"
157 " second nearest neighbors (-NN2)\n"
158 " third nearest neighbors (-NN3)\n"
159 " \n"
160 " Note: by default, this program clusters on the absolute values \n"
161 " of the voxels \n"
162 "----------------------------------------------------------------------- \n"
163 "Usage:\n"
164 " \n"
165 " 3dclust [editing options] [other options] rmm vmul dset ... \n"
166 " \n"
167 " *OR*\n"
168 " \n"
169 " 3dclust [editing options] -NNx dset ...\n"
170 " where '-NNx' is one of '-NN1' or '-NN2' or '-NN3':\n"
171 " -NN1 == 1st nearest-neighbor (faces touching) clustering\n"
172 " -NN2 == 2nd nearest-neighbor (edges touching) clustering\n"
173 " -NN2 == 3rd nearest-neighbor (corners touching) clustering\n"
174 " Optionally, you can put an integer after the '-NNx' option, to\n"
175 " indicate the minimum number of voxels to allow in a cluster;\n"
176 " for example: -NN2 60\n"
177 "----------------------------------------------------------------------- \n"
178 "Examples: \n"
179 "--------- \n"
180 " \n"
181 " 3dclust -1clip 0.3 5 2000 func+orig'[1]' \n"
182 " 3dclust -1noneg -1thresh 0.3 5 2000 func+orig'[1]' \n"
183 " 3dclust -1noneg -1thresh 0.3 5 2000 func+orig'[1]' func+orig'[3] \n"
184 " \n"
185 " 3dclust -noabs -1clip 0.5 -dxyz=1 1 10 func+orig'[1]' \n"
186 " 3dclust -noabs -1clip 0.5 5 700 func+orig'[1]' \n"
187 " \n"
188 " 3dclust -noabs -2clip 0 999 -dxyz=1 1 10 func+orig'[1]' \n"
189 " \n"
190 " 3dclust -1clip 0.3 5 3000 func+orig'[1]' \n"
191 " 3dclust -quiet -1clip 0.3 5 3000 func+orig'[1]' \n"
192 " 3dclust -summarize -quiet -1clip 0.3 5 3000 func+orig'[1]' \n"
193 " 3dclust -1Dformat -1clip 0.3 5 3000 func+orig'[1]' > out.1D\n"
194 "----------------------------------------------------------------------- \n"
195 " \n"
196 "Arguments (must be included on command line): \n"
197 "--------- \n"
198 " \n"
199 "THE OLD WAY TO SPECIFY THE TYPE OF CLUSTERING\n"
200 "\n"
201 " rmm : cluster connection radius (in millimeters). \n"
202 " All nonzero voxels closer than rmm millimeters \n"
203 " (center-to-center distance) to the given voxel are \n"
204 " included in the cluster. \n"
205 " * If rmm = 0, then clusters are defined by nearest-\n"
206 " neighbor connectivity \n"
207 " \n"
208 " vmul : minimum cluster volume (micro-liters) \n"
209 " i.e., determines the size of the volume cluster. \n"
210 " * If vmul = 0, then all clusters are kept. \n"
211 " * If vmul < 0, then the absolute vmul is the minimum\n"
212 " number of voxels allowed in a cluster. \n"
213 "\n"
214 " If you do not use one of the '-NNx' options, you must give the\n"
215 " numbers for rmm and vmul just before the input dataset name(s)\n"
216 "\n"
217 "THE NEW WAY TO SPECIFY TYPE OF CLUSTERING [13 Jul 2017]\n"
218 "\n"
219 " -NN1 or -NN2 or -NN3\n"
220 "\n"
221 " If you use one of these '-NNx' options, you do NOT give the rmm\n"
222 " and vmul values. Instead, after all the options that start with '-',\n"
223 " you just give the input dataset name(s).\n"
224 " If you want to set a minimum cluster size using '-NNx', put the minimum\n"
225 " voxel count immediately after, as in '-NN3 100'.\n"
226 "\n"
227 "FOLLOWED BY ONE (or more) DATASETS\n"
228 " \n"
229 " dset : input dataset (more than one allowed, but only the \n"
230 " first sub-brick of the dataset) \n"
231 " \n"
232 " The results are sent to standard output (i.e., the screen): \n"
233 " if you want to save them in a file, then use redirection, as in\n"
234 "\n"
235 " 3dclust -1thresh 0.4 -NN2 Elvis.nii'[1]' > Elvis.clust.txt\n"
236 " \n"
237 "----------------------------------------------------------------------- \n"
238 " \n"
239 "Options: \n"
240 "------- \n"
241 " \n"
242 " Editing options are as in 3dmerge (see 3dmerge -help) \n"
243 " (including -1thresh, -1dindex, -1tindex, -dxyz=1 options) \n"
244 "\n"
245 " -NN1 => described earlier;\n"
246 " -NN2 => replaces the use of 'rmm' to specify the\n"
247 " -NN3 => clustering method (vmul is set to 2 voxels)\n"
248 " \n"
249 " -noabs => Use the signed voxel intensities (not the absolute \n"
250 " value) for calculation of the mean and Standard \n"
251 " Error of the Mean (SEM) \n"
252 " \n"
253 " -summarize => Write out only the total nonzero voxel \n"
254 " count and volume for each dataset \n"
255 " \n"
256 " -nosum => Suppress printout of the totals \n"
257 " \n"
258 " -verb => Print out a progress report (to stderr) \n"
259 " as the computations proceed \n"
260 " \n"
261 " -1Dformat => Write output in 1D format (now default). You can \n"
262 " redirect the output to a .1D file and use the file \n"
263 " as input to whereami for obtaining Atlas-based \n"
264 " information on cluster locations. \n"
265 " See whereami -help for more info. \n"
266 " -no_1Dformat => Do not write output in 1D format. \n"
267 " \n"
268 " -quiet => Suppress all non-essential output \n"
269 " \n"
270 " -mni => If the input dataset has the +tlrc view, this option \n"
271 " will transform the output xyz-coordinates from TLRC to \n"
272 " MNI space.\n"
273 " \n"
274 " N.B.0: Only use this option if the dataset is in Talairach \n"
275 " space, NOT when it is already in MNI space. \n"
276 " N.B.1: The MNI template brain is about 5 mm higher (in S), \n"
277 " 10 mm lower (in I), 5 mm longer (in PA), and tilted \n"
278 " about 3 degrees backwards, relative to the Talairach- \n"
279 " Tournoux Atlas brain. For more details, see, e.g.: \n"
280 " https://imaging.mrc-cbu.cam.ac.uk/imaging/MniTalairach\n"
281 " N.B.2: If the input dataset does not have the +tlrc view, \n"
282 " then the only effect is to flip the output coordinates\n"
283 " to the 'LPI' (neuroscience) orientation, as if you \n"
284 " gave the '-orient LPI' option.) \n"
285 " \n"
286 " -isovalue => Clusters will be formed only from contiguous (in the \n"
287 " rmm sense) voxels that also have the same value. \n"
288 " \n"
289 " N.B.: The normal method is to cluster all contiguous \n"
290 " nonzero voxels together. \n"
291 " \n"
292 " -isomerge => Clusters will be formed from each distinct value \n"
293 " in the dataset; spatial contiguity will not be \n"
294 " used (but you still have to supply rmm and vmul \n"
295 " on the command line). \n"
296 " \n"
297 " N.B.: 'Clusters' formed this way may well have components \n"
298 " that are widely separated! \n"
299 "\n"
300 " -inmask => If 3dClustSim put an internal attribute into the \n"
301 " input dataset that describes a mask, 3dclust will \n"
302 " use this mask to eliminate voxels before clustering, \n"
303 " if you give this option. '-inmask' is how the AFNI \n"
304 " AFNI Clusterize GUI works by default. \n"
305 " [If there is no internal mask in the dataset] \n"
306 " [header, then '-inmask' doesn't do anything.] \n"
307 "\n"
308 " N.B.: The usual way for 3dClustSim to have put this internal \n"
309 " mask into a functional dataset is via afni_proc.py. \n"
310 " \n"
311 " -prefix ppp => Write a new dataset that is a copy of the \n"
312 " input, but with all voxels not in a cluster \n"
313 " set to zero; the new dataset's prefix is 'ppp' \n"
314 " \n"
315 " N.B.: Use of the -prefix option only affects the \n"
316 " first input dataset. \n"
317 "\n"
318 " -savemask q => Write a new dataset that is an ordered mask, such \n"
319 " that the largest cluster is labeled '1', the next \n"
320 " largest '2' and so forth. Should be the same as \n"
321 " '3dmerge -1clust_order' or Clusterize 'SaveMsk'. \n"
322 " -binary => This turns the output of '-savemask' into a binary \n"
323 " (0 or 1) mask, rather than a cluster-index mask. \n"
324 " **-->> If no clusters are found, the mask is not written! \n"
325 "\n"
326 "----------------------------------------------------------------------- \n"
327 " N.B.: 'N.B.' is short for 'Nota Bene', Latin for 'Note Well'; \n"
328 " also see http://en.wikipedia.org/wiki/Nota_bene \n"
329 "----------------------------------------------------------------------- \n"
330 " \n"
331 "E.g., 3dclust -1clip 0.3 5 3000 func+orig'[1]' \n"
332 " \n"
333 " The above command tells 3dclust to find potential cluster volumes for \n"
334 " dataset func+orig, sub-brick #1, where the threshold has been set \n"
335 " to 0.3 (i.e., ignore voxels with activation threshold >0.3 or <-0.3). \n"
336 " Voxels must be no more than 5 mm apart, and the cluster volume \n"
337 " must be at least 3000 micro-liters in size. \n"
338 " \n"
339 "Explanation of 3dclust Output: \n"
340 "----------------------------- \n"
341 " \n"
342 " Volume : Volume that makes up the cluster, in microliters (mm^3)\n"
343 " (or the number of voxels, if -dxyz=1 is given) \n"
344 " \n"
345 " CM RL : Center of mass (CM) for the cluster in the Right-Left \n"
346 " direction (i.e., the coordinates for the CM) \n"
347 " \n"
348 " CM AP : Center of mass for the cluster in the \n"
349 " Anterior-Posterior direction \n"
350 " \n"
351 " CM IS : Center of mass for the cluster in the \n"
352 " Inferior-Superior direction \n"
353 " \n"
354 " minRL, maxRL : Bounding box for the cluster, min and max \n"
355 " coordinates in the Right-Left direction \n"
356 " \n"
357 " minAP, maxAP : Min and max coordinates in the Anterior-Posterior \n"
358 " direction of the volume cluster \n"
359 " \n"
360 " minIS, max IS: Min and max coordinates in the Inferior-Superior \n"
361 " direction of the volume cluster \n"
362 " \n"
363 " Mean : Mean value for the volume cluster \n"
364 " \n"
365 " SEM : Standard Error of the Mean for the volume cluster \n"
366 " \n"
367 " Max Int : Maximum Intensity value for the volume cluster \n"
368 " \n"
369 " MI RL : Coordinate of the Maximum Intensity value in the \n"
370 " Right-Left direction of the volume cluster \n"
371 " \n"
372 " MI AP : Coordinate of the Maximum Intensity value in the \n"
373 " Anterior-Posterior direction of the volume cluster \n"
374 " \n"
375 " MI IS : Coordinate of the Maximum Intensity value in the \n"
376 " Inferior-Superior direction of the volume cluster \n"
377 "----------------------------------------------------------------------- \n"
378 " \n"
379 "Nota Bene: \n"
380 " \n"
381 " * The program does not work on complex- or rgb-valued datasets! \n"
382 " \n"
383 " * Using the -1noneg option is strongly recommended! \n"
384 " \n"
385 " * 3D+time datasets are allowed, but only if you use the \n"
386 " -1tindex and -1dindex options. \n"
387 " \n"
388 " * Bucket datasets are allowed, but you will almost certainly \n"
389 " want to use the -1tindex and -1dindex options with these. \n"
390 " \n"
391 " * SEM values are not realistic for interpolated data sets! \n"
392 " A ROUGH correction is to multiply the SEM of the interpolated \n"
393 " data set by the square root of the number of interpolated \n"
394 " voxels per original voxel. \n"
395 " \n"
396 " * If you use -dxyz=1, then rmm should be given in terms of \n"
397 " voxel edges (not mm) and vmul should be given in terms of \n"
398 " voxel counts (not microliters). Thus, to connect to only \n"
399 " 3D nearest neighbors and keep clusters of 10 voxels or more, \n"
400 " use something like '3dclust -dxyz=1 1.01 10 dset+orig'. \n"
401 " In the report, 'Volume' will be voxel count, but the rest of \n"
402 " the coordinate dependent information will be in actual xyz \n"
403 " millimeters. \n"
404 " \n"
405 " * The default coordinate output order is DICOM. If you prefer \n"
406 " the SPM coordinate order, use the option '-orient LPI' or \n"
407 " set the environment variable AFNI_ORIENT to 'LPI'. For more \n"
408 " information, see file README.environment. \n"
409 ) ;
410 PRINT_COMPILE_DATE ; exit(0) ;
411 }
412
413 mainENTRY("3dclust main"); machdep(); AFNI_logger("3dclust",argc,argv);
414 PRINT_VERSION("3dclust") ; AUTHOR(PROGRAM_AUTHOR) ;
415 INFO_message("*** Consider using program 3dClusterize instead of 3dclust ***") ;
416
417 THD_coorder_fill( my_getenv("AFNI_ORIENT") , &CL_cord ) ; /* July 1997 */
418 CL_read_opts( argc , argv ) ;
419 nopt = CL_nopt ;
420
421 if (CL_1Dform) {
422 sprintf(c1d, "#");
423 sprintf(c1dn, " ");
424 } else {
425 c1d[0] = '\0';
426 c1dn[0] = '\0';
427 }
428 if( CL_do_mni )
429 THD_coorder_fill( "LPI" , &CL_cord ) ; /* 30 Apr 2002 */
430
431 /*----- Identify software -----*/
432 #if 0
433 if( !CL_quiet ){
434 printf ("%s\n%s\n", c1d, c1d);
435 printf ("%sProgram: %s \n", c1d, PROGRAM_NAME);
436 printf ("%sAuthor: %s \n", c1d, PROGRAM_AUTHOR);
437 printf ("%sDate: %s \n", c1d, PROGRAM_DATE);
438 printf ("%s\n", c1d);
439 }
440 #endif
441
442 if( do_NN ){ /* 12 Jul 2017 */
443 CL_edopt.fake_dxyz = 1 ;
444 switch( do_NN ){
445 default:
446 case 1: rmm = 1.11f ; break ;
447 case 2: rmm = 1.44f ; break ;
448 case 3: rmm = 1.77f ; break ;
449 }
450 if( NNvox > 0 ) vmul = (float)NNvox ;
451 else vmul = 2.0f ;
452 } else { /* the OLDE way (with rmm and vmul) */
453 if( nopt+3 > argc )
454 ERROR_exit("No rmm or vmul arguments?") ;
455
456 rmm = strtod( argv[nopt++] , NULL ) ;
457 vmul = strtod( argv[nopt++] , NULL ) ;
458 if( rmm < 0.0 )
459 ERROR_exit("Illegal rmm=%f",rmm) ;
460 else if ( rmm == 0.0f ){
461 CL_edopt.fake_dxyz = 1 ; /* 26 Dec 2007 */
462 rmm = 1.11f ;
463 }
464 }
465
466 /* BDW 26 March 1999 */
467
468 if( CL_edopt.clust_rmm >= 0.0 ){ /* 01 Nov 1999 */
469 WARNING_message("-1clust can't be used in 3dclust") ;
470 CL_edopt.clust_rmm = -1.0 ;
471 }
472
473 /**-- loop over datasets --**/
474
475 dset = NULL ;
476 for( iarg=nopt ; iarg < argc ; iarg++ ){
477 if( dset != NULL ) THD_delete_3dim_dataset(dset,False) ; /* flush old */
478
479 dset = THD_open_dataset( argv[iarg] ) ; /* open new */
480
481 if( dset == NULL ){ /* failed? */
482 ERROR_message("Can't open dataset %s -- skipping it",argv[iarg]) ;
483 continue ;
484 }
485 if( DSET_NUM_TIMES(dset) > 1 &&
486 ( CL_edopt.iv_fim < 0 || CL_edopt.iv_thr < 0 ) ){ /* no time */
487 ERROR_message( /* dependence! */
488 "Cannot use time-dependent dataset %s",argv[iarg]) ;
489 continue ;
490 }
491
492 THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ; /* no mmap */
493 if( CL_verbose )
494 INFO_message("Loading dataset %s",argv[iarg]) ;
495 DSET_load(dset); CHECK_LOAD_ERROR(dset); /* read in */
496 EDIT_one_dataset( dset , &CL_edopt ) ; /* editing? */
497
498 /* search for ASCII mask string [02 Aug 2011] */
499 /* results: nmask == 0 ==> no mask found at all
500 mask != NULL ==> this is the mask
501 nmask > 0 && mask == NULL ==> mask found but not used
502 nmask < 0 ==> mask found but not usable */
503
504 if( mask != NULL ){ free(mask); mask = NULL; } /* clear mask stuff */
505 nmask = 0 ;
506
507 { ATR_string *atr = THD_find_string_atr(dset->dblk,"AFNI_CLUSTSIM_MASK") ;
508 if( atr != NULL ){ /* mask stored as B64 string */
509 nmask = mask_b64string_nvox(atr->ch) ; /* length of mask */
510 if( nmask != DSET_NVOX(dset) ) /* must match dataset */
511 nmask = -1 ;
512 else if( !no_inmask )
513 mask = mask_from_b64string(atr->ch,&nmask) ;
514 } else { /* mask name stored in NN1 NIML header */
515 atr = THD_find_string_atr(dset->dblk,"AFNI_CLUSTSIM_NN1") ;
516 if( atr != NULL ){ /* mask stored as reference to a dataset */
517 NI_element *nel = NI_read_element_fromstring(atr->ch) ;
518 char *nnn ;
519 nnn = NI_get_attribute(nel,"mask_dset_name") ; /* dataset name? */
520 if( nnn == NULL ){
521 nnn = NI_get_attribute(nel,"mask_dset_idcode") ; /* try idcode */
522 if( nnn != NULL ) nmask = -1 ; /* can't use idcode */
523 } else {
524 THD_3dim_dataset *mset = THD_open_dataset(nnn) ; /* try to read */
525 if( mset != NULL ){ /* the dataset */
526 nmask = DSET_NVOX(mset) ;
527 if( nmask != DSET_NVOX(dset) ) /* must match dataset */
528 nmask = -1 ;
529 else if( !no_inmask )
530 mask = THD_makemask(mset,0,1.0f,0.0f) ;
531 DSET_delete(mset) ;
532 } else {
533 nmask = -1 ; /* can't read dataset */
534 }
535 }
536 NI_free_element(nel) ;
537 }
538 }
539 } /* end of trying to load internal mask */
540
541 /* 30 Apr 2002: check if -mni should be used here */
542
543 do_mni = (CL_do_mni && dset->view_type == VIEW_TALAIRACH_TYPE) ;
544
545 if( CL_ivfim < 0 )
546 ivfim = DSET_PRINCIPAL_VALUE(dset) ; /* useful data */
547 else
548 ivfim = CL_ivfim ; /* 16 Sep 1999 */
549
550 /* 02 Aug 2011: mask the data, maybe */
551
552 mri_maskify( DSET_BRICK(dset,ivfim) , mask ) ; /* does nada if mask==NULL */
553
554 /* and get a pointer to the data */
555
556 vfim = DSET_ARRAY(dset,ivfim) ; /* ptr to data */
557 fimfac = DSET_BRICK_FACTOR(dset,ivfim) ; /* scl factor */
558 if( vfim == NULL ){
559 ERROR_message("Cannot access data brick[%d] in dataset %s",ivfim,argv[iarg]) ;
560 continue ;
561 }
562
563 if( !AFNI_GOOD_FUNC_DTYPE( DSET_BRICK_TYPE(dset,ivfim) ) ||
564 DSET_BRICK_TYPE(dset,ivfim) == MRI_rgb ){
565
566 ERROR_message("Illegal datum type in dataset %s",argv[iarg]) ;
567 continue ;
568 }
569
570 nx = dset->daxes->nxx ; dx = fabs(dset->daxes->xxdel) ;
571 ny = dset->daxes->nyy ; dy = fabs(dset->daxes->yydel) ;
572 nz = dset->daxes->nzz ; dz = fabs(dset->daxes->zzdel) ;
573 nxy = nx * ny ; nxyz = nxy * nz ;
574
575 if( CL_edopt.fake_dxyz ){ dxf = dyf = dzf = 1.0 ; } /* 24 Jan 2001 */
576 else { dxf = dx ; dyf = dy ; dzf = dz ; }
577
578 if( vmul >= 0.0 )
579 ptmin = (int) (vmul / (dxf*dyf*dzf) + 0.99) ;
580 else
581 ptmin = (int) fabs(vmul) ; /* 30 Apr 2002 */
582
583 /*-- print report header --*/
584 if( !CL_quiet ){
585
586 if( CL_summarize != 1 ){
587 printf( "%s\n"
588 "%sCluster report for file %s %s\n"
589 #if 0
590 "%s[3D Dataset Name: %s ]\n"
591 "%s[ Short Label: %s ]\n"
592 #endif
593 "%s[Connectivity radius = %.2f mm Volume threshold = %.2f ]\n"
594 "%s[Single voxel volume = %.1f (microliters) ]\n"
595 "%s[Voxel datum type = %s ]\n"
596 "%s[Voxel dimensions = %.3f mm X %.3f mm X %.3f mm ]\n"
597 "%s[Coordinates Order = %s ]\n",
598 c1d,
599 c1d, argv[iarg] , do_mni ? "[MNI coords]" : "" , /* 30 Apr 2002 */
600 #if 0
601 c1d, dset->self_name ,
602 c1d, dset->label1 ,
603 #endif
604 c1d, rmm , ptmin*dx*dy*dz ,
605 c1d, dx*dy*dz ,
606 c1d, MRI_TYPE_name[ DSET_BRICK_TYPE(dset,ivfim) ] ,
607 c1d, dx,dy,dz,
608 c1d, CL_cord.orcode );
609
610 if( CL_edopt.fake_dxyz ) /* 24 Jan 2001 */
611 printf("%s[Fake voxel dimen = %.3f mm X %.3f mm X %.3f mm ]\n",
612 c1d, dxf,dyf,dzf) ;
613
614 if( nmask > 0 && mask != NULL ) /* 02 Aug 2011 */
615 printf("%s[Using internal mask]\n",c1d) ;
616 else if( nmask > 0 )
617 printf("%s[Skipping internal mask]\n",c1d) ;
618 else if( nmask < 0 )
619 printf("%s[Un-usable internal mask]\n",c1d) ; /* should not happen */
620
621 if (CL_noabs) /* BDW 19 Jan 1999 */
622 printf ("%sMean and SEM based on Signed voxel intensities: \n%s\n", c1d, c1d);
623 else
624 printf ("%sMean and SEM based on Absolute Value "
625 "of voxel intensities: \n%s\n", c1d, c1d);
626
627 printf (
628 "%sVolume CM %s CM %s CM %s min%s max%s min%s max%s min%s max%s Mean SEM Max Int MI %s MI %s MI %s\n"
629 "%s------ ----- ----- ----- ----- ----- ----- ----- ----- ----- ------- ------- ------- ----- ----- -----\n",
630
631 c1d,
632 ORIENT_tinystr[ CL_cord.xxor ] ,
633 ORIENT_tinystr[ CL_cord.yyor ] ,
634 ORIENT_tinystr[ CL_cord.zzor ] ,
635 ORIENT_tinystr[ CL_cord.xxor ] , ORIENT_tinystr[ CL_cord.xxor ] ,
636 ORIENT_tinystr[ CL_cord.yyor ] , ORIENT_tinystr[ CL_cord.yyor ] ,
637 ORIENT_tinystr[ CL_cord.zzor ] , ORIENT_tinystr[ CL_cord.zzor ] ,
638 ORIENT_tinystr[ CL_cord.xxor ] ,
639 ORIENT_tinystr[ CL_cord.yyor ] ,
640 ORIENT_tinystr[ CL_cord.zzor ] ,
641 c1d
642 ) ;
643
644 } else {
645 if (CL_noabs) /* BDW 19 Jan 1999 */
646 printf ("%sMean and SEM based on Signed voxel intensities: \n", c1d);
647 else
648 printf ("%sMean and SEM based on Absolute Value "
649 "of voxel intensities: \n", c1d);
650 printf("%sCluster summary for file %s %s\n" ,
651 c1d, argv[iarg] , do_mni ? "[MNI coords]" : "");
652 printf("%sVolume CM %s CM %s CM %s Mean SEM \n", c1d, ORIENT_tinystr[ CL_cord.xxor ],ORIENT_tinystr[ CL_cord.yyor ] ,ORIENT_tinystr[ CL_cord.zzor ]);
653 }
654 } /* end of report header */
655
656 /*-- actually find the clusters in the dataset */
657
658 clar = NIH_find_clusters( nx,ny,nz , dxf,dyf,dzf ,
659 DSET_BRICK_TYPE(dset,ivfim) , vfim , rmm ,
660 CL_edopt.isomode ) ;
661
662 /*-- don't need dataset data any more --*/
663
664 PURGE_DSET( dset ) ;
665
666 if( clar == NULL || clar->num_clu == 0 ){
667 printf("%s** NO CLUSTERS FOUND ***\n", c1d) ;
668 if( AFNI_yesenv("AFNI_3dclust_report_zero") ) printf(" 0\n") ;
669 if( clar != NULL ) DESTROY_CLARR(clar) ;
670 continue ; /* next dataset */
671 }
672
673 /** edit for volume (June 1995) **/
674
675 INIT_CLARR(clbig) ;
676 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
677 cl = clar->clar[iclu] ;
678 if( cl != NULL && cl->num_pt >= ptmin ){ /* big enough */
679 ADDTO_CLARR(clbig,cl) ; /* copy pointer */
680 clar->clar[iclu] = NULL ; /* null out original */
681 }
682 }
683 DESTROY_CLARR(clar) ;
684 clar = clbig ;
685 if( clar == NULL || clar->num_clu == 0 ){
686 printf("%s** NO CLUSTERS FOUND ***\n", c1d) ;
687 if( AFNI_yesenv("AFNI_3dclust_report_zero") ) printf(" 0\n") ;
688 if( clar != NULL ) DESTROY_CLARR(clar) ;
689 continue ;
690 }
691
692 /** end of June 1995 addition **/
693
694 /** sort clusters by size, to make a nice report **/
695
696 if( clar->num_clu < 3333 ){
697 SORT_CLARR(clar) ;
698 } else if( CL_summarize != 1 ){
699 printf("%s** TOO MANY CLUSTERS TO SORT BY VOLUME ***\n", c1d) ;
700 }
701
702 /*-- 29 Nov 2001: write out an edited dataset? --*/
703 if( CL_prefix != NULL || CL_savemask != NULL ){
704
705 if (iarg == nopt) {
706 int qv ; short *mmm=NULL ;
707
708 /* make a mask of voxels to keep */
709
710 mmm = (short *) calloc(sizeof(short),nxyz) ;
711 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
712 cl = clar->clar[iclu] ; if( cl == NULL ) continue ;
713 for( ipt=0 ; ipt < cl->num_pt ; ipt++ ){
714 ii = cl->i[ipt] ; jj = cl->j[ipt] ; kk = cl->k[ipt] ;
715 mmm[ii+jj*nx+kk*nxy] = (do_binary) ? 1 : (iclu+1) ;
716 }
717 }
718
719 if( CL_prefix != NULL ){
720 DSET_load( dset ) ; /* reload data from disk */
721 PREP_LOADED_DSET_4_REWRITE(dset, CL_prefix); /* ZSS Dec 2011 */
722
723 tross_Make_History( "3dclust" , argc , argv , dset ) ;
724
725 /* mask out each sub-brick */
726
727 for( qv=0 ; qv < DSET_NVALS(dset) ; qv++ ){
728
729 switch( DSET_BRICK_TYPE(dset,qv) ){
730
731 case MRI_short:{
732 short *bar = (short *) DSET_ARRAY(dset,qv) ;
733 for( ii=0 ; ii < nxyz ; ii++ )
734 if( mmm[ii] == 0 ) bar[ii] = 0 ;
735 }
736 break ;
737
738 case MRI_byte:{
739 byte *bar = (byte *) DSET_ARRAY(dset,qv) ;
740 for( ii=0 ; ii < nxyz ; ii++ )
741 if( mmm[ii] == 0 ) bar[ii] = 0 ;
742 }
743 break ;
744
745 case MRI_int:{
746 int *bar = (int *) DSET_ARRAY(dset,qv) ;
747 for( ii=0 ; ii < nxyz ; ii++ )
748 if( mmm[ii] == 0 ) bar[ii] = 0 ;
749 }
750 break ;
751
752 case MRI_float:{
753 float *bar = (float *) DSET_ARRAY(dset,qv) ;
754 for( ii=0 ; ii < nxyz ; ii++ )
755 if( mmm[ii] == 0 ) bar[ii] = 0.0 ;
756 }
757 break ;
758
759 case MRI_double:{
760 double *bar = (double *) DSET_ARRAY(dset,qv) ;
761 for( ii=0 ; ii < nxyz ; ii++ )
762 if( mmm[ii] == 0 ) bar[ii] = 0.0 ;
763 }
764 break ;
765
766 case MRI_complex:{
767 complex *bar = (complex *) DSET_ARRAY(dset,qv) ;
768 for( ii=0 ; ii < nxyz ; ii++ )
769 if( mmm[ii] == 0 ) bar[ii].r = bar[ii].i = 0.0 ;
770 }
771 break ;
772
773 case MRI_rgb:{
774 byte *bar = (byte *) DSET_ARRAY(dset,qv) ;
775 for( ii=0 ; ii < nxyz ; ii++ )
776 if( mmm[ii] == 0 ) bar[3*ii] = bar[3*ii+1] = bar[3*ii+2] = 0 ;
777 }
778 break ;
779 } /* end of switch over sub-brick type */
780 } /* end of loop over sub-bricks */
781
782 /* write dataset out */
783
784 DSET_write(dset) ; WROTE_DSET(dset) ; PURGE_DSET(dset) ;
785 }
786
787 if( CL_savemask != NULL ){ /* 26 Aug 2011 */
788 THD_3dim_dataset *qset ;
789 qset = EDIT_empty_copy(dset) ;
790 EDIT_dset_items( qset ,
791 ADN_prefix , CL_savemask ,
792 ADN_nvals , 1 ,
793 ADN_none ) ;
794 EDIT_substitute_brick(qset,0,MRI_short,mmm) ; mmm = NULL ;
795 tross_Copy_History( dset , qset ) ;
796 tross_Make_History( "3dclust" , argc , argv , qset ) ;
797 DSET_write(qset) ; WROTE_DSET(qset) ; DSET_delete(qset) ;
798 }
799
800 if( mmm != NULL ) free(mmm) ;
801
802 } else { /** Bad news **/
803 WARNING_message(
804 "Output volume not written for input %s . You either \n"
805 "have bad datasets on the command line (check output warnings),\n"
806 "or multiple valid datasets as input. In the latter case, \n"
807 "-prefix does not work.\n",
808 DSET_BRIKNAME(dset));
809 }
810
811 } /* end of saving a dataset */
812
813 ndet = 0 ;
814
815 vol_total = nvox_total = 0 ;
816 glmmsum = glmssum = glsqsum = glxxsum = glyysum = glzzsum = 0;
817
818 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
819 cl = clar->clar[iclu] ;
820 if( cl == NULL || cl->num_pt < ptmin ) continue ; /* no good */
821
822 volsum = cl->num_pt * dxf*dyf*dzf ;
823 xxsum = yysum = zzsum = mmsum = mssum = 0.0 ;
824 xxmax = yymax = zzmax = mmmax = msmax = 0.0 ;
825 sqsum = sem = 0;
826
827 /* These should be pegged at whatever actual max/min values are */
828 RLmax = APmax = ISmax = -1000;
829 RLmin = APmin = ISmin = 1000;
830
831 for( ipt=0 ; ipt < cl->num_pt ; ipt++ ){
832
833 #if 0
834 /** this is obsolete and nonfunctional code **/
835 IJK_TO_THREE( cl->ijk[ipt] , ii,jj,kk , nx,nxy ) ;
836 #endif
837 ii = cl->i[ipt] ; jj = cl->j[ipt] ; kk = cl->k[ipt] ;
838
839 fv = THD_3dind_to_3dmm( dset , TEMP_IVEC3(ii,jj,kk) ) ;
840 fv = THD_3dmm_to_dicomm( dset , fv ) ;
841 xx = fv.xyz[0] ; yy = fv.xyz[1] ; zz = fv.xyz[2] ;
842 if( !do_mni )
843 THD_dicom_to_coorder( &CL_cord , &xx,&yy,&zz ) ; /* July 1997 */
844 else
845 THD_3tta_to_3mni( &xx , &yy , &zz ) ; /* 30 Apr 2002 */
846
847 ms = cl->mag[ipt]; /* BDW 18 Jan 1999 */
848 mm = fabs(ms);
849
850 mssum += ms;
851 mmsum += mm;
852
853 sqsum += mm * mm;
854 xxsum += mm * xx ; yysum += mm * yy ; zzsum += mm * zz ;
855 if( mm > mmmax ){
856 xxmax = xx ; yymax = yy ; zzmax = zz ;
857 mmmax = mm ; msmax = ms ;
858 }
859
860 /* Dimensions: */
861 if ( xx > RLmax )
862 RLmax = xx;
863 if ( xx < RLmin )
864 RLmin = xx;
865 if ( yy > APmax )
866 APmax = yy;
867 if ( yy < APmin )
868 APmin = yy;
869 if ( zz > ISmax )
870 ISmax = zz;
871 if ( zz < ISmin )
872 ISmin = zz;
873
874 }
875 if( mmsum == 0.0 ) continue ;
876
877 glmssum += mssum;
878 glmmsum += mmsum;
879 glsqsum += sqsum ;
880 glxxsum += xxsum;
881 glyysum += yysum;
882 glzzsum += zzsum;
883
884 ndet++ ;
885 xxsum /= mmsum ; yysum /= mmsum ; zzsum /= mmsum ;
886
887 if (CL_noabs) mean = mssum / cl->num_pt; /* BDW 19 Jan 1999 */
888 else mean = mmsum / cl->num_pt;
889
890 if( fimfac != 0.0 )
891 { mean *= fimfac; msmax *= fimfac;
892 sqsum *= fimfac*fimfac; } /* BDW 12/27/96 */
893
894 /* MSB 11/1/96 Calculate SEM using SEM^2=s^2/N,
895 where s^2 = (SUM Y^2)/N - (Ymean)^2
896 where sqsum = (SUM Y^2 ) */
897
898 if (cl->num_pt > 1)
899 {
900 sem = (sqsum - (cl->num_pt * mean * mean)) / (cl->num_pt - 1);
901 if (sem > 0.0) sem = sqrt( sem / cl->num_pt ); else sem = 0.0;
902 }
903 else
904 sem = 0.0;
905
906 if( CL_summarize != 1 ){
907 MCW_fc7(mean, buf1) ;
908 MCW_fc7(msmax,buf2) ;
909 MCW_fc7(sem ,buf3) ;
910
911 printf("%s%6.0f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %7s %7s %7s %5.1f %5.1f %5.1f \n",
912 c1dn, volsum, xxsum, yysum, zzsum, RLmin, RLmax, APmin, APmax, ISmin, ISmax, buf1, buf3, buf2, xxmax, yymax, zzmax ) ;
913 }
914
915 nvox_total += cl->num_pt ;
916 vol_total += volsum ;
917
918 }
919
920 DESTROY_CLARR(clar) ;
921 if( ndet == 0 ){
922 printf("%s** NO CLUSTERS FOUND ABOVE THRESHOLD VOLUME ***\n", c1d) ;
923 if( AFNI_yesenv("AFNI_3dclust_report_zero") ) printf(" 0\n") ;
924 }
925
926
927 /* MSB 11/1/96 Calculate global SEM */
928
929 if (CL_noabs) glmean = glmssum / nvox_total; /* BDW 19 Jan 1999 */
930 else glmean = glmmsum / nvox_total;
931
932 /* BDW 12/27/96 */
933 if( fimfac != 0.0 ){ glsqsum *= fimfac*fimfac ; glmean *= fimfac ; }
934 if (nvox_total > 1)
935 {
936 sem = (glsqsum - (nvox_total*glmean*glmean)) / (nvox_total - 1);
937 if (sem > 0.0) sem = sqrt( sem / nvox_total ); else sem = 0.0;
938 }
939 else
940 sem = 0.0;
941
942 glxxsum /= glmmsum ; glyysum /= glmmsum ; glzzsum /= glmmsum ;
943
944 /* MSB 11/1/96 Modified so that mean and SEM would print in correct column */
945 if( CL_summarize == 1 )
946 { if( !CL_quiet )
947 printf( "%s------ ----- ----- ----- -------- -------- \n", c1d);
948 printf("%s%6.0f %5.1f %5.1f %5.1f %8.1f %6.3f\n" , c1d, vol_total, glxxsum, glyysum, glzzsum, glmean, sem ) ;
949 }
950 else if( ndet > 1 && CL_summarize != -1 )
951 {
952 MCW_fc7(glmean ,buf1) ;
953 MCW_fc7(sem ,buf3) ;
954 if( !CL_quiet )
955 printf ("%s------ ----- ----- ----- ----- ----- ----- ----- ----- ----- ------- ------- ------- ----- ----- -----\n", c1d);
956 printf ("%s%6.0f %5.1f %5.1f %5.1f %7s %7s \n",
957 c1d, vol_total, glxxsum, glyysum, glzzsum, buf1, buf3 ) ;
958 }
959
960 }
961
962 exit(0) ;
963 }
964
965
966 /*---------------------------------------------------------------------------*/
967 /*
968 read the arguments, and load the global variables
969 */
970
971 #ifdef CLDEBUG
972 # define DUMP1 fprintf(stderr,"ARG: %s\n",argv[nopt])
973 # define DUMP2 fprintf(stderr,"ARG: %s %s\n",argv[nopt],argv[nopt+1])
974 # define DUMP3 fprintf(stderr,"ARG: %s %s %s\n",argv[nopt],argv[nopt+1],argv[nopt+2])
975 #else
976 # define DUMP1
977 # define DUMP2
978 # define DUMP3
979 #endif
980
CL_read_opts(int argc,char * argv[])981 void CL_read_opts( int argc , char * argv[] )
982 {
983 int nopt = 1 ;
984 int ival ;
985
986 INIT_EDOPT( &CL_edopt ) ;
987
988 while( nopt < argc && argv[nopt][0] == '-' ){
989
990 /**** check for editing options ****/
991
992 ival = EDIT_check_argv( argc , argv , nopt , &CL_edopt ) ;
993 if( ival > 0 ){
994 nopt += ival ;
995 continue ;
996 }
997
998 #if 0 /* These two are now captured in EDIT_check_argv,
999 remove this block next time you see it. ZSS March 2010 */
1000 /**** 30 Apr 2002: -isovalue and -isomerge ****/
1001 if( strcmp(argv[nopt],"-isovalue") == 0 ){
1002 CL_isomode = ISOVALUE_MODE ;
1003 nopt++ ; continue ;
1004 }
1005
1006 if( strcmp(argv[nopt],"-isomerge") == 0 ){
1007 CL_isomode = ISOMERGE_MODE ;
1008 nopt++ ; continue ;
1009 }
1010 #endif
1011
1012 if( strcmp(argv[nopt],"-NN1") == 0 ){ /* 12 Jul 2017 */
1013 do_NN = 1 ; nopt++ ;
1014 if( nopt < argc && isdigit(argv[nopt][0]) )
1015 NNvox = (int)strtod(argv[nopt++],NULL) ;
1016 else
1017 NNvox = 0 ;
1018 continue ;
1019 }
1020 if( strcmp(argv[nopt],"-NN2") == 0 ){
1021 do_NN = 2 ; nopt++ ;
1022 if( nopt < argc && isdigit(argv[nopt][0]) )
1023 NNvox = (int)strtod(argv[nopt++],NULL) ;
1024 else
1025 NNvox = 0 ;
1026 continue ;
1027 }
1028 if( strcmp(argv[nopt],"-NN3") == 0 ){
1029 do_NN = 3 ; nopt++ ;
1030 if( nopt < argc && isdigit(argv[nopt][0]) )
1031 NNvox = (int)strtod(argv[nopt++],NULL) ;
1032 else
1033 NNvox = 0 ;
1034 continue ;
1035 }
1036 if( strcmp(argv[nopt],"-NN") == 0 ){
1037 nopt++ ;
1038 if( nopt >= argc ) ERROR_exit("need argument after '-NN'") ;
1039 do_NN = (int)strtod(argv[nopt],NULL) ;
1040 if( do_NN < 1 || do_NN > 3 )
1041 ERROR_exit("Illegal value '%s' after '-NN'",argv[nopt]) ;
1042 nopt++ ;
1043 if( nopt < argc && isdigit(argv[nopt][0]) )
1044 NNvox = (int)strtod(argv[nopt++],NULL) ;
1045 else
1046 NNvox = 0 ;
1047 continue ;
1048 }
1049
1050 if( strcmp(argv[nopt],"-no_inmask") == 0 ){ /* 02 Aug 2011 */
1051 no_inmask = 1 ; nopt++ ; continue ;
1052 }
1053 if( strcmp(argv[nopt],"-inmask") == 0 ){
1054 no_inmask = 0 ; nopt++ ; continue ;
1055 }
1056
1057 /**** 30 Apr 2002: -mni ****/
1058
1059 if( strcmp(argv[nopt],"-mni") == 0 ){
1060 CL_do_mni = 1 ;
1061 nopt++ ; continue ;
1062 }
1063
1064 /**** 29 Nov 2001: -prefix ****/
1065
1066 if( strcmp(argv[nopt],"-prefix") == 0 ){
1067 if( ++nopt >= argc ){
1068 ERROR_exit("need an argument after -prefix!") ;
1069 }
1070 CL_prefix = argv[nopt] ;
1071 if( !THD_filename_ok(CL_prefix) ){
1072 ERROR_exit("-prefix string is illegal: %s\n",CL_prefix);
1073 }
1074 nopt++ ; continue ;
1075 }
1076
1077 /**** 26 Aug 2011: -savemask ****/
1078
1079 if( strcmp(argv[nopt],"-savemask") == 0 ){
1080 if( ++nopt >= argc ){
1081 ERROR_exit("need an argument after -savemask!\n") ;
1082 }
1083 CL_savemask = argv[nopt] ;
1084 if( !THD_filename_ok(CL_savemask) ){
1085 ERROR_exit("-savemask string is illegal: %s\n",CL_savemask);
1086 }
1087 nopt++ ; continue ;
1088 }
1089
1090 if( strcmp(argv[nopt],"-binary") == 0 ){ /* 26 May 2017 */
1091 do_binary = 1 ; nopt++ ; continue ;
1092 }
1093
1094 /**** Sep 16 1999: -1tindex and -1dindex ****/
1095
1096 if( strncmp(argv[nopt],"-1dindex",5) == 0 ){
1097 if( ++nopt >= argc ){
1098 ERROR_exit("need an argument after -1dindex!\n") ;
1099 }
1100 CL_ivfim = CL_edopt.iv_fim = (int) strtod( argv[nopt++] , NULL ) ;
1101 continue ;
1102 }
1103
1104 if( strncmp(argv[nopt],"-1tindex",5) == 0 ){
1105 if( ++nopt >= argc ){
1106 ERROR_exit("need an argument after -1tindex!\n") ;
1107 }
1108 CL_ivthr = CL_edopt.iv_thr = (int) strtod( argv[nopt++] , NULL ) ;
1109 continue ;
1110 }
1111
1112 /**** -summarize ****/
1113
1114 if( strncmp(argv[nopt],"-summarize",5) == 0 ){
1115 CL_summarize = 1 ;
1116 nopt++ ; continue ;
1117 }
1118
1119 /**** -nosum [05 Jul 2001] ****/
1120
1121 if( strncmp(argv[nopt],"-nosum",6) == 0 ){
1122 CL_summarize = -1 ;
1123 nopt++ ; continue ;
1124 }
1125
1126 /**** -verbose ****/
1127
1128 if( strncmp(argv[nopt],"-verbose",5) == 0 ){
1129 CL_verbose = CL_edopt.verbose = 1 ;
1130 nopt++ ; continue ;
1131 }
1132
1133 /**** -quiet ****/
1134
1135 if( strncmp(argv[nopt],"-quiet",5) == 0 ){
1136 CL_quiet = 1 ;
1137 nopt++ ; continue ;
1138 }
1139
1140 /**** -1Dformat ****/
1141
1142 if( strncmp(argv[nopt],"-1Dformat",5) == 0 ){
1143 CL_1Dform = 1 ;
1144 nopt++ ; continue ;
1145 }
1146
1147 if( strncmp(argv[nopt],"-no_1Dformat",7) == 0 ){
1148 CL_1Dform = 0 ;
1149 nopt++ ; continue ;
1150 }
1151
1152 /**** -orient code ****/
1153
1154 if( strncmp(argv[nopt],"-orient",5) == 0 ){
1155 THD_coorder_fill( argv[++nopt] , &CL_cord ) ; /* July 1997 */
1156 nopt++ ; continue ;
1157 }
1158
1159 /**** -noabs ****/ /* BDW 19 Jan 1999 */
1160
1161 if( strncmp(argv[nopt],"-noabs",6) == 0 ){
1162 CL_noabs = 1 ;
1163 nopt++ ; continue ;
1164 }
1165
1166 /**** unknown switch ****/
1167
1168 ERROR_message("Unrecognized option %s",argv[nopt]) ;
1169 suggest_best_prog_option(argv[0], argv[nopt]);
1170 exit(1) ;
1171
1172 } /* end of loop over options */
1173
1174 #ifdef CLDEBUG
1175 printf("** Finished with options\n") ;
1176 #endif
1177
1178 CL_nopt = nopt ;
1179 return ;
1180 }
1181
1182 /*---------------------------------------------------------------------------*/
1183
MCW_fc7(float qval,char * buf)1184 void MCW_fc7( float qval , char * buf )
1185 {
1186 float aval = fabs(qval) ;
1187 int lv , il ;
1188 char lbuf[16] ;
1189
1190 /* special case if the value is an integer */
1191
1192 lv = (int) qval ;
1193
1194 if( qval == lv && abs(lv) < 100000 ){
1195 if( lv >= 0 ) sprintf( buf , " %d" , lv ) ;
1196 else sprintf( buf , "%d" , lv ) ;
1197 return ;
1198 }
1199
1200 /* macro to strip trailing zeros from output */
1201
1202 #define BSTRIP \
1203 for( il=6 ; il>1 && lbuf[il]=='0' ; il-- ) lbuf[il] = '\0'
1204
1205 /* noninteger: choose floating format based on magnitude */
1206
1207 lv = (int) (10.0001 + log10(aval)) ;
1208
1209 switch( lv ){
1210
1211 default:
1212 if( qval > 0.0 ) sprintf( lbuf , "%7.1e" , qval ) ;
1213 else sprintf( lbuf , "%7.0e" , qval ) ;
1214 break ;
1215
1216 case 7: /* 0.001 -0.01 */
1217 case 8: /* 0.01 -0.1 */
1218 case 9: /* 0.1 -1 */
1219 case 10: /* 1 -9.99 */
1220 sprintf( lbuf , "%7.4f" , qval ) ; BSTRIP ; break ;
1221
1222 case 11: /* 10-99.9 */
1223 sprintf( lbuf , "%7.3f" , qval ) ; BSTRIP ; break ;
1224
1225 case 12: /* 100-999.9 */
1226 sprintf( lbuf , "%7.2f" , qval ) ; BSTRIP ; break ;
1227
1228 case 13: /* 1000-9999.9 */
1229 sprintf( lbuf , "%7.1f" , qval ) ; BSTRIP ; break ;
1230
1231 case 14: /* 10000-99999.9 */
1232 sprintf( lbuf , "%7.0f" , qval ) ; break ;
1233
1234 }
1235 strcpy(buf,lbuf) ;
1236 return ;
1237 }
1238 /*---------------------------------------------------------------------------*/
1239