00471 {
00472   int nx, ny, nz, nxy, nxyz, ixyz;       /* voxel counters */
00473 
00474 
00475   /*----- Initialize local variables -----*/
00476   nx = DSET_NX(anat);   ny = DSET_NY(anat);   nz = DSET_NZ(anat);
00477   nxy = nx*ny;   nxyz = nxy*nz;
00478 
00479 
00480   /*----- Initialize voxel classification indicators -----*/
00481   *cv = (short *) malloc (nxyz * sizeof(short));
00482   MTEST (*cv);
00483   for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00484     (*cv)[ixyz] = 0;
00485 
00486 
00487   /*----- Initialize random number generator -----*/
00488   srand48 (1234567);
00489 
00490 
00491   return (1);
00492 }

00351 {
00352   char message[MAX_STRING_LENGTH];    /* error message */
00353   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
00354   int ierror;                         /* number of errors in editing data */
00355 
00356   
00357   /*----- make an empty copy of input dataset -----*/
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   /*----- deallocate memory -----*/   
00385   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00386   
00387 }

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 }

00170 {
00171   int nopt = 1;                     /* input option argument counter */
00172   int ival;                         /* integer input */
00173   float fval;                       /* float input */
00174   char message[MAX_STRING_LENGTH];  /* error message */
00175 
00176 
00177   /*----- does user request help menu? -----*/
00178   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();  
00179    
00180   
00181   /*----- add to program log -----*/
00182   AFNI_logger (PROGRAM_NAME,argc,argv); 
00183 
00184 
00185   /*----- main loop over input options -----*/
00186   while (nopt < argc )
00187     {
00188 
00189       /*-----   -anat filename   -----*/
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           /** RWCox [05 Dec 2002]
00212               If input is a byte dataset, make a short copy of it. **/
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       /*-----   -min_val a  -----*/
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       /*-----   -max_val b  -----*/
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       /*-----   -min_conn m  -----*/
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       /*-----   -max_conn n  -----*/
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       /*-----   -nosmooth -----*/
00292       if (strcmp(argv[nopt], "-nosmooth") == 0)
00293         {
00294           nosmooth = TRUE;
00295           nopt++;
00296           continue;
00297         }
00298 
00299 
00300       /*-----   -mask  -----*/
00301       if (strncmp(argv[nopt], "-mask", 5) == 0)
00302         {
00303           write_mask = TRUE;
00304           nopt++;
00305           continue;
00306         }
00307 
00308 
00309       /*-----   -quiet  -----*/
00310       if (strncmp(argv[nopt], "-quiet", 6) == 0)
00311         {
00312           quiet = TRUE;
00313           nopt++;
00314           continue;
00315         }
00316 
00317 
00318       /*-----   -prefix prefixname   -----*/
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       /*----- unknown command -----*/
00332       sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
00333       SI_error (message);
00334       
00335     }
00336 
00337   
00338 }

00401 {
00402   float parameters [DIMENSION];    /* parameters for PDF estimation */
00403   int nxyz;
00404   short * sfim;
00405   int ixyz, icount;
00406   short * rfim;
00407   int lower_cutoff = 25;
00408 
00409 
00410   /*----- save command line for history notes -----*/
00411   commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
00412 
00413 
00414   /*----- Get operator inputs -----*/
00415   get_options (argc, argv);
00416 
00417 
00418   /*----- Verify that inputs are acceptable -----*/
00419   verify_inputs();
00420 
00421 
00422   /*----- Initialize local variables -----*/
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   /*----- Just use voxels whose intensity is above the lower cutoff -----*/
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   /*----- Get PDF estimate and set voxel intensity limits -----*/
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 }

00705 {
00706   int ii ;
00707 
00708   /*----- Identify software -----*/
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   /*----- Program initialization -----*/
00725   initialize_program (argc, argv);
00726 
00727   
00728   /*----- Perform automatic segmentation -----*/
00729   SEGMENT_auto();
00730 
00731   exit(0);
00732 }

00667 {
00668   short * cv = NULL;           /* volume with 1's at target voxel locations */
00669   int ok;                      /* Boolean for successful operation */
00670 
00671 
00672   /*----- Perform initialization for automatic segmentation algorithm -----*/
00673   ok = auto_initialize (&cv);
00674   if (! ok)  return;
00675 
00676 
00677   /*----- Segment intracranial voxels -----*/
00678   segment_volume (cv);
00679 
00680   /*----- Perform voxel connectivity tests -----*/
00681   connectivity_tests (cv);
00682 
00683   /*----- Put estimated target structure into output dataset -----*/
00684   target_into_dataset (cv);
00685 
00686   /*----- Write out the segmented dataset -----*/
00687   write_afni_data (cv);
00688 
00689   return ;
00690 
00691 }

00079 {
00080   fprintf (stderr, "\n");
00081   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00082   exit(1);
00083 }

00505 {
00506   short * anat_data  = NULL;               /* data from anatomical image */
00507   int nx, ny, nz, nxy, nxyz, ixyz;         /* voxel counters */
00508  
00509 
00510   /*----- Initialize local variables -----*/
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       /*----- Reset to zero those voxels which lie outside the brain -----*/
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   /*----- deallocate memory -----*/   
00527   THD_delete_3dim_dataset (anat, False);   anat = NULL;
00528 
00529 
00530   return;
00531 }

00548 {
00549   int nxyz;                           /* number of voxels */
00550   int ii;                             /* voxel index */
00551   THD_3dim_dataset * dset=NULL;       /* input afni data set pointer */
00552   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
00553   int ierror;                         /* number of errors in editing data */
00554   int ibuf[32];                       /* integer buffer */
00555   float fbuf[MAX_STAT_AUX];           /* float buffer */
00556   float fimfac;                       /* scale factor for short data */
00557   int output_datum;                   /* data type for output data */
00558   char * filename;                    /* prefix filename for output */
00559 
00560 
00561   /*----- initialize local variables -----*/
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   /*-- make an empty copy of this dataset, for eventual output --*/
00568   new_dset = EDIT_empty_copy( dset ) ;
00569 
00570 
00571   /*----- Record history of dataset -----*/
00572   tross_Copy_History( dset , new_dset ) ;
00573   if( commandline != NULL )
00574      tross_Append_History( new_dset , commandline ) ;
00575 
00576   
00577   /*----- deallocate memory -----*/   
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   /*----- attach bricks to new data set -----*/
00628   mri_fix_data_pointer (cv, DSET_BRICK(new_dset,0)); 
00629   fimfac = 1.0;
00630   
00631 
00632   /*----- write afni data set -----*/
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   /*----- deallocate memory -----*/   
00656   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00657   
00658 }