1 #include "mrilib.h"
2
3 #undef DEBUG
4
5 #undef ASSIF
6 #define ASSIF(p,v) if( p!= NULL ) *p = v
7 #define NOOP
8
9 /* local prototypes */
10 static int find_connected_set(byte *, int, int, int, int,
11 THD_ivec3 *, int_list *);
12 static int set_mask_vals (byte *, int_list *, byte);
13
14 static int dall = 1024 ;
15
16 # define DALL 1024 /* Allocation size for cluster arrays */
17
18 /*---------------------------------------------------------------------*/
19
20 static int verb = 0 ; /* 28 Oct 2003 */
THD_automask_verbose(int v)21 void THD_automask_verbose( int v ){ verb = v ; }
22
23 /*---------------------------------------------------------------------*/
24
25 static int exterior_clip = 0 ;
THD_automask_extclip(int e)26 void THD_automask_extclip( int e ){ exterior_clip = e ; }
27
28 /*---------------------------------------------------------------------*/
29
30 static float clfrac = 0.5f ; /* 20 Mar 2006 */
THD_automask_set_clipfrac(float f)31 void THD_automask_set_clipfrac( float f )
32 {
33 clfrac = (f >= 0.1f && f <= 0.9f) ? f : 0.5f ;
34 }
35
36 /*---------------------------------------------------------------------*/
37
38 /* parameters for erode/restore peeling */
39
40 static int peelcount = 1 ; /* 24 Oct 2006 */
41 static int peelthr = 17 ;
THD_automask_set_peelcounts(int p,int t)42 void THD_automask_set_peelcounts( int p , int t )
43 {
44 peelcount = (p > 0) ? p : 1 ;
45 peelthr = (t >= 6 && t <= 26) ? t : 17 ;
46 }
47
48 /*---------------------------------------------------------------------*/
49
50 static int gradualize = 1 ;
THD_automask_set_gradualize(int n)51 void THD_automask_set_gradualize( int n ){ gradualize = n; }
52
53 /*---------------------------------------------------------------------*/
54
55 static int cheapo = 0 ;
THD_automask_set_cheapo(int n)56 void THD_automask_set_cheapo( int n ){ cheapo = n; } /* 13 Aug 2007 */
57
58 /*---------------------------------------------------------------------*/
59
60 static int onlypos = 0 ;
THD_automask_set_onlypos(int n)61 void THD_automask_set_onlypos( int n ){ onlypos = n; } /* 09 Nov 2020 */
62
63 /*---------------------------------------------------------------------*/
64
mask_count(int nvox,byte * mmm)65 INLINE int mask_count( int nvox , byte *mmm )
66 {
67 register int ii , nn ;
68 if( nvox <= 0 || mmm == NULL ) return 0 ;
69 for( nn=ii=0 ; ii < nvox ; ii++ ) nn += (mmm[ii] != 0) ;
70 return nn ;
71 }
72
73 /*---------------------------------------------------------------------*/
74
mask_intersect_count(int nvox,byte * mmm,byte * nnn)75 int mask_intersect_count( int nvox , byte *mmm , byte *nnn )
76 {
77 register int nint , ii ;
78 if( nvox <= 0 || mmm == NULL || nnn == NULL ) return 0 ;
79 for( nint=ii=0 ; ii < nvox ; ii++ ) nint += (mmm[ii] && nnn[ii]) ;
80 return nint ;
81 }
82
83 /*---------------------------------------------------------------------*/
84
mask_union_count(int nvox,byte * mmm,byte * nnn)85 int mask_union_count( int nvox , byte *mmm , byte *nnn )
86 {
87 register int nint , ii ;
88 if( nvox <= 0 ) return 0 ;
89 if( mmm == NULL && nnn != NULL ) return mask_count( nvox , nnn ) ;
90 if( mmm != NULL && nnn == NULL ) return mask_count( nvox , mmm ) ;
91 for( nint=ii=0 ; ii < nvox ; ii++ ) nint += (mmm[ii] || nnn[ii]) ;
92 return nint ;
93 }
94
95 /*---------------------------------------------------------------------*/
96
mask_rgyrate(int nx,int ny,int nz,byte * mmm)97 float_triple mask_rgyrate( int nx, int ny, int nz , byte *mmm )
98 {
99 float_triple xyz={0.0f,0.0f,0.0f} ;
100 float xx,yy,zz , kq,jq , xc,yc,zc ; int ii,jj,kk , vv , nmmm ;
101
102 if( nx < 1 || ny < 1 || nz < 1 || mmm == NULL ) return xyz ;
103
104 xc = yc = zc = 0.0f ; nmmm = 0 ;
105 for( vv=kk=0 ; kk < nz ; kk++ ){
106 for( jj=0 ; jj < ny ; jj++ ){
107 for( ii=0 ; ii < nx ; ii++,vv++ ){
108 if( mmm[vv] ){ xc += ii ; yc += jj ; zc += kk ; nmmm++ ; }
109 }}}
110 if( nmmm <= 1 ) return xyz ;
111 xc /= nmmm ; yc /= nmmm ; zc /= nmmm ;
112
113 xx = yy = zz = 0.0f ;
114 for( vv=kk=0 ; kk < nz ; kk++ ){
115 kq = (kk-zc)*(kk-zc) ;
116 for( jj=0 ; jj < ny ; jj++ ){
117 jq = (jj-yc)*(jj-yc) ;
118 for( ii=0 ; ii < nx ; ii++,vv++ ){
119 if( mmm[vv] ){ xx += (ii-xc)*(ii-xc) ; yy += jq ; zz += kq ; }
120 }}}
121
122 xyz.a = xx/nmmm ; xyz.b = yy/nmmm ; xyz.c = zz/nmmm ;
123 return xyz ;
124 }
125
126
127
128
129 /*
130 Apply a mask dataset to all voxels and across all sub-bricks of a
131 dataset (migrated from the mothership, 3dAutomask.c.
132 */
thd_apply_mask(THD_3dim_dataset * dset,byte * mask,char * prefix)133 THD_3dim_dataset * thd_apply_mask( THD_3dim_dataset *dset, byte *mask,
134 char *prefix )
135 {
136 THD_3dim_dataset *out_dset;
137
138 int i,j, nbriks, nvox;
139 float *data_fptr, *out_fptr;
140 byte *data_bptr, *out_bptr, *mask_ptr;
141 short *data_iptr, *out_iptr;
142 MRI_IMARR *fim_array;
143 MRI_IMAGE *fim, *data_im, *outdata_im;
144
145 nvox = DSET_NVOX(dset);
146 nbriks = dset->dblk->nvals;
147 out_dset = EDIT_empty_copy(dset) ;
148 tross_Copy_History (dset, out_dset);
149 EDIT_dset_items( out_dset ,
150 ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
151 ADN_prefix , prefix ,
152 ADN_label1 , prefix ,
153 ADN_datum_all , DSET_BRICK_TYPE(dset,0) ,
154 ADN_none ) ;
155 /* make new Image Array */
156 INIT_IMARR(fim_array);
157 for(i=0;i<nbriks;i++) {
158 fim = mri_new_conforming( DSET_BRICK(dset,i) , DSET_BRICK_TYPE(dset,i) ) ;
159 ADDTO_IMARR(fim_array, fim);
160 }
161 out_dset->dblk->brick = fim_array; /* update pointer to data */
162
163
164 for(i=0;i<nbriks;i++) {
165 data_im = DSET_BRICK(dset, i);
166 outdata_im = DSET_BRICK(out_dset, i);
167 mask_ptr = mask;
168 switch(DSET_BRICK_TYPE(dset,i) ){
169 default:
170 return NULL;
171 case MRI_short:{
172 data_iptr = mri_data_pointer(data_im);
173 out_iptr = mri_data_pointer(outdata_im);
174 for(j=0;j<nvox;j++) {
175 if(*mask_ptr++) {
176 * out_iptr++ = *data_iptr++;
177 }
178 else {
179 *out_iptr++ = 0;
180 data_iptr++;
181 }
182 }
183 }
184 break;
185
186 case MRI_float:{
187 data_fptr = (float *) mri_data_pointer(data_im);
188 out_fptr = (float *) mri_data_pointer(outdata_im);
189 for(j=0;j<nvox;j++) {
190 if(*mask_ptr++) {
191 *out_fptr++ = *data_fptr++;
192 }
193 else {
194 *out_fptr++ = 0.0;
195 data_fptr++;
196 }
197 }
198 }
199 break;
200
201 case MRI_byte:{
202 data_bptr = (byte *) mri_data_pointer(data_im);
203 out_bptr = (byte *) mri_data_pointer(outdata_im);
204 for(j=0;j<nvox;j++) {
205 if(*mask_ptr++) {
206 *out_bptr++ = *data_bptr++;
207 }
208 else {
209 *out_bptr++ = 0;
210 data_bptr++;
211 }
212 }
213 }
214 break;
215
216 }
217
218 DSET_BRICK_FACTOR(out_dset, i) = DSET_BRICK_FACTOR(dset,i) ;
219 }
220
221 THD_load_statistics( out_dset );
222 return(out_dset);
223
224 }
225
226 /*---------------------------------------------------------------------*/
227 /*! Make a byte mask for a 3D+time dataset -- 13 Aug 2001 - RWCox.
228 - compare to thd_makemask.c
229 - 05 Mar 2003: modified to put most code into mri_automask_image().
230 -----------------------------------------------------------------------*/
231
THD_automask(THD_3dim_dataset * dset)232 byte * THD_automask( THD_3dim_dataset *dset )
233 {
234 MRI_IMAGE *medim ;
235 byte *mmm ;
236
237 ENTRY("THD_automask") ;
238
239 /* average at each voxel */
240
241 if( onlypos ) medim = THD_avepos_brick(dset) ;
242 else medim = THD_aveabs_brick(dset) ;
243 if( medim == NULL ) RETURN(NULL);
244
245 mmm = mri_automask_image( medim ) ;
246
247 mri_free(medim) ; RETURN(mmm) ;
248 }
249
250 /*---------------------------------------------------------------------*/
251 /*! Make a byte mask from the average of an array of 3D images.
252 We assume that they all have the same (nx,ny,nz) dimensions.
253 -----------------------------------------------------------------------*/
254
mri_automask_imarr(MRI_IMARR * imar)255 byte * mri_automask_imarr( MRI_IMARR *imar ) /* 18 Nov 2004 */
256 {
257 MRI_IMAGE *avim , *tim , *qim ;
258 byte *mmm ;
259 int ii , jj , nvox,nim ;
260 float fac , *avar , *qar ;
261
262 ENTRY("mri_automask_imarr") ;
263
264 if( imar == NULL || IMARR_COUNT(imar) < 1 ) RETURN(NULL) ;
265
266 nim = IMARR_COUNT(imar) ;
267 if( nim == 1 ){
268 mmm = mri_automask_image( IMARR_SUBIMAGE(imar,0) ) ;
269 RETURN(mmm) ;
270 }
271
272 avim = mri_new_conforming( IMARR_SUBIMAGE(imar,0) , MRI_float ) ;
273 avar = MRI_FLOAT_PTR(avim) ;
274 nvox = avim->nvox ;
275 for( jj=0 ; jj < nim ; jj++ ){
276 tim = IMARR_SUBIMAGE(imar,jj) ;
277 if( tim->kind != MRI_float ) qim = mri_to_float(tim) ;
278 else qim = tim ;
279 qar = MRI_FLOAT_PTR(qim) ;
280 for( ii=0 ; ii < nvox ; ii++ ) avar[ii] += qar[ii] ;
281 if( qim != tim ) mri_free(qim) ;
282 }
283 fac = 1.0f / (float)nim ;
284 for( ii=0 ; ii < nvox ; ii++ ) avar[ii] *= fac ;
285 mmm = mri_automask_image( avim ) ;
286 mri_free(avim) ;
287 RETURN(mmm) ;
288 }
289
290 /*---------------------------------------------------------------------*/
291 /*! Make a byte mask from an image (3D). Adapted from THD_automask()
292 to make it possible to do this on an image directly.
293 -----------------------------------------------------------------------*/
294
mri_automask_image(MRI_IMAGE * im)295 byte * mri_automask_image( MRI_IMAGE *im )
296 {
297 float clip_val , *mar ;
298 byte *mmm = NULL ;
299 int nvox , ii,jj , nmm , nx,ny,nz ;
300 MRI_IMAGE *medim ;
301
302 ENTRY("mri_automask_image") ;
303
304 if( im == NULL ) RETURN(NULL) ;
305
306 if( im->nz == 1 ) RETURN( mri_automask_image2D(im) ) ; /* 12 Jun 2010 */
307
308 if( im->kind != MRI_float ) medim = mri_to_float(im) ;
309 else medim = im ;
310
311 /* find clip value to excise small stuff */
312
313 clip_val = THD_cliplevel(medim,clfrac) ;
314
315 if( verb ) ININFO_message("Fixed clip level = %f\n",clip_val) ;
316
317 /* create mask of values above clip value */
318
319 nvox = medim->nvox ;
320 mar = MRI_FLOAT_PTR(medim) ;
321 mmm = (byte *) calloc( sizeof(byte), nvox ) ;
322
323 if( !gradualize ){
324 for( nmm=ii=0 ; ii < nvox ; ii++ )
325 if( mar[ii] >= clip_val ){ mmm[ii] = 1; nmm++; }
326 } else {
327 MRI_IMAGE *cim; float *car, cbot=1.e+38,ctop=-1.e+38 ;
328 cim = THD_cliplevel_gradual(medim,clfrac); car = MRI_FLOAT_PTR(cim);
329 for( nmm=ii=0 ; ii < nvox ; ii++ ){
330 if( mar[ii] >= car[ii] ){ mmm[ii] = 1; nmm++; }
331 if( car[ii] < cbot ) cbot = car[ii] ;
332 if( car[ii] > ctop ) ctop = car[ii] ;
333 }
334 if( verb ) ININFO_message("Used gradual clip level = %f .. %f",cbot,ctop) ;
335 mri_free(cim) ;
336 }
337
338 if( verb ) ININFO_message("Number voxels above clip level = %d\n",nmm) ;
339 if( im != medim && (!exterior_clip || nmm==0) ){ mri_free(medim); medim=NULL; }
340 if( nmm == 0 ){ cheapo=0; RETURN(mmm); } /* should not happen */
341
342 /*-- 6 Mar 2009: if we don't have volume data, stop here [rickr] --*/
343 if( im->nx < 2 || im->ny < 2 || im->nz < 2 ) RETURN(mmm);
344
345 /*-- 10 Apr 2002: only keep the largest connected component --*/
346
347 nx = im->nx ; ny = im->ny ; nz = im->nz ;
348 dall = (nx*ny*nz)/128 ; /* allocation delta for clustering */
349
350 THD_mask_clust( nx,ny,nz, mmm ) ;
351
352 /* 18 Apr 2002: now erode the resulting volume
353 (to break off any thinly attached pieces) */
354
355 #if 0
356 if( peelcount > 1 ){ /* 25 Oct 2006 */
357 THD_mask_erodemany( nx,ny,nz, mmm, 1 ) ;
358 THD_mask_clust( nx,ny,nz, mmm ) ;
359 THD_mask_fillin_once( nx,ny,nz , mmm , 1 ) ;
360 }
361 #endif
362
363 THD_mask_erodemany( nx,ny,nz, mmm, peelcount ) ; /* 24 Oct 2006: multiple layers */
364
365 /* now recluster it, and again keep only the largest survivor */
366
367 THD_mask_clust( nx,ny,nz, mmm ) ;
368
369 /* 19 Apr 2002: fill in small holes */
370
371 jj = ii = THD_mask_fillin_once( nx,ny,nz , mmm , 1 ) ;
372 if( ii > 0 ){
373 jj += ii = THD_mask_fillin_once( nx,ny,nz , mmm , 1 ) ;
374 if( ii > 0 ){
375 jj += ii = THD_mask_fillin_once( nx,ny,nz , mmm , 1 ) ;
376 }
377 }
378
379 if( cheapo ){
380 if( medim != im ) mri_free(medim) ; /* 13 Aug 2007 */
381 cheapo = 0 ; RETURN(mmm) ;
382 }
383
384 if( jj > 0 && verb )
385 ININFO_message("Filled %5d voxels in small holes; now have %d voxels\n",
386 jj , mask_count(nvox,mmm) ) ;
387
388 nmm = 1 ;
389 jj = rint(0.016*nx) ; nmm = MAX(nmm,jj) ;
390 jj = rint(0.016*ny) ; nmm = MAX(nmm,jj) ;
391 jj = rint(0.016*nz) ; nmm = MAX(nmm,jj) ;
392
393 if( nmm > 1 || jj > 0 ){
394 for( jj=0,ii=2 ; ii < nmm ; ii++ )
395 jj += THD_mask_fillin_once( nx,ny,nz , mmm , ii ) ;
396 jj += THD_mask_fillin_completely( nx,ny,nz, mmm , nmm ) ;
397 if( jj > 0 && verb )
398 ININFO_message("Filled %5d voxels in large holes; now have %d voxels\n",
399 jj , mask_count(nvox,mmm) ) ;
400 }
401
402 THD_mask_erodemany( nx,ny,nz, mmm, 1 ) ;
403 THD_mask_clust( nx,ny,nz, mmm ) ;
404
405 /* 28 May 2002:
406 invert the mask, then find the largest cluster of 1's;
407 this will be the outside of the brain;
408 put this back into the mask, then invert again;
409 the effect will be to fill any holes left inside the brain */
410
411 for( ii=0 ; ii < nvox ; ii++ ) mmm[ii] = !mmm[ii] ;
412
413 if( verb ) ININFO_message("Clustering non-brain voxels ...\n") ;
414 THD_mask_clust( nx,ny,nz, mmm ) ;
415
416 /* mask is now 1 for non-brain voxels;
417 if we want to clip off voxels neighboring the non-brain
418 mask AND whose values are below clip_val, do so now */
419
420 if( exterior_clip ){
421 float tclip=9999.9*clip_val ;
422 jj = THD_mask_clip_neighbors( nx,ny,nz , mmm , clip_val,tclip,mar ) ;
423 if( im != medim ) mri_free(medim) ;
424 if( jj > 0 && verb )
425 ININFO_message("Removed %d exterior voxels below clip level\n",jj);
426 } else {
427 jj = 0 ;
428 }
429
430 /* and re-invert mask to get brain voxels */
431
432 for( ii=0 ; ii < nvox ; ii++ ) mmm[ii] = !mmm[ii] ;
433 if( verb ) ININFO_message("Mask now has %d voxels\n",mask_count(nvox,mmm)) ;
434
435 if( exterior_clip && jj > 0 ){
436 THD_mask_erodemany( nx,ny,nz, mmm, 1 ) ;
437 THD_mask_clust( nx,ny,nz, mmm ) ;
438 }
439
440 RETURN(mmm) ;
441 }
442
443 /*---------------------------------------------------------------------*/
444 /*! Find voxels not in the mask, but neighboring it, that are also
445 below the clip threshold in mar[]. Add these to the mask.
446 Repeat until done. Return value is number of voxels added.
447 -----------------------------------------------------------------------*/
448
THD_mask_clip_neighbors(int nx,int ny,int nz,byte * mmm,float clip_val,float tclip,float * mar)449 int THD_mask_clip_neighbors( int nx, int ny, int nz ,
450 byte *mmm, float clip_val, float tclip, float *mar )
451 {
452 /* int jm,jp , km,kp , im,ip ; */
453 int ii,jj,kk , ntot=0,nnew , j3,k3,i3 , nxy=nx*ny ;
454
455 if( mmm == NULL || mar == NULL ) return 0 ;
456
457 do{
458 nnew = 0 ;
459 for( kk=1 ; kk < nz-1 ; kk++ ){
460 k3 = kk*nxy ;
461 for( jj=1 ; jj < ny-1 ; jj++ ){
462 j3 = k3 + jj*nx ;
463 for( ii=1 ; ii < nx-1 ; ii++ ){
464 i3 = ii+j3 ;
465 if( mmm[i3] || /* in mask */
466 (mar[i3] >= clip_val && mar[i3] <= tclip) ) continue ; /* or is OK */
467
468 /* If here, voxel IS NOT in mask, and IS below threshold.
469 If any neighbors are also in mask, then add it to mask. */
470
471 if( mmm[i3-1] || mmm[i3+1] ||
472 mmm[i3-nx] || mmm[i3+nx] ||
473 mmm[i3-nxy] || mmm[i3+nxy] ){ mmm[i3] = 1; nnew++; }
474 }}}
475 ntot += nnew ;
476 } while( nnew > 0 ) ;
477
478 return ntot ;
479 }
480
481 /*---------------------------------------------------------------------*/
482 /*! Fill in a byte mask. Filling is done by looking to each side
483 (plus/minus) of a non-filled voxel, and seeing if there is a
484 filled voxel on both sides. This looking is done parallel to
485 the x-, y-, and z-axes, out to distance nside voxels.
486 - nx,ny,nz = dimensions of mask
487 - mmm = mask itself (will be altered)
488 - nside = width of fill in look to each side
489 - Return value is number of filled in voxels
490 -----------------------------------------------------------------------*/
491
THD_mask_fillin_once(int nx,int ny,int nz,byte * mmm,int nside)492 int THD_mask_fillin_once( int nx, int ny, int nz, byte *mmm, int nside )
493 {
494 int ii,jj,kk , nsx,nsy,nsz , nxy,nxyz , iv,jv,kv,ll , nfill ;
495 byte *nnn ;
496 int nx2,nx3,nx4 , nxy2,nxy3,nxy4 ;
497
498 ENTRY("THD_mask_fillin_once") ;
499
500 if( mmm == NULL || nside <= 0 ) RETURN(0) ;
501
502 nsx = (nx-1)/2 ; if( nsx > nside ) nsx = nside ;
503 nsy = (ny-1)/2 ; if( nsy > nside ) nsy = nside ;
504 nsz = (nz-1)/2 ; if( nsz > nside ) nsz = nside ;
505
506 if( nsx == 0 && nsy == 0 && nsz == 0 ) RETURN(0) ;
507
508 #ifdef DEBUG
509 ININFO_message("THD_mask_fillin_once: nsx=%d nsy=%d nsz=%d\n",nsx,nsy,nsz);
510 #else
511 STATUSi(" :: nsx",nsx); STATUSi(" :: nsy",nsy); STATUSi(" :: nsz",nsz);
512 #endif
513
514 nxy = nx*ny ; nxyz = nxy*nz ; nfill = 0 ;
515
516 nx2 = 2*nx ; nx3 = 3*nx ; nx4 = 4*nx ;
517 nxy2 = 2*nxy ; nxy3 = 3*nxy ; nxy4 = 4*nxy ;
518
519 nnn = AFMALL(byte, nxyz) ; /* stores filled in values */
520
521 /* loop over voxels */
522
523 #define FILLVOX \
524 do{ nnn[iv] = 1; nfill++; goto NextVox; } while(0)
525
526 for( kk=nsz ; kk < nz-nsz ; kk++ ){
527 kv = kk*nxy ;
528 for( jj=nsy ; jj < ny-nsy ; jj++ ){
529 jv = jj*nx + kv ;
530 for( ii=nsx ; ii < nx-nsx ; ii++ ){
531 iv = ii+jv ;
532 if( mmm[iv] ) continue ; /* already filled */
533
534 /* check in +x direction, then -x if +x hits */
535
536 switch( nsx ){
537 case 0: break ;
538 case 1:
539 if( mmm[iv+1] && mmm[iv-1] ) FILLVOX;
540 break ;
541
542 case 2:
543 if( (mmm[iv+1]||mmm[iv+2]) &&
544 (mmm[iv-1]||mmm[iv-2]) ) FILLVOX;
545 break ;
546
547 case 3:
548 if( (mmm[iv+1]||mmm[iv+2]||mmm[iv+3]) &&
549 (mmm[iv-1]||mmm[iv-2]||mmm[iv-3]) ) FILLVOX;
550 break ;
551
552 /*** This stuff is for eventual expansion to hard-coding larger steps ***/
553 #if 0
554 #define XPLU4(q) mmm[q+1]||mmm[q+2]||mmm[q+3]||mmm[q+4]
555 #define XMIN4(q) mmm[q-1]||mmm[q-2]||mmm[q-3]||mmm[q-4]
556
557 #define XPLU5(q) XPLU4(q) ||mmm[q+5]
558 #define XPLU6(q) XPLU5(q) ||mmm[q+6]
559 #define XPLU7(q) XPLU6(q) ||mmm[q+7]
560 #define XPLU8(q) XPLU7(q) ||mmm[q+8]
561 #define XPLU9(q) XPLU8(q) ||mmm[q+9]
562 #define XPLU10(q) XPLU9(q) ||mmm[q+10]
563 #define XPLU11(q) XPLU10(q)||mmm[q+11]
564
565 #define XMIN5(q) XMIN4(q) ||mmm[q-5]
566 #define XMIN6(q) XMIN5(q) ||mmm[q-6]
567 #define XMIN7(q) XMIN6(q) ||mmm[q-7]
568 #define XMIN8(q) XMIN7(q) ||mmm[q-8]
569 #define XMIN9(q) XMIN8(q) ||mmm[q-9]
570 #define XMIN10(q) XMIN9(q) ||mmm[q-10]
571 #define XMIN11(q) XMIN10(q)||mmm[q-11]
572 #endif
573
574 case 4:
575 if( (mmm[iv+1]||mmm[iv+2]||mmm[iv+3]||mmm[iv+4]) &&
576 (mmm[iv-1]||mmm[iv-2]||mmm[iv-3]||mmm[iv-4]) ) FILLVOX;
577 break ;
578
579 default:
580 for( ll=1 ; ll <= nsx ; ll++ ) if( mmm[iv+ll] ) break ;
581 if( ll <= nsx ){
582 for( ll=1 ; ll <= nsx ; ll++ ) if( mmm[iv-ll] ) break ;
583 if( ll <= nsx ) FILLVOX;
584 }
585 break ;
586 }
587
588 /* check in +y direction, then -y if +y hits */
589
590 switch( nsy ){
591 case 0: break ;
592 case 1:
593 if( mmm[iv+nx] && mmm[iv-nx] ) FILLVOX;
594 break ;
595
596 case 2:
597 if( (mmm[iv+nx]||mmm[iv+nx2]) &&
598 (mmm[iv-nx]||mmm[iv-nx2]) ) FILLVOX;
599 break ;
600
601 case 3:
602 if( (mmm[iv+nx]||mmm[iv+nx2]||mmm[iv+nx3]) &&
603 (mmm[iv-nx]||mmm[iv-nx2]||mmm[iv-nx3]) ) FILLVOX;
604 break ;
605
606 case 4:
607 if( (mmm[iv+nx]||mmm[iv+nx2]||mmm[iv+nx3]||mmm[iv+nx4]) &&
608 (mmm[iv-nx]||mmm[iv-nx2]||mmm[iv-nx3]||mmm[iv-nx4]) ) FILLVOX;
609 break ;
610
611 default:
612 if( (mmm[iv+nx]||mmm[iv+nx2]||mmm[iv+nx3]||mmm[iv+nx4]) == 0 ){
613 for( ll=5 ; ll <= nsy ; ll++ ) if( mmm[iv+ll*nx] ) break ;
614 } else {
615 ll = 1 ;
616 }
617 if( ll <= nsy ){
618 if( (mmm[iv-nx]||mmm[iv-nx2]||mmm[iv-nx3]||mmm[iv-nx4]) == 0 ){
619 for( ll=5 ; ll <= nsy ; ll++ ) if( mmm[iv-ll*nx] ) break ;
620 } else {
621 ll = 1 ;
622 }
623 if( ll <= nsy ) FILLVOX;
624 }
625 break ;
626 }
627
628 /* check in +z direction, then -z if +z hits */
629
630 switch( nsz ){
631 case 0: break ;
632 case 1:
633 if( mmm[iv+nxy] && mmm[iv-nxy] ) FILLVOX;
634 break ;
635
636 case 2:
637 if( (mmm[iv+nxy]||mmm[iv+nxy2]) &&
638 (mmm[iv-nxy]||mmm[iv-nxy2]) ) FILLVOX;
639 break ;
640
641 case 3:
642 if( (mmm[iv+nxy]||mmm[iv+nxy2]||mmm[iv+nxy3]) &&
643 (mmm[iv-nxy]||mmm[iv-nxy2]||mmm[iv-nxy3]) ) FILLVOX;
644 break ;
645
646 case 4:
647 if( (mmm[iv+nxy]||mmm[iv+nxy2]||mmm[iv+nxy3]||mmm[iv+nxy4]) &&
648 (mmm[iv-nxy]||mmm[iv-nxy2]||mmm[iv-nxy3]||mmm[iv-nxy4]) ) FILLVOX;
649 break ;
650
651 default:
652 if( (mmm[iv+nxy]||mmm[iv+nxy2]||mmm[iv+nxy3]||mmm[iv+nxy4]) == 0 ){
653 for( ll=5 ; ll <= nsz ; ll++ ) if( mmm[iv+ll*nxy] ) break ;
654 } else {
655 ll = 1 ;
656 }
657 if( ll <= nsz ){
658 if( (mmm[iv-nxy]||mmm[iv-nxy2]||mmm[iv-nxy3]||mmm[iv-nxy4]) == 0 ){
659 for( ll=5 ; ll <= nsz ; ll++ ) if( mmm[iv-ll*nxy] ) break ;
660 } else {
661 ll = 1 ;
662 }
663 if( ll <= nsz ) FILLVOX;
664 }
665 break ;
666 }
667
668 NextVox: ; /* end of loop over ii */
669 } } }
670
671 /* copy fills back into mmm */
672
673 if( nfill > 0 ){
674 for( iv=0 ; iv < nxyz ; iv++ ) if( nnn[iv] ) mmm[iv] = 1 ;
675 }
676
677 #ifdef DEBUG
678 ININFO_message("THD_mask_fillin_once: nfill=%d\n",nfill) ;
679 #else
680 STATUSi(" >> nfill",nfill) ;
681 #endif
682
683 free(nnn) ; RETURN(nfill) ;
684 }
685
686 /*----------------------------------------------------------------------*/
687 /*! Fill in a byte mask, repeatedly until it doesn't fill any more.
688 Return value is number of voxels filled.
689 ------------------------------------------------------------------------*/
690
THD_mask_fillin_completely(int nx,int ny,int nz,byte * mmm,int nside)691 int THD_mask_fillin_completely( int nx, int ny, int nz, byte *mmm, int nside )
692 {
693 int nfill=0 , kfill ;
694
695 ENTRY("THD_mask_fillin_completely") ;
696
697 do{
698 kfill = THD_mask_fillin_once( nx,ny,nz , mmm , nside ) ;
699 nfill += kfill ;
700 } while( kfill > 0 ) ;
701
702 RETURN(nfill) ;
703 }
704
705 /*----------------------------------------------------------------------*/
706 /*! Fill any holes in a byte mask. 26 Apr 2012 [rickr]
707
708 A hole is defined as a connected set of zero voxels that does not
709 reach a volume edge.
710
711 method: use less than two volumes of extra memory to be fast
712 (one volume is to make 0/1 copy of first, so values are known)
713 make 0/1 copy of volume
714 for each voxel
715 if set, continue
716 find connected set of zero voxels, setting to 2
717 if this does not contain an edge voxel, reset to 1
718 clear all edge voxels (set any 2 to 0)
719
720 if dirs, fill holes restricted to those axis directions
721
722 return -1 on error, else number of voxels filled
723 ------------------------------------------------------------------------*/
724
THD_mask_fill_holes(int nx,int ny,int nz,byte * mmm,THD_ivec3 * dirs,int verb)725 int THD_mask_fill_holes( int nx, int ny, int nz, byte *mmm,
726 THD_ivec3 * dirs, int verb )
727 {
728 int_list Cset; /* current connected set of voxels */
729 byte * bmask=NULL;
730 int nholes=0, nfilled=0, voxel;
731 int nvox, has_edge;
732
733 // rcr - pass ivec3 check_dir (data axis directions)
734
735 ENTRY("THD_mask_fill_holes");
736
737 nvox = nx*ny*nz;
738
739 if( verb > 1 ) INFO_message("filling holes");
740
741 /* make a 0/1 copy that we can edit along the way */
742 bmask = (byte *)malloc(nvox*sizeof(byte));
743 if( !bmask ) {
744 ERROR_message("TFMH: alloc failure of %d bytes", nvox);
745 RETURN(-1);
746 }
747 for( voxel = 0; voxel < nvox; voxel++ )
748 if( mmm[voxel] ) bmask[voxel] = 1;
749 else bmask[voxel] = 0;
750
751 /* initialize empty connected set list */
752 if( init_int_list(&Cset, 100) != 100 ) {
753 ERROR_message("TFMH: failed init_int_list of size %d", 100);
754 RETURN(-1);
755 }
756
757 /* start at any unset voxel and fill as either a hole or an edge set */
758 for( voxel = 0; voxel < nvox; voxel++ ) {
759 if( bmask[voxel] ) continue; /* skip any set voxel */
760
761 /* current Cset of voxels will be initialized to 2 */
762 has_edge = find_connected_set(bmask, nx, ny, nz, voxel, dirs, &Cset);
763 if( has_edge < 0 ) RETURN(-1);
764 if( Cset.num < 1 ){ ERROR_message("TFMH: FCS failure"); RETURN(-1); }
765
766 if( ! has_edge ) {
767 /* then this was a hole, count and reset new values to 1 */
768 if(verb>1) INFO_message("found hole of %d voxels", Cset.num);
769 nholes++; /* track filled interior holes */
770 nfilled += Cset.num;
771 set_mask_vals(bmask, &Cset, 1); /* modify these 2's to 1's */
772 } else if (verb>1) INFO_message("found edge set of %d voxels", Cset.num);
773 }
774
775 /* now all voxels should be set, apply bmask to mmm :
776 * set all values where bmask == 1 (ignore bmask == 2) */
777 for( voxel = 0; voxel < nvox; voxel++ )
778 if( bmask[voxel] == 1 && !mmm[voxel] ) mmm[voxel] = 1;
779
780 if(verb) INFO_message("filled %d holes (%d voxels)", nholes, nfilled);
781
782 RETURN(nfilled);
783 }
784
785
786 /*------------------------------------------------------------------*/
787 /*! set list of mask voxels to the given value */
set_mask_vals(byte * bmask,int_list * L,byte val)788 static int set_mask_vals(byte * bmask, int_list * L, byte val)
789 {
790 int ind;
791
792 ENTRY("set_mask_vals");
793 for( ind = 0; ind < L->num; ind++ )
794 bmask[L->list[ind]] = val;
795 RETURN(0);
796 }
797
798
799 /*------------------------------------------------------------------*/
800 /*! find connected set of zero voxels, and note whether this is an
801 * edge set (whether any voxel is at a volume edge)
802
803 modify new set of bmask values to 2
804
805 dirs is a restriction of connected sets - only search in those dirs
806 (basically to allow for planar hole filling)
807
808 method:
809 set voxel
810 init as search list
811 while search list is not empty
812 init new search list
813 for each voxel in list
814 for each zero neighbor (6, 18 or 26? start with 6)
815 (restrict neighbors to axis dirs)
816 set and add to new list
817 current list = new list
818
819 return whether there was an edge voxel, or -1 on error
820 *------------------------------------------------------------------*/
find_connected_set(byte * bmask,int nx,int ny,int nz,int voxel,THD_ivec3 * dirs,int_list * IL)821 static int find_connected_set(byte * bmask, int nx, int ny, int nz, int voxel,
822 THD_ivec3 * dirs, int_list * IL)
823 {
824 int_list SL1, SL2; /* search lists */
825 int_list * spa, *spb, *spt; /* changing pointers to search lists */
826 int index, nxy = nx*ny, newvox, testvox;
827 int has_edge = 0; /* is there an edge voxel */
828 int newval = 2; /* easy to change */
829 int i, j, k;
830 int doi=1, doj=1, dok=1; /* directions to search, from dirs */
831
832 ENTRY("find_connected_set");
833
834 if( !bmask || nx<0 || ny<0 || nz<0 || voxel<0 || !IL ) {
835 ERROR_message("FCS: invalid inputs"); RETURN(-1);
836 }
837
838 IL->num = 0; /* set current list as empty */
839
840 if( bmask[voxel] ) RETURN(0); /* nothing to do */
841
842 bmask[voxel] = newval; /* set current voxel and init search lists */
843 if( init_int_list(&SL1, 1000) != 1000 ) RETURN(-1);
844 if( init_int_list(&SL2, 1000) != 1000 ) RETURN(-1);
845
846 /* for a little more speed, do not move list contents around */
847 spa = &SL1; spb = &SL2;
848
849 /* using 1000 as the list memory increment size */
850 if( add_to_int_list(spa, voxel, 1000) < 0 ) RETURN(-1);
851
852 /* note directions to apply */
853 if( dirs ) {
854 doi = dirs->ijk[0];
855 doj = dirs->ijk[1];
856 dok = dirs->ijk[2];
857 }
858
859 while( spa->num > 0 ) {
860 spb->num = 0; /* init as empty */
861 extend_int_list(IL, spa); /* track all connected voxels */
862
863 /* for each voxel in list, add any unset neighbors to new list */
864 for( index = 0; index < spa->num; index++ ) {
865 newvox = spa->list[index];
866 IJK_TO_THREE(newvox, i,j,k, nx,nxy);
867
868 /* check whether this is an edge voxel */
869 if( ! has_edge ) {
870 if( dirs ) {
871 if ( doi && ( i==0 || i == nx-1 ) ) has_edge = 1;
872 else if( doj && ( j==0 || j == ny-1 ) ) has_edge = 1;
873 else if( dok && ( k==0 || k == nz-1 ) ) has_edge = 1;
874 } else if( i==0 || j==0 || k==0 || i == nx-1 || j==ny-1 || k==nz-1 )
875 has_edge = 1;
876 }
877
878 /* in each of 6 directions:
879 * if neighbor is not off edge and is not set
880 * set and add to next search list
881 */
882 if( doi ) {
883 testvox = newvox-1;
884 if( i > 0 && !bmask[testvox] ) { /* should never happen */
885 bmask[testvox] = newval;
886 add_to_int_list(spb, testvox, 1000);
887 }
888 testvox = newvox+1;
889 if( i < nx-1 && !bmask[testvox] ) {
890 bmask[testvox] = newval;
891 add_to_int_list(spb, testvox, 1000);
892 }
893 }
894
895 if( doj ) {
896 testvox = newvox-nx;
897 if( j > 0 && !bmask[testvox] ) {
898 bmask[testvox] = newval;
899 add_to_int_list(spb, testvox, 1000);
900 }
901 testvox = newvox+nx;
902 if( j < ny-1 && !bmask[testvox] ) {
903 bmask[testvox] = newval;
904 add_to_int_list(spb, testvox, 1000);
905 }
906 }
907
908 if( dok ) {
909 testvox = newvox-nxy;
910 if( k > 0 && !bmask[testvox] ) {
911 bmask[testvox] = newval;
912 add_to_int_list(spb, testvox, 1000);
913 }
914 testvox = newvox+nxy;
915 if( k < nz-1 && !bmask[testvox] ) {
916 bmask[testvox] = newval;
917 add_to_int_list(spb, testvox, 1000);
918 }
919 }
920 }
921
922 /* now switch search list pointers */
923 spt = spa; spa = spb; spb = spt;
924 }
925
926 RETURN(has_edge);
927 }
928
929
930 /*------------------------------------------------------------------*/
931
932 /*! Put (i,j,k) into the current cluster, if it is nonzero. */
933
934 # undef CPUT
935 # define CPUT(i,j,k) \
936 do{ ijk = THREE_TO_IJK(i,j,k,nx,nxy) ; \
937 if( mmm[ijk] ){ \
938 if( nnow == nall ){ /* increase array lengths */ \
939 nall += dall + nall/4 ; \
940 inow = (short *) realloc(inow,sizeof(short)*nall) ; \
941 jnow = (short *) realloc(jnow,sizeof(short)*nall) ; \
942 know = (short *) realloc(know,sizeof(short)*nall) ; \
943 } \
944 inow[nnow] = i ; jnow[nnow] = j ; know[nnow] = k ; \
945 nnow++ ; mmm[ijk] = 0 ; \
946 } } while(0)
947
948 /*------------------------------------------------------------------*/
949 /*! Find the biggest cluster of nonzeros in the byte mask mmm.
950 --------------------------------------------------------------------*/
951
THD_mask_clust(int nx,int ny,int nz,byte * mmm)952 void THD_mask_clust( int nx, int ny, int nz, byte *mmm )
953 {
954 int ii,jj,kk, icl , nxy,nxyz , ijk , ijk_last ;
955 int ip,jp,kp , im,jm,km ;
956 int nbest ; short *ibest, *jbest , *kbest ;
957 int nnow ; short *inow , *jnow , *know ; int nall ;
958
959 ENTRY("THD_mask_clust") ;
960
961 if( mmm == NULL ) EXRETURN ;
962
963 nxy = nx*ny ; nxyz = nxy * nz ;
964
965 nbest = 0 ;
966 ibest = AFMALL(short, sizeof(short)) ;
967 jbest = AFMALL(short, sizeof(short)) ;
968 kbest = AFMALL(short, sizeof(short)) ;
969
970 /*--- scan through array, find nonzero point, build a cluster, ... ---*/
971 if(verb) ININFO_message("Clustering voxels ...");
972
973 ijk_last = 0 ; if( dall < DALL ) dall = DALL ;
974 while(1) {
975 /* find next nonzero point */
976
977 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( mmm[ijk] ) break ;
978 if( ijk == nxyz ) break ; /* didn't find any! */
979
980 ijk_last = ijk+1 ; /* start here next time */
981
982 /* init current cluster list with this point */
983
984 mmm[ijk] = 0 ; /* clear found point */
985 nall = DALL ; /* # allocated pts */
986 nnow = 1 ; /* # pts in cluster */
987 inow = (short *) malloc(sizeof(short)*DALL) ; /* coords of pts */
988 jnow = (short *) malloc(sizeof(short)*DALL) ;
989 know = (short *) malloc(sizeof(short)*DALL) ;
990 IJK_TO_THREE(ijk, inow[0],jnow[0],know[0] , nx,nxy) ;
991
992 /*--
993 for each point in cluster:
994 check neighboring points for nonzero entries in mmm
995 enter those into cluster (and clear them in mmm)
996 continue until end of cluster is reached
997 (note that cluster is expanding as we progress)
998 --*/
999
1000 for( icl=0 ; icl < nnow ; icl++ ){
1001 ii = inow[icl] ; jj = jnow[icl] ; kk = know[icl] ;
1002 im = ii-1 ; jm = jj-1 ; km = kk-1 ;
1003 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ;
1004
1005 if( im >= 0 ) CPUT(im,jj,kk) ;
1006 if( ip < nx ) CPUT(ip,jj,kk) ;
1007 if( jm >= 0 ) CPUT(ii,jm,kk) ;
1008 if( jp < ny ) CPUT(ii,jp,kk) ;
1009 if( km >= 0 ) CPUT(ii,jj,km) ;
1010 if( kp < nz ) CPUT(ii,jj,kp) ;
1011 }
1012
1013 /* see if now cluster is larger than best yet */
1014
1015 if( nnow > nbest ){ /* now is bigger; */
1016 free(ibest) ; free(jbest) ; free(kbest) ; /* replace best */
1017 nbest = nnow ; ibest = inow ; /* with now */
1018 jbest = jnow ; kbest = know ;
1019 } else { /* old is bigger */
1020 free(inow) ; free(jnow) ; free(know) ; /* toss now */
1021 }
1022
1023 } /* loop ends when all nonzero points are clustered */
1024
1025 /* put 1's back in at all points in best cluster */
1026 for( icl=0 ; icl < nbest ; icl++ ){
1027 ijk = THREE_TO_IJK(ibest[icl],jbest[icl],kbest[icl],nx,nxy) ;
1028 mmm[ijk] = 1 ;
1029 }
1030 free(ibest) ; free(jbest) ; free(kbest) ;
1031
1032 if( verb ) ININFO_message("Largest cluster has %d voxels\n",nbest) ;
1033
1034 EXRETURN ;
1035 }
1036
1037 /*--------------------------------------------------------------------------*/
1038 /*! Erode away nonzero voxels that aren't neighbored by mostly other
1039 nonzero voxels. Then (maybe) restore those that were eroded that
1040 are neighbors of survivors.
1041 The neighbors are the 6,18,26 voxels closest
1042 in 3D (facing, edge and corner neighbors) by specifying NN=1,2 or 3.
1043 The byte array passed here is returned with voxel elements zeroed
1044 ----------------------------------------------------------------------------*/
THD_mask_erode(int nx,int ny,int nz,byte * mmm,int redilate,byte NN)1045 void THD_mask_erode( int nx, int ny, int nz, byte *mmm, int redilate, byte NN )
1046 {
1047 int ii,jj,kk , jy,kz, im,jm,km , ip,jp,kp , num ;
1048 int victim ; /* do we apply operation to current voxel? */
1049 int nxy=nx*ny , nxyz=nxy*nz ;
1050 int jmkm,jykm,jpkm, jmkz,jykz,jpkz, jmkp,jykp,jpkp;
1051 byte *nnn ;
1052 int nxm1, nym1, nzm1;
1053 ENTRY("THD_mask_erode") ;
1054
1055 if( mmm == NULL ) EXRETURN ;
1056
1057 nnn = (byte *)calloc(sizeof(byte),nxyz) ; /* mask of eroded voxels */
1058 if( nnn == NULL ) EXRETURN ; /* WTF? */
1059
1060 /* mark interior voxels that don't have 17 out of 18 nonzero nbhrs */
1061 STATUS("marking to erode") ;
1062 for( kk=0 ; kk < nz ; kk++ ){
1063 kz = kk*nxy ; km = kz-nxy ; kp = kz+nxy ;
1064 if( kk == 0 ) km = kz ;
1065 if( kk == nz-1 ) kp = kz ;
1066
1067 for( jj=0 ; jj < ny ; jj++ ){
1068 jy = jj*nx ; jm = jy-nx ; jp = jy+nx ;
1069 if( jj == 0 ) jm = jy ;
1070 if( jj == ny-1 ) jp = jy ;
1071 jykz = jy+kz;
1072
1073 for( ii=0 ; ii < nx ; ii++ ){
1074 if( mmm[ii+jykz] ){ /* count nonzero nbhrs */
1075 im = ii-1 ; ip = ii+1 ;
1076 if(( ii == 0 )|| (jj == 0) || (kk==0) ||
1077 ( ii == nx-1 ) || ( jj == ny-1 ) || ( kk == nz-1 ))
1078 nnn[ii+jykz] = 1;
1079 else {
1080
1081 // count neighbors
1082 // note slices are semi-graphically shown
1083 // the text order seems to indicate ijk as row,col,slice
1084 // but they're really col,row,slice. It doesn't really matter
1085 // because all combinations are used similarly
1086
1087 /* Complicated logic is confusing both compilers (warnings) */
1088 /* and coders, so simplify. Use 'victim' to decide whether */
1089 /* to apply the current operation. Then set nnn later. */
1090 /* First check NN1 voxels, then if needed, NN2 and NN3. */
1091 /* 16 Mar 2021 [rickr] */
1092 victim = 0 ;
1093 num = mmm[ii+jy+km]
1094 + mmm[im+jy+kz]
1095 + mmm[ii+jm+kz] + mmm[ii+jp+kz]
1096 + mmm[ip+jykz]
1097 + mmm[ii+jy+kp];
1098 if(num<6) victim = 1 ; /* mark to erode */
1099
1100 /* check NN2 cases (18 neighbors), if not already a victim */
1101 if(NN>=2 && !victim) {
1102 num += mmm[im+jy+km]
1103 + mmm[ii+jm+km] + mmm[ii+jp+km]
1104 + mmm[ip+jy+km]
1105 + mmm[im+jm+kz] + mmm[im+jp+kz]
1106
1107 + mmm[ip+jm+kz] + mmm[ip+jp+kz]
1108 + mmm[im+jy+kp]
1109 + mmm[ii+jm+kp] + mmm[ii+jp+kp]
1110 + mmm[ip+jy+kp] ;
1111 // if not enough neighbors, erode
1112 if( num < 18 ) victim = 1 ; /* mark to erode */
1113 }
1114
1115 /* check NN3 cases (26 neighbors), if not already a victim */
1116 if(NN==3 && !victim) {
1117 num += mmm[im+jm+km] + mmm[im+jp+km]
1118 + mmm[ip+jm+km] + mmm[ip+jp+km]
1119 + mmm[im+jm+kp] + mmm[im+jp+kp]
1120 + mmm[ip+jm+kp] + mmm[ip+jp+kp];
1121
1122 if( num < 26 ) victim = 1 ; /* mark to erode */
1123 }
1124
1125 /* mark victims for erosion */
1126 nnn[ii+jykz] = victim ;
1127
1128 } // end non-edge box (ii,jj,kk=0,...)
1129 } // end if in mask
1130 } } }
1131
1132 STATUS("eroding") ;
1133 for( jj=ii=0 ; ii < nxyz ; ii++ ) /* actually erode */
1134 if( nnn[ii] ){ mmm[ii] = 0 ; jj++ ; }
1135
1136 if( verb && jj > 0 ) ININFO_message("Eroded %d voxels\n",jj) ;
1137
1138 /* optionally re-dilate eroded voxels that are next to survivors */
1139 /* useful for cleaning up 1 voxel deep necks that are eroded in at least
1140 two directions */
1141 if(redilate) {
1142 STATUS("marking to redilate") ;
1143 for( kk=0 ; kk < nz ; kk++ ){
1144 kz = kk*nxy ; km = kz-nxy ; kp = kz+nxy ;
1145 if( kk == 0 ) km = kz ;
1146 if( kk == nz-1 ) kp = kz ;
1147
1148 for( jj=0 ; jj < ny ; jj++ ){
1149 jy = jj*nx ; jm = jy-nx ; jp = jy+nx ;
1150 if( jj == 0 ) jm = jy ;
1151 if( jj == ny-1 ) jp = jy ;
1152
1153 for( ii=0 ; ii < nx ; ii++ ){
1154 if( nnn[ii+jy+kz] ){ /* was eroded */
1155
1156 im = ii-1 ; ip = ii+1 ;
1157 // borders make more sense here as placeholders than for
1158 // erosion case
1159 if( ii == 0 ) im = 0 ;
1160 if( ii == nx-1 ) ip = ii ;
1161
1162 /* Complicated logic is confusing both compilers (warnings) */
1163 /* and coders, so simplify. */
1164 /* Set victim if any NN1 voxel is set (else go to NN2,3). */
1165 /* (dglen indentation provides a "visual" of neighbors) */
1166 /* (nnn[] will later be set to victim) 16 Mar 2021 [rickr] */
1167
1168 /* NN1 case, check for set first neighbors */
1169 victim = mmm[ii+jy+km]
1170 || mmm[im+jy+kz]
1171 || mmm[ii+jm+kz] || mmm[ii+jp+kz]
1172 || mmm[ip+jy+kz]
1173 || mmm[ii+jy+kp];
1174
1175 /* similar NN2 case */
1176 if(NN>=2 && !victim)
1177 victim = mmm[im+jy+km]
1178 || mmm[ii+jm+km] || mmm[ii+jp+km]
1179 || mmm[ip+jy+km]
1180 || mmm[im+jm+kz] || mmm[im+jp+kz]
1181
1182 || mmm[ip+jm+kz] || mmm[ip+jp+kz]
1183 || mmm[im+jy+kp]
1184 || mmm[ii+jm+kp] || mmm[ii+jp+kp]
1185 || mmm[ip+jy+kp];
1186
1187 /* similar NN3 case */
1188 if(NN>=3 && !victim)
1189 victim = mmm[im+jm+km] || mmm[im+jp+km]
1190 || mmm[ip+jm+km] || mmm[ip+jp+km]
1191 || mmm[im+jm+kp] || mmm[im+jp+kp]
1192 || mmm[ip+jm+kp] || mmm[ip+jp+kp];
1193
1194 /* if we found a victim (neighbor), mark for redilation */
1195 nnn[ii+jy+kz] = victim;
1196
1197 } // end if in mask
1198 }}} /* end of ii,jj,kk loops */
1199
1200 /* actually do the dilation */
1201
1202 STATUS("redilating") ;
1203 for( jj=ii=0 ; ii < nxyz ; ii++ )
1204 if( nnn[ii] ){ mmm[ii] = 1 ; jj++ ; }
1205
1206 if( verb && jj > 0 ) ININFO_message("Restored %d eroded voxels\n",jj) ;
1207
1208 } /* end of redilate */
1209
1210 free(nnn) ; EXRETURN ;
1211 }
1212
1213
1214 /*--------------------------------------------------------------------------*/
1215 /*! Erode away nonzero voxels that aren't neighbored by mostly other nonzero
1216 voxels.
1217
1218 Pass nerode to match ndil from THD_mask_dilate.
1219 For each masked voxel, count unmasked neighbors. If there are
1220 at least nerode, then erode the voxel.
1221
1222 Pass NN to specify NN=1,2,3, again to match THD_mask_dilate.
1223 The nerode count will be restricted by the given NN level.
1224 19 May 2020 [rickr]
1225
1226 This allows for symmetric dilate/erode operations.
1227 It is symmetric with THD_mask_dilate. 7 May 2012 [rickr]
1228 ----------------------------------------------------------------------------*/
1229
THD_mask_erode_sym(int nx,int ny,int nz,byte * mmm,int nerode,int NN)1230 void THD_mask_erode_sym( int nx, int ny, int nz, byte *mmm, int nerode,
1231 int NN )
1232 {
1233 /* use ix as well, for consistency at a tiny cost of speed */
1234 int ii,jj,kk , ix,jy,kz, im,jm,km , ip,jp,kp , num ;
1235 int nxy=nx*ny , nxyz=nxy*nz, nmask ;
1236 int nmax; /* depends on NN */
1237 byte *nnn ;
1238
1239 ENTRY("THD_mask_erode_sym") ;
1240
1241 if( mmm == NULL ) EXRETURN ;
1242
1243 /* set nmax to be the full neighbor count at the given NN level */
1244 if( NN >= 3 ) nmax = 26; /* 3x3x3 box minus center */
1245 else if( NN == 2 ) nmax = 18; /* minus the 8 corners */
1246 else nmax = 6; /* count only face neighbors */
1247
1248 /* for clarity, first restrict zero neighbor count to what is appropriate */
1249 if ( nerode < 1 ) nerode = 1 ;
1250 else if( nerode > nmax ) nerode = nmax ;
1251
1252 /* then switch logic from (>= nerode zero neighbors)
1253 * to (<= (nmax-nerode) set ones) */
1254 nmask = nmax-nerode;
1255
1256 nnn = (byte *)calloc(sizeof(byte),nxyz) ; /* mask of eroded voxels */
1257 if( ! nnn ) EXRETURN ; /* death to Ming! */
1258
1259 /* mark for erosion voxels that do not have sufficient nonzero nhbrs */
1260 STATUS("marking to erode") ;
1261
1262 /* check for edge voxels first, use planar duplication */
1263 for( kk=0 ; kk < nz ; kk++ ){
1264 kz = kk*nxy ; km = kz-nxy ; kp = kz+nxy ;
1265 if( kk == 0 ) km = kz ;
1266 if( kk == nz-1 ) kp = kz ;
1267
1268 for( jj=0 ; jj < ny ; jj++ ){
1269 jy = jj*nx ; jm = jy-nx ; jp = jy+nx ;
1270 if( jj == 0 ) jm = jy ;
1271 if( jj == ny-1 ) jp = jy ;
1272
1273 for( ii=0 ; ii < nx ; ii++ ){
1274 /* if this voxel is not set, skip it */
1275 if( ! mmm[ii+jy+kz] ) continue;
1276
1277 ix = ii ; im = ii-1 ; ip = ii+1 ;
1278 if( ii == 0 ) im = 0 ;
1279 if( ii == nx-1 ) ip = ii ;
1280
1281 /* now count nonzero nbhrs, one NN level at a time */
1282
1283 /* NN1, where exactly one axis is not at the center */
1284 /* (and that can either plus or minus) */
1285 /* mod: i j k */
1286 num = mmm[im+jy+kz] + mmm[ix+jm+kz] + mmm[ix+jy+km]
1287 + mmm[ip+jy+kz] + mmm[ix+jp+kz] + mmm[ix+jy+kp];
1288
1289 /* NN2, where exactly two axes are not at the center */
1290 if( NN >= 2 )
1291 /* fix: i j k */
1292 num += mmm[ix+jm+km] + mmm[im+jy+km] + mmm[im+jm+kz]
1293 + mmm[ix+jm+kp] + mmm[im+jy+kp] + mmm[im+jp+kz]
1294 + mmm[ix+jp+km] + mmm[ip+jy+km] + mmm[ip+jm+kz]
1295 + mmm[ix+jp+kp] + mmm[ip+jy+kp] + mmm[ip+jp+kz];
1296
1297 /* NN3, where all three axes are not at the center */
1298 if( NN >= 3 )
1299 /* im column ip column */
1300 num += mmm[im+jm+km] + mmm[ip+jm+km]
1301 + mmm[im+jm+kp] + mmm[ip+jm+kp]
1302 + mmm[im+jp+km] + mmm[ip+jp+km]
1303 + mmm[im+jp+kp] + mmm[ip+jp+kp];
1304
1305 /* if not enough set neighbors, mark to erode */
1306 if( num <= nmask ) nnn[ii+jy+kz] = 1 ;
1307 } } } /* for kk,jj,ii loops */
1308
1309 STATUS("eroding") ;
1310 for( jj=ii=0 ; ii < nxyz ; ii++ ) /* actually erode */
1311 if( nnn[ii] ){ mmm[ii] = 0 ; jj++ ; }
1312
1313 if( verb && jj > 0 ) ININFO_message("Eroded %d voxels\n",jj) ;
1314
1315 free(nnn) ; EXRETURN ;
1316 }
1317
1318 /*--------------------------------------------------------------------------*/
1319 /*! Generalization of THD_mask_erode(), to peel away multiple layers and
1320 then redilate.
1321 ----------------------------------------------------------------------------*/
1322
THD_mask_erodemany(int nx,int ny,int nz,byte * mmm,int npeel)1323 void THD_mask_erodemany( int nx, int ny, int nz, byte *mmm, int npeel )
1324 {
1325 int ii,jj,kk , jy,kz, im,jm,km , ip,jp,kp , num , pp ;
1326 int nxy=nx*ny , nxyz=nxy*nz ;
1327 byte *nnn,*qqq , bpp , bth ;
1328 int realpeelthr;
1329
1330 ENTRY("THD_mask_erodemany") ;
1331
1332 if( mmm == NULL || npeel < 1 || nxyz < 27 ) EXRETURN ;
1333
1334 /* until allowing for NN1,3 below, allow for NN2 max */
1335 realpeelthr = MIN(18,peelthr);
1336
1337 nnn = (byte *)calloc(sizeof(byte),nxyz) ; /* mask of eroded voxels */
1338 if( nnn == NULL ) EXRETURN ; /* WTF? */
1339 qqq = (byte *)malloc(sizeof(byte)*nxyz) ; /* another copy */
1340 if( qqq == NULL ){ free(nnn); EXRETURN; } /* WTF? */
1341
1342 /* mark interior voxels that don't have 'peelthr' out of 18 nonzero nbhrs */
1343
1344 STATUS("peelings, nothing more than peelings") ;
1345 for( pp=1 ; pp <= npeel ; pp++ ){ /* pp = peel layer index */
1346 bpp = (byte)pp ;
1347 for( kk=0 ; kk < nz ; kk++ ){
1348 kz = kk*nxy ; km = kz-nxy ; kp = kz+nxy ;
1349 if( kk == 0 ) km = kz ;
1350 if( kk == nz-1 ) kp = kz ;
1351
1352 for( jj=0 ; jj < ny ; jj++ ){
1353 jy = jj*nx ; jm = jy-nx ; jp = jy+nx ;
1354 if( jj == 0 ) jm = jy ;
1355 if( jj == ny-1 ) jp = jy ;
1356
1357 for( ii=0 ; ii < nx ; ii++ ){
1358 if( mmm[ii+jy+kz] ){ /* count nonzero nbhrs */
1359 im = ii-1 ; ip = ii+1 ;
1360 if( ii == 0 ) im = 0 ;
1361 if( ii == nx-1 ) ip = ii ;
1362 num = mmm[im+jy+km]
1363 + mmm[ii+jm+km] + mmm[ii+jy+km] + mmm[ii+jp+km]
1364 + mmm[ip+jy+km]
1365 + mmm[im+jm+kz] + mmm[im+jy+kz] + mmm[im+jp+kz]
1366 + mmm[ii+jm+kz] + mmm[ii+jp+kz]
1367 + mmm[ip+jm+kz] + mmm[ip+jy+kz] + mmm[ip+jp+kz]
1368 + mmm[im+jy+kp]
1369 + mmm[ii+jm+kp] + mmm[ii+jy+kp] + mmm[ii+jp+kp]
1370 + mmm[ip+jy+kp] ;
1371 if( num < realpeelthr ) nnn[ii+jy+kz] = bpp ; /* mark to erode */
1372 }
1373 }}}
1374 for( ii=0 ; ii < nxyz ; ii++ ) /* actually erode */
1375 if( nnn[ii] ) mmm[ii] = 0 ;
1376
1377 } /* end of loop over peeling layers */
1378
1379 /* re-dilate eroded voxels that are next to survivors */
1380
1381 STATUS("unpeelings") ;
1382 for( pp=npeel ; pp >= 1 ; pp-- ){ /* loop from innermost peel to outer */
1383 bpp = (byte)pp ; memset(qqq,0,sizeof(byte)*nxyz) ;
1384 bth = (pp==npeel) ? 0 : 1 ;
1385 for( kk=0 ; kk < nz ; kk++ ){
1386 kz = kk*nxy ; km = kz-nxy ; kp = kz+nxy ;
1387 if( kk == 0 ) km = kz ;
1388 if( kk == nz-1 ) kp = kz ;
1389
1390 for( jj=0 ; jj < ny ; jj++ ){
1391 jy = jj*nx ; jm = jy-nx ; jp = jy+nx ;
1392 if( jj == 0 ) jm = jy ;
1393 if( jj == ny-1 ) jp = jy ;
1394
1395 for( ii=0 ; ii < nx ; ii++ ){
1396 if( nnn[ii+jy+kz] >= bpp && !mmm[ii+jy+kz] ){ /* was eroded before */
1397 im = ii-1 ; ip = ii+1 ;
1398 if( ii == 0 ) im = 0 ;
1399 if( ii == nx-1 ) ip = ii ;
1400 qqq[ii+jy+kz] = /* count any surviving nbhrs */
1401 mmm[im+jy+km]
1402 + mmm[ii+jm+km] + mmm[ii+jy+km] + mmm[ii+jp+km]
1403 + mmm[ip+jy+km]
1404 + mmm[im+jm+kz] + mmm[im+jy+kz] + mmm[im+jp+kz]
1405 + mmm[ii+jm+kz] + mmm[ii+jp+kz]
1406 + mmm[ip+jm+kz] + mmm[ip+jy+kz] + mmm[ip+jp+kz]
1407 + mmm[im+jy+kp]
1408 + mmm[ii+jm+kp] + mmm[ii+jy+kp] + mmm[ii+jp+kp]
1409 + mmm[ip+jy+kp] ;
1410 }
1411 }}} /* end of ii,jj,kk loops */
1412
1413 /* actually do the dilation */
1414
1415 for( ii=0 ; ii < nxyz ; ii++ )
1416 if( qqq[ii] > bth ) mmm[ii] = 1 ;
1417
1418 } /* end of redilate loop */
1419
1420 free(qqq); free(nnn); EXRETURN;
1421 }
1422
1423 /*--------------------------------------------------------------------------*/
1424 /*! Dilate a mask - that is, fill in zero voxels that have at least ndil
1425 neighbors in the mask. The neighbors are the 8/18/26 voxels closest
1426 in 3D defined by the NN Nearest Neighbor mode NN1=facing,NN2=edge,NN3=corners
1427 Past default has been NN2.
1428 * not sure about usefulness of ndil threshold, but keeping it
1429 because it's used across many programs now [DRG]
1430 Return value is number of voxels added to the mask.
1431 ----------------------------------------------------------------------------*/
1432
THD_mask_dilate(int nx,int ny,int nz,byte * mmm,int ndil,byte NN)1433 int THD_mask_dilate( int nx, int ny, int nz, byte *mmm , int ndil, byte NN )
1434 {
1435 int ii,jj,kk , jy,kz, im,jm,km , ip,jp,kp , num ;
1436 int nxy=nx*ny , nxyz=nxy*nz , nadd ;
1437 byte *nnn ;
1438
1439 if( mmm == NULL ) return 0 ;
1440 if( ndil < 1 ) ndil = 1 ;
1441
1442 nnn = (byte*)calloc(sizeof(byte),nxyz) ; /* mask of dilated voxels */
1443
1444 /* mark exterior voxels neighboring enough interior voxels */
1445
1446 for( kk=0 ; kk < nz ; kk++ ){
1447 kz = kk*nxy ; km = kz-nxy ; kp = kz+nxy ;
1448 if( kk == 0 ) km = kz ;
1449 if( kk == nz-1 ) kp = kz ;
1450
1451 for( jj=0 ; jj < ny ; jj++ ){
1452 jy = jj*nx ; jm = jy-nx ; jp = jy+nx ;
1453 if( jj == 0 ) jm = jy ;
1454 if( jj == ny-1 ) jp = jy ;
1455
1456 for( ii=0 ; ii < nx ; ii++ ){
1457 if( mmm[ii+jy+kz] == 0 ){ /* consider 0 voxels, count nonzero nbhrs */
1458 im = ii-1 ; ip = ii+1 ;
1459 if( ii == 0 ) im = 0 ;
1460 if( ii == nx-1 ) ip = ii ;
1461 num = /* count non-zero nbhrs */
1462 mmm[ii+jy+km]
1463 + mmm[im+jy+kz]
1464 + mmm[ii+jm+kz] + mmm[ii+jp+kz]
1465 + mmm[ip+jy+kz]
1466 + mmm[ii+jy+kp];
1467
1468 if(NN>=2)
1469 num += // count nn2 neighbors
1470 mmm[im+jy+km]
1471 + mmm[ii+jm+km] + mmm[ii+jp+km]
1472 + mmm[ip+jy+km]
1473 + mmm[im+jm+kz] + mmm[im+jp+kz]
1474
1475 + mmm[ip+jm+kz] + mmm[ip+jp+kz]
1476 + mmm[im+jy+kp]
1477 + mmm[ii+jm+kp] + mmm[ii+jp+kp]
1478 + mmm[ip+jy+kp];
1479
1480 if(NN==3) // consider corners (corners of -1,+1 slices)
1481 num +=
1482 mmm[im+jm+km] + mmm[im+jp+km]
1483 + mmm[ip+jm+km] + mmm[ip+jp+km]
1484 + mmm[im+jm+kp] + mmm[im+jp+kp]
1485 + mmm[ip+jm+kp] + mmm[ip+jp+kp];
1486
1487 if( num >= ndil ) nnn[ii+jy+kz] = 1 ; /* mark to dilate */
1488 } // end if voxel not already in mask
1489
1490 } } } // end for kk,jj,ii
1491
1492 for( nadd=ii=0 ; ii < nxyz ; ii++ ) /* copy it to */
1493 if( nnn[ii] && !mmm[ii] ){ mmm[ii] = 1; nadd++; } /* input mask */
1494
1495 free(nnn) ; return nadd ;
1496 }
1497
1498 /*---------------------------------------------------------------------*/
1499 /* clip in autobox by default but allow turning it off */
1500 /* thanks Judd */
1501
1502 static int bbox_clip=1 ;
THD_autobbox_clip(int c)1503 void THD_autobbox_clip( int c ){ bbox_clip = c; }
1504
1505 /* extra padding? (should be non-negative) */
1506
1507 static int bbox_npad=0 ;
THD_autobbox_npad(int c)1508 void THD_autobbox_npad( int c ){ bbox_npad = c; }
1509
1510 /* don't allow padding to enlarge dataset? [08 Jan 2019] */
1511
1512 static int bbox_noexpand=0 ;
THD_autobbox_noexpand(int c)1513 void THD_autobbox_noexpand( int c ){ bbox_noexpand = c; }
1514
1515 /*---------------------------------------------------------------------*/
1516 /*! Find a bounding box for the main cluster of large-ish voxels.
1517 [xm..xp] will be box for x index, etc.
1518 Send in a prefix and you shall get a cropped volume back. Else
1519 you get NULL and you will like it!
1520 -----------------------------------------------------------------------*/
1521
THD_autobbox(THD_3dim_dataset * dset,int * xm,int * xp,int * ym,int * yp,int * zm,int * zp,char * prefix)1522 THD_3dim_dataset * THD_autobbox( THD_3dim_dataset *dset ,
1523 int *xm, int *xp , int *ym, int *yp , int *zm, int *zp,
1524 char *prefix)
1525 {
1526 MRI_IMAGE *medim ;
1527 THD_3dim_dataset *cropped = NULL;
1528 float clip_val , *mar ;
1529 int nvox , ii ;
1530
1531 ENTRY("THD_autobbox") ;
1532
1533 /* get brick of max absolute value at each voxel */
1534
1535 medim = THD_maxabs_brick(dset) ; if( medim == NULL ) RETURN(NULL) ;
1536
1537 mar = MRI_FLOAT_PTR(medim) ;
1538 nvox = medim->nvox ;
1539 for( ii=0 ; ii < nvox ; ii++ ) mar[ii] = fabs(mar[ii]) ;
1540 if( bbox_clip ){
1541 clip_val = THD_cliplevel(medim,clfrac) ;
1542 for( ii=0 ; ii < nvox ; ii++ )
1543 if( mar[ii] < clip_val ) mar[ii] = 0.0 ;
1544 }
1545
1546 MRI_autobbox( medim , xm,xp , ym,yp , zm,zp ) ;
1547 mri_free(medim) ;
1548
1549 if (prefix) {
1550 int nx=DSET_NX(dset), ny=DSET_NY(dset), nz=DSET_NZ(dset) ;
1551 int xb,xt , yb,yt , zb,zt ;
1552
1553 xb = *xm - bbox_npad ; if( bbox_noexpand && xb < 0 ) xb = 0 ;
1554 xt = *xp + bbox_npad ; if( bbox_noexpand && xt >= nx ) xt = nx-1 ;
1555 yb = *ym - bbox_npad ; if( bbox_noexpand && yb < 0 ) yb = 0 ;
1556 yt = *yp + bbox_npad ; if( bbox_noexpand && yt >= ny ) yt = ny-1 ;
1557 zb = *zm - bbox_npad ; if( bbox_noexpand && zb < 0 ) zb = 0 ;
1558 zt = *zp + bbox_npad ; if( bbox_noexpand && zt >= nz ) zt = nz-1 ;
1559
1560 if( ! (bbox_noexpand &&
1561 xb == 0 && xt == nx-1 &&
1562 yb == 0 && yt == ny-1 && zb == 0 && zt == nz-1) ){
1563 cropped = THD_zeropad( dset ,
1564 -xb , xt-(nx-1) ,
1565 -yb , yt-(ny-1) ,
1566 -zb , zt-(nz-1) ,
1567 prefix , ZPAD_IJK ) ;
1568 if( cropped == NULL ) /* should not transpire */
1569 ERROR_message("THD_autobbox: could not create cropped volume :(") ;
1570 }
1571 }
1572 RETURN(cropped) ;
1573 }
1574
1575 /*------------------------------------------------------------------------*/
1576
1577 static int bbox_clust=1 ;
MRI_autobbox_clust(int c)1578 void MRI_autobbox_clust( int c ){ bbox_clust = c; }
1579
MRI_autobbox(MRI_IMAGE * qim,int * xm,int * xp,int * ym,int * yp,int * zm,int * zp)1580 void MRI_autobbox( MRI_IMAGE *qim ,
1581 int *xm, int *xp , int *ym, int *yp , int *zm, int *zp )
1582 {
1583 MRI_IMAGE *fim ;
1584 float *mar ;
1585 byte *mmm = NULL ;
1586 int nvox , ii,jj,kk , nmm , nx,ny,nz,nxy ;
1587
1588 ENTRY("MRI_autobbox") ;
1589
1590 ASSIF(xm,0) ; ASSIF(xp,0) ; /* initialize output params */
1591 ASSIF(ym,0) ; ASSIF(yp,0) ;
1592 ASSIF(zm,0) ; ASSIF(zp,0) ;
1593
1594 /* find largest component as in first part of THD_automask() */
1595
1596 if( qim->kind != MRI_float ) fim = mri_to_float(qim) ;
1597 else fim = qim ;
1598
1599 nvox = fim->nvox ;
1600 mar = MRI_FLOAT_PTR(fim) ;
1601 mmm = (byte *) calloc( sizeof(byte) , nvox ) ;
1602 for( nmm=ii=0 ; ii < nvox ; ii++ )
1603 if( mar[ii] != 0.0 ){ mmm[ii] = 1; nmm++; }
1604
1605 if( fim != qim ) mri_free(fim) ;
1606
1607 if( nmm == 0 ){ free(mmm); EXRETURN; }
1608
1609 nx = qim->nx; ny = qim->ny; nz = qim->nz; nxy = nx*ny;
1610
1611 if( bbox_clust ){
1612 THD_mask_clust( nx,ny,nz, mmm ) ;
1613 THD_mask_erodemany( nx,ny,nz, mmm, peelcount ) ;
1614 THD_mask_clust( nx,ny,nz, mmm ) ;
1615 }
1616
1617 /* For each plane direction,
1618 find the first and last index that have nonzero voxels in that plane */
1619
1620 if( xm != NULL ){
1621 for( ii=0 ; ii < nx ; ii++ )
1622 for( kk=0 ; kk < nz ; kk++ )
1623 for( jj=0 ; jj < ny ; jj++ )
1624 if( mmm[ii+jj*nx+kk*nxy] ) goto CP5 ;
1625 CP5: ASSIF(xm,ii) ;
1626 }
1627
1628 if( xp != NULL ){
1629 for( ii=nx-1 ; ii >= 0 ; ii-- )
1630 for( kk=0 ; kk < nz ; kk++ )
1631 for( jj=0 ; jj < ny ; jj++ )
1632 if( mmm[ii+jj*nx+kk*nxy] ) goto CP6 ;
1633 CP6: ASSIF(xp,ii) ;
1634 }
1635
1636 if( ym != NULL ){
1637 for( jj=0 ; jj < ny ; jj++ )
1638 for( kk=0 ; kk < nz ; kk++ )
1639 for( ii=0 ; ii < nx ; ii++ )
1640 if( mmm[ii+jj*nx+kk*nxy] ) goto CP3 ;
1641 CP3: ASSIF(ym,jj) ;
1642 }
1643
1644 if( yp != NULL ){
1645 for( jj=ny-1 ; jj >= 0 ; jj-- )
1646 for( kk=0 ; kk < nz ; kk++ )
1647 for( ii=0 ; ii < nx ; ii++ )
1648 if( mmm[ii+jj*nx+kk*nxy] ) goto CP4 ;
1649 CP4: ASSIF(yp,jj) ;
1650 }
1651
1652 if( zm != NULL ){
1653 for( kk=0 ; kk < nz ; kk++ )
1654 for( jj=0 ; jj < ny ; jj++ )
1655 for( ii=0 ; ii < nx ; ii++ )
1656 if( mmm[ii+jj*nx+kk*nxy] ) goto CP1 ;
1657 CP1: ASSIF(zm,kk) ;
1658 }
1659
1660 if( zp != NULL ){
1661 for( kk=nz-1 ; kk >= 0 ; kk-- )
1662 for( jj=0 ; jj < ny ; jj++ )
1663 for( ii=0 ; ii < nx ; ii++ )
1664 if( mmm[ii+jj*nx+kk*nxy] ) goto CP2 ;
1665 CP2: ASSIF(zp,kk) ;
1666 }
1667
1668 free(mmm) ; EXRETURN ;
1669 }
1670
1671 /*------------------------------------------------------------------------*/
1672
MRI_autobbox_byte(MRI_IMAGE * qim,int * xm,int * xp,int * ym,int * yp,int * zm,int * zp)1673 void MRI_autobbox_byte( MRI_IMAGE *qim ,
1674 int *xm, int *xp , int *ym, int *yp , int *zm, int *zp )
1675 {
1676 byte *mmm ;
1677 int nvox , ii,jj,kk , nmm , nx,ny,nz,nxy ;
1678
1679 ENTRY("MRI_autobbox_byte") ;
1680
1681 ASSIF(xm,0) ; ASSIF(xp,0) ; /* initialize output params */
1682 ASSIF(ym,0) ; ASSIF(yp,0) ;
1683 ASSIF(zm,0) ; ASSIF(zp,0) ;
1684
1685 if( qim->kind != MRI_byte ) EXRETURN ; /* bad input */
1686
1687 mmm = MRI_BYTE_PTR(qim) ; if( qim == NULL ) EXRETURN ;
1688
1689 nx = qim->nx; ny = qim->ny; nz = qim->nz; nxy = nx*ny; nvox = nxy*nz ;
1690
1691 /* For each plane direction,
1692 find the first and last index that have nonzero voxels in that plane */
1693
1694 if( xm != NULL ){
1695 for( ii=0 ; ii < nx ; ii++ )
1696 for( kk=0 ; kk < nz ; kk++ )
1697 for( jj=0 ; jj < ny ; jj++ )
1698 if( mmm[ii+jj*nx+kk*nxy] ) goto CP5 ;
1699 CP5: ASSIF(xm,ii) ;
1700 }
1701
1702 if( xp != NULL ){
1703 for( ii=nx-1 ; ii >= 0 ; ii-- )
1704 for( kk=0 ; kk < nz ; kk++ )
1705 for( jj=0 ; jj < ny ; jj++ )
1706 if( mmm[ii+jj*nx+kk*nxy] ) goto CP6 ;
1707 CP6: ASSIF(xp,ii) ;
1708 }
1709
1710 if( ym != NULL ){
1711 for( jj=0 ; jj < ny ; jj++ )
1712 for( kk=0 ; kk < nz ; kk++ )
1713 for( ii=0 ; ii < nx ; ii++ )
1714 if( mmm[ii+jj*nx+kk*nxy] ) goto CP3 ;
1715 CP3: ASSIF(ym,jj) ;
1716 }
1717
1718 if( yp != NULL ){
1719 for( jj=ny-1 ; jj >= 0 ; jj-- )
1720 for( kk=0 ; kk < nz ; kk++ )
1721 for( ii=0 ; ii < nx ; ii++ )
1722 if( mmm[ii+jj*nx+kk*nxy] ) goto CP4 ;
1723 CP4: ASSIF(yp,jj) ;
1724 }
1725
1726 if( zm != NULL ){
1727 for( kk=0 ; kk < nz ; kk++ )
1728 for( jj=0 ; jj < ny ; jj++ )
1729 for( ii=0 ; ii < nx ; ii++ )
1730 if( mmm[ii+jj*nx+kk*nxy] ) goto CP1 ;
1731 CP1: ASSIF(zm,kk) ;
1732 }
1733
1734 if( zp != NULL ){
1735 for( kk=nz-1 ; kk >= 0 ; kk-- )
1736 for( jj=0 ; jj < ny ; jj++ )
1737 for( ii=0 ; ii < nx ; ii++ )
1738 if( mmm[ii+jj*nx+kk*nxy] ) goto CP2 ;
1739 CP2: ASSIF(zp,kk) ;
1740 }
1741
1742 EXRETURN ;
1743 }
1744
1745 /*------------------------------------------------------------------------*/
1746
THD_peel_mask(int nx,int ny,int nz,byte * mmm,int pdepth)1747 int THD_peel_mask( int nx, int ny, int nz , byte *mmm, int pdepth )
1748 {
1749 int nxy=nx*ny , ii,jj,kk , ijk , bot,top , pd=pdepth ;
1750 /* int nxp=nx-pd , nyp=ny-pd , nzp=nz-pd ; */
1751 int num=0 , dnum , nite ;
1752
1753 for( nite=0 ; nite < pd ; nite++ ){
1754 dnum = 0 ;
1755
1756 for( kk=0 ; kk < nz ; kk++ ){
1757 for( jj=0 ; jj < ny ; jj++ ){
1758 ijk = jj*nx + kk*nxy ;
1759 for( bot=0 ; bot < nx && !mmm[bot+ijk]; bot++ ) ;
1760 top = bot+pd ; if( top >= nx ) continue ;
1761 for( ii=bot+1 ; ii <= top && mmm[ii+ijk] ; ii++ ) ;
1762 if( ii <= top ){ mmm[bot+ijk] = 0; dnum++; }
1763 }}
1764
1765 for( kk=0 ; kk < nz ; kk++ ){
1766 for( jj=0 ; jj < ny ; jj++ ){
1767 ijk = jj*nx + kk*nxy ;
1768 for( top=nx-1 ; top >= 0 && !mmm[top+ijk]; top-- ) ;
1769 bot = top-pd ; if( bot < 0 ) continue ;
1770 for( ii=top-1 ; ii >= bot && mmm[ii+ijk] ; ii-- ) ;
1771 if( ii >= bot ){ mmm[top+ijk] = 0; dnum++; }
1772 }}
1773
1774 for( kk=0 ; kk < nz ; kk++ ){
1775 for( ii=0 ; ii < nx ; ii++ ){
1776 ijk = ii + kk*nxy ;
1777 for( bot=0 ; bot < ny && !mmm[bot*nx+ijk] ; bot++ ) ;
1778 top = bot+pd ;
1779 if( top >= ny ) continue ;
1780 for( jj=bot+1 ; jj <= top && mmm[jj*nx+ijk] ; jj++ ) ;
1781 if( jj <= top ){ mmm[bot*nx+ijk] = 0; dnum++; }
1782 }}
1783
1784 for( kk=0 ; kk < nz ; kk++ ){
1785 for( ii=0 ; ii < nx ; ii++ ){
1786 ijk = ii + kk*nxy ;
1787 for( top=ny-1 ; top >= 0 && !mmm[top*nx+ijk] ; top-- ) ;
1788 bot = top-pd ; if( bot < 0 ) continue ;
1789 for( jj=top-1 ; jj >= bot && mmm[jj*nx+ijk] ; jj-- ) ;
1790 if( jj >= bot ){ mmm[top*nx+ijk] = 0; dnum++; }
1791 }}
1792
1793 for( jj=0 ; jj < ny ; jj++ ){
1794 for( ii=0 ; ii < nx ; ii++ ){
1795 ijk = ii + jj*nx ;
1796 for( top=nz-1 ; top >= 0 && !mmm[top*nxy+ijk] ; top-- ) ;
1797 bot = top-pd ; if( bot < 0 ) continue ;
1798 for( kk=top-1 ; kk >= bot && mmm[kk*nxy+ijk] ; kk-- ) ;
1799 if( kk >= bot ){ mmm[top*nxy+ijk] = 0; dnum++; }
1800 }}
1801
1802 num += dnum ;
1803 if( dnum == 0 ) break ;
1804 }
1805
1806 return num ;
1807 }
1808
1809 /*!
1810 Return a vector representing the number of layers peeled
1811 to reach a voxel in mask. ZSS March 02 2010
1812 nn is 1,2 or 3 for nearest neighbor mode NN1, NN2, NN3
1813 */
THD_mask_depth(int nx,int ny,int nz,byte * mask,byte preservemask,short * usethisdepth,byte NN)1814 short *THD_mask_depth ( int nx, int ny, int nz, byte *mask,
1815 byte preservemask, short *usethisdepth, byte NN)
1816 {
1817 int ii, np, niter, ncpmask, nxyz;
1818 byte *cpmask=NULL;
1819 short *depth = NULL;
1820
1821 if ((nxyz = nx * ny * nz) < 0) {
1822 if (verb) ERROR_message("Bad dims");
1823 return(NULL);
1824 }
1825
1826 if (preservemask) {
1827 cpmask = (byte *)malloc(nxyz*sizeof(byte));
1828 memcpy(cpmask, mask, nxyz*sizeof(byte));
1829 } else {
1830 cpmask = mask;
1831 }
1832
1833 if (!cpmask) {
1834 if (verb) ERROR_message("NULL mask (or mask copy) pointer");
1835 return(NULL);
1836 }
1837
1838 if (usethisdepth) {
1839 depth = usethisdepth;
1840 } else {
1841 if (!(depth = (short *)calloc(nxyz, sizeof(short)))) {
1842 if (verb) ERROR_message("Failed to allocate for %d shorts", nxyz);
1843 return(NULL);
1844 }
1845 }
1846 if (!depth) {
1847 if (verb) ERROR_message("NULL depth vector");
1848 return(NULL);
1849 }
1850
1851 ncpmask = THD_countmask( nxyz , cpmask );
1852
1853 niter=0;
1854 while ( ncpmask > 0) {
1855 for (ii=0 ; ii < nxyz; ++ii) {
1856 if (cpmask[ii]) ++depth[ii];
1857 }
1858 /* peel it */
1859 THD_mask_erode( nx, ny, nz, cpmask, 0, NN ) ;
1860 np = ncpmask - THD_countmask( nxyz , cpmask );
1861 if (verb) INFO_message("Peeled %d voxels from mask of %d (now %d)\n",
1862 np, ncpmask,
1863 ncpmask - np) ;
1864
1865 ncpmask -= np;
1866 ++niter;
1867 if (!np && ncpmask) {
1868 WARNING_message("Nothing left to peel, after %d iterations.\n"
1869 " however %d voxels remain in cpmask!\n"
1870 " Jumping ship.\n",
1871 niter, ncpmask);
1872 break;
1873 }
1874 if (ncpmask < 0) {
1875 ERROR_message("Behavioral problems. ncpmask is < 0!\n"
1876 "Hiding head in sand.");
1877 break;
1878 }
1879 }
1880
1881 if (cpmask != mask) free(cpmask); cpmask=NULL;
1882
1883 return(depth);
1884 }
1885
1886 /*------------------------------------------------------------------*/
1887
1888 #undef DALL
1889 #define DALL 128
1890
1891 /*! Put (i,j,k) into the current cluster, if it is nonzero. */
1892
1893 # undef CPUT
1894 # define CPUT(i,j) \
1895 do{ ijk = (i) + (j)*nx ; \
1896 if( mmm[ijk] ){ \
1897 if( nnow == nall ){ /* increase array lengths */ \
1898 nall += DALL + nall/4 ; \
1899 inow = (short *) realloc(inow,sizeof(short)*nall) ; \
1900 jnow = (short *) realloc(jnow,sizeof(short)*nall) ; \
1901 } \
1902 inow[nnow] = i; jnow[nnow] = j; nnow++; mmm[ijk] = 0; \
1903 } } while(0)
1904
1905 /*------------------------------------------------------------------*/
1906 /*! Find the biggest cluster of nonzeros in the byte mask mmm,
1907 and keep all clusters at least kfrac times that size
1908 (0 < kfrac <= 1, duh).
1909 *//*----------------------------------------------------------------*/
1910
THD_mask_clust2D(int nx,int ny,float kfrac,byte * mmm)1911 void THD_mask_clust2D( int nx, int ny, float kfrac, byte *mmm )
1912 {
1913 int ii,jj, icl , nxy, ijk , ijk_last ;
1914 int ip,jp, im,jm , nbest ;
1915 short *inow , *jnow ; int nall , nnow , nkeep ;
1916
1917 int clust_num=0 ;
1918 int *clust_len=NULL ;
1919 short **clust_iii=NULL , **clust_jjj=NULL ;
1920
1921 ENTRY("THD_mask_clust2D") ;
1922
1923 if( mmm == NULL ) EXRETURN ;
1924
1925 nxy = nx*ny ;
1926 nbest = 0 ;
1927
1928 /*--- scan through array, find nonzero point, build a cluster, ... ---*/
1929
1930 ijk_last = 0 ;
1931 while(1) {
1932 /* find next nonzero point */
1933
1934 for( ijk=ijk_last ; ijk < nxy ; ijk++ ) if( mmm[ijk] ) break ;
1935 if( ijk == nxy ) break ; /* didn't find any! */
1936
1937 ijk_last = ijk+1 ; /* start here next time */
1938
1939 /* init current cluster list with this point */
1940
1941 mmm[ijk] = 0 ; /* clear found point */
1942 nall = 16 ; /* # allocated pts */
1943 nnow = 1 ; /* # pts in cluster */
1944 inow = (short *) malloc(sizeof(short)*16) ; /* coords of pts */
1945 jnow = (short *) malloc(sizeof(short)*16) ;
1946 inow[0] = ijk % nx ; jnow[0] = ijk / nx ;
1947
1948 /*--
1949 for each point in cluster:
1950 check neighboring points for nonzero entries in mmm
1951 enter those into cluster (and clear them in mmm)
1952 continue until end of cluster is reached
1953 (note that cluster is expanding as we progress)
1954 --*/
1955
1956 for( icl=0 ; icl < nnow ; icl++ ){
1957 ii = inow[icl] ; jj = jnow[icl] ;
1958 im = ii-1 ; jm = jj-1 ;
1959 ip = ii+1 ; jp = jj+1 ;
1960
1961 if( im >= 0 ) CPUT(im,jj) ;
1962 if( ip < nx ) CPUT(ip,jj) ;
1963 if( jm >= 0 ) CPUT(ii,jm) ;
1964 if( jp < ny ) CPUT(ii,jp) ;
1965 }
1966
1967 /* save all clusters for later processing */
1968
1969 clust_num++ ;
1970 clust_len = (int * )realloc(clust_len,sizeof(int )*clust_num) ;
1971 clust_iii = (short **)realloc(clust_iii,sizeof(short *)*clust_num) ;
1972 clust_jjj = (short **)realloc(clust_jjj,sizeof(short *)*clust_num) ;
1973 clust_len[clust_num-1] = nnow ;
1974 clust_iii[clust_num-1] = inow ;
1975 clust_jjj[clust_num-1] = jnow ;
1976
1977 if( nnow > nbest ) nbest = nnow ; /* nbest = size of biggest cluster */
1978
1979 } /* loop ends when all nonzero points are clustered */
1980
1981 /* put 1's back in at all points from biggest clusters */
1982
1983 nkeep = (int)(kfrac * nbest) ;
1984 if( nkeep <= 0 || nkeep > nbest ) nkeep = nbest ;
1985
1986 for( icl=0 ; icl < clust_num ; icl++ ){
1987 nnow = clust_len[icl] ;
1988 if( nnow >= nkeep ){
1989 inow = clust_iii[icl] ; jnow = clust_jjj[icl] ;
1990 for( ii=0 ; ii < nnow ; ii++ ) mmm[ inow[ii] + jnow[ii]*nx ] = 1 ;
1991 }
1992 free(clust_iii[icl]) ; free(clust_jjj[icl]) ;
1993 }
1994 free(clust_iii); free(clust_jjj); free(clust_len) ;
1995
1996 EXRETURN ;
1997 }
1998
1999 /*--------------------------------------------------------------------------*/
2000 /* 2D version of peeling/restoring code [14 Sep 2020 - birthday of SL Cox!] */
2001 /*--------------------------------------------------------------------------*/
2002
THD_mask_erodemany2D(int nx,int ny,byte * mmm,int npeel)2003 void THD_mask_erodemany2D( int nx, int ny, byte *mmm, int npeel )
2004 {
2005 int ii,jj , jy, im,jm , ip,jp , num , pp ;
2006 int nxy=nx*ny ;
2007 byte *nnn,*qqq , bpp , bth ;
2008 const int realpeelthr = 6 ;
2009
2010 ENTRY("THD_mask_erodemany2D") ;
2011
2012 if( mmm == NULL || npeel < 1 || nxy < 16 ) EXRETURN ;
2013
2014 nnn = (byte *)calloc(sizeof(byte),nxy) ; /* mask of eroded voxels */
2015 if( nnn == NULL ) EXRETURN ; /* WTF? */
2016 qqq = (byte *)malloc(sizeof(byte)*nxy) ; /* another copy */
2017 if( qqq == NULL ){ free(nnn); EXRETURN; } /* WTF-squared? */
2018
2019 /* mark interior voxels that don't have 'peelthr' out of 8 nonzero nbhrs */
2020
2021 STATUS("peelings, nothing more than peelings") ;
2022 for( pp=1 ; pp <= npeel ; pp++ ){ /* pp = peel layer index */
2023
2024 bpp = (byte)pp ;
2025
2026 for( jj=0 ; jj < ny ; jj++ ){
2027 jy = jj*nx ; jm = jy-nx ; jp = jy+nx ;
2028 if( jj == 0 ) jm = jy ;
2029 if( jj == ny-1 ) jp = jy ;
2030
2031 for( ii=0 ; ii < nx ; ii++ ){
2032 if( mmm[ii+jy] ){ /* count nonzero nbhrs */
2033 im = ii-1 ; ip = ii+1 ;
2034 if( ii == 0 ) im = 0 ;
2035 if( ii == nx-1 ) ip = ii ;
2036 num = mmm[im+jm] + mmm[ii+jm] + mmm[ip+jm]
2037 + mmm[im+jy] + mmm[ip+jy]
2038 + mmm[im+jp] + mmm[ii+jp] + mmm[ip+jp] ;
2039 if( num < realpeelthr ) nnn[ii+jy] = bpp ; /* mark to erode */
2040 }
2041 }} /* end of ii,jj loops */
2042
2043 for( ii=0 ; ii < nxy ; ii++ ) /* actually erode */
2044 if( nnn[ii] ) mmm[ii] = 0 ;
2045
2046 } /* end of loop over peeling layers */
2047
2048 /* re-dilate eroded voxels that are next to survivors */
2049
2050 STATUS("unpeelings") ;
2051 for( pp=npeel ; pp >= 1 ; pp-- ){ /* loop from innermost peel to outer */
2052
2053 bpp = (byte)pp ; memset(qqq,0,sizeof(byte)*nxy) ;
2054 bth = (pp==npeel) ? 0 : 1 ; /* threshold */
2055
2056 for( jj=0 ; jj < ny ; jj++ ){
2057 jy = jj*nx ; jm = jy-nx ; jp = jy+nx ;
2058 if( jj == 0 ) jm = jy ;
2059 if( jj == ny-1 ) jp = jy ;
2060
2061 for( ii=0 ; ii < nx ; ii++ ){
2062 if( nnn[ii+jy] >= bpp && !mmm[ii+jy] ){ /* was eroded before */
2063 im = ii-1 ; ip = ii+1 ;
2064 if( ii == 0 ) im = 0 ;
2065 if( ii == nx-1 ) ip = ii ;
2066 qqq[ii+jy] = /* count any surviving nbhrs */
2067 mmm[im+jm] + mmm[ii+jm] + mmm[ip+jm]
2068 + mmm[im+jy] + mmm[ip+jy]
2069 + mmm[im+jp] + mmm[ii+jp] + mmm[ip+jp] ;
2070 }
2071 }} /* end of ii,jj loops */
2072
2073 /* actually do the dilation */
2074
2075 for( ii=0 ; ii < nxy ; ii++ )
2076 if( qqq[ii] > bth ) mmm[ii] = 1 ;
2077
2078 } /* end of redilate loop */
2079
2080 free(qqq); free(nnn); EXRETURN;
2081 }
2082
2083 /*---------------------------------------------------------------------*/
2084 /*! Make a byte mask from an image (2D).
2085 -----------------------------------------------------------------------*/
2086
mri_automask_image2D(MRI_IMAGE * im)2087 byte * mri_automask_image2D( MRI_IMAGE *im ) /* 12 Jun 2010 */
2088 {
2089 float clip_val , *mar ;
2090 byte *mmm = NULL ;
2091 int nvox, ii, nmm ;
2092 MRI_IMAGE *medim ;
2093
2094 ENTRY("mri_automask_image2D") ;
2095
2096 if( im == NULL || im->nx < 16 || im->ny < 16 ) RETURN(NULL) ;
2097
2098 medim = mri_to_float(im) ;
2099 mar = MRI_FLOAT_PTR(medim) ;
2100 nvox = medim->nvox ;
2101 for( ii=0 ; ii < nvox ; ii++ ) mar[ii] = fabsf(mar[ii]) ;
2102
2103 clip_val = THD_cliplevel(medim,clfrac) ;
2104 mmm = (byte *) calloc( sizeof(byte), nvox ) ;
2105
2106 for( nmm=ii=0 ; ii < nvox ; ii++ )
2107 if( mar[ii] >= clip_val ){ mmm[ii] = 1; nmm++; }
2108
2109 mri_free(medim) ;
2110 if( nmm == 0 ){ free(mmm) ; RETURN(NULL) ; } /* should not happen */
2111 if( nmm <= 2 || nmm == nvox ) RETURN(mmm) ; /* very unlikely */
2112
2113 THD_mask_erodemany2D( im->nx,im->ny, mmm, 3 ) ; /* 14 Sep 2020 */
2114
2115 THD_mask_clust2D( im->nx,im->ny,0.5f, mmm ) ;
2116
2117 for( ii=0 ; ii < nvox ; ii++ ) mmm[ii] = !mmm[ii] ;
2118
2119 THD_mask_clust2D( im->nx,im->ny,0.9f, mmm ) ;
2120
2121 for( ii=0 ; ii < nvox ; ii++ ) mmm[ii] = !mmm[ii] ;
2122
2123 RETURN(mmm) ;
2124 }
2125