00001 #include "mrilib.h"
00002 
00003 /*----------------------------------------------------------------*/
00004 /*! Modified MCW_find_clusters():
00005     - added mode variable
00006     - ISOVALUE_MODE => make clusters that are both contiguous
00007                        and have the same numerical value
00008     - ISOMERGE_MODE => make clusters that have the same numerical
00009                        value, regardless of contiguity
00010 ------------------------------------------------------------------*/
00011 
00012 MCW_cluster_array * NIH_find_clusters(
00013                        int nx, int ny, int nz,
00014                        float dx, float dy, float dz,
00015                        int ftype , void * fim ,
00016                        float max_dist , int mode )
00017 {
00018    MCW_cluster_array *clust_arr ;
00019    MCW_cluster       *clust , *mask=NULL ;
00020    int ii,jj,kk ,  nxy,nxyz , ijk , ijk_last , mnum ;
00021    int icl , jma , ijkcl , ijkma , did_one ;
00022    float fimv ;
00023    short *sfar ;
00024    float *ffar ;
00025    byte  *bfar ;
00026    short ic, jc, kc;
00027    short im, jm, km;
00028 
00029 ENTRY("NIH_find_clusters") ;
00030 
00031    if( fim == NULL ) RETURN(NULL) ;
00032 
00033    switch( ftype ){
00034       default: RETURN(NULL) ;
00035       case MRI_short:  sfar = (short *) fim ; break ;
00036       case MRI_byte :  bfar = (byte  *) fim ; break ;
00037       case MRI_float:  ffar = (float *) fim ; break ;
00038    }
00039 
00040    /* default => use older code (in edt_clust.c) */
00041 
00042    if( mode <= 0 || mode > ISOMERGE_MODE ){
00043      RETURN( MCW_find_clusters( nx,ny,nz , dx,dy,dz ,
00044                                ftype,fim , max_dist ) ) ;
00045    }
00046 
00047    /*--- make a cluster that is a mask of points closer than max_dist ---*/
00048 
00049    if( mode == ISOVALUE_MODE ){
00050      mask = MCW_build_mask (nx, ny, nz, dx, dy, dz, max_dist);
00051      if (mask == NULL)
00052      {
00053         fprintf(stderr, "Unable to build mask in NIH_find_clusters");
00054         RETURN(NULL);
00055      }
00056      mnum = mask->num_pt ;
00057    }
00058 
00059    nxy = nx*ny ; nxyz = nxy * nz ;
00060 
00061    /*--- scan through array, find nonzero point, build a cluster, ... ---*/
00062 
00063    INIT_CLARR(clust_arr) ;
00064 
00065    ijk_last = 0 ;
00066    do {
00067 
00068       /* find nonzero point in 3D array, starting at ijk_last */
00069 
00070       switch( ftype ){
00071          case MRI_short:
00072             for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( sfar[ijk] != 0 ) break ;
00073             if( ijk < nxyz ){
00074                fimv = sfar[ijk] ; sfar[ijk] = 0 ;  /* save found point */
00075             }
00076          break ;
00077 
00078          case MRI_byte:
00079             for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( bfar[ijk] != 0 ) break ;
00080             if( ijk < nxyz ){
00081                fimv = bfar[ijk] ; bfar[ijk] = 0 ;  /* save found point */
00082             }
00083          break ;
00084 
00085          case MRI_float:
00086             for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( ffar[ijk] != 0.0 ) break ;
00087             if( ijk < nxyz ){
00088                fimv = ffar[ijk] ; ffar[ijk] = 0.0 ;  /* save found point */
00089             }
00090          break ;
00091       }
00092       if( ijk == nxyz ) break ;  /* didn't find any nonzero point! */
00093 
00094       ijk_last = ijk+1 ;         /* start here next time */
00095 
00096       INIT_CLUSTER(clust) ;                        /* make a new cluster */
00097       IJK_TO_THREE(ijk,ic,jc,kc,nx,nxy) ;          /* find 3D index */
00098       ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ;  /* start cluster off */
00099 
00100       switch( mode ){
00101 
00102         case ISOVALUE_MODE:{ /* for each point in cluster:
00103                                 check points offset by mask for nonzero entries in fim
00104                                 enter those into cluster
00105                                 continue until end of cluster is reached
00106                                 (note that cluster is expanding as we progress) */
00107 
00108          switch( ftype ){
00109             case MRI_short:
00110                for( icl=0 ; icl < clust->num_pt ; icl++ ){
00111                   ic = clust->i[icl];  /* check around this point */
00112                   jc = clust->j[icl];
00113                   kc = clust->k[icl];
00114 
00115                   for( jma=0 ; jma < mnum ; jma++ ){
00116                      im = ic + mask->i[jma];  /* offset by mask */
00117                      jm = jc + mask->j[jma];
00118                      km = kc + mask->k[jma];
00119                      if( im < 0 || im >= nx ||
00120                          jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
00121 
00122                      ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
00123                      if( ijkma < ijk_last || ijkma >= nxyz || sfar[ijkma] != fimv ) continue ;
00124 
00125                      ADDTO_CLUSTER( clust , im, jm, km, sfar[ijkma] ) ;
00126                      sfar[ijkma] = 0 ;
00127                   }
00128                }
00129             break ;
00130 
00131             case MRI_byte:
00132                for( icl=0 ; icl < clust->num_pt ; icl++ ){
00133                   ic = clust->i[icl];
00134                   jc = clust->j[icl];
00135                   kc = clust->k[icl];
00136 
00137                   for( jma=0 ; jma < mnum ; jma++ ){
00138                      im = ic + mask->i[jma];
00139                      jm = jc + mask->j[jma];
00140                      km = kc + mask->k[jma];
00141                      if( im < 0 || im >= nx ||
00142                          jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
00143 
00144                      ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
00145                      if( ijkma < ijk_last || ijkma >= nxyz || bfar[ijkma] != fimv ) continue ;
00146 
00147                      ADDTO_CLUSTER( clust , im, jm, km, bfar[ijkma] ) ;
00148                      bfar[ijkma] = 0 ;
00149                   }
00150                }
00151             break ;
00152 
00153             case MRI_float:
00154                for( icl=0 ; icl < clust->num_pt ; icl++ ){
00155                   ic = clust->i[icl];
00156                   jc = clust->j[icl];
00157                   kc = clust->k[icl];
00158 
00159                   for( jma=0 ; jma < mnum ; jma++ ){
00160                      im = ic + mask->i[jma];
00161                      jm = jc + mask->j[jma];
00162                      km = kc + mask->k[jma];
00163                      if( im < 0 || im >= nx ||
00164                          jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
00165 
00166                      ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
00167                      if( ijkma < ijk_last || ijkma >= nxyz || ffar[ijkma] != fimv ) continue ;
00168 
00169                      ADDTO_CLUSTER( clust , im, jm, km, ffar[ijkma] ) ;
00170                      ffar[ijkma] = 0.0 ;
00171                   }
00172                }
00173             break ;
00174          } /* end of switch on array type */
00175         }
00176         break ; /* end of ISOVALUE_MODE */
00177 
00178         /*.................................*/
00179 
00180         case ISOMERGE_MODE:{  /* find other points with the same value */
00181 
00182          switch( ftype ){
00183             case MRI_short:
00184               for( ijk=ijk_last ; ijk < nxyz ; ijk++ )
00185                 if( sfar[ijk] == fimv ){
00186                   IJK_TO_THREE(ijk,ic,jc,kc,nx,nxy) ;          /* find 3D index */
00187                   ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ;  /* start cluster off */
00188                   sfar[ijk] = 0 ;
00189                 }
00190             break ;
00191 
00192             case MRI_byte:
00193               for( ijk=ijk_last ; ijk < nxyz ; ijk++ )
00194                 if( bfar[ijk] == fimv ){
00195                   IJK_TO_THREE(ijk,ic,jc,kc,nx,nxy) ;          /* find 3D index */
00196                   ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ;  /* start cluster off */
00197                   bfar[ijk] = 0 ;
00198                 }
00199             break ;
00200 
00201             case MRI_float:
00202               for( ijk=ijk_last ; ijk < nxyz ; ijk++ )
00203                 if( ffar[ijk] == fimv ){
00204                   IJK_TO_THREE(ijk,ic,jc,kc,nx,nxy) ;          /* find 3D index */
00205                   ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ;  /* start cluster off */
00206                   ffar[ijk] = 0 ;
00207                 }
00208             break ;
00209          }
00210         }
00211         break ; /* end of ISOMERGE_MODE */
00212 
00213       }
00214 
00215       ADDTO_CLARR(clust_arr,clust) ;
00216    } while( 1 ) ;
00217 
00218    if( mask != NULL ) KILL_CLUSTER(mask) ;
00219 
00220    if( clust_arr->num_clu <= 0 ){ DESTROY_CLARR(clust_arr) ; }
00221 
00222    RETURN(clust_arr) ;
00223 }