Doxygen Source Code Documentation
thd_brainormalize.c File Reference
#include "mrilib.h"
#include "thd_brainormalize.h"
Go to the source code of this file.
Data Structures | |
struct | basin |
struct | clipvec |
struct | shortvox |
Defines | |
#define | IJK(i, j, k) ((i)+(j)*nx+(k)*nxy) |
#define | DALL 4096 |
#define | DPUT(i, j, k, d) |
#define | CPUT(i, j, k) |
#define | DBALL 32768 |
#define | BDEP(i) (baslist[i]->depth) |
#define | INIT_BASIN(iv) |
#define | KILL_BASIN(ib) |
#define | ADDTO_BASIN(ib, iv) |
#define | MERGE_BASIN(ib, ic) |
#define | BASECHECK(a, b, c) |
Functions | |
void | mri_brainormalize_verbose (int v) |
void | mri_brainormalize_initialize (float dx, float dy, float dz) |
float | THD_BN_dxyz () |
int | THD_BN_nx () |
int | THD_BN_ny () |
int | THD_BN_nz () |
int | mask_count (int nvox, byte *mmm) |
short * | THD_mask_distize (int nx, int ny, int nz, byte *mmm, byte *ccc) |
void | clustedit3D (int nx, int ny, int nz, byte *mmm, int csize) |
float | partial_cliplevel (MRI_IMAGE *im, float mfrac, int ibot, int itop, int jbot, int jtop, int kbot, int ktop) |
clipvec | get_octant_clips (MRI_IMAGE *im, float mfrac) |
INLINE float | pointclip (int ii, int jj, int kk, clipvec *cv) |
byte * | mri_short2mask (MRI_IMAGE *im) |
sort_shortvox (int n, shortvox *ar, int dec, float botperc, float topperc) | |
MRI_IMAGE * | mri_watershedize (MRI_IMAGE *sim, float prefac) |
byte * | make_peel_mask (int nx, int ny, int nz, byte *mmm, int pdepth) |
void | zedit_mask (int nx, int ny, int nz, byte *mmm, int zdepth, int zbot) |
void | ijkwarp (float i, float j, float k, float *x, float *y, float *z) |
void | ijk_invwarp (float x, float y, float z, float *i, float *j, float *k) |
void | brainnormalize_coord (float ispat, float jspat, float kspat, float *iorig, float *jorig, float *korig, THD_3dim_dataset *origset, float *xrai_orig, float *yrai_orig, float *zrai_orig) |
takes in voxel indices into the Spat Normed volume (RAI) and returns voxel indices and coordinates in the original volume. Used to figure out shift to apply to surface model to align it with original volume. | |
MRI_IMAGE * | mri_brainormalize (MRI_IMAGE *im, int xxor, int yyor, int zzor, MRI_IMAGE **imout_origp, MRI_IMAGE **imout_edge) |
Variables | |
int | verb = 0 |
float | thd_bn_dxyz = 1.0 |
int | thd_bn_nx = 167 |
int | thd_bn_ny = 212 |
int | thd_bn_nz = 175 |
float | ai |
float | bi |
float | aj |
float | bj |
float | ak |
float | bk |
Define Documentation
|
Value: { register basin *bb = baslist[ib] ; \ if( bb->num == bb->nall ){ \ bb->nall = (int)(1.2*bb->nall)+32 ; \ bb->ivox = (int *)realloc( (void *)bb->ivox , \ sizeof(int)*bb->nall ) ; \ } \ bb->ivox[bb->num] = (iv) ; bb->num++ ; \ svox[iv].basin = (ib) ; } Definition at line 649 of file thd_brainormalize.c. Referenced by mri_watershedize(). |
|
Value: |
|
Definition at line 622 of file thd_brainormalize.c. Referenced by mri_watershedize(). |
|
Value: do{ ijk = (i)+(j)*nx+(k)*nxy ; \ if( mmm[ijk] ){ \ if( nnow == nall ){ \ nall += DALL ; \ inow = (short *) realloc((void *)inow,sizeof(short)*nall) ; \ jnow = (short *) realloc((void *)jnow,sizeof(short)*nall) ; \ know = (short *) realloc((void *)know,sizeof(short)*nall) ; \ } \ inow[nnow] = (i); jnow[nnow] = (j); know[nnow] = (k); \ nnow++ ; mmm[ijk] = 0 ; \ } } while(0) Definition at line 134 of file thd_brainormalize.c. |
|
Definition at line 53 of file thd_brainormalize.c. Referenced by THD_mask_distize(). |
|
Definition at line 620 of file thd_brainormalize.c. Referenced by mri_watershedize(). |
|
Value: do{ ijk = (i)+(j)*nx+(k)*nxy ; \ if( mmm[ijk] == 0 ) break ; \ if( nnow == nall ){ \ nall += DALL ; \ inow = (short *) realloc((void *)inow,sizeof(short)*nall) ; \ jnow = (short *) realloc((void *)jnow,sizeof(short)*nall) ; \ know = (short *) realloc((void *)know,sizeof(short)*nall) ; \ } \ inow[nnow] = (i); jnow[nnow] = (j); know[nnow] = (k); \ nnow++ ; mmm[ijk] = 0 ; ddd[ijk] = (d) ; \ } while(0) Definition at line 59 of file thd_brainormalize.c. Referenced by THD_mask_distize(). |
|
Definition at line 5 of file thd_brainormalize.c. Referenced by clustedit3D(), mri_watershedize(), partial_cliplevel(), and SUMA_process_NIML_data(). |
|
Value: { register int qb=nbtop; \ if( qb >= nball ){ \ register int qqb=nball+DBALL,zb ; \ baslist = (basin **)realloc((void *)baslist, \ sizeof(basin *)*qqb) ; \ for( zb=nball ; zb < qqb ; zb++ ) baslist[zb] = NULL ; \ nball = qqb ; \ } \ baslist[qb] = (basin *) malloc(sizeof(basin)) ; \ baslist[qb]->num = 1 ; \ baslist[qb]->nall = 1 ; \ baslist[qb]->depth = svox[iv].val ; \ baslist[qb]->ivox = (int *)malloc(sizeof(int)) ; \ baslist[qb]->ivox[0] = (iv) ; \ svox[iv].basin = qb ; nbtop++ ; \ } Definition at line 624 of file thd_brainormalize.c. Referenced by mri_watershedize(). |
|
Value: { if( baslist[ib] != NULL ){ \ free((void *)baslist[ib]->ivox) ; \ free((void *)baslist[ib]) ; \ baslist[ib] = NULL ; } \ } Definition at line 642 of file thd_brainormalize.c. Referenced by mri_watershedize(). |
|
Value: { register basin *bb = baslist[ib], *cc = baslist[ic] ; \ int zz = bb->num + cc->num ; \ if( bb->nall < zz ){ \ bb->nall = zz+1 ; \ bb->ivox = (int *)realloc( (void *)bb->ivox , \ sizeof(int)*bb->nall ) ; \ } \ memcpy( bb->ivox + bb->num , \ cc->ivox , sizeof(int) * cc->num ) ; \ bb->num = zz ; \ for( zz=0 ; zz < cc->num ; zz++ ) \ svox[cc->ivox[zz]].basin = (ib) ; \ KILL_BASIN(ic) ; } Definition at line 659 of file thd_brainormalize.c. Referenced by mri_watershedize(). |
Function Documentation
|
takes in voxel indices into the Spat Normed volume (RAI) and returns voxel indices and coordinates in the original volume. Used to figure out shift to apply to surface model to align it with original volume.
Definition at line 1023 of file thd_brainormalize.c. References CURRENT_DAXES, THD_3dim_dataset::daxes, ijkwarp(), THD_dataxes::nxx, THD_dataxes::nyy, THD_dataxes::nzz, ORI_A2P_TYPE, ORI_I2S_TYPE, ORI_L2R_TYPE, ORI_P2A_TYPE, ORI_R2L_TYPE, ORI_S2I_TYPE, THD_3dmm_to_dicomm(), THD_BN_XORG, THD_BN_YORG, THD_BN_ZORG, THD_dataxes::xxdel, THD_dataxes::xxorg, THD_dataxes::xxorient, THD_fvec3::xyz, THD_dataxes::yydel, THD_dataxes::yyorg, THD_dataxes::yyorient, THD_dataxes::zzdel, THD_dataxes::zzorg, and THD_dataxes::zzorient.
01027 { 01028 THD_dataxes * daxes ; 01029 THD_fvec3 fv, fvdic ; 01030 float irai, jrai, krai; 01031 01032 /* find out corresponding indices in original dset */ 01033 ijkwarp(ispat, jspat, kspat , &irai, &jrai, &krai); 01034 01035 01036 01037 /* These indices assume an RAI dset orientation */ 01038 /* Find out what these indices should be in origset's orientation */ 01039 switch( origset->daxes->xxorient ){ 01040 case ORI_R2L_TYPE: *iorig = irai ; break ; 01041 case ORI_L2R_TYPE: *iorig = origset->daxes->nxx - irai ; break ; 01042 case ORI_P2A_TYPE: *iorig = origset->daxes->nxx - jrai ; break ; 01043 case ORI_A2P_TYPE: *iorig = jrai ; break ; 01044 case ORI_I2S_TYPE: *iorig = krai ; break ; 01045 case ORI_S2I_TYPE: *iorig = origset->daxes->nxx - krai ; break ; 01046 } 01047 switch( origset->daxes->yyorient ){ 01048 case ORI_R2L_TYPE: *jorig = irai ; break ; 01049 case ORI_L2R_TYPE: *jorig = origset->daxes->nyy - irai ; break ; 01050 case ORI_P2A_TYPE: *jorig = origset->daxes->nyy - jrai ; break ; 01051 case ORI_A2P_TYPE: *jorig = jrai ; break ; 01052 case ORI_I2S_TYPE: *jorig = krai ; break ; 01053 case ORI_S2I_TYPE: *jorig = origset->daxes->nyy - krai ; break ; 01054 } 01055 switch( origset->daxes->zzorient ){ 01056 case ORI_R2L_TYPE: *korig = irai ; break ; 01057 case ORI_L2R_TYPE: *korig = origset->daxes->nzz - irai ; break ; 01058 case ORI_P2A_TYPE: *korig = origset->daxes->nzz - jrai ; break ; 01059 case ORI_A2P_TYPE: *korig = jrai ; break ; 01060 case ORI_I2S_TYPE: *korig = krai ; break ; 01061 case ORI_S2I_TYPE: *korig = origset->daxes->nzz - krai ; break ; 01062 } 01063 01064 01065 /* change indices into mm coords in orig dset*/ 01066 daxes = CURRENT_DAXES(origset) ; 01067 01068 fv.xyz[0] = daxes->xxorg + *iorig * daxes->xxdel ; /* 3dfind_to_3dmm */ 01069 fv.xyz[1] = daxes->yyorg + *jorig * daxes->yydel ; 01070 fv.xyz[2] = daxes->zzorg + *korig * daxes->zzdel ; 01071 01072 fvdic = THD_3dmm_to_dicomm(origset, fv ); /* 3dmm_to_dicomm */ 01073 *xrai_orig = fvdic.xyz[0]; 01074 *yrai_orig = fvdic.xyz[1]; 01075 *zrai_orig = fvdic.xyz[2]; 01076 01077 /* report for sanity */ 01078 if (0) { 01079 fprintf(stderr, "brainnormalize_coord:\n" 01080 " ijk_spat_rai = [%f %f %f]\n" 01081 " ijk_orig_rai = [%f %f %f] (in rai order, not native to iset!)\n" 01082 " ijk_orig = [%f %f %f] (in native order)\n" 01083 " XYZ_orig = [%f %f %f]\n" 01084 " Origin spat = [%f %f %f]\n", 01085 ispat, jspat, kspat, 01086 irai, jrai, krai , 01087 *iorig, *jorig, *korig , 01088 *xrai_orig, *yrai_orig, *zrai_orig, 01089 THD_BN_XORG, THD_BN_YORG, THD_BN_ZORG); 01090 } 01091 return; 01092 } |
|
Remove clusters below csize from the binary mask mmm. Definition at line 150 of file thd_brainormalize.c. References CPUT, free, IJK, malloc, mmm, nz, and realloc. Referenced by mri_short2mask().
00151 { 00152 int ii,jj,kk, icl , nxy,nxyz , ijk , ijk_last ; 00153 int ip,jp,kp , im,jm,km ; 00154 int nnow , nall , nsav , nkill ; 00155 short *inow , *jnow , *know ; 00156 short *isav , *jsav , *ksav ; 00157 00158 if( nx < 1 || ny < 1 || nz < 1 || mmm == NULL || csize < 2 ) return ; 00159 00160 nxy = nx*ny ; nxyz = nxy*nz ; 00161 00162 nall = 8 ; /* # allocated pts */ 00163 inow = (short *) malloc(sizeof(short)*nall) ; /* coords of pts */ 00164 jnow = (short *) malloc(sizeof(short)*nall) ; 00165 know = (short *) malloc(sizeof(short)*nall) ; 00166 00167 nsav = nkill = 0 ; 00168 isav = (short *) malloc(sizeof(short)) ; 00169 jsav = (short *) malloc(sizeof(short)) ; 00170 ksav = (short *) malloc(sizeof(short)) ; 00171 00172 /*--- scan through array, find nonzero point, build a cluster, ... ---*/ 00173 00174 #if 0 00175 if( verb ) fprintf(stderr," + clustedit3D: threshold size=%d voxels\n",csize) ; 00176 #endif 00177 00178 ijk_last = 0 ; 00179 while(1) { 00180 /* find next nonzero point */ 00181 00182 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( mmm[ijk] ) break ; 00183 if( ijk == nxyz ) break ; /* didn't find any! */ 00184 00185 ijk_last = ijk+1 ; /* start here next time */ 00186 00187 /* init current cluster list with this point */ 00188 00189 mmm[ijk] = 0 ; /* clear found point */ 00190 nnow = 1 ; /* # pts in cluster */ 00191 inow[0] = ijk % nx ; 00192 jnow[0] = (ijk%nxy)/nx ; 00193 know[0] = ijk / nxy ; 00194 00195 /*-- 00196 for each point in cluster: 00197 check 6 neighboring points for nonzero entries in mmm 00198 enter those into cluster (and clear them in mmm) 00199 continue until end of cluster is reached 00200 (note that cluster size nnow is expanding as we progress) 00201 --*/ 00202 00203 for( icl=0 ; icl < nnow ; icl++ ){ 00204 ii = inow[icl] ; jj = jnow[icl] ; kk = know[icl] ; 00205 im = ii-1 ; jm = jj-1 ; km = kk-1 ; 00206 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ; 00207 00208 if( im >= 0 ) CPUT(im,jj,kk) ; 00209 if( ip < nx ) CPUT(ip,jj,kk) ; 00210 if( jm >= 0 ) CPUT(ii,jm,kk) ; 00211 if( jp < ny ) CPUT(ii,jp,kk) ; 00212 if( km >= 0 ) CPUT(ii,jj,km) ; 00213 if( kp < nz ) CPUT(ii,jj,kp) ; 00214 } 00215 00216 /* see if now cluster is large enough to live */ 00217 00218 if( nnow >= csize ){ 00219 kk = nsav+nnow ; 00220 isav = (short *)realloc((void *)isav,sizeof(short)*kk) ; 00221 jsav = (short *)realloc((void *)jsav,sizeof(short)*kk) ; 00222 ksav = (short *)realloc((void *)ksav,sizeof(short)*kk) ; 00223 memcpy(isav+nsav,inow,sizeof(short)*nnow) ; 00224 memcpy(jsav+nsav,jnow,sizeof(short)*nnow) ; 00225 memcpy(ksav+nsav,know,sizeof(short)*nnow) ; 00226 nsav = kk ; 00227 #if 0 00228 if( verb ) 00229 fprintf(stderr," + clustedit3D: saved cluster with %d voxels\n",nnow); 00230 #endif 00231 } else { 00232 nkill += nnow ; 00233 } 00234 00235 } /* loop ends when all nonzero points are clustered */ 00236 00237 free((void *)inow); free((void *)jnow); free((void *)know); 00238 00239 /* copy saved points back into mmm */ 00240 00241 for( ii=0 ; ii < nsav ; ii++ ) 00242 mmm[ IJK(isav[ii],jsav[ii],ksav[ii]) ] = 1 ; 00243 00244 free((void *)isav); free((void *)jsav); free((void *)ksav) ; 00245 00246 #if 0 00247 if( verb ) 00248 fprintf(stderr," + clustedit3D totals:" 00249 " saved=%d killed=%d nxyz=%d\n",nsav,nkill,nxyz) ; 00250 #endif 00251 return ; 00252 } |
|
Get the cliplevel for each octant about the center-of-mass. Definition at line 342 of file thd_brainormalize.c. References clipvec::clip_000, clipvec::clip_001, clipvec::clip_010, clipvec::clip_011, clipvec::clip_100, clipvec::clip_101, clipvec::clip_110, clipvec::clip_111, clipvec::clip_max, clipvec::clip_min, clipvec::dxi, clipvec::dyi, clipvec::dzi, ENTRY, MRI_IMAGE::kind, MAX, MIN, MRI_SHORT_PTR, MRI_IMAGE::nx, MRI_IMAGE::ny, MRI_IMAGE::nz, nz, partial_cliplevel(), RETURN, clipvec::x0, clipvec::x1, clipvec::y0, clipvec::y1, clipvec::z0, and clipvec::z1. Referenced by mri_brainormalize(), and mri_short2mask().
00343 { 00344 float xcm,ycm,zcm , sum,val , clip_min , clip_max ; 00345 int ii,jj,kk , nx,ny,nz,nxy , ic,jc,kc , it,jt,kt , ijk ; 00346 short *sar ; 00347 clipvec cv ; 00348 00349 ENTRY("get_octant_clips") ; 00350 00351 cv.clip_000 = -1 ; /* flags error return */ 00352 00353 if( im == NULL || im->kind != MRI_short ) RETURN(cv) ; 00354 00355 nx = im->nx ; ny = im->ny ; nz = im->nz ; nxy = nx*ny ; 00356 it = nx-1 ; jt = ny-1 ; kt = nz-1 ; 00357 00358 /* compute CM of image */ 00359 00360 sar = MRI_SHORT_PTR(im) ; if( sar == NULL ) RETURN(cv) ; 00361 00362 xcm = ycm = zcm = sum = 0.0 ; 00363 for( ijk=kk=0 ; kk < nz ; kk++ ){ 00364 for( jj=0 ; jj < ny ; jj++ ){ 00365 for( ii=0 ; ii < nx ; ii++,ijk++ ){ 00366 val = (float)sar[ijk]; if( val <= 0.0 ) continue; 00367 sum += val ; 00368 xcm += val * ii ; 00369 ycm += val * jj ; 00370 zcm += val * kk ; 00371 }}} 00372 if( sum == 0.0 ) RETURN(cv) ; 00373 ic = (int)rint(xcm/sum); jc = (int)rint(ycm/sum); kc = (int)rint(zcm/sum); 00374 00375 /* compute cliplevel in each octant about the CM */ 00376 00377 val = 0.5 * partial_cliplevel( im,mfrac , 0,it , 0,jt , 0,kt ) ; 00378 00379 cv.clip_000 = partial_cliplevel( im,mfrac, 0 ,ic+2, 0 ,jc+2, 0 ,kc+2 ); 00380 cv.clip_100 = partial_cliplevel( im,mfrac, ic-2,it , 0 ,jc+2, 0 ,kc+2 ); 00381 cv.clip_010 = partial_cliplevel( im,mfrac, 0 ,ic+2, jc-2,jt , 0 ,kc+2 ); 00382 cv.clip_110 = partial_cliplevel( im,mfrac, ic-2,it , jc-2,jt , 0 ,kc+2 ); 00383 cv.clip_001 = partial_cliplevel( im,mfrac, 0 ,ic+2, 0 ,jc+2, kc-2,kt ); 00384 cv.clip_101 = partial_cliplevel( im,mfrac, ic-2,it , 0 ,jc+2, kc-2,kt ); 00385 cv.clip_011 = partial_cliplevel( im,mfrac, 0 ,ic+2, jc-2,jt , kc-2,kt ); 00386 cv.clip_111 = partial_cliplevel( im,mfrac, ic-2,it , jc-2,jt , kc-2,kt ); 00387 00388 if( cv.clip_000 < val ) cv.clip_000 = val ; /* don't let */ 00389 if( cv.clip_100 < val ) cv.clip_100 = val ; /* clip levels */ 00390 if( cv.clip_010 < val ) cv.clip_010 = val ; /* get too */ 00391 if( cv.clip_110 < val ) cv.clip_110 = val ; /* small! */ 00392 if( cv.clip_001 < val ) cv.clip_001 = val ; 00393 if( cv.clip_101 < val ) cv.clip_101 = val ; 00394 if( cv.clip_011 < val ) cv.clip_011 = val ; 00395 if( cv.clip_111 < val ) cv.clip_111 = val ; 00396 00397 clip_min = cv.clip_000 ; 00398 clip_min = MIN(clip_min,cv.clip_100) ; 00399 clip_min = MIN(clip_min,cv.clip_010) ; 00400 clip_min = MIN(clip_min,cv.clip_110) ; 00401 clip_min = MIN(clip_min,cv.clip_001) ; 00402 clip_min = MIN(clip_min,cv.clip_101) ; 00403 clip_min = MIN(clip_min,cv.clip_011) ; 00404 clip_min = MIN(clip_min,cv.clip_111) ; cv.clip_min = clip_min ; 00405 00406 clip_max = cv.clip_000 ; 00407 clip_max = MAX(clip_max,cv.clip_100) ; 00408 clip_max = MAX(clip_max,cv.clip_010) ; 00409 clip_max = MAX(clip_max,cv.clip_110) ; 00410 clip_max = MAX(clip_max,cv.clip_001) ; 00411 clip_max = MAX(clip_max,cv.clip_101) ; 00412 clip_max = MAX(clip_max,cv.clip_011) ; 00413 clip_max = MAX(clip_max,cv.clip_111) ; cv.clip_max = clip_max ; 00414 00415 /* (x0,y0,z0) = center of lowest octant 00416 (x1,y1,z1) = center of highest octant */ 00417 00418 cv.x0 = 0.5*ic ; cv.x1 = 0.5*(ic+it) ; 00419 cv.y0 = 0.5*jc ; cv.y1 = 0.5*(jc+jt) ; 00420 cv.z0 = 0.5*kc ; cv.z1 = 0.5*(kc+kt) ; 00421 cv.dxi = (cv.x1 > cv.x0) ? 1.0/(cv.x1-cv.x0) : 0.0 ; 00422 cv.dyi = (cv.y1 > cv.y0) ? 1.0/(cv.y1-cv.y0) : 0.0 ; 00423 cv.dzi = (cv.z1 > cv.z0) ? 1.0/(cv.z1-cv.z0) : 0.0 ; 00424 00425 #if 0 00426 if( verb ) 00427 fprintf(stderr," + get_octant_clips: min clip=%.1f\n" 00428 " clip_000=%.1f clip_100=%.1f clip_010=%.1f clip_110=%.1f\n" 00429 " clip_001=%.1f clip_101=%.1f clip_011=%.1f clip_111=%.1f\n" 00430 " (x0,y0,z0)=(%.1f,%.1f,%.1f) (x1,y1,z1)=(%.1f,%.1f,%.1f)\n" , 00431 val , 00432 cv.clip_000 , cv.clip_100 , cv.clip_010 , cv.clip_110 , 00433 cv.clip_001 , cv.clip_101 , cv.clip_011 , cv.clip_111 , 00434 cv.x0 , cv.y0 , cv.z0 , cv.x1 , cv.y1 , cv.z1 ) ; 00435 #endif 00436 00437 RETURN(cv) ; 00438 } |
|
Definition at line 994 of file thd_brainormalize.c. References ai, aj, ak, bi, bj, bk, and i. Referenced by mri_brainormalize().
|
|
Definition at line 986 of file thd_brainormalize.c. References ai, aj, ak, bi, bj, bk, and i. Referenced by brainnormalize_coord(), and mri_brainormalize().
|
|
Definition at line 864 of file thd_brainormalize.c. References calloc, free, mask_count(), mmm, nz, THD_mask_clust(), THD_mask_erode(), and top.
00865 { 00866 int nxy=nx*ny,nxyz=nxy*nz , ii,jj,kk , ijk , bot,top , pd=pdepth ; 00867 byte *ppp ; 00868 00869 if( mmm == NULL || pdepth <= 1 ) return NULL ; 00870 00871 ppp = (byte *)calloc(sizeof(byte),nxyz) ; 00872 00873 for( kk=0 ; kk < nz ; kk++ ){ 00874 for( jj=0 ; jj < ny ; jj++ ){ 00875 ijk = jj*nx + kk*nxy ; 00876 for( bot=0 ; bot < nx && !mmm[bot+ijk]; bot++ ) ; 00877 top = bot+pd ; if( top >= nx ) continue ; 00878 for( ii=bot+1 ; ii <= top && mmm[ii+ijk] ; ii++ ) ; 00879 if( ii <= top ){ 00880 top = ii; for( ii=bot ; ii <= top ; ii++ ) ppp[ii+ijk] = 1; 00881 } 00882 }} 00883 00884 for( kk=0 ; kk < nz ; kk++ ){ 00885 for( jj=0 ; jj < ny ; jj++ ){ 00886 ijk = jj*nx + kk*nxy ; 00887 for( top=nx-1 ; top >= 0 && !mmm[top+ijk]; top-- ) ; 00888 bot = top-pd ; if( bot < 0 ) continue ; 00889 for( ii=top-1 ; ii >= bot && mmm[ii+ijk] ; ii-- ) ; 00890 if( ii >= bot ){ 00891 bot = ii; for( ii=top ; ii >= bot ; ii-- ) ppp[ii+ijk] = 1; 00892 } 00893 }} 00894 00895 for( kk=0 ; kk < nz ; kk++ ){ 00896 for( ii=0 ; ii < nx ; ii++ ){ 00897 ijk = ii + kk*nxy ; 00898 for( bot=0 ; bot < ny && !mmm[bot*nx+ijk] ; bot++ ) ; 00899 top = bot+pd ; 00900 if( top >= ny ) continue ; 00901 for( jj=bot+1 ; jj <= top && mmm[jj*nx+ijk] ; jj++ ) ; 00902 if( jj <= top ){ 00903 top = jj; for( jj=bot ; jj <= top ; jj++ ) ppp[jj*nx+ijk] = 1; 00904 } 00905 }} 00906 00907 for( kk=0 ; kk < nz ; kk++ ){ 00908 for( ii=0 ; ii < nx ; ii++ ){ 00909 ijk = ii + kk*nxy ; 00910 for( top=ny-1 ; top >= 0 && !mmm[top*nx+ijk] ; top-- ) ; 00911 bot = top-pd ; if( bot < 0 ) continue ; 00912 for( jj=top-1 ; jj >= bot && mmm[jj*nx+ijk] ; jj-- ) ; 00913 if( jj >= bot ){ 00914 bot = jj; for( jj=top ; jj >= bot ; jj-- ) ppp[jj*nx+ijk] = 1; 00915 } 00916 }} 00917 00918 #if 1 00919 for( jj=0 ; jj < ny ; jj++ ){ 00920 for( ii=0 ; ii < nx ; ii++ ){ 00921 ijk = ii + jj*nx ; 00922 for( top=nz-1 ; top >= 0 && !mmm[top*nxy+ijk] ; top-- ) ; 00923 bot = top-pd ; if( bot < 0 ) continue ; 00924 for( kk=top-1 ; kk >= bot && mmm[kk*nxy+ijk] ; kk-- ) ; 00925 if( kk >= bot ){ 00926 bot = kk; for( kk=top ; kk >= bot ; kk-- ) ppp[kk*nxy+ijk] = 1; 00927 } 00928 }} 00929 #endif 00930 00931 kk = mask_count(nxyz,ppp) ; 00932 if( kk == 0 ){ free((void *)ppp) ; return NULL ; } 00933 if( verb ) fprintf(stderr," + Initial peel mask has %d voxels\n",kk ) ; 00934 THD_mask_erode( nx,ny,nz, ppp ) ; 00935 THD_mask_clust( nx,ny,nz, ppp ) ; 00936 kk = mask_count(nxyz,ppp) ; 00937 if( kk == 0 ){ free((void *)ppp) ; return NULL ; } 00938 if( verb ) fprintf(stderr," + Final peel mask has %d voxels\n",kk ) ; 00939 00940 return ppp ; 00941 } |
|
Count number of nonzeros in mask array Definition at line 43 of file thd_brainormalize.c. References mmm. Referenced by make_peel_mask(), mri_automask_image(), mri_brainormalize(), mri_short2mask(), and zedit_mask().
00044 { 00045 int ii , nn ; 00046 for( nn=ii=0 ; ii < nvox ; ii++ ) nn += (mmm[ii] != 0) ; 00047 return nn ; 00048 } |
|
Definition at line 1106 of file thd_brainormalize.c. References abs, AFNI_noenv(), AFNI_yesenv(), ai, aj, ak, bi, bj, bk, calloc, clipvec::clip_000, MRI_IMAGE::dx, MRI_IMAGE::dy, MRI_IMAGE::dz, ENTRY, free, get_octant_clips(), ijk_invwarp(), ijkwarp(), MRI_IMAGE::kind, LAST_ORIENT_TYPE, malloc, mask_count(), MRI_BYTE_PTR, MRI_COPY_AUX, MRI_CUBIC, mri_flip3D(), mri_free(), mri_maxabs(), mri_new_conforming, MRI_NN, mri_short2mask(), MRI_SHORT_PTR, mri_to_short(), mri_warp3D(), mri_warp3D_method(), MRI_IMAGE::nx, MRI_IMAGE::ny, MRI_IMAGE::nz, nz, ORI_A2P_TYPE, ORI_I2S_TYPE, ORI_L2R_TYPE, ORI_P2A_TYPE, ORI_R2L_TYPE, ORI_S2I_TYPE, pointclip(), RETURN, SHORTIZE, thd_bn_dxyz, thd_bn_nx, thd_bn_ny, thd_bn_nz, THD_BN_XCM, THD_BN_XORG, THD_BN_YCM, THD_BN_YORG, THD_BN_ZCM, THD_BN_ZHEIGHT, THD_BN_ZORG, THD_mask_clust(), THD_mask_dilate(), THD_mask_distize(), THD_mask_erode(), THD_mask_fillin_once(), verb, MRI_IMAGE::xo, MRI_IMAGE::yo, z1, and MRI_IMAGE::zo. Referenced by main().
01107 { 01108 MRI_IMAGE *sim , *tim , *bim ; 01109 short *sar , sval ; 01110 int ii,jj,kk,ijk,ktop,kbot , nx,ny,nz,nxy,nxyz ; 01111 float val , icm,jcm,kcm,sum , dx,dy,dz ; 01112 byte *mask , *bar ; 01113 float sim_dx, sim_dy, sim_dz, sim_xo, sim_yo, sim_zo; 01114 int sim_nx, sim_ny, sim_nz; 01115 int *zcount , *hist,*gist , z1,z2,z3 ; 01116 MRI_IMAGE *imout_orig = NULL; 01117 01118 ENTRY("mri_brainormalize") ; 01119 01120 if( im == NULL || xxor < 0 || xxor > LAST_ORIENT_TYPE || 01121 yyor < 0 || yyor > LAST_ORIENT_TYPE || 01122 zzor < 0 || zzor > LAST_ORIENT_TYPE ) RETURN(NULL) ; 01123 01124 if( im->nx < 16 || im->ny < 16 || im->nz < 16 ) RETURN(NULL) ; 01125 01126 val = mri_maxabs(im) ; if( val <= 0.0 ) RETURN(NULL) ; 01127 01128 /* make a short copy */ 01129 01130 if( verb ) fprintf(stderr,"++mri_brainormalize: copying input\n") ; 01131 01132 if( im->kind == MRI_short || im->kind == MRI_byte || im->kind == MRI_rgb ) 01133 tim = mri_to_short( 1.0 , im ) ; 01134 else 01135 tim = mri_to_short( 32767.0/val , im ) ; 01136 01137 /* flip to RAI orientation */ 01138 01139 ii = jj = kk = 0 ; 01140 switch( xxor ){ 01141 case ORI_R2L_TYPE: ii = 1 ; break ; 01142 case ORI_L2R_TYPE: ii = -1 ; break ; 01143 case ORI_P2A_TYPE: jj = -1 ; break ; 01144 case ORI_A2P_TYPE: jj = 1 ; break ; 01145 case ORI_I2S_TYPE: kk = 1 ; break ; 01146 case ORI_S2I_TYPE: kk = -1 ; break ; 01147 } 01148 switch( yyor ){ 01149 case ORI_R2L_TYPE: ii = 2 ; break ; 01150 case ORI_L2R_TYPE: ii = -2 ; break ; 01151 case ORI_P2A_TYPE: jj = -2 ; break ; 01152 case ORI_A2P_TYPE: jj = 2 ; break ; 01153 case ORI_I2S_TYPE: kk = 2 ; break ; 01154 case ORI_S2I_TYPE: kk = -2 ; break ; 01155 } 01156 switch( zzor ){ 01157 case ORI_R2L_TYPE: ii = 3 ; break ; 01158 case ORI_L2R_TYPE: ii = -3 ; break ; 01159 case ORI_P2A_TYPE: jj = -3 ; break ; 01160 case ORI_A2P_TYPE: jj = 3 ; break ; 01161 case ORI_I2S_TYPE: kk = 3 ; break ; 01162 case ORI_S2I_TYPE: kk = -3 ; break ; 01163 } 01164 01165 if( ii==1 && jj==2 && kk==3 ){ /* no flip needed */ 01166 sim = tim ; 01167 } else { /* flipization */ 01168 if( verb ) 01169 fprintf(stderr,"++mri_brainormalize: flipping to RAI orientation\n") ; 01170 sim = mri_flip3D( ii,jj,kk , tim ) ; 01171 mri_free(tim) ; 01172 if( sim == NULL ) RETURN(NULL) ; /* bad orientation codes? */ 01173 } 01174 01175 sar = MRI_SHORT_PTR(sim) ; 01176 if( sar == NULL ){ mri_free(sim); RETURN(NULL); } /* bad image? */ 01177 01178 /* make a binary mask */ 01179 01180 nx = sim->nx ; ny = sim->ny ; nz = sim->nz ; nxy = nx*ny ; nxyz = nxy*nz ; 01181 dx = fabs(sim->dx) ; if( dx == 0.0 ) dx = 1.0 ; 01182 dy = fabs(sim->dy) ; if( dy == 0.0 ) dy = 1.0 ; 01183 dz = fabs(sim->dz) ; if( dz == 0.0 ) dz = 1.0 ; 01184 01185 /* save some info to create an output image with the same number of slices as original image*/ 01186 sim_dx = sim->dx; sim_dy = sim->dy; sim_dz = sim->dz; 01187 sim_xo = 0.0; sim_yo = 0.0; sim_zo = 0.0; /* origins are added after this function returns.*/ 01188 sim_nx = sim->nx; sim_ny = sim->ny; sim_nz = sim->nz; 01189 01190 if( verb ) fprintf(stderr,"++mri_brainormalize: making mask\n") ; 01191 mask = mri_short2mask( sim ) ; 01192 01193 if( mask == NULL ){ mri_free(sim); RETURN(NULL); } 01194 01195 /* fill in any isolated holes in mask */ 01196 01197 (void) THD_mask_fillin_once( nx,ny,nz , mask , 2 ) ; /* thd_automask.c */ 01198 THD_mask_dilate ( nx,ny,nz , mask , 5 ) ; 01199 THD_mask_dilate ( nx,ny,nz , mask , 5 ) ; 01200 01201 kk = mask_count(nxyz,mask) ; 01202 if( verb ) 01203 fprintf(stderr,"++mri_brainormalize: filled in mask now has %d voxels\n",kk) ; 01204 01205 if( kk <= 999 ){ free((void *)mask); mri_free(sim); RETURN(NULL); } 01206 01207 /* descending from Superior: 01208 count biggest blob in each slice 01209 find Superiormost location that has 3 01210 slices in a row with "a lot" of stuff 01211 zero out all stuff out above that slice */ 01212 01213 zcount = (int *) malloc( sizeof(int) *nz ) ; /* slice counts */ 01214 for( kk=nz-1 ; kk >= 0 ; kk-- ){ 01215 zcount[kk] = mask_count( nxy , mask+kk*nxy ) ; 01216 } 01217 01218 #if 0 01219 if( verb ){ 01220 fprintf(stderr,"++mri_brainormalize: zcount from top slice #%d\n",nz-1) ; 01221 for( kk=nz-1 ; kk >= 0 ; kk-- ){ 01222 fprintf(stderr," %.3f",((double)(zcount[kk]))/((double)nxy) ) ; 01223 if( (nz-kk)%10 == 0 && kk > 0 ) fprintf(stderr,"\n") ; 01224 } 01225 fprintf(stderr,"\n") ; 01226 } 01227 #endif 01228 01229 /* search down for topmost slice that meets the criterion */ 01230 01231 z1 = (int)(0.010*nxy) ; 01232 z2 = (int)(0.015*nxy) ; 01233 z3 = (int)(0.020*nxy) ; 01234 for( kk=nz-1 ; kk > 2 ; kk-- ) 01235 if( zcount[kk] >= z1 && zcount[kk-1] >= z2 && zcount[kk-2] >= z3 ) break ; 01236 01237 free((void *)zcount) ; 01238 if( kk <= 2 ){ free((void *)mask); mri_free(sim); RETURN(NULL); } 01239 01240 /* zero out all above the slice we just found */ 01241 01242 ktop = kk ; 01243 if( ktop < nz-1 ){ 01244 if( verb ) 01245 fprintf(stderr,"++mri_brainormalize: top clip above slice %d\n",ktop) ; 01246 memset( mask+(ktop+1)*nxy , 0 , nxy*(nz-1-ktop)*sizeof(byte) ) ; 01247 } 01248 01249 /* find slice index THD_BN_ZHEIGHT mm below that top slice */ 01250 01251 jj = (int)( ktop-THD_BN_ZHEIGHT/dz ) ; 01252 if( jj >= 0 ){ 01253 if( verb ) 01254 fprintf(stderr,"++mri_brainormalize: bot clip below slice %d\n",jj) ; 01255 memset( mask , 0 , nxy*(jj+1)*sizeof(byte) ) ; 01256 } 01257 01258 kk = mask_count(nxyz,mask) ; 01259 if( kk <= 999 ){ free((void *)mask); mri_free(sim); RETURN(NULL); } 01260 01261 /* apply mask to image (will also remove any negative values) */ 01262 01263 if( verb ) 01264 fprintf(stderr,"++mri_brainormalize: applying mask to image; %d voxels\n",kk) ; 01265 for( ii=0 ; ii < nxyz ; ii++ ) if( !mask[ii] ) sar[ii] = 0 ; 01266 01267 free((void *)mask) ; /* done with this mask */ 01268 01269 /* compute CM of masked image (indexes, not mm) */ 01270 01271 icm = jcm = kcm = sum = 0.0 ; 01272 #ifndef THD_BN_CMTOP 01273 kbot = 0 ; 01274 ktop = nz-1 ; 01275 #else 01276 kbot = (int)rint( ktop-110.0/dz ); if( kbot < 0 ) kbot = 0; 01277 #endif 01278 for( ijk=kbot*nxy,kk=kbot ; kk <= ktop ; kk++ ){ 01279 for( jj=0 ; jj < ny ; jj++ ){ 01280 for( ii=0 ; ii < nx ; ii++,ijk++ ){ 01281 val = (float)sar[ijk] ; 01282 sum += val ; 01283 icm += val * ii ; 01284 jcm += val * jj ; 01285 kcm += val * kk ; 01286 }}} 01287 if( sum == 0.0 ){ mri_free(sim); RETURN(NULL); } /* huh? */ 01288 01289 ai = thd_bn_dxyz/dx ; bi = icm/sum - ai*(THD_BN_XCM-THD_BN_XORG)/thd_bn_dxyz ; 01290 aj = thd_bn_dxyz/dy ; bj = jcm/sum - aj*(THD_BN_YCM-THD_BN_YORG)/thd_bn_dxyz ; 01291 ak = thd_bn_dxyz/dz ; bk = kcm/sum - ak*(THD_BN_ZCM-THD_BN_ZORG)/thd_bn_dxyz ; 01292 01293 if( verb ) fprintf(stderr,"++mri_brainormalize: warping to standard grid\n a = [%f %f %f], b = [%f %f %f]\n", ai, aj, ak, bi, bj, bk) ; 01294 01295 mri_warp3D_method( MRI_CUBIC ) ; 01296 tim = mri_warp3D( sim , thd_bn_nx,thd_bn_ny,thd_bn_nz , ijkwarp ) ; 01297 mri_free(sim) ; 01298 01299 tim->dx = tim->dy = tim->dz = thd_bn_dxyz ; 01300 tim->xo = THD_BN_XORG ; 01301 tim->yo = THD_BN_YORG ; 01302 tim->zo = THD_BN_ZORG ; 01303 01304 nx = tim->nx ; ny = tim->ny ; nz = tim->nz ; nxy = nx*ny ; nxyz = nxy*nz ; 01305 sar = MRI_SHORT_PTR(tim) ; 01306 01307 /*-- rescale to partially uniformize --*/ 01308 01309 { clipvec bvec ; float bval , sv ; 01310 bvec = get_octant_clips( tim , 0.40f ) ; 01311 if( bvec.clip_000 > 0.0f ){ 01312 for( ijk=kk=0 ; kk < nz ; kk++ ){ 01313 for( jj=0 ; jj < ny ; jj++ ){ 01314 for( ii=0 ; ii < nx ; ii++,ijk++ ){ 01315 bval = pointclip( ii,jj,kk , &bvec ) ; /* cliplevel here */ 01316 sv = 1000.0f * sar[ijk] / bval ; 01317 sar[ijk] = SHORTIZE(sv) ; 01318 } 01319 }}} 01320 } 01321 01322 /*-- build another mask now --*/ 01323 if( !AFNI_noenv("REMASK") ){ 01324 int sbot,stop , nwid , cbot,ctop , ibot,itop ; 01325 float dsum , ws , *wt ; 01326 int pval[128] , wval[128] , npk , tval , nmask,nhalf ; 01327 float pdif ; 01328 short mbot,mtop ; 01329 01330 /* build histogram */ 01331 01332 hist = (int *) calloc(sizeof(int),32768) ; 01333 gist = (int *) calloc(sizeof(int),32768) ; 01334 01335 memset( hist , 0 , sizeof(int)*32768 ) ; 01336 for( ii=0 ; ii < nxyz ; ii++ ) hist[sar[ii]]++ ; 01337 for( sbot=1 ; sbot < 32768 && hist[sbot]==0 ; sbot++ ) ; /* nada */ 01338 if( sbot == 32768 ) goto Remask_Done ; 01339 for( stop=32768-1 ; stop > sbot && hist[stop]==0 ; stop-- ) ; /* nada */ 01340 if( stop == sbot ) goto Remask_Done ; 01341 01342 /* find median */ 01343 01344 nmask = 0 ; 01345 for( ii=sbot ; ii <= stop ; ii++ ) nmask += hist[ii] ; 01346 nhalf = nmask / 2 ; nmask = 0 ; 01347 for( ii=sbot ; ii <= stop && nmask < nhalf ; ii++ ) nmask += hist[ii] ; 01348 cbot = 0.40 * ii ; 01349 ctop = 1.60 * ii ; 01350 01351 #if 0 01352 /* smooth histogram */ 01353 01354 nwid = rint(0.10*cbot) ; 01355 01356 if( nwid <= 0 ){ 01357 memcpy( gist , hist , sizeof(int)*32768 ) ; 01358 } else { 01359 ws = 0.0f ; 01360 wt = (float *)malloc(sizeof(float)*(2*nwid+1)) ; 01361 for( ii=0 ; ii <= 2*nwid ; ii++ ){ 01362 wt[ii] = nwid-abs(nwid-ii) + 0.5f ; 01363 ws += wt[ii] ; 01364 } 01365 for( ii=0 ; ii <= 2*nwid ; ii++ ) wt[ii] /= ws ; 01366 for( jj=cbot ; jj <= ctop ; jj++ ){ 01367 ibot = jj-nwid ; if( ibot < sbot ) ibot = sbot ; 01368 itop = jj+nwid ; if( itop > stop ) itop = stop ; 01369 ws = 0.0 ; 01370 for( ii=ibot ; ii <= itop ; ii++ ) 01371 ws += wt[nwid-jj+ii] * hist[ii] ; 01372 gist[jj] = rint(ws) ; 01373 } 01374 free(wt) ; 01375 } 01376 01377 /* scan for peaks */ 01378 01379 npk = 0 ; 01380 for( ii=cbot+2 ; ii <= ctop-2 ; ii++ ){ 01381 if( gist[ii] > gist[ii-1] && 01382 gist[ii] > gist[ii-2] && 01383 gist[ii] > gist[ii+1] && 01384 gist[ii] > gist[ii+2] ){ 01385 pval[npk]=ii; wval[npk++] = gist[ii]; 01386 } 01387 01388 else if( gist[ii] == gist[ii+1] && /* very special case */ 01389 gist[ii] > gist[ii-1] && 01390 gist[ii] > gist[ii-2] && 01391 gist[ii] > gist[ii+2] ){ 01392 pval[npk]=ii+0.5; wval[npk++] = gist[ii]; 01393 } 01394 01395 else if( gist[ii] == gist[ii+1] && /* super special case */ 01396 gist[ii] == gist[ii-1] && 01397 gist[ii] > gist[ii-2] && 01398 gist[ii] > gist[ii+2] ){ 01399 pval[npk]=ii; wval[npk++] = gist[ii]; 01400 } 01401 } 01402 01403 if( npk > 2 ){ /* find largest two peaks and keep only them */ 01404 float pval_top, pval_2nd, wval_top, wval_2nd , www; int iii,itop ; 01405 www = wval[0] ; iii = 0 ; 01406 for( ii=1 ; ii < npk ; ii++ ){ 01407 if( wval[ii] > www ){ www = wval[ii] ; iii = ii ; } 01408 } 01409 pval_top = pval[iii] ; wval_top = www ; itop = iii ; 01410 www = -1 ; iii = -1 ; 01411 for( ii=0 ; ii < npk ; ii++ ){ 01412 if( ii != itop && wval[ii] > www ){ www = wval[ii] ; iii = ii ; } 01413 } 01414 pval_2nd = pval[iii] ; wval_2nd = www ; 01415 01416 /* make sure peaks are increasing in pval */ 01417 01418 if( pval_top < pval_2nd ){ 01419 pval[0] = pval_top ; wval[0] = wval_top ; 01420 pval[1] = pval_2nd ; wval[1] = wval_2nd ; 01421 } else { 01422 pval[0] = pval_2nd ; wval[0] = wval_2nd ; 01423 pval[1] = pval_top ; wval[1] = wval_top ; 01424 } 01425 npk = 2 ; 01426 } 01427 01428 if( npk == 2 ){ 01429 jj = gist[pval[0]] ; tval = pval[0] ; 01430 for( ii=pval[0]+1 ; ii < pval[1] ; ii++ ){ 01431 if( gist[ii] < jj ){ tval = ii ; jj = gist[ii] ; } 01432 } 01433 01434 pdif = 1.5f * (tval-pval[0]) ; 01435 if( pdif < pval[1]-pval[0] ) pdif = pval[1]-pval[0] ; 01436 mbot = (short)(pval[0]-pdif) ; 01437 if( mbot < cbot ) mbot = cbot ; 01438 01439 pdif = 1.5f * (pval[1]-tval) ; 01440 if( pdif < pval[1]-pval[0] ) pdif = pval[1]-pval[0] ; 01441 mtop = (short)(pval[1]+pdif) ; 01442 if( mtop > ctop ) mtop = ctop ; 01443 } else { 01444 mbot = cbot ; mtop = ctop ; 01445 } 01446 mtop = stop+1 ; /* effectively, no threshold here */ 01447 #endif 01448 01449 mbot = cbot ; mtop = 32767 ; 01450 01451 if( verb ) 01452 fprintf(stderr,"++mri_brainormalize: masking standard image %d..%d\n",mbot,mtop) ; 01453 01454 mask = (byte *) malloc( sizeof(byte)*nxyz ) ; 01455 for( ii=0 ; ii < nxyz ; ii++ ) 01456 mask[ii] = (sar[ii] > mbot) && (sar[ii] < mtop) ; 01457 01458 THD_mask_erode( nx,ny,nz, mask ) ; 01459 THD_mask_clust( nx,ny,nz, mask ) ; 01460 for( ii=0 ; ii < nxyz ; ii++ ) mask[ii] = !mask[ii] ; 01461 THD_mask_clust( nx,ny,nz, mask ) ; 01462 for( ii=0 ; ii < nxyz ; ii++ ) mask[ii] = !mask[ii] ; 01463 01464 for( ii=0 ; ii < nxyz ; ii++ ) if( !mask[ii] ) sar[ii] = 0 ; 01465 free((void *)mask) ; 01466 01467 Remask_Done: 01468 free((void *)hist) ; free((void *)gist) ; 01469 01470 } 01471 #if 0 01472 else 01473 #endif 01474 { 01475 /*-- clip top 1% of values that have survived --*/ 01476 01477 hist = (int *) calloc(sizeof(int),32768) ; 01478 for( ii=0 ; ii < nxyz ; ii++ ) hist[sar[ii]]++ ; 01479 for( ii=kk=0 ; ii < 32767 ; ii++ ) kk += hist[ii] ; 01480 kk = (int)(0.01*kk) ; ktop = 0 ; 01481 for( jj=0,ii=32767 ; ii > 0 && jj < kk ; ii-- ){ 01482 jj += hist[ii] ; if( hist[ii] > 0 && ktop == 0 ) ktop = ii ; 01483 } 01484 jj = ii ; 01485 if( verb ) fprintf(stderr," + 99%% clipping at %d (from %d)\n",jj,ktop) ; 01486 for( ii=0 ; ii < nxyz ; ii++ ) if( sar[ii] > jj ) sar[ii] = jj ; 01487 01488 free((void *)hist) ; 01489 } 01490 01491 /* distize? */ 01492 01493 if( AFNI_yesenv("DISTIZE") ){ 01494 byte *ccc = (byte *)calloc(sizeof(byte),nxyz); 01495 short *ddd ; 01496 int kbot=(int)rint(0.45*nz) , ktop=(int)rint(0.65*nz) , 01497 jbot=(int)rint(0.30*ny) , jtop=(int)rint(0.70*ny) , 01498 ibot=(int)rint(0.30*nx) , itop=(int)rint(0.70*nx) ; 01499 01500 mask = (byte *)malloc( sizeof(byte)*nxyz ) ; 01501 for( ii=0 ; ii < nxyz ; ii++ ) mask[ii] = (sar[ii] > 0) ; 01502 for( kk=kbot ; kk <= ktop ; kk++ ){ 01503 for( jj=jbot ; jj <= jtop ; jj++ ){ 01504 for( ii=ibot ; ii <= itop ; ii++ ){ 01505 ijk = ii + jj*nx + kk*nxy ; 01506 ccc[ijk] = mask[ijk] ; 01507 }}} 01508 if( verb ) fprintf(stderr," + distizing\n") ; 01509 ddd = THD_mask_distize( nx,ny,nz , mask , ccc ) ; 01510 if( ddd != NULL ){ 01511 int id,jd,kd , ijk , dijk ; float ff ; 01512 for( ijk=0 ; ijk < nxyz ; ijk++ ){ 01513 if( ddd[ijk] > 0 ){ 01514 ii = ijk % nx ; jj = (ijk%nxy)/nx ; kk = ijk / nxy ; 01515 if( ii < ibot ) id = ibot-ii ; 01516 else if( ii > itop ) id = ii-itop ; else id = 0 ; 01517 if( jj < jbot ) jd = jbot-jj ; 01518 else if( jj > jtop ) jd = jj-jtop ; else jd = 0 ; 01519 if( kk < kbot ) kd = kbot-kk ; 01520 else if( kk > ktop ) kd = kk-ktop ; else kd = 0 ; 01521 dijk = id+jd+kd+1 ; 01522 ff = (100.0f * ddd[ijk]) / (float)dijk - 98.9f ; 01523 if( ff > 255.0f ) ff = 255.0f ; 01524 sar[ijk] = (short)ff ; 01525 } else { 01526 sar[ijk] = 0 ; 01527 } 01528 } 01529 free((void *)ddd) ; 01530 } 01531 free((void *)mask); free((void *)ccc); 01532 } 01533 01534 /* create a spat norm of the original volume ZSS */ 01535 if (imout_origp) { 01536 mri_warp3D_method( MRI_NN ) ; 01537 if (1 && verb) fprintf(stderr,"thd_brainormalize (ZSS):\n n: %d %d %d\n d: %f %f %f\n o: %f %f %f\n ", sim_nx, sim_ny, sim_nz, sim_dx, sim_dy, sim_dz, sim_xo, sim_yo, sim_zo); 01538 imout_orig = mri_warp3D( tim, sim_nx, sim_ny, sim_nz, ijk_invwarp ); 01539 imout_orig->dx = sim_dx; imout_orig->dy = sim_dy; imout_orig->dz = sim_dz; 01540 imout_orig->xo = sim_xo; imout_orig->yo = sim_yo; imout_orig->zo = sim_zo; 01541 *imout_origp = imout_orig; 01542 } 01543 01544 if (imout_edge) { /* create an edge version NOT EXISTING YET */ 01545 *imout_edge = NULL; 01546 } 01547 01548 /*-- convert output to bytes --*/ 01549 01550 bim = mri_new_conforming( tim , MRI_byte ) ; 01551 MRI_COPY_AUX(bim,tim) ; 01552 bar = MRI_BYTE_PTR(bim) ; 01553 01554 jj = 0 ; 01555 for( ii=0 ; ii < nxyz ; ii++ ) if( sar[ii] > jj ) jj = sar[ii] ; 01556 01557 if( jj > 255 ){ 01558 float fac = 255.0 / jj ; 01559 if( verb ) fprintf(stderr," + scaling by fac=%g\n",fac) ; 01560 for( ii=0 ; ii < nxyz ; ii++ ) bar[ii] = (byte)(fac*sar[ii]+0.49) ; 01561 } else { 01562 for( ii=0 ; ii < nxyz ; ii++ ) bar[ii] = (byte)sar[ii] ; 01563 } 01564 mri_free(tim) ; 01565 01566 /*-- done!!! --*/ 01567 01568 RETURN(bim) ; 01569 } |
|
Definition at line 13 of file thd_brainormalize.c. References MIN, thd_bn_dxyz, thd_bn_nx, thd_bn_ny, and thd_bn_nz. Referenced by main().
00014 { 00015 /* Set the reinterpolation resolution to the smallest delta */ 00016 thd_bn_dxyz = MIN(fabs(dx), fabs(dy)); thd_bn_dxyz = MIN(thd_bn_dxyz, fabs(dz)); 00017 thd_bn_nx = (int)( (float)thd_bn_nx / thd_bn_dxyz ); 00018 thd_bn_ny = (int)( (float)thd_bn_ny / thd_bn_dxyz ); 00019 thd_bn_nz = (int)( (float)thd_bn_nz / thd_bn_dxyz ); 00020 return; 00021 } |
|
Definition at line 8 of file thd_brainormalize.c. Referenced by main().
|
|
Make a mask from the 3D image, based on gradually changing cliplevels. This feature is to allow for gradual shading in the MRI volume. Definition at line 467 of file thd_brainormalize.c. References clipvec::clip_000, clipvec::clip_001, clipvec::clip_010, clipvec::clip_011, clipvec::clip_100, clipvec::clip_101, clipvec::clip_110, clipvec::clip_111, clustedit3D(), ENTRY, get_octant_clips(), MRI_IMAGE::kind, malloc, mask_count(), MRI_SHORT_PTR, MRI_IMAGE::nx, MRI_IMAGE::ny, MRI_IMAGE::nz, nz, pointclip(), and RETURN. Referenced by mri_brainormalize().
00468 { 00469 int ii,jj,kk,ijk , nx,ny,nz,nxy,nxyz ; 00470 clipvec bvec , tvec ; 00471 short *sar ; 00472 byte *mask ; 00473 float bval , tval ; 00474 00475 ENTRY("mri_short2mask") ; 00476 if( im == NULL || im->kind != MRI_short ) RETURN(NULL) ; 00477 sar = MRI_SHORT_PTR(im) ; if( sar == NULL ) RETURN(NULL) ; 00478 00479 nx = im->nx ; ny = im->ny ; nz = im->nz ; nxy = nx*ny ; nxyz = nxy*nz ; 00480 00481 bvec = get_octant_clips( im , 0.40f ) ; 00482 if( bvec.clip_000 < 0.0 ) RETURN(NULL) ; 00483 00484 tvec = bvec ; 00485 tvec.clip_000 *= 9.91 ; 00486 tvec.clip_100 *= 9.91 ; 00487 tvec.clip_010 *= 9.91 ; 00488 tvec.clip_110 *= 9.91 ; 00489 tvec.clip_001 *= 9.91 ; 00490 tvec.clip_101 *= 9.91 ; 00491 tvec.clip_011 *= 9.91 ; 00492 tvec.clip_111 *= 9.91 ; 00493 00494 /* create mask, clipping at a level that varies spatially */ 00495 00496 #if 0 00497 if( verb ) fprintf(stderr," + mri_short2mask: clipping\n") ; 00498 #endif 00499 00500 mask = (byte *) malloc( sizeof(byte)*nxyz ) ; 00501 for( ijk=kk=0 ; kk < nz ; kk++ ){ 00502 for( jj=0 ; jj < ny ; jj++ ){ 00503 for( ii=0 ; ii < nx ; ii++,ijk++ ){ 00504 bval = pointclip( ii,jj,kk , &bvec ) ; /* cliplevel here */ 00505 #if 0 00506 tval = pointclip( ii,jj,kk , &tvec ) ; /* cliplevel here */ 00507 mask[ijk] = (sar[ijk] >= bval && sar[ijk] <= tval) ; /* binarize */ 00508 #else 00509 mask[ijk] = (sar[ijk] >= bval) ; /* binarize */ 00510 #endif 00511 }}} 00512 00513 /* remove small clusters */ 00514 00515 clustedit3D( nx,ny,nz , mask , (int)rint(0.02*nxyz) ) ; 00516 00517 if( verb ) fprintf(stderr," + mri_short2mask: %d voxels survive\n", 00518 mask_count(nxyz,mask) ) ; 00519 00520 RETURN(mask) ; 00521 } |
|
Definition at line 676 of file thd_brainormalize.c. References ADDTO_BASIN, shortvox::basin, BDEP, calloc, DBALL, ENTRY, free, i, shortvox::i, IJK, INIT_BASIN, shortvox::j, shortvox::k, KILL_BASIN, MRI_IMAGE::kind, malloc, MERGE_BASIN, MRI_COPY_AUX, mri_new_conforming, MRI_SHORT_PTR, MRI_IMAGE::nx, MRI_IMAGE::ny, MRI_IMAGE::nz, nz, qsort_intint(), RETURN, sort_shortvox(), shortvox::val, and verb. Referenced by main().
00677 { 00678 MRI_IMAGE *tim ; 00679 int ii,jj,kk , pp,qq , nx,ny,nz,nxy,nxyz , nvox ; 00680 int ip,jp,kp , im,jm,km ; 00681 short *sar , *tar ; 00682 shortvox *svox ; 00683 int *isvox , *bcount,*bname ; 00684 int nb,vb,mb,m,mu,mq,mz , bp[6] , hpf ; 00685 00686 basin **baslist ; 00687 int nball , nbtop ; 00688 00689 ENTRY("watershedize") ; 00690 00691 if( sim == NULL || sim->kind != MRI_short ) RETURN(NULL) ; 00692 sar = MRI_SHORT_PTR(sim) ; if( sar == NULL ) RETURN(NULL) ; 00693 00694 nx = sim->nx; ny = sim->ny; nz = sim->nz; nxy = nx*ny; nxyz = nxy*nz; 00695 00696 /* count number of voxels > 0 */ 00697 00698 for( nvox=0,pp=0 ; pp < nxyz ; pp++ ) if( sar[pp] > 0 ) nvox++ ; 00699 if( nvox <= 999 ) RETURN(NULL) ; 00700 00701 if( verb ) fprintf(stderr," + mri_watershedize: %d voxels input\n",nvox) ; 00702 00703 /* create voxel lists */ 00704 00705 svox = (shortvox *) malloc( sizeof(shortvox)* nvox ) ; 00706 isvox = (int *) malloc( sizeof(int) * nxyz ) ; 00707 for( qq=pp=0 ; pp < nxyz ; pp++ ){ 00708 if( sar[pp] > 0 ){ /* save this one: */ 00709 ii = pp % nx ; /* spatial indexes */ 00710 jj = (pp%nxy) / nx ; 00711 kk = pp / nxy ; 00712 svox[qq].i = ii ; 00713 svox[qq].j = jj ; 00714 svox[qq].k = kk ; 00715 svox[qq].val = sar[pp] ; /* value */ 00716 svox[qq].basin = -1 ; /* which basin */ 00717 qq++ ; 00718 isvox[pp] = qq ; /* where in list */ 00719 } else { 00720 isvox[pp] = -1 ; /* voxel not in list */ 00721 } 00722 } 00723 00724 /* sort voxel list into descending order */ 00725 00726 if( verb ) fprintf(stderr," + mri_watershedize: sorting voxels\n") ; 00727 00728 sort_shortvox( nvox , svox , 1 , 0.00 , 0.02 ) ; 00729 00730 /* create basin for first (deepest) voxel */ 00731 00732 nball = DBALL ; 00733 nbtop = 0 ; 00734 baslist = (basin **) calloc(sizeof(basin *),nball) ; 00735 00736 INIT_BASIN(0) ; 00737 00738 hpf = (int)rint(prefac*svox[0].val) ; /* preflood */ 00739 00740 /* scan voxels as they get shallower, and basinate them */ 00741 00742 if( verb ){ 00743 fprintf(stderr," + mri_watershedize: basinating voxels\n") ; 00744 fprintf(stderr," data range: %d..%d preflood_height=%d\n", 00745 svox[nvox-1].val , svox[0].val , hpf ) ; 00746 } 00747 00748 for( pp=1 ; pp < nvox ; pp++ ){ 00749 00750 ii = svox[pp].i; jj = svox[pp].j; kk = svox[pp].k; /* where */ 00751 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ; /* nbhrs */ 00752 im = ii-1 ; jm = jj-1 ; km = kk-1 ; 00753 00754 if( verb && pp%100000 == 0 ) fprintf(stderr, (pp%1000000)?".":"!") ; 00755 00756 /* macro checks if (a,b,c) voxel is in the list; 00757 if so and it is already in a basin, then 00758 make a list of basins encountered: 00759 nb = number of unique basins encountered (0..6) 00760 mb = index of deepest basin encountered (0..nb-1) 00761 vb = value (depth) of deepest basin encountered 00762 bp[m] = index of m-th basin encountered (m=0..nb-1) */ 00763 00764 #undef BASECHECK 00765 #define BASECHECK(a,b,c) \ 00766 { qq = isvox[IJK(a,b,c)] ; \ 00767 if( qq >= 0 && svox[qq].basin >= 0 ){ \ 00768 qq = svox[qq].basin ; \ 00769 for( m=0 ; m < nb && bp[m] != qq ; m++ ) ; \ 00770 if( m == nb ){ \ 00771 bp[nb] = qq ; \ 00772 if( BDEP(qq) > vb ){ mb = nb; vb = BDEP(qq); } \ 00773 nb++ ; \ 00774 } \ 00775 } \ 00776 } 00777 00778 nb = 0 ; vb = -1 ; mb = -1 ; /* initialize counters */ 00779 if( ip < nx ) BASECHECK(ip,jj,kk) ; /* check each neighbor */ 00780 if( im >= 0 ) BASECHECK(im,jj,kk) ; /* for basin-ositiness */ 00781 if( jp < ny ) BASECHECK(ii,jp,kk) ; 00782 if( jm >= 0 ) BASECHECK(ii,jm,kk) ; 00783 if( kp < nz ) BASECHECK(ii,jj,kp) ; 00784 if( km >= 0 ) BASECHECK(ii,jj,km) ; 00785 00786 if( nb == 0 ){ /*** this voxel is isolated ==> create new basin ****/ 00787 00788 INIT_BASIN(pp) ; 00789 00790 } else { /*** this voxel has deeper neighbors ***/ 00791 00792 mq = bp[mb] ; /* assign voxel to best basin */ 00793 ADDTO_BASIN( mq , pp ) ; 00794 00795 /* if have more than one neighbor, other */ 00796 if( nb > 1 ){ /* basins could be merged with the best */ 00797 mz = svox[pp].val ; /* depth of this voxel */ 00798 for( m=0 ; m < nb ; m++ ){ 00799 if( m == mb ) continue ; /* can't merge with itself */ 00800 mu = bp[m] ; 00801 if( BDEP(mu)-mz <= hpf ){ /* basin not TOO much deeper */ 00802 MERGE_BASIN(mq,mu) ; 00803 } 00804 } 00805 } 00806 } 00807 } /* end of loop over voxels */ 00808 00809 /* at this point, all voxels in svox are assigned to a basin */ 00810 00811 free((void *)isvox) ; 00812 00813 /* count number of basines left */ 00814 00815 for( mu=m=0 ; m < nbtop ; m++ ) 00816 if( baslist[m] != NULL ) mu++ ; 00817 00818 if( verb ) fprintf(stderr,"\n++ %d active basins left, out of %d\n",mu,nbtop) ; 00819 00820 bcount = (int *) calloc(sizeof(int),mu) ; /* number in each basin */ 00821 bname = (int *) calloc(sizeof(int),mu) ; 00822 isvox = (int *) calloc(sizeof(int),nbtop) ; /* new index */ 00823 00824 for( m=ii=0 ; m < nbtop ; m++ ) 00825 if( baslist[m] != NULL ){ isvox[m] = ii; bname[ii] = ii; ii++; KILL_BASIN(m); } 00826 free((void *)baslist) ; 00827 00828 for( pp=0 ; pp < nvox ; pp++ ){ 00829 m = svox[pp].basin ; /* old basin name for this voxel */ 00830 ii = isvox[m] ; /* new basin name for this voxel */ 00831 svox[pp].basin = ii ; /* reassign name in this voxel */ 00832 bcount[ii]++ ; /* count number in this basin */ 00833 } 00834 00835 tim = mri_new_conforming( sim , MRI_short ) ; /* output image */ 00836 MRI_COPY_AUX(tim,sim) ; 00837 tar = MRI_SHORT_PTR(tim) ; 00838 00839 for( ii=0 ; ii < mu ; ii++ ) bcount[ii] = -bcount[ii] ; 00840 qsort_intint( mu , bcount , bname ) ; /* sort into decreasing order */ 00841 for( ii=0 ; ii < mu ; ii++ ) bcount[ii] = -bcount[ii] ; 00842 00843 if( verb ) 00844 fprintf(stderr," + top 9 basin counts: %d %d %d %d %d %d %d %d %d\n", 00845 bcount[0] , bcount[1] , bcount[2] , bcount[3] , 00846 bcount[4] , bcount[5] , bcount[6] , bcount[7] , bcount[8] ) ; 00847 00848 for( ii=0 ; ii < mu ; ii++ ) isvox[ii] = ii ; 00849 qsort_intint( mu , bname , isvox ) ; 00850 00851 for( pp=0 ; pp < nvox ; pp++ ){ 00852 m = svox[pp].basin ; jj = isvox[m]+1 ; if( jj > 32767 ) jj = 32767 ; 00853 tar[IJK(svox[pp].i,svox[pp].j,svox[pp].k)] = jj ; 00854 } 00855 00856 free((void *)isvox) ; free((void *)svox ); 00857 free((void *)bcount); free((void *)bname); 00858 00859 return tim ; 00860 } |
|
Find the cliplevel in just part of the image, with i=ibot..itop (inclusive), et cetera. Definition at line 258 of file thd_brainormalize.c. References calloc, ENTRY, free, IJK, MRI_IMAGE::kind, MRI_SHORT_PTR, MRI_IMAGE::nx, MRI_IMAGE::ny, MRI_IMAGE::nz, nz, and RETURN. Referenced by get_octant_clips().
00261 { 00262 int ncut,ib , qq,nold ; 00263 int ii,jj,kk , nx,ny,nz,nxy ; 00264 int *hist ; 00265 short *sar ; 00266 byte *bar ; 00267 int nhist , npos , nhalf , val ; 00268 double dsum ; 00269 00270 ENTRY("partial_cliplevel") ; 00271 if( im == NULL || im->kind != MRI_short ) RETURN(1.0) ; 00272 00273 if( mfrac <= 0.0 || mfrac >= 0.99 ) mfrac = 0.50 ; 00274 00275 nx = im->nx ; ny = im->ny ; nz = im->nz ; nxy = nx*ny ; 00276 00277 if( ibot < 0 ) ibot = 0 ; 00278 if( jbot < 0 ) jbot = 0 ; 00279 if( kbot < 0 ) kbot = 0 ; 00280 if( itop >= nx ) itop = nx-1 ; 00281 if( jtop >= ny ) jtop = ny-1 ; 00282 if( ktop >= nz ) ktop = nz-1 ; 00283 00284 if( itop < ibot || jtop < jbot || ktop < kbot ) RETURN(1.0) ; 00285 00286 /*-- allocate histogram --*/ 00287 00288 nhist = 32767 ; 00289 hist = (int *) calloc(sizeof(int),nhist+1) ; 00290 00291 /*-- make histogram of positive entries --*/ 00292 00293 dsum = 0.0 ; /* will be sum of squares */ 00294 npos = 0 ; /* will be number of positive values */ 00295 sar = MRI_SHORT_PTR(im) ; 00296 for( kk=kbot ; kk <= ktop ; kk++ ){ 00297 for( jj=jbot ; jj <= jtop ; jj++ ){ 00298 for( ii=ibot ; ii <= itop ; ii++ ){ 00299 val = sar[IJK(ii,jj,kk)] ; 00300 if( val > 0 && val <= nhist ){ 00301 hist[val]++ ; 00302 dsum += (double)(val)*(double)(val); npos++; 00303 } 00304 }}} 00305 00306 if( npos <= 999 ){ free((void *)hist); RETURN(1.0); } 00307 00308 /*-- initialize cut position to include upper 65% of positive voxels --*/ 00309 00310 qq = 0.65 * npos ; ib = rint(0.5*sqrt(dsum/npos)) ; 00311 for( kk=0,ii=nhist-1 ; ii >= ib && kk < qq ; ii-- ) kk += hist[ii] ; 00312 00313 /*-- algorithm --*/ 00314 00315 ncut = ii ; /* initial cut position */ 00316 qq = 0 ; /* iteration count */ 00317 do{ 00318 for( npos=0,ii=ncut ; ii < nhist ; ii++ ) npos += hist[ii]; /* number >= cut */ 00319 nhalf = npos/2 ; 00320 for( kk=0,ii=ncut ; ii < nhist && kk < nhalf ; ii++ ) /* find median of */ 00321 kk += hist[ii] ; /* valuess >= cut */ 00322 nold = ncut ; /* last cut */ 00323 ncut = mfrac * ii ; /* new cut */ 00324 qq++ ; 00325 } while( qq < 20 && ncut != nold ) ; /* iterate until done, or at most 20 times */ 00326 00327 free((void *)hist) ; 00328 RETURN( (float)(ncut) ) ; 00329 } |
|
Return the cliplevel at any point, tri-linearly interpolated between octant centers. Definition at line 444 of file thd_brainormalize.c. References clipvec::clip_000, clipvec::clip_001, clipvec::clip_010, clipvec::clip_011, clipvec::clip_100, clipvec::clip_101, clipvec::clip_110, clipvec::clip_111, clipvec::dxi, clipvec::dyi, clipvec::dzi, clipvec::x0, x0, clipvec::y0, y0, y1, clipvec::z0, z0, and z1. Referenced by mri_brainormalize(), and mri_short2mask().
00445 { 00446 float x1,y1,z1 , x0,y0,z0 , val ; 00447 00448 /* get relative position in box defined by octant centers */ 00449 00450 x1 = (ii-cv->x0)*cv->dxi; if(x1 < 0.0) x1=0.0; else if(x1 > 1.0) x1=1.0; 00451 y1 = (jj-cv->y0)*cv->dyi; if(y1 < 0.0) y1=0.0; else if(y1 > 1.0) y1=1.0; 00452 z1 = (kk-cv->z0)*cv->dzi; if(z1 < 0.0) z1=0.0; else if(z1 > 1.0) z1=1.0; 00453 00454 x0 = 1.0-x1 ; y0 = 1.0-y1 ; z0 = 1.0-z1 ; 00455 00456 val = cv->clip_000 * x0*y0*z0 + cv->clip_100 * x1*y0*z0 00457 + cv->clip_010 * x0*y1*z0 + cv->clip_110 * x1*y1*z0 00458 + cv->clip_001 * x0*y0*z1 + cv->clip_101 * x1*y0*z1 00459 + cv->clip_011 * x0*y1*z1 + cv->clip_111 * x1*y1*z1 ; 00460 return val ; 00461 } |
|
Sort array of shortvox into increasing order (decreasing if dec != 0). Definition at line 530 of file thd_brainormalize.c. References calloc, free, and shortvox::val. Referenced by mri_watershedize().
00531 { 00532 int ii,jj , sbot,stop,nsv , sval , pbot,ptop ; 00533 int *hsv , *csv ; 00534 shortvox *tar ; 00535 00536 if( n < 2 || ar == NULL ) return ; 00537 00538 /* decreasing order desired? flip values */ 00539 00540 if( dec ){ 00541 float tmp ; 00542 for( ii=0 ; ii < n ; ii++ ) ar[ii].val = -ar[ii].val ; 00543 tmp = botperc ; botperc = topperc ; topperc = tmp ; 00544 } 00545 00546 /* find range of values */ 00547 00548 sbot = stop = ar[0].val ; 00549 for( ii=1 ; ii < n ; ii++ ){ 00550 sval = ar[ii].val ; 00551 if( sval < sbot ) sbot = sval ; 00552 else if( sval > stop ) stop = sval ; 00553 } 00554 nsv = stop-sbot+1 ; /* number of distinct values */ 00555 if( nsv <= 1 ) return ; 00556 00557 /* build hsv[i] = how many have value = sbot+i 00558 csv[i] = how many have value < sbot+i, i=0..nsv-1 */ 00559 00560 hsv = (int *)calloc(sizeof(int),nsv) ; 00561 csv = (int *)calloc(sizeof(int),nsv+1) ; 00562 for( ii=0 ; ii < n ; ii++ ) hsv[ar[ii].val-sbot]++ ; 00563 for( ii=1 ; ii <= nsv ; ii++ ) csv[ii] = csv[ii-1]+hsv[ii-1] ; 00564 free((void *)hsv) ; 00565 00566 if( botperc > 0.0 && botperc < 50.0 ){ 00567 jj = (int)rint(0.01*botperc*n) ; 00568 for( ii=0 ; ii < nsv && csv[ii] <= jj ; ii++ ) ; 00569 pbot = ii+sbot ; 00570 if( verb ) fprintf(stderr," + sort_shortvox: sbot=%d pbot=%d\n",sbot,pbot) ; 00571 } else { 00572 pbot = sbot ; 00573 } 00574 00575 if( topperc > 0.0 && topperc < 50.0 ){ 00576 jj = (int)rint(0.01*(100.0-topperc)*n) ; 00577 for( ii=0 ; ii < nsv && csv[ii] <= jj ; ii++ ) ; 00578 ptop = ii+sbot ; 00579 if( verb ) fprintf(stderr," + sort_shortvox: stop=%d ptop=%d\n",stop,ptop) ; 00580 } else { 00581 ptop = stop ; 00582 } 00583 00584 /* copy from ar into temp array tar, 00585 putting each one into its place as given by csv */ 00586 00587 tar = (shortvox *)calloc(sizeof(shortvox),n) ; 00588 for( ii=0 ; ii < n ; ii++ ){ 00589 sval = ar[ii].val - sbot ; /* sval is in 0..nsv-1 now */ 00590 tar[ csv[sval] ] = ar[ii] ; 00591 csv[sval]++ ; 00592 } 00593 00594 if( pbot > sbot ){ 00595 for( ii=0 ; ii < n ; ii++ ) 00596 if( tar[ii].val < pbot ) tar[ii].val = pbot ; 00597 } 00598 if( ptop < stop ){ 00599 for( ii=0 ; ii < n ; ii++ ) 00600 if( tar[ii].val > ptop ) tar[ii].val = ptop ; 00601 } 00602 00603 /* copy back into ar */ 00604 00605 memcpy( ar , tar , sizeof(shortvox)*n ) ; 00606 free((void *)tar) ; free((void *)csv) ; 00607 00608 /* unflip? */ 00609 00610 if( dec ) 00611 for( ii=0 ; ii < n ; ii++ ) ar[ii].val = -ar[ii].val ; 00612 00613 return ; 00614 } |
|
Definition at line 22 of file thd_brainormalize.c. References thd_bn_dxyz.
00023 { 00024 return thd_bn_dxyz; 00025 } |
|
Definition at line 26 of file thd_brainormalize.c. References thd_bn_nx.
00027 { 00028 return thd_bn_nx; 00029 } |
|
Definition at line 30 of file thd_brainormalize.c. References thd_bn_ny.
00031 { 00032 return thd_bn_ny; 00033 } |
|
Definition at line 34 of file thd_brainormalize.c. References thd_bn_nz.
00035 { 00036 return thd_bn_nz; 00037 } |
|
Definition at line 74 of file thd_brainormalize.c. References DALL, DPUT, free, malloc, mmm, and nz. Referenced by mri_brainormalize().
00075 { 00076 short *ddd , dnow ; 00077 int ii,jj,kk , nxy=nx*ny , nxyz=nx*ny*nz , ijk ; 00078 int ip,jp,kp , im,jm,km , icl ; 00079 int nccc,nmmm , nall,nnow ; 00080 short *inow , *jnow , *know ; 00081 float drat ; 00082 00083 if( mmm == NULL || ccc == NULL ) return NULL ; 00084 00085 ddd = (short *)malloc( sizeof(short)*nxyz ) ; 00086 nccc = nmmm = 0 ; 00087 for( ii=0 ; ii < nxyz ; ii++ ){ 00088 if( ccc[ii] ){ ddd[ii] = 1; nccc++; nmmm++; } 00089 else if( mmm[ii] ){ ddd[ii] = -1; nmmm++; } 00090 else { ddd[ii] = 0; } 00091 } 00092 if( nccc == 0 ){ free((void *)ddd); return NULL; } 00093 00094 nall = nccc+DALL ; /* # allocated pts */ 00095 inow = (short *) malloc(sizeof(short)*nall) ; /* coords of pts */ 00096 jnow = (short *) malloc(sizeof(short)*nall) ; 00097 know = (short *) malloc(sizeof(short)*nall) ; 00098 nnow = 0 ; 00099 00100 for( ii=0 ; ii < nxyz ; ii++ ){ 00101 if( ccc[ii] ){ 00102 inow[nnow] = ii % nx ; 00103 jnow[nnow] = (ii%nxy)/nx ; 00104 know[nnow] = ii / nxy ; 00105 mmm[ii] = 0 ; 00106 nnow++ ; 00107 } 00108 } 00109 00110 for( icl=0 ; icl < nnow ; icl++ ){ 00111 ii = inow[icl] ; jj = jnow[icl] ; kk = know[icl] ; ijk = ii+jj*nx+kk*nxy ; 00112 im = ii-1 ; jm = jj-1 ; km = kk-1 ; 00113 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ; dnow = ddd[ijk]+1 ; 00114 00115 if( im >= 0 ) DPUT(im,jj,kk,dnow) ; 00116 if( ip < nx ) DPUT(ip,jj,kk,dnow) ; 00117 if( jm >= 0 ) DPUT(ii,jm,kk,dnow) ; 00118 if( jp < ny ) DPUT(ii,jp,kk,dnow) ; 00119 if( km >= 0 ) DPUT(ii,jj,km,dnow) ; 00120 if( kp < nz ) DPUT(ii,jj,kp,dnow) ; 00121 } 00122 00123 for( ii=0 ; ii < nxyz ; ii++ ) mmm[ii] = (ddd[ii] > 0) ; 00124 00125 free((void *)inow); free((void *)jnow); free((void *)know); 00126 00127 return ddd ; 00128 } |
|
Definition at line 945 of file thd_brainormalize.c. References calloc, free, mask_count(), mmm, nz, THD_mask_clust(), THD_mask_erode(), and top.
00946 { 00947 int nxy=nx*ny,nxyz=nxy*nz , ii,jj,kk , ijk , bot,top ; 00948 int zd=zdepth , zb=zbot , zt , zslab , zz ; 00949 byte *ppp , *zzz ; 00950 00951 if( mmm == NULL ) return ; 00952 00953 if( zd < 1 ) zd = 1 ; 00954 if( zb < zd ) zb = zd ; 00955 zslab = 2*zd+1 ; 00956 00957 for( kk=nz-1 ; kk >= zb ; kk-- ){ 00958 jj = mask_count( nxy , mmm+kk*nxy ) ; 00959 if( jj > 0.005*nxy ) break ; 00960 } 00961 zt = kk-zd ; if( zt < zb ) return ; 00962 00963 ppp = (byte *)calloc(sizeof(byte),nxyz) ; 00964 zzz = (byte *)calloc(sizeof(byte),nxy*zslab) ; 00965 00966 for( zz=zb ; zz <= zt ; zz++ ){ 00967 memcpy( zzz , mmm+(zz-zd)*nxy , nxy*zslab ) ; 00968 THD_mask_erode( nx,ny,zslab, zzz ) ; 00969 THD_mask_clust( nx,ny,zslab, zzz ) ; 00970 memcpy( ppp+zz*nxy , zzz+zd*nxy , nxy ) ; 00971 } 00972 free((void *)zzz) ; 00973 memcpy( mmm+zb*nxy , ppp+zb*nxy , (zt-zb+1)*nxy ) ; 00974 free((void *)ppp) ; 00975 THD_mask_erode( nx,ny,nz, mmm ) ; 00976 THD_mask_clust( nx,ny,nz, mmm ) ; 00977 } |
Variable Documentation
|
Index warping function for mri_warp3D() call. Definition at line 984 of file thd_brainormalize.c. Referenced by ijk_invwarp(), ijkwarp(), and mri_brainormalize(). |
|
Index warping function for mri_warp3D() call. Definition at line 984 of file thd_brainormalize.c. Referenced by ijk_invwarp(), ijkwarp(), and mri_brainormalize(). |
|
Index warping function for mri_warp3D() call. Definition at line 984 of file thd_brainormalize.c. Referenced by ijk_invwarp(), ijkwarp(), and mri_brainormalize(). |
|
Index warping function for mri_warp3D() call. Definition at line 984 of file thd_brainormalize.c. Referenced by ijk_invwarp(), ijkwarp(), and mri_brainormalize(). |
|
Index warping function for mri_warp3D() call. Definition at line 984 of file thd_brainormalize.c. Referenced by ijk_invwarp(), ijkwarp(), and mri_brainormalize(). |
|
Index warping function for mri_warp3D() call. Definition at line 984 of file thd_brainormalize.c. Referenced by ijk_invwarp(), ijkwarp(), and mri_brainormalize(). |
|
Definition at line 9 of file thd_brainormalize.c. Referenced by mri_brainormalize(), mri_brainormalize_initialize(), and THD_BN_dxyz(). |
|
Definition at line 10 of file thd_brainormalize.c. Referenced by mri_brainormalize(), mri_brainormalize_initialize(), and THD_BN_nx(). |
|
Definition at line 11 of file thd_brainormalize.c. Referenced by mri_brainormalize(), mri_brainormalize_initialize(), and THD_BN_ny(). |
|
Definition at line 12 of file thd_brainormalize.c. Referenced by mri_brainormalize(), mri_brainormalize_initialize(), and THD_BN_nz(). |
|
Definition at line 7 of file thd_brainormalize.c. Referenced by mri_brainormalize(), mri_brainormalize_verbose(), and mri_watershedize(). |