1 /*
2
3 ver = 1.0; date = Dec 1, 2021
4 + [PT] start of this program
5
6 ver = 1.2; date = Dec 3, 2021
7 + [PT] getting there, for basic functionality. Defaults work pretty well.
8
9 ver = 1.3; date = Dec 6, 2021
10 + [PT] start to add scaling of edge values...
11
12 ver = 1.31; date = Dec 6, 2021
13 + [PT] ... and scaling now activated: using DOG value at each edge vox
14
15 ver = 1.4; date = Dec 9, 2021
16 + [PT] add -only2D opt, for DRG
17
18 ver = 1.5; date = Dec 23, 2021
19 + [PT] run with binarized/optimized EulerDist form. shd be faster
20
21 ver = 1.51; date = Dec 26, 2021
22 + [PT] all output dsets get their history full made
23
24 ver = 1.52; date = Dec 26, 2021
25 + [PT] start adding in pieces for masking---have NOT applied yet
26
27 ver = 1.6; date = Dec 26, 2021
28 + [PT] add in -automask, automask+X, and '-mask ..' option behavior
29
30
31 *** still might add:
32 - multi-spatial scale requests
33
34 */
35
36 #include <stdio.h>
37 #include <stdlib.h>
38 #include <math.h>
39 #include <unistd.h>
40 #include "debugtrace.h"
41 #include "mrilib.h"
42
43 int run_edge_dog( int comline, PARAMS_edge_dog opts,
44 int argc, char *argv[] );
45
46 // --------------------------------------------------------------------------
47
usage_3dedgedog()48 int usage_3dedgedog()
49 {
50 char *author = "PA Taylor and DR Glen (SSCC, NIMH, NIH)";
51 PARAMS_edge_dog opts;
52
53 opts = set_edge_dog_defaults();
54
55 printf(
56 "\n"
57 "Overview ~1~ \n"
58 "\n"
59 "This program calculates edges in an image using the Difference of Gaussians\n"
60 "(DOG) method by Wilson and Giese (1977) and later combined with work by\n"
61 "Marr and Hildreth (1980) to provide a computationally efficient\n"
62 "approximation to their Lagrangian of Gaussian (LOG) method for calculating\n"
63 "edges in an image. This is a fascinating set of papers to read. But you\n"
64 "don't have to take *my* word for it!...\n"
65 "\n"
66 "Generating edges in this way has some interesting properties, such as\n"
67 "numerical efficiency and edges that are closed loops/surfaces. The edges\n"
68 "can be tuned to focus on structures of a particular size, too, which can be\n"
69 "particularly useful in some applications.\n"
70 "\n"
71 "written by: %s\n"
72 "\n"
73 "Description ~2~ \n"
74 "\n"
75 "The primary papers for learning more about the DOG and LOG methods are:\n"
76 "\n"
77 " Wilson HR, Giese SC (1977). Threshold visibility of frequency\n"
78 " gradient patterns. Vision Res. 17(10):1177-90.\n"
79 " doi: 10.1016/0042-6989(77)90152-3. PMID: 595381.\n"
80 "\n"
81 " Marr D, Hildreth E (1980). Theory of edge detection. Proc R Soc\n"
82 " Lond B Biol Sci. 207(1167):187-217.\n"
83 " doi: 10.1098/rspb.1980.0020. PMID: 6102765.\n"
84 "\n"
85 "Thanks to C. Rorden for pointing these papers out and discussing them.\n"
86 "\n"
87 "The current code here extends/tweaks the MH1980 algorithm a bit. It runs\n"
88 "in 3D by default (a straightforward extension), it also employs the\n"
89 "Eulerian Distance Transform (EDT) to pick out the actual edges from the\n"
90 "DOG step---see 3dEulerDist for more information about the EDT.\n"
91 "\n"
92 "The DOG-based edges require specifying a couple parameters, the main\n"
93 "one being interpretable as a minimal 'scale size' for structures. In this\n"
94 "code, this is the 'sigma_rad' (or 'sigma_nvox', if you want to specify it\n"
95 "in terms of the number of voxels along a given axis), which is the 'inner\n"
96 "Gaussian' sigma value, if you are following MH1980. The default for this\n"
97 "sigma_rad parameter is set based on the expected average thickness of adult\n"
98 "human GM, but it is easily alterable at the command line for any other\n"
99 "values.\n"
100 "\n"
101 "==========================================================================\n"
102 "\n"
103 "Command usage and option list ~1~ \n"
104 "\n"
105 " 3dedgedog [options] -prefix PREF -input DSET\n"
106 "\n"
107 "where: \n"
108 "\n"
109 " -input DSET :(req) input dataset\n"
110 "\n"
111 " -prefix PREF :(req) output prefix name\n"
112 "\n"
113 " -mask MASK :mask dataset. NB: this mask is only applied *after*\n"
114 " the EDT has been calculated. Therefore, the boundaries\n"
115 " of this mask have no affect on the calculated distance\n"
116 " values, except for potentially zeroing some out at the\n"
117 " end. Mask only gets made from [0]th vol.\n"
118 "\n"
119 " -automask :alternative to '-mask ..', for automatic internal\n"
120 " calculation of a mask in the usual AFNI way. Again, this\n"
121 " mask is only applied after all calcs (so using this does\n"
122 " not speed up the calc or affect distance values).\n"
123 " ** Special note: you can also write '-automask+X', where\n"
124 " X is some integer; this will dilate the initial automask\n"
125 " X number of times (as in 3dAllineate); must have X>0.\n"
126 "\n"
127 " -sigma_rad RRR :radius for 'inner' Gaussian, in units of mm; RRR must\n"
128 " by greater than zero (def: %.2f). Default is chosen to\n"
129 " approximate typical GM thickness in human adults.\n"
130 "\n"
131 " -sigma_nvox NNN :define radius for 'inner' Gaussian by providing a\n"
132 " multiplicative factor for voxel edge length, which will\n"
133 " be applied in each direction; NNN can be any float\n"
134 " greater than zero. This is an alternative to the\n"
135 " '-sigma_rad ..' opt (def: use '-sigma_rad' and its\n"
136 " default value).\n"
137 "\n"
138 " -ratio_sigma RS :the ratio of inner and outer Gaussian sigma values.\n"
139 " That is, RS defines the size of the outer Gaussian,\n"
140 " by scaling up the inner value. RS can be any float\n"
141 " greater than 1 (def: %.2f). See 'Notes' for more about\n"
142 " this parameter.\n"
143 "\n"
144 " -output_intermed :use this option flag if you would like to output some\n"
145 " intermediate dataset(s):\n"
146 " + DOG (difference of Gaussian)\n"
147 " + EDT2 (Euler Distance Transform, dist**2 values),\n"
148 " [0]th vol only\n"
149 " (def: not output). Output names will be user-entered\n"
150 " prefix with a representative suffix appended.\n"
151 "\n"
152 " -edge_bnd_NN EBN :specify the 'nearest neighbor' (NN) value for the\n"
153 " connectedness of the drawn boundaries. EBN must be\n"
154 " one of the following integer values:\n"
155 " 1 -> for face only\n"
156 " 2 -> for face+edge\n"
157 " 3 -> for face+edge+node\n"
158 " (def: %d).\n"
159 "\n"
160 " -edge_bnd_side EBS :specify which boundary layer around the zero-layer\n"
161 " to use in the algorithm. EBS must be one of the\n"
162 " following keywords:\n"
163 " \"NEG\" -> for negative (inner) boundary\n"
164 " \"POS\" -> for positive (outer) boundary\n"
165 " \"BOTH\" -> for both (inner+outer) boundary\n"
166 " \"BOTH_SIGN\" -> for both (inner+outer) boundary,\n"
167 " with pos/neg sides keeping sign\n"
168 " (def: \"%s\").\n"
169 "\n"
170 " -edge_bnd_scale :by default, this program outputs a mask of edges, so\n"
171 " edge locations have value=1, and everything else is 0.\n"
172 " Using this option means the edges will have values\n"
173 " scaled to have a relative magnitude between 0 and 100\n"
174 " (NB: the output dset will still be datum=short)\n"
175 " depending on the gradient value at the edge.\n"
176 "\n"
177 " -only2D SLI :instead of estimating edges in full 3D volume, calculate\n"
178 " edges just in 2D, per plane. Provide the slice plane\n"
179 " you want to run along as the single argument SLI:\n"
180 " \"axi\" -> for axial slice\n"
181 " \"cor\" -> for coronal slice\n"
182 " \"sag\" -> for sagittal slice\n"
183 "\n"
184 "==========================================================================\n"
185 "\n"
186 "Notes ~1~\n"
187 "\n"
188 "The value of sigma_rad ~2~\n"
189 "\n"
190 "(... which sounds like the title of a great story, no? Anyways...)\n"
191 "This parameter represents the ratio of the width of the two Gaussians that\n"
192 "are blurred in the first stage of the DOG estimation. In the limit that\n"
193 "sigma_rad approaches 1, the DOG -> LOG. So, we want to keep the value of\n"
194 "this parameter in the general vicinity of 1 (and it can't be less than 1,\n"
195 "because the ratio is of the outer-to-the-inner Gaussian). MH1980 suggested\n"
196 "that sigma_rad=1.6 was optimal 'on engineering grounds' of bandwidth\n"
197 "sensitivity of filters. This is *very approximate* reasoning, but provides\n"
198 "another reference datum for selection.\n"
199 "\n"
200 "Because the DOG approximation used here is for visual purposes of MRI\n"
201 "datasets, often even more specifically for alignment purposes, we have\n"
202 "chosen a default value that seemed visually appropriate to real data.\n"
203 "Values of sigma_rad close to one show much noisier, scattered images---that\n"
204 "is, they pick up *lots* of contrast differences, probably too many for most\n"
205 "visualization purposes. Edge images smoothen as sigma_rad increases, but\n"
206 "as it gets larger, it can also blend together edges of features---such as\n"
207 "gyri of the brain with dura. So, long story short, the default value here\n"
208 "tries to pick a reasonable middle ground.\n"
209 "\n"
210 "==========================================================================\n"
211 "\n"
212 "Examples ~1~\n"
213 "\n"
214 "1) Basic case:\n"
215 " 3dedgedog \\\n"
216 " -input anat+orig.HEAD \\\n"
217 " -prefix anat_EDGE.nii.gz \n"
218 "\n"
219 "2) Same as above, but output both edges from the DOG+EDT steps, keeping\n"
220 " the sign of each side:\n"
221 " 3dedgedog \\\n"
222 " -edge_bnd_side BOTH_SIGN \\\n"
223 " -input anat+orig.HEAD \\\n"
224 " -prefix anat_EDGE_BOTHS.nii.gz \n"
225 "\n"
226 "3) Output both sides of edges, and scale the edge values (by DOG value):\n"
227 " 3dedgedog \\\n"
228 " -edge_bnd_side BOTH_SIGN \\\n"
229 " -edge_bnd_scale \\\n"
230 " -input anat+orig.HEAD \\\n"
231 " -prefix anat_EDGE_BOTHS_SCALE.nii.gz \n"
232 "\n"
233 "4) Increase scale size of edged shapes to 2.7mm:\n"
234 " 3dedgedog \\\n"
235 " -sigma_rad 2.7 \\\n"
236 " -edge_bnd_scale \\\n"
237 " -input anat+orig.HEAD \\\n"
238 " -prefix anat_EDGE_BOTHS_SCALE.nii.gz \n"
239 "\n"
240 "5) Apply automasking, with a bit of mask dilation so outer boundary is\n"
241 " included:\n"
242 " 3dedgedog \\\n"
243 " -automask+2 \\\n"
244 " -input anat+orig.HEAD \\\n"
245 " -prefix anat_EDGE_AMASK.nii.gz \n"
246 "\n"
247 "==========================================================================\n"
248 "\n",
249 author, opts.sigma_rad[0], opts.ratio_sigma, opts.edge_bnd_NN,
250 opts.edge_bnd_side_user );
251
252 return 0;
253 }
254
main(int argc,char * argv[])255 int main(int argc, char *argv[]) {
256
257 int ii = 0;
258 int iarg;
259 PARAMS_edge_dog InOpts;
260 float tmp;
261 int itmp;
262 char *auto_tstring = NULL;
263
264 mainENTRY("3dedgedog"); machdep();
265
266 // fill option struct with defaults
267 InOpts = set_edge_dog_defaults();
268
269 // ****************************************************************
270 // parse command line arguments
271 // ****************************************************************
272
273 if (argc == 1) { usage_3dedgedog(); exit(0); }
274
275 /* scan through args */
276 iarg = 1;
277 while( iarg < argc && argv[iarg][0] == '-' ){
278
279 if( strcmp(argv[iarg],"-help") == 0 ||
280 strcmp(argv[iarg],"-h") == 0 ) {
281 usage_3dedgedog();
282 exit(0);
283 }
284
285 if( strcmp(argv[iarg],"-input") == 0 ){
286 if( ++iarg >= argc )
287 ERROR_exit("Need argument after '%s'", argv[iarg-1]);
288
289 InOpts.input_name = strdup(argv[iarg]) ;
290 iarg++ ; continue ;
291 }
292
293 if( strcmp(argv[iarg],"-mask") == 0 ){
294 if( ++iarg >= argc )
295 ERROR_exit("Need argument after '%s'", argv[iarg-1]);
296
297 InOpts.mask_name = strdup(argv[iarg]) ;
298 iarg++ ; continue ;
299 }
300
301 if( strncmp(argv[iarg],"-automask", 9) == 0 ) {
302 InOpts.do_automask = 1;
303
304 auto_tstring = argv[iarg];
305 if ( auto_tstring[9] == '+' && auto_tstring[10] != '\0' ) {
306 InOpts.amask_ndil = (int) strtod(auto_tstring+10, NULL);
307 if ( InOpts.amask_ndil <= 0 )
308 ERROR_exit("Cannot have number of dilations (%d) <=0.",
309 InOpts.amask_ndil);
310 }
311
312 iarg++ ; continue ;
313 }
314
315 if( strcmp(argv[iarg],"-prefix") == 0 ){
316 if( ++iarg >= argc )
317 ERROR_exit("Need argument after '%s'", argv[iarg-1]);
318
319 InOpts.prefix = strdup(argv[iarg]) ;
320
321 if( !THD_filename_ok(InOpts.prefix) )
322 ERROR_exit("Illegal name after '%s'", argv[iarg-1]);
323 iarg++ ; continue ;
324 }
325
326 if( strcmp(argv[iarg],"-sigma_rad") == 0) {
327 if( ++iarg >= argc )
328 ERROR_exit("Need argument after '%s'", argv[iarg-1]);
329
330 tmp = atof(argv[iarg]);
331 if( tmp <= 0 )
332 ERROR_exit("Need positive value after '%s'", argv[iarg-1]);
333
334 InOpts.sigma_rad[0] = tmp;
335 InOpts.sigma_rad[1] = tmp;
336 InOpts.sigma_rad[2] = tmp;
337
338 iarg++ ; continue ;
339 }
340
341 if( strcmp(argv[iarg],"-sigma_nvox") == 0) {
342 if( ++iarg >= argc )
343 ERROR_exit("Need argument after '%s'", argv[iarg-1]);
344
345 tmp = atof(argv[iarg]);
346 if( tmp <= 0 )
347 ERROR_exit("Need positive value after '%s'", argv[iarg-1]);
348
349 InOpts.sigma_nvox[0] = tmp;
350 InOpts.sigma_nvox[1] = tmp;
351 InOpts.sigma_nvox[2] = tmp;
352
353 iarg++ ; continue ;
354 }
355
356 if( strcmp(argv[iarg],"-edge_bnd_NN") == 0) {
357 if( ++iarg >= argc )
358 ERROR_exit("Need argument after '%s'", argv[iarg-1]);
359
360 itmp = atoi(argv[iarg]);
361 if( 1 <= itmp && itmp <=3 )
362 InOpts.edge_bnd_NN = itmp;
363 else
364 ERROR_exit("Need either 1, 2 or 3 after '%s'", argv[iarg-1]);
365
366 iarg++ ; continue ;
367 }
368
369 if( strcmp(argv[iarg],"-edge_bnd_side") == 0) {
370 if( ++iarg >= argc )
371 ERROR_exit("Need argument after '%s'", argv[iarg-1]);
372
373 if( strcmp(argv[iarg],"NEG") == 0)
374 InOpts.edge_bnd_side = -1;
375 else if( strcmp(argv[iarg],"POS") == 0)
376 InOpts.edge_bnd_side = 1;
377 else if( strcmp(argv[iarg],"BOTH") == 0)
378 InOpts.edge_bnd_side = 2;
379 else if( strcmp(argv[iarg],"BOTH_SIGN") == 0)
380 InOpts.edge_bnd_side = 3;
381 else
382 ERROR_exit("Need either \"NEG\", \"POS\" or \"BOTH\" "
383 "after '%s'", argv[iarg-1]);
384
385 iarg++ ; continue ;
386 }
387
388 if( strcmp(argv[iarg],"-edge_bnd_scale") == 0 ){
389 InOpts.edge_bnd_scale = 1;
390 iarg++ ; continue ;
391 }
392
393 if( strcmp(argv[iarg],"-ratio_sigma") == 0) {
394 if( ++iarg >= argc )
395 ERROR_exit("Need argument after '%s'", argv[iarg-1]);
396
397 tmp = atof(argv[iarg]);
398 if( tmp <= 1.0 )
399 ERROR_exit("Need value >1 after '%s'", argv[iarg-1]);
400
401 InOpts.ratio_sigma = tmp;
402
403 iarg++ ; continue ;
404 }
405
406 if( strcmp(argv[iarg],"-output_intermed") == 0 ){
407 InOpts.do_output_intermed = 1;
408 iarg++ ; continue ;
409 }
410
411 // same as in 3dEulerDist
412 if( strcmp(argv[iarg],"-only2D") == 0) {
413 if( ++iarg >= argc )
414 ERROR_exit("Need argument after '%s'", argv[iarg-1]);
415
416 if( strcmp(argv[iarg],"cor") == 0 || \
417 strcmp(argv[iarg],"axi") == 0 || \
418 strcmp(argv[iarg],"sag") == 0 )
419 InOpts.only2D = strdup(argv[iarg]);
420 else
421 ERROR_exit("Need either \"cor\", \"axi\" or \"sag\" "
422 "after '%s'", argv[iarg-1]);
423
424 iarg++ ; continue ;
425 }
426
427 ERROR_message("Bad option '%s'\n",argv[iarg]);
428 suggest_best_prog_option(argv[0], argv[iarg]);
429 exit(1);
430 }
431
432 // ****************************************************************
433 // verify presence+behavior of inputs
434 // ****************************************************************
435
436 INFO_message("3dedgedog: verify inputs");
437
438 if ( !InOpts.input_name ) {
439 ERROR_message("You need to provide an input dset with '-input ..'");
440 exit(1);
441 }
442
443 if ( !InOpts.prefix )
444 ERROR_exit("Need an output name via '-prefix ..'\n");
445
446 if ( InOpts.mask_name && InOpts.do_automask )
447 ERROR_exit("Cannot combine '-mask ..' and '-automask'. "
448 "You must choose which you *really* want.\n");
449
450 // DONE FILLING, now do the work
451 ii = run_edge_dog(1, InOpts, argc, argv);
452
453 return 0;
454 }
455
456
run_edge_dog(int comline,PARAMS_edge_dog opts,int argc,char * argv[])457 int run_edge_dog( int comline, PARAMS_edge_dog opts,
458 int argc, char *argv[] )
459 {
460 int i, j, k, idx;
461 int nn;
462 int nx, ny, nz, nxy, nvox, nvals;
463
464 THD_3dim_dataset *dset_input = NULL; // input
465 THD_3dim_dataset *dset_mask = NULL; // mask
466 THD_3dim_dataset *dset_dog = NULL; // intermed/out
467 THD_3dim_dataset *dset_bnd = NULL; // output
468
469 char prefix_dog[THD_MAX_PREFIX];
470
471 byte *mask_arr = NULL; // what mask becomes
472 int nmask = -1;
473
474 ENTRY("run_edge_dog");
475
476 dset_input = THD_open_dataset(opts.input_name);
477 if( dset_input == NULL )
478 ERROR_exit("Can't open dataset '%s'", opts.input_name);
479 DSET_load(dset_input); CHECK_LOAD_ERROR(dset_input);
480
481 nx = DSET_NX(dset_input);
482 ny = DSET_NY(dset_input);
483 nz = DSET_NZ(dset_input);
484 nxy = nx*ny;
485 nvox = DSET_NVOX(dset_input);
486 nvals = DSET_NVALS(dset_input);
487
488 // get/make mask array, if user asks
489 if( opts.mask_name ) {
490 dset_mask = THD_open_dataset(opts.mask_name);
491 if( dset_mask == NULL )
492 ERROR_exit("Can't open dataset '%s'", opts.mask_name);
493 DSET_load(dset_mask); CHECK_LOAD_ERROR(dset_mask);
494
495 if( THD_dataset_mismatch( dset_input, dset_mask ) )
496 ERROR_exit("Mismatch between input and mask dsets!\n");
497
498 // mask comes from just [0]th vol
499 mask_arr = THD_makemask( dset_mask, 0, 0, -1 );
500 nmask = THD_countmask( nvox, mask_arr );
501 if( opts.verb )
502 INFO_message("Number of voxels in automask = %d / %d", nmask, nvox);
503 if( nmask < 1 )
504 ERROR_exit("Automask is too small to process");
505
506 DSET_delete(dset_mask);
507 if(dset_mask)
508 free(dset_mask);
509 }
510 else if( opts.do_automask ) {
511 mask_arr = THD_automask( dset_input );
512 if( mask_arr == NULL )
513 ERROR_exit("Can't create -automask from input dataset?");
514
515 if( opts.amask_ndil ){
516 if( opts.verb )
517 INFO_message("Inflate automask %d times", opts.amask_ndil);
518
519 // a very odd recipe for inflating a mask... mirroring
520 // default 3dAllineate behavior (slightly different than that
521 // of 3dAutomask; commented out)
522 for( i=0 ; i < opts.amask_ndil ; i++ ){
523 j = THD_mask_dilate( nx, ny, nz, mask_arr, 3, 2 );
524 j = THD_mask_fillin_once( nx, ny, nz, mask_arr, 2 );
525 //j = THD_mask_fillin_completely( nx, ny, nz , mask_arr, 1 );
526 }
527 }
528
529 nmask = THD_countmask( nvox, mask_arr );
530 if( opts.verb )
531 INFO_message("Number of voxels in automask = %d / %d", nmask, nvox);
532 if( nmask < 1 )
533 ERROR_exit("Automask is too small to process");
534 }
535
536 // if only running in 2D, figure out which slice that is;
537 // same as in 3dEulerDist
538 if ( opts.only2D )
539 i = choose_axes_for_plane( dset_input, opts.only2D,
540 opts.axes_to_proc, opts.verb );
541
542 /* Prepare header for output by copying that of input, and then
543 changing items as necessary */
544 dset_dog = EDIT_empty_copy( dset_input );
545 i = build_edge_dog_suppl_prefix( &opts, prefix_dog, "_DOG" );
546 EDIT_dset_items(dset_dog,
547 ADN_nvals, nvals,
548 ADN_datum_all, MRI_float,
549 ADN_prefix, prefix_dog,
550 ADN_none );
551
552 // calculate DOG
553 INFO_message("Calculate DOG");
554 for( nn=0 ; nn<nvals ; nn++ )
555 i = calc_edge_dog_DOG(dset_dog, opts, dset_input, nn);
556
557 // Output DOG data, if asked
558 if ( opts.do_output_intermed ){
559 INFO_message("Output intermediate dset: %s", prefix_dog);
560
561 THD_load_statistics( dset_dog );
562 if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(dset_dog)) )
563 ERROR_exit("Can't overwrite existing dataset '%s'",
564 DSET_HEADNAME(dset_dog));
565 tross_Make_History("3dedgedog", argc, argv, dset_dog);
566
567 // write dset
568 THD_write_3dim_dataset(NULL, NULL, dset_dog, True);
569 }
570
571 // ------------------------ calc edge/bnd ---------------------------
572
573 dset_bnd = EDIT_empty_copy( dset_input );
574 EDIT_dset_items(dset_bnd,
575 ADN_nvals, nvals,
576 ADN_datum_all, MRI_short,
577 ADN_prefix, opts.prefix,
578 ADN_none );
579
580 INFO_message("Calculate boundaries");
581
582 // calculate edges/bnds: might be several ways to these from DOG data
583 for( nn=0 ; nn<nvals ; nn++ ){
584 i = calc_edge_dog_BND(dset_bnd, opts, dset_dog, nn, argc, argv);
585
586 if( opts.edge_bnd_scale ){
587 if( !nn )
588 INFO_message("Scale boundaries");
589 i = scale_edge_dog_BND(dset_bnd, opts, dset_dog, nn);
590 }
591 }
592
593 if ( mask_arr ) {
594 if( opts.verb )
595 INFO_message("Apply mask");
596 THD_3dim_dataset *dset_tmp = NULL;
597
598 dset_tmp = thd_apply_mask( dset_bnd, mask_arr,
599 opts.prefix );
600 DSET_delete( dset_bnd );
601 dset_bnd = dset_tmp;
602 dset_tmp = NULL;
603 }
604
605 INFO_message("Output main dset: %s", opts.prefix);
606
607 THD_load_statistics( dset_bnd );
608 if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(dset_bnd)) )
609 ERROR_exit("Can't overwrite existing dataset '%s'",
610 DSET_HEADNAME(dset_bnd));
611 tross_Make_History("3dedgedog", argc, argv, dset_bnd);
612
613 // write edge/bnd
614 THD_write_3dim_dataset(NULL, NULL, dset_bnd, True);
615
616 // ------ write out automask, if used and intermediates are output ------
617 if( opts.do_output_intermed && opts.do_automask ){
618 char prefix_amask[THD_MAX_PREFIX];
619 THD_3dim_dataset *dset_amask = NULL;
620
621 dset_amask = EDIT_empty_copy( dset_input );
622 i = build_edge_dog_suppl_prefix( &opts, prefix_amask, "_AMASK" );
623
624 if( opts.verb )
625 INFO_message("Write out automask dset: %s", prefix_amask);
626
627 EDIT_dset_items(dset_amask,
628 ADN_nvals, 1,
629 ADN_datum_all, MRI_byte,
630 ADN_prefix, prefix_amask,
631 ADN_none );
632
633 EDIT_substitute_brick(dset_amask, 0, MRI_byte, mask_arr);
634
635 THD_load_statistics( dset_amask );
636 if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(dset_amask)) )
637 ERROR_exit("Can't overwrite existing dataset '%s'",
638 DSET_HEADNAME(dset_amask));
639 tross_Make_History("3dedgedog", argc, argv, dset_amask);
640
641 // write and free dset
642 THD_write_3dim_dataset(NULL, NULL, dset_amask, True);
643 DSET_delete(dset_amask);
644 free(dset_amask);
645 mask_arr = NULL;
646 }
647
648
649 // ------------------------- free stuff ---------------------------
650
651 if( dset_input ){
652 DSET_delete(dset_input);
653 free(dset_input);
654 }
655
656 if( dset_dog ){
657 DSET_delete(dset_dog);
658 free(dset_dog);
659 }
660
661 if( dset_bnd ){
662 DSET_delete(dset_bnd);
663 free(dset_bnd);
664 }
665
666 if( mask_arr )
667 free(mask_arr);
668
669 return 0;
670 }
671
672