1 #include "mrilib.h"
2
3 #undef ALLOW_FILLIN /* 28 May 2002 */
4
5
main(int argc,char * argv[])6 int main( int argc , char * argv[] )
7 {
8 THD_3dim_dataset *dset , *mset=NULL, *masked_dset ;
9 char *amprefix = "automask" ;
10 char *prefix = NULL;
11 byte *mask ;
12 int iarg=1 , fillin=0 , nmask,nfill , dilate=0 , dd , erode = 0;
13 int dilate_flag = 0, erode_flag = 0;
14 float SIhh=0.0 ; /* 06 Mar 2003 */
15 int SIax=0 , SIbot=0,SItop=0 ;
16 int verb=1 ;
17 float clfrac=0.5 ; /* 20 Mar 2006 */
18 int peels=1, nbhrs=17 ; /* 24 Oct 2006 */
19 byte NN = 2; /* edge neighbors 23 Dec 2019 DRG */
20 byte new_NN = 0; /* flag user settings for neighbors and NNx*/
21 int new_nbhrs = 0;
22 int max_nbhrs = 18; /* maximum neighbors for NN2 */
23 int apply_mask = 0; /* 17 Nov 2009 */
24 char *apply_prefix = NULL, *depthprefix=NULL;
25 short *depth=NULL, dodepth=0; /* 02 March 2010 ZSS */
26
27
28 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
29 printf(
30 "Usage: 3dAutomask [options] dataset\n"
31 "Input dataset is EPI 3D+time, or a skull-stripped anatomical.\n"
32 "Output dataset is a brain-only mask dataset.\n"
33 "\n"
34 "This program by itself does NOT do 'skull-stripping'. Use\n"
35 "program 3dSkullStrip for that purpose!\n"
36 "\n"
37 "Method:\n"
38 " + Uses 3dClipLevel algorithm to find clipping level.\n"
39 " + Keeps only the largest connected component of the\n"
40 " supra-threshold voxels, after an erosion/dilation step.\n"
41 " + Writes result as a 'fim' type of functional dataset,\n"
42 " which will be 1 inside the mask and 0 outside the mask.\n"
43 "\n"
44 "Options:\n"
45 "--------\n"
46 " -prefix ppp = Write mask into dataset with prefix 'ppp'.\n"
47 " [Default == 'automask']\n"
48 "\n"
49 " -apply_prefix ppp = Apply mask to input dataset and save\n"
50 " masked dataset. If an apply_prefix is given\n"
51 " and not the usual prefix, the only output\n"
52 " will be the applied dataset\n"
53 "\n"
54 " -clfrac cc = Set the 'clip level fraction' to 'cc', which\n"
55 " must be a number between 0.1 and 0.9.\n"
56 " A small 'cc' means to make the initial threshold\n"
57 " for clipping (a la 3dClipLevel) smaller, which\n"
58 " will tend to make the mask larger. [default=0.5]\n"
59 "\n"
60 " -nograd = The program uses a 'gradual' clip level by default.\n"
61 " To use a fixed clip level, use '-nograd'.\n"
62 " [Change to gradual clip level made 24 Oct 2006.]\n"
63 "\n"
64 " -peels pp = Peel (erode) the mask 'pp' times, \n"
65 " then unpeel (dilate). Using NN2 neighborhoods,\n"
66 " clips off protuberances less than 2*pp voxels\n"
67 " thick. Turn off by setting to 0. [Default == 1]\n"
68 "\n"
69 " -NN1 -NN2 -NN3 = Erode and dilate using different neighbor definitions\n"
70 " NN1=faces, NN2=edges, NN3= corners [Default=NN2]\n"
71 " Applies to erode and dilate options, if present.\n"
72 " Note the default peeling processes still use NN2\n"
73 " unless the peels are set to 0\n"
74 "\n"
75 " -nbhrs nn = Define the number of neighbors needed for a voxel\n"
76 " NOT to be eroded. The 18 nearest neighbors in\n"
77 " the 3D lattice are used, so 'nn' should be between\n"
78 " 6 and 26. [Default == 17]\n"
79 "\n"
80 " -q = Don't write progress messages (i.e., be quiet).\n"
81 "\n"
82 " -eclip = After creating the mask, remove exterior\n"
83 " voxels below the clip threshold.\n"
84 "\n"
85 " -dilate nd = Dilate the mask outwards 'nd' times.\n"
86 "\n"
87 " -erode ne = Erode the mask inwards 'ne' times.\n"
88 #ifdef ALLOW_FILLIN
89 "\n"
90 " -fillin nnn = Fill in holes inside the mask of width up\n"
91 " to 'nnn' voxels. [Default == 0 == no fillin]\n"
92 #endif
93 "\n"
94 " -SI hh = After creating the mask, find the most superior\n"
95 " voxel, then zero out everything more than 'hh'\n"
96 " millimeters inferior to that. hh=130 seems to\n"
97 " be decent (i.e., for Homo Sapiens brains).\n"
98 "\n"
99 " -depth DEP = Produce a dataset (DEP) that shows how many peel \n"
100 " operations it takes to get to a voxel in the mask.\n"
101 " The higher the number, the deeper a voxel is located \n"
102 " in the mask. Note this uses the NN1,2,3 neighborhoods\n"
103 " above with a default of 2 for edge-sharing neighbors\n"
104 #ifdef ALLOW_FILLIN
105 " None of -peels, -dilate, -fillin, or -erode affect this option.\n"
106 #else
107 " None of -peels, -dilate, or -erode affect this option.\n"
108 #endif
109 ) ;
110
111 printf(
112 "--------------------------------------------------------------------\n"
113 "How to make an edge-of-brain mask from an anatomical volume:\n"
114 "* 3dSkullStrip to create a brain-only dataset; say, Astrip+orig\n"
115 "* 3dAutomask -prefix Amask Astrip+orig\n"
116 "* Create a mask of edge-only voxels via\n"
117 " 3dcalc -a Amask+orig -b a+i -c a-i -d a+j -e a-j -f a+k -g a-k \\\n"
118 " -expr 'ispositive(a)*amongst(0,b,c,d,e,f,g)' -prefix Aedge\n"
119 " which will be 1 at all voxels in the brain mask that have a\n"
120 " nearest neighbor that is NOT in the brain mask.\n"
121 "* cf. '3dcalc -help' DIFFERENTIAL SUBSCRIPTS for information\n"
122 " on the 'a+i' et cetera inputs used above.\n"
123 "* In regions where the brain mask is 'stair-stepping', then the\n"
124 " voxels buried inside the corner of the steps probably won't\n"
125 " show up in this edge mask:\n"
126 " ...00000000...\n"
127 " ...aaa00000...\n"
128 " ...bbbaa000...\n"
129 " ...bbbbbaa0...\n"
130 " Only the 'a' voxels are in this edge mask, and the 'b' voxels\n"
131 " down in the corners won't show up, because they only touch a\n"
132 " 0 voxel on a corner, not face-on. Depending on your use for\n"
133 " the edge mask, this effect may or may not be a problem.\n"
134 "--------------------------------------------------------------------\n"
135 ) ;
136
137 PRINT_COMPILE_DATE ; exit(0) ;
138 }
139
140 mainENTRY("3dAutomask main"); machdep(); AFNI_logger("3dAutomask",argc,argv);
141 PRINT_VERSION("3dAutomask") ; AUTHOR("Emperor Zhark") ;
142
143 /*-- options --*/
144
145 while( iarg < argc && argv[iarg][0] == '-' ){
146
147 if( strncmp(argv[iarg],"-peel",5) == 0 ){ /* 24 Oct 2006 */
148 peels = (int)strtod( argv[++iarg] , NULL ) ;
149 iarg++ ; continue ;
150 }
151
152 if( strcmp(argv[iarg],"-NN1") == 0 ){ /* 23 Dec 2019 */
153 new_NN = 1 ;
154 iarg++ ; continue ;
155 }
156
157 if( strcmp(argv[iarg],"-NN2") == 0 ){
158 new_NN = 2 ;
159 iarg++ ; continue ;
160 }
161
162 if( strcmp(argv[iarg],"-NN3") == 0 ){
163 new_NN = 3 ;
164 iarg++ ; continue ;
165 }
166
167 if( strncmp(argv[iarg],"-nbhr",5) == 0 ){ /* 24 Oct 2006 */
168 nbhrs = (int)strtod( argv[++iarg] , NULL ) ;
169 iarg++ ; continue ;
170 }
171
172 if( strncmp(argv[iarg],"-nograd",5) == 0 ){
173 THD_automask_set_gradualize(0) ;
174 iarg++ ; continue ;
175 }
176
177 if( strcmp(argv[iarg],"-clfrac") == 0 || strcmp(argv[iarg],"-mfrac") == 0 ){ /* 20 Mar 2006 */
178 clfrac = strtod( argv[++iarg] , NULL ) ;
179 if( clfrac < 0.1f || clfrac > 0.9f )
180 ERROR_exit("-clfrac value %f is illegal!",clfrac) ;
181 THD_automask_set_clipfrac(clfrac) ;
182 iarg++ ; continue ;
183 }
184
185 if( strcmp(argv[iarg],"-SI") == 0 ){ /* 06 Mar 2003 */
186 SIhh = strtod( argv[++iarg] , NULL ) ;
187 if( SIhh <= 0.0 )
188 ERROR_exit("-SI value %f is illegal!\n",SIhh) ;
189 iarg++ ; continue ;
190 }
191
192 if( strcmp(argv[iarg],"-eclip") == 0 ){ /* 28 Oct 2003 */
193 THD_automask_extclip(1) ;
194 iarg++ ; continue ;
195 }
196
197 if( strcmp(argv[iarg],"-depth") == 0 ){
198 depthprefix = argv[++iarg] ;
199 if( !THD_filename_ok(depthprefix) )
200 ERROR_exit("-depth %s is illegal!\n",depthprefix) ;
201 dodepth = 1 ;
202 iarg++ ; continue ;
203 }
204
205 if( strcmp(argv[iarg],"-q") == 0 ){ /* 28 Oct 2003 */
206 verb = 0 ; iarg++ ; continue ;
207 }
208
209 if( strcmp(argv[iarg],"-prefix") == 0 ){
210 prefix = argv[++iarg] ;
211 if( !THD_filename_ok(prefix) )
212 ERROR_exit("-prefix %s is illegal!\n",prefix) ;
213 iarg++ ; continue ;
214 }
215
216 #ifdef ALLOW_FILLIN
217 if( strcmp(argv[iarg],"-fillin") == 0 ){
218 fillin = strtol( argv[++iarg] , NULL , 10 ) ;
219 if( fillin < 0 )
220 ERROR_exit("-fillin %s is illegal!\n",argv[iarg]) ;
221 else if( fillin > 0 )
222 fillin = (fillin+2) / 2 ;
223 iarg++ ; continue ;
224 }
225 #endif
226
227 if( strcmp(argv[iarg],"-dilate") == 0 ){
228 dilate = strtol( argv[++iarg] , NULL , 10 ) ;
229 dilate_flag = 1;
230 if( dilate < 0 )
231 ERROR_exit("-dilate %s is illegal!\n",argv[iarg]);
232 iarg++ ; continue ;
233 }
234
235 if( strcmp(argv[iarg],"-erode") == 0 ){
236 erode = strtol( argv[++iarg] , NULL , 10 ) ;
237 erode_flag = 1;
238 if( erode < 0 )
239 ERROR_exit("-erode %s is illegal!\n",argv[iarg]);
240 iarg++ ; continue ;
241 }
242
243 if( strcmp(argv[iarg],"-apply_prefix") == 0 ){
244 apply_prefix = argv[++iarg] ;
245 if( !THD_filename_ok(apply_prefix) )
246 ERROR_exit("-apply_prefix %s is illegal!\n",apply_prefix) ;
247 apply_mask = 1;
248 iarg++ ; continue ;
249 }
250
251 ERROR_exit("ILLEGAL option: %s\n",argv[iarg]) ;
252 }
253
254 // nearest neighbor mode changes override neighbor choices
255 if(new_NN)
256 NN = new_NN;
257
258 switch(NN){
259 case(1):
260 max_nbhrs = 6;
261 break;
262 case(3):
263 max_nbhrs = 26;
264 break;
265
266 case(2):
267 default:
268 max_nbhrs = 18;
269 }
270
271 if(new_NN) /* keep historic default of 17 for NN2
272 unless changed via NNx or nbhrs options */
273 nbhrs = max_nbhrs;
274
275 /* check for valid number of neighbors */
276 if(new_nbhrs){
277 if((new_nbhrs< max_nbhrs) && (new_nbhrs>0))
278 nbhrs = new_nbhrs;
279 else
280 ERROR_exit("-nbhrs %d is illegal (choose between 1 and %d)!\n",
281 new_nbhrs, max_nbhrs);
282 }
283
284 /* set default number of peels (erosions/dilations combo, neighbors) */
285 THD_automask_set_peelcounts(peels,nbhrs) ;
286
287 if((dilate_flag+erode_flag)>1)
288 WARNING_message("Combining dilate,erode options is like peels option");
289
290 /*-- read data --*/
291
292 dset = THD_open_dataset(argv[iarg]); CHECK_OPEN_ERROR(dset,argv[iarg]);
293 if( DSET_BRICK_TYPE(dset,0) != MRI_short &&
294 DSET_BRICK_TYPE(dset,0) != MRI_byte &&
295 DSET_BRICK_TYPE(dset,0) != MRI_float ){
296 ERROR_exit("Illegal dataset datum type\n") ;
297 }
298 if( verb ) INFO_message("Loading dataset %s\n",argv[iarg]) ;
299 DSET_load(dset) ; CHECK_LOAD_ERROR(dset) ;
300
301 /*** do all the real work now ***/
302
303 if( verb ) INFO_message("Forming automask\n") ;
304 if( verb ) THD_automask_verbose(1) ;
305 mask = THD_automask( dset ) ;
306 if( mask == NULL )
307 ERROR_exit("Mask creation fails for unknown reasons!\n");
308
309 if (dodepth) { /* ZSS March 02 2010 */
310 if (!(depth = THD_mask_depth( DSET_NX(dset),
311 DSET_NY(dset),
312 DSET_NZ(dset), mask, 1, NULL, NN))) {
313 ERROR_exit("Failed to get depth vector!\n");
314 }
315 }
316
317 /* 30 Aug 2002 (modified 05 Mar 2003 to do fillin, etc, after dilation) */
318 if(dilate||erode) {
319 int ii,nx,ny,nz , nmm ;
320
321 nx = DSET_NX(dset) ; ny = DSET_NY(dset) ; nz = DSET_NZ(dset) ;
322 nmm = 1 ;
323 ii = rint(0.032*nx) ; nmm = MAX(nmm,ii) ;
324 ii = rint(0.032*ny) ; nmm = MAX(nmm,ii) ;
325 ii = rint(0.032*nz) ; nmm = MAX(nmm,ii) ;
326
327 if( verb && dilate) INFO_message("Dilating automask\n") ;
328 for( dd=0 ; dd < dilate ; dd++ ){
329 THD_mask_dilate ( nx,ny,nz , mask, 3, NN ) ;
330 THD_mask_fillin_completely( nx,ny,nz , mask, nmm ) ;
331 }
332
333 /* 3 May 2006 - drg- eroding option added */
334 if( verb && erode) INFO_message("Eroding automask\n") ;
335 for( dd=0 ; dd < erode ; dd++ ){
336 THD_mask_erode ( nx,ny,nz , mask, 0, NN) ;
337 THD_mask_fillin_completely( nx,ny,nz , mask, nmm ) ;
338 }
339
340 if(dilate||erode) {
341 nmm = nx*ny*nz ;
342 for( ii=0 ; ii < nmm ; ii++ ) mask[ii] = !mask[ii] ;
343 THD_mask_clust( nx,ny,nz, mask ) ;
344 for( ii=0 ; ii < nmm ; ii++ ) mask[ii] = !mask[ii] ;
345 }
346 }
347
348 /* 18 Apr 2002: print voxel count */
349
350 nmask = THD_countmask( DSET_NVOX(dset) , mask ) ;
351 if( verb ) INFO_message("%d voxels in the mask [out of %d: %.2f%%]\n",
352 nmask,DSET_NVOX(dset), (100.0*nmask)/DSET_NVOX(dset) ) ;
353 if( nmask == 0 )
354 ERROR_exit("No voxels? Quitting without saving mask\n");
355
356 /* 18 Apr 2002: maybe fill in voxels */
357
358 #ifdef ALLOW_FILLIN
359 if( fillin > 0 ){
360 nfill = THD_mask_fillin_completely(
361 DSET_NX(dset),DSET_NY(dset),DSET_NZ(dset), mask, fillin ) ;
362 if( verb ) INFO_message("%d voxels filled in; %d voxels total\n",
363 nfill,nfill+nmask ) ;
364 }
365 #endif
366
367 /** 04 Jun 2002: print cut plane report **/
368
369 { int nx=DSET_NX(dset), ny=DSET_NY(dset), nz=DSET_NZ(dset), nxy=nx*ny ;
370 int ii,jj,kk ;
371
372 #if 0
373 { int xm=-1,xp=-1,ym=-1,yp=-1,zm=-1,zp=-1 ;
374 THD_autobbox( dset , &xm,&xp , &ym,&yp , &zm,&zp, NULL ) ;
375 INFO_message("Auto bbox: x=%d..%d y=%d..%d z=%d..%d\n",
376 xm,xp,ym,yp,zm,zp ) ;
377 }
378 #endif
379
380 for( ii=0 ; ii < nx ; ii++ )
381 for( kk=0 ; kk < nz ; kk++ )
382 for( jj=0 ; jj < ny ; jj++ )
383 if( mask[ii+jj*nx+kk*nxy] ) goto CP5 ;
384 CP5: if( verb )
385 INFO_message("first %3d x-planes are zero [from %c]\n",
386 ii,ORIENT_tinystr[dset->daxes->xxorient][0]) ;
387 if( SIhh > 0.0 && ORIENT_tinystr[dset->daxes->xxorient][0] == 'S' ){
388 SIax = 1 ; SIbot = ii + (int)(SIhh/fabs(DSET_DX(dset))+0.5) ; SItop = nx-1 ;
389 }
390
391 for( ii=nx-1 ; ii >= 0 ; ii-- )
392 for( kk=0 ; kk < nz ; kk++ )
393 for( jj=0 ; jj < ny ; jj++ )
394 if( mask[ii+jj*nx+kk*nxy] ) goto CP6 ;
395 CP6: if( verb )
396 INFO_message("last %3d x-planes are zero [from %c]\n",
397 nx-1-ii,ORIENT_tinystr[dset->daxes->xxorient][1]) ;
398 if( SIhh > 0.0 && ORIENT_tinystr[dset->daxes->xxorient][1] == 'S' ){
399 SIax = 1 ; SIbot = 0 ; SItop = ii - (int)(SIhh/fabs(DSET_DX(dset))+0.5) ;
400 }
401
402 for( jj=0 ; jj < ny ; jj++ )
403 for( kk=0 ; kk < nz ; kk++ )
404 for( ii=0 ; ii < nx ; ii++ )
405 if( mask[ii+jj*nx+kk*nxy] ) goto CP3 ;
406 CP3: if( verb )
407 INFO_message("first %3d y-planes are zero [from %c]\n",
408 jj,ORIENT_tinystr[dset->daxes->yyorient][0]) ;
409 if( SIhh > 0.0 && ORIENT_tinystr[dset->daxes->yyorient][0] == 'S' ){
410 SIax = 2 ; SIbot = jj + (int)(SIhh/fabs(DSET_DY(dset))+0.5) ; SItop = ny-1 ;
411 }
412
413 for( jj=ny-1 ; jj >= 0 ; jj-- )
414 for( kk=0 ; kk < nz ; kk++ )
415 for( ii=0 ; ii < nx ; ii++ )
416 if( mask[ii+jj*nx+kk*nxy] ) goto CP4 ;
417 CP4: if( verb )
418 INFO_message("last %3d y-planes are zero [from %c]\n",
419 ny-1-jj,ORIENT_tinystr[dset->daxes->yyorient][1]) ;
420 if( SIhh > 0.0 && ORIENT_tinystr[dset->daxes->yyorient][1] == 'S' ){
421 SIax = 2 ; SIbot = 0 ; SItop = jj - (int)(SIhh/fabs(DSET_DY(dset))+0.5) ;
422 }
423
424 for( kk=0 ; kk < nz ; kk++ )
425 for( jj=0 ; jj < ny ; jj++ )
426 for( ii=0 ; ii < nx ; ii++ )
427 if( mask[ii+jj*nx+kk*nxy] ) goto CP1 ;
428 CP1: if( verb )
429 INFO_message("first %3d z-planes are zero [from %c]\n",
430 kk,ORIENT_tinystr[dset->daxes->zzorient][0]) ;
431 if( SIhh > 0.0 && ORIENT_tinystr[dset->daxes->zzorient][0] == 'S' ){
432 SIax = 3 ; SIbot = kk + (int)(SIhh/fabs(DSET_DZ(dset))+0.5) ; SItop = nz-1 ;
433 }
434
435 for( kk=nz-1 ; kk >= 0 ; kk-- )
436 for( jj=0 ; jj < ny ; jj++ )
437 for( ii=0 ; ii < nx ; ii++ )
438 if( mask[ii+jj*nx+kk*nxy] ) goto CP2 ;
439 CP2: if( verb )
440 INFO_message("last %3d z-planes are zero [from %c]\n",
441 nz-1-kk,ORIENT_tinystr[dset->daxes->zzorient][1]) ;
442 if( SIhh > 0.0 && ORIENT_tinystr[dset->daxes->zzorient][1] == 'S' ){
443 SIax = 3 ; SIbot = 0 ; SItop = kk - (int)(SIhh/fabs(DSET_DZ(dset))+0.5) ;
444 }
445
446 /* 06 Mar 2003: cut off stuff below SIhh mm from most Superior point */
447
448 if( SIax > 0 && SIbot <= SItop ){
449 char *cax="xyz" ;
450 if( verb )
451 INFO_message("SI clipping mask along axis %c from %d..%d\n" ,
452 cax[SIax-1] , SIbot,SItop ) ;
453 switch( SIax ){
454 case 1:
455 for( ii=SIbot ; ii <= SItop ; ii++ )
456 for( kk=0 ; kk < nz ; kk++ )
457 for( jj=0 ; jj < ny ; jj++ ) mask[ii+jj*nx+kk*nxy] = 0 ;
458 break ;
459 case 2:
460 for( jj=SIbot ; jj <= SItop ; jj++ )
461 for( kk=0 ; kk < nz ; kk++ )
462 for( ii=0 ; ii < nx ; ii++ ) mask[ii+jj*nx+kk*nxy] = 0 ;
463 break ;
464 case 3:
465 for( kk=SIbot ; kk <= SItop ; kk++ )
466 for( jj=0 ; jj < ny ; jj++ )
467 for( ii=0 ; ii < nx ; ii++ ) mask[ii+jj*nx+kk*nxy] = 0 ;
468 break ;
469 }
470 nmask = THD_countmask( DSET_NVOX(dset) , mask ) ;
471 if( verb )
472 INFO_message("%d voxels left [out of %d]\n",nmask,DSET_NVOX(dset)) ;
473 }
474 }
475
476
477 /* create output dataset */
478 if(prefix==NULL) {
479 if(apply_prefix == NULL)
480 prefix = amprefix;
481 }
482
483 if(prefix) {
484 mset = EDIT_empty_copy( dset ) ;
485 EDIT_dset_items( mset ,
486 ADN_prefix , prefix ,
487 ADN_datum_all , MRI_byte ,
488 ADN_nvals , 1 ,
489 ADN_ntt , 0 ,
490 ADN_type , HEAD_FUNC_TYPE ,
491 ADN_func_type , FUNC_FIM_TYPE ,
492 ADN_none ) ;
493 EDIT_substitute_brick( mset , 0 , MRI_byte , mask ) ;
494
495 /* 16 Apr 2002: make history */
496
497 tross_Copy_History( dset , mset ) ;
498 tross_Make_History( "3dAutomask", argc,argv, mset ) ;
499
500 DSET_write( mset ) ;
501 if( verb ) WROTE_DSET(mset) ;
502 }
503
504 if (apply_mask) {
505 if(verb) INFO_message("applying mask to original data\n");
506 masked_dset = thd_apply_mask(dset, mask, apply_prefix);
507 if(masked_dset){
508 if(verb) INFO_message("Writing masked data\n");
509 tross_Copy_History( dset , masked_dset ) ;
510 tross_Make_History( "3dAutomask", argc,argv, masked_dset ) ;
511 DSET_write(masked_dset);
512 if( verb ) WROTE_DSET(masked_dset) ;
513 DSET_unload( masked_dset ) ; /* don't need data any more */
514 }
515 else {
516 ERROR_exit("Could not apply mask to dataset");
517 }
518 }
519
520 if (dodepth) { /* ZSS March 02 2010 */
521 if(verb) INFO_message("Writing depth dataset\n");
522 if (mset) DSET_delete(mset); mset = NULL;
523 mset = EDIT_empty_copy( dset ) ;
524 EDIT_dset_items( mset ,
525 ADN_prefix , depthprefix ,
526 ADN_datum_all , MRI_short ,
527 ADN_nvals , 1 ,
528 ADN_ntt , 0 ,
529 ADN_type , HEAD_FUNC_TYPE ,
530 ADN_func_type , FUNC_FIM_TYPE ,
531 ADN_none ) ;
532 EDIT_substitute_brick( mset , 0 , MRI_short , depth ) ;
533
534 /* 16 Apr 2002: make history */
535
536 tross_Copy_History( dset , mset ) ;
537 tross_Make_History( "3dAutomask", argc,argv, mset ) ;
538
539 DSET_write( mset ) ;
540 if( verb ) WROTE_DSET(mset) ;
541
542 DSET_delete(mset); mset = NULL;
543 }
544
545 DSET_unload( dset ) ; /* don't need data any more */
546
547 if( verb ) INFO_message("CPU time = %f sec\n",COX_cpu_time()) ;
548
549 exit(0) ;
550 }
551