1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 extern int *z_rand_order(int bot, int top, long int seed);
10
11 /*----------------
12 Another quickie.
13 ------------------*/
14
main(int argc,char * argv[])15 int main( int argc , char * argv[] )
16 {
17 int narg , nvox=0 , ii,jj,kk,vv , mcount , iv , mc , ndset,ndval , i,j,k ;
18 THD_3dim_dataset *mask_dset=NULL , **input_dset=NULL ;
19 float mask_bot=666.0 , mask_top=-666.0 ;
20 byte *mmm = NULL ;
21 char *oname = NULL , *obuf , *otemp ;
22 FILE *ofile ;
23 MRI_IMAGE *flim ;
24 float *flar ;
25 int no_ijk=0 , yes_xyz=0, lpi_xyz = 0 ;
26 int yes_index=0 ; /*-- 09 May 2003 [rickr] --*/
27 byte *cmask=NULL ; int ncmask=0 ;
28 int verb=1 ;
29 int yes_niml=0 , no_zero=0 , numz ; /* 04 Feb 2008 */
30 int nout ; char **bout , *niml_name="maskdump" ;
31 int nrand = -1;
32 byte *bmask = NULL ; /*-- box+ball mask: moved here 09 Sep 2009 --*/
33
34 int box_num=0 ; float *box_dat=NULL ; /* 09 May 2003 - RWCox */
35 int ball_num=0; float *ball_dat=NULL; /* 09 Sep 2009 - RWCox */
36 int nx=0,ny=0,nz=0,nxy=0,nxyz ;
37 unsigned int nrandseed = 1234u;
38 float dx, dy, dz; /* scale mm to voxels [24 Nov 2021 rickr] */
39
40 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
41 printf(
42 "Usage: 3dmaskdump [options] dataset dataset ...\n"
43 "Writes to an ASCII file values from the input datasets\n"
44 "which satisfy the mask criteria given in the options.\n"
45 "If no options are given, then all voxels are included.\n"
46 "This might result in a GIGANTIC output file.\n"
47 "Options:\n"
48 " -mask mset Means to use the dataset 'mset' as a mask:\n"
49 " Only voxels with nonzero values in 'mset'\n"
50 " will be printed from 'dataset'. Note\n"
51 " that the mask dataset and the input dataset\n"
52 " must have the same number of voxels.\n"
53 " -mrange a b Means to further restrict the voxels from\n"
54 " 'mset' so that only those mask values\n"
55 " between 'a' and 'b' (inclusive) will\n"
56 " be used. If this option is not given,\n"
57 " all nonzero values from 'mset' are used.\n"
58 " Note that if a voxel is zero in 'mset', then\n"
59 " it won't be included, even if a < 0 < b.\n"
60 /*-- 09 May 2003: add -index option [rickr] */
61 " -index Means to write out the dataset index values.\n"
62 " -noijk Means not to write out the i,j,k values.\n"
63 " -xyz Means to write the x,y,z coordinates from\n"
64 " the 1st input dataset at the start of each\n"
65 " output line. These coordinates are in\n"
66 " the 'RAI' (DICOM) order.\n"
67 " -o fname Means to write output to file 'fname'.\n"
68 " [default = stdout, which you won't like]\n"
69 "\n"
70 " -cmask 'opts' Means to execute the options enclosed in single\n"
71 " quotes as a 3dcalc-like program, and produce\n"
72 " produce a mask from the resulting 3D brick.\n"
73 " Examples:\n"
74 " -cmask '-a fred+orig[7] -b zork+orig[3] -expr step(a-b)'\n"
75 " produces a mask that is nonzero only where\n"
76 " the 7th sub-brick of fred+orig is larger than\n"
77 " the 3rd sub-brick of zork+orig.\n"
78 " -cmask '-a fred+orig -expr 1-bool(k-7)'\n"
79 " produces a mask that is nonzero only in the\n"
80 " 7th slice (k=7); combined with -mask, you\n"
81 " could use this to extract just selected voxels\n"
82 " from particular slice(s).\n"
83 " Notes: * You can use both -mask and -cmask in the same\n"
84 " run - in this case, only voxels present in\n"
85 " both masks will be dumped.\n"
86 " * Only single sub-brick calculations can be\n"
87 " used in the 3dcalc-like calculations -\n"
88 " if you input a multi-brick dataset here,\n"
89 " without using a sub-brick index, then only\n"
90 " its 0th sub-brick will be used.\n"
91 " * Do not use quotes inside the 'opts' string!\n"
92 "\n"
93 " -xbox x y z Means to put a 'mask' down at the dataset (not DICOM)\n"
94 " coordinates of 'x y z' mm. By default, this box is\n"
95 " 1 voxel wide in each direction. You can specify\n"
96 " instead a range of coordinates using a colon ':'\n"
97 " after the coordinates; for example:\n"
98 " -xbox 22:27 31:33 44\n"
99 " means a box from (x,y,z)=(22,31,44) to (27,33,44).\n"
100 " NOTE: dataset coordinates are NOT the coordinates you\n"
101 " typically see in AFNI's main controller top left corner.\n"
102 " Those coordinates are typically in either RAI/DICOM order\n"
103 " or in LPI/SPM order and should be used with -dbox and\n"
104 " -nbox, respectively.\n"
105 "\n"
106 " -dbox x y z Means the same as -xbox, but the coordinates are in\n"
107 " RAI/DICOM order (+x=Left, +y=Posterior, +z=Superior).\n"
108 " If your AFNI environment variable AFNI_ORIENT is set to\n"
109 " RAI, these coordinates correspond to those you'd enter\n"
110 " into the 'Jump to (xyz)' control in AFNI, and to\n"
111 " those output by 3dclust.\n"
112 " NOTE: It is possible to make AFNI and/or 3dclust output \n"
113 " coordinates in an order different from the one specified \n"
114 " by AFNI_ORIENT, but you'd have to work hard on that. \n"
115 " In any case, the order is almost always specified along \n"
116 " with the coordinates. If you see RAI/DICOM, then use \n"
117 " -dbox. If you see LPI/SPM then use -nbox. \n"
118 "\n"
119 " -nbox x y z Means the same as -xbox, but the coordinates are in\n"
120 " LPI/SPM or 'neuroscience' order where the signs of the\n"
121 " x and y coordinates are reversed relative to RAI/DICOM.\n"
122 " (+x=Right, +y=Anterior, +z=Superior)\n"
123 "\n"
124 " -ibox i j k Means to put a 'mask' down at the voxel indexes\n"
125 " given by 'i j k'. By default, this picks out\n"
126 " just 1 voxel. Again, you can use a ':' to specify\n"
127 " a range (now in voxels) of locations.\n"
128 " Notes: * Boxes are cumulative; that is, if you specify more\n"
129 " than 1 box, you'll get more than one region.\n"
130 " * If a -mask and/or -cmask option is used, then\n"
131 " the INTERSECTION of the boxes with these masks\n"
132 " determines which voxels are output; that is,\n"
133 " a voxel must be inside some box AND inside the\n"
134 " mask in order to be selected for output.\n"
135 " * If boxes select more than 1 voxel, the output lines\n"
136 " are NOT necessarily in the order of the options on\n"
137 " the command line.\n"
138 " * Coordinates (for -xbox, -dbox, and -nbox) are relative\n"
139 " to the first dataset on the command line.\n"
140 " * It may be helpful to slightly pad boxes, to be sure they\n"
141 " contain the desired voxel centers.\n"
142 "\n"
143 " -xball x y z r Means to put a ball (sphere) mask down at dataset\n"
144 " coordinates (x,y,z) with radius r.\n"
145 " -dball x y z r Same, but (x,y,z) are in RAI/DICOM order.\n"
146 " -nball x y z r Same, but (x,y,z) are in LPI/SPM order.\n"
147 " Notes: * The combined (set UNION) of all ball and/or box masks\n"
148 " is created first. Then, if a -mask and/or -cmask\n"
149 " option was used, then the ball+box mask will be\n"
150 " INTERSECTED with the existing mask.\n"
151 " * Balls not centered over voxels, or those applied to\n"
152 " anisotropic volumes may not appear symmetric.\n"
153 " * Consider slight padding to handle truncation.\n"
154 "\n"
155 " -nozero Means to skip output of any voxel where all the\n"
156 " data values are zero.\n"
157 "\n"
158 " -n_rand N_RAND Means to keep only N_RAND randomly selected\n"
159 " voxels from what would have been the output.\n"
160 "\n"
161 " -n_randseed SEED Seed the random number generator with SEED,\n"
162 " instead of the default seed of %u\n"
163 "\n"
164 " -niml name Means to output data in the XML/NIML format that\n"
165 " is compatible with input back to AFNI via\n"
166 " the READ_NIML_FILE command.\n"
167 " * 'name' is the 'target_name' for the NIML header\n"
168 " field, which is the name that will be assigned\n"
169 " to the dataset when it is sent into AFNI.\n"
170 " * Also implies '-noijk' and '-xyz' and '-nozero'.\n"
171 "\n"
172 " -quiet Means not to print progress messages to stderr.\n"
173 "\n"
174 "Inputs after the last option are datasets whose values you\n"
175 "want to be dumped out. These datasets (and the mask) can\n"
176 "use the sub-brick selection mechanism (described in the\n"
177 "output of '3dcalc -help') to choose which values you get.\n"
178 "\n"
179 "Each selected voxel gets one line of output:\n"
180 " i j k val val val ....\n"
181 "where (i,j,k) = 3D index of voxel in the dataset arrays,\n"
182 "and val = the actual voxel value. Note that if you want\n"
183 "the mask value to be output, you have to include that\n"
184 "dataset in the dataset input list again, after you use\n"
185 "it in the '-mask' option.\n"
186 "\n"
187 "* To eliminate the 'i j k' columns, use the '-noijk' option.\n"
188 "* To add spatial coordinate columns, use the '-xyz' option.\n"
189 "\n"
190 "N.B.: This program doesn't work with complex-valued datasets!\n"
191 , nrandseed) ;
192
193 printf("\n" MASTER_SHORTHELP_STRING ) ;
194
195 PRINT_COMPILE_DATE ; exit(0) ;
196 }
197
198 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
199
200 mainENTRY("3dmaskdump main"); machdep() ;
201
202
203 { int new_argc ; char ** new_argv ;
204 addto_args( argc , argv , &new_argc , &new_argv ) ;
205 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
206 }
207
208 AFNI_logger("3dmaskdump",argc,argv) ;
209
210 /* scan argument list */
211
212 narg = 1 ;
213 while( narg < argc && argv[narg][0] == '-' ){
214
215 if( strcasecmp(argv[narg],"-niml") == 0 ){ /* 04 Feb 2008 */
216 yes_niml = 1 ; narg++ ;
217 if( narg >= argc ) ERROR_exit("-niml must be followed by a name") ;
218 niml_name = strdup(argv[narg]) ;
219 if( !THD_filename_pure(niml_name) )
220 ERROR_exit("Illegal 'name' after '-niml'") ;
221 narg++ ; continue ;
222 }
223
224 if( strcasecmp(argv[narg],"-nozero") == 0 ){ /* 04 Feb 2008 */
225 no_zero = 1 ; narg++ ; continue ;
226 }
227
228
229 if( strcmp(argv[narg],"-quiet") == 0 ){ /* 09 May 2003 - RWC */
230 verb = 0 ; narg++ ; continue ;
231 }
232
233 if( strcmp(argv[narg]+2,"box") == 0 ){ /* 09 May 2003 - RWC */
234 float xbot,xtop , ybot,ytop , zbot,ztop , btyp=0.0 ;
235 int nn ;
236 char code = *(argv[narg]+1) ; /* should be 'x', 'd' , 'n', or 'i' */
237 switch( code ){
238 case 'x': btyp = BOX_XYZ ; break ;
239 case 'd': btyp = BOX_DIC ; break ;
240 case 'n': btyp = BOX_NEU ; break ;
241 case 'i': btyp = BOX_IJK ; break ;
242 default: ERROR_exit("Unknown 'box' option %s\n",argv[narg]) ;
243 }
244 if( narg+3 >= argc )
245 ERROR_exit("need 3 arguments after %s\n",argv[narg]);
246 nn = sscanf( argv[narg+1] , "%f:%f" , &xbot , &xtop ) ;
247 if( nn < 1 )
248 ERROR_exit("Can't decode %s after %s\n",argv[narg+1],argv[narg]);
249 else if( nn == 1 )
250 xtop=xbot ;
251 nn = sscanf( argv[narg+2] , "%f:%f" , &ybot , &ytop ) ;
252 if( nn < 1 )
253 ERROR_exit("Can't decode %s after %s\n",argv[narg+2],argv[narg]);
254 else if( nn == 1 )
255 ytop=ybot ;
256 nn = sscanf( argv[narg+3] , "%f:%f" , &zbot , &ztop ) ;
257 if( nn < 1 )
258 ERROR_exit("Can't decode %s after %s\n",argv[narg+3],argv[narg]);
259 else if( nn == 1 )
260 ztop=zbot ;
261 box_dat = (float *) realloc( box_dat , sizeof(float)*BOXLEN*(box_num+1) ) ;
262 box_dat[0+BOXLEN*box_num] = btyp ;
263 box_dat[1+BOXLEN*box_num] = xbot ;
264 box_dat[2+BOXLEN*box_num] = xtop ;
265 box_dat[3+BOXLEN*box_num] = ybot ;
266 box_dat[4+BOXLEN*box_num] = ytop ;
267 box_dat[5+BOXLEN*box_num] = zbot ;
268 box_dat[6+BOXLEN*box_num] = ztop ;
269 box_num++ ; narg += 4 ; continue ;
270 }
271
272 if( strcmp(argv[narg]+2,"ball") == 0 ){ /* 09 Sep 2009 - RWC */
273 float xcen,ycen,zcen,rad , btyp=0.0 ;
274 char code = *(argv[narg]+1) ; /* should be 'x', 'd' , or 'n' */
275 switch( code ){
276 case 'x': btyp = BALL_XYZ ; break ;
277 case 'd': btyp = BALL_DIC ; break ;
278 case 'n': btyp = BALL_NEU ; break ;
279 default: ERROR_exit("Unknown 'ball' option %s",argv[narg]) ;
280 }
281 if( narg+4 >= argc )
282 ERROR_exit("need 4 arguments after %s\n",argv[narg]);
283 xcen = strtod( argv[++narg] , NULL ) ;
284 ycen = strtod( argv[++narg] , NULL ) ;
285 zcen = strtod( argv[++narg] , NULL ) ;
286 rad = strtod( argv[++narg] , NULL ) ;
287 if( rad <= 0.0f ){
288 WARNING_message("%s radius=%s !?",argv[narg-4],argv[narg]) ; rad = 0.0f;
289 }
290
291 ball_dat = (float *) realloc( ball_dat , sizeof(float)*BOXLEN*(ball_num+1) ) ;
292 ball_dat[0+BOXLEN*ball_num] = btyp ;
293 ball_dat[1+BOXLEN*ball_num] = xcen ;
294 ball_dat[2+BOXLEN*ball_num] = ycen ;
295 ball_dat[3+BOXLEN*ball_num] = zcen ;
296 ball_dat[4+BOXLEN*ball_num] = rad ;
297 ball_num++ ; narg++ ;
298 continue ;
299 }
300
301 /*-- 09 May 2003: option to output index value (for Mike B) [rickr] --*/
302
303 if( strcmp(argv[narg],"-index") == 0 ){
304 yes_index = 1 ;
305 narg++ ; continue ;
306 }
307
308 if( strcmp(argv[narg],"-noijk") == 0 ){
309 no_ijk = 1 ;
310 narg++ ; continue ;
311 }
312
313 if( strcmp(argv[narg],"-xyz") == 0 ){ /* 23 Mar 2003 */
314 yes_xyz = 1 ;
315 narg++ ; continue ;
316 }
317
318 if( strcmp(argv[narg],"-lpi") == 0 ){ /* 19 Apr 2013 */
319 lpi_xyz = 1 ;
320 narg++ ; continue ;
321 }
322
323 if( strcmp(argv[narg],"-cmask") == 0 ){ /* 16 Mar 2000 */
324 if( narg+1 >= argc )
325 ERROR_exit("-cmask option requires a following argument!\n");
326 cmask = EDT_calcmask( argv[++narg] , &ncmask, 0 ) ;
327 if( cmask == NULL ) ERROR_exit("Can't compute -cmask!\n");
328 narg++ ; continue ;
329 }
330
331 if( strncmp(argv[narg],"-mask",5) == 0 ){
332 if( mask_dset != NULL )
333 ERROR_exit("Cannot have two -mask options!\n") ;
334 if( narg+1 >= argc )
335 ERROR_exit("-mask option requires a following argument!\n");
336 mask_dset = THD_open_dataset( argv[++narg] ) ;
337 if( mask_dset == NULL )
338 ERROR_exit("Cannot open mask dataset!\n") ;
339 if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex )
340 ERROR_exit("Cannot deal with complex-valued mask dataset!\n");
341 narg++ ; continue ;
342 }
343
344 if( strncmp(argv[narg],"-mrange",5) == 0 ){
345 if( narg+2 >= argc )
346 ERROR_exit("-mrange option requires 2 following arguments!\n");
347 mask_bot = strtod( argv[++narg] , NULL ) ;
348 mask_top = strtod( argv[++narg] , NULL ) ;
349 if( mask_top < mask_top )
350 ERROR_exit("-mrange inputs are illegal!\n") ;
351 narg++ ; continue ;
352 }
353
354 if( strcmp(argv[narg],"-n_rand") == 0 ){
355 if( narg+1 >= argc )
356 ERROR_exit("-n_rand option requires 1 argument!\n");
357 nrand = (int)strtod( argv[++narg] , NULL ) ;
358 narg++ ; continue ;
359 }
360
361 if( strcmp(argv[narg],"-n_randseed") == 0 || strcmp(argv[narg],"-nrandseed") == 0 ){
362 double temp ;
363 if( narg+1 >= argc )
364 ERROR_exit("-n_randseed option requires 1 argument!\n");
365 temp = strtod( argv[++narg] , NULL ) ;
366 nrandseed = temp ;
367 if( nrandseed == 0 ){ /* 29 Mar 2011 */
368 nrandseed = ((unsigned int)time(NULL))%666667u + 17u*(unsigned int)getpid() ;
369 INFO_message("nrandseed initialized to %u",nrandseed) ;
370 }
371 narg++ ; continue ;
372 }
373
374 if( strcmp(argv[narg],"-o") == 0 ){
375 if( narg+1 >= argc )
376 ERROR_exit("-o needs an argument after it!\n");
377 oname = argv[++narg] ;
378 if( ! THD_filename_ok(oname) )
379 ERROR_exit("name after -o is illegal!\n");
380 if( THD_is_file(oname) )
381 ERROR_exit("file %s already exists!\n",oname);
382 narg++ ; continue ;
383 }
384
385 ERROR_message("Unknown option: %s\n",argv[narg]) ;
386 suggest_best_prog_option(argv[0], argv[narg]);
387 exit(1);
388 }
389
390 if( yes_niml ){
391 no_ijk = 1; yes_xyz = 1; yes_index = 0; no_zero = 1; /* 04 Feb 2008 */
392 }
393
394 /* keep -quiet quiet 27 Mar 2006 [rickr] */
395 if( verb > 0 ) PRINT_VERSION("3dmaskdump") ;
396
397 /* should have at least one more argument */
398
399 ndset = argc - narg ;
400 if( ndset <= 0 ) ERROR_exit("No input dataset!?\n") ;
401
402 /* open all input datasets */
403
404 input_dset = (THD_3dim_dataset **) malloc( sizeof(THD_3dim_dataset *)*ndset ) ;
405 for( ndval=ii=0 ; ii < ndset ; ii++ ){
406 input_dset[ii] = THD_open_dataset( argv[narg+ii] ) ;
407 if( input_dset[ii] == NULL )
408 ERROR_exit("Can't open dataset %s!\n",argv[narg+ii]) ;
409
410 if( DSET_BRICK_TYPE(input_dset[ii],0) == MRI_complex )
411 ERROR_exit("Cannot deal with complex-valued input dataset!\n");
412
413 if( ii == 0 ){
414 nvox = DSET_NVOX(input_dset[0]) ;
415 nx = DSET_NX(input_dset[0]) ;
416 ny = DSET_NY(input_dset[0]) ; nxy = nx*ny ;
417 nz = DSET_NZ(input_dset[0]) ; nxyz = nxy*nz ;
418 } else if( DSET_NX(input_dset[ii]) != nx ||
419 DSET_NY(input_dset[ii]) != ny || DSET_NZ(input_dset[ii]) != nz ){
420 ERROR_exit("Dataset %s does not match in size!\n",argv[narg+ii]);
421 }
422
423 ndval += DSET_NVALS(input_dset[ii]) ;
424 }
425
426 /* make a byte mask from mask dataset */
427
428 if( mask_dset == NULL ){
429 mmm = NULL ;
430 if( verb ) INFO_message("%d voxels in the entire dataset (no mask)\n",nvox) ;
431 } else {
432 if( DSET_NVOX(mask_dset) != nvox )
433 ERROR_exit("Input and mask datasets are not same dimensions!\n");
434 mmm = THD_makemask( mask_dset , 0 , mask_bot,mask_top ) ;
435 mcount = THD_countmask( nvox , mmm ) ;
436 if( mcount <= 0 ) ERROR_exit("No voxels in the mask!\n") ;
437 if( verb ) INFO_message("%d voxels in the mask\n",mcount) ;
438 DSET_delete(mask_dset) ;
439 }
440
441 /* 16 Mar 2000: deal with the cmask */
442
443 if( cmask != NULL ){
444 if( ncmask != nvox )
445 ERROR_exit("Input and cmask datasets are not same dimensions!\n");
446 if( mmm != NULL ){
447 for( ii=0 ; ii < nvox ; ii++ ) mmm[ii] = (mmm[ii] && cmask[ii]) ;
448 free(cmask) ;
449 mcount = THD_countmask( nvox , mmm ) ;
450 if( mcount <= 0 ) ERROR_exit("No voxels in the mask+cmask!\n") ;
451 if( verb ) INFO_message("%d voxels in the mask+cmask\n",mcount) ;
452 } else {
453 mmm = cmask ;
454 mcount = THD_countmask( nvox , mmm ) ;
455 if( mcount <= 0 ) ERROR_exit("No voxels in the cmask!\n") ;
456 if( verb ) INFO_message("%d voxels in the cmask\n",mcount) ;
457 }
458 }
459
460 /*-- 09 May 2003: make a mask corresponding to the boxen - RWC --*/
461
462 if( box_num > 0 ){
463 int bb, ibot,itop, jbot,jtop, kbot,ktop , btyp , ii,jj,kk ;
464 float xbot,xtop, ybot,ytop, zbot,ztop ;
465 THD_fvec3 dv,xv ;
466 THD_3dim_dataset *dset = input_dset[0] ;
467 float xmin=dset->daxes->xxmin , xmax=dset->daxes->xxmax ;
468 float ymin=dset->daxes->yymin , ymax=dset->daxes->yymax ;
469 float zmin=dset->daxes->zzmin , zmax=dset->daxes->zzmax ;
470
471 if( bmask == NULL ) bmask = calloc(1,nvox) ;
472
473 for( bb=0 ; bb < box_num ; bb++ ){
474 btyp = box_dat[0+BOXLEN*bb];
475 xbot = box_dat[1+BOXLEN*bb]; xtop = box_dat[2+BOXLEN*bb];
476 ybot = box_dat[3+BOXLEN*bb]; ytop = box_dat[4+BOXLEN*bb];
477 zbot = box_dat[5+BOXLEN*bb]; ztop = box_dat[6+BOXLEN*bb];
478
479 if( btyp != BOX_IJK ){ /* convert coords to indexes */
480
481 if( btyp == BOX_NEU ){ /* coords from Neuroscience to DICOM */
482 xbot = -xbot; xtop = -xtop; ybot = -ybot; ytop = -ytop; btyp = BOX_DIC;
483 }
484 if( btyp == BOX_DIC ){ /* coords from DICOM to dataset */
485 LOAD_FVEC3(dv,xbot,ybot,zbot) ;
486 xv = THD_dicomm_to_3dmm( dset , dv ) ;
487 UNLOAD_FVEC3(xv,xbot,ybot,zbot) ;
488 LOAD_FVEC3(dv,xtop,ytop,ztop) ;
489 xv = THD_dicomm_to_3dmm( dset , dv ) ;
490 UNLOAD_FVEC3(xv,xtop,ytop,ztop) ;
491 }
492 if( xbot < xmin && xtop < xmin ) continue ; /* skip box if outside dataset */
493 if( xbot > xmax && xtop > xmax ) continue ;
494 if( ybot < ymin && ytop < ymin ) continue ;
495 if( ybot > ymax && ytop > ymax ) continue ;
496 if( zbot < zmin && ztop < zmin ) continue ;
497 if( zbot > zmax && ztop > zmax ) continue ;
498 LOAD_FVEC3(dv,xbot,ybot,zbot) ;
499 xv = THD_3dmm_to_3dfind( dset , dv ) ; /* coords from dataset to index */
500 UNLOAD_FVEC3(xv,xbot,ybot,zbot) ;
501 LOAD_FVEC3(dv,xtop,ytop,ztop) ;
502 xv = THD_3dmm_to_3dfind( dset , dv ) ;
503 UNLOAD_FVEC3(xv,xtop,ytop,ztop) ;
504 }
505
506 /* flip as float indices, before truncating [24 Nov 2021 rickr] */
507 if( xbot > xtop ){ dx = xbot; xbot = xtop; xtop = dx; }
508 if( ybot > ytop ){ dx = ybot; ybot = ytop; ytop = dx; }
509 if( zbot > ztop ){ dx = zbot; zbot = ztop; ztop = dx; }
510
511 /* do not round, it could add unrequested voxels [24 Nov 2021 rickr] */
512 ibot = ceilf(xbot) ; jbot = ceilf(ybot) ; kbot = ceilf(zbot) ;
513 itop = floorf(xtop) ; jtop = floorf(ytop) ; ktop = floorf(ztop) ;
514
515 /* skip box if outside dataset */
516 if ( itop < 0 || ibot >= nx ) continue;
517 if ( jtop < 0 || jbot >= ny ) continue;
518 if ( ktop < 0 || kbot >= nz ) continue;
519
520 /* constrain values to dataset dimensions */
521 if ( ibot < 0 ) ibot = 0; if ( itop >= nx ) itop = nx-1;
522 if ( jbot < 0 ) jbot = 0; if ( jtop >= ny ) jtop = ny-1;
523 if ( kbot < 0 ) kbot = 0; if ( ktop >= nz ) ktop = nz-1;
524
525 for( kk=kbot ; kk <= ktop ; kk++ )
526 for( jj=jbot ; jj <= jtop ; jj++ )
527 for( ii=ibot ; ii <= itop ; ii++ ) bmask[ii+jj*nx+kk*nxy] = 1 ;
528
529 } /* end of loop over list of boxes */
530
531 } /* end of box mask */
532
533 /*-- 09 Sep 2009: make a mask corresponding to the balls - RWC --*/
534
535 if( ball_num > 0 ){
536 int bb, ibot,itop, jbot,jtop, kbot,ktop , btyp , ii,jj,kk ;
537 float xcen,ycen,zcen , rad , xx,yy,zz , dist , icen,jcen,kcen ;
538 THD_fvec3 dv,xv ;
539 THD_3dim_dataset *dset = input_dset[0] ;
540 float xmin=dset->daxes->xxmin , xmax=dset->daxes->xxmax ;
541 float ymin=dset->daxes->yymin , ymax=dset->daxes->yymax ;
542 float zmin=dset->daxes->zzmin , zmax=dset->daxes->zzmax ;
543
544 if( bmask == NULL ) bmask = calloc(1,nvox) ; /* will include boxes */
545
546 /* loop over list of balls */
547
548 for( bb=0 ; bb < ball_num ; bb++ ){
549 btyp = ball_dat[0+BOXLEN*bb] ;
550 xcen = ball_dat[1+BOXLEN*bb] ; ycen = ball_dat[2+BOXLEN*bb] ;
551 zcen = ball_dat[3+BOXLEN*bb] ; rad = ball_dat[4+BOXLEN*bb] ;
552
553 /* convert center coords to dataset indexes */
554
555 if( btyp == BALL_NEU ){ /* coords from Neuroscience to DICOM */
556 xcen = -xcen; ycen = -ycen; btyp = BALL_DIC;
557 }
558 if( btyp == BALL_DIC ){ /* coords from DICOM to dataset */
559 LOAD_FVEC3(dv,xcen,ycen,zcen) ;
560 xv = THD_dicomm_to_3dmm( dset , dv ) ;
561 UNLOAD_FVEC3(xv,xcen,ycen,zcen) ;
562 }
563 if( xcen < xmin || xcen > xmax ) continue ; /* skip ball if outside */
564 if( ycen < ymin || ycen > ymax ) continue ;
565 if( zcen < zmin || zcen > zmax ) continue ;
566 LOAD_FVEC3(dv,xcen,ycen,zcen) ;
567 xv = THD_3dmm_to_3dfind( dset , dv ) ; /* coords from dataset to index */
568 UNLOAD_FVEC3(xv,icen,jcen,kcen) ;
569
570 /* create a bounding box around the ball to add */
571 /* scale radius to voxels, and do not round [24 Nov 2021 rickr] */
572 dx = fabs(DSET_DX(dset));
573 dy = fabs(DSET_DY(dset));
574 dz = fabs(DSET_DZ(dset));
575 ibot = ceilf(icen-rad/dx) ; itop = floorf(icen+rad/dx) ;
576 jbot = ceilf(jcen-rad/dy) ; jtop = floorf(jcen+rad/dy) ;
577 kbot = ceilf(kcen-rad/dz) ; ktop = floorf(kcen+rad/dz) ;
578
579 /* skip test box if outside dataset [9 Dec 2021 rickr] */
580 if ( itop < 0 || ibot >= nx ) continue;
581 if ( jtop < 0 || jbot >= ny ) continue;
582 if ( ktop < 0 || kbot >= nz ) continue;
583
584 /* constrain values to dataset dimensions [9 Dec 2021 rickr] */
585 if ( ibot < 0 ) ibot = 0; if ( itop >= nx ) itop = nx-1;
586 if ( jbot < 0 ) jbot = 0; if ( jtop >= ny ) jtop = ny-1;
587 if ( kbot < 0 ) kbot = 0; if ( ktop >= nz ) ktop = nz-1;
588
589 rad = rad*rad ;
590
591 for( kk=kbot ; kk <= ktop ; kk++ ){
592 for( jj=jbot ; jj <= jtop ; jj++ ){
593 for( ii=ibot ; ii <= itop ; ii++ ){
594 LOAD_FVEC3( dv , ii,jj,kk ) ; /* convert to xyz coords */
595 xv = THD_3dfind_to_3dmm( dset , dv ) ; /* then test distance^2 */
596 UNLOAD_FVEC3( xv , xx,yy,zz ) ; /* xyz of ball center. */
597 dist = SQR(xx-xcen) + SQR(yy-ycen) + SQR(zz-zcen) ;
598 if( dist <= rad ) bmask[ii+jj*nx+kk*nxy] = 1 ;
599 }}}
600 } /* end of loop over list of balls */
601
602 } /* end of ball maskology */
603
604 /*--- intersect boxes+balls mask with any other mask we have now ---*/
605
606 if( bmask != NULL ){
607 mcount = THD_countmask( nvox , bmask ) ;
608 if( verb ) INFO_message("%d voxels in the boxes and/or balls\n",mcount) ;
609 if( mcount == 0 )
610 ERROR_exit("Can't continue with no voxels insides boxes+balls!\n");
611 if( mmm != NULL ){
612 for( ii=0 ; ii < nvox ; ii++ )
613 mmm[ii] = (mmm[ii] && bmask[ii]) ; /* intersection */
614 free(bmask) ;
615 mcount = THD_countmask( nvox , mmm ) ;
616 if( mcount <= 0 ) ERROR_exit("No voxels in the mask INTERSECT boxes+balls!\n") ;
617 if( verb ) INFO_message("%d voxels in the mask INTERSECT boxes+balls\n",mcount) ;
618 } else {
619 mmm = bmask ;
620 if( verb ) INFO_message("Using only the boxes+balls mask") ;
621 }
622 }
623
624 /*----- read in input dataset bricks -----*/
625
626 for( ii=0 ; ii < ndset ; ii++ ){
627 DSET_load(input_dset[ii]) ; CHECK_LOAD_ERROR(input_dset[ii]) ;
628 }
629
630 /*--- open the output file ---*/
631
632 if( oname != NULL ){
633 ofile = fopen( oname , "w" ) ;
634 if( ofile == NULL ) ERROR_exit("Can't open output file %s\n",oname);
635 } else {
636 ofile = stdout ;
637 }
638
639 /*--- output string buffers ---*/
640
641 /*-- 09 May 2003: add room for the index (9 -> 10) [rickr] --*/
642 obuf = (char *) malloc( sizeof(char) * (ndval+10) * 16 ) ;
643
644 /* do we want a random set? */
645 if (nrand > 0) {
646 int *iok = (int *)calloc(nvox, sizeof(int));
647 int *ranorder=NULL;
648 int cnt=0 , ccc ;
649
650 cnt = 0;
651 for( ii=0 ; ii < nvox ; ii++ ){
652 if( mmm == NULL || mmm[ii]) {
653 iok[cnt] = ii;
654 ++cnt;
655 }
656 }
657 if( verb )
658 INFO_message("Dumping %d randomly selected voxels \n"
659 "out of %d originally in the mask.\n",
660 ( (cnt<nrand) ? cnt : nrand ), cnt);
661 ccc = cnt-1 ;
662 /* z_rand_order line causes a compiler error on macosx_10.3_G5
663 *
664 * 3dmaskdump.c: In function `main':
665 * 3dmaskdump.c:723: error: unrecognizable insn:
666 * ...
667 * 3dmaskdump.c:723: internal compiler error: in extract_insn,
668 * at recog.c:2175
669 * Old OS and system, no bug report filed. 15 Oct 2009 [rickr]
670 */
671 ranorder = z_rand_order(0, ccc, nrandseed);
672 if (mmm) free(mmm); mmm = (byte *)calloc(nvox, sizeof(byte));
673 for (ii=0; ii < ( (cnt<nrand) ? cnt : nrand ); ++ii) {
674 mmm[iok[ranorder[ii]]] = 1;
675 }
676
677 free(iok); iok = NULL;
678 free(ranorder); ranorder = NULL;
679 }
680
681 /*------- loop over voxels -------*/
682
683 nout = 0 ; bout = NULL ;
684 for( ii=0 ; ii < nvox ; ii++ ){
685 if( mmm != NULL && mmm[ii] == 0 ) continue ; /* skip this voxel */
686
687 obuf[0] = '\0' ; /* start it off with nothing */
688
689 /*-- 09 May 2003: optionally display index [rickr] --*/
690 if ( yes_index ){
691 otemp = MV_format_fval((float)ii);
692 strcat(obuf,otemp); strcat(obuf," ");
693 }
694
695 if( ! no_ijk ){
696 i = DSET_index_to_ix( input_dset[0] , ii ) ; /* voxel indexes */
697 j = DSET_index_to_jy( input_dset[0] , ii ) ;
698 k = DSET_index_to_kz( input_dset[0] , ii ) ;
699
700 otemp = MV_format_fval((float)i); strcat(obuf,otemp); strcat(obuf," ");
701 otemp = MV_format_fval((float)j); strcat(obuf,otemp); strcat(obuf," ");
702 otemp = MV_format_fval((float)k); strcat(obuf,otemp); strcat(obuf," ");
703 }
704
705 if( yes_xyz ){ /* 23 Mar 2003 */
706 THD_ivec3 ind ; THD_fvec3 vec ;
707 i = DSET_index_to_ix( input_dset[0] , ii ) ;
708 j = DSET_index_to_jy( input_dset[0] , ii ) ;
709 k = DSET_index_to_kz( input_dset[0] , ii ) ;
710 LOAD_IVEC3(ind,i,j,k) ;
711 vec = THD_3dind_to_3dmm ( input_dset[0] , ind ) ;
712 vec = THD_3dmm_to_dicomm( input_dset[0] , vec ) ;
713 if(lpi_xyz) { /* for LPI output, negate RAI's x and y */
714 vec.xyz[0] = - vec.xyz[0]; vec.xyz[1] = -vec.xyz[1];
715 }
716 otemp = MV_format_fval(vec.xyz[0]); strcat(obuf,otemp); strcat(obuf," ");
717 otemp = MV_format_fval(vec.xyz[1]); strcat(obuf,otemp); strcat(obuf," ");
718 otemp = MV_format_fval(vec.xyz[2]); strcat(obuf,otemp); strcat(obuf," ");
719 }
720
721 for( numz=kk=0 ; kk < ndset ; kk++ ){
722 flim = THD_extract_series( ii , input_dset[kk] , 0 ) ;
723 flar = MRI_FLOAT_PTR(flim) ;
724 for( vv=0 ; vv < flim->nx ; vv++ ){
725 if( flar[vv] == 0.0f ) numz++ ; /* 04 Feb 2008 */
726 otemp = MV_format_fval(flar[vv]) ;
727 strcat(obuf,otemp) ; strcat(obuf," ") ;
728 }
729 mri_free(flim) ;
730 }
731 if( no_zero && numz == ndval ) continue ; /* 04 Feb 2008 */
732
733 jj = strlen(obuf) ; obuf[jj-1] = '\0' ; /* kill last blank */
734
735 if( yes_niml ){ /* 04 Feb 2008: save output for NIML-ization */
736 nout++ ;
737 bout = (char **)realloc(bout,sizeof(char *)*nout) ;
738 bout[nout-1] = strdup(obuf) ;
739 } else {
740 fprintf(ofile,"%s\n",obuf) ; /* regular output */
741 }
742
743 } /*--- end of loop over voxels ---*/
744
745 if( yes_niml ){ /*----- 04 Feb 2008: NIML formatted output -----*/
746 if( nout > 0 ){
747 char *gstr = EDIT_get_geometry_string(input_dset[0]) ;
748 fprintf(ofile,"<VOLUME_DATA_SPARSE\n") ;
749 if( gstr != NULL && *gstr != '\0' )
750 fprintf(ofile," geometry_string='%s'\n",gstr) ;
751 fprintf(ofile," target_name='%s'\n",niml_name) ;
752 fprintf(ofile," ni_type='%d*float'\n",ndval+3) ; /* +3 allows for xyz */
753 fprintf(ofile," ni_dimen='%d'\n",nout) ;
754 fprintf(ofile,">\n") ;
755 for( ii=0 ; ii < nout ; ii++ ){
756 fprintf(ofile,"%s\n",bout[ii]); free(bout[ii]);
757 }
758 fprintf(ofile,"</VOLUME_DATA_SPARSE>\n") ;
759 free(bout) ;
760 if( verb ) INFO_message("-niml: output %d voxels",nout) ;
761 } else {
762 ERROR_exit("-niml output requested, but no data to output!") ;
763 }
764 }
765
766 /*----- Free at last! -----*/
767
768 fflush(ofile); fsync(fileno(ofile)); fclose(ofile); exit(0) ;
769 }
770