1 #include "mrilib.h"
2
3 /*----------------------------------------------------------------*/
4 /*! Modified MCW_find_clusters():
5 - added mode variable
6 - ISOVALUE_MODE => make clusters that are both contiguous
7 and have the same numerical value
8 - ISOMERGE_MODE => make clusters that have the same numerical
9 value, regardless of contiguity
10 ------------------------------------------------------------------*/
11
NIH_find_clusters(int nx,int ny,int nz,float dx,float dy,float dz,int ftype,void * fim,float max_dist,int mode)12 MCW_cluster_array * NIH_find_clusters(
13 int nx, int ny, int nz,
14 float dx, float dy, float dz,
15 int ftype , void * fim ,
16 float max_dist , int mode )
17 {
18 MCW_cluster_array *clust_arr ;
19 MCW_cluster *clust , *mask=NULL ;
20 int ii,jj,kk , nxy,nxyz , ijk=0 , ijk_last , mnum=0 ;
21 int icl , jma , ijkcl , ijkma , did_one ;
22 float fimv=0.0 ;
23 short *sfar=NULL ;
24 float *ffar=NULL ;
25 byte *bfar=NULL ;
26 short ic, jc, kc;
27 short im, jm, km;
28
29 ENTRY("NIH_find_clusters") ;
30
31 if( fim == NULL ) RETURN(NULL) ;
32
33 switch( ftype ){
34 default: RETURN(NULL) ;
35 case MRI_short: sfar = (short *) fim ; break ;
36 case MRI_byte : bfar = (byte *) fim ; break ;
37 case MRI_float: ffar = (float *) fim ; break ;
38 }
39
40 /* default => use older code (in edt_clust.c) */
41
42 if( mode <= 0 || mode > ISOMERGE_MODE ){
43 RETURN( MCW_find_clusters( nx,ny,nz , dx,dy,dz ,
44 ftype,fim , max_dist ) ) ;
45 }
46
47 /*--- make a cluster that is a mask of points closer than max_dist ---*/
48
49 if( mode == ISOVALUE_MODE ){
50 mask = MCW_build_mask (dx, dy, dz, max_dist);
51 if (mask == NULL)
52 {
53 fprintf(stderr, "Unable to build mask in NIH_find_clusters");
54 RETURN(NULL);
55 }
56 mnum = mask->num_pt ;
57 }
58
59 nxy = nx*ny ; nxyz = nxy * nz ;
60
61 /*--- scan through array, find nonzero point, build a cluster, ... ---*/
62
63 INIT_CLARR(clust_arr) ;
64
65 ijk_last = 0 ;
66 do {
67
68 /* find nonzero point in 3D array, starting at ijk_last */
69
70 switch( ftype ){
71 case MRI_short:
72 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( sfar[ijk] != 0 ) break ;
73 if( ijk < nxyz ){
74 fimv = sfar[ijk] ; sfar[ijk] = 0 ; /* save found point */
75 }
76 break ;
77
78 case MRI_byte:
79 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( bfar[ijk] != 0 ) break ;
80 if( ijk < nxyz ){
81 fimv = bfar[ijk] ; bfar[ijk] = 0 ; /* save found point */
82 }
83 break ;
84
85 case MRI_float:
86 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( ffar[ijk] != 0.0 ) break ;
87 if( ijk < nxyz ){
88 fimv = ffar[ijk] ; ffar[ijk] = 0.0 ; /* save found point */
89 }
90 break ;
91 }
92 if( ijk == nxyz ) break ; /* didn't find any nonzero point! */
93
94 ijk_last = ijk+1 ; /* start here next time */
95
96 INIT_CLUSTER(clust) ; /* make a new cluster */
97 IJK_TO_THREE(ijk,ic,jc,kc,nx,nxy) ; /* find 3D index */
98 ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ; /* start cluster off */
99
100 switch( mode ){
101
102 case ISOVALUE_MODE:{ /* for each point in cluster:
103 check points offset by mask for nonzero entries in fim
104 enter those into cluster
105 continue until end of cluster is reached
106 (note that cluster is expanding as we progress) */
107
108 switch( ftype ){
109 case MRI_short:
110 for( icl=0 ; icl < clust->num_pt ; icl++ ){
111 ic = clust->i[icl]; /* check around this point */
112 jc = clust->j[icl];
113 kc = clust->k[icl];
114
115 for( jma=0 ; jma < mnum ; jma++ ){
116 im = ic + mask->i[jma]; /* offset by mask */
117 jm = jc + mask->j[jma];
118 km = kc + mask->k[jma];
119 if( im < 0 || im >= nx ||
120 jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
121
122 ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
123 if( ijkma < ijk_last || ijkma >= nxyz || sfar[ijkma] != fimv ) continue ;
124
125 ADDTO_CLUSTER( clust , im, jm, km, sfar[ijkma] ) ;
126 sfar[ijkma] = 0 ;
127 }
128 }
129 break ;
130
131 case MRI_byte:
132 for( icl=0 ; icl < clust->num_pt ; icl++ ){
133 ic = clust->i[icl];
134 jc = clust->j[icl];
135 kc = clust->k[icl];
136
137 for( jma=0 ; jma < mnum ; jma++ ){
138 im = ic + mask->i[jma];
139 jm = jc + mask->j[jma];
140 km = kc + mask->k[jma];
141 if( im < 0 || im >= nx ||
142 jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
143
144 ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
145 if( ijkma < ijk_last || ijkma >= nxyz || bfar[ijkma] != fimv ) continue ;
146
147 ADDTO_CLUSTER( clust , im, jm, km, bfar[ijkma] ) ;
148 bfar[ijkma] = 0 ;
149 }
150 }
151 break ;
152
153 case MRI_float:
154 for( icl=0 ; icl < clust->num_pt ; icl++ ){
155 ic = clust->i[icl];
156 jc = clust->j[icl];
157 kc = clust->k[icl];
158
159 for( jma=0 ; jma < mnum ; jma++ ){
160 im = ic + mask->i[jma];
161 jm = jc + mask->j[jma];
162 km = kc + mask->k[jma];
163 if( im < 0 || im >= nx ||
164 jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
165
166 ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
167 if( ijkma < ijk_last || ijkma >= nxyz || ffar[ijkma] != fimv ) continue ;
168
169 ADDTO_CLUSTER( clust , im, jm, km, ffar[ijkma] ) ;
170 ffar[ijkma] = 0.0 ;
171 }
172 }
173 break ;
174 } /* end of switch on array type */
175 }
176 break ; /* end of ISOVALUE_MODE */
177
178 /*.................................*/
179
180 case ISOMERGE_MODE:{ /* find other points with the same value */
181
182 switch( ftype ){
183 case MRI_short:
184 for( ijk=ijk_last ; ijk < nxyz ; ijk++ )
185 if( sfar[ijk] == fimv ){
186 IJK_TO_THREE(ijk,ic,jc,kc,nx,nxy) ; /* find 3D index */
187 ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ; /* start cluster off */
188 sfar[ijk] = 0 ;
189 }
190 break ;
191
192 case MRI_byte:
193 for( ijk=ijk_last ; ijk < nxyz ; ijk++ )
194 if( bfar[ijk] == fimv ){
195 IJK_TO_THREE(ijk,ic,jc,kc,nx,nxy) ; /* find 3D index */
196 ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ; /* start cluster off */
197 bfar[ijk] = 0 ;
198 }
199 break ;
200
201 case MRI_float:
202 for( ijk=ijk_last ; ijk < nxyz ; ijk++ )
203 if( ffar[ijk] == fimv ){
204 IJK_TO_THREE(ijk,ic,jc,kc,nx,nxy) ; /* find 3D index */
205 ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ; /* start cluster off */
206 ffar[ijk] = 0 ;
207 }
208 break ;
209 }
210 }
211 break ; /* end of ISOMERGE_MODE */
212
213 }
214
215 ADDTO_CLARR(clust_arr,clust) ;
216 } while( 1 ) ;
217
218 if( mask != NULL ) KILL_CLUSTER(mask) ;
219
220 if( clust_arr->num_clu <= 0 ){ DESTROY_CLARR(clust_arr) ; }
221
222 RETURN(clust_arr) ;
223 }
224