1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1999-2001, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*---------------------------------------------------------------------------*/
8 /*
9 Program to automatically segment the intracranial region.
10
11 File: 3dIntracranial.c
12 Author: B. Douglas Ward
13 Date: 04 June 1999
14
15
16 Mod: Corrected initialization of PDF estimation flag.
17 Date: 06 August 1999
18
19 Mod: Added changes for incorporating History notes.
20 Date: 09 September 1999
21
22 Mod: Interface changes for new estpdf.c
23 Date: 28 January 2000
24
25 Mod: Added call to AFNI_logger.
26 Date: 15 August 2001
27
28 Mod: Added option to suppress spatial smoothing of segmentation mask.
29 Date: 03 December 2001
30
31 Mod: Set MAX_STRING_LENGTH equal to THD_MAX_NAME.
32 Date: 02 December 2002
33
34 Mod: Convert input from MRI_byte to MRI_short if needed.
35 Date: 05 December 2002
36
37 Mod: Modified the -help menu -- P Christidis
38 Date: 21 July 2005
39 */
40
41 /*---------------------------------------------------------------------------*/
42
43 #define PROGRAM_NAME "3dIntracranial" /* name of this program */
44 #define PROGRAM_AUTHOR "B. D. Ward" /* program author */
45 #define PROGRAM_INITIAL "04 June 1999" /* date of initial program release */
46 #define PROGRAM_LATEST "21 July 2005" /* date of latest program revision */
47
48 /*---------------------------------------------------------------------------*/
49 /*
50 Include header files.
51 */
52
53
54 #include "mrilib.h"
55 #include "Intracranial.h"
56
57
58 /*---------------------------------------------------------------------------*/
59 /*
60 Global variables and constants.
61 */
62
63 #define USE_QUIET
64
65 static char * anat_filename = NULL; /* file name for input anat dataset */
66 static char * prefix_filename = NULL; /* prefix name for output dataset */
67 static RwcBoolean write_mask = FALSE; /* flag for generate 'fim' mask */
68 static RwcBoolean quiet = FALSE; /* flag for suppress screen output */
69 static char * commandline = NULL ; /* command line for history notes */
70 static RwcBoolean nosmooth = FALSE; /* flag to delete spatial smoothing */
71
72
73 /*---------------------------------------------------------------------------*/
74 /*
75 Print error message and stop.
76 */
77
SI_error(char * message)78 void SI_error (char * message)
79 {
80 fprintf (stderr, "\n");
81 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
82 exit(1);
83 }
84
85
86 /*---------------------------------------------------------------------------*/
87 /*
88 Include source code.
89 */
90
91 #include "estpdf3.c" /* code for PDF estimation */
92 #include "Intracranial.c" /* code for intracranial segmentation */
93
94
95 /*---------------------------------------------------------------------------*/
96 /*
97 Routine to display 3dIntracranial help menu.
98 */
99
display_help_menu()100 void display_help_menu()
101 {
102 printf
103 (
104 "3dIntracranial - performs automatic segmentation of intracranial region.\n"
105 " \n"
106 " This program will strip the scalp and other non-brain tissue from a \n"
107 " high-resolution T1 weighted anatomical dataset. \n"
108 " \n"
109 "** Nota Bene: the newer program 3dSkullStrip should also be considered \n"
110 "** for this functionality -- it usually works better. \n"
111 " \n"
112 "----------------------------------------------------------------------- \n"
113 " \n"
114 "Usage: \n"
115 "----- \n"
116 " \n"
117 "3dIntracranial \n"
118 " -anat filename => Filename of anat dataset to be segmented \n"
119 " \n"
120 " [-min_val a] => Minimum voxel intensity limit \n"
121 " Default: Internal PDF estimate for lower bound \n"
122 " \n"
123 " [-max_val b] => Maximum voxel intensity limit \n"
124 " Default: Internal PDF estimate for upper bound \n"
125 " \n"
126 " [-min_conn m] => Minimum voxel connectivity to enter \n"
127 " Default: m=4 \n"
128 " \n"
129 " [-max_conn n] => Maximum voxel connectivity to leave \n"
130 " Default: n=2 \n"
131 " \n"
132 " [-nosmooth] => Suppress spatial smoothing of segmentation mask \n"
133 " \n"
134 " [-mask] => Generate functional image mask (complement) \n"
135 " Default: Generate anatomical image \n"
136 " \n"
137 " [-quiet] => Suppress output to screen \n"
138 " \n"
139 " -prefix pname => Prefix name for file to contain segmented image \n"
140 " \n"
141 " ** NOTE **: The newer program 3dSkullStrip will probably give \n"
142 " better segmentation results than 3dIntracranial! \n"
143
144 "----------------------------------------------------------------------- \n"
145 " \n"
146 "Examples: \n"
147 "-------- \n"
148 " \n"
149 " 3dIntracranial -anat elvis+orig -prefix elvis_strip \n"
150 " \n"
151 " 3dIntracranial -min_val 30 -max_val 350 -anat elvis+orig -prefix strip\n"
152 " \n"
153 " 3dIntracranial -nosmooth -quiet -anat elvis+orig -prefix elvis_strip \n"
154 " \n"
155 "----------------------------------------------------------------------- \n"
156 );
157
158 PRINT_COMPILE_DATE ; exit(0);
159 }
160
161
162 /*---------------------------------------------------------------------------*/
163 /*
164 Routine to get user specified input options.
165 */
166
get_options(int argc,char ** argv)167 void get_options
168 (
169 int argc, /* number of input arguments */
170 char ** argv /* array of input arguments */
171 )
172
173 {
174 int nopt = 1; /* input option argument counter */
175 int ival; /* integer input */
176 float fval; /* float input */
177 char message[MAX_STRING_LENGTH]; /* error message */
178
179
180 /*----- does user request help menu? -----*/
181 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
182
183
184 /*----- add to program log -----*/
185 AFNI_logger (PROGRAM_NAME,argc,argv);
186
187
188 /*----- main loop over input options -----*/
189 while (nopt < argc )
190 {
191
192 /*----- -anat filename -----*/
193 if (strncmp(argv[nopt], "-anat", 5) == 0)
194 {
195 nopt++;
196 if (nopt >= argc) SI_error ("need argument after -anat ");
197 anat_filename = malloc (sizeof(char) * MAX_STRING_LENGTH);
198 MTEST (anat_filename);
199 strcpy (anat_filename, argv[nopt]);
200
201 anat = THD_open_dataset (anat_filename);
202 CHECK_OPEN_ERROR(anat,anat_filename) ;
203 DSET_load(anat) ; CHECK_LOAD_ERROR(anat) ;
204
205 /** RWCox [05 Dec 2002]
206 If input is a byte dataset, make a short copy of it. **/
207
208 if( DSET_BRICK_TYPE(anat,0) == MRI_byte ){
209 THD_3dim_dataset *qset ;
210 register byte *bar ; register short *sar ;
211 register int ii,nvox ;
212 fprintf(stderr,"++ WARNING: converting input dataset from byte to short\n") ;
213 qset = EDIT_empty_copy(anat) ;
214 nvox = DSET_NVOX(anat) ;
215 bar = (byte *) DSET_ARRAY(anat,0) ;
216 sar = (short *)malloc(sizeof(short)*nvox) ;
217 for( ii=0 ; ii < nvox ; ii++ ) sar[ii] = (short) bar[ii] ;
218 EDIT_substitute_brick( qset , 0 , MRI_short , sar ) ;
219 DSET_delete(anat) ; anat = qset ;
220 }
221
222 nopt++;
223 continue;
224 }
225
226
227 /*----- -min_val a -----*/
228 if (strncmp(argv[nopt], "-min_val", 8) == 0)
229 {
230 nopt++;
231 if (nopt >= argc) SI_error ("need argument after -min_val ");
232 sscanf (argv[nopt], "%f", &fval);
233 if (fval < 0.0)
234 SI_error ("illegal argument after -min_val ");
235 min_val_float = fval;
236 min_val_int = 0;
237 nopt++;
238 continue;
239 }
240
241
242 /*----- -max_val b -----*/
243 if (strncmp(argv[nopt], "-max_val", 8) == 0)
244 {
245 nopt++;
246 if (nopt >= argc) SI_error ("need argument after -max_val ");
247 sscanf (argv[nopt], "%f", &fval);
248 if (fval < 0.0)
249 SI_error ("illegal argument after -max_val ");
250 max_val_float = fval;
251 max_val_int = 0;
252 nopt++;
253 continue;
254 }
255
256
257 /*----- -min_conn m -----*/
258 if (strncmp(argv[nopt], "-min_conn", 9) == 0)
259 {
260 nopt++;
261 if (nopt >= argc) SI_error ("need argument after -min_conn ");
262 sscanf (argv[nopt], "%d", &ival);
263 if ((ival < 1) || (ival > 7))
264 SI_error ("illegal argument after -min_conn ");
265 min_conn_int = ival;
266 nopt++;
267 continue;
268 }
269
270
271 /*----- -max_conn n -----*/
272 if (strncmp(argv[nopt], "-max_conn", 9) == 0)
273 {
274 nopt++;
275 if (nopt >= argc) SI_error ("need argument after -max_conn ");
276 sscanf (argv[nopt], "%d", &ival);
277 if ((ival < -1) || (ival > 5))
278 SI_error ("illegal argument after -max_conn ");
279 max_conn_int = ival;
280 nopt++;
281 continue;
282 }
283
284
285 /*----- -nosmooth -----*/
286 if (strcmp(argv[nopt], "-nosmooth") == 0)
287 {
288 nosmooth = TRUE;
289 nopt++;
290 continue;
291 }
292
293
294 /*----- -mask -----*/
295 if (strncmp(argv[nopt], "-mask", 5) == 0)
296 {
297 write_mask = TRUE;
298 nopt++;
299 continue;
300 }
301
302
303 /*----- -quiet -----*/
304 if (strncmp(argv[nopt], "-quiet", 6) == 0)
305 {
306 quiet = TRUE;
307 nopt++;
308 continue;
309 }
310
311
312 /*----- -prefix prefixname -----*/
313 if (strncmp(argv[nopt], "-prefix", 7) == 0)
314 {
315 nopt++;
316 if (nopt >= argc) SI_error ("need argument after -prefix ");
317 prefix_filename = malloc (sizeof(char) * MAX_STRING_LENGTH);
318 MTEST (prefix_filename);
319 strcpy (prefix_filename, argv[nopt]);
320 nopt++;
321 continue;
322 }
323
324
325 /*----- unknown command -----*/
326 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
327 SI_error (message);
328
329 }
330
331
332 }
333
334
335 /*---------------------------------------------------------------------------*/
336 /*
337 Routine to check whether one output file already exists.
338 */
339
check_one_output_file(char * filename)340 void check_one_output_file
341 (
342 char * filename /* name of output file */
343 )
344
345 {
346 char message[MAX_STRING_LENGTH]; /* error message */
347 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
348 int ierror; /* number of errors in editing data */
349
350
351 /*----- make an empty copy of input dataset -----*/
352 new_dset = EDIT_empty_copy ( anat );
353
354
355 ierror = EDIT_dset_items( new_dset ,
356 ADN_prefix , filename ,
357 ADN_label1 , filename ,
358 ADN_self_name , filename ,
359 ADN_none ) ;
360
361 if( ierror > 0 )
362 {
363 sprintf (message,
364 "*** %d errors in attempting to create output dataset!\n",
365 ierror);
366 SI_error (message);
367 }
368
369 if( THD_is_file(new_dset->dblk->diskptr->header_name) )
370 {
371 sprintf (message,
372 "Output dataset file %s already exists"
373 "--cannot continue!\a\n",
374 new_dset->dblk->diskptr->header_name);
375 SI_error (message);
376 }
377
378 /*----- deallocate memory -----*/
379 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
380
381 }
382
383
384 /*---------------------------------------------------------------------------*/
385 /*
386 Program initialization.
387 */
388
initialize_program(int argc,char ** argv)389 void initialize_program
390 (
391 int argc, /* number of input arguments */
392 char ** argv /* array of input arguments */
393 )
394
395 {
396 float parameters [DIMENSION]; /* parameters for PDF estimation */
397 int nxyz;
398 short * sfim;
399 int ixyz, icount;
400 short * rfim;
401 int lower_cutoff = 25;
402
403
404 /*----- save command line for history notes -----*/
405 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
406
407
408 /*----- Get operator inputs -----*/
409 get_options (argc, argv);
410
411
412 /*----- Verify that inputs are acceptable -----*/
413 verify_inputs();
414
415
416 /*----- Initialize local variables -----*/
417 if (anat == NULL) SI_error ("Unable to read anat dataset");
418 nxyz = DSET_NX(anat) * DSET_NY(anat) * DSET_NZ(anat);
419 sfim = (short *) DSET_BRICK_ARRAY(anat,0) ;
420 if (sfim == NULL) SI_error ("Unable to read anat dataset");
421 rfim = (short *) malloc (sizeof(short) * nxyz); MTEST (rfim);
422
423
424 /*----- Just use voxels whose intensity is above the lower cutoff -----*/
425 icount = 0;
426 for (ixyz = 0; ixyz < nxyz; ixyz++)
427 if (sfim[ixyz] > lower_cutoff)
428 {
429 rfim[icount] = sfim[ixyz];
430 icount++;
431 }
432 if (! quiet)
433 printf ("%d voxels above lower cutoff = %d \n", icount, lower_cutoff);
434
435
436 /*----- Get PDF estimate and set voxel intensity limits -----*/
437 if (min_val_int || max_val_int) estpdf_short (icount, rfim, parameters);
438 if (min_val_int) min_val_float = parameters[4] - 2.0*parameters[5];
439 if (max_val_int) max_val_float = parameters[7] + 2.0*parameters[8];
440
441
442 if (! quiet)
443 {
444 printf ("\n");
445 printf ("Control inputs: \n");
446 printf ("anat filename = %s \n", anat_filename);
447 printf ("min value = %f \n", min_val_float);
448 printf ("max value = %f \n", max_val_float);
449 printf ("min conn = %d \n", min_conn_int);
450 printf ("max conn = %d \n", max_conn_int);
451 printf ("prefix filename = %s \n", prefix_filename);
452 }
453
454
455 }
456
457
458 /*---------------------------------------------------------------------------*/
459 /*
460 Perform initialization for automatic segmentation algorithm.
461 */
462
auto_initialize(short ** cv)463 int auto_initialize (short ** cv)
464
465 {
466 int nx, ny, nz, nxy, nxyz, ixyz; /* voxel counters */
467
468
469 /*----- Initialize local variables -----*/
470 nx = DSET_NX(anat); ny = DSET_NY(anat); nz = DSET_NZ(anat);
471 nxy = nx*ny; nxyz = nxy*nz;
472
473
474 /*----- Initialize voxel classification indicators -----*/
475 *cv = (short *) malloc (nxyz * sizeof(short));
476 MTEST (*cv);
477 for (ixyz = 0; ixyz < nxyz; ixyz++)
478 (*cv)[ixyz] = 0;
479
480
481 /*----- Initialize random number generator -----*/
482 srand48 (1234567);
483
484
485 return (1);
486 }
487
488
489 /*---------------------------------------------------------------------------*/
490 /*
491 Put estimated target structure into output dataset.
492 */
493
target_into_dataset(short * cv)494 void target_into_dataset
495 (
496 short * cv /* volume with 1's at target voxel locations */
497 )
498
499 {
500 short * anat_data = NULL; /* data from anatomical image */
501 int nx, ny, nz, nxy, nxyz, ixyz; /* voxel counters */
502
503
504 /*----- Initialize local variables -----*/
505 anat_data = (short *) DSET_BRICK_ARRAY(anat,0) ;
506 nx = DSET_NX(anat); ny = DSET_NY(anat); nz = DSET_NZ(anat);
507 nxy = nx*ny; nxyz = nxy*nz;
508
509 if (! write_mask)
510 {
511 /*----- Reset to zero those voxels which lie outside the brain -----*/
512 for (ixyz = 0; ixyz < nxyz; ixyz++)
513 {
514 if (cv[ixyz]) cv[ixyz] = 0;
515 else cv[ixyz] = anat_data[ixyz];
516 }
517 }
518
519
520 /*----- deallocate memory -----*/
521 THD_delete_3dim_dataset (anat, False); anat = NULL;
522
523
524 return;
525 }
526
527
528 /*---------------------------------------------------------------------------*/
529 /*
530 Routine to write one AFNI dataset.
531
532 The output is either a segmented anatomical dataset,
533 or a mask 'fim' type dataset.
534
535 */
536
write_afni_data(short * cv)537 void write_afni_data
538 (
539 short * cv /* volume with 1's at target voxel locations */
540 )
541
542 {
543 int nxyz; /* number of voxels */
544 int ii; /* voxel index */
545 THD_3dim_dataset * dset=NULL; /* input afni data set pointer */
546 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
547 int ierror; /* number of errors in editing data */
548 int ibuf[32]; /* integer buffer */
549 float fbuf[MAX_STAT_AUX]; /* float buffer */
550 float fimfac; /* scale factor for short data */
551 int output_datum; /* data type for output data */
552 char * filename; /* prefix filename for output */
553
554
555 /*----- initialize local variables -----*/
556 filename = prefix_filename;
557 dset = THD_open_dataset (anat_filename); CHECK_OPEN_ERROR(dset,anat_filename);
558 nxyz = DSET_NX(dset) * DSET_NY(dset) * DSET_NZ(dset);
559
560
561 /*-- make an empty copy of this dataset, for eventual output --*/
562 new_dset = EDIT_empty_copy( dset ) ;
563
564
565 /*----- Record history of dataset -----*/
566 tross_Copy_History( dset , new_dset ) ;
567 if( commandline != NULL )
568 tross_Append_History( new_dset , commandline ) ;
569
570
571 /*----- deallocate memory -----*/
572 THD_delete_3dim_dataset (dset, False); dset = NULL ;
573
574
575 output_datum = MRI_short ;
576
577 ibuf[0] = output_datum ;
578
579 if (write_mask)
580 {
581 int func_type = FUNC_FIM_TYPE;
582 ierror = EDIT_dset_items( new_dset ,
583 ADN_prefix , filename ,
584 ADN_label1 , filename ,
585 ADN_self_name , filename ,
586 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
587 ADN_func_type , func_type ,
588 ADN_nvals , FUNC_nvals[func_type] ,
589 ADN_datum_array , ibuf ,
590 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
591 ADN_none ) ;
592 }
593 else
594 {
595 ierror = EDIT_dset_items( new_dset ,
596 ADN_prefix , filename ,
597 ADN_label1 , filename ,
598 ADN_self_name , filename ,
599 ADN_datum_array , ibuf ,
600 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
601 ADN_none ) ;
602 }
603
604
605 if( ierror > 0 ){
606 fprintf(stderr,
607 "*** %d errors in attempting to create output dataset!\n",
608 ierror ) ;
609 exit(1) ;
610 }
611
612
613 if( THD_deathcon() && THD_is_file(new_dset->dblk->diskptr->header_name) ){
614 fprintf(stderr,
615 "*** Output dataset file %s already exists--cannot continue!\a\n",
616 new_dset->dblk->diskptr->header_name ) ;
617 exit(1) ;
618 }
619
620
621 /*----- attach bricks to new data set -----*/
622 mri_fix_data_pointer (cv, DSET_BRICK(new_dset,0));
623 fimfac = 1.0;
624
625
626 /*----- write afni data set -----*/
627 if (! quiet)
628 {
629 if (write_mask) printf("\nWriting functional (mask) dataset: ");
630 else printf ("\nWriting anatomical dataset: ");
631 printf("%s\n", DSET_BRIKNAME(new_dset) ) ;
632 }
633
634 printf("NOTE: You will probably get better results by using\n"
635 " 3dSkullStrip -input %s -prefix %s\n" ,
636 anat_filename , prefix_filename ) ;
637
638 for( ii=0 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ;
639 (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ;
640
641 fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ;
642 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
643
644 THD_load_statistics( new_dset ) ;
645 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
646
647
648 /*----- deallocate memory -----*/
649 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
650
651 }
652
653
654 /*---------------------------------------------------------------------------*/
655 /*
656 Perform automatic image segmentation.
657 */
658
SEGMENT_auto()659 void SEGMENT_auto ()
660 {
661 short * cv = NULL; /* volume with 1's at target voxel locations */
662 int ok; /* binary, for successful operation */
663
664
665 /*----- Perform initialization for automatic segmentation algorithm -----*/
666 ok = auto_initialize (&cv);
667 if (! ok) return;
668
669
670 /*----- Segment intracranial voxels -----*/
671 segment_volume (cv);
672
673 /*----- Perform voxel connectivity tests -----*/
674 connectivity_tests (cv);
675
676 /*----- Put estimated target structure into output dataset -----*/
677 target_into_dataset (cv);
678
679 /*----- Write out the segmented dataset -----*/
680 write_afni_data (cv);
681
682 return ;
683
684 }
685
686
687 /*---------------------------------------------------------------------------*/
688 /*
689 This is the main routine for program 3dIntracranial.
690 */
691
main(int argc,char ** argv)692 int main
693 (
694 int argc, /* number of input arguments */
695 char ** argv /* array of input arguments */
696 )
697
698 {
699 int ii ;
700
701 /*----- Identify software -----*/
702
703 #if 0
704 for( ii=1 ; ii < argc ; ii++ )
705 if( strncmp(argv[ii],"-quiet",6) == 0 ) break ;
706 if( ii == argc ){
707 printf ("\n\n");
708 printf ("Program: %s \n", PROGRAM_NAME);
709 printf ("Author: %s \n", PROGRAM_AUTHOR);
710 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
711 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
712 printf ("\n");
713 }
714 #endif
715
716 mainENTRY("3dIntracranial:main") ; machdep() ; PRINT_VERSION("3dIntracranial") ;
717 AUTHOR(PROGRAM_AUTHOR) ;
718 WARNING_message("This program (3dIntracranial) is old, obsolete, and not maintained!") ;
719 INFO_message("3dSkullStrip is almost always superior to 3dIntracranial :)") ;
720
721
722 /*----- Program initialization -----*/
723 initialize_program (argc, argv);
724
725
726 /*----- Perform automatic segmentation -----*/
727 SEGMENT_auto();
728
729 exit(0);
730 }
731
732 /*---------------------------------------------------------------------------*/
733
734