Doxygen Source Code Documentation
thd_automask.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | USE_DILATE |
#define | USE_FILLIN |
#define | FILLVOX do{ nnn[iv] = 1; nfill++; goto NextVox; } while(0) |
#define | DALL 1024 |
#define | CPUT(i, j, k) |
Functions | |
void | THD_automask_verbose (int v) |
void | THD_automask_extclip (int e) |
int | mask_count (int nvox, byte *mmm) |
byte * | THD_automask (THD_3dim_dataset *dset) |
byte * | mri_automask_imarr (MRI_IMARR *imar) |
byte * | mri_automask_image (MRI_IMAGE *im) |
int | THD_mask_clip_neighbors (int nx, int ny, int nz, byte *mmm, float clip_val, float tclip, float *mar) |
int | THD_mask_fillin_once (int nx, int ny, int nz, byte *mmm, int nside) |
int | THD_mask_fillin_completely (int nx, int ny, int nz, byte *mmm, int nside) |
void | THD_mask_clust (int nx, int ny, int nz, byte *mmm) |
void | THD_mask_erode (int nx, int ny, int nz, byte *mmm) |
void | THD_mask_dilate (int nx, int ny, int nz, byte *mmm, int ndil) |
void | THD_autobbox (THD_3dim_dataset *dset, int *xm, int *xp, int *ym, int *yp, int *zm, int *zp) |
void | MRI_autobbox (MRI_IMAGE *qim, int *xm, int *xp, int *ym, int *yp, int *zm, int *zp) |
int | THD_peel_mask (int nx, int ny, int nz, byte *mmm, int pdepth) |
Variables | |
int | verb = 0 |
int | exterior_clip = 0 |
int | dall = 0 |
Define Documentation
|
Value: do{ ijk = THREE_TO_IJK(i,j,k,nx,nxy) ; \ if( mmm[ijk] ){ \ if( nnow == nall ){ \ nall += dall ; \ inow = (short *) realloc(inow,sizeof(short)*nall) ; \ jnow = (short *) realloc(jnow,sizeof(short)*nall) ; \ know = (short *) realloc(know,sizeof(short)*nall) ; \ } \ inow[nnow] = i ; jnow[nnow] = j ; know[nnow] = k ; \ nnow++ ; mmm[ijk] = 0 ; \ } } while(0) Definition at line 451 of file thd_automask.c. Referenced by clustedit3D(), and THD_mask_clust(). |
|
Definition at line 447 of file thd_automask.c. Referenced by THD_mask_clust(). |
|
|
|
Definition at line 3 of file thd_automask.c. |
|
Definition at line 4 of file thd_automask.c. |
Function Documentation
|
Count number of nonzeros in mask array Definition at line 18 of file thd_automask.c. References mmm.
00019 { 00020 int ii , nn ; 00021 for( nn=ii=0 ; ii < nvox ; ii++ ) nn += (mmm[ii] != 0) ; 00022 return nn ; 00023 } |
|
Definition at line 730 of file thd_automask.c. References calloc, ENTRY, fim, free, MRI_IMAGE::kind, mmm, MRI_FLOAT_PTR, mri_free(), mri_to_float(), MRI_IMAGE::nvox, MRI_IMAGE::nx, MRI_IMAGE::ny, MRI_IMAGE::nz, nz, THD_mask_clust(), and THD_mask_erode(). Referenced by mri_3dalign_setup(), and THD_autobbox().
00732 { 00733 MRI_IMAGE *fim ; 00734 float *mar ; 00735 byte *mmm = NULL ; 00736 int nvox , ii,jj,kk , nmm , nx,ny,nz,nxy ; 00737 00738 ENTRY("MRI_autobbox") ; 00739 00740 /* find largest component as in first part of THD_automask() */ 00741 00742 if( qim->kind != MRI_float ) fim = mri_to_float(qim) ; 00743 else fim = qim ; 00744 00745 nvox = fim->nvox ; 00746 mar = MRI_FLOAT_PTR(fim) ; 00747 mmm = (byte *) calloc( sizeof(byte) , nvox ) ; 00748 for( nmm=ii=0 ; ii < nvox ; ii++ ) 00749 if( mar[ii] != 0.0 ){ mmm[ii] = 1; nmm++; } 00750 00751 if( fim != qim ) mri_free(fim) ; 00752 00753 if( nmm == 0 ){ free(mmm); EXRETURN; } 00754 00755 nx = qim->nx; ny = qim->ny; nz = qim->nz; nxy = nx*ny; 00756 00757 THD_mask_clust( nx,ny,nz, mmm ) ; 00758 THD_mask_erode( nx,ny,nz, mmm ) ; 00759 THD_mask_clust( nx,ny,nz, mmm ) ; 00760 00761 /* For each plane direction, 00762 find the first and last index that have nonzero voxels in that plane */ 00763 00764 if( xm != NULL ){ 00765 for( ii=0 ; ii < nx ; ii++ ) 00766 for( kk=0 ; kk < nz ; kk++ ) 00767 for( jj=0 ; jj < ny ; jj++ ) 00768 if( mmm[ii+jj*nx+kk*nxy] ) goto CP5 ; 00769 CP5: *xm = ii ; 00770 } 00771 00772 if( xp != NULL ){ 00773 for( ii=nx-1 ; ii >= 0 ; ii-- ) 00774 for( kk=0 ; kk < nz ; kk++ ) 00775 for( jj=0 ; jj < ny ; jj++ ) 00776 if( mmm[ii+jj*nx+kk*nxy] ) goto CP6 ; 00777 CP6: *xp = ii ; 00778 } 00779 00780 if( ym != NULL ){ 00781 for( jj=0 ; jj < ny ; jj++ ) 00782 for( kk=0 ; kk < nz ; kk++ ) 00783 for( ii=0 ; ii < nx ; ii++ ) 00784 if( mmm[ii+jj*nx+kk*nxy] ) goto CP3 ; 00785 CP3: *ym = jj ; 00786 } 00787 00788 if( yp != NULL ){ 00789 for( jj=ny-1 ; jj >= 0 ; jj-- ) 00790 for( kk=0 ; kk < nz ; kk++ ) 00791 for( ii=0 ; ii < nx ; ii++ ) 00792 if( mmm[ii+jj*nx+kk*nxy] ) goto CP4 ; 00793 CP4: *yp = jj ; 00794 } 00795 00796 if( zm != NULL ){ 00797 for( kk=0 ; kk < nz ; kk++ ) 00798 for( jj=0 ; jj < ny ; jj++ ) 00799 for( ii=0 ; ii < nx ; ii++ ) 00800 if( mmm[ii+jj*nx+kk*nxy] ) goto CP1 ; 00801 CP1: *zm = kk ; 00802 } 00803 00804 if( zp != NULL ){ 00805 for( kk=nz-1 ; kk >= 0 ; kk-- ) 00806 for( jj=0 ; jj < ny ; jj++ ) 00807 for( ii=0 ; ii < nx ; ii++ ) 00808 if( mmm[ii+jj*nx+kk*nxy] ) goto CP2 ; 00809 CP2: *zp = kk ; 00810 } 00811 00812 free(mmm) ; EXRETURN ; 00813 } |
|
Make a byte mask from an image (3D). Adapted from THD_automask() to make it possible to do this on an image directly. ----------------------------------------------------------------------- Definition at line 92 of file thd_automask.c. References AFNI_yesenv(), calloc, dall, ENTRY, exterior_clip, ININFO_message(), MRI_IMAGE::kind, mask_count(), MAX, mmm, MRI_FLOAT_PTR, mri_free(), mri_to_float(), MRI_IMAGE::nvox, MRI_IMAGE::nx, MRI_IMAGE::ny, MRI_IMAGE::nz, nz, RETURN, THD_cliplevel(), THD_mask_clip_neighbors(), THD_mask_clust(), THD_mask_erode(), THD_mask_fillin_completely(), THD_mask_fillin_once(), THD_peel_mask(), and verb. Referenced by main(), mri_automask_imarr(), read_input_data(), and THD_automask().
00093 { 00094 float clip_val , *mar ; 00095 byte *mmm = NULL ; 00096 int nvox , ii,jj , nmm , nx,ny,nz ; 00097 MRI_IMAGE *medim ; 00098 00099 ENTRY("mri_automask_image") ; 00100 00101 if( im == NULL ) return NULL ; 00102 00103 if( im->kind != MRI_float ) medim = mri_to_float(im) ; 00104 else medim = im ; 00105 00106 /* find clip value to excise small stuff */ 00107 00108 clip_val = THD_cliplevel(medim,0.5) ; 00109 00110 if( verb ) ININFO_message("Clip level = %f\n",clip_val) ; 00111 00112 /* create mask of values above clip value */ 00113 00114 nvox = medim->nvox ; 00115 mar = MRI_FLOAT_PTR(medim) ; 00116 mmm = (byte *) calloc( sizeof(byte), nvox ) ; 00117 for( nmm=ii=0 ; ii < nvox ; ii++ ) 00118 if( mar[ii] >= clip_val ){ mmm[ii] = 1; nmm++; } 00119 00120 if( AFNI_yesenv("TOPCLIP") ){ 00121 float tclip = 3.1*clip_val ; 00122 if( verb ) ININFO_message("Top clip = %f\n",tclip) ; 00123 for( ii=0 ; ii < nvox ; ii++ ) 00124 if( mar[ii] > tclip ) mmm[ii] = 0 ; 00125 for( nmm=ii=0 ; ii < nvox ; ii++ ) if( mmm[ii] ) nmm++ ; 00126 } 00127 00128 if( verb ) ININFO_message("Number voxels above clip level = %d\n",nmm) ; 00129 if( im != medim && (!exterior_clip || nmm==0) ){ mri_free(medim); medim=NULL; } 00130 if( nmm == 0 ) RETURN(mmm) ; /* should not happen */ 00131 00132 /*-- 10 Apr 2002: only keep the largest connected component --*/ 00133 00134 nx = im->nx ; ny = im->ny ; nz = im->nz ; 00135 dall = (nx*ny*nz)/128 ; /* allocation delta for clustering */ 00136 00137 if( verb ) ININFO_message("Clustering voxels above clip level ...\n") ; 00138 THD_mask_clust( nx,ny,nz, mmm ) ; 00139 00140 /* 18 Apr 2002: now erode the resulting volume 00141 (to break off any thinly attached pieces) */ 00142 00143 THD_mask_erode( nx,ny,nz, mmm ) ; 00144 00145 /* now recluster it, and again keep only the largest survivor */ 00146 00147 THD_mask_clust( nx,ny,nz, mmm ) ; 00148 00149 #ifdef USE_FILLIN 00150 /* 19 Apr 2002: fill in small holes */ 00151 00152 jj = ii = THD_mask_fillin_once( nx,ny,nz , mmm , 1 ) ; 00153 if( ii > 0 ){ 00154 jj += ii = THD_mask_fillin_once( nx,ny,nz , mmm , 1 ) ; 00155 if( ii > 0 ){ 00156 jj += ii = THD_mask_fillin_once( nx,ny,nz , mmm , 1 ) ; 00157 } 00158 } 00159 if( jj > 0 && verb ) 00160 ININFO_message("Filled %d voxels in small holes; now have %d voxels\n", 00161 jj , mask_count(nvox,mmm) ) ; 00162 00163 if( AFNI_yesenv("PEEL") ){ 00164 jj = THD_peel_mask( nx,ny,nz , mmm , 7 ) ; 00165 if( jj > 0 ){ 00166 ININFO_message("Peeled %d voxels from surface\n",jj) ; 00167 THD_mask_erode( nx,ny,nz, mmm ) ; 00168 THD_mask_clust( nx,ny,nz, mmm ) ; 00169 } 00170 } 00171 00172 nmm = 1 ; 00173 jj = rint(0.016*nx) ; nmm = MAX(nmm,jj) ; 00174 jj = rint(0.016*ny) ; nmm = MAX(nmm,jj) ; 00175 jj = rint(0.016*nz) ; nmm = MAX(nmm,jj) ; 00176 00177 if( nmm > 1 || jj > 0 ){ 00178 for( jj=0,ii=2 ; ii < nmm ; ii++ ) 00179 jj += THD_mask_fillin_once( nx,ny,nz , mmm , ii ) ; 00180 jj += THD_mask_fillin_completely( nx,ny,nz, mmm , nmm ) ; 00181 if( jj > 0 && verb ) 00182 ININFO_message("Filled %d voxels in large holes; now have %d voxels\n", 00183 jj , mask_count(nvox,mmm) ) ; 00184 } 00185 00186 THD_mask_erode( nx,ny,nz, mmm ) ; 00187 THD_mask_clust( nx,ny,nz, mmm ) ; 00188 #endif 00189 00190 /* 28 May 2002: 00191 invert the mask, then find the largest cluster of 1's; 00192 this will be the outside of the brain; 00193 put this back into the mask, then invert again; 00194 the effect will be to fill any holes left inside the brain */ 00195 00196 for( ii=0 ; ii < nvox ; ii++ ) mmm[ii] = !mmm[ii] ; 00197 00198 if( verb ) ININFO_message("Clustering non-brain voxels ...\n") ; 00199 THD_mask_clust( nx,ny,nz, mmm ) ; 00200 00201 /* mask is now 1 for non-brain voxels; 00202 if we want to clip off voxels neighboring the non-brain 00203 mask AND whose values are below clip_val, do so now */ 00204 00205 if( exterior_clip ){ 00206 float tclip ; 00207 tclip = AFNI_yesenv("TOPCLIP") ? 3.1*clip_val : 9999.9*clip_val ;; 00208 jj = THD_mask_clip_neighbors( nx,ny,nz , mmm , clip_val,tclip,mar ) ; 00209 if( im != medim ) mri_free(medim) ; 00210 if( jj > 0 && verb ) 00211 ININFO_message("Removed %d exterior voxels below clip level\n",jj); 00212 } else { 00213 jj = 0 ; 00214 } 00215 00216 /* and re-invert mask to get brain voxels */ 00217 00218 for( ii=0 ; ii < nvox ; ii++ ) mmm[ii] = !mmm[ii] ; 00219 if( verb ) ININFO_message("Mask now has %d voxels\n",mask_count(nvox,mmm)) ; 00220 00221 if( exterior_clip && jj > 0 ){ 00222 THD_mask_erode( nx,ny,nz, mmm ) ; 00223 THD_mask_clust( nx,ny,nz, mmm ) ; 00224 } 00225 00226 RETURN(mmm) ; 00227 } |
|
Make a byte mask from the average of an array of 3D images. We assume that they all have the same (nx,ny,nz) dimensions. ----------------------------------------------------------------------- Definition at line 52 of file thd_automask.c. References ENTRY, IMARR_COUNT, IMARR_SUBIMAGE, MRI_IMAGE::kind, mmm, mri_automask_image(), MRI_FLOAT_PTR, mri_free(), mri_new_conforming, mri_to_float(), MRI_IMAGE::nvox, and RETURN.
00053 { 00054 MRI_IMAGE *avim , *tim , *qim ; 00055 byte *mmm ; 00056 int ii , jj , nvox,nim ; 00057 float fac , *avar , *qar ; 00058 00059 ENTRY("mri_automask_imarr") ; 00060 00061 if( imar == NULL || IMARR_COUNT(imar) < 1 ) RETURN(NULL) ; 00062 00063 nim = IMARR_COUNT(imar) ; 00064 if( nim == 1 ){ 00065 mmm = mri_automask_image( IMARR_SUBIMAGE(imar,0) ) ; 00066 RETURN(mmm) ; 00067 } 00068 00069 avim = mri_new_conforming( IMARR_SUBIMAGE(imar,0) , MRI_float ) ; 00070 avar = MRI_FLOAT_PTR(avim) ; 00071 nvox = avim->nvox ; 00072 for( jj=0 ; jj < nim ; jj++ ){ 00073 tim = IMARR_SUBIMAGE(imar,jj) ; 00074 if( tim->kind != MRI_float ) qim = mri_to_float(tim) ; 00075 else qim = tim ; 00076 qar = MRI_FLOAT_PTR(qim) ; 00077 for( ii=0 ; ii < nvox ; ii++ ) avar[ii] += qar[ii] ; 00078 if( qim != tim ) mri_free(qim) ; 00079 } 00080 fac = 1.0f / (float)nim ; 00081 for( ii=0 ; ii < nvox ; ii++ ) avar[ii] *= fac ; 00082 mmm = mri_automask_image( avim ) ; 00083 mri_free(avim) ; 00084 RETURN(mmm) ; 00085 } |
|
Find a bounding box for the main cluster of large-ish voxels. [xm..xp] will be box for x index, etc. ----------------------------------------------------------------------- Definition at line 704 of file thd_automask.c. References ENTRY, MRI_autobbox(), MRI_FLOAT_PTR, mri_free(), MRI_IMAGE::nvox, THD_cliplevel(), and THD_median_brick(). Referenced by main().
00706 { 00707 MRI_IMAGE *medim ; 00708 float clip_val , *mar ; 00709 int nvox , ii ; 00710 00711 ENTRY("THD_autobbox") ; 00712 00713 medim = THD_median_brick(dset) ; if( medim == NULL ) EXRETURN ; 00714 00715 mar = MRI_FLOAT_PTR(medim) ; 00716 nvox = medim->nvox ; 00717 for( ii=0 ; ii < nvox ; ii++ ) mar[ii] = fabs(mar[ii]) ; 00718 00719 clip_val = THD_cliplevel(medim,0.5) ; 00720 for( ii=0 ; ii < nvox ; ii++ ) 00721 if( mar[ii] < clip_val ) mar[ii] = 0.0 ; 00722 00723 MRI_autobbox( medim , xm,xp , ym,yp , zm,zp ) ; 00724 00725 mri_free(medim) ; EXRETURN ; 00726 } |
|
Make a byte mask for a 3D+time dataset -- 13 Aug 2001 - RWCox.
Definition at line 31 of file thd_automask.c. References ENTRY, mmm, mri_automask_image(), mri_free(), RETURN, and THD_median_brick(). Referenced by main().
00032 { 00033 MRI_IMAGE *medim ; 00034 byte *mmm ; 00035 00036 ENTRY("THD_automask") ; 00037 00038 /* median at each voxel */ 00039 00040 medim = THD_median_brick(dset) ; if( medim == NULL ) RETURN(NULL); 00041 00042 mmm = mri_automask_image( medim ) ; 00043 00044 mri_free(medim) ; RETURN(mmm) ; 00045 } |
|
Definition at line 12 of file thd_automask.c. References exterior_clip. Referenced by main().
00012 { exterior_clip = e ; } |
|
Definition at line 9 of file thd_automask.c. Referenced by main().
|
|
Find voxels not in the mask, but neighboring it, that are also below the clip threshold in mar[]. Add these to the mask. Repeat until done. Return value is number of voxels added. ----------------------------------------------------------------------- Definition at line 235 of file thd_automask.c. Referenced by mri_automask_image().
00237 { 00238 int ii,jj,kk , ntot=0,nnew , jm,jp,j3 , km,kp,k3 , im,ip,i3 , nxy=nx*ny ; 00239 00240 if( mmm == NULL || mar == NULL ) return 0 ; 00241 00242 do{ 00243 nnew = 0 ; 00244 for( kk=1 ; kk < nz-1 ; kk++ ){ 00245 k3 = kk*nxy ; 00246 for( jj=1 ; jj < ny-1 ; jj++ ){ 00247 j3 = k3 + jj*nx ; 00248 for( ii=1 ; ii < nx-1 ; ii++ ){ 00249 i3 = ii+j3 ; 00250 if( mmm[i3] || /* in mask */ 00251 (mar[i3] >= clip_val && mar[i3] <= tclip) ) continue ; /* or is OK */ 00252 00253 /* If here, voxel IS NOT in mask, and IS below threshold. 00254 If any neighbors are also in mask, then add it to mask. */ 00255 00256 if( mmm[i3-1] || mmm[i3+1] || 00257 mmm[i3-nx] || mmm[i3+nx] || 00258 mmm[i3-nxy] || mmm[i3+nxy] ){ mmm[i3] = 1; nnew++; } 00259 }}} 00260 ntot += nnew ; 00261 } while( nnew > 0 ) ; 00262 00263 return ntot ; 00264 } |
|
Find the biggest cluster of nonzeros in the byte mask mmm. -------------------------------------------------------------------- Definition at line 468 of file thd_automask.c. References AFMALL, CPUT, DALL, dall, free, IJK_TO_THREE, ININFO_message(), malloc, mmm, nz, and THREE_TO_IJK. Referenced by main(), make_peel_mask(), MRI_autobbox(), mri_automask_image(), mri_brainormalize(), mri_warp3D_align_setup(), THD_orient_guess(), and zedit_mask().
00469 { 00470 int ii,jj,kk, icl , nxy,nxyz , ijk , ijk_last , mnum ; 00471 int ip,jp,kp , im,jm,km ; 00472 int nbest ; short *ibest, *jbest , *kbest ; 00473 int nnow ; short *inow , *jnow , *know ; int nall ; 00474 00475 if( mmm == NULL ) return ; 00476 00477 nxy = nx*ny ; nxyz = nxy * nz ; 00478 00479 nbest = 0 ; 00480 ibest = AFMALL(short, sizeof(short)) ; 00481 jbest = AFMALL(short, sizeof(short)) ; 00482 kbest = AFMALL(short, sizeof(short)) ; 00483 00484 /*--- scan through array, find nonzero point, build a cluster, ... ---*/ 00485 00486 ijk_last = 0 ; if( dall < DALL ) dall = DALL ; 00487 while(1) { 00488 /* find next nonzero point */ 00489 00490 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( mmm[ijk] ) break ; 00491 if( ijk == nxyz ) break ; /* didn't find any! */ 00492 00493 ijk_last = ijk+1 ; /* start here next time */ 00494 00495 /* init current cluster list with this point */ 00496 00497 mmm[ijk] = 0 ; /* clear found point */ 00498 nall = DALL ; /* # allocated pts */ 00499 nnow = 1 ; /* # pts in cluster */ 00500 inow = (short *) malloc(sizeof(short)*DALL) ; /* coords of pts */ 00501 jnow = (short *) malloc(sizeof(short)*DALL) ; 00502 know = (short *) malloc(sizeof(short)*DALL) ; 00503 IJK_TO_THREE(ijk, inow[0],jnow[0],know[0] , nx,nxy) ; 00504 00505 /*-- 00506 for each point in cluster: 00507 check neighboring points for nonzero entries in mmm 00508 enter those into cluster (and clear them in mmm) 00509 continue until end of cluster is reached 00510 (note that cluster is expanding as we progress) 00511 --*/ 00512 00513 for( icl=0 ; icl < nnow ; icl++ ){ 00514 ii = inow[icl] ; jj = jnow[icl] ; kk = know[icl] ; 00515 im = ii-1 ; jm = jj-1 ; km = kk-1 ; 00516 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ; 00517 00518 if( im >= 0 ) CPUT(im,jj,kk) ; 00519 if( ip < nx ) CPUT(ip,jj,kk) ; 00520 if( jm >= 0 ) CPUT(ii,jm,kk) ; 00521 if( jp < ny ) CPUT(ii,jp,kk) ; 00522 if( km >= 0 ) CPUT(ii,jj,km) ; 00523 if( kp < nz ) CPUT(ii,jj,kp) ; 00524 } 00525 00526 /* see if now cluster is larger than best yet */ 00527 00528 if( nnow > nbest ){ /* now is bigger; */ 00529 free(ibest) ; free(jbest) ; free(kbest) ; /* replace best */ 00530 nbest = nnow ; ibest = inow ; /* with now */ 00531 jbest = jnow ; kbest = know ; 00532 } else { /* old is bigger */ 00533 free(inow) ; free(jnow) ; free(know) ; /* toss now */ 00534 } 00535 00536 } /* loop ends when all nonzero points are clustered */ 00537 00538 /* put 1's back in at all points in best cluster */ 00539 00540 for( icl=0 ; icl < nbest ; icl++ ){ 00541 ijk = THREE_TO_IJK(ibest[icl],jbest[icl],kbest[icl],nx,nxy) ; 00542 mmm[ijk] = 1 ; 00543 } 00544 free(ibest) ; free(jbest) ; free(kbest) ; 00545 00546 if( verb ) ININFO_message("Largest cluster has %d voxels\n",nbest) ; 00547 00548 return ; 00549 } |
|
Dilate a mask - that is, fill in zero voxels that have at least ndil neighbors in the mask. The neighbors are the 18 voxels closest in 3D (nearest and next-nearest neighbors). ---------------------------------------------------------------------------- Definition at line 651 of file thd_automask.c. References calloc, free, mmm, and nz. Referenced by main(), and mri_brainormalize().
00652 { 00653 int ii,jj,kk , jy,kz, im,jm,km , ip,jp,kp , num ; 00654 int nxy=nx*ny , nxyz=nxy*nz ; 00655 byte *nnn ; 00656 00657 if( mmm == NULL ) return ; 00658 if( ndil < 1 ) ndil = 1 ; 00659 else if( ndil > 17 ) ndil = 17 ; 00660 00661 nnn = (byte*)calloc(sizeof(byte),nxyz) ; /* mask of dilated voxels */ 00662 00663 /* mark exterior voxels neighboring enough interior voxels */ 00664 00665 for( kk=0 ; kk < nz ; kk++ ){ 00666 kz = kk*nxy ; km = kz-nxy ; kp = kz+nxy ; 00667 if( kk == 0 ) km = kz ; 00668 else if( kk == nz-1 ) kp = kz ; 00669 00670 for( jj=0 ; jj < ny ; jj++ ){ 00671 jy = jj*nx ; jm = jy-nx ; jp = jy+nx ; 00672 if( jj == 0 ) jm = jy ; 00673 else if( jj == ny-1 ) jp = jy ; 00674 00675 for( ii=0 ; ii < nx ; ii++ ){ 00676 if( mmm[ii+jy+kz] == 0 ){ /* count nonzero nbhrs */ 00677 im = ii-1 ; ip = ii+1 ; 00678 if( ii == 0 ) im = 0 ; 00679 else if( ii == nx-1 ) ip = ii ; 00680 num = mmm[im+jy+km] 00681 + mmm[ii+jm+km] + mmm[ii+jy+km] + mmm[ii+jp+km] 00682 + mmm[ip+jy+km] 00683 + mmm[im+jm+kz] + mmm[im+jy+kz] + mmm[im+jp+kz] 00684 + mmm[ii+jm+kz] + mmm[ii+jp+kz] 00685 + mmm[ip+jm+kz] + mmm[ip+jy+kz] + mmm[ip+jp+kz] 00686 + mmm[im+jy+kp] 00687 + mmm[ii+jm+kp] + mmm[ii+jy+kp] + mmm[ii+jp+kp] 00688 + mmm[ip+jy+kp] ; 00689 if( num >= ndil ) nnn[ii+jy+kz] = 1 ; /* mark to dilate */ 00690 } 00691 } } } 00692 00693 for( ii=0 ; ii < nxyz ; ii++ ) /* actually dilate */ 00694 if( nnn[ii] ) mmm[ii] = 1 ; 00695 00696 free(nnn) ; return ; 00697 } |
|
Erode away nonzero voxels that aren't neighbored by mostly other nonzero voxels. Then restore those that were eroded that are neighbors of survivors. The neighbors are the 18 voxels closest in 3D (nearest and next-nearest neighbors). ---------------------------------------------------------------------------- Definition at line 558 of file thd_automask.c. References calloc, free, ININFO_message(), mmm, nz, and verb. Referenced by make_peel_mask(), MRI_autobbox(), mri_automask_image(), mri_brainormalize(), mri_warp3D_align_setup(), and zedit_mask().
00559 { 00560 int ii,jj,kk , jy,kz, im,jm,km , ip,jp,kp , num ; 00561 int nxy=nx*ny , nxyz=nxy*nz ; 00562 byte *nnn ; 00563 00564 if( mmm == NULL ) return ; 00565 00566 nnn = (byte*)calloc(sizeof(byte),nxyz) ; /* mask of eroded voxels */ 00567 00568 /* mark interior voxels that don't have 17 out of 18 nonzero nbhrs */ 00569 00570 for( kk=0 ; kk < nz ; kk++ ){ 00571 kz = kk*nxy ; km = kz-nxy ; kp = kz+nxy ; 00572 if( kk == 0 ) km = kz ; 00573 else if( kk == nz-1 ) kp = kz ; 00574 00575 for( jj=0 ; jj < ny ; jj++ ){ 00576 jy = jj*nx ; jm = jy-nx ; jp = jy+nx ; 00577 if( jj == 0 ) jm = jy ; 00578 else if( jj == ny-1 ) jp = jy ; 00579 00580 for( ii=0 ; ii < nx ; ii++ ){ 00581 if( mmm[ii+jy+kz] ){ /* count nonzero nbhrs */ 00582 im = ii-1 ; ip = ii+1 ; 00583 if( ii == 0 ) im = 0 ; 00584 else if( ii == nx-1 ) ip = ii ; 00585 num = mmm[im+jy+km] 00586 + mmm[ii+jm+km] + mmm[ii+jy+km] + mmm[ii+jp+km] 00587 + mmm[ip+jy+km] 00588 + mmm[im+jm+kz] + mmm[im+jy+kz] + mmm[im+jp+kz] 00589 + mmm[ii+jm+kz] + mmm[ii+jp+kz] 00590 + mmm[ip+jm+kz] + mmm[ip+jy+kz] + mmm[ip+jp+kz] 00591 + mmm[im+jy+kp] 00592 + mmm[ii+jm+kp] + mmm[ii+jy+kp] + mmm[ii+jp+kp] 00593 + mmm[ip+jy+kp] ; 00594 if( num < 17 ) nnn[ii+jy+kz] = 1 ; /* mark to erode */ 00595 } 00596 } } } 00597 00598 for( jj=ii=0 ; ii < nxyz ; ii++ ) /* actually erode */ 00599 if( nnn[ii] ){ mmm[ii] = 0 ; jj++ ; } 00600 00601 if( verb && jj > 0 ) ININFO_message("Eroded %d voxels\n",jj) ; 00602 00603 /* re-dilate eroded voxels that are next to survivors */ 00604 00605 #ifdef USE_DILATE 00606 for( kk=0 ; kk < nz ; kk++ ){ 00607 kz = kk*nxy ; km = kz-nxy ; kp = kz+nxy ; 00608 if( kk == 0 ) km = kz ; 00609 else if( kk == nz-1 ) kp = kz ; 00610 00611 for( jj=0 ; jj < ny ; jj++ ){ 00612 jy = jj*nx ; jm = jy-nx ; jp = jy+nx ; 00613 if( jj == 0 ) jm = jy ; 00614 else if( jj == ny-1 ) jp = jy ; 00615 00616 for( ii=0 ; ii < nx ; ii++ ){ 00617 if( nnn[ii+jy+kz] ){ /* was eroded */ 00618 im = ii-1 ; ip = ii+1 ; 00619 if( ii == 0 ) im = 0 ; 00620 else if( ii == nx-1 ) ip = ii ; 00621 nnn[ii+jy+kz] = /* see if has any nbhrs */ 00622 mmm[im+jy+km] 00623 || mmm[ii+jm+km] || mmm[ii+jy+km] || mmm[ii+jp+km] 00624 || mmm[ip+jy+km] 00625 || mmm[im+jm+kz] || mmm[im+jy+kz] || mmm[im+jp+kz] 00626 || mmm[ii+jm+kz] || mmm[ii+jp+kz] 00627 || mmm[ip+jm+kz] || mmm[ip+jy+kz] || mmm[ip+jp+kz] 00628 || mmm[im+jy+kp] 00629 || mmm[ii+jm+kp] || mmm[ii+jy+kp] || mmm[ii+jp+kp] 00630 || mmm[ip+jy+kp] ; 00631 } 00632 } } } 00633 00634 /* actually do the dilation */ 00635 00636 for( jj=ii=0 ; ii < nxyz ; ii++ ) 00637 if( nnn[ii] ){ mmm[ii] = 1 ; jj++ ; } 00638 00639 if( verb && jj > 0 ) ININFO_message("Restored %d eroded voxels\n",jj) ; 00640 #endif 00641 00642 free(nnn) ; return ; 00643 } |
|
Fill in a byte mask, repeatedly until it doesn't fill any more. Return value is number of voxels filled. ------------------------------------------------------------------------ Definition at line 431 of file thd_automask.c. References ENTRY, mmm, nz, RETURN, and THD_mask_fillin_once(). Referenced by main(), and mri_automask_image().
00432 { 00433 int nfill=0 , kfill ; 00434 00435 ENTRY("THD_mask_fillin_completely") ; 00436 00437 do{ 00438 kfill = THD_mask_fillin_once( nx,ny,nz , mmm , nside ) ; 00439 nfill += kfill ; 00440 } while( kfill > 0 ) ; 00441 00442 RETURN(nfill) ; 00443 } |
|
Fill in a byte mask. Filling is done by looking to each side (plus/minus) of a non-filled voxel, and seeing if there is a filled voxel on both sides. This looking is done parallel to the x-, y-, and z-axes, out to distance nside voxels.
Definition at line 277 of file thd_automask.c. References AFMALL, ENTRY, free, mmm, nz, and RETURN. Referenced by mri_automask_image(), mri_brainormalize(), SUMA_SurfGridIntersect(), and THD_mask_fillin_completely().
00278 { 00279 int ii,jj,kk , nsx,nsy,nsz , nxy,nxyz , iv,jv,kv,ll , nfill ; 00280 byte *nnn ; 00281 int nx2,nx3,nx4 , nxy2,nxy3,nxy4 ; 00282 00283 ENTRY("THD_mask_fillin_once") ; 00284 00285 if( mmm == NULL || nside <= 0 ) RETURN(0) ; 00286 00287 nsx = (nx-1)/2 ; if( nsx > nside ) nsx = nside ; 00288 nsy = (ny-1)/2 ; if( nsy > nside ) nsy = nside ; 00289 nsz = (nz-1)/2 ; if( nsz > nside ) nsz = nside ; 00290 00291 if( nsx == 0 && nsy == 0 && nsz == 0 ) RETURN(0) ; 00292 00293 #ifdef DEBUG 00294 fprintf(stderr,"THD_mask_fillin_once: nsx=%d nsy=%d nsz=%d\n",nsx,nsy,nsz); 00295 #endif 00296 00297 nxy = nx*ny ; nxyz = nxy*nz ; nfill = 0 ; 00298 00299 nx2 = 2*nx ; nx3 = 3*nx ; nx4 = 4*nx ; 00300 nxy2 = 2*nxy ; nxy3 = 3*nxy ; nxy4 = 4*nxy ; 00301 00302 nnn = AFMALL(byte, nxyz) ; /* stores filled in values */ 00303 00304 /* loop over voxels */ 00305 00306 #define FILLVOX \ 00307 do{ nnn[iv] = 1; nfill++; goto NextVox; } while(0) 00308 00309 for( kk=nsz ; kk < nz-nsz ; kk++ ){ 00310 kv = kk*nxy ; 00311 for( jj=nsy ; jj < ny-nsy ; jj++ ){ 00312 jv = jj*nx + kv ; 00313 for( ii=nsx ; ii < nx-nsx ; ii++ ){ 00314 iv = ii+jv ; 00315 if( mmm[iv] ) continue ; /* already filled */ 00316 00317 /* check in +x direction, then -x if +x hits */ 00318 00319 switch( nsx ){ 00320 case 1: 00321 if( mmm[iv+1] && mmm[iv-1] ) FILLVOX; 00322 break ; 00323 00324 case 2: 00325 if( (mmm[iv+1]||mmm[iv+2]) && 00326 (mmm[iv-1]||mmm[iv-2]) ) FILLVOX; 00327 break ; 00328 00329 case 3: 00330 if( (mmm[iv+1]||mmm[iv+2]||mmm[iv+3]) && 00331 (mmm[iv-1]||mmm[iv-2]||mmm[iv-3]) ) FILLVOX; 00332 break ; 00333 00334 case 4: 00335 if( (mmm[iv+1]||mmm[iv+2]||mmm[iv+3]||mmm[iv+4]) && 00336 (mmm[iv-1]||mmm[iv-2]||mmm[iv-3]||mmm[iv-4]) ) FILLVOX; 00337 break ; 00338 00339 default: 00340 for( ll=1 ; ll <= nsx ; ll++ ) if( mmm[iv+ll] ) break ; 00341 if( ll <= nsx ){ 00342 for( ll=1 ; ll <= nsx ; ll++ ) if( mmm[iv-ll] ) break ; 00343 if( ll <= nsx ) FILLVOX; 00344 } 00345 break ; 00346 } 00347 00348 /* check in +y direction, then -y if +y hits */ 00349 00350 switch( nsy ){ 00351 case 1: 00352 if( mmm[iv+nx] && mmm[iv-nx] ) FILLVOX; 00353 break ; 00354 00355 case 2: 00356 if( (mmm[iv+nx]||mmm[iv+nx2]) && 00357 (mmm[iv-nx]||mmm[iv-nx2]) ) FILLVOX; 00358 break ; 00359 00360 case 3: 00361 if( (mmm[iv+nx]||mmm[iv+nx2]||mmm[iv+nx3]) && 00362 (mmm[iv-nx]||mmm[iv-nx2]||mmm[iv-nx3]) ) FILLVOX; 00363 break ; 00364 00365 case 4: 00366 if( (mmm[iv+nx]||mmm[iv+nx2]||mmm[iv+nx3]||mmm[iv+nx4]) && 00367 (mmm[iv-nx]||mmm[iv-nx2]||mmm[iv-nx3]||mmm[iv-nx4]) ) FILLVOX; 00368 break ; 00369 00370 default: 00371 for( ll=1 ; ll <= nsy ; ll++ ) if( mmm[iv+ll*nx] ) break ; 00372 if( ll <= nsy ){ 00373 for( ll=1 ; ll <= nsy ; ll++ ) if( mmm[iv-ll*nx] ) break ; 00374 if( ll <= nsy ) FILLVOX; 00375 } 00376 break ; 00377 } 00378 00379 /* check in +z direction, then -z if +z hits */ 00380 00381 switch( nsz ){ 00382 case 1: 00383 if( mmm[iv+nxy] && mmm[iv-nxy] ) FILLVOX; 00384 break ; 00385 00386 case 2: 00387 if( (mmm[iv+nxy]||mmm[iv+nxy2]) && 00388 (mmm[iv-nxy]||mmm[iv-nxy2]) ) FILLVOX; 00389 break ; 00390 00391 case 3: 00392 if( (mmm[iv+nxy]||mmm[iv+nxy2]||mmm[iv+nxy3]) && 00393 (mmm[iv-nxy]||mmm[iv-nxy2]||mmm[iv-nxy3]) ) FILLVOX; 00394 break ; 00395 00396 case 4: 00397 if( (mmm[iv+nxy]||mmm[iv+nxy2]||mmm[iv+nxy3]||mmm[iv+nxy4]) && 00398 (mmm[iv-nxy]||mmm[iv-nxy2]||mmm[iv-nxy3]||mmm[iv-nxy4]) ) FILLVOX; 00399 break ; 00400 00401 default: 00402 for( ll=1 ; ll <= nsz ; ll++ ) if( mmm[iv+ll*nxy] ) break ; 00403 if( ll <= nsz ){ 00404 for( ll=1 ; ll <= nsz ; ll++ ) if( mmm[iv-ll*nxy] ) break ; 00405 if( ll <= nsz ) FILLVOX; 00406 } 00407 break ; 00408 } 00409 00410 NextVox: ; /* end of loop over ii */ 00411 } } } 00412 00413 /* copy fills back into mmm */ 00414 00415 if( nfill > 0 ){ 00416 for( iv=0 ; iv < nxyz ; iv++ ) if( nnn[iv] ) mmm[iv] = 1 ; 00417 } 00418 00419 #ifdef DEBUG 00420 fprintf(stderr,"THD_mask_fillin_once: nfill=%d\n",nfill) ; 00421 #endif 00422 00423 free(nnn) ; RETURN(nfill) ; 00424 } |
|
Definition at line 817 of file thd_automask.c. Referenced by mri_automask_image().
00818 { 00819 int nxy=nx*ny , ii,jj,kk , ijk , bot,top , pd=pdepth ; 00820 int nxp=nx-pd , nyp=ny-pd , nzp=nz-pd ; 00821 int num=0 , dnum , nite ; 00822 00823 for( nite=0 ; nite < pd ; nite++ ){ 00824 dnum = 0 ; 00825 00826 for( kk=0 ; kk < nz ; kk++ ){ 00827 for( jj=0 ; jj < ny ; jj++ ){ 00828 ijk = jj*nx + kk*nxy ; 00829 for( bot=0 ; bot < nx && !mmm[bot+ijk]; bot++ ) ; 00830 top = bot+pd ; if( top >= nx ) continue ; 00831 for( ii=bot+1 ; ii <= top && mmm[ii+ijk] ; ii++ ) ; 00832 if( ii <= top ){ mmm[bot+ijk] = 0; dnum++; } 00833 }} 00834 00835 for( kk=0 ; kk < nz ; kk++ ){ 00836 for( jj=0 ; jj < ny ; jj++ ){ 00837 ijk = jj*nx + kk*nxy ; 00838 for( top=nx-1 ; top >= 0 && !mmm[top+ijk]; top-- ) ; 00839 bot = top-pd ; if( bot < 0 ) continue ; 00840 for( ii=top-1 ; ii >= bot && mmm[ii+ijk] ; ii-- ) ; 00841 if( ii >= bot ){ mmm[top+ijk] = 0; dnum++; } 00842 }} 00843 00844 for( kk=0 ; kk < nz ; kk++ ){ 00845 for( ii=0 ; ii < nx ; ii++ ){ 00846 ijk = ii + kk*nxy ; 00847 for( bot=0 ; bot < ny && !mmm[bot*nx+ijk] ; bot++ ) ; 00848 top = bot+pd ; 00849 if( top >= ny ) continue ; 00850 for( jj=bot+1 ; jj <= top && mmm[jj*nx+ijk] ; jj++ ) ; 00851 if( jj <= top ){ mmm[bot*nx+ijk] = 0; dnum++; } 00852 }} 00853 00854 for( kk=0 ; kk < nz ; kk++ ){ 00855 for( ii=0 ; ii < nx ; ii++ ){ 00856 ijk = ii + kk*nxy ; 00857 for( top=ny-1 ; top >= 0 && !mmm[top*nx+ijk] ; top-- ) ; 00858 bot = top-pd ; if( bot < 0 ) continue ; 00859 for( jj=top-1 ; jj >= bot && mmm[jj*nx+ijk] ; jj-- ) ; 00860 if( jj >= bot ){ mmm[top*nx+ijk] = 0; dnum++; } 00861 }} 00862 00863 for( jj=0 ; jj < ny ; jj++ ){ 00864 for( ii=0 ; ii < nx ; ii++ ){ 00865 ijk = ii + jj*nx ; 00866 for( top=nz-1 ; top >= 0 && !mmm[top*nxy+ijk] ; top-- ) ; 00867 bot = top-pd ; if( bot < 0 ) continue ; 00868 for( kk=top-1 ; kk >= bot && mmm[kk*nxy+ijk] ; kk-- ) ; 00869 if( kk >= bot ){ mmm[top*nxy+ijk] = 0; dnum++; } 00870 }} 00871 00872 num += dnum ; 00873 if( dnum == 0 ) break ; 00874 } 00875 00876 return num ; 00877 } |
Variable Documentation
|
Definition at line 14 of file thd_automask.c. Referenced by mri_automask_image(), and THD_mask_clust(). |
|
Definition at line 11 of file thd_automask.c. Referenced by mri_automask_image(), and THD_automask_extclip(). |
|
Definition at line 8 of file thd_automask.c. Referenced by mri_automask_image(), THD_automask_verbose(), and THD_mask_erode(). |