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 #include "mrilib.h"
00008 
00009 void find_base_value( int nxyz , short * sfim , int * base , int * peak ) ;
00010 void remove_isolated_stuff( THD_3dim_dataset * qset ) ;
00011 void xyz_to_ijk( THD_3dim_dataset * ds , float x , float y , float z ,
00012                                          int * i , int * j , int * k  ) ;
00013 void DRAW_2dfiller( int nx , int ny , int ix , int jy , byte * ar ) ;
00014 
00015 /*---------------------------------------------------------------------------*/
00016 
00017 #define DD(i,j,k) dar[(i)+(j)*nx+(k)*nxy]
00018 #define QQ(i,j,k) qar[(i)+(j)*nx+(k)*nxy]
00019 #define BB(i,j)   bar[(i)+(j)*nx]
00020 
00021 static int verbose = 0 ;
00022 
00023 int main( int argc , char * argv[] )
00024 {
00025    THD_3dim_dataset * dset , * qset ;
00026    short * dar , * qar ;
00027 
00028    char prefix[128] = "strip" ;
00029    int base = 0 , hyval = 0 , hyper = 0 ;
00030    double metric_thresh = 1.05 , area_thresh = 9999.0 ;
00031 
00032    int iarg , nx,ny,nz,nxy,nxyz , dx,dy,dz , ii,jj,kk,  ic,jc,kc , ktop , cc,pp ;
00033    byte * bar ;
00034    double qsum , sum , rr , metric,area ;
00035    int   nsum , didit ;
00036 
00037    /* help? */
00038 
00039    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00040       printf(
00041              "Usage: 3dstrip [options] dataset\n"
00042              "Attempts to strip the scalp tissue from the input anatomical dataset.\n"
00043              "This dataset must be in AC-PC or Talairach coordinates, must be stored\n"
00044              "as shorts, and must have cubical voxels.\n"
00045              "\n"
00046              "Options:\n"
00047              "  -base   bb = Set all voxels at or below value 'bb' to zero.\n"
00048              "                 [default = computed from the data histogram]\n"
00049              "  -metric mm = Set the metric threshold to 'mm'; the metric is\n"
00050              "                 used to determine the annularity of the removed\n"
00051              "                 tissue in each axial slice.\n"
00052              "                 [default = 1.05 (unitless); must be > 1.0]\n"
00053              "  -area   aa = Set the area threshold to 'aa' square millimeters;\n"
00054              "                 if the amount of removed tissue would be larger \n"
00055              "                 than 'aa', or the metric is larger than 'mm',\n"
00056              "                 then nothing will be removed in or below the slice.\n"
00057              "                 [default = 9999]\n"
00058              "  -hyper     = Flags the elimination of hyperintensity voxels after\n"
00059              "                 the stripping procedure.\n"
00060              "                 [default = do not remove these voxels]\n"
00061              "  -hyclip    = Flags the setting of hyperintensity voxels back to\n"
00062              "                 the hyperintensity threshold (rather than to zero).\n"
00063              "  -hyval  hh = Set the hyperintensity threshold to 'hh'; all voxels\n"
00064              "                 with intensity 'hh' or above will be set to zero\n"
00065              "                 (for -hyper) or set to this value (for -hyclip).\n"
00066              "                 [default = computed from the data histogram]\n"
00067              "  -prefix pp = Use 'pp' for the prefix of the output dataset.\n"
00068              "                 [default = 'strip']\n"
00069              "  -verbose   = Print out progress reports, including the area\n"
00070              "                 and metric for each axial slice\n"
00071              "\n"
00072              "The metric is average(r**2)/average(r)**2 for all candidate removal\n"
00073              "voxels, where r is measured from the center of the slice.  This ratio\n"
00074              "is 1 for a circle, and will be nearly 1 for a thin annular region.\n"
00075              "The default threshold of 1.05 is chosen from experience.  The default\n"
00076              "area threshold is also chosen from experience.  The goal of the thresholds\n"
00077              "is to avoid removing brain tissue.\n"
00078            ) ;
00079       exit(0) ;
00080    }
00081 
00082    /* arguments */
00083 
00084    iarg = 1 ;
00085    while( iarg < argc && argv[iarg][0] == '-' ){
00086 
00087       if( strncmp(argv[iarg],"-verbose",5) == 0 ){
00088          verbose = 1 ;
00089          iarg++ ; continue ;
00090       }
00091 
00092       if( strcmp(argv[iarg],"-metric") == 0 ){
00093          metric_thresh = strtod( argv[++iarg] , NULL ) ;
00094          if( metric_thresh <= 1.0 ){fprintf(stderr,"** Illegal value after -metric\n");exit(1);}
00095          iarg++ ; continue ;
00096       }
00097 
00098       if( strcmp(argv[iarg],"-area") == 0 ){
00099          area_thresh = strtod( argv[++iarg] , NULL ) ;
00100          if( area_thresh <= 1.0 ){fprintf(stderr,"** Illegal value after -area\n");exit(1);}
00101          iarg++ ; continue ;
00102       }
00103 
00104       if( strcmp(argv[iarg],"-base") == 0 ){
00105          base = strtol( argv[++iarg] , NULL , 10 ) ;
00106          if( base <= 0 ){fprintf(stderr,"** Illegal value after -base\n");exit(1);}
00107          iarg++ ; continue ;
00108       }
00109 
00110       if( strcmp(argv[iarg],"-prefix") == 0 ){
00111          strcpy( prefix , argv[++iarg] ) ;
00112          iarg++ ; continue ;
00113       }
00114 
00115       if( strcmp(argv[iarg],"-hyper") == 0 ){
00116          hyper = 1 ;
00117          iarg++ ; continue ;
00118       }
00119 
00120       if( strcmp(argv[iarg],"-hyclip") == 0 ){
00121          hyper = 2 ;
00122          iarg++ ; continue ;
00123       }
00124 
00125       if( strcmp(argv[iarg],"-hyval") == 0 ){
00126          hyval = strtol( argv[++iarg] , NULL , 10 ) ;
00127          if( hyval <= 0 ){fprintf(stderr,"** Illegal value after -hyval\n");exit(1);}
00128          iarg++ ; continue ;
00129       }
00130 
00131       fprintf(stderr,"** Illegal option: %s\n",argv[iarg]) ; exit(1) ;
00132    }
00133 
00134    /* open dataset */
00135 
00136    if( iarg >= argc ){fprintf(stderr,"** No dataset name on command line?\n");exit(1);}
00137 
00138    dset = THD_open_dataset( argv[iarg] ) ;
00139    if( dset == NULL ){fprintf(stderr,"** Can't open dataset %s\n",argv[iarg]);exit(1);}
00140 
00141    dx = DSET_DX(dset) ; dy = DSET_DY(dset) ; dz = DSET_DZ(dset) ;
00142    nx = DSET_NX(dset) ; ny = DSET_NY(dset) ; nz = DSET_NZ(dset) ;
00143 
00144    nxy = nx*ny ; nxyz = nx*ny*nz ;
00145 
00146    if( dx <= 0.0 || dx != dy || dx != dz ){
00147       fprintf(stderr,"** Dataset voxel shape is illegal!\n") ; exit(1) ;
00148    }
00149 
00150    if( dset->view_type == VIEW_ORIGINAL_TYPE ){
00151       fprintf(stderr,"** Dataset is in the +orig view, which is illegal!\n") ; exit(1) ;
00152    }
00153 
00154    DSET_mallocize(dset) ; DSET_load(dset) ;
00155    if( !DSET_LOADED(dset) ){fprintf(stderr,"** Can't load dataset %s\n",argv[iarg]);exit(1);}
00156 
00157    dar = DSET_ARRAY(dset,0) ;
00158 
00159    /* find base value, if needed */
00160 
00161    if( base <= 0 || (hyval <= 0 && hyper) ){
00162       int bb , hh ;
00163       find_base_value( nxyz , dar , &bb , &hh ) ;
00164       if( base <= 0 ){
00165          base = bb ;
00166          if( verbose ) printf("-- Base value set to %d\n",base) ;
00167       }
00168       if( hyval <= 0 && hyper && hh > base+1 ){
00169          hyval = hh ;
00170          if( verbose ) printf("-- Hyperintensity value set to %d\n",hyval) ;
00171       }
00172    }
00173 
00174    /* make new dataset */
00175 
00176    qset = EDIT_empty_copy( dset ) ;
00177    EDIT_dset_items( qset,
00178                        ADN_prefix    , prefix ,
00179                        ADN_nvals     , 1 ,
00180                        ADN_ntt       , 0 ,
00181                     ADN_none ) ;
00182 
00183    /* copy dar to qar, omiting values <= base */
00184 
00185    qar = (short *) malloc( sizeof(short) * nxyz ) ; /* another array */
00186    if( qar == NULL ){fprintf(stderr,"** Can't malloc workspace!\n");exit(1);}
00187 
00188    EDIT_substitute_brick( qset , 0 , MRI_short , qar ) ;
00189 
00190    for( ii=pp=0 ; ii < nxyz ; ii++ ){
00191       if( dar[ii] > base ){
00192          qar[ii] = dar[ii] ;
00193       } else {
00194          qar[ii] = 0 ; pp++ ;
00195       }
00196    }
00197    if( verbose ) printf("-- Blasted %d voxels at or below base\n",pp) ;
00198 
00199    /* set edges to zero */
00200 
00201    for( jj=0 ; jj < ny ; jj++ )
00202       for( ii=0 ; ii < nx ; ii++ ){ QQ(ii,jj,0) = QQ(ii,jj,nz-1) = 0 ; }
00203 
00204    for( kk=0 ; kk < nz ; kk++ )
00205       for( jj=0 ; jj < ny ; jj++ ){ QQ(0,jj,kk) = QQ(nx-1,jj,kk) = 0 ; }
00206 
00207    for( kk=0 ; kk < nz ; kk++ )
00208       for( ii=0 ; ii < nx ; ii++ ){ QQ(ii,0,kk) = QQ(ii,ny-1,kk) = 0 ; }
00209 
00210    /* if in Talairach view, remove stuff above the brain */
00211    /* [coordinate order is x=R-L, y=A-P, z=I-S]          */
00212 
00213    if( dset->view_type == VIEW_TALAIRACH_TYPE ){
00214       xyz_to_ijk( dset , 0.0 , 0.0 , 75.0 , NULL,NULL , &kc ) ;
00215       for( kk=kc ; kk < nz ; kk++ )
00216          for( jj=0 ; jj < ny ; jj++ )
00217             for( ii=0 ; ii < nx ; ii++ ) QQ(ii,jj,kk) = 0 ;
00218       ktop = kc-1 ;
00219    } else {
00220       ktop = nz-1 ;
00221    }
00222 
00223    /* now prune the dataset of the little isolated junk */
00224 
00225    remove_isolated_stuff( qset ) ;
00226 
00227    /* find center of brain */
00228 
00229    xyz_to_ijk( dset , 0.0,15.0,0.0 , &ic , &jc , NULL ) ;  /* a decent place */
00230 
00231    /* for each slice:
00232         build a mask of blocking points (zero voxels in the slice)
00233         find some outermost points
00234         flood fill from the outermost points
00235         compute a metric of "annularity" of the filled points
00236         if the slice passes the metric, blast out the filled points
00237    */
00238 
00239    bar = (byte *) malloc( sizeof(byte) * nxy ) ;
00240    didit = 0 ;
00241 
00242    for( kk=ktop ; kk >= 0 ; kk-- ){
00243 
00244       memset( bar , 0 , sizeof(byte) * nxy ) ;       /* build the mask */
00245       cc = 0 ;
00246       for( jj=0 ; jj < ny ; jj++ )
00247          for( ii=0 ; ii < nx ; ii++ )
00248             if( QQ(ii,jj,kk) == 0 ){ BB(ii,jj) = 1 ; cc++ ; }
00249 
00250       if( cc < 0.1*nxy ){
00251          if( verbose ) printf("-- Slice %3d: skipping because not enough zeros\n",kk) ;
00252          continue ;
00253       }
00254 
00255       /* scan in from right, fill from the first nonmasked point we hit */
00256 
00257       for( ii=0 ; ii < ic ; ii++ )
00258          if( BB(ii,jc) != 1 ) break ;
00259       if( ii < ic ) DRAW_2dfiller( nx,ny , ii,jc , bar ) ;
00260 
00261       /* scan in from left */
00262 
00263       for( ii=nx-1 ; ii > ic ; ii-- )
00264          if( BB(ii,jc) != 1 ) break ;
00265       if( ii > ic && BB(ii,jc) == 0 ) DRAW_2dfiller( nx,ny , ii,jc , bar ) ;
00266 
00267       /* scan in from diagonals */
00268 
00269       pp = MIN(ic,jc) ;
00270       for( ii=ic-pp,jj=jc-pp ; ii < ic && jj < jc ; ii++,jj++ )
00271          if( BB(ii,jj) != 1 ) break ;
00272       if( ii < ic && jj < jc && BB(ii,jj) == 0 ) DRAW_2dfiller( nx,ny , ii,jj , bar ) ;
00273 
00274       pp = MIN(ic,ny-1-jc) ;
00275       for( ii=ic-pp,jj=jc+pp ; ii < ic && jj > jc ; ii++,jj-- )
00276          if( BB(ii,jj) != 1 ) break ;
00277       if( ii < ic && jj > jc && BB(ii,jj) == 0 ) DRAW_2dfiller( nx,ny , ii,jj , bar ) ;
00278 
00279       pp = MIN(nx-1-ic,jc) ;
00280       for( ii=ic+pp,jj=jc-pp ; ii > ic && jj < jc ; ii--,jj++ )
00281          if( BB(ii,jj) != 1 ) break ;
00282       if( ii > ic && jj < jc && BB(ii,jj) == 0 ) DRAW_2dfiller( nx,ny , ii,jj , bar ) ;
00283 
00284       pp = MIN(nx-1-ic,ny-1-jc) ;
00285       for( ii=ic+pp,jj=jc+pp ; ii > ic && jj > jc ; ii--,jj-- )
00286          if( BB(ii,jj) != 1 ) break ;
00287       if( ii > ic && jj > jc && BB(ii,jj) == 0 ) DRAW_2dfiller( nx,ny , ii,jj , bar ) ;
00288 
00289       /* we don't scan in from top/bottom since those might be
00290          clipped off in Talairach coordinates, which would be disastrous! */
00291 
00292       /* compute metric for marked points:
00293            metric = average( r**2 ) / average(r)**2
00294            this should be nearly 1 for an annular region */
00295 
00296       qsum = sum = 0.0 ; nsum = 0 ;
00297       for( jj=0 ; jj < ny ; jj++ ){
00298          for( ii=0 ; ii < nx ; ii++ ){
00299             if( BB(ii,jj) == 2 ){
00300                rr    = (ii-ic)*(ii-ic) + (jj-jc)*(jj-jc) ;
00301                qsum += rr ;
00302                sum  += sqrt(rr) ;
00303                nsum++ ;
00304             }
00305          }
00306       }
00307 
00308       if( nsum < 10 ){
00309          if( verbose ) printf("-- Slice %3d: skipping because not enough boundary fillin\n",kk) ;
00310          continue ;
00311       }
00312 
00313       qsum   = qsum/nsum ;
00314       sum    = (sum*sum)/(nsum*nsum) ;
00315       metric = qsum / sum ;
00316       area   = nsum * dx*dx ;
00317       if( verbose ) printf("-- Slice %3d: metric =%8.4f  area =%8.0f\n",kk,metric,area) ;
00318 
00319       if( metric <= metric_thresh && area <= area_thresh ){
00320          for( jj=0 ; jj < ny ; jj++ ){
00321             for( ii=0 ; ii < nx ; ii++ ){
00322                if( BB(ii,jj) == 2 ) QQ(ii,jj,kk) = 0 ;
00323             }
00324          }
00325          didit++ ;
00326       } else if( didit > 0 ) {
00327          break ;  /* don't go below first slice that fails, after some slice worked */
00328       }
00329    }
00330 
00331    if( ! didit ){fprintf(stderr,"** No slices were trimmed -- end of run!\n");exit(1);}
00332 
00333    /* if in Talairach view, remove stuff to the sides of the brain   */
00334    /* (done after the slice stuff so the scalp annulus isn't broken) */
00335 
00336    if( dset->view_type == VIEW_TALAIRACH_TYPE ){
00337 
00338       xyz_to_ijk( dset ,  0.0, -71.0, 0.0 , NULL , &jc , NULL ) ;   /* anterior */
00339       for( kk=0 ; kk < nz ; kk++ )
00340          for( jj=0 ; jj <= jc ; jj++ )
00341             for( ii=0 ; ii < nx ; ii++ ) QQ(ii,jj,kk) = 0 ;
00342 
00343       xyz_to_ijk( dset ,  0.0, 103.0, 0.0 , NULL , &jc , NULL ) ;   /* posterior */
00344       for( kk=0 ; kk < nz ; kk++ )
00345          for( jj=jc ; jj < ny ; jj++ )
00346             for( ii=0 ; ii < nx ; ii++ ) QQ(ii,jj,kk) = 0 ;
00347 
00348       xyz_to_ijk( dset ,  -69.0, 0.0, 0.0 , &ic , NULL , NULL ) ;   /* right */
00349       for( kk=0 ; kk < nz ; kk++ )
00350          for( jj=0 ; jj < ny ; jj++ )
00351             for( ii=0 ; ii <= ic ; ii++ ) QQ(ii,jj,kk) = 0 ;
00352 
00353       xyz_to_ijk( dset ,  69.0, 0.0, 0.0 , &ic , NULL , NULL ) ;    /* left */
00354       for( kk=0 ; kk < nz ; kk++ )
00355          for( jj=0 ; jj < ny ; jj++ )
00356             for( ii=ic ; ii < nx ; ii++ ) QQ(ii,jj,kk) = 0 ;
00357 
00358       xyz_to_ijk( dset , 0.0, 0.0, -43.0 , NULL , &jc , &kc ) ;     /* below */
00359       for( kk=0 ; kk <= kc ; kk++ )
00360          for( jj=0 ; jj <= jc ; jj++ )
00361             for( ii=0 ; ii < nx ; ii++ ) QQ(ii,jj,kk) = 0 ;
00362    }
00363 
00364    /* get rid of the hyperintensity stuff */
00365 
00366    if( hyper && hyval > base+1 ){
00367       short val ;
00368 
00369            if( hyper == 1 ) val = 0 ;
00370       else if( hyper == 2 ) val = hyval ;
00371 
00372       for( ii=pp=0 ; ii < nxyz ; ii++ ){
00373          if( qar[ii] >= hyval ){ qar[ii] = val ; pp++ ; }
00374       }
00375       if( verbose ) printf("-- Processed %d voxels at or above hyval\n",pp) ;
00376    }
00377 
00378    /* re-prune the dataset of the little isolated junk */
00379 
00380    remove_isolated_stuff( qset ) ;
00381 
00382    /* write output */
00383 
00384    DSET_write(qset) ;
00385    exit(0) ;
00386 }
00387 
00388 /*==================================================================================*/
00389 
00390 #define NPMAX 128
00391 
00392 void find_base_value( int nxyz , short * sfim , int * nbase , int * hyval )
00393 {
00394    float bper = 60.0 , bmin = 1 ;
00395 
00396    int ii , kk , nbin , sval , sum , nbot , a,b,c , npeak,ntop , nvox ;
00397    int * fbin ;
00398    int   kmin[NPMAX] , kmax[NPMAX] ;
00399    int   nmin        , nmax        ;
00400 
00401    /*-- make histogram of shorts --*/
00402 
00403    fbin = (int *) malloc( sizeof(int) * 32768 ) ;
00404    for( kk=0 ; kk < 32768 ; kk++ ) fbin[kk] = 0 ;
00405 
00406    nvox = 0 ;
00407 
00408    for( ii=0 ; ii < nxyz ; ii++ ){
00409       kk = sfim[ii] ; if( kk >= 0 ){ fbin[kk]++ ; nvox++ ; }
00410    }
00411 
00412    /*-- find largest value --*/
00413 
00414    for( kk=32767 ; kk > 0 ; kk-- ) if( fbin[kk] > 0 ) break ;
00415    if( kk == 0 ){fprintf(stderr,"** find_base_value: All voxels are zero!\n");exit(1);}
00416    nbin = kk+1 ;
00417 
00418    /*-- find bper point in cumulative distribution --*/
00419 
00420    sval = 0.01 * bper * nvox ;
00421    sum  = 0 ;
00422    for( kk=0 ; kk < nbin ; kk++ ){
00423       sum += fbin[kk] ; if( sum >= sval ) break ;
00424    }
00425    nbot = kk ; if( nbot == 0 ) nbot = 1 ; if( bmin > nbot ) nbot = bmin ;
00426    if( nbot >= nbin-9 ){
00427       fprintf(stderr,"** find_base_value: Base point on histogram too high\n");
00428       exit(1);
00429    }
00430 
00431    /*-- smooth histogram --*/
00432 
00433    b = fbin[nbot-1] ; c = fbin[nbot] ;
00434    for( kk=nbot ; kk < nbin ; kk++ ){
00435       a = b ; b = c ; c = fbin[kk+1] ; fbin[kk] = 0.25*(a+c+2*b) ;
00436    }
00437 
00438    /*-- find minima and maxima above bper point --*/
00439 
00440    nmin = nmax = 0 ;
00441    for( kk=nbot+1 ; kk < nbin ; kk++ ){
00442       if( fbin[kk] < fbin[kk-1] && fbin[kk] < fbin[kk+1] && nmin < NPMAX ){
00443          kmin[nmin++] = kk ;
00444       } else if( fbin[kk] > fbin[kk-1] && fbin[kk] > fbin[kk+1] && nmax < NPMAX ){
00445          kmax[nmax++] = kk ;
00446       }
00447    }
00448 
00449    /*-- find the largest two maxima --*/
00450 
00451    if( nmax == 0 ){
00452       fprintf(stderr,"** find_base_value: No histogram maxima above base point\n");
00453       exit(1);
00454    }
00455 
00456    if( nmax == 1 ){
00457       npeak = kmax[0] ; ntop = 0 ;
00458    } else {
00459       int f1,f2 , k1,k2 , fk , klow,kup ;
00460 
00461       k1 = 0 ; f1 = fbin[kmax[0]] ;
00462       k2 = 1 ; f2 = fbin[kmax[1]] ;
00463       if( f1 < f2 ){
00464          k1 = 1 ; f1 = fbin[kmax[1]] ;
00465          k2 = 0 ; f2 = fbin[kmax[0]] ;
00466       }
00467 
00468       for( kk=2 ; kk < nmax ; kk++ ){
00469          fk = fbin[kmax[kk]] ;
00470          if( fk > f1 ){
00471             f2 = f1 ; k2 = k1 ;
00472             f1 = fk ; k1 = kk ;
00473          } else if( fk > f2 ){
00474             f2 = fk ; k2 = kk ;
00475          }
00476       }
00477       npeak = MIN( kmax[k1] , kmax[k2] ) ;  /* smaller bin of the 2 top peaks */
00478 
00479       /* find valley between 2 peaks */
00480 
00481       ntop  = MAX( kmax[k1] , kmax[k2] ) ;
00482 
00483       fk = fbin[ntop] ; klow = ntop ;
00484       for( kk=ntop-1 ; kk >= npeak ; kk-- ){
00485          if( fbin[kk] < fk ){ fk = fbin[kk] ; klow = kk ; }
00486       }
00487       fk  = MAX( 0.10*fk , 0.05*fbin[ntop] ) ;
00488       kup = MIN( nbin-1 , ntop+3*(ntop-klow+2) ) ;
00489       for( kk=ntop+1 ; kk <= kup ; kk++ ) if( fbin[kk] < fk ) break ;
00490 
00491       ntop = kk ;
00492    }
00493 
00494    for( kk=npeak-1 ; kk > 0 ; kk-- )
00495       if( fbin[kk] < fbin[kk-1] && fbin[kk] < fbin[kk+1] ) break ;
00496 
00497    if( ntop == 0 ) ntop = npeak + (npeak-kk) ;
00498 
00499    *nbase = kk ;
00500    *hyval = ntop ;
00501 
00502    free(fbin) ; return ;
00503 }
00504 
00505 /*==================================================================================*/
00506 
00507 void remove_isolated_stuff( THD_3dim_dataset * qset )
00508 {
00509    short * qar = DSET_ARRAY(qset,0) ;
00510    short nb[27] ;
00511    int   nblast , nx,ny,nz , nxy,nxyz , ii,jj,kk,ll,cc ;
00512    float dx = DSET_DX(qset) ;
00513    EDIT_options edopt ;
00514 
00515    nx = DSET_NX(qset) ; ny = DSET_NY(qset) ; nz = DSET_NZ(qset) ;
00516    nxy = nx*ny ; nxyz = nx*ny*nz ;
00517 
00518    do{ nblast = 0 ;
00519 
00520        /* edit 3x3x3 volumes */
00521 
00522        for( kk=1 ; kk < nz-1 ; kk++ ){
00523          for( jj=1 ; jj < ny-1 ; jj++ ){
00524             for( ii=1 ; ii < nx-1 ; ii++ ){
00525                if( QQ(ii,jj,kk) > 0 ){           /* must have at least 3 nonzero   */
00526                   nb[ 0] = QQ(ii-1,jj-1,kk-1) ;  /* voxels in a 3x3x3 neighborhood */
00527                   nb[ 1] = QQ(ii  ,jj-1,kk-1) ;
00528                   nb[ 2] = QQ(ii+1,jj-1,kk-1) ;
00529                   nb[ 3] = QQ(ii-1,jj  ,kk-1) ;
00530                   nb[ 4] = QQ(ii  ,jj  ,kk-1) ;
00531                   nb[ 5] = QQ(ii+1,jj  ,kk-1) ;
00532                   nb[ 6] = QQ(ii-1,jj+1,kk-1) ;
00533                   nb[ 7] = QQ(ii  ,jj+1,kk-1) ;
00534                   nb[ 8] = QQ(ii+1,jj+1,kk-1) ;
00535                   nb[ 9] = QQ(ii-1,jj-1,kk  ) ;
00536                   nb[10] = QQ(ii  ,jj-1,kk  ) ;
00537                   nb[11] = QQ(ii+1,jj-1,kk  ) ;
00538                   nb[12] = QQ(ii-1,jj  ,kk  ) ;
00539                   nb[13] = QQ(ii  ,jj  ,kk  ) ;
00540                   nb[14] = QQ(ii+1,jj  ,kk  ) ;
00541                   nb[15] = QQ(ii-1,jj-1,kk  ) ;
00542                   nb[16] = QQ(ii  ,jj-1,kk  ) ;
00543                   nb[17] = QQ(ii+1,jj-1,kk  ) ;
00544                   nb[18] = QQ(ii-1,jj+1,kk+1) ;
00545                   nb[19] = QQ(ii  ,jj+1,kk+1) ;
00546                   nb[20] = QQ(ii+1,jj+1,kk+1) ;
00547                   nb[21] = QQ(ii-1,jj  ,kk+1) ;
00548                   nb[22] = QQ(ii  ,jj  ,kk+1) ;
00549                   nb[23] = QQ(ii+1,jj  ,kk+1) ;
00550                   nb[24] = QQ(ii-1,jj+1,kk+1) ;
00551                   nb[25] = QQ(ii  ,jj+1,kk+1) ;
00552                   nb[26] = QQ(ii+1,jj+1,kk+1) ;
00553 
00554                   for( ll=cc=0 ; ll < 27 ; ll++ ) if( nb[ll] > 0 ) cc++ ;
00555                   if( cc < 4 ){ QQ(ii,jj,kk) = 0 ; nblast++ ; }
00556                }
00557             }
00558          }
00559        }
00560 
00561        /* edit 3x3 areas in each slice orientation */
00562 
00563        for( kk=1 ; kk < nz-1 ; kk++ ){
00564          for( jj=1 ; jj < ny-1 ; jj++ ){
00565             for( ii=1 ; ii < nx-1 ; ii++ ){
00566                if( QQ(ii,jj,kk) > 0 ){         /* must have at least 2 nonzero */
00567                   nb[ 0] = QQ(ii-1,jj-1,kk) ;  /* voxels in a 3x3 neighborhood */
00568                   nb[ 1] = QQ(ii  ,jj-1,kk) ;
00569                   nb[ 2] = QQ(ii+1,jj-1,kk) ;
00570                   nb[ 3] = QQ(ii-1,jj  ,kk) ;
00571                   nb[ 4] = QQ(ii  ,jj  ,kk) ;
00572                   nb[ 5] = QQ(ii+1,jj  ,kk) ;
00573                   nb[ 6] = QQ(ii-1,jj+1,kk) ;
00574                   nb[ 7] = QQ(ii  ,jj+1,kk) ;
00575                   nb[ 8] = QQ(ii+1,jj+1,kk) ;
00576                   for( ll=cc=0 ; ll < 9 ; ll++ ) if( nb[ll] > 0 ) cc++ ;
00577                   if( cc < 2 ){ QQ(ii,jj,kk) = 0 ; nblast++ ; }
00578                }
00579             }
00580          }
00581        }
00582 
00583        for( jj=1 ; jj < ny-1 ; jj++ ){
00584           for( kk=1 ; kk < nz-1 ; kk++ ){
00585             for( ii=1 ; ii < nx-1 ; ii++ ){
00586                if( QQ(ii,jj,kk) > 0 ){         /* must have at least 2 nonzero */
00587                   nb[ 0] = QQ(ii-1,jj,kk-1) ;  /* voxels in a 3x3 neighborhood */
00588                   nb[ 1] = QQ(ii  ,jj,kk-1) ;
00589                   nb[ 2] = QQ(ii+1,jj,kk-1) ;
00590                   nb[ 3] = QQ(ii-1,jj,kk  ) ;
00591                   nb[ 4] = QQ(ii  ,jj,kk  ) ;
00592                   nb[ 5] = QQ(ii+1,jj,kk  ) ;
00593                   nb[ 6] = QQ(ii-1,jj,kk+1) ;
00594                   nb[ 7] = QQ(ii  ,jj,kk+1) ;
00595                   nb[ 8] = QQ(ii+1,jj,kk+1) ;
00596                   for( ll=cc=0 ; ll < 9 ; ll++ ) if( nb[ll] > 0 ) cc++ ;
00597                   if( cc < 2 ){ QQ(ii,jj,kk) = 0 ; nblast++ ; }
00598                }
00599             }
00600          }
00601        }
00602 
00603        for( ii=1 ; ii < nx-1 ; ii++ ){
00604          for( jj=1 ; jj < ny-1 ; jj++ ){
00605             for( kk=1 ; kk < nz-1 ; kk++ ){
00606                if( QQ(ii,jj,kk) > 0 ){         /* must have at least 2 nonzero */
00607                   nb[ 0] = QQ(ii,jj-1,kk-1) ;  /* voxels in a 3x3 neighborhood */
00608                   nb[ 1] = QQ(ii,jj  ,kk-1) ;
00609                   nb[ 2] = QQ(ii,jj+1,kk-1) ;
00610                   nb[ 3] = QQ(ii,jj-1,kk  ) ;
00611                   nb[ 4] = QQ(ii,jj  ,kk  ) ;
00612                   nb[ 5] = QQ(ii,jj+1,kk  ) ;
00613                   nb[ 6] = QQ(ii,jj-1,kk+1) ;
00614                   nb[ 7] = QQ(ii,jj  ,kk+1) ;
00615                   nb[ 8] = QQ(ii,jj+1,kk+1) ;
00616                   for( ll=cc=0 ; ll < 9 ; ll++ ) if( nb[ll] > 0 ) cc++ ;
00617                   if( cc < 2 ){ QQ(ii,jj,kk) = 0 ; nblast++ ; }
00618                }
00619             }
00620          }
00621        }
00622 
00623       if( verbose ) printf("-- Blasted %d voxels for being too isolated\n",nblast) ;
00624    } while( nblast > 0 ) ;
00625 
00626    /* cluster data, remove small clusters */
00627 
00628    INIT_EDOPT( &edopt ) ;
00629    edopt.edit_clust = ECFLAG_SAME;
00630    edopt.clust_rmm  = 1.01*dx ;
00631    edopt.clust_vmul = 25000.0 ;
00632    EDIT_one_dataset( qset , &edopt ) ;
00633 
00634    return ;
00635 }
00636 
00637 /*---------------------------------------------------------------------------*/
00638 
00639 void xyz_to_ijk( THD_3dim_dataset * ds , float x , float y , float z ,
00640                                          int * i , int * j , int * k  )
00641 {
00642    THD_fvec3 fv ;
00643    THD_ivec3 iv ;
00644 
00645    LOAD_FVEC3( fv , x,y,z ) ;
00646    fv = THD_dicomm_to_3dmm( ds , fv ) ;
00647    iv = THD_3dmm_to_3dind ( ds , fv ) ;
00648    if( i != NULL ) *i = iv.ijk[0] ;
00649    if( j != NULL ) *j = iv.ijk[1] ;
00650    if( k != NULL ) *k = iv.ijk[2] ;
00651    return ;
00652 }
00653 
00654 /*---------------------------------------------------------------------------
00655    Flood filling a byte array:
00656      nx = 1st dimension
00657      ny = 2nd dimension
00658      ix = start point
00659      jy = end point
00660      ar = array, with 0's everwhere except 1's as barriers to flooding
00661 
00662    All filled points (starting with ix,jy) will get the value 2.
00663 -----------------------------------------------------------------------------*/
00664 
00665 void DRAW_2dfiller( int nx , int ny , int ix , int jy , byte * ar )
00666 {
00667    int ii,jj , ip,jp , num ;
00668 
00669 #define AR(i,j) ar[(i)+(j)*nx]
00670 
00671    /* fill out in cross from 1st point */
00672 
00673    ip = ix ; jp = jy ; AR(ip,jp) = 2 ;
00674 
00675    for( ii=ip+1; ii < nx && AR(ii,jp) == 0; ii++ ) AR(ii,jp) = 2;
00676    for( ii=ip-1; ii >= 0 && AR(ii,jp) == 0; ii-- ) AR(ii,jp) = 2;
00677    for( jj=jp+1; jj < ny && AR(ip,jj) == 0; jj++ ) AR(ip,jj) = 2;
00678    for( jj=jp-1; jj >= 0 && AR(ip,jj) == 0; jj-- ) AR(ip,jj) = 2;
00679 
00680    /* brute force repetition of the cross technique */
00681 
00682    do {
00683       num = 0 ;
00684       for( jp=0 ; jp < ny ; jp++ ){
00685          for( ip=0 ; ip < nx ; ip++ ){
00686             if( AR(ip,jp) == 2 ){
00687                for( ii=ip+1; ii < nx && AR(ii,jp) == 0; ii++ ){ AR(ii,jp) = 2; num++; }
00688                for( ii=ip-1; ii >= 0 && AR(ii,jp) == 0; ii-- ){ AR(ii,jp) = 2; num++; }
00689                for( jj=jp+1; jj < ny && AR(ip,jj) == 0; jj++ ){ AR(ip,jj) = 2; num++; }
00690                for( jj=jp-1; jj >= 0 && AR(ip,jj) == 0; jj-- ){ AR(ip,jj) = 2; num++; }
00691             }
00692          }
00693       }
00694    } while( num > 0 ) ;
00695 
00696    return ;
00697 }