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