00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   This file contains routines used by plug_intracranial and 3dIntracranial.
00010 
00011   File:    Intracranial.c
00012   Author:  B. Douglas Ward
00013   Date:    04 June 1999
00014 
00015   Mod:     Correction to initialization in center of mass calculation.
00016   Date:    11 February 2000
00017 
00018   Mod:     Added option to suppress spatial smoothing of segmentation mask.
00019   Date:    03 December 2001
00020 */
00021 
00022 
00023 /*---------------------------------------------------------------------------*/
00024 /*
00025   Verify that inputs are acceptable.
00026 */
00027 
00028 int verify_inputs()
00029 {
00030   char message[MAX_STRING_LENGTH];    /* error message */
00031 
00032 
00033   /*----- Check for required datasets -----*/
00034   if (anat == NULL)
00035     { 
00036       sprintf (message, "No anatomical image");
00037       SI_error (message); 
00038       return (0);
00039     }
00040   if (! ISANAT(anat))
00041     { 
00042       sprintf (message, "Input dataset must be anatomical image");
00043       SI_error (message); 
00044       return (0);
00045     }
00046   if (DSET_BRICK_TYPE(anat,0) != MRI_short)
00047     {
00048       sprintf (message, "Anatomical image must have data type: short integer");
00049       SI_error (message);
00050       return (0);
00051     }
00052 
00053 #ifdef PLUG_INTRACRANIAL
00054   if (dset == NULL)
00055     { 
00056       sprintf (message, "No output dataset");
00057       SI_error (message);
00058       return (0);
00059     }
00060   if (DSET_BRICK_TYPE(dset,0) != MRI_short)
00061     {
00062       sprintf (message, "Output dataset must have data type: short integer");
00063       SI_error (message);
00064       return (0);
00065     }
00066 #else
00067   if (prefix_filename == NULL)
00068     {
00069       sprintf (message, "Must specify output prefix filename");
00070       SI_error (message);
00071       return (0);
00072     }
00073 
00074   /*----- Check whether output file already exists -----*/
00075   check_one_output_file (prefix_filename);
00076 #endif
00077 
00078 
00079   /*----- Check compatibility of datasets -----*/
00080 #ifdef PLUG_INTRACRANIAL
00081   if (! EQUIV_DATAXES (anat->daxes, dset->daxes))
00082     {
00083       sprintf (message, 
00084                "Anatomical image and output dataset are incompatible");
00085       SI_error (message);
00086       return (0);
00087     }
00088 #endif
00089     
00090 
00091   /*----- Check allowed range of intensity values -----*/
00092   if (min_val_float >= max_val_float)
00093     {
00094       sprintf (message, "min_val >= max_val ???");
00095       SI_error (message);
00096       return (0);
00097     }
00098 
00099   return (1);
00100 
00101 }
00102 
00103 
00104 /*---------------------------------------------------------------------------*/
00105 /*
00106    Flood filling a byte array:
00107      nx = 1st dimension
00108      ny = 2nd dimension
00109      ix = start point
00110      jy = end point
00111      ar = array, with 0's everwhere except 1's as barriers to flooding
00112 
00113    All filled points (starting with ix,jy) will get the value 2.
00114 
00115    This routine was adapted from:  plug_drawdset.
00116 */
00117 
00118 void draw_2dfiller( int nx , int ny , int ix , int jy , byte * ar )
00119 {
00120    int ii,jj , ip,jp , num ;
00121 
00122 #define AR(i,j) ar[(i)+(j)*nx]
00123 
00124 
00125    /*----- Test for early termination -----*/
00126    if (AR(ix,jy) != 0)  return;
00127 
00128 
00129    /* fill out in cross from 1st point */
00130 
00131    ip = ix ; jp = jy ; AR(ip,jp) = 2 ;
00132 
00133    for( ii=ip+1; ii < nx && AR(ii,jp) == 0; ii++ ) AR(ii,jp) = 2;
00134    for( ii=ip-1; ii >= 0 && AR(ii,jp) == 0; ii-- ) AR(ii,jp) = 2;
00135    for( jj=jp+1; jj < ny && AR(ip,jj) == 0; jj++ ) AR(ip,jj) = 2;
00136    for( jj=jp-1; jj >= 0 && AR(ip,jj) == 0; jj-- ) AR(ip,jj) = 2;
00137 
00138    /* brute force repetition of the cross technique */
00139 
00140    do {
00141       num = 0 ;
00142       for( jp=0 ; jp < ny ; jp++ ){
00143          for( ip=0 ; ip < nx ; ip++ ){
00144             if( AR(ip,jp) == 2 ){
00145                for( ii=ip+1; ii < nx && AR(ii,jp) == 0; ii++ ){ AR(ii,jp) = 2; num++; }
00146                for( ii=ip-1; ii >= 0 && AR(ii,jp) == 0; ii-- ){ AR(ii,jp) = 2; num++; }
00147                for( jj=jp+1; jj < ny && AR(ip,jj) == 0; jj++ ){ AR(ip,jj) = 2; num++; }
00148                for( jj=jp-1; jj >= 0 && AR(ip,jj) == 0; jj-- ){ AR(ip,jj) = 2; num++; }
00149             }
00150          }
00151       }
00152    } while( num > 0 ) ;
00153 
00154    return ;
00155 }
00156 
00157 
00158 /*---------------------------------------------------------------------------*/
00159 /*
00160   Find the center of mass of the brain image.
00161 */
00162 
00163 void center_of_mass (int * cx, int * cy, int * cz)
00164 
00165 {
00166   int nx, ny, nz, nxy, nxyz;               /* dataset dimensions in voxels */
00167   int ix, jy, kz, ixyz;                    /* voxel indices */
00168   short * anat_data  = NULL;               /* data from anatomical image */
00169   float f, sum, sumx, sumy, sumz;          /* sums for calculating means */
00170 
00171 
00172   /*----- Progress report -----*/
00173   if (! quiet)  printf ("\nCalculating center of mass \n");
00174 
00175 
00176   /*----- Initialize local variables -----*/
00177   anat_data  = (short *) DSET_BRICK_ARRAY(anat,0) ;
00178   nx = DSET_NX(anat);   ny = DSET_NY(anat);   nz = DSET_NZ(anat);
00179   nxy = nx * ny;   nxyz = nxy * nz;
00180 
00181 
00182   /*----- Sum over all voxels -----*/
00183   sum = 0.0;   sumx = 0.0;   sumy = 0.0;   sumz = 0.0;
00184   for (kz = 0;  kz < nz;  kz++)
00185     for (jy = 0;  jy < ny;  jy++)
00186       for (ix = 0;  ix < nx;  ix++)
00187         {
00188           ixyz = ix + jy*nx + kz*nxy;
00189           f = anat_data[ixyz];
00190           sum += f;   sumx += ix*f;   sumy += jy*f;   sumz += kz*f;
00191         }
00192 
00193 
00194   /*----- Calculate center of mass -----*/
00195   *cx = sumx / sum;
00196   *cy = sumy / sum;
00197   *cz = sumz / sum;
00198 
00199   if (! quiet)  printf ("Center of mass:   ix = %d   jy = %d   kz = %d  \n", 
00200                         *cx, *cy, *cz);
00201  
00202 }
00203 
00204 
00205 /*---------------------------------------------------------------------------*/
00206 /*
00207   Segment the slices perpendicular to the x-axis.
00208 */
00209 
00210 void segment_x_slices
00211 (
00212   short * cv,                  /* volume with 1's at non-brain locations */
00213   int cx,                      /* center slice location */
00214   Boolean axial_slice          /* true if x-slices are axial slices */
00215 )
00216 
00217 {
00218   const int MAXNPTS = 1000;    /* max. number of random initial test points */
00219   const int MAXFOUND = 10;     /* max. number of initial search points */
00220 
00221   int nx, ny, nz, nxy, nxz, nyz, nxyz;     /* dataset dimensions in voxels */
00222   int ix, jy, kz, ixy, ixz, iyz, ixyz;     /* voxel indices */
00223 
00224   int m;                       /* slice index */
00225   short * anat_data  = NULL;   /* data from anatomical image */
00226   float z;                     /* voxel gray-scale intensity */
00227   int nmpts, found;            /* random search counters */
00228   byte * slice = NULL;         /* current slice brain mask */
00229   byte * pslice = NULL;        /* previous slice non-brain mask */
00230 
00231   int count;                   /* count of brain voxels in current slice */
00232   int maxcount = 0;            /* maximum number of brain voxels in a slice */
00233   float threshold = 0.05;      /* stop if count falls below 5% of maximum */
00234 
00235 
00236   /*----- Progress report -----*/
00237   if (! quiet)  printf ("\nSegmenting slices perpendicular to x-axis \n");
00238 
00239 
00240   /*----- Initialize local variables -----*/
00241   anat_data  = (short *) DSET_BRICK_ARRAY(anat,0) ;
00242   nx = DSET_NX(anat);   ny = DSET_NY(anat);   nz = DSET_NZ(anat);
00243   nxy = nx*ny;   nxz = nx*nz;   nyz = ny*nz;   nxyz = nxy*nz;
00244 
00245 
00246   /*----- Allocate memory -----*/
00247   slice = (byte *)  malloc (nyz * sizeof(byte));
00248   pslice = (byte *) malloc (nyz * sizeof(byte));
00249 
00250 
00251   /*----- Loop over x-slices -----*/
00252   m = 0; 
00253   while (m <= nx)
00254     {
00255       /*----- Start from center slice -----*/
00256       if (m <= cx)  ix = cx - m;
00257       else          ix = m-1;
00258       
00259 
00260       /*----- Initialize mask  -----*/
00261       if (ix == cx)  
00262         {
00263           if (axial_slice)
00264             for (iyz = 0;  iyz < nyz;  iyz++)   pslice[iyz] = 0;
00265           else
00266             {
00267               for (kz = 0;  kz < nz;  kz++)
00268                 for (jy = 0;  jy < ny;  jy++)
00269                   {
00270                     iyz = jy + kz*ny;
00271                     ixyz = ix + jy*nx + kz*nxy;
00272                     if (cv[ixyz])   pslice[iyz] = 2;
00273                     else            pslice[iyz] = 0;
00274                   }
00275             }
00276         }
00277 
00278 
00279       /*----- Examine each voxel for possible inclusion in brain -----*/
00280       for (kz = 0;  kz < nz;  kz++)
00281         for (jy = 0;  jy < ny;  jy++)
00282           {
00283             iyz = jy + kz*ny;
00284             ixyz = ix + jy*nx + kz*nxy;
00285             z = anat_data[ixyz];
00286 
00287             if (pslice[iyz] == 2)        slice[iyz] = 1;
00288             else if (z < min_val_float)  slice[iyz] = 1;
00289             else if (z > max_val_float)  slice[iyz] = 1;
00290             else                         slice[iyz] = 0;
00291           }
00292 
00293 
00294       /*----- Fill brain voxels from center out -----*/
00295       draw_2dfiller (ny, nz, ny/2, nz/2, slice);
00296 
00297 
00298       /*----- Use additional random starting points for robustness -----*/
00299       nmpts = 0;   found = 0;
00300       while ((nmpts < MAXNPTS) && (found < MAXFOUND))
00301         {
00302           nmpts++;
00303           if ((ix == cx) && axial_slice)
00304             {
00305               jy = (3*ny/8) + (int) (rand_uniform(0.0,1.0) * ny/4);
00306               kz = (3*nz/8) + (int) (rand_uniform(0.0,1.0) * nz/4);
00307             }
00308           else
00309             {
00310               jy = (int) (rand_uniform(0.0,1.0) * ny);
00311               kz = (int) (rand_uniform(0.0,1.0) * nz);
00312             }
00313           if (slice[jy+kz*ny] != 1)  found++;
00314           draw_2dfiller (ny, nz, jy, kz, slice);
00315         }
00316 
00317 
00318       /*----- Filled brain voxels become barriers -----*/
00319       for (iyz = 0;  iyz < nyz;  iyz++)
00320         if (slice[iyz] == 2)  pslice[iyz] = 1;
00321         else                  pslice[iyz] = 0;
00322 
00323 
00324       /*----- Now fill non-brain voxels from outside in -----*/
00325       draw_2dfiller (ny, nz, 0,    0,    pslice);
00326       draw_2dfiller (ny, nz, ny-1, 0,    pslice);
00327       draw_2dfiller (ny, nz, 0,    nz-1, pslice);
00328       draw_2dfiller (ny, nz, ny-1, nz-1, pslice);
00329 
00330 
00331       /*----- Save results for this slice -----*/
00332       count = 0;
00333       for (kz = 0;  kz < nz;  kz++)
00334         for (jy = 0;  jy < ny;  jy++)
00335           {
00336             iyz = jy + kz*ny;
00337             ixyz = ix + jy*nx + kz*nxy;
00338             if (pslice[iyz] != 2)
00339               {
00340                 cv[ixyz] = 0;
00341                 count++;
00342               }
00343           }
00344 
00345 
00346       /*----- Slightly increase size of brain mask for next slice -----*/
00347       for (kz = 1;  kz < nz-1;  kz++)
00348         for (jy = 1;  jy < ny-1;  jy++)
00349           {
00350             iyz = jy + kz*ny;
00351             ixyz = ix + jy*nx + kz*nxy;
00352             if (!cv[ixyz])  
00353               {
00354                 pslice[iyz]    = 0;
00355                 pslice[iyz+1]  = 0;
00356                 pslice[iyz-1]  = 0;
00357                 pslice[iyz-ny] = 0;
00358                 pslice[iyz+ny] = 0;
00359               }
00360           }
00361 
00362 
00363       /*----- Examine count of brain voxels for this slice -----*/
00364       if (count > maxcount)  maxcount = count;
00365       if (count < threshold * maxcount)
00366         {
00367           if (m <= cx)  m = cx;
00368           else          m = nx;
00369         }
00370 
00371 
00372       /*----- Increment slice index -----*/
00373       m++;
00374 
00375     }  /* Loop over x-slices */
00376 
00377 
00378   /*----- Release memory -----*/
00379   free (pslice);   pslice = NULL;
00380   free (slice);    slice = NULL;
00381 
00382 
00383   return;
00384 
00385 }
00386 
00387 
00388 /*---------------------------------------------------------------------------*/
00389 /*
00390   Segment the slices perpendicular to the y-axis.
00391 */
00392 
00393 void segment_y_slices
00394 (
00395   short * cv,                  /* volume with 1's at non-brain locations */
00396   int cy,                      /* center slice location */
00397   Boolean axial_slice          /* true if y-slices are axial slices */
00398 )
00399 
00400 {
00401   const int MAXNPTS = 1000;    /* max. number of random initial test points */
00402   const int MAXFOUND = 10;     /* max. number of initial search points */
00403 
00404   int nx, ny, nz, nxy, nxz, nyz, nxyz;     /* dataset dimensions in voxels */
00405   int ix, jy, kz, ixy, ixz, iyz, ixyz;     /* voxel indices */
00406 
00407   int m;                       /* slice index */
00408   short * anat_data  = NULL;   /* data from anatomical image */
00409   float z;                     /* voxel gray-scale intensity */
00410   int nmpts, found;            /* random search counters */
00411   byte * slice = NULL;         /* current slice brain mask */
00412   byte * pslice = NULL;        /* previous slice non-brain mask */
00413 
00414   int count;                   /* count of brain voxels in current slice */
00415   int maxcount = 0;            /* maximum number of brain voxels in a slice */
00416   float threshold = 0.05;      /* stop if count falls below 5% of maximum */
00417 
00418 
00419   /*----- Progress report -----*/
00420   if (! quiet)  printf ("\nSegmenting slices perpendicular to y-axis \n");
00421 
00422 
00423   /*----- Initialize local variables -----*/
00424   anat_data  = (short *) DSET_BRICK_ARRAY(anat,0) ;
00425   nx = DSET_NX(anat);   ny = DSET_NY(anat);   nz = DSET_NZ(anat);
00426   nxy = nx*ny;   nxz = nx*nz;   nyz = ny*nz;   nxyz = nxy*nz;
00427 
00428 
00429   /*----- Allocate memory -----*/
00430   slice = (byte *)  malloc (nxz * sizeof(byte));
00431   pslice = (byte *) malloc (nxz * sizeof(byte));
00432 
00433 
00434   /*----- Loop over y-slices -----*/
00435   m = 0;
00436   while (m <= ny)
00437     {
00438       /*----- Start from center slice -----*/
00439       if (m <= cy)  jy = cy - m;
00440       else          jy = m-1;
00441       
00442 
00443       /*----- Initialize mask  -----*/
00444       if (jy == cy) 
00445         { 
00446           if (axial_slice)
00447             for (ixz = 0;  ixz < nxz;  ixz++)   pslice[ixz] = 0;
00448           else 
00449             { 
00450               for (kz = 0;  kz < nz;  kz++)
00451                 for (ix = 0;  ix < nx;  ix++)
00452                   {
00453                     ixz = ix + kz*nx;
00454                     ixyz = ix + jy*nx + kz*nxy;
00455                     if (cv[ixyz])   pslice[ixz] = 2;
00456                     else            pslice[ixz] = 0;
00457                   }
00458             }
00459         }
00460 
00461 
00462       /*----- Examine each voxel for possible inclusion in brain -----*/
00463       for (kz = 0;  kz < nz;  kz++)
00464         for (ix = 0;  ix < nx;  ix++)
00465           {
00466             ixz = ix + kz*nx;
00467             ixyz = ix + jy*nx + kz*nxy;
00468             z = anat_data[ixyz];
00469 
00470             if (pslice[ixz] == 2)        slice[ixz] = 1;
00471             else if (z < min_val_float)  slice[ixz] = 1;
00472             else if (z > max_val_float)  slice[ixz] = 1;
00473             else                         slice[ixz] = 0;
00474           }
00475 
00476 
00477       /*----- Fill brain voxels from center out -----*/
00478       draw_2dfiller (nx, nz, nx/2, nz/2, slice);
00479 
00480 
00481       /*----- Use additional random starting points for robustness -----*/
00482       nmpts = 0;   found = 0;
00483       while ((nmpts < MAXNPTS) && (found < MAXFOUND))
00484         {
00485           nmpts++;
00486           if ((jy == cy) && axial_slice)
00487             {
00488               ix = (3*nx/8) + (int) (rand_uniform(0.0,1.0) * nx/4);
00489               kz = (3*nz/8) + (int) (rand_uniform(0.0,1.0) * nz/4);
00490             }
00491           else
00492             {
00493               ix = (int) (rand_uniform(0.0,1.0) * nx);
00494               kz = (int) (rand_uniform(0.0,1.0) * nz);
00495             }
00496           if (slice[ix+kz*nx] != 1)  found++;
00497           draw_2dfiller (nx, nz, ix, kz, slice);
00498         }
00499 
00500 
00501       /*----- Filled brain voxels become barriers -----*/
00502       for (ixz = 0;  ixz < nxz;  ixz++)
00503         if (slice[ixz] == 2)  pslice[ixz] = 1;
00504         else                  pslice[ixz] = 0;
00505 
00506 
00507       /*----- Now fill non-brain voxels from outside in -----*/
00508       draw_2dfiller (nx, nz, 0,    0,    pslice);
00509       draw_2dfiller (nx, nz, nx-1, 0,    pslice);
00510       draw_2dfiller (nx, nz, 0,    nz-1, pslice);
00511       draw_2dfiller (nx, nz, nx-1, nz-1, pslice);
00512 
00513 
00514       /*----- Save results for this slice -----*/
00515       count = 0;
00516       for (kz = 0;  kz < nz;  kz++)
00517         for (ix = 0;  ix < nx;  ix++)
00518           {
00519             ixz = ix + kz*nx;
00520             ixyz = ix + jy*nx + kz*nxy;
00521             if (pslice[ixz] != 2)  
00522               {
00523                 cv[ixyz] = 0;
00524                 count++;
00525               }
00526           }
00527 
00528 
00529       /*----- Slightly increase size of brain mask for next slice -----*/
00530       for (kz = 1;  kz < nz-1;  kz++)
00531         for (ix = 1;  ix < nx-1;  ix++)
00532           {
00533             ixz = ix + kz*nx;
00534             ixyz = ix + jy*nx + kz*nxy;
00535             if (!cv[ixyz])  
00536               {
00537                 pslice[ixz]    = 0;
00538                 pslice[ixz+1]  = 0;
00539                 pslice[ixz-1]  = 0;
00540                 pslice[ixz-nx] = 0;
00541                 pslice[ixz+nx] = 0;
00542               }
00543           }
00544 
00545 
00546       /*----- Examine count of brain voxels for this slice -----*/
00547       if (count > maxcount)  maxcount = count;
00548       if (count < threshold * maxcount)
00549         {
00550           if (m <= cy)  m = cy;
00551           else          m = ny;
00552         }
00553       
00554 
00555       /*----- Increment slice index -----*/
00556       m++;
00557 
00558     }  /* Loop over y-slices */
00559 
00560 
00561   /*----- Release memory -----*/
00562   free (pslice);   pslice = NULL;
00563   free (slice);    slice = NULL;
00564 
00565 
00566   return;
00567 
00568 }
00569 
00570 
00571 /*---------------------------------------------------------------------------*/
00572 /*
00573   Segment the slices perpendicular to the z-axis.
00574 */
00575 
00576 void segment_z_slices
00577 (
00578   short * cv,                  /* volume with 1's at non-brain locations */
00579   int cz,                      /* center slice location */
00580   Boolean axial_slice          /* true if z-slices are axial slices */
00581 )
00582 
00583 {
00584   const int MAXNPTS = 1000;    /* max. number of random initial test points */
00585   const int MAXFOUND = 10;     /* max. number of initial search points */
00586 
00587   int nx, ny, nz, nxy, nxz, nyz, nxyz;     /* dataset dimensions in voxels */
00588   int ix, jy, kz, ixy, ixz, iyz, ixyz;     /* voxel indices */
00589 
00590   int m;                       /* slice index */
00591   short * anat_data  = NULL;   /* data from anatomical image */
00592   float z;                     /* voxel gray-scale intensity */
00593   int nmpts, found;            /* random search counters */
00594   byte * slice = NULL;         /* current slice brain mask */
00595   byte * pslice = NULL;        /* previous slice non-brain mask */
00596 
00597   int count;                   /* count of brain voxels in current slice */
00598   int maxcount = 0;            /* maximum number of brain voxels in a slice */
00599   float threshold = 0.05;      /* stop if count falls below 5% of maximum */
00600 
00601 
00602   /*----- Progress report -----*/
00603   if (! quiet)  printf ("\nSegmenting slices perpendicular to z-axis \n");
00604 
00605 
00606   /*----- Initialize local variables -----*/
00607   anat_data  = (short *) DSET_BRICK_ARRAY(anat,0) ;
00608   nx = DSET_NX(anat);   ny = DSET_NY(anat);   nz = DSET_NZ(anat);
00609   nxy = nx*ny;   nxz = nx*nz;   nyz = ny*nz;   nxyz = nxy*nz;
00610 
00611 
00612   /*----- Allocate memory -----*/
00613   slice = (byte *)  malloc (nxy * sizeof(byte));
00614   pslice = (byte *) malloc (nxy * sizeof(byte));
00615 
00616 
00617   /*----- Loop over z-slices -----*/
00618   m = 0;
00619   while (m <= nz)
00620     {
00621       /*----- Start from center slice -----*/
00622       if (m <= cz)  kz = cz - m;
00623       else          kz = m-1;
00624 
00625       
00626       /*----- Initialize mask -----*/
00627       if (kz == cz)  
00628         {
00629           if (axial_slice)
00630             for (ixy = 0;  ixy < nxy;  ixy++)   pslice[ixy] = 0;
00631           else
00632             {
00633               for (ix = 0;  ix < nx;  ix++)
00634                 for (jy = 0;  jy < ny;  jy++)
00635                   {
00636                     ixy = ix + jy*nx;
00637                     ixyz = ix + jy*nx + kz*nxy;
00638                     if (cv[ixyz])   pslice[ixy] = 2;
00639                     else            pslice[ixy] = 0;
00640                   }
00641             }
00642         }
00643 
00644 
00645       /*----- Examine each voxel for possible inclusion in brain -----*/
00646       for (jy = 0;  jy < ny;  jy++)
00647         for (ix = 0;  ix < nx;  ix++)
00648           {
00649             ixy = ix + jy*nx;
00650             ixyz = ix + jy*nx + kz*nxy;
00651             z = anat_data[ixyz];
00652 
00653             if (pslice[ixy] == 2)        slice[ixy] = 1;
00654             else if (z < min_val_float)  slice[ixy] = 1;
00655             else if (z > max_val_float)  slice[ixy] = 1;
00656             else                         slice[ixy] = 0;
00657           }
00658 
00659 
00660       /*----- Fill brain voxels from center out -----*/
00661       draw_2dfiller (nx, ny, nx/2, ny/2, slice);
00662 
00663 
00664       /*----- Use additional random starting points for robustness -----*/
00665       nmpts = 0;   found = 0;
00666       while ((nmpts < MAXNPTS) && (found < MAXFOUND))
00667         {
00668           nmpts++;
00669           if ((kz == cz) && axial_slice)
00670             {
00671               ix = (3*nx/8) + (int) (rand_uniform(0.0,1.0) * nx/4);
00672               jy = (3*ny/8) + (int) (rand_uniform(0.0,1.0) * ny/4);
00673             }
00674           else
00675             {
00676               ix = (int) (rand_uniform(0.0,1.0) * nx);
00677               jy = (int) (rand_uniform(0.0,1.0) * ny);
00678             }
00679           if (slice[ix+jy*nx] != 1)  found++;
00680           draw_2dfiller (nx, ny, ix, jy, slice);
00681         }
00682 
00683 
00684       /*----- Filled brain voxels become barriers -----*/
00685       for (ixy = 0;  ixy < nxy;  ixy++)
00686         if (slice[ixy] == 2)  pslice[ixy] = 1;
00687         else                  pslice[ixy] = 0;
00688 
00689 
00690       /*----- Now fill non-brain voxels from outside in -----*/
00691       draw_2dfiller (nx, ny, 0,    0,    pslice);
00692       draw_2dfiller (nx, ny, nx-1, 0,    pslice);
00693       draw_2dfiller (nx, ny, 0,    ny-1, pslice);
00694       draw_2dfiller (nx, ny, nx-1, ny-1, pslice);
00695 
00696 
00697       /*----- Save results for this slice -----*/
00698       count = 0;
00699       for (jy = 0;  jy < ny;  jy++)
00700         for (ix = 0;  ix < nx;  ix++)
00701           {
00702             ixy = ix + jy*nx;
00703             ixyz = ix + jy*nx + kz*nxy;
00704             if (pslice[ixy] != 2)   
00705               {
00706                 cv[ixyz] = 0;
00707                 count++;
00708               }
00709           }
00710 
00711 
00712       /*----- Slightly increase size of brain mask for next slice -----*/
00713       for (jy = 1;  jy < ny-1;  jy++)
00714         for (ix = 1;  ix < nx-1;  ix++)
00715           {
00716             ixy = ix + jy*nx;
00717             ixyz = ix + jy*nx + kz*nxy;
00718             if (!cv[ixyz])  
00719               {
00720                 pslice[ixy]    = 0;
00721                 pslice[ixy+1]  = 0;
00722                 pslice[ixy-1]  = 0;
00723                 pslice[ixy-nx] = 0;
00724                 pslice[ixy+nx] = 0;
00725               }
00726           }
00727 
00728 
00729       /*----- Examine count of brain voxels for this slice -----*/
00730       if (count > maxcount)  maxcount = count;
00731       if (count < threshold * maxcount)
00732         {
00733           if (m <= cz)  m = cz;
00734           else          m = nz;
00735         }
00736 
00737 
00738       /*----- Increment slice index -----*/
00739       m++;
00740 
00741     }  /* Loop over z-slices */
00742 
00743 
00744   /*----- Release memory -----*/
00745   free (pslice);   pslice = NULL;
00746   free (slice);    slice = NULL;
00747 
00748 
00749   return;
00750 
00751 }
00752 
00753 
00754 /*---------------------------------------------------------------------------*/
00755 /*
00756   Create envelope for segmented image.
00757 */
00758 
00759 void segment_envelope
00760 (
00761   short * cv,                  /* volume with 1's at non-brain locations */
00762   int cx, int cy, int cz       /* location of center of mass */
00763 )
00764 
00765 {
00766 #define MAXLAT 90
00767 #define MAXLNG 180
00768 
00769   int nx, ny, nz, nxy, nxz, nyz, nxyz;     /* dataset dimensions in voxels */
00770   int ix, jy, kz, ixy, ixz, iyz, ixyz;     /* voxel indices */
00771 
00772   /* 24 Mar 2004: on Mac OS X, the compiler generates incorrect code
00773                   for these large arrays when optimization is turned on;
00774                   instead, use malloc() to get them if AIZE is not defined */
00775 
00776 #undef AIZE
00777 #ifndef DARWIN
00778 # define AIZE
00779 #endif
00780 
00781 #ifdef  AIZE
00782   float radius[MAXLAT][MAXLNG];      /* max. radius vector to brain voxel */
00783   float smradius[MAXLAT][MAXLNG];    /* array of smoothed radius vectors */
00784 #else
00785   float **radius , **smradius ;
00786 #endif
00787 
00788   int ilat, jlat;              /* radius array latitude index */
00789   int ilng, jlng;              /* radius array longitude index */
00790   float deltalat, deltalng;    /* increments in latitude and longitude */
00791   float lat, clat, slat;       /* voxel latitude */
00792   float lng, clng, slng;       /* voxel longitude */
00793   int ir, i, j;                /* indices */
00794   float rxy;                   /* radius to voxel in xy-plane */
00795   float r;                     /* radius to voxel from center of brain */
00796   int nr;                      /* maximum test radius */
00797   float maxr;                  /* maximum radius to voxel inside the brain */
00798 
00799 
00800   /*----- Progress report -----*/
00801   if (! quiet)  printf ("\nCreating envelope for segmented image \n");
00802 
00803 #ifndef AIZE
00804   radius   = (float **) malloc(sizeof(float)*MAXLAT) ;
00805   smradius = (float **) malloc(sizeof(float)*MAXLAT) ;
00806   for (ilat = 0;  ilat < MAXLAT;  ilat++){
00807       radius[ilat] = (float *)malloc(sizeof(float)*MAXLNG) ;
00808     smradius[ilat] = (float *)malloc(sizeof(float)*MAXLNG) ;
00809   }
00810 #endif
00811 
00812 
00813   /*----- Initialize local variables -----*/
00814   nx = DSET_NX(anat);   ny = DSET_NY(anat);   nz = DSET_NZ(anat);
00815   nxy = nx*ny;   nxz = nx*nz;   nyz = ny*nz;   nxyz = nxy*nz;
00816   deltalat = PI / MAXLAT;
00817   deltalng = 2.0 * PI / MAXLNG;
00818   nr = nx;  
00819   if (ny > nr)  nr = ny;
00820   if (nz > nr)  nr = nz;
00821 
00822     
00823   
00824   /*----- Calculate the radius vector -----*/
00825   for (ilat = 0;  ilat < MAXLAT;  ilat++)
00826     {
00827       lat = ilat * deltalat - PI/2;
00828       slat = sin(lat);
00829       clat = cos(lat);
00830       
00831       for (ilng = 0;  ilng < MAXLNG;  ilng++)
00832         {
00833           lng = ilng * deltalng;
00834           clng = cos(lng);
00835           slng = sin(lng);
00836           
00837           radius[ilat][ilng] = 0.0;
00838           for (ir = 0;  ir < nr;  ir++)
00839             {
00840               ix = cx + ir*clat*clng;
00841               jy = cy + ir*clat*slng;
00842               kz = cz + ir * slat;
00843 
00844               if ((ix >= 0) && (ix < nx) && (jy >= 0) && (jy < ny)
00845                   && (kz >= 0) && (kz < nz))
00846                 {
00847                   ixyz = ix + jy*nx + kz*nxy;
00848                   if (! cv[ixyz])  radius[ilat][ilng] = (float) ir; 
00849                 }
00850             }
00851         }
00852     }
00853   
00854   /*----- Smooth the radius vectors -----*/
00855   for (ilat = 0;  ilat < MAXLAT;  ilat++)
00856     {
00857       for (ilng = 0;  ilng < MAXLNG;  ilng++)
00858         {
00859           smradius[ilat][ilng] = 0.0;
00860           
00861           for (i = -1;  i <= 1;  i++)
00862             {
00863               jlat = (ilat + i + MAXLAT) % MAXLAT;
00864               
00865               if (ilat == 0)  jlat = 0;
00866               if (ilat == MAXLAT-1)  jlat = MAXLAT-1;
00867               
00868               for (j = -2;  j <= 2;  j++)
00869                 {
00870                   jlng = (ilng + j + MAXLNG) % MAXLNG;
00871                   smradius[ilat][ilng] += radius[jlat][jlng] / 15.0;
00872                 }
00873             }
00874         }
00875     }
00876 
00877   /*----- Find maximum radius to a brain voxel -----*/
00878   maxr = 0.0;
00879   for (ilat = 0;  ilat < MAXLAT;  ilat++)
00880     for (ilng = 0;  ilng < MAXLNG;  ilng++)    
00881       if (maxr < smradius[ilat][ilng])  maxr = smradius[ilat][ilng];
00882           
00883   /*----- Smooth the brain mask -----*/
00884   for (ix = 0;  ix < nx;  ix++){
00885     for (jy = 0;  jy < ny;  jy++)
00886       {
00887         rxy = hypot ((float) (ix-cx), (float) (jy-cy));
00888         lng = atan2 (jy-cy,  ix-cx);
00889         ilng = (int) ((lng/deltalng) + MAXLNG) % MAXLNG;
00890 
00891         for (kz = 0;  kz < nz;  kz++)
00892           {
00893             ixyz = ix + jy*nx + kz*nxy;
00894 
00895             r = hypot (rxy, (float) (kz-cz));     
00896 
00897             if (r < maxr)
00898               {
00899                 lat = atan2 (kz-cz, rxy);
00900                 ilat = (int) (((lat+PI/2)/deltalat) + MAXLAT) % MAXLAT;
00901                 
00902                 if (r < smradius[ilat][ilng])  cv[ixyz] = 0;
00903               }
00904             
00905           }
00906       }
00907    }
00908 
00909 #ifndef AIZE
00910   for (ilat = 0;  ilat < MAXLAT;  ilat++){
00911     free(  radius[ilat]) ; free(smradius[ilat]) ;
00912   }
00913   free(smradius); free(radius);
00914 #endif
00915 }
00916 
00917 
00918 /*---------------------------------------------------------------------------*/
00919 /*
00920   Segment the intracranial voxels
00921 */
00922 
00923 void segment_volume 
00924 (
00925   short * cv                    /* volume with 1's at non-brain locations */
00926 )
00927 
00928 {
00929   THD_dataxes * daxes;                  /* dataset axes */
00930   char xxorient, yyorient, zzorient;    /* dataset axes orientations */ 
00931   int nxyz;                             /* dataset dimension in voxels */
00932   int ixyz;                             /* voxel indices */
00933   int cx, cy, cz;                       /* location of center of mass */
00934 
00935 
00936   /*----- Initialize local variables -----*/
00937   daxes = anat->daxes;
00938   xxorient = ORIENT_typestr[daxes->xxorient][0];
00939   yyorient = ORIENT_typestr[daxes->yyorient][0];
00940   zzorient = ORIENT_typestr[daxes->zzorient][0];
00941   nxyz = DSET_NX(anat) * DSET_NY(anat) * DSET_NZ(anat);
00942 
00943 
00944   /*----- Calculate center of mass of brain image -----*/
00945   center_of_mass (&cx, &cy, &cz);
00946 
00947 
00948   /*----- Initialize mask for non-brain voxels -----*/
00949   for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00950     cv[ixyz] = 1;
00951 
00952 
00953   /*----- Segment the axial slices -----*/
00954        if ((xxorient == 'S') || (xxorient == 'I'))  segment_x_slices (cv,cx,1);
00955   else if ((yyorient == 'S') || (yyorient == 'I'))  segment_y_slices (cv,cy,1);
00956   else if ((zzorient == 'S') || (zzorient == 'I'))  segment_z_slices (cv,cz,1);
00957   else SI_error ("Unable to determine dataset orientation");
00958  
00959 
00960   /*----- Segment the sagittal slices -----*/
00961        if ((xxorient == 'L') || (xxorient == 'R'))  segment_x_slices (cv,cx,0);
00962   else if ((yyorient == 'L') || (yyorient == 'R'))  segment_y_slices (cv,cy,0);
00963   else if ((zzorient == 'L') || (zzorient == 'R'))  segment_z_slices (cv,cz,0);
00964   else SI_error ("Unable to determine dataset orientation");
00965   
00966  
00967   /*----- Segment the coronal slices -----*/       
00968        if ((xxorient == 'A') || (xxorient == 'P'))  segment_x_slices (cv,cx,0);
00969   else if ((yyorient == 'A') || (yyorient == 'P'))  segment_y_slices (cv,cy,0);
00970   else if ((zzorient == 'A') || (zzorient == 'P'))  segment_z_slices (cv,cz,0);
00971   else SI_error ("Unable to determine dataset orientation");
00972 
00973 
00974   /*----- Create envelope for segmented image -----*/
00975   if (! nosmooth)
00976     segment_envelope (cv, cx, cy, cz);
00977  
00978 
00979   return;
00980 }
00981 
00982 
00983 /*---------------------------------------------------------------------------*/
00984 /*
00985   Perform voxel connectivity tests.
00986 */
00987 
00988 int connectivity_tests (short * cv)
00989 
00990 {
00991   int nx, ny, nz, nxy, nxyz;        /* voxel counters */
00992   int ix, jy, kz, ixyz;             /* voxel indices */
00993   byte * dv = NULL;                 /* volume for intermediate data */
00994 
00995 
00996   /*----- Progress report -----*/
00997   if (! quiet)  printf ("\nPerforming voxel connectivity tests \n");
00998 
00999 
01000   /*----- Initialize local variables -----*/
01001   nx = DSET_NX(anat);   ny = DSET_NY(anat);   nz = DSET_NZ(anat);
01002   nxy = nx*ny;   nxyz = nxy*nz;
01003 
01004 
01005   /*----- Initialize voxel classification indicators -----*/
01006   dv = (byte *) malloc (nxyz * sizeof(byte));
01007   MTEST (dv);   if (dv == NULL)  return (0);
01008   for (ixyz = 0;  ixyz < nxyz;  ixyz++)  
01009     dv[ixyz] = 0;
01010 
01011 
01012   /*----- Determine which voxels will leave the target structure -----*/
01013   if (max_conn_int > -1)
01014     {
01015       for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01016         {
01017           if (! cv[ixyz])
01018             {
01019               IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);
01020               if (ix > 1)     dv[THREE_TO_IJK(ix-1,jy,kz,nx,nxy)]++;
01021               if (ix < nx-1)  dv[THREE_TO_IJK(ix+1,jy,kz,nx,nxy)]++;
01022               if (jy > 1)     dv[THREE_TO_IJK(ix,jy-1,kz,nx,nxy)]++;
01023               if (jy < ny-1)  dv[THREE_TO_IJK(ix,jy+1,kz,nx,nxy)]++;
01024               if (kz > 1)     dv[THREE_TO_IJK(ix,jy,kz-1,nx,nxy)]++;
01025               if (kz < nz-1)  dv[THREE_TO_IJK(ix,jy,kz+1,nx,nxy)]++;
01026             }
01027         }
01028 
01029       for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01030         {
01031           if (dv[ixyz] <= max_conn_int)  cv[ixyz] = 1;
01032           dv[ixyz] = 0;
01033         }
01034     }
01035 
01036 
01037   /*----- Determine which voxels will enter the target structure -----*/
01038   if (min_conn_int < 7)
01039     {
01040       for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01041         {
01042           if (! cv[ixyz])
01043             {
01044               IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);
01045               if (ix > 1)     dv[THREE_TO_IJK(ix-1,jy,kz,nx,nxy)]++;
01046               if (ix < nx-1)  dv[THREE_TO_IJK(ix+1,jy,kz,nx,nxy)]++;
01047               if (jy > 1)     dv[THREE_TO_IJK(ix,jy-1,kz,nx,nxy)]++;
01048               if (jy < ny-1)  dv[THREE_TO_IJK(ix,jy+1,kz,nx,nxy)]++;
01049               if (kz > 1)     dv[THREE_TO_IJK(ix,jy,kz-1,nx,nxy)]++;          
01050               if (kz < nz-1)  dv[THREE_TO_IJK(ix,jy,kz+1,nx,nxy)]++;
01051             }
01052         }
01053       
01054       for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01055         if (dv[ixyz] >=  min_conn_int)  cv[ixyz] = 0;
01056     }
01057 
01058 
01059   /*----- Deallocate memory -----*/
01060   free (dv);   dv = NULL;
01061 
01062   
01063   return (1);
01064 }
01065 
01066 
01067 /*---------------------------------------------------------------------------*/