1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 /*-----------------------------------------------------------------------------*/
10 /*! Routine to make a cluster that is a mask of points closer than max_dist.
11 - dx, dy, dz = voxel dimensions
12 - max_dist = maximum distance for a point to be included in the mask
13 Date : 11 September 1996
14
15 To correct error due to ambiguity in identification of clusters,
16 voxel coordinates are now stored as 3 separate short integers.
17 BDW 06 March 1997
18
19 30 Apr 2002: max_dist input as <= 0 now gives NN connectivity
20
21 N.B.: The cluster does NOT contain the (0,0,0) point!
22 N.B.: The cluster is not sorted by radius from the (0,0,0) point!
23 To do this, see MCW_radsort_cluster().
24 -----------------------------------------------------------------------------*/
25
MCW_build_mask(float dx,float dy,float dz,float max_dist)26 MCW_cluster * MCW_build_mask( float dx, float dy, float dz, float max_dist )
27 {
28 int ii, jj, kk, idx, jdy, kdz, mnum;
29 float xq, yq, zq, dist_q;
30 MCW_cluster *mask;
31
32 ENTRY("MCW_build_mask") ;
33
34 if( max_dist <= 0.0 ){ /* 30 Apr 2002 */
35 dx = dy = dz = 1.0f ; max_dist = 1.01f ;
36 } else {
37 if( dx <= 0.0f ) dx = 1.0f ; /* something sensible */
38 if( dy <= 0.0f ) dy = 1.0f ;
39 if( dz <= 0.0f ) dz = 1.0f ;
40 }
41
42 idx = max_dist/dx ; jdy = max_dist/dy ; kdz = max_dist/dz ;
43
44 if( idx < 1 && jdy < 1 && kdz < 1 ){
45 WARNING_message("Illegal input to MCW_build_mask:"
46 " dx=%g dy=%g dz=%g max_dist=%g" ,
47 dx,dy,dz,max_dist ) ;
48 RETURN( NULL );
49 }
50
51 INIT_CLUSTER(mask) ;
52
53 dist_q = max_dist * max_dist ;
54
55 for( kk=-kdz ; kk <= kdz ; kk++ ){
56 zq = (kk*dz) * (kk*dz) ;
57 for( jj=-jdy ; jj <= jdy ; jj++ ){
58 yq = zq + (jj*dy) * (jj*dy) ;
59 for( ii=-idx ; ii <= idx ; ii++ ){
60 xq = yq + (ii*dx)*(ii*dx) ;
61 if( xq <= dist_q && xq > 0.0f ){
62 ADDTO_CLUSTER( mask , ii, jj, kk, 0 ) ;
63 }
64 }
65 }
66 }
67
68 mnum = mask->num_pt ;
69 if( mnum < 1 ){
70 KILL_CLUSTER(mask) ;
71 WARNING_message("MCW_build_mask error: mask has only %d elements!",mnum);
72 RETURN( NULL );
73 }
74
75 RETURN (mask);
76 }
77
78 /*----------------------------------------------------------------------*/
79 /*! Like MCW_build_mask(), but adds the (0,0,0) point to the cluster,
80 and then sorts by radius (so the points closest to the origin
81 are first in the array). The radius from the origin is stored
82 in the 'mag' element of the cluster struct.
83 ------------------------------------------------------------------------*/
84
MCW_spheremask(float dx,float dy,float dz,float radius)85 MCW_cluster * MCW_spheremask( float dx, float dy, float dz, float radius )
86 {
87 MCW_cluster *mask=NULL;
88 int ii , nn ;
89 float x,y,z ;
90
91 if( radius != 0.0f )
92 mask = MCW_build_mask(dx,dy,dz,radius) ;
93
94 if( mask == NULL ){ INIT_CLUSTER(mask) ; }
95 ADDTO_CLUSTER(mask,0,0,0,0) ;
96
97 /** sorting stuff added 20 Oct 2006 **/
98
99 if( dx <= 0.0f ) dx = 1.0f ;
100 if( dy <= 0.0f ) dy = 1.0f ;
101 if( dz <= 0.0f ) dz = 1.0f ;
102 nn = mask->num_pt ;
103 for( ii=0 ; ii < nn ; ii++ ){
104 x = mask->i[ii]*dx; y = mask->j[ii]*dy; z = mask->k[ii]*dz;
105
106 /* if((x==0.0)&&(y==0.0)&&(z==0.0))
107 mask->mag[ii] = 1.0;
108 else
109 */
110
111 mask->mag[ii] = sqrt(x*x+y*y+z*z) ;
112 }
113 MCW_sort_cluster( mask ) ;
114 return mask ;
115 }
116
117 /*----------------------------------------------------------------------*/
118 /*! Like MCW_spheremask(), but builds a rectangular parallelopiped. */
119
MCW_rectmask(float dx,float dy,float dz,float xh,float yh,float zh)120 MCW_cluster * MCW_rectmask( float dx, float dy, float dz,
121 float xh, float yh, float zh )
122 {
123 int ii, jj, kk, idx, jdy, kdz ;
124 MCW_cluster *mask;
125
126 if( dx <= 0.0f ) dx = 1.0f ;
127 if( dy <= 0.0f ) dy = 1.0f ;
128 if( dz <= 0.0f ) dz = 1.0f ;
129 if( xh < 0.0f ) xh = 0.0f ;
130 if( yh < 0.0f ) yh = 0.0f ;
131 if( zh < 0.0f ) zh = 0.0f ;
132
133 idx = (int)(xh/dx) ;
134 jdy = (int)(yh/dy) ;
135 kdz = (int)(zh/dz) ;
136
137 INIT_CLUSTER(mask) ;
138
139 ADDTO_CLUSTER(mask,0,0,0,0) ; /* always keep central point first */
140 for( kk=-kdz ; kk <= kdz ; kk++ ){
141 for( jj=-jdy ; jj <= jdy ; jj++ ){
142 for( ii=-idx ; ii <= idx ; ii++ ){
143 if (ii || jj || kk)
144 ADDTO_CLUSTER( mask , ii,jj,kk , 0 ) ;
145 }}}
146
147 return mask ;
148 }
149
150 /*----------------------------------------------------------------------*/
151 /*! Like MCW_spheremask(), but builds a rhombic dodecahedron.
152 Volume = 2 * radius**3.
153 *//*--------------------------------------------------------------------*/
154
MCW_rhddmask(float dx,float dy,float dz,float radius)155 MCW_cluster * MCW_rhddmask( float dx, float dy, float dz, float radius )
156 {
157 int ii, jj, kk, idx, jdy, kdz ;
158 float a,b,c ;
159 MCW_cluster *mask;
160
161 if( radius <= 0.0 ){ /* 30 Apr 2002 */
162 dx = dy = dz = 1.0f ; radius = 1.01f ;
163 } else {
164 if( dx <= 0.0f ) dx = 1.0f ;
165 if( dy <= 0.0f ) dy = 1.0f ;
166 if( dz <= 0.0f ) dz = 1.0f ;
167 }
168
169 idx = (int)(radius/dx) ;
170 jdy = (int)(radius/dy) ;
171 kdz = (int)(radius/dz) ;
172
173 INIT_CLUSTER(mask) ;
174
175 ADDTO_CLUSTER(mask,0,0,0,0) ; /* always keep central point first */
176 for( kk=-kdz ; kk <= kdz ; kk++ ){
177 c = kk*dz ;
178 for( jj=-jdy ; jj <= jdy ; jj++ ){
179 b = jj*dy ;
180 for( ii=-idx ; ii <= idx ; ii++ ){
181 if (ii || jj || kk) {
182 a = ii*dx ;
183 if( fabsf(a+b) <= radius &&
184 fabsf(a-b) <= radius &&
185 fabsf(a+c) <= radius &&
186 fabsf(a-c) <= radius &&
187 fabsf(b+c) <= radius &&
188 fabsf(b-c) <= radius ) ADDTO_CLUSTER( mask , ii,jj,kk , 0 ) ;
189 }
190 }}}
191
192 return mask ;
193 }
194
195 /*----------------------------------------------------------------------*/
196 /*! Like MCW_spheremask(), but builds a truncated octahedron.
197 Volume = 4 * radius**3.
198 *//*--------------------------------------------------------------------*/
199
200 #undef FAS
201 #undef TOHD_inside
202 #define FAS(a,b) (fabsf(a) <= (b))
203 #define TOHD_inside(a,b,c,siz) \
204 ( FAS((a),(siz)) && FAS((b),(siz)) && FAS((c),(siz)) && \
205 FAS((a)+(b)+(c),1.5f*(siz)) && FAS((a)-(b)+(c),1.5f*(siz)) && \
206 FAS((a)+(b)-(c),1.5f*(siz)) && FAS((a)-(b)-(c),1.5f*(siz)) )
207
MCW_tohdmask(float dx,float dy,float dz,float radius)208 MCW_cluster * MCW_tohdmask( float dx, float dy, float dz, float radius )
209 {
210 int ii, jj, kk, idx, jdy, kdz ;
211 float a,b,c ;
212 MCW_cluster *mask;
213
214 if( radius <= 0.0 ){ /* 30 Apr 2002 */
215 dx = dy = dz = 1.0f ; radius = 1.01f ;
216 } else {
217 if( dx <= 0.0f ) dx = 1.0f ;
218 if( dy <= 0.0f ) dy = 1.0f ;
219 if( dz <= 0.0f ) dz = 1.0f ;
220 }
221
222 idx = (int)(radius/dx) ;
223 jdy = (int)(radius/dy) ;
224 kdz = (int)(radius/dz) ;
225
226 INIT_CLUSTER(mask) ;
227
228 ADDTO_CLUSTER(mask,0,0,0,0) ; /* always keep central point first */
229 for( kk=-kdz ; kk <= kdz ; kk++ ){
230 c = kk*dz ;
231 for( jj=-jdy ; jj <= jdy ; jj++ ){
232 b = jj*dy ;
233 for( ii=-idx ; ii <= idx ; ii++ ){
234 a = ii*dx ;
235 if( (ii || jj || kk) && TOHD_inside(a,b,c,radius) )
236 ADDTO_CLUSTER( mask , ii,jj,kk , 0 ) ;
237 }}}
238
239 return mask ;
240 }
241
MCW_showmask(MCW_cluster * nbhd,char * opening,char * closing,FILE * fout)242 void MCW_showmask (MCW_cluster *nbhd, char *opening, char *closing, FILE *fout)
243 {
244 int ii;
245 if (!fout) fout = stdout;
246 if (opening) fprintf(fout, "%s", opening);
247 if (!nbhd) {
248 fprintf(fout, "NULL nbhd\n");
249 } else {
250 fprintf(fout, "Neighborhood of %d voxels (%d allocated), %s mag.\n",
251 nbhd->num_pt, nbhd->num_all, nbhd->mag?"with":"without");
252 if (nbhd->mag) {
253 for (ii=0; ii<nbhd->num_pt; ++ii) {
254 fprintf (fout, "Offset[I J K]: %+03d %+03d %+03d, Mag: %f\n",
255 nbhd->i[ii], nbhd->j[ii], nbhd->k[ii], nbhd->mag[ii]);
256 }
257 } else {
258 for (ii=0; ii<nbhd->num_pt; ++ii) {
259 fprintf (fout, "Offset[I J K]: %+03d %+03d %+03d\n",
260 nbhd->i[ii], nbhd->j[ii], nbhd->k[ii]);
261 }
262 }
263 }
264 if (closing) fprintf(fout, "%s", closing);
265 return;
266 }
267