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(). |
