00001 #include "mrilib.h"
00002
00003
00004
00005
00006
00007
00008
00009
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
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
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
00062
00063 INIT_CLARR(clust_arr) ;
00064
00065 ijk_last = 0 ;
00066 do {
00067
00068
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 ;
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 ;
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 ;
00089 }
00090 break ;
00091 }
00092 if( ijk == nxyz ) break ;
00093
00094 ijk_last = ijk+1 ;
00095
00096 INIT_CLUSTER(clust) ;
00097 IJK_TO_THREE(ijk,ic,jc,kc,nx,nxy) ;
00098 ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ;
00099
00100 switch( mode ){
00101
00102 case ISOVALUE_MODE:{
00103
00104
00105
00106
00107
00108 switch( ftype ){
00109 case MRI_short:
00110 for( icl=0 ; icl < clust->num_pt ; icl++ ){
00111 ic = clust->i[icl];
00112 jc = clust->j[icl];
00113 kc = clust->k[icl];
00114
00115 for( jma=0 ; jma < mnum ; jma++ ){
00116 im = ic + mask->i[jma];
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 }
00175 }
00176 break ;
00177
00178
00179
00180 case ISOMERGE_MODE:{
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) ;
00187 ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ;
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) ;
00196 ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ;
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) ;
00205 ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ;
00206 ffar[ijk] = 0 ;
00207 }
00208 break ;
00209 }
00210 }
00211 break ;
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 }