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 find clusters of points !=0; return an array of clusters
11 [the input array fim is destroyed in this
12 computation; it may be recreated by MCW_cluster_to_vol]
13
14 Nov 1995: fim array may now be short, byte, or float,
15 signified by new ftype argument.
16
17 Coordinates of voxels in clusters are now stored as 3 separate
18 short integers, to correct error due to abiguity in
19 identification of clusters.
20 BDW 06 March 1997
21 ------------------------------------------------------------------*/
22
MCW_find_clusters(int nx,int ny,int nz,float dx,float dy,float dz,int ftype,void * fim,float max_dist)23 MCW_cluster_array * MCW_find_clusters(
24 int nx, int ny, int nz,
25 float dx, float dy, float dz,
26 int ftype , void * fim ,
27 float max_dist )
28 {
29 MCW_cluster_array * clust_arr ;
30 MCW_cluster * clust , * mask ;
31 int ii,jj,kk , nxy,nxyz , ijk=0 , ijk_last , mnum ;
32 int icl , jma , ijkcl , ijkma , did_one ;
33 float fimv=0.0 ;
34 short * sfar=NULL ;
35 float * ffar=NULL ;
36 byte * bfar=NULL ;
37 short ic, jc, kc;
38 short im, jm, km;
39
40 ENTRY("MCW_find_clusters") ;
41
42 if( fim == NULL ) RETURN(NULL) ;
43
44 switch( ftype ){
45 default: RETURN(NULL) ;
46 case MRI_short: sfar = (short *) fim ; break ;
47 case MRI_byte : bfar = (byte *) fim ; break ;
48 case MRI_float: ffar = (float *) fim ; break ;
49 }
50
51 /*--- make a cluster that is a mask of points closer than max_dist ---*/
52
53 mask = MCW_build_mask (dx, dy, dz, max_dist);
54 if (mask == NULL)
55 {
56 fprintf (stderr, "Unable to build mask in MCW_find_clusters");
57 RETURN(NULL) ;
58 }
59
60 nxy = nx*ny ; nxyz = nxy * nz ;
61
62 mnum = mask->num_pt ;
63
64 /*--- scan through array, find next nonzero point, build a cluster, ... ---*/
65
66 INIT_CLARR(clust_arr) ;
67
68 ijk_last = 0 ; /* starting point for scan */
69 do {
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 points! */
93
94 #ifdef CLUST_DEBUG
95 printf(" starting cluster at ijk=%d\n",ijk) ;
96 #endif
97
98 ijk_last = ijk+1 ; /* start scanning here next time */
99
100 INIT_CLUSTER(clust) ; /* make a new (empty) cluster */
101 IJK_TO_THREE(ijk,ic,jc,kc,nx,nxy) ;
102 ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ; /* start it off! */
103
104 /*--
105 for each point in cluster:
106 check points offset by the mask for nonzero entries in fim
107 enter those into cluster
108 continue until end of cluster is reached
109 (note that cluster is expanding as we progress)
110 --*/
111
112 switch( ftype ){
113 case MRI_short: /* clust->num_pt increases with each ADDTO_CLUSTER */
114 for( icl=0 ; icl < clust->num_pt ; icl++ ){
115 ic = clust->i[icl]; /* check around this point */
116 jc = clust->j[icl];
117 kc = clust->k[icl];
118
119 for( jma=0 ; jma < mnum ; jma++ ){ /* scanning mask */
120 im = ic + mask->i[jma]; /* offsets from (ic,jc,kc) pt */
121 jm = jc + mask->j[jma];
122 km = kc + mask->k[jma];
123 if( im < 0 || im >= nx || /* outside the 3D volume? */
124 jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
125
126 ijkma = THREE_TO_IJK (im, jm, km, nx, nxy); /* 3D index */
127
128 /* is ijkma: a previous point?
129 outside the box?
130 not "active"? -- if any of these, skip it */
131
132 if( ijkma < ijk_last || ijkma >= nxyz || sfar[ijkma] == 0 ) continue ;
133
134 ADDTO_CLUSTER( clust , im, jm, km, sfar[ijkma] ) ; /* new! */
135 sfar[ijkma] = 0 ; /* mark it as used up */
136 }
137 }
138 break ;
139
140 case MRI_byte:
141 for( icl=0 ; icl < clust->num_pt ; icl++ ){
142 ic = clust->i[icl];
143 jc = clust->j[icl];
144 kc = clust->k[icl];
145
146 for( jma=0 ; jma < mnum ; jma++ ){
147 im = ic + mask->i[jma];
148 jm = jc + mask->j[jma];
149 km = kc + mask->k[jma];
150 if( im < 0 || im >= nx ||
151 jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
152
153 ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
154 if( ijkma < ijk_last || ijkma >= nxyz || bfar[ijkma] == 0 ) continue ;
155
156 ADDTO_CLUSTER( clust , im, jm, km, bfar[ijkma] ) ;
157 bfar[ijkma] = 0 ;
158 }
159 }
160 break ;
161
162 case MRI_float:
163 for( icl=0 ; icl < clust->num_pt ; icl++ ){
164 ic = clust->i[icl];
165 jc = clust->j[icl];
166 kc = clust->k[icl];
167
168 for( jma=0 ; jma < mnum ; jma++ ){
169 im = ic + mask->i[jma];
170 jm = jc + mask->j[jma];
171 km = kc + mask->k[jma];
172 if( im < 0 || im >= nx ||
173 jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
174
175 ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
176 if( ijkma < ijk_last || ijkma >= nxyz || ffar[ijkma] == 0.0 ) continue ;
177
178 ADDTO_CLUSTER( clust , im, jm, km, ffar[ijkma] ) ;
179 ffar[ijkma] = 0.0 ;
180 }
181 }
182 break ;
183 }
184
185 /* Cluster 'clust' is complete, add it to the cluster array */
186
187 ADDTO_CLARR(clust_arr,clust) ;
188 } while( 1 ) ;
189
190 /* the mask is now obsolete */
191
192 KILL_CLUSTER(mask) ;
193
194 /* Didn't find anything at all? Weird. */
195
196 if( clust_arr->num_clu <= 0 ){ DESTROY_CLARR(clust_arr) ; }
197
198 RETURN(clust_arr) ;
199 }
200
201 /*---------------------------------------------------------------
202 Write the points stored in a cluster back into a volume
203
204 Coordinates of voxels in clusters are now stored as 3 separate
205 short integers, to correct error due to abiguity in
206 identification of clusters.
207 BDW 06 March 1997
208 -----------------------------------------------------------------*/
209
MCW_cluster_to_vol(int nx,int ny,int nz,int ftype,void * fim,MCW_cluster * clust)210 void MCW_cluster_to_vol( int nx , int ny , int nz ,
211 int ftype , void *fim , MCW_cluster *clust )
212 {
213 int icl, ijk ;
214 int nxy ;
215 short *sfar ;
216 float *ffar ;
217 byte *bfar ;
218
219 ENTRY("MCW_cluster_to_vol") ;
220
221 if( fim == NULL || clust == NULL ) EXRETURN ;
222
223 nxy = nx * ny;
224
225 switch( ftype ){
226 case MRI_short:
227 sfar = (short *) fim ;
228 for( icl=0 ; icl < clust->num_pt ; icl++ )
229 {
230 ijk = THREE_TO_IJK (clust->i[icl], clust->j[icl], clust->k[icl],
231 nx, nxy);
232 sfar[ijk] = clust->mag[icl] ;
233 }
234 EXRETURN ;
235
236 case MRI_byte:
237 bfar = (byte *) fim ;
238 for( icl=0 ; icl < clust->num_pt ; icl++ )
239 {
240 ijk = THREE_TO_IJK (clust->i[icl], clust->j[icl], clust->k[icl],
241 nx, nxy);
242 bfar[ijk] = clust->mag[icl] ;
243 }
244 EXRETURN ;
245
246 case MRI_float:
247 ffar = (float *) fim ;
248 for( icl=0 ; icl < clust->num_pt ; icl++ )
249 {
250 ijk = THREE_TO_IJK (clust->i[icl], clust->j[icl], clust->k[icl],
251 nx, nxy);
252 ffar[ijk] = clust->mag[icl] ;
253 }
254 EXRETURN ;
255 }
256
257 EXRETURN ; /* should not be reached */
258 }
259
260 /*-------------------------------------------------------------------*/
261 /* Put the values from the dataset at the cluster locations into
262 the 'mag' component of the cluster, and zero out the dataset
263 at those location. Can use MCW_cluster_to_vol() to restore later.
264 This pair of functions is used in Clusterize for Flash-ing.
265 ---------------------------------------------------------------------*/
266
MCW_vol_to_cluster(int nx,int ny,int nz,int ftype,void * fim,MCW_cluster * clust)267 void MCW_vol_to_cluster( int nx , int ny , int nz ,
268 int ftype , void *fim , MCW_cluster *clust )
269 {
270 int icl, ijk ;
271 int nxy ;
272 short *sfar ;
273 float *ffar ;
274 byte *bfar ;
275
276 ENTRY("MCW_vol_to_cluster") ;
277
278 if( fim == NULL || clust == NULL ) EXRETURN ;
279
280 nxy = nx * ny;
281
282 switch( ftype ){
283 case MRI_short:
284 sfar = (short *) fim ;
285 for( icl=0 ; icl < clust->num_pt ; icl++ )
286 {
287 ijk = THREE_TO_IJK (clust->i[icl], clust->j[icl], clust->k[icl],
288 nx, nxy);
289 clust->mag[icl] = sfar[ijk] ; sfar[ijk] = 0 ;
290 }
291 EXRETURN ;
292
293 case MRI_byte:
294 bfar = (byte *) fim ;
295 for( icl=0 ; icl < clust->num_pt ; icl++ )
296 {
297 ijk = THREE_TO_IJK (clust->i[icl], clust->j[icl], clust->k[icl],
298 nx, nxy);
299 clust->mag[icl] = bfar[ijk] ; bfar[ijk] = 0 ;
300 }
301 EXRETURN ;
302
303 case MRI_float:
304 ffar = (float *) fim ;
305 for( icl=0 ; icl < clust->num_pt ; icl++ )
306 {
307 ijk = THREE_TO_IJK (clust->i[icl], clust->j[icl], clust->k[icl],
308 nx, nxy);
309 clust->mag[icl] = ffar[ijk] ; ffar[ijk] = 0.0f ;
310 }
311 EXRETURN ;
312 }
313
314 EXRETURN ; /* should not be reached */
315 }
316
317
318 /*----------------------------------------------------------------------------
319 Erosion and dilation of 3d clusters.
320
321 Purpose: Sometimes a cluster consists of several main bodies connected by
322 a narrow path. The objective is to sever the connecting path,
323 while keeping the main bodies intact, thereby producing separate
324 clusters.
325
326 Method: Erosion is applied to eliminate the outer layer of voxels. This
327 should cut the connecting path. The original main bodies are
328 then restored by dilating the result.
329
330 Author: B. Douglas Ward
331
332 Date: 16 June 1998
333
334
335 -----------------------------------------------------------------------------*/
336
MCW_erode_clusters(int nx,int ny,int nz,float dx,float dy,float dz,int ftype,void * fim,float max_dist,float pv,int dilate)337 void MCW_erode_clusters
338 (
339 int nx, int ny, int nz, /* dimensions of volume fim */
340 float dx, float dy, float dz, /* voxel dimensions */
341 int ftype, /* data type */
342 void * fim, /* volume data */
343 float max_dist, /* voxel connectivity radius */
344 float pv, /* fraction pv of voxels within max_dist must
345 be in the cluster */
346 int dilate /* boolean for perform dilation phase */
347 )
348
349 {
350 MCW_cluster * mask = NULL; /* mask determines nbhd membership */
351 int minimum; /* minimum number of voxels in nbhd */
352 int count; /* count of voxels in neighborhood */
353 int nxy, nxyz; /* numbers of voxels */
354 int ijk, iv, jv, kv; /* voxel indices */
355 int ijkm, im, jm, km; /* voxel indices */
356 int imask, nmask; /* mask indices */
357 short * sfar=NULL; /* pointer to short data */
358 byte * bfar=NULL; /* pointer to byte data */
359 float * ffar=NULL; /* pointer to float data */
360 float * efim = NULL; /* copy of eroded voxels */
361
362 ENTRY("MCW_erode_clusters") ;
363
364
365 /*----- Just in case -----*/
366 if ( fim == NULL ) EXRETURN;
367
368
369 /*----- Initialize local variables -----*/
370 nxy = nx * ny; nxyz = nxy * nz;
371
372
373 /*----- Set pointer to input data -----*/
374 switch (ftype)
375 {
376 default: EXRETURN;
377 case MRI_short: sfar = (short *) fim; break;
378 case MRI_byte : bfar = (byte *) fim; break;
379 case MRI_float: ffar = (float *) fim; break;
380 }
381
382
383 /*----- Initialization for copy of eroded voxels -----*/
384 efim = (float *) malloc (sizeof(float) * nxyz);
385 if (efim == NULL)
386 {
387 fprintf (stderr, "Unable to allocate memory in MCW_erode_clusters");
388 EXRETURN;
389 }
390 for (ijk = 0; ijk < nxyz; ijk++)
391 efim[ijk] = 0.0;
392
393
394 /*--- Make a cluster that is a mask of points closer than max_dist ---*/
395 mask = MCW_build_mask (dx, dy, dz, max_dist);
396 if (mask == NULL)
397 {
398 fprintf (stderr, "Unable to build mask in MCW_erode_clusters");
399 EXRETURN;
400 }
401
402
403 /*----- Calculate minimum number of voxels in nbhd. for non-erosion -----*/
404 nmask = mask->num_pt ;
405 minimum = floor(pv*nmask + 0.99);
406 if (minimum <= 0) EXRETURN; /*----- Nothing will be eroded -----*/
407
408
409 /*----- Step 1: Identify voxels to be eroded -----*/
410 switch (ftype)
411 {
412 case MRI_short:
413 for (ijk = 0; ijk < nxyz; ijk++)
414 {
415 if (sfar[ijk] == 0) continue;
416 IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
417
418 /*----- Count number of active voxels in the neighborhood -----*/
419 count = 0;
420 for (imask = 0; imask < nmask; imask++)
421 {
422 im = iv + mask->i[imask];
423 jm = jv + mask->j[imask];
424 km = kv + mask->k[imask];
425 if ( im < 0 || jm < 0 || km < 0 ||
426 im >= nx || jm >= ny || km >= nz ) continue;
427 ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
428 if (sfar[ijkm] != 0) count++;
429 }
430
431 /*----- Record voxel to be eroded -----*/
432 if (count < minimum) efim[ijk] = (float) sfar[ijk];
433 }
434 break;
435
436 case MRI_byte:
437 for (ijk = 0; ijk < nxyz; ijk++)
438 {
439 if (bfar[ijk] == 0) continue;
440 IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
441
442 /*----- Count number of active voxels in the neighborhood -----*/
443 count = 0;
444 for (imask = 0; imask < nmask; imask++)
445 {
446 im = iv + mask->i[imask];
447 jm = jv + mask->j[imask];
448 km = kv + mask->k[imask];
449 if ( im < 0 || jm < 0 || km < 0 ||
450 im >= nx || jm >= ny || km >= nz ) continue;
451 ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
452 if (bfar[ijkm] != 0) count++;
453 }
454
455 /*----- Record voxel to be eroded -----*/
456 if (count < minimum) efim[ijk] = (float) bfar[ijk];
457 }
458 break;
459
460 case MRI_float:
461 for (ijk = 0; ijk < nxyz; ijk++)
462 {
463 if (ffar[ijk] == 0.0) continue;
464 IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
465
466 /*----- Count number of active voxels in the neighborhood -----*/
467 count = 0;
468 for (imask = 0; imask < nmask; imask++)
469 {
470 im = iv + mask->i[imask];
471 jm = jv + mask->j[imask];
472 km = kv + mask->k[imask];
473 if ( im < 0 || jm < 0 || km < 0 ||
474 im >= nx || jm >= ny || km >= nz ) continue;
475 ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
476 if (ffar[ijkm] != 0.0) count++;
477 }
478
479 /*----- Record voxel to be eroded -----*/
480 if (count < minimum) efim[ijk] = ffar[ijk];
481 }
482 break;
483
484 }
485
486
487 /*----- Step 2: Erode voxels -----*/
488 switch (ftype)
489 {
490 case MRI_short:
491 for (ijk = 0; ijk < nxyz; ijk++)
492 if (efim[ijk] != 0.0) sfar[ijk] = 0;
493 break;
494
495 case MRI_byte:
496 for (ijk = 0; ijk < nxyz; ijk++)
497 if (efim[ijk] != 0.0) bfar[ijk] = 0;
498 break;
499
500 case MRI_float:
501 for (ijk = 0; ijk < nxyz; ijk++)
502 if (efim[ijk] != 0.0) ffar[ijk] = 0.0;
503 break;
504 }
505
506
507 /*----- Proceed with dilation phase? -----*/
508 if (dilate)
509 {
510
511
512 /*----- Step 3: Identify voxels to be restored -----*/
513 switch (ftype)
514 {
515 case MRI_short:
516 for (ijk = 0; ijk < nxyz; ijk++)
517 {
518 if (efim[ijk] == 0.0) continue;
519 IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
520
521 /*---- Determine if any active voxels in the neighborhood ----*/
522 for (imask = 0; imask < nmask; imask++)
523 {
524 im = iv + mask->i[imask];
525 jm = jv + mask->j[imask];
526 km = kv + mask->k[imask];
527 if ( im < 0 || jm < 0 || km < 0 ||
528 im >= nx || jm >= ny || km >= nz ) continue;
529 ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
530 if (sfar[ijkm] != 0) break;
531 }
532
533 /*----- Reset voxel not to be restored -----*/
534 if (imask == nmask) efim[ijk] = 0.0;
535 }
536 break;
537
538 case MRI_byte:
539 for (ijk = 0; ijk < nxyz; ijk++)
540 {
541 if (efim[ijk] == 0.0) continue;
542 IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
543
544 /*---- Determine if any active voxels in the neighborhood ----*/
545 for (imask = 0; imask < nmask; imask++)
546 {
547 im = iv + mask->i[imask];
548 jm = jv + mask->j[imask];
549 km = kv + mask->k[imask];
550 if ( im < 0 || jm < 0 || km < 0 ||
551 im >= nx || jm >= ny || km >= nz ) continue;
552 ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
553 if (bfar[ijkm] != 0) break;
554 }
555
556 /*----- Reset voxel not to be restored -----*/
557 if (imask == nmask) efim[ijk] = 0.0;
558 }
559 break;
560
561 case MRI_float:
562 for (ijk = 0; ijk < nxyz; ijk++)
563 {
564 if (efim[ijk] == 0.0) continue;
565 IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
566
567 /*---- Determine if any active voxels in the neighborhood ----*/
568 for (imask = 0; imask < nmask; imask++)
569 {
570 im = iv + mask->i[imask];
571 jm = jv + mask->j[imask];
572 km = kv + mask->k[imask];
573 if ( im < 0 || jm < 0 || km < 0 ||
574 im >= nx || jm >= ny || km >= nz ) continue;
575 ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
576 if (ffar[ijkm] != 0.0) break;
577 }
578
579 /*----- Reset voxel not to be restored -----*/
580 if (imask == nmask) efim[ijk] = 0.0;
581 }
582 break;
583 }
584
585
586 /*----- Step 4: Restore voxels -----*/
587 switch (ftype)
588 {
589 case MRI_short:
590 for (ijk = 0; ijk < nxyz; ijk++)
591 if (efim[ijk] != 0.0) sfar[ijk] = (short) efim[ijk];
592 break;
593
594 case MRI_byte:
595 for (ijk = 0; ijk < nxyz; ijk++)
596 if (efim[ijk] != 0.0) bfar[ijk] = (byte) efim[ijk];
597 break;
598
599 case MRI_float:
600 for (ijk = 0; ijk < nxyz; ijk++)
601 if (efim[ijk] != 0.0) ffar[ijk] = efim[ijk];
602 break;
603 }
604
605 } /* if (dilate) */
606
607
608 /*----- Release memory -----*/
609 KILL_CLUSTER(mask) ;
610 free (efim); efim = NULL;
611 EXRETURN ;
612 }
613