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