Skip to content

AFNI/NIfTI Server

Sections
Personal tools
You are here: Home » AFNI » Documentation
 

Doxygen Source Code Documentation

Main Page   Alphabetical List   Data Structures   File List   Data Fields   Globals   Search  

thd_automask.c

Go to the documentation of this file.
00001 #include "mrilib.h"
00002 
00003 #define USE_DILATE
00004 #define USE_FILLIN
00005 
00006 #undef DEBUG
00007 
00008 static int verb = 0 ;                            /* 28 Oct 2003 */
00009 void THD_automask_verbose( int v ){ verb = v ; }
00010 
00011 static int exterior_clip = 0 ;
00012 void THD_automask_extclip( int e ){ exterior_clip = e ; }
00013 
00014 static int dall = 0 ;
00015 
00016 /*---------------------------------------------------------------------*/
00017 
00018 static int mask_count( int nvox , byte *mmm )
00019 {
00020    int ii , nn ;
00021    for( nn=ii=0 ; ii < nvox ; ii++ ) nn += (mmm[ii] != 0) ;
00022    return nn ;
00023 }
00024 
00025 /*---------------------------------------------------------------------*/
00026 /*! Make a byte mask for a 3D+time dataset -- 13 Aug 2001 - RWCox.
00027      - compare to thd_makemask.c
00028      - 05 Mar 2003: modified to put most code into mri_automask_image().
00029 -----------------------------------------------------------------------*/
00030 
00031 byte * THD_automask( THD_3dim_dataset *dset )
00032 {
00033    MRI_IMAGE *medim ;
00034    byte *mmm ;
00035 
00036 ENTRY("THD_automask") ;
00037 
00038    /* median at each voxel */
00039 
00040    medim = THD_median_brick(dset) ; if( medim == NULL ) RETURN(NULL);
00041 
00042    mmm = mri_automask_image( medim ) ;
00043 
00044    mri_free(medim) ; RETURN(mmm) ;
00045 }
00046 
00047 /*---------------------------------------------------------------------*/
00048 /*! Make a byte mask from the average of an array of 3D images.
00049     We assume that they all have the same (nx,ny,nz) dimensions.
00050 -----------------------------------------------------------------------*/
00051 
00052 byte * mri_automask_imarr( MRI_IMARR *imar )  /* 18 Nov 2004 */
00053 {
00054    MRI_IMAGE *avim , *tim , *qim ;
00055    byte *mmm ;
00056    int ii , jj , nvox,nim ;
00057    float fac , *avar , *qar ;
00058 
00059 ENTRY("mri_automask_imarr") ;
00060 
00061    if( imar == NULL || IMARR_COUNT(imar) < 1 ) RETURN(NULL) ;
00062 
00063    nim = IMARR_COUNT(imar) ;
00064    if( nim == 1 ){
00065      mmm = mri_automask_image( IMARR_SUBIMAGE(imar,0) ) ;
00066      RETURN(mmm) ;
00067    }
00068 
00069    avim = mri_new_conforming( IMARR_SUBIMAGE(imar,0) , MRI_float ) ;
00070    avar = MRI_FLOAT_PTR(avim) ;
00071    nvox = avim->nvox ;
00072    for( jj=0 ; jj < nim ; jj++ ){
00073      tim = IMARR_SUBIMAGE(imar,jj) ;
00074      if( tim->kind != MRI_float ) qim = mri_to_float(tim) ;
00075      else                         qim = tim ;
00076      qar = MRI_FLOAT_PTR(qim) ;
00077      for( ii=0 ; ii < nvox ; ii++ ) avar[ii] += qar[ii] ;
00078      if( qim != tim ) mri_free(qim) ;
00079    }
00080    fac = 1.0f / (float)nim ;
00081    for( ii=0 ; ii < nvox ; ii++ ) avar[ii] *= fac ;
00082    mmm = mri_automask_image( avim ) ;
00083    mri_free(avim) ;
00084    RETURN(mmm) ;
00085 }
00086 
00087 /*---------------------------------------------------------------------*/
00088 /*! Make a byte mask from an image (3D).  Adapted from THD_automask()
00089     to make it possible to do this on an image directly.
00090 -----------------------------------------------------------------------*/
00091 
00092 byte * mri_automask_image( MRI_IMAGE *im )
00093 {
00094    float clip_val , *mar ;
00095    byte *mmm = NULL ;
00096    int nvox , ii,jj , nmm , nx,ny,nz ;
00097    MRI_IMAGE *medim ;
00098 
00099 ENTRY("mri_automask_image") ;
00100 
00101    if( im == NULL ) return NULL ;
00102 
00103    if( im->kind != MRI_float ) medim = mri_to_float(im) ;
00104    else                        medim = im ;
00105 
00106    /* find clip value to excise small stuff */
00107 
00108    clip_val = THD_cliplevel(medim,0.5) ;
00109 
00110    if( verb ) ININFO_message("Clip level = %f\n",clip_val) ;
00111 
00112    /* create mask of values above clip value */
00113 
00114    nvox = medim->nvox ;
00115    mar  = MRI_FLOAT_PTR(medim) ;
00116    mmm  = (byte *) calloc( sizeof(byte), nvox ) ;
00117    for( nmm=ii=0 ; ii < nvox ; ii++ )
00118      if( mar[ii] >= clip_val ){ mmm[ii] = 1; nmm++; }
00119 
00120    if( AFNI_yesenv("TOPCLIP") ){
00121      float tclip = 3.1*clip_val ;
00122      if( verb ) ININFO_message("Top clip = %f\n",tclip) ;
00123      for( ii=0 ; ii < nvox ; ii++ )
00124        if( mar[ii] > tclip ) mmm[ii] = 0 ;
00125      for( nmm=ii=0 ; ii < nvox ; ii++ ) if( mmm[ii] ) nmm++ ;
00126    }
00127 
00128    if( verb ) ININFO_message("Number voxels above clip level = %d\n",nmm) ;
00129    if( im != medim && (!exterior_clip || nmm==0) ){ mri_free(medim); medim=NULL; }
00130    if( nmm == 0 ) RETURN(mmm) ;  /* should not happen */
00131 
00132    /*-- 10 Apr 2002: only keep the largest connected component --*/
00133 
00134    nx = im->nx ; ny = im->ny ; nz = im->nz ;
00135    dall = (nx*ny*nz)/128 ;  /* allocation delta for clustering */
00136 
00137    if( verb ) ININFO_message("Clustering voxels above clip level ...\n") ;
00138    THD_mask_clust( nx,ny,nz, mmm ) ;
00139 
00140    /* 18 Apr 2002: now erode the resulting volume
00141                    (to break off any thinly attached pieces) */
00142 
00143    THD_mask_erode( nx,ny,nz, mmm ) ;
00144 
00145    /* now recluster it, and again keep only the largest survivor */
00146 
00147    THD_mask_clust( nx,ny,nz, mmm ) ;
00148 
00149 #ifdef USE_FILLIN
00150    /* 19 Apr 2002: fill in small holes */
00151 
00152    jj = ii = THD_mask_fillin_once( nx,ny,nz , mmm , 1 ) ;
00153    if( ii > 0 ){
00154      jj += ii = THD_mask_fillin_once( nx,ny,nz , mmm , 1 ) ;
00155      if( ii > 0 ){
00156        jj += ii = THD_mask_fillin_once( nx,ny,nz , mmm , 1 ) ;
00157      }
00158    }
00159    if( jj > 0 && verb )
00160     ININFO_message("Filled   %d voxels in small holes; now have %d voxels\n",
00161             jj , mask_count(nvox,mmm) ) ;
00162 
00163    if( AFNI_yesenv("PEEL") ){
00164      jj = THD_peel_mask( nx,ny,nz , mmm , 7 ) ;
00165      if( jj > 0 ){
00166        ININFO_message("Peeled %d voxels from surface\n",jj) ;
00167        THD_mask_erode( nx,ny,nz, mmm ) ;
00168        THD_mask_clust( nx,ny,nz, mmm ) ;
00169      }
00170    }
00171 
00172    nmm = 1 ;
00173    jj  = rint(0.016*nx) ; nmm = MAX(nmm,jj) ;
00174    jj  = rint(0.016*ny) ; nmm = MAX(nmm,jj) ;
00175    jj  = rint(0.016*nz) ; nmm = MAX(nmm,jj) ;
00176 
00177    if( nmm > 1 || jj > 0 ){
00178      for( jj=0,ii=2 ; ii < nmm ; ii++ )
00179        jj += THD_mask_fillin_once( nx,ny,nz , mmm , ii ) ;
00180      jj += THD_mask_fillin_completely( nx,ny,nz, mmm , nmm ) ;
00181      if( jj > 0 && verb )
00182       ININFO_message("Filled   %d voxels in large holes; now have %d voxels\n",
00183               jj , mask_count(nvox,mmm) ) ;
00184    }
00185 
00186    THD_mask_erode( nx,ny,nz, mmm ) ;
00187    THD_mask_clust( nx,ny,nz, mmm ) ;
00188 #endif
00189 
00190    /* 28 May 2002:
00191       invert the mask, then find the largest cluster of 1's;
00192       this will be the outside of the brain;
00193       put this back into the mask, then invert again;
00194       the effect will be to fill any holes left inside the brain */
00195 
00196    for( ii=0 ; ii < nvox ; ii++ ) mmm[ii] = !mmm[ii] ;
00197 
00198    if( verb ) ININFO_message("Clustering non-brain voxels ...\n") ;
00199    THD_mask_clust( nx,ny,nz, mmm ) ;
00200 
00201    /* mask is now 1 for non-brain voxels;
00202       if we want to clip off voxels neighboring the non-brain
00203       mask AND whose values are below clip_val, do so now     */
00204 
00205    if( exterior_clip ){
00206      float tclip ;
00207      tclip = AFNI_yesenv("TOPCLIP") ? 3.1*clip_val : 9999.9*clip_val ;;
00208      jj = THD_mask_clip_neighbors( nx,ny,nz , mmm , clip_val,tclip,mar ) ;
00209      if( im != medim ) mri_free(medim) ;
00210      if( jj > 0 && verb )
00211        ININFO_message("Removed  %d exterior voxels below clip level\n",jj);
00212    } else {
00213      jj = 0 ;
00214    }
00215 
00216    /* and re-invert mask to get brain voxels */
00217 
00218    for( ii=0 ; ii < nvox ; ii++ ) mmm[ii] = !mmm[ii] ;
00219    if( verb ) ININFO_message("Mask now has %d voxels\n",mask_count(nvox,mmm)) ;
00220 
00221    if( exterior_clip && jj > 0 ){
00222      THD_mask_erode( nx,ny,nz, mmm ) ;
00223      THD_mask_clust( nx,ny,nz, mmm ) ;
00224    }
00225 
00226    RETURN(mmm) ;
00227 }
00228 
00229 /*---------------------------------------------------------------------*/
00230 /*! Find voxels not in the mask, but neighboring it, that are also
00231     below the clip threshold in mar[].  Add these to the mask.
00232     Repeat until done.  Return value is number of voxels added.
00233 -----------------------------------------------------------------------*/
00234 
00235 int THD_mask_clip_neighbors( int nx, int ny, int nz ,
00236                             byte *mmm, float clip_val, float tclip, float *mar )
00237 {
00238    int ii,jj,kk , ntot=0,nnew , jm,jp,j3 , km,kp,k3 , im,ip,i3 , nxy=nx*ny ;
00239 
00240    if( mmm == NULL || mar == NULL ) return 0 ;
00241 
00242    do{
00243     nnew = 0 ;
00244     for( kk=1 ; kk < nz-1 ; kk++ ){
00245      k3 = kk*nxy ;
00246      for( jj=1 ; jj < ny-1 ; jj++ ){
00247       j3 = k3 + jj*nx ;
00248       for( ii=1 ; ii < nx-1 ; ii++ ){
00249        i3 = ii+j3 ;
00250        if( mmm[i3] ||                                             /* in mask */
00251            (mar[i3] >= clip_val && mar[i3] <= tclip) ) continue ; /* or is OK */
00252 
00253        /* If here, voxel IS NOT in mask, and IS below threshold.
00254           If any neighbors are also in mask, then add it to mask. */
00255 
00256        if( mmm[i3-1]   || mmm[i3+1]   ||
00257            mmm[i3-nx]  || mmm[i3+nx]  ||
00258            mmm[i3-nxy] || mmm[i3+nxy]   ){ mmm[i3] = 1; nnew++; }
00259     }}}
00260     ntot += nnew ;
00261    } while( nnew > 0 ) ;
00262 
00263    return ntot ;
00264 }
00265 
00266 /*---------------------------------------------------------------------*/
00267 /*! Fill in a byte mask.  Filling is done by looking to each side
00268     (plus/minus) of a non-filled voxel, and seeing if there is a
00269     filled voxel on both sides.  This looking is done parallel to
00270     the x-, y-, and z-axes, out to distance nside voxels.
00271     - nx,ny,nz = dimensions of mask
00272     - mmm      = mask itself (will be altered)
00273     - nside    = width of fill in look to each side
00274     - Return value is number of filled in voxels
00275 -----------------------------------------------------------------------*/
00276 
00277 int THD_mask_fillin_once( int nx, int ny, int nz, byte *mmm, int nside )
00278 {
00279    int ii,jj,kk , nsx,nsy,nsz , nxy,nxyz , iv,jv,kv,ll , nfill ;
00280    byte *nnn ;
00281    int nx2,nx3,nx4 , nxy2,nxy3,nxy4 ;
00282 
00283 ENTRY("THD_mask_fillin_once") ;
00284 
00285    if( mmm == NULL || nside <= 0 ) RETURN(0) ;
00286 
00287    nsx = (nx-1)/2 ; if( nsx > nside ) nsx = nside ;
00288    nsy = (ny-1)/2 ; if( nsy > nside ) nsy = nside ;
00289    nsz = (nz-1)/2 ; if( nsz > nside ) nsz = nside ;
00290 
00291    if( nsx == 0 && nsy == 0 && nsz == 0 ) RETURN(0) ;
00292 
00293 #ifdef DEBUG
00294    fprintf(stderr,"THD_mask_fillin_once: nsx=%d nsy=%d nsz=%d\n",nsx,nsy,nsz);
00295 #endif
00296 
00297    nxy = nx*ny ; nxyz = nxy*nz ; nfill = 0 ;
00298 
00299    nx2  = 2*nx  ; nx3  = 3*nx  ; nx4  = 4*nx  ;
00300    nxy2 = 2*nxy ; nxy3 = 3*nxy ; nxy4 = 4*nxy ;
00301 
00302    nnn = AFMALL(byte, nxyz) ;  /* stores filled in values */
00303 
00304    /* loop over voxels */
00305 
00306 #define FILLVOX                                     \
00307  do{ nnn[iv] = 1; nfill++; goto NextVox; } while(0)
00308 
00309    for( kk=nsz ; kk < nz-nsz ; kk++ ){
00310      kv = kk*nxy ;
00311      for( jj=nsy ; jj < ny-nsy ; jj++ ){
00312        jv = jj*nx + kv ;
00313        for( ii=nsx ; ii < nx-nsx ; ii++ ){
00314          iv = ii+jv ;
00315          if( mmm[iv] ) continue ;     /* already filled */
00316 
00317          /* check in +x direction, then -x if +x hits */
00318 
00319          switch( nsx ){
00320            case 1:
00321              if( mmm[iv+1] && mmm[iv-1] ) FILLVOX;
00322            break ;
00323 
00324            case 2:
00325              if( (mmm[iv+1]||mmm[iv+2]) &&
00326                  (mmm[iv-1]||mmm[iv-2])   ) FILLVOX;
00327            break ;
00328 
00329            case 3:
00330              if( (mmm[iv+1]||mmm[iv+2]||mmm[iv+3]) &&
00331                  (mmm[iv-1]||mmm[iv-2]||mmm[iv-3])   ) FILLVOX;
00332            break ;
00333 
00334            case 4:
00335              if( (mmm[iv+1]||mmm[iv+2]||mmm[iv+3]||mmm[iv+4]) &&
00336                  (mmm[iv-1]||mmm[iv-2]||mmm[iv-3]||mmm[iv-4])   ) FILLVOX;
00337            break ;
00338 
00339            default:
00340              for( ll=1 ; ll <= nsx ; ll++ ) if( mmm[iv+ll] ) break ;
00341              if( ll <= nsx ){
00342                for( ll=1 ; ll <= nsx ; ll++ ) if( mmm[iv-ll] ) break ;
00343                if( ll <= nsx ) FILLVOX;
00344              }
00345            break ;
00346          }
00347 
00348          /* check in +y direction, then -y if +y hits */
00349 
00350          switch( nsy ){
00351            case 1:
00352              if( mmm[iv+nx] && mmm[iv-nx] ) FILLVOX;
00353            break ;
00354 
00355            case 2:
00356              if( (mmm[iv+nx]||mmm[iv+nx2]) &&
00357                  (mmm[iv-nx]||mmm[iv-nx2])   ) FILLVOX;
00358            break ;
00359 
00360            case 3:
00361              if( (mmm[iv+nx]||mmm[iv+nx2]||mmm[iv+nx3]) &&
00362                  (mmm[iv-nx]||mmm[iv-nx2]||mmm[iv-nx3])   ) FILLVOX;
00363            break ;
00364 
00365            case 4:
00366              if( (mmm[iv+nx]||mmm[iv+nx2]||mmm[iv+nx3]||mmm[iv+nx4]) &&
00367                  (mmm[iv-nx]||mmm[iv-nx2]||mmm[iv-nx3]||mmm[iv-nx4])   ) FILLVOX;
00368            break ;
00369 
00370            default:
00371              for( ll=1 ; ll <= nsy ; ll++ ) if( mmm[iv+ll*nx] ) break ;
00372              if( ll <= nsy ){
00373                for( ll=1 ; ll <= nsy ; ll++ ) if( mmm[iv-ll*nx] ) break ;
00374                if( ll <= nsy ) FILLVOX;
00375              }
00376            break ;
00377          }
00378 
00379          /* check in +z direction, then -z if +z hits */
00380 
00381          switch( nsz ){
00382            case 1:
00383              if( mmm[iv+nxy] && mmm[iv-nxy] ) FILLVOX;
00384            break ;
00385 
00386            case 2:
00387              if( (mmm[iv+nxy]||mmm[iv+nxy2]) &&
00388                  (mmm[iv-nxy]||mmm[iv-nxy2])   ) FILLVOX;
00389            break ;
00390 
00391            case 3:
00392              if( (mmm[iv+nxy]||mmm[iv+nxy2]||mmm[iv+nxy3]) &&
00393                  (mmm[iv-nxy]||mmm[iv-nxy2]||mmm[iv-nxy3])   ) FILLVOX;
00394            break ;
00395 
00396            case 4:
00397              if( (mmm[iv+nxy]||mmm[iv+nxy2]||mmm[iv+nxy3]||mmm[iv+nxy4]) &&
00398                  (mmm[iv-nxy]||mmm[iv-nxy2]||mmm[iv-nxy3]||mmm[iv-nxy4])   ) FILLVOX;
00399            break ;
00400 
00401            default:
00402              for( ll=1 ; ll <= nsz ; ll++ ) if( mmm[iv+ll*nxy] ) break ;
00403              if( ll <= nsz ){
00404                for( ll=1 ; ll <= nsz ; ll++ ) if( mmm[iv-ll*nxy] ) break ;
00405                if( ll <= nsz ) FILLVOX;
00406              }
00407            break ;
00408          }
00409 
00410          NextVox: ; /* end of loop over ii */
00411    } } }
00412 
00413    /* copy fills back into mmm */
00414 
00415    if( nfill > 0 ){
00416      for( iv=0 ; iv < nxyz ; iv++ ) if( nnn[iv] ) mmm[iv] = 1 ;
00417    }
00418 
00419 #ifdef DEBUG
00420    fprintf(stderr,"THD_mask_fillin_once: nfill=%d\n",nfill) ;
00421 #endif
00422 
00423    free(nnn) ; RETURN(nfill) ;
00424 }
00425 
00426 /*----------------------------------------------------------------------*/
00427 /*! Fill in a byte mask, repeatedly until it doesn't fill any more.
00428     Return value is number of voxels filled.
00429 ------------------------------------------------------------------------*/
00430 
00431 int THD_mask_fillin_completely( int nx, int ny, int nz, byte *mmm, int nside )
00432 {
00433    int nfill=0 , kfill ;
00434 
00435 ENTRY("THD_mask_fillin_completely") ;
00436 
00437    do{
00438       kfill = THD_mask_fillin_once( nx,ny,nz , mmm , nside ) ;
00439       nfill += kfill ;
00440    } while( kfill > 0 ) ;
00441 
00442    RETURN(nfill) ;
00443 }
00444 
00445 /*------------------------------------------------------------------*/
00446 
00447 # define DALL 1024  /* Allocation size for cluster arrays */
00448 
00449 /*! Put (i,j,k) into the current cluster, if it is nonzero. */
00450 
00451 # define CPUT(i,j,k)                                            \
00452   do{ ijk = THREE_TO_IJK(i,j,k,nx,nxy) ;                        \
00453       if( mmm[ijk] ){                                           \
00454         if( nnow == nall ){ /* increase array lengths */        \
00455           nall += dall ;                                        \
00456           inow = (short *) realloc(inow,sizeof(short)*nall) ;   \
00457           jnow = (short *) realloc(jnow,sizeof(short)*nall) ;   \
00458           know = (short *) realloc(know,sizeof(short)*nall) ;   \
00459         }                                                       \
00460         inow[nnow] = i ; jnow[nnow] = j ; know[nnow] = k ;      \
00461         nnow++ ; mmm[ijk] = 0 ;                                 \
00462       } } while(0)
00463 
00464 /*------------------------------------------------------------------*/
00465 /*! Find the biggest cluster of nonzeros in the byte mask mmm.
00466 --------------------------------------------------------------------*/
00467 
00468 void THD_mask_clust( int nx, int ny, int nz, byte *mmm )
00469 {
00470    int ii,jj,kk, icl ,  nxy,nxyz , ijk , ijk_last , mnum ;
00471    int ip,jp,kp , im,jm,km ;
00472    int nbest ; short *ibest, *jbest , *kbest ;
00473    int nnow  ; short *inow , *jnow  , *know  ; int nall ;
00474 
00475    if( mmm == NULL ) return ;
00476 
00477    nxy = nx*ny ; nxyz = nxy * nz ;
00478 
00479    nbest = 0 ;
00480    ibest = AFMALL(short, sizeof(short)) ;
00481    jbest = AFMALL(short, sizeof(short)) ;
00482    kbest = AFMALL(short, sizeof(short)) ;
00483 
00484    /*--- scan through array, find nonzero point, build a cluster, ... ---*/
00485 
00486    ijk_last = 0 ; if( dall < DALL ) dall = DALL ;
00487    while(1) {
00488      /* find next nonzero point */
00489 
00490      for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( mmm[ijk] ) break ;
00491      if( ijk == nxyz ) break ;  /* didn't find any! */
00492 
00493      ijk_last = ijk+1 ;         /* start here next time */
00494 
00495      /* init current cluster list with this point */
00496 
00497      mmm[ijk] = 0 ;                                /* clear found point */
00498      nall = DALL ;                                 /* # allocated pts */
00499      nnow = 1 ;                                    /* # pts in cluster */
00500      inow = (short *) malloc(sizeof(short)*DALL) ; /* coords of pts */
00501      jnow = (short *) malloc(sizeof(short)*DALL) ;
00502      know = (short *) malloc(sizeof(short)*DALL) ;
00503      IJK_TO_THREE(ijk, inow[0],jnow[0],know[0] , nx,nxy) ;
00504 
00505      /*--
00506         for each point in cluster:
00507            check neighboring points for nonzero entries in mmm
00508            enter those into cluster (and clear them in mmm)
00509            continue until end of cluster is reached
00510              (note that cluster is expanding as we progress)
00511      --*/
00512 
00513      for( icl=0 ; icl < nnow ; icl++ ){
00514        ii = inow[icl] ; jj = jnow[icl] ; kk = know[icl] ;
00515        im = ii-1      ; jm = jj-1      ; km = kk-1 ;
00516        ip = ii+1      ; jp = jj+1      ; kp = kk+1 ;
00517 
00518        if( im >= 0 ) CPUT(im,jj,kk) ;
00519        if( ip < nx ) CPUT(ip,jj,kk) ;
00520        if( jm >= 0 ) CPUT(ii,jm,kk) ;
00521        if( jp < ny ) CPUT(ii,jp,kk) ;
00522        if( km >= 0 ) CPUT(ii,jj,km) ;
00523        if( kp < nz ) CPUT(ii,jj,kp) ;
00524      }
00525 
00526      /* see if now cluster is larger than best yet */
00527 
00528      if( nnow > nbest ){                         /* now is bigger; */
00529        free(ibest) ; free(jbest) ; free(kbest) ; /* replace best  */
00530        nbest = nnow ; ibest = inow ;             /* with now     */
00531        jbest = jnow ; kbest = know ;
00532      } else {                                    /* old is bigger */
00533        free(inow) ; free(jnow) ; free(know) ;    /* toss now     */
00534      }
00535 
00536    } /* loop ends when all nonzero points are clustered */
00537 
00538    /* put 1's back in at all points in best cluster */
00539 
00540    for( icl=0 ; icl < nbest ; icl++ ){
00541       ijk = THREE_TO_IJK(ibest[icl],jbest[icl],kbest[icl],nx,nxy) ;
00542       mmm[ijk] = 1 ;
00543    }
00544    free(ibest) ; free(jbest) ; free(kbest) ;
00545 
00546    if( verb ) ININFO_message("Largest cluster has %d voxels\n",nbest) ;
00547 
00548    return ;
00549 }
00550 
00551 /*--------------------------------------------------------------------------*/
00552 /*! Erode away nonzero voxels that aren't neighbored by mostly other
00553     nonzero voxels.  Then restore those that were eroded that are
00554     neighbors of survivors.  The neighbors are the 18 voxels closest
00555     in 3D (nearest and next-nearest neighbors).
00556 ----------------------------------------------------------------------------*/
00557 
00558 void THD_mask_erode( int nx, int ny, int nz, byte *mmm )
00559 {
00560    int ii,jj,kk , jy,kz, im,jm,km , ip,jp,kp , num ;
00561    int nxy=nx*ny , nxyz=nxy*nz ;
00562    byte *nnn ;
00563 
00564    if( mmm == NULL ) return ;
00565 
00566    nnn = (byte*)calloc(sizeof(byte),nxyz) ;  /* mask of eroded voxels */
00567 
00568    /* mark interior voxels that don't have 17 out of 18 nonzero nbhrs */
00569 
00570    for( kk=0 ; kk < nz ; kk++ ){
00571     kz = kk*nxy ; km = kz-nxy ; kp = kz+nxy ;
00572          if( kk == 0    ) km = kz ;
00573     else if( kk == nz-1 ) kp = kz ;
00574 
00575     for( jj=0 ; jj < ny ; jj++ ){
00576      jy = jj*nx ; jm = jy-nx ; jp = jy+nx ;
00577           if( jj == 0    ) jm = jy ;
00578      else if( jj == ny-1 ) jp = jy ;
00579 
00580      for( ii=0 ; ii < nx ; ii++ ){
00581        if( mmm[ii+jy+kz] ){           /* count nonzero nbhrs */
00582          im = ii-1 ; ip = ii+1 ;
00583               if( ii == 0    ) im = 0 ;
00584          else if( ii == nx-1 ) ip = ii ;
00585          num =  mmm[im+jy+km]
00586               + mmm[ii+jm+km] + mmm[ii+jy+km] + mmm[ii+jp+km]
00587               + mmm[ip+jy+km]
00588               + mmm[im+jm+kz] + mmm[im+jy+kz] + mmm[im+jp+kz]
00589               + mmm[ii+jm+kz]                 + mmm[ii+jp+kz]
00590               + mmm[ip+jm+kz] + mmm[ip+jy+kz] + mmm[ip+jp+kz]
00591               + mmm[im+jy+kp]
00592               + mmm[ii+jm+kp] + mmm[ii+jy+kp] + mmm[ii+jp+kp]
00593               + mmm[ip+jy+kp] ;
00594          if( num < 17 ) nnn[ii+jy+kz] = 1 ;  /* mark to erode */
00595        }
00596    } } }
00597 
00598    for( jj=ii=0 ; ii < nxyz ; ii++ )            /* actually erode */
00599      if( nnn[ii] ){ mmm[ii] = 0 ; jj++ ; }
00600 
00601    if( verb && jj > 0 ) ININFO_message("Eroded   %d voxels\n",jj) ;
00602 
00603    /* re-dilate eroded voxels that are next to survivors */
00604 
00605 #ifdef USE_DILATE
00606    for( kk=0 ; kk < nz ; kk++ ){
00607     kz = kk*nxy ; km = kz-nxy ; kp = kz+nxy ;
00608          if( kk == 0    ) km = kz ;
00609     else if( kk == nz-1 ) kp = kz ;
00610 
00611     for( jj=0 ; jj < ny ; jj++ ){
00612      jy = jj*nx ; jm = jy-nx ; jp = jy+nx ;
00613           if( jj == 0    ) jm = jy ;
00614      else if( jj == ny-1 ) jp = jy ;
00615 
00616      for( ii=0 ; ii < nx ; ii++ ){
00617        if( nnn[ii+jy+kz] ){           /* was eroded */
00618          im = ii-1 ; ip = ii+1 ;
00619               if( ii == 0    ) im = 0 ;
00620          else if( ii == nx-1 ) ip = ii ;
00621          nnn[ii+jy+kz] =              /* see if has any nbhrs */
00622                mmm[im+jy+km]
00623             || mmm[ii+jm+km] || mmm[ii+jy+km] || mmm[ii+jp+km]
00624             || mmm[ip+jy+km]
00625             || mmm[im+jm+kz] || mmm[im+jy+kz] || mmm[im+jp+kz]
00626             || mmm[ii+jm+kz]                  || mmm[ii+jp+kz]
00627             || mmm[ip+jm+kz] || mmm[ip+jy+kz] || mmm[ip+jp+kz]
00628             || mmm[im+jy+kp]
00629             || mmm[ii+jm+kp] || mmm[ii+jy+kp] || mmm[ii+jp+kp]
00630             || mmm[ip+jy+kp] ;
00631        }
00632    } } }
00633 
00634    /* actually do the dilation */
00635 
00636    for( jj=ii=0 ; ii < nxyz ; ii++ )
00637      if( nnn[ii] ){ mmm[ii] = 1 ; jj++ ; }
00638 
00639    if( verb && jj > 0 ) ININFO_message("Restored %d eroded voxels\n",jj) ;
00640 #endif
00641 
00642    free(nnn) ; return ;
00643 }
00644 
00645 /*--------------------------------------------------------------------------*/
00646 /*! Dilate a mask - that is, fill in zero voxels that have at least ndil
00647     neighbors in the mask.  The neighbors are the 18 voxels closest
00648     in 3D (nearest and next-nearest neighbors).
00649 ----------------------------------------------------------------------------*/
00650 
00651 void THD_mask_dilate( int nx, int ny, int nz, byte *mmm , int ndil )
00652 {
00653    int ii,jj,kk , jy,kz, im,jm,km , ip,jp,kp , num ;
00654    int nxy=nx*ny , nxyz=nxy*nz ;
00655    byte *nnn ;
00656 
00657    if( mmm == NULL ) return ;
00658         if( ndil < 1  ) ndil =  1 ;
00659    else if( ndil > 17 ) ndil = 17 ;
00660 
00661    nnn = (byte*)calloc(sizeof(byte),nxyz) ;  /* mask of dilated voxels */
00662 
00663    /* mark exterior voxels neighboring enough interior voxels */
00664 
00665    for( kk=0 ; kk < nz ; kk++ ){
00666     kz = kk*nxy ; km = kz-nxy ; kp = kz+nxy ;
00667          if( kk == 0    ) km = kz ;
00668     else if( kk == nz-1 ) kp = kz ;
00669 
00670     for( jj=0 ; jj < ny ; jj++ ){
00671      jy = jj*nx ; jm = jy-nx ; jp = jy+nx ;
00672           if( jj == 0    ) jm = jy ;
00673      else if( jj == ny-1 ) jp = jy ;
00674 
00675      for( ii=0 ; ii < nx ; ii++ ){
00676        if( mmm[ii+jy+kz] == 0 ){           /* count nonzero nbhrs */
00677          im = ii-1 ; ip = ii+1 ;
00678               if( ii == 0    ) im = 0 ;
00679          else if( ii == nx-1 ) ip = ii ;
00680          num =  mmm[im+jy+km]
00681               + mmm[ii+jm+km] + mmm[ii+jy+km] + mmm[ii+jp+km]
00682               + mmm[ip+jy+km]
00683               + mmm[im+jm+kz] + mmm[im+jy+kz] + mmm[im+jp+kz]
00684               + mmm[ii+jm+kz]                 + mmm[ii+jp+kz]
00685               + mmm[ip+jm+kz] + mmm[ip+jy+kz] + mmm[ip+jp+kz]
00686               + mmm[im+jy+kp]
00687               + mmm[ii+jm+kp] + mmm[ii+jy+kp] + mmm[ii+jp+kp]
00688               + mmm[ip+jy+kp] ;
00689          if( num >= ndil ) nnn[ii+jy+kz] = 1 ;  /* mark to dilate */
00690        }
00691    } } }
00692 
00693    for( ii=0 ; ii < nxyz ; ii++ )            /* actually dilate */
00694      if( nnn[ii] ) mmm[ii] = 1 ;
00695 
00696    free(nnn) ; return ;
00697 }
00698 
00699 /*---------------------------------------------------------------------*/
00700 /*! Find a bounding box for the main cluster of large-ish voxels.
00701     [xm..xp] will be box for x index, etc.
00702 -----------------------------------------------------------------------*/
00703 
00704 void THD_autobbox( THD_3dim_dataset *dset ,
00705                    int *xm, int *xp , int *ym, int *yp , int *zm, int *zp )
00706 {
00707    MRI_IMAGE *medim ;
00708    float clip_val , *mar ;
00709    int nvox , ii ;
00710 
00711 ENTRY("THD_autobbox") ;
00712 
00713    medim = THD_median_brick(dset) ; if( medim == NULL ) EXRETURN ;
00714 
00715    mar  = MRI_FLOAT_PTR(medim) ;
00716    nvox = medim->nvox ;
00717    for( ii=0 ; ii < nvox ; ii++ ) mar[ii] = fabs(mar[ii]) ;
00718 
00719    clip_val = THD_cliplevel(medim,0.5) ;
00720    for( ii=0 ; ii < nvox ; ii++ )
00721      if( mar[ii] < clip_val ) mar[ii] = 0.0 ;
00722 
00723    MRI_autobbox( medim , xm,xp , ym,yp , zm,zp ) ;
00724 
00725    mri_free(medim) ; EXRETURN ;
00726 }
00727 
00728 /*------------------------------------------------------------------------*/
00729 
00730 void MRI_autobbox( MRI_IMAGE *qim ,
00731                    int *xm, int *xp , int *ym, int *yp , int *zm, int *zp )
00732 {
00733    MRI_IMAGE *fim ;
00734    float *mar ;
00735    byte *mmm = NULL ;
00736    int nvox , ii,jj,kk , nmm , nx,ny,nz,nxy ;
00737 
00738 ENTRY("MRI_autobbox") ;
00739 
00740    /* find largest component as in first part of THD_automask() */
00741 
00742    if( qim->kind != MRI_float ) fim = mri_to_float(qim) ;
00743    else                         fim = qim ;
00744 
00745    nvox = fim->nvox ;
00746    mar  = MRI_FLOAT_PTR(fim) ;
00747    mmm  = (byte *) calloc( sizeof(byte) , nvox ) ;
00748    for( nmm=ii=0 ; ii < nvox ; ii++ )
00749      if( mar[ii] != 0.0 ){ mmm[ii] = 1; nmm++; }
00750 
00751    if( fim != qim ) mri_free(fim) ;
00752 
00753    if( nmm == 0 ){ free(mmm); EXRETURN; }
00754 
00755    nx = qim->nx; ny = qim->ny; nz = qim->nz; nxy = nx*ny;
00756 
00757    THD_mask_clust( nx,ny,nz, mmm ) ;
00758    THD_mask_erode( nx,ny,nz, mmm ) ;
00759    THD_mask_clust( nx,ny,nz, mmm ) ;
00760 
00761    /* For each plane direction,
00762       find the first and last index that have nonzero voxels in that plane */
00763 
00764    if( xm != NULL ){
00765      for( ii=0 ; ii < nx ; ii++ )
00766        for( kk=0 ; kk < nz ; kk++ )
00767          for( jj=0 ; jj < ny ; jj++ )
00768            if( mmm[ii+jj*nx+kk*nxy] ) goto CP5 ;
00769      CP5: *xm = ii ;
00770    }
00771 
00772    if( xp != NULL ){
00773      for( ii=nx-1 ; ii >= 0 ; ii-- )
00774        for( kk=0 ; kk < nz ; kk++ )
00775          for( jj=0 ; jj < ny ; jj++ )
00776            if( mmm[ii+jj*nx+kk*nxy] ) goto CP6 ;
00777      CP6: *xp = ii ;
00778    }
00779 
00780    if( ym != NULL ){
00781      for( jj=0 ; jj < ny ; jj++ )
00782        for( kk=0 ; kk < nz ; kk++ )
00783          for( ii=0 ; ii < nx ; ii++ )
00784            if( mmm[ii+jj*nx+kk*nxy] ) goto CP3 ;
00785      CP3: *ym = jj ;
00786    }
00787 
00788    if( yp != NULL ){
00789      for( jj=ny-1 ; jj >= 0 ; jj-- )
00790        for( kk=0 ; kk < nz ; kk++ )
00791          for( ii=0 ; ii < nx ; ii++ )
00792            if( mmm[ii+jj*nx+kk*nxy] ) goto CP4 ;
00793      CP4: *yp = jj ;
00794    }
00795 
00796    if( zm != NULL ){
00797      for( kk=0 ; kk < nz ; kk++ )
00798        for( jj=0 ; jj < ny ; jj++ )
00799          for( ii=0 ; ii < nx ; ii++ )
00800            if( mmm[ii+jj*nx+kk*nxy] ) goto CP1 ;
00801      CP1: *zm = kk ;
00802    }
00803 
00804    if( zp != NULL ){
00805      for( kk=nz-1 ; kk >= 0 ; kk-- )
00806        for( jj=0 ; jj < ny ; jj++ )
00807          for( ii=0 ; ii < nx ; ii++ )
00808            if( mmm[ii+jj*nx+kk*nxy] ) goto CP2 ;
00809      CP2: *zp = kk ;
00810    }
00811 
00812    free(mmm) ; EXRETURN ;
00813 }
00814 
00815 /*------------------------------------------------------------------------*/
00816 
00817 int THD_peel_mask( int nx, int ny, int nz , byte *mmm, int pdepth )
00818 {
00819    int nxy=nx*ny , ii,jj,kk , ijk , bot,top , pd=pdepth ;
00820    int nxp=nx-pd , nyp=ny-pd , nzp=nz-pd ;
00821    int num=0 , dnum , nite ;
00822 
00823    for( nite=0 ; nite < pd ; nite++ ){
00824     dnum = 0 ;
00825 
00826     for( kk=0 ; kk < nz ; kk++ ){
00827      for( jj=0 ; jj < ny ; jj++ ){
00828        ijk = jj*nx + kk*nxy ;
00829        for( bot=0 ; bot < nx && !mmm[bot+ijk]; bot++ ) ;
00830        top = bot+pd ; if( top >= nx ) continue ;
00831        for( ii=bot+1 ; ii <= top && mmm[ii+ijk] ; ii++ ) ;
00832        if( ii <= top ){ mmm[bot+ijk] = 0; dnum++; }
00833     }}
00834 
00835     for( kk=0 ; kk < nz ; kk++ ){
00836      for( jj=0 ; jj < ny ; jj++ ){
00837        ijk = jj*nx + kk*nxy ;
00838        for( top=nx-1 ; top >= 0 && !mmm[top+ijk]; top-- ) ;
00839        bot = top-pd ; if( bot < 0 ) continue ;
00840        for( ii=top-1 ; ii >= bot && mmm[ii+ijk] ; ii-- ) ;
00841        if( ii >= bot ){ mmm[top+ijk] = 0; dnum++; }
00842     }}
00843 
00844     for( kk=0 ; kk < nz ; kk++ ){
00845      for( ii=0 ; ii < nx ; ii++ ){
00846        ijk = ii + kk*nxy ;
00847        for( bot=0 ; bot < ny && !mmm[bot*nx+ijk] ; bot++ ) ;
00848        top = bot+pd ;
00849        if( top >= ny ) continue ;
00850        for( jj=bot+1 ; jj <= top && mmm[jj*nx+ijk] ; jj++ ) ;
00851        if( jj <= top ){ mmm[bot*nx+ijk] = 0; dnum++; }
00852     }}
00853 
00854     for( kk=0 ; kk < nz ; kk++ ){
00855      for( ii=0 ; ii < nx ; ii++ ){
00856        ijk = ii + kk*nxy ;
00857        for( top=ny-1 ; top >= 0 && !mmm[top*nx+ijk] ; top-- ) ;
00858        bot = top-pd ; if( bot < 0 ) continue ;
00859        for( jj=top-1 ; jj >= bot && mmm[jj*nx+ijk] ; jj-- ) ;
00860        if( jj >= bot ){ mmm[top*nx+ijk] = 0; dnum++; }
00861     }}
00862 
00863     for( jj=0 ; jj < ny ; jj++ ){
00864      for( ii=0 ; ii < nx ; ii++ ){
00865        ijk = ii + jj*nx ;
00866        for( top=nz-1 ; top >= 0 && !mmm[top*nxy+ijk] ; top-- ) ;
00867        bot = top-pd ; if( bot < 0 ) continue ;
00868        for( kk=top-1 ; kk >= bot && mmm[kk*nxy+ijk] ; kk-- ) ;
00869        if( kk >= bot ){ mmm[top*nxy+ijk] = 0; dnum++; }
00870     }}
00871 
00872     num += dnum ;
00873     if( dnum == 0 ) break ;
00874    }
00875 
00876    return num ;
00877 }
 
Powered by Plone

This site conforms to the following standards: