1 #include "mrilib.h"
2
3 /*--------------------------------------------------------------------------*/
4
5 #undef IJK
6 #define IJK(p,q,r) ((p)+(q)*nx+(r)*nxy)
7
8 #undef NRMAX
9 #define NRMAX 19
10
mri_mask_localcount(MRI_IMAGE * maskim,int iv,int jv,int kv,int nrad,float * rad,short * ccount,short * vcount)11 int mri_mask_localcount( MRI_IMAGE *maskim , int iv, int jv, int kv ,
12 int nrad, float *rad, short *ccount, short *vcount )
13 {
14 byte *mmm ;
15 MCW_cluster *cl ;
16 int ii,jj,kk , ncl , nx,ny,nz,nxy,nxyz , ijk ;
17 int ip,jp,kp , im,jm,km , icl , ic,jc,kc , rr ;
18 float dx,dy,dz , val , di,dj,dk , rsq , rq[NRMAX];
19
20 ENTRY("mri_mask_localcount") ;
21
22 nx = maskim->nx; ny = maskim->ny; nz = maskim->nz; nxy = nx*ny; nxyz = nxy*nz;
23 mmm = MRI_BYTE_PTR(maskim) ;
24 ijk = IJK(iv,jv,kv) ; if( mmm[ijk] == 0 ) RETURN(0) ;
25
26 ii = iv ; jj = jv ; kk = kv ;
27 dx = maskim->dx; dy = maskim->dy; dz = maskim->dz;
28
29 #undef CCLUST
30 #define CCLUST(p,q,r) \
31 do{ int pqr=IJK(p,q,r) ; \
32 if( mmm[pqr] ){ \
33 di = (ii-(p))*dx; dj = (jj-(q))*dy; dk = (kk-(r))*dz; \
34 val = di*di+dj*dj+dk*dk ; \
35 if( val <= rsq ){ ADDTO_CLUSTER(cl,p,q,r,val); mmm[pqr]=0; } \
36 } \
37 } while(0)
38
39 INIT_CLUSTER(cl) ; ADDTO_CLUSTER(cl,ii,jj,kk,0) ; mmm[ijk] = 0 ;
40 for( rr=0 ; rr < nrad ; rr++ ){
41 rq[rr] = rsq = SQR(rad[rr]) ;
42 for( icl=0 ; icl < cl->num_pt ; icl++ ){
43 ic = cl->i[icl]; jc = cl->j[icl]; kc = cl->k[icl];
44 im = ic-1 ; jm = jc-1 ; km = kc-1 ;
45 ip = ic+1 ; jp = jc+1 ; kp = kc+1 ;
46 if( im >= 0 && im < nx ) CCLUST(im,jc,kc) ;
47 if( ip >= 0 && ip < nx ) CCLUST(ip,jc,kc) ;
48 if( jm >= 0 && jm < ny ) CCLUST(ic,jm,kc) ;
49 if( jp >= 0 && jp < ny ) CCLUST(ic,jp,kc) ;
50 if( km >= 0 && km < nz ) CCLUST(ic,jc,km) ;
51 if( kp >= 0 && kp < nz ) CCLUST(ic,jc,kp) ;
52 }
53 ccount[rr] = (short)cl->num_pt ;
54 }
55 for( icl=0 ; icl < cl->num_pt ; icl++ )
56 mmm[IJK(cl->i[icl],cl->j[icl],cl->k[icl])] = 1 ; /* restore mmm */
57 KILL_CLUSTER(cl) ;
58
59 val = rad[nrad-1] / dx ;
60 im = (int)(ii-val) ; if( im < 0 ) im = 0 ;
61 ip = (int)(ii+val+1.0) ; if( ip >= nx ) ip = nx-1 ;
62 val = rad[nrad-1] / dy ;
63 jm = (int)(jj-val) ; if( jm < 0 ) jm = 0 ;
64 jp = (int)(jj+val+1.0) ; if( jp >= ny ) jp = ny-1 ;
65 val = rad[nrad-1] / dz ;
66 km = (int)(kk-val) ; if( km < 0 ) km = 0 ;
67 kp = (int)(kk+val+1.0) ; if( kp >= nz ) kp = nz-1 ;
68 for( rr=0 ; rr < nrad ; rr++ ) vcount[rr] = 0 ;
69 icl = 0 ;
70 for( kc=km ; kc <= kp ; kc++ ){
71 dk = (kk-kc)*dz; dk = dk*dk ;
72 for( jc=jm ; jc <= jp ; jc++ ){
73 dj = (jj-jc)*dy; dj = dj*dj + dk ;
74 for( ic=im ; ic <= ip ; ic++ ){
75 di = (ii-ic)*dx; val = di*di+dj ;
76 for( rr=0 ; rr < nrad ; rr++ ) if( val <= rq[rr] ) vcount[rr]++ ;
77 }}}
78
79 RETURN(1) ;
80 }
81
82 /*--------------------------------------------------------------------------*/
83
main(int argc,char * argv[])84 int main( int argc , char *argv[] )
85 {
86 THD_3dim_dataset *inset , *outset ;
87 MRI_IMAGE *medim , *mim ;
88 byte *mmm ; float *mar ;
89 short *dar05,*dar10,*dar15,*dar20 ;
90 float *dim10,*dim15,*dim20 ;
91 int nx,ny,nz , nxy,nxyz , ii,jj,kk , nmm , icm,jcm,kcm , ijk ;
92 int ibot,itop , jbot,jtop , kbot,ktop ;
93 float isum,jsum,ksum , clip_val , lg10,lg15,lg20 ;
94 int id,jd,kd , dijk , di,dj,dk ; float ff ;
95 char *prefix = "dimize" ;
96 #define NRAD 6
97 float rad[NRAD] = { 5.0f , 8.0f , 11.0f , 14.0f , 17.0f , 20.0f } ;
98 float *bb[NRAD] ;
99 short ccc[NRAD] , vvv[NRAD] ; int rr ;
100
101 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
102 printf("3dDimize dataset\n"); PRINT_COMPILE_DATE ; exit(0);
103 }
104
105 mainENTRY("3dDimize"); machdep();
106
107 inset = THD_open_dataset(argv[1]); CHECK_OPEN_ERROR((inset,argv[1]);
108 DSET_load(inset); CHECK_LOAD_ERROR(inset) ;
109 medim = THD_median_brick(inset); if( medim == NULL ) ERROR_exit("Can't median");
110 DSET_unload(inset) ;
111 clip_val = THD_cliplevel(medim,0.50f) ;
112
113 /* create mask of values above clip value */
114
115 nx = medim->nx; ny = medim->ny; nz = medim->nz; nxy = nx*ny; nxyz = nxy*nz;
116 mar = MRI_FLOAT_PTR(medim) ;
117 mim = mri_new_conforming(medim,MRI_byte); mmm = MRI_BYTE_PTR(mim);
118 for( ii=0 ; ii < nxyz ; ii++ ) mmm[ii] = (mar[ii] >= clip_val) ;
119 mri_free(medim) ;
120 MRI_COPY_AUX(mim,DSET_BRICK(inset,0)) ;
121
122 THD_mask_clust( nx,ny,nz, mmm ) ;
123 THD_mask_erode( nx,ny,nz, mmm, 1 ) ;
124 THD_mask_clust( nx,ny,nz, mmm ) ;
125
126 outset = EDIT_empty_copy(inset) ;
127 EDIT_dset_items( outset ,
128 ADN_nvals , NRAD-1 ,
129 ADN_ntt , NRAD-1 ,
130 ADN_prefix , prefix ,
131 ADN_none ) ;
132 for( rr=0 ; rr < NRAD-1 ; rr++ ){
133 bb[rr] = (float *)calloc(sizeof(float),nxyz) ;
134 EDIT_substitute_brick( outset , rr , MRI_float , bb[rr] ) ;
135 }
136
137 for( kk=0 ; kk < nz ; kk++ ){
138 for( jj=0 ; jj < ny ; jj++ ){
139 for( ii=0 ; ii < nx ; ii++ ){
140 ijk = IJK(ii,jj,kk) ; if( mmm[ijk] == 0 ) continue ;
141 rr = mri_mask_localcount( mim , ii,jj,kk , NRAD,rad , ccc,vvv ) ;
142 if( rr > 0 ){
143 for( rr=0 ; rr < NRAD-1 ; rr++ )
144 bb[rr][ijk] = 3.0f * logf( ccc[rr+1]/(float)ccc[rr] )
145 / logf( vvv[rr+1]/(float)vvv[rr] ) ;
146 }
147 }}}
148
149 DSET_write( outset ) ;
150 WROTE_DSET( outset ) ;
151 exit(0) ;
152 }
153