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