1 /*
2 Working from P. Lauren's distanceField.c, which was made in
3 connection with PA Taylor's Python library for this.
4
5 ver = 2.0; date = Nov 29, 2021
6 + [PT] this program has been a longtime coming. This version merges
7 P Taylor's Python version in lib_EDT.py with P Lauren's concurrent
8 work on a C version (which had been compared/developed in part with
9 the aformentioned lib_EDT.py).
10
11 ver = 2.1; date = Dec 1, 2021
12 + [PT] Rearrange how main run_EDT_3D works, and how calc_EDT_3D() is
13 called. Mostly, we want to make it easier for other programs to use
14 the latter func
15
16 ver = 2.2; date = Dec 1, 2021
17 + [PT] Bug fix: had delta set incorrectly in 2 out of 3
18 calc_EDT_3d_dim?() funcs, because was using *DZ* everywhere.
19
20 ver = 2.3; date = Dec 8, 2021
21 + [PT] new opt, '-only2D ..' so that EDT can be calced in only 2D, if
22 desired. For DRG, may he use it well.
23
24 ver = 2.4; date = Dec 9, 2021
25 + [PT] internal tweak: allocate some tmp arrays less, doesn't really
26 make much difference in speed :(
27
28 ver = 2.5; date = Dec 9, 2021
29 + [PT] fix 2D selection---was only correct for some dset orientations
30
31 ver = 2.6; date = Dec 21, 2021
32 + [PT] change run_EDTD_per_line to use a working array
33 - also start adding in -binary_only opt (not being used yet)
34
35 ver = 2.7; date = Dec 23, 2021
36 + [PT] -binary_only opt now working well (faster code running for special
37 case of binary input mask)
38
39 ver = 2.71; date = Dec 24, 2021
40 + [PT] minor efficiency increase in some index calcs
41
42 ver = 2.72; date = Dec 26, 2021
43 + [PT] fix help to give correct option name: -bounds_are_not_zero
44 - who knew that missing out 'not' could change meaning so much?
45
46 */
47
48 #include <stdio.h>
49 #include <stdlib.h>
50 #include <math.h>
51 #include <unistd.h>
52 #include "debugtrace.h"
53 #include "mrilib.h"
54
55 int run_EDT_3D( int comline, PARAMS_euler_dist opts,
56 int argc, char *argv[] );
57
58 // --------------------------------------------------------------------------
59
usage_3dEulerDist()60 int usage_3dEulerDist()
61 {
62 char *author = "PA Taylor and P Lauren (SSCC, NIMH, NIH)";
63
64 printf(
65 "\n"
66 "Overview ~1~ \n"
67 "\n"
68 "This program calculates the Eulerian Distance Transform (EDT).\n"
69 "\n"
70 "Basically, this means calculating the Euclidean distance of each\n"
71 "voxel's centroid to the nearest boundary with a separate ROI (well, to be\n"
72 "brutally technical, to centroid of the nearest voxel in a neighboring ROI.\n"
73 "The input dataset should be a map of ROIs (so, integer-valued). The\n"
74 "EDT values are calculated throughout the entire FOV by default,\n"
75 "even in the zero/background regions (there is an option to control this).\n"
76 "\n"
77 "written by: %s\n"
78 "\n"
79 "Description ~2~ \n"
80 "\n"
81 "This code calculates the Euclidean Distance Transform (EDT) for 3D\n"
82 "volumes following this nice, efficient algorithm, by Felzenszwalb\n"
83 "and Huttenlocher (2012; FH2012):\n"
84 "\n"
85 " Felzenszwalb PF, Huttenlocher DP (2012). Distance Transforms of\n"
86 " Sampled Functions. Theory of Computing 8:415-428.\n"
87 " https://cs.brown.edu/people/pfelzens/papers/dt-final.pdf\n"
88 "\n"
89 "Thanks to C. Rorden for pointing this paper out and discussing it.\n"
90 "\n"
91 "The current code here extends/tweaks the FH2012 algorithm to a more\n"
92 "general case of having several different ROIs present, for running\n"
93 "in 3D (trivial extension), and for having voxels of non-unity and\n"
94 "non-isotropic lengths. It does this by utilizing the fact that at\n"
95 "its very heart, the FH2012 algorithm works line by line and can even\n"
96 "be thought of as working boundary-by-boundary.\n"
97 "\n"
98 "Here, the zero-valued 'background' is also just treated like an ROI,\n"
99 "with one difference. At a FOV boundary, the zero-valued\n"
100 "ROI/backgroud is treated as open, so that the EDT value at each\n"
101 "'zero' voxel is always to one of the shapes within the FOV. For\n"
102 "nonzero ROIs, one can treat the FOV boundary *either* as an ROI edge\n"
103 "(EDT value there will be 1 edge length) *or* as being open.\n"
104 "\n"
105 "==========================================================================\n"
106 "\n"
107 "Command usage and option list ~1~ \n"
108 "\n"
109 " 3dEulerDist [options] -prefix PREF -input DSET\n"
110 "\n"
111 "where: \n"
112 "\n"
113 " -input DSET :(req) input dataset\n"
114 "\n"
115 " -prefix PREF :(req) output prefix name\n"
116 "\n"
117 " -mask MASK :mask dataset. NB: this mask is only applied *after*\n"
118 " the EDT has been calculated. Therefore, the boundaries\n"
119 " of this mask have no affect on the calculated distance\n"
120 " values, except for potentially zeroing some out at the\n"
121 " end.\n"
122 "\n"
123 " -dist_sq :by default, the output EDT volume contains distance\n"
124 " values. By using this option, the output values are\n"
125 " distance**2.\n"
126 "\n"
127 " -zeros_are_zero :by default, EDT values are output for the full FOV,\n"
128 " even zero-valued regions. If this option is used, EDT\n"
129 " values are only reported within the nonzero locations\n"
130 " of the input dataset.\n"
131 "\n"
132 " -zeros_are_neg :if this option is used, EDT in the zero/background\n"
133 " of the input will be negative (def: they are positive).\n"
134 " This opt cannot be used if '-zeros_are_zero' is.\n"
135 "\n"
136 " -nz_are_neg :if this option is used, EDT in the nonzero ROI regions\n"
137 " of the input will be negative (def: they are positive).\n"
138 "\n"
139 " -bounds_are_not_zero :this flag affects how FOV boundaries are treated for\n"
140 " nonzero ROIs: by default, they are viewed as ROI\n"
141 " boundaries (so the FOV is a closed boundary for an ROI,\n"
142 " as if the FOV were padded by an extra layer of zeros);\n"
143 " but when this option is used, the ROI behaves as if it\n"
144 " continued 'infinitely' at the FOV boundary (so it is\n"
145 " an open boundary). Zero-valued ROIs (= background)\n"
146 " are not affected by this option.\n"
147 "\n"
148 " -ignore_voxdims :this EDT algorithm works in terms of physical distance\n"
149 " and uses the voxel dimension info in each direction, by\n"
150 " default. However, using this option will ignore voxel\n"
151 " size, producing outputs as if each voxel dimension was\n"
152 " unity.\n"
153 "\n"
154 " -only2D SLI :instead of running full 3D EDT, run just in 2D, per.\n"
155 " plane. Provide the slice plane you want to run along\n"
156 " as the single argument SLI:\n"
157 " \"axi\" -> for axial slice\n"
158 " \"cor\" -> for coronal slice\n"
159 " \"sag\" -> for sagittal slice\n"
160 "\n"
161 " -binary_only :if the input is a binary mask or should be treated as\n"
162 " one (all nonzero voxels -> 1; all zeros stay 0), then\n"
163 " using this option will speed up the calculation. See\n"
164 " Notes below for more explanation of this. NOT ON YET!\n"
165 "\n"
166 " -verb V :manage verbosity when running code (def: 1).\n"
167 " Providing a V of 0 means to run quietly.\n"
168 "\n"
169 "==========================================================================\n"
170 "\n"
171 "Notes ~1~\n"
172 "\n"
173 "The original EDT algorithm of FH2012 was developed for a simple binary\n"
174 "mask input (and actually for homogeneous data grids of spacing=1). This\n"
175 "program, however, was built to handle more generalized cases of inputs,\n"
176 "namely ROI maps (and arbitrary voxel dimensions).\n"
177 "\n"
178 "The tradeoff of the expansion to handling ROI maps is an increase in\n"
179 "processing time---the original binary-mask algorithm is *very* efficient,\n"
180 "and the generalized one is still pretty quick but less so.\n"
181 "\n"
182 "So, if you know that your input should be treated as a binary mask, then\n"
183 "you can use the '-binary_only' option to utilize the more efficient\n"
184 "(and less generalized) algorithm. The output dataset should be the same\n"
185 "in either case---this option flag is purely about speed of computation.\n"
186 "\n"
187 "All other options about outputting dist**2 or negative values/etc. can be\n"
188 "used in conjunction with the '-binary_only', too.\n"
189 "\n"
190 "==========================================================================\n"
191 "\n"
192 "Examples ~1~\n"
193 "\n"
194 "1) Basic case:\n"
195 " 3dEulerDist \\\n"
196 " -input roi_map.nii.gz \\\n"
197 " -prefix roi_map_EDT.nii.gz \n"
198 "\n"
199 "2) Same as above, but only output distances within nonzero regions/ROIs:\n"
200 " 3dEulerDist \\\n"
201 " -zeros_are_zero \\\n"
202 " -input roi_map.nii.gz \\\n"
203 " -prefix roi_map_EDT_ZZ.nii.gz \n"
204 "\n"
205 "3) Output distance-squared at each voxel:\n"
206 " 3dEulerDist \\\n"
207 " -dist_sq \\\n"
208 " -input mask.nii.gz \\\n"
209 " -prefix mask_EDT_SQ.nii.gz \n"
210 "\n"
211 "4) Distinguish ROIs from nonzero background by making the former have\n"
212 " negative distance values in output:\n"
213 " 3dEulerDist \\\n"
214 " -nz_are_neg \\\n"
215 " -input roi_map.nii.gz \\\n"
216 " -prefix roi_map_EDT_NZNEG.nii.gz \n"
217 "\n"
218 "5) Have output voxel values represent (number of vox)**2 from a boundary;\n"
219 " voxel dimensions are ignored here:\n"
220 " 3dEulerDist \\\n"
221 " -ignore_voxdims \\\n"
222 " -dist_sq \\\n"
223 " -input roi_map.nii.gz \\\n"
224 " -prefix roi_map_EDT_SQ_VOX.nii.gz \n"
225 "\n"
226 "6) Basic case, with option for speed-up because the input is a binary mask\n"
227 " (i.e., only ones and zeros); any of the other above options can\n"
228 " be combined with this, too:\n"
229 " 3dEulerDist \\\n"
230 " -binary_only \\\n"
231 " -input roi_mask.nii.gz \\\n"
232 " -prefix roi_mask_EDT.nii.gz \n"
233 "\n"
234 "==========================================================================\n"
235 "\n",
236 author );
237
238 return 0;
239 }
240
main(int argc,char * argv[])241 int main(int argc, char *argv[]) {
242
243 int ii = 0;
244 int iarg;
245 PARAMS_euler_dist InOpts;
246 int itmp;
247
248 mainENTRY("3dEulerDist"); machdep();
249
250 // fill option struct with defaults
251 InOpts = set_euler_dist_defaults();
252
253 // ****************************************************************
254 // parse command line arguments
255 // ****************************************************************
256
257 if (argc == 1) { usage_3dEulerDist(); exit(0); }
258
259 /* scan through args */
260 iarg = 1;
261 while( iarg < argc && argv[iarg][0] == '-' ){
262
263 if( strcmp(argv[iarg],"-help") == 0 ||
264 strcmp(argv[iarg],"-h") == 0 ) {
265 usage_3dEulerDist();
266 exit(0);
267 }
268
269 if( strcmp(argv[iarg],"-input") == 0 ){
270 if( ++iarg >= argc )
271 ERROR_exit("Need argument after '%s'", argv[iarg-1]);
272
273 InOpts.input_name = strdup(argv[iarg]) ;
274 iarg++ ; continue ;
275 }
276
277 if( strcmp(argv[iarg],"-mask") == 0 ){
278 if( ++iarg >= argc )
279 ERROR_exit("Need argument after '%s'", argv[iarg-1]);
280
281 InOpts.mask_name = strdup(argv[iarg]) ;
282 iarg++ ; continue ;
283 }
284
285 if( strcmp(argv[iarg],"-prefix") == 0 ){
286 if( ++iarg >= argc )
287 ERROR_exit("Need argument after '%s'", argv[iarg-1]);
288
289 InOpts.prefix = strdup(argv[iarg]) ;
290
291 if( !THD_filename_ok(InOpts.prefix) )
292 ERROR_exit("Illegal name after '%s'", argv[iarg-1]);
293 iarg++ ; continue ;
294 }
295
296 if( strcmp(argv[iarg],"-zeros_are_zero") == 0) {
297 InOpts.zeros_are_zeroed = 1;
298 iarg++ ; continue ;
299 }
300
301 if( strcmp(argv[iarg],"-zeros_are_neg") == 0) {
302 InOpts.zero_region_sign = -1;
303 iarg++ ; continue ;
304 }
305
306 if( strcmp(argv[iarg],"-nz_are_neg") == 0) {
307 InOpts.nz_region_sign = -1;
308 iarg++ ; continue ;
309 }
310
311 if( strcmp(argv[iarg],"-bounds_are_not_zero") == 0) {
312 InOpts.bounds_are_zero = 0;
313 iarg++ ; continue ;
314 }
315
316 if( strcmp(argv[iarg],"-dist_sq") == 0) {
317 InOpts.dist_sq = 1;
318 iarg++ ; continue ;
319 }
320
321 if( strcmp(argv[iarg],"-ignore_voxdims") == 0) {
322 InOpts.ignore_voxdims = 1;
323 iarg++ ; continue ;
324 }
325
326 if( strcmp(argv[iarg],"-only2D") == 0) {
327 if( ++iarg >= argc )
328 ERROR_exit("Need argument after '%s'", argv[iarg-1]);
329
330 if( strcmp(argv[iarg],"cor") == 0 || \
331 strcmp(argv[iarg],"axi") == 0 || \
332 strcmp(argv[iarg],"sag") == 0 )
333 InOpts.only2D = strdup(argv[iarg]);
334 else
335 ERROR_exit("Need either \"cor\", \"axi\" or \"sag\" "
336 "after '%s'", argv[iarg-1]);
337
338 iarg++ ; continue ;
339 }
340
341 if( strcmp(argv[iarg],"-binary_only") == 0) {
342 InOpts.binary_only = 1;
343 iarg++ ; continue ;
344 }
345
346 if( strcmp(argv[iarg],"-verb") == 0) {
347 if( ++iarg >= argc )
348 ERROR_exit("Need int>0 argument after '%s'", argv[iarg-1]);
349
350 itmp = atoi(argv[iarg]);
351 if( itmp <= 0 )
352 InOpts.verb = 0;
353 else
354 InOpts.verb = itmp;
355
356 iarg++ ; continue ;
357 }
358
359 ERROR_message("Bad option '%s'\n",argv[iarg]);
360 suggest_best_prog_option(argv[0], argv[iarg]);
361 exit(1);
362 }
363
364 // ****************************************************************
365 // verify presence+behavior of inputs
366 // ****************************************************************
367
368 if ( InOpts.verb )
369 INFO_message("EDT: starting to check inputs...");
370
371 if ( !InOpts.input_name ) {
372 ERROR_message("You need to provide an input dset with '-input ..'");
373 exit(1);
374 }
375
376 if ( !InOpts.prefix )
377 ERROR_exit("Need an output name via '-prefix ..'\n");
378
379 if ( InOpts.zeros_are_zeroed && InOpts.zero_region_sign==-1 )
380 ERROR_exit("Cannot combine '-zeros_are_zero' and '-zeros_are_neg'. "
381 "You must choose which you *really* want.\n");
382
383 // DONE FILLING, now call
384 ii = run_EDT_3D(1, InOpts, argc, argv);
385
386 return 0;
387 }
388
389
run_EDT_3D(int comline,PARAMS_euler_dist opts,int argc,char * argv[])390 int run_EDT_3D( int comline, PARAMS_euler_dist opts,
391 int argc, char *argv[] )
392 {
393 int i, j, k, idx;
394 int nn;
395 int nx, ny, nz, nxy, nvox, nvals;
396 THD_3dim_dataset *dset_roi = NULL; // input
397 THD_3dim_dataset *dset_mask = NULL; // mask
398 THD_3dim_dataset *dset_edt = NULL; // output
399
400 ENTRY("run_EDT_3D");
401
402 dset_roi = THD_open_dataset(opts.input_name);
403 if( (dset_roi == NULL ))
404 ERROR_exit("Can't open dataset '%s'", opts.input_name);
405 DSET_load(dset_roi); CHECK_LOAD_ERROR(dset_roi);
406
407
408
409 if( opts.mask_name ) {
410 dset_mask = THD_open_dataset(opts.mask_name);
411 if( dset_mask == NULL )
412 ERROR_exit("Can't open dataset '%s'", opts.mask_name);
413 DSET_load(dset_mask); CHECK_LOAD_ERROR(dset_mask);
414
415 if( THD_dataset_mismatch( dset_roi , dset_mask ) )
416 ERROR_exit("Mismatch between input and mask dsets!\n");
417 }
418
419 // if only running in 2D, figure out which slice that is
420 if ( opts.only2D )
421 i = choose_axes_for_plane( dset_roi, opts.only2D,
422 opts.axes_to_proc, opts.verb );
423
424 nx = DSET_NX(dset_roi);
425 ny = DSET_NY(dset_roi);
426 nz = DSET_NZ(dset_roi);
427 nxy = nx*ny;
428 nvox = DSET_NVOX(dset_roi);
429 nvals = DSET_NVALS(dset_roi);
430
431 /* Prepare header for output by copying that of input, and then
432 changing items as necessary */
433 dset_edt = EDIT_empty_copy( dset_roi );
434 EDIT_dset_items(dset_edt,
435 ADN_datum_all, MRI_float ,
436 ADN_prefix, opts.prefix,
437 ADN_none );
438
439 if( opts.binary_only ){
440 for( nn=0 ; nn<nvals ; nn++ )
441 i = calc_EDT_3D_BIN(dset_edt, opts, dset_roi, dset_mask, nn);
442 }
443 else{
444 for( nn=0 ; nn<nvals ; nn++ )
445 i = calc_EDT_3D(dset_edt, opts, dset_roi, dset_mask, nn);
446 }
447
448 // free input dset
449 DSET_delete(dset_roi);
450 free(dset_roi);
451
452 // prepare to write output
453 THD_load_statistics( dset_edt );
454 if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(dset_edt)) )
455 ERROR_exit("Can't overwrite existing dataset '%s'",
456 DSET_HEADNAME(dset_edt));
457 tross_Make_History("3dEulerDist", argc, argv, dset_edt);
458
459 // write and free dset
460 THD_write_3dim_dataset(NULL, NULL, dset_edt, True);
461 DSET_delete(dset_edt);
462 free(dset_edt);
463
464 // free more
465 if( dset_mask ){
466 DSET_delete(dset_mask);
467 free(dset_mask);
468 }
469
470 return 0;
471 }
472
473