1 #include "mrilib.h"
2
3 /*--------------------------------------------------------------------------*/
4 static float hbot1 = 1.0f ;
5 static float htop1 = -1.0f ;
6 static float hbot2 = 1.0f ;
7 static float htop2 = -1.0f ;
8
mri_nbistat_setclip(float hb1,float ht1,float hb2,float ht2)9 void mri_nbistat_setclip( float hb1, float ht1 , float hb2, float ht2 )
10 {
11 hbot1 = hb1 ; htop1 = ht1 ; hbot2 = hb2 ; htop2 = ht2 ;
12 }
13
14 /*--------------------------------------------------------------------------*/
15 static MRI_IMAGE *wim = NULL ;
16 static MRI_IMAGE *wnim = NULL ;
17
mri_bistat_setweight(MRI_IMAGE * wm)18 void mri_bistat_setweight( MRI_IMAGE *wm ) /* 14 Aug 2007 */
19 {
20 if( wim != NULL ){ mri_free(wim); wim = NULL; }
21 if( wm != NULL ){ wim = mri_to_float(wm); }
22 return ;
23 }
24
25 /*--------------------------------------------------------------------------*/
26 /*! Input = 2 1D images, and an NBISTAT_ code to compute some statistic.
27 Output = statistic's value.
28 *//*------------------------------------------------------------------------*/
29
mri_nbistat(int code,MRI_IMAGE * im,MRI_IMAGE * jm)30 float mri_nbistat( int code , MRI_IMAGE *im , MRI_IMAGE *jm )
31 {
32 MRI_IMAGE *fim , *gim ;
33 float *far , *gar ; float outval=0.0f ;
34 int npt , ii ;
35 static int lref=0 ; static float **ref=NULL ; /* 26 Apr 2012 */
36
37 if( im == NULL || jm == NULL || im->nvox == 0 || im->nvox != jm->nvox )
38 return outval ;
39
40 /* convert input to float format, if not already there */
41
42 fim = (im->kind == MRI_float) ? im : mri_to_float(im) ;
43 gim = (jm->kind == MRI_float) ? jm : mri_to_float(jm) ;
44 far = MRI_FLOAT_PTR(fim) ; /* array of values to statisticate */
45 gar = MRI_FLOAT_PTR(gim) ;
46 npt = fim->nvox ; /* number of values */
47
48 if( hbot1 < htop1 ){
49 for( ii=0 ; ii < npt ; ii++ )
50 if( far[ii] < hbot1 ) far[ii] = hbot1 ;
51 else if( far[ii] > htop1 ) far[ii] = htop1 ;
52 }
53 if( hbot2 < htop2 ){
54 for( ii=0 ; ii < npt ; ii++ )
55 if( gar[ii] < hbot2 ) gar[ii] = hbot2 ;
56 else if( gar[ii] > htop2 ) gar[ii] = htop2 ;
57 }
58
59 #pragma omp critical
60 { if( (code == NBISTAT_L2SLOPE || code == NBISTAT_L1SLOPE) && lref < npt ){
61 if( ref == NULL ){
62 ref = (float **)malloc(sizeof(float *)*2) ; ref[0] = NULL ;
63 }
64 ref[0] = (float *)realloc(ref[0],sizeof(float)*npt) ; lref = npt ;
65 for( ii=0 ; ii < npt ; ii++ ) ref[0][ii] = 1.0f ;
66 }
67 }
68
69 switch( code ){
70
71 case NBISTAT_NUM: outval = (float)npt ; break ; /* quite easy */
72
73 case NBISTAT_SPEARMAN_CORR:
74 outval = THD_spearman_corr( npt , far , gar ) ; break ;
75
76 case NBISTAT_QUADRANT_CORR:
77 outval = THD_quadrant_corr( npt , far , gar ) ; break ;
78
79 case NBISTAT_PEARSON_CORR:
80 if( wnim == NULL )
81 outval = THD_pearson_corr( npt , far , gar ) ;
82 else
83 outval = THD_pearson_corr_wt( npt , far , gar , MRI_FLOAT_PTR(wnim) ) ;
84 break ;
85
86 case NBISTAT_MUTUAL_INFO:
87 outval = THD_mutual_info( npt , far , gar ) ; break ;
88
89 case NBISTAT_NORMUT_INFO:
90 outval = THD_norm_mutinf( npt , far , gar ) ;
91 if( outval != 0.0f ) outval = 1.0f / outval ;
92 break ;
93
94 case NBISTAT_JOINT_ENTROPY:
95 outval = THD_jointentrop( npt , far , gar ) ; break ;
96
97 case NBISTAT_HELLINGER:
98 outval = THD_hellinger( npt , far , gar ) ; break ;
99
100 case NBISTAT_CORR_RATIO_M:
101 THD_corr_ratio_mode(1) ;
102 outval = THD_corr_ratio( npt , far , gar ) ; break ;
103
104 case NBISTAT_CORR_RATIO_A:
105 THD_corr_ratio_mode(2) ;
106 outval = THD_corr_ratio( npt , far , gar ) ; break ;
107
108 case NBISTAT_CORR_RATIO_U:
109 THD_corr_ratio_mode(0) ;
110 outval = THD_corr_ratio( npt , far , gar ) ; break ;
111
112 case NBISTAT_NCD:
113 outval = THD_ncdfloat( npt , far , gar ) ; break ;
114
115 case NBISTAT_L2SLOPE:{
116 float *qfit ;
117 ref[1] = far ;
118 qfit = lsqfit( npt , gar , NULL , 2 , ref ) ;
119 if( qfit != NULL ){ outval = qfit[1] ; free(qfit) ; }
120 }
121 break ;
122
123 case NBISTAT_L1SLOPE:{
124 float qfit[2] , val ;
125 ref[1] = far ;
126 val = cl1_solve( npt , 2 , gar , ref , qfit , 0 ) ;
127 if( val >= 0.0f ) outval = qfit[1] ;
128 }
129 break ;
130
131 case NBISTAT_EUCLIDIAN_DIST: /* May 4 2012 ZSS */
132 outval = THD_distance( npt , far , gar , 0) ; break ;
133
134 case NBISTAT_CITYBLOCK_DIST: /* May 4 2012 ZSS */
135 outval = THD_distance( npt , far , gar , 1) ; break ;
136
137 }
138
139 /* cleanup and exit */
140
141 if( fim != im ) mri_free(fim) ;
142 if( gim != jm ) mri_free(gim) ;
143 return outval ;
144 }
145
146 /*--------------------------------------------------------------------------*/
147 /*! Compute a local statistic at each voxel of an image pair, possibly with
148 a mask; 'local' is defined with a neighborhood; 'statistic' is defined
149 by an NBISTAT_ code.
150 *//*------------------------------------------------------------------------*/
151
mri_localbistat(MRI_IMAGE * im,MRI_IMAGE * jm,byte * mask,MCW_cluster * nbhd,int code)152 MRI_IMAGE * mri_localbistat( MRI_IMAGE *im, MRI_IMAGE *jm ,
153 byte *mask, MCW_cluster *nbhd, int code )
154 {
155 MRI_IMAGE *outim , *nbim , *nbjm ;
156 float *outar ;
157 int ii,jj,kk , nx,ny,nz , ijk ;
158
159 ENTRY("mri_localbistat") ;
160
161 if( im == NULL || nbhd == NULL ) RETURN(NULL) ;
162
163 outim = mri_new_conforming( im , MRI_float ) ;
164 outar = MRI_FLOAT_PTR(outim) ;
165 nx = outim->nx ; ny = outim->ny ; nz = outim->nz ;
166
167 ijk = (int)cbrt((double)nbhd->num_pt) ; /* for entropy, etc. */
168 set_2Dhist_hbin( ijk ) ;
169
170 for( ijk=kk=0 ; kk < nz ; kk++ ){
171 for( jj=0 ; jj < ny ; jj++ ){
172 for( ii=0 ; ii < nx ; ii++ ){
173 nbim = mri_get_nbhd( im , mask , ii,jj,kk , nbhd ) ;
174 nbjm = mri_get_nbhd( jm , mask , ii,jj,kk , nbhd ) ;
175 outar[ijk++] = mri_nbistat( code , nbim , nbjm ) ;
176 mri_free(nbim) ; mri_free(nbjm) ;
177 }}}
178
179 RETURN(outim) ;
180 }
181
182 /*--------------------------------------------------------------------------*/
183
184 static int verb=0 , vn=0 ;
THD_localbistat_verb(int i)185 void THD_localbistat_verb(int i){ verb=i; vn=0; }
186
vstep_print(void)187 static void vstep_print(void)
188 {
189 static char xx[10] = "0123456789" ;
190 fprintf(stderr , "%c" , xx[vn%10] ) ;
191 if( vn%10 == 9) fprintf(stderr,".") ;
192 vn++ ;
193 }
194
195 /*--------------------------------------------------------------------------*/
196
THD_localbistat(THD_3dim_dataset * dset,THD_3dim_dataset * eset,byte * mask,MCW_cluster * nbhd,int ncode,int * code)197 THD_3dim_dataset * THD_localbistat( THD_3dim_dataset *dset ,
198 THD_3dim_dataset *eset , byte *mask ,
199 MCW_cluster *nbhd , int ncode, int *code )
200 {
201 THD_3dim_dataset *oset ;
202 MRI_IMAGE *nbim , *nbjm ;
203 int iv,cc , nvin,nvout , nx,ny,nz,nxyz , ii,jj,kk,ijk ;
204 float **aar ;
205 int vstep , iiv , jjv , nvd,nve ;
206
207 ENTRY("THD_localbistat") ;
208
209 if( dset == NULL || eset == NULL ||
210 nbhd == NULL || ncode < 1 || code == NULL ) RETURN(NULL);
211
212 if( DSET_NVOX(dset) != DSET_NVOX(eset) ) RETURN(NULL) ;
213 DSET_load(dset) ; if( !DSET_LOADED(dset) ) RETURN(NULL) ;
214 DSET_load(eset) ; if( !DSET_LOADED(eset) ) RETURN(NULL) ;
215
216 oset = EDIT_empty_copy( dset ) ;
217 nvd = DSET_NVALS( dset ) ;
218 nve = DSET_NVALS( eset ) ;
219 nvin = MAX(nvd,nve) ;
220 nvout = nvin * ncode ;
221 EDIT_dset_items( oset ,
222 ADN_nvals , nvout ,
223 ADN_datum_all , MRI_float ,
224 ADN_ntt , nvout ,
225 ADN_nsl , 0 ,
226 ADN_brick_fac , NULL ,
227 ADN_prefix , "localbistat" ,
228 ADN_none ) ;
229
230 nx = DSET_NX(dset) ;
231 ny = DSET_NY(dset) ;
232 nz = DSET_NZ(dset) ; nxyz = nx*ny*nz ;
233
234 vstep = (verb && nxyz > 66666) ? nxyz/50 : 0 ;
235 if( vstep ) fprintf(stderr,"++ voxel loop:") ;
236
237 aar = (float **)malloc(sizeof(float *)*ncode) ;
238
239 for( iv=0 ; iv < nvin ; iv++ ){
240 iiv = iv ; if( iiv >= nvd ) iiv = nvd-1 ; /* sub-brick for dset */
241 jjv = iv ; if( jjv >= nve ) jjv = nve-1 ; /* sub-brick for eset */
242 for( cc=0 ; cc < ncode ; cc++ ){
243 aar[cc] = (float *)malloc(sizeof(float)*nxyz) ;
244 if( aar[cc] == NULL )
245 ERROR_exit("THD_localbistat: out of memory at iv=%d cc=%d",iv,cc);
246 }
247
248 for( ijk=kk=0 ; kk < nz ; kk++ ){
249 for( jj=0 ; jj < ny ; jj++ ){
250 for( ii=0 ; ii < nx ; ii++,ijk++ ){
251 if( vstep && ijk%vstep==vstep-1 ) vstep_print() ;
252 nbim = THD_get_dset_nbhd( dset,iiv , mask,ii,jj,kk , nbhd ) ;
253 nbjm = THD_get_dset_nbhd( eset,jjv , mask,ii,jj,kk , nbhd ) ;
254 if( wim != NULL ) wnim = mri_get_nbhd( wim , mask,ii,jj,kk , nbhd ) ;
255 for( cc=0 ; cc < ncode ; cc++ )
256 aar[cc][ijk] = mri_nbistat( code[cc] , nbim,nbjm ) ;
257 mri_free(nbim) ; mri_free(nbjm) ;
258 if( wnim != NULL ){ mri_free(wnim); wnim = NULL; }
259 }}}
260
261 if( iiv < nvd-1 ) DSET_unload_one(dset,iiv) ;
262 if( jjv < nve-1 ) DSET_unload_one(eset,jjv) ;
263
264 for( cc=0 ; cc < ncode ; cc++ )
265 EDIT_substitute_brick( oset , iv*ncode+cc , MRI_float , aar[cc] ) ;
266 }
267
268 if( vstep ) fprintf(stderr,"\n") ;
269 free((void *)aar) ;
270 RETURN(oset) ;
271 }
272