1 #include "mrilib.h"
2 #include "parser.h"
3 
4 floatvecvec * symmetric_semi_rCDF( THD_3dim_dataset *dset ,
5                                    byte *mask , float top , int nbin ) ;
6 
main(int argc,char * argv[])7 int main( int argc , char *argv[] )
8 {
9    int nopt=1 , nvox ;
10    char *prefix="Normalizer" ;
11    THD_3dim_dataset *samset=NULL, *inset=NULL, *outset=NULL ;
12    byte *mask=NULL ; int nmask,nmask_hits ;
13    float apar,bpar,cpar,dpar ;
14    floatvec *rfv , *pfv ;
15 
16    /*-----------------------------------------------------------------------*/
17 
18    if( argc < 3 || strcasecmp(argv[1],"-help") == 0 ){
19      printf("\n"
20       "Usage: 3dNormalizer [options]\n"
21       "\n"
22       "This program builds the CDF of the '-sample' dataset, then\n"
23       "uses that to normalize the distribution of the '-input' dataset.\n"
24       "\n"
25       "Options\n"
26       "-------\n"
27       " -sample sss  = Read dataset 'sss' as the sample.  MANDATORY\n"
28       " -input  iii  = Read dataset 'iii' as the input.   MANDATORY\n"
29       " -mask   mmm  = Read dataset 'iii' as the mask.\n"
30       " -prefix ppp  = Use 'ppp' for the output datset prefix.\n"
31       "\n"
32       "Currently highly experimental! -- Zhark the nonGaussian\n"
33       "\n"
34      ) ;
35      exit(0) ;
36    }
37 
38    /*-----------------------------------------------------------------------*/
39 
40    while( nopt < argc && argv[nopt][0] == '-' ){
41 
42      /*---*/
43 
44      if( strcasecmp(argv[nopt],"-prefix") == 0 ){
45        if( ++nopt >= argc )
46          ERROR_exit("Need argument after '%s'",argv[nopt-1]) ;
47        prefix = strdup(argv[nopt]) ;
48        if( !THD_filename_ok(prefix) )
49          ERROR_exit("Prefix '%s' is not acceptable",prefix) ;
50        nopt++ ; continue ;
51      }
52 
53      /*---*/
54 
55      if( strcasecmp(argv[nopt],"-sample") == 0 ){
56        if( ++nopt >= argc )
57          ERROR_exit("Need argument after '%s'",argv[nopt-1]) ;
58        if( samset != NULL )
59          ERROR_exit("Can't use '-sample' more than once!") ;
60        samset = THD_open_dataset(argv[nopt]) ;
61        CHECK_OPEN_ERROR(samset,argv[nopt]) ;
62        nopt++ ; continue ;
63      }
64 
65      /*---*/
66 
67      if( strcasecmp(argv[nopt],"-input") == 0 ){
68        if( ++nopt >= argc )
69          ERROR_exit("Need argument after '%s'",argv[nopt-1]) ;
70        if( inset != NULL )
71          ERROR_exit("Can't use '-input' more than once!") ;
72        inset = THD_open_dataset(argv[nopt]) ;
73        CHECK_OPEN_ERROR(inset,argv[nopt]) ;
74        nopt++ ; continue ;
75      }
76 
77      /*---*/
78 
79      if( strcasecmp(argv[nopt],"-mask") == 0 ){
80        bytevec *bvec ;
81        if( mask != NULL )
82          ERROR_exit("Can't use '-mask' twice!") ;
83        if( ++nopt >= argc )
84          ERROR_exit("Need argument after '%s'",argv[nopt-1]) ;
85        bvec = THD_create_mask_from_string(argv[nopt]) ;
86        if( bvec == NULL )
87          ERROR_exit("Can't create mask from '-mask' option") ;
88        mask = bvec->ar ; nmask = bvec->nar ;
89        nmask_hits = THD_countmask( nmask , mask ) ;
90        if( nmask_hits > 0 )
91          INFO_message("%d voxels in -mask dataset",nmask_hits) ;
92        else
93          ERROR_exit("no nonzero voxels in -mask dataset") ;
94        nopt++ ; continue ;
95      }
96 
97      /*---*/
98 
99      ERROR_exit("Unknown option '%s'",argv[nopt]) ;
100 
101    }
102 
103    /*-----------------------------------------------------------------------*/
104 
105    if( samset == NULL ) ERROR_exit("'-sample' is mandatory") ;
106    if( inset  == NULL ) ERROR_exit("'-input' is mandatory") ;
107 
108    nvox = DSET_NVOX(samset) ;
109 
110    if( DSET_NVOX(samset) != DSET_NVOX(inset) )
111      ERROR_exit("'-sample' and '-inset' datasets don't match in voxel count!") ;
112 
113    if( mask != NULL && DSET_NVOX(samset) != nmask )
114      ERROR_exit("'-sample' and '-mask' datasets don't match in voxel count!") ;
115 
116    if( mask == NULL ){
117      mask = (byte *)malloc(sizeof(byte)*nvox) ;
118      memset( mask , 1 , sizeof(byte)*nvox ) ;
119      nmask = nmask_hits = nvox ;
120      INFO_message("no -mask option ==> using all %d voxels",nmask) ;
121    }
122 
123    /*-----------------------------------------------------------------------*/
124 
125    { floatvecvec *ovv ;
126      DSET_load(samset) ; CHECK_LOAD_ERROR(samset) ;
127      ovv = symmetric_semi_rCDF( samset , mask , 5.0f , 100 ) ;
128      rfv = ovv->fvar + 0 ;
129      pfv = ovv->fvar + 1 ;
130 #if 1
131      mri_write_floatvec( modify_afni_prefix(prefix,NULL,".cdf.1D") , rfv ) ;
132      mri_write_floatvec( modify_afni_prefix(prefix,NULL,".pdf.1D") , pfv ) ;
133 #endif
134      DSET_unload(samset) ;
135    }
136 
137 #undef  LNCOSH
138 #define LNCOSH(x) (fabsf(x)+logf(0.5f+0.5f*expf(-2.0f*fabsf(x))))
139 #undef  HFUNC
140 #define HFUNC(x)  (bpar*(x)+apar*(LNCOSH(dpar*(x)-cpar)-LNCOSH(cpar)))
141 
142    { float *qv, *wv, *xv , *fitv, parbot[26],partop[26],parout[26] ;
143      int ii, nval=rfv->nar ; float dx=rfv->dx ;
144      qv = (float *)malloc(sizeof(float)*nval) ;
145      wv = (float *)malloc(sizeof(float)*nval) ;
146      xv = (float *)malloc(sizeof(float)*nval) ;
147      for( ii=0 ; ii < nval ; ii++ ){
148        qv[ii] = qginv(0.5*rfv->ar[ii]) ;
149             if( qv[ii] <  1.5f ) wv[ii] = 0.7f ;
150        else if( qv[ii] <= 4.0f ) wv[ii] = 1.0f ;
151        else                      wv[ii] = 0.02f ;
152        xv[ii] = ii*dx ;
153        qv[ii] = qv[ii] - xv[ii] ;
154      }
155      parbot[0] =  0.0f ; partop[0] = 2.0f ; /* limits on a */
156      parbot[1] = -0.5f ; partop[1] = 0.5f ; /* limits on b */
157      parbot[2] =  0.1f ; partop[2] = 2.9f ; /* limits on c */
158      parbot[3] =  0.2f ; partop[3] = 2.2f ; /* limits on c */
159 #if 0
160      powell_set_verbose(2) ;
161 #endif
162      fitv = PARSER_fitter( nval , xv , qv ,
163                            "b*x+a*(logcosh(d*x-c)-logcosh(c))" , "x" ,
164                            parbot , partop , parout , 1 , wv ) ;
165      if( fitv == NULL )
166        ERROR_exit("PARSER_fitter() fails :-(") ;
167      apar = parout[0] ;
168      bpar = parout[1] ;
169      cpar = parout[2] ;
170      dpar = parout[3] ;
171 #if 1
172      INFO_message("apar=%g  bpar=%g  cpar=%g  dpar=%g",apar,bpar,cpar,dpar) ;
173      { MRI_IMAGE *qim = mri_new(nval,2,MRI_float) ;
174        float *qar = MRI_FLOAT_PTR(qim) ;
175        for( ii=0 ; ii < nval ; ii++ ){
176          qar[ii] = qv[ii] ; qar[ii+nval] = fitv[ii] ;
177        }
178        mri_write_1D( modify_afni_prefix(prefix,NULL,".qfit.1D") , qim ) ;
179        mri_free(qim) ;
180      }
181 #endif
182    }
183 
184    /*-----------------------------------------------------------------------*/
185 
186    { MRI_IMAGE *bim ; int nv=DSET_NVALS(inset),iv,jj ; float *bar ;
187 
188      DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;
189      bpar = bpar + 1.0f ;
190 
191      outset = EDIT_empty_copy(inset) ;
192      EDIT_dset_items( outset ,
193                         ADN_prefix    , prefix ,
194                         ADN_datum_all , MRI_float ,
195                       ADN_none ) ;
196      tross_Copy_History(inset,outset) ;
197      tross_Make_History("3dNormalizer",argc,argv,outset) ;
198 
199      for( iv=0 ; iv < nv ; iv++ ){
200        bim = THD_extract_float_brick(iv,inset) ; bar = MRI_FLOAT_PTR(bim) ;
201        if( DSET_BRICK_STATCODE(inset,iv) == FUNC_ZT_TYPE ){
202          for( jj=0 ; jj < nvox ; jj++ ){
203            bar[jj] = HFUNC(bar[jj]) ;
204          }
205        }
206        EDIT_substitute_brick( outset , iv , MRI_float , bar ) ;
207      }
208 
209      DSET_write(outset) ; WROTE_DSET(outset) ;
210    }
211 
212    exit(0) ;
213 }
214 
215 
216 /*--------------------------------------------------------------------------*/
217 
symmetric_semi_rCDF(THD_3dim_dataset * dset,byte * mask,float top,int nbin)218 floatvecvec * symmetric_semi_rCDF( THD_3dim_dataset *dset ,
219                                    byte *mask , float top , int nbin )
220 {
221    floatvec *rfv=NULL , *pfv=NULL ; floatvecvec *ovv=NULL ;
222    int64vec *riv=NULL ;
223    float dx=top/nbin , dxinv=nbin/top , val , scl ;
224    int ival , ii,jj , nvox,nval ; int64_t ntot , nsum ;
225    MRI_IMAGE *bim ; float *bar ;
226 
227 ENTRY("symmetric_semi_rCDF") ;
228 
229    MAKE_int64vec( riv , nbin+1 ) ;
230 
231    nvox = DSET_NVOX(dset) ; nval = DSET_NVALS(dset) ;
232 
233    ntot=0 ;
234    for( ival=0 ; ival < nval ; ival++ ){
235      bim = THD_extract_float_brick( ival , dset ) ;
236      if( bim == NULL ) continue ;
237      bar = MRI_FLOAT_PTR(bim) ;
238      for( ii=0 ; ii < nvox ; ii++ ){
239        val = fabsf(bar[ii]) ;
240        if( val > 0.0f && val < top && mask[ii] ){
241          jj = (int)(val*dxinv) ; riv->ar[jj]++ ; ntot++ ;
242        }
243      }
244    }
245 
246    if( ntot == 0 ){
247      KILL_int64vec(riv) ; RETURN(NULL) ;
248    }
249 
250    scl = 1.0f / ntot ;
251 
252    ovv = (floatvecvec *)malloc(sizeof(floatvecvec)) ;
253    ovv->nvec = 2 ;
254    ovv->fvar = (floatvec *)malloc(sizeof(floatvec)*2) ;
255    ovv->fvar[0].nar = nbin+1 ;
256    ovv->fvar[0].ar  = (float *)calloc(sizeof(float),(nbin+1)) ;
257    ovv->fvar[1].nar = nbin+1 ;
258    ovv->fvar[1].ar  = (float *)calloc(sizeof(float),(nbin+1)) ;
259    rfv = ovv->fvar + 0 ; rfv->dx = dx ; rfv->x0 = 0.0f ;
260    pfv = ovv->fvar + 1 ; pfv->dx = dx ; pfv->x0 = 0.0f ;
261 
262    nsum = 0 ;
263    for( jj=nbin ; jj >=0 ; jj-- ){
264      nsum += riv->ar[jj] ;
265      rfv->ar[jj] = scl * nsum ;
266      pfv->ar[jj] = riv->ar[jj] * scl ;
267    }
268 
269    KILL_int64vec(riv) ;
270    RETURN(ovv) ;
271 }
272