00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #define PROGRAM_NAME "3dIntracranial"
00044 #define PROGRAM_AUTHOR "B. D. Ward"
00045 #define PROGRAM_INITIAL "04 June 1999"
00046 #define PROGRAM_LATEST "21 July 2005"
00047
00048
00049
00050
00051
00052
00053
00054 #include "mrilib.h"
00055 #include "Intracranial.h"
00056
00057
00058
00059
00060
00061
00062
00063 #define USE_QUIET
00064
00065 static char * anat_filename = NULL;
00066 static char * prefix_filename = NULL;
00067 static Boolean write_mask = FALSE;
00068 static Boolean quiet = FALSE;
00069 static char * commandline = NULL ;
00070 static Boolean nosmooth = FALSE;
00071
00072
00073
00074
00075
00076
00077
00078 void SI_error (char * message)
00079 {
00080 fprintf (stderr, "\n");
00081 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00082 exit(1);
00083 }
00084
00085
00086
00087
00088
00089
00090
00091 #include "estpdf3.c"
00092 #include "Intracranial.c"
00093
00094
00095
00096
00097
00098
00099
00100 void display_help_menu()
00101 {
00102 printf
00103 (
00104 "3dIntracranial - performs automatic segmentation of intracranial region.\n"
00105 " \n"
00106 " This program will strip the scalp and other non-brain tissue from a \n"
00107 " high-resolution T1 weighted anatomical dataset. \n"
00108 " \n"
00109 "----------------------------------------------------------------------- \n"
00110 " \n"
00111 "Usage: \n"
00112 "----- \n"
00113 " \n"
00114 "3dIntracranial \n"
00115 " -anat filename => Filename of anat dataset to be segmented \n"
00116 " \n"
00117 " [-min_val a] => Minimum voxel intensity limit \n"
00118 " Default: Internal PDF estimate for lower bound \n"
00119 " \n"
00120 " [-max_val b] => Maximum voxel intensity limit \n"
00121 " Default: Internal PDF estimate for upper bound \n"
00122 " \n"
00123 " [-min_conn m] => Minimum voxel connectivity to enter \n"
00124 " Default: m=4 \n"
00125 " \n"
00126 " [-max_conn n] => Maximum voxel connectivity to leave \n"
00127 " Default: n=2 \n"
00128 " \n"
00129 " [-nosmooth] => Suppress spatial smoothing of segmentation mask \n"
00130 " \n"
00131 " [-mask] => Generate functional image mask (complement) \n"
00132 " Default: Generate anatomical image \n"
00133 " \n"
00134 " [-quiet] => Suppress output to screen \n"
00135 " \n"
00136 " -prefix pname => Prefix name for file to contain segmented image \n"
00137 " \n"
00138 " ** NOTE **: The newer program 3dSkullStrip will probably give \n"
00139 " better segmentation results than 3dIntracranial! \n"
00140
00141 "----------------------------------------------------------------------- \n"
00142 " \n"
00143 "Examples: \n"
00144 "-------- \n"
00145 " \n"
00146 " 3dIntracranial -anat elvis+orig -prefix elvis_strip \n"
00147 " \n"
00148 " 3dIntracranial -min_val 30 -max_val 350 -anat elvis+orig -prefix strip\n"
00149 " \n"
00150 " 3dIntracranial -nosmooth -quiet -anat elvis+orig -prefix elvis_strip \n"
00151 " \n"
00152 "----------------------------------------------------------------------- \n"
00153 );
00154
00155 exit(0);
00156 }
00157
00158
00159
00160
00161
00162
00163
00164 void get_options
00165 (
00166 int argc,
00167 char ** argv
00168 )
00169
00170 {
00171 int nopt = 1;
00172 int ival;
00173 float fval;
00174 char message[MAX_STRING_LENGTH];
00175
00176
00177
00178 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
00179
00180
00181
00182 AFNI_logger (PROGRAM_NAME,argc,argv);
00183
00184
00185
00186 while (nopt < argc )
00187 {
00188
00189
00190 if (strncmp(argv[nopt], "-anat", 5) == 0)
00191 {
00192 nopt++;
00193 if (nopt >= argc) SI_error ("need argument after -anat ");
00194 anat_filename = malloc (sizeof(char) * MAX_STRING_LENGTH);
00195 MTEST (anat_filename);
00196 strcpy (anat_filename, argv[nopt]);
00197
00198 anat = THD_open_dataset (anat_filename);
00199 if (!ISVALID_3DIM_DATASET (anat))
00200 {
00201 sprintf (message, "Can't open dataset: %s\n", anat_filename);
00202 SI_error (message);
00203 }
00204 THD_load_datablock (anat->dblk);
00205 if (DSET_ARRAY(anat,0) == NULL)
00206 {
00207 sprintf (message, "Can't access data: %s\n", anat_filename);
00208 SI_error (message);
00209 }
00210
00211
00212
00213
00214 if( DSET_BRICK_TYPE(anat,0) == MRI_byte ){
00215 THD_3dim_dataset *qset ;
00216 register byte *bar ; register short *sar ;
00217 register int ii,nvox ;
00218 fprintf(stderr,"++ WARNING: converting input dataset from byte to short\n") ;
00219 qset = EDIT_empty_copy(anat) ;
00220 nvox = DSET_NVOX(anat) ;
00221 bar = (byte *) DSET_ARRAY(anat,0) ;
00222 sar = (short *)malloc(sizeof(short)*nvox) ;
00223 for( ii=0 ; ii < nvox ; ii++ ) sar[ii] = (short) bar[ii] ;
00224 EDIT_substitute_brick( qset , 0 , MRI_short , sar ) ;
00225 DSET_delete(anat) ; anat = qset ;
00226 }
00227
00228 nopt++;
00229 continue;
00230 }
00231
00232
00233
00234 if (strncmp(argv[nopt], "-min_val", 8) == 0)
00235 {
00236 nopt++;
00237 if (nopt >= argc) SI_error ("need argument after -min_val ");
00238 sscanf (argv[nopt], "%f", &fval);
00239 if (fval < 0.0)
00240 SI_error ("illegal argument after -min_val ");
00241 min_val_float = fval;
00242 min_val_int = 0;
00243 nopt++;
00244 continue;
00245 }
00246
00247
00248
00249 if (strncmp(argv[nopt], "-max_val", 8) == 0)
00250 {
00251 nopt++;
00252 if (nopt >= argc) SI_error ("need argument after -max_val ");
00253 sscanf (argv[nopt], "%f", &fval);
00254 if (fval < 0.0)
00255 SI_error ("illegal argument after -max_val ");
00256 max_val_float = fval;
00257 max_val_int = 0;
00258 nopt++;
00259 continue;
00260 }
00261
00262
00263
00264 if (strncmp(argv[nopt], "-min_conn", 9) == 0)
00265 {
00266 nopt++;
00267 if (nopt >= argc) SI_error ("need argument after -min_conn ");
00268 sscanf (argv[nopt], "%d", &ival);
00269 if ((ival < 1) || (ival > 7))
00270 SI_error ("illegal argument after -min_conn ");
00271 min_conn_int = ival;
00272 nopt++;
00273 continue;
00274 }
00275
00276
00277
00278 if (strncmp(argv[nopt], "-max_conn", 9) == 0)
00279 {
00280 nopt++;
00281 if (nopt >= argc) SI_error ("need argument after -max_conn ");
00282 sscanf (argv[nopt], "%d", &ival);
00283 if ((ival < -1) || (ival > 5))
00284 SI_error ("illegal argument after -max_conn ");
00285 max_conn_int = ival;
00286 nopt++;
00287 continue;
00288 }
00289
00290
00291
00292 if (strcmp(argv[nopt], "-nosmooth") == 0)
00293 {
00294 nosmooth = TRUE;
00295 nopt++;
00296 continue;
00297 }
00298
00299
00300
00301 if (strncmp(argv[nopt], "-mask", 5) == 0)
00302 {
00303 write_mask = TRUE;
00304 nopt++;
00305 continue;
00306 }
00307
00308
00309
00310 if (strncmp(argv[nopt], "-quiet", 6) == 0)
00311 {
00312 quiet = TRUE;
00313 nopt++;
00314 continue;
00315 }
00316
00317
00318
00319 if (strncmp(argv[nopt], "-prefix", 7) == 0)
00320 {
00321 nopt++;
00322 if (nopt >= argc) SI_error ("need argument after -prefix ");
00323 prefix_filename = malloc (sizeof(char) * MAX_STRING_LENGTH);
00324 MTEST (prefix_filename);
00325 strcpy (prefix_filename, argv[nopt]);
00326 nopt++;
00327 continue;
00328 }
00329
00330
00331
00332 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
00333 SI_error (message);
00334
00335 }
00336
00337
00338 }
00339
00340
00341
00342
00343
00344
00345
00346 void check_one_output_file
00347 (
00348 char * filename
00349 )
00350
00351 {
00352 char message[MAX_STRING_LENGTH];
00353 THD_3dim_dataset * new_dset=NULL;
00354 int ierror;
00355
00356
00357
00358 new_dset = EDIT_empty_copy ( anat );
00359
00360
00361 ierror = EDIT_dset_items( new_dset ,
00362 ADN_prefix , filename ,
00363 ADN_label1 , filename ,
00364 ADN_self_name , filename ,
00365 ADN_none ) ;
00366
00367 if( ierror > 0 )
00368 {
00369 sprintf (message,
00370 "*** %d errors in attempting to create output dataset!\n",
00371 ierror);
00372 SI_error (message);
00373 }
00374
00375 if( THD_is_file(new_dset->dblk->diskptr->header_name) )
00376 {
00377 sprintf (message,
00378 "Output dataset file %s already exists"
00379 "--cannot continue!\a\n",
00380 new_dset->dblk->diskptr->header_name);
00381 SI_error (message);
00382 }
00383
00384
00385 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00386
00387 }
00388
00389
00390
00391
00392
00393
00394
00395 void initialize_program
00396 (
00397 int argc,
00398 char ** argv
00399 )
00400
00401 {
00402 float parameters [DIMENSION];
00403 int nxyz;
00404 short * sfim;
00405 int ixyz, icount;
00406 short * rfim;
00407 int lower_cutoff = 25;
00408
00409
00410
00411 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
00412
00413
00414
00415 get_options (argc, argv);
00416
00417
00418
00419 verify_inputs();
00420
00421
00422
00423 if (anat == NULL) SI_error ("Unable to read anat dataset");
00424 nxyz = DSET_NX(anat) * DSET_NY(anat) * DSET_NZ(anat);
00425 sfim = (short *) DSET_BRICK_ARRAY(anat,0) ;
00426 if (sfim == NULL) SI_error ("Unable to read anat dataset");
00427 rfim = (short *) malloc (sizeof(short) * nxyz); MTEST (rfim);
00428
00429
00430
00431 icount = 0;
00432 for (ixyz = 0; ixyz < nxyz; ixyz++)
00433 if (sfim[ixyz] > lower_cutoff)
00434 {
00435 rfim[icount] = sfim[ixyz];
00436 icount++;
00437 }
00438 if (! quiet)
00439 printf ("%d voxels above lower cutoff = %d \n", icount, lower_cutoff);
00440
00441
00442
00443 if (min_val_int || max_val_int) estpdf_short (icount, rfim, parameters);
00444 if (min_val_int) min_val_float = parameters[4] - 2.0*parameters[5];
00445 if (max_val_int) max_val_float = parameters[7] + 2.0*parameters[8];
00446
00447
00448 if (! quiet)
00449 {
00450 printf ("\n");
00451 printf ("Control inputs: \n");
00452 printf ("anat filename = %s \n", anat_filename);
00453 printf ("min value = %f \n", min_val_float);
00454 printf ("max value = %f \n", max_val_float);
00455 printf ("min conn = %d \n", min_conn_int);
00456 printf ("max conn = %d \n", max_conn_int);
00457 printf ("prefix filename = %s \n", prefix_filename);
00458 }
00459
00460
00461 }
00462
00463
00464
00465
00466
00467
00468
00469 int auto_initialize (short ** cv)
00470
00471 {
00472 int nx, ny, nz, nxy, nxyz, ixyz;
00473
00474
00475
00476 nx = DSET_NX(anat); ny = DSET_NY(anat); nz = DSET_NZ(anat);
00477 nxy = nx*ny; nxyz = nxy*nz;
00478
00479
00480
00481 *cv = (short *) malloc (nxyz * sizeof(short));
00482 MTEST (*cv);
00483 for (ixyz = 0; ixyz < nxyz; ixyz++)
00484 (*cv)[ixyz] = 0;
00485
00486
00487
00488 srand48 (1234567);
00489
00490
00491 return (1);
00492 }
00493
00494
00495
00496
00497
00498
00499
00500 void target_into_dataset
00501 (
00502 short * cv
00503 )
00504
00505 {
00506 short * anat_data = NULL;
00507 int nx, ny, nz, nxy, nxyz, ixyz;
00508
00509
00510
00511 anat_data = (short *) DSET_BRICK_ARRAY(anat,0) ;
00512 nx = DSET_NX(anat); ny = DSET_NY(anat); nz = DSET_NZ(anat);
00513 nxy = nx*ny; nxyz = nxy*nz;
00514
00515 if (! write_mask)
00516 {
00517
00518 for (ixyz = 0; ixyz < nxyz; ixyz++)
00519 {
00520 if (cv[ixyz]) cv[ixyz] = 0;
00521 else cv[ixyz] = anat_data[ixyz];
00522 }
00523 }
00524
00525
00526
00527 THD_delete_3dim_dataset (anat, False); anat = NULL;
00528
00529
00530 return;
00531 }
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543 void write_afni_data
00544 (
00545 short * cv
00546 )
00547
00548 {
00549 int nxyz;
00550 int ii;
00551 THD_3dim_dataset * dset=NULL;
00552 THD_3dim_dataset * new_dset=NULL;
00553 int ierror;
00554 int ibuf[32];
00555 float fbuf[MAX_STAT_AUX];
00556 float fimfac;
00557 int output_datum;
00558 char * filename;
00559
00560
00561
00562 filename = prefix_filename;
00563 dset = THD_open_dataset (anat_filename);
00564 nxyz = DSET_NX(dset) * DSET_NY(dset) * DSET_NZ(dset);
00565
00566
00567
00568 new_dset = EDIT_empty_copy( dset ) ;
00569
00570
00571
00572 tross_Copy_History( dset , new_dset ) ;
00573 if( commandline != NULL )
00574 tross_Append_History( new_dset , commandline ) ;
00575
00576
00577
00578 THD_delete_3dim_dataset (dset, False); dset = NULL ;
00579
00580
00581 output_datum = MRI_short ;
00582
00583 ibuf[0] = output_datum ;
00584
00585 if (write_mask)
00586 {
00587 int func_type = FUNC_FIM_TYPE;
00588 ierror = EDIT_dset_items( new_dset ,
00589 ADN_prefix , filename ,
00590 ADN_label1 , filename ,
00591 ADN_self_name , filename ,
00592 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
00593 ADN_func_type , func_type ,
00594 ADN_nvals , FUNC_nvals[func_type] ,
00595 ADN_datum_array , ibuf ,
00596 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
00597 ADN_none ) ;
00598 }
00599 else
00600 {
00601 ierror = EDIT_dset_items( new_dset ,
00602 ADN_prefix , filename ,
00603 ADN_label1 , filename ,
00604 ADN_self_name , filename ,
00605 ADN_datum_array , ibuf ,
00606 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
00607 ADN_none ) ;
00608 }
00609
00610
00611 if( ierror > 0 ){
00612 fprintf(stderr,
00613 "*** %d errors in attempting to create output dataset!\n",
00614 ierror ) ;
00615 exit(1) ;
00616 }
00617
00618
00619 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
00620 fprintf(stderr,
00621 "*** Output dataset file %s already exists--cannot continue!\a\n",
00622 new_dset->dblk->diskptr->header_name ) ;
00623 exit(1) ;
00624 }
00625
00626
00627
00628 mri_fix_data_pointer (cv, DSET_BRICK(new_dset,0));
00629 fimfac = 1.0;
00630
00631
00632
00633 if (! quiet)
00634 {
00635 if (write_mask) printf("\nWriting functional (mask) dataset: ");
00636 else printf ("\nWriting anatomical dataset: ");
00637 printf("%s\n", DSET_BRIKNAME(new_dset) ) ;
00638
00639 printf("NOTE: You may get better results by trying\n"
00640 " 3dSkullStrip -input %s -prefix %s\n" ,
00641 anat_filename , prefix_filename ) ;
00642 }
00643
00644
00645 for( ii=0 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ;
00646 (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ;
00647
00648 fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ;
00649 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
00650
00651 THD_load_statistics( new_dset ) ;
00652 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
00653
00654
00655
00656 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00657
00658 }
00659
00660
00661
00662
00663
00664
00665
00666 void SEGMENT_auto ()
00667 {
00668 short * cv = NULL;
00669 int ok;
00670
00671
00672
00673 ok = auto_initialize (&cv);
00674 if (! ok) return;
00675
00676
00677
00678 segment_volume (cv);
00679
00680
00681 connectivity_tests (cv);
00682
00683
00684 target_into_dataset (cv);
00685
00686
00687 write_afni_data (cv);
00688
00689 return ;
00690
00691 }
00692
00693
00694
00695
00696
00697
00698
00699 int main
00700 (
00701 int argc,
00702 char ** argv
00703 )
00704
00705 {
00706 int ii ;
00707
00708
00709
00710 for( ii=1 ; ii < argc ; ii++ )
00711 if( strncmp(argv[ii],"-quiet",6) == 0 ) break ;
00712 if( ii == argc ){
00713 printf ("\n\n");
00714 printf ("Program: %s \n", PROGRAM_NAME);
00715 printf ("Author: %s \n", PROGRAM_AUTHOR);
00716 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
00717 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
00718 printf ("\n");
00719 }
00720
00721 mainENTRY("3dIntracranial:main") ; machdep() ;
00722
00723
00724
00725 initialize_program (argc, argv);
00726
00727
00728
00729 SEGMENT_auto();
00730
00731 exit(0);
00732 }
00733
00734
00735