1 #include "mrilib.h"
2
3 /***---------------------------------------------------------***/
4 /*** This file is intended to be #include-d into mri_nwarp.c ***/
5 /***---------------------------------------------------------***/
6
7 /***
8 #ifdef USE_OMP
9 #include <omp.h>
10 #endif
11 ***/
12
13 #ifndef USE_OMP
omp_get_thread_num(void)14 static int omp_get_thread_num(void) { return 0 ; }
omp_get_max_threads(void)15 static int omp_get_max_threads(void){ return 1 ; }
16 #endif
17
18 typedef struct {
19 int meth ;
20 } INCOR_generic ;
21
22 typedef struct {
23 int meth ;
24 int npt ;
25 double sx , sxx , sy , syy , sxy , sw ;
26 } INCOR_pearson ;
27
28 typedef struct {
29 int meth ;
30 int npt ;
31 double sx , sxx , sy , syy , sxy , sw ;
32 double xcbot , xctop , ycbot , yctop ;
33 double xdbot , xdtop , ydbot , ydtop ;
34 } INCOR_pearclp ;
35
36 typedef struct {
37 int meth ;
38 int nbin ;
39 float *xc , *yc , *xyc , nww ;
40 float xxbot , xxtop , yybot , yytop ;
41 float xcbot , xctop , ycbot , yctop ;
42 } INCOR_2Dhist ;
43
44 typedef struct {
45 int numblok , *nelm ;
46 int nx,ny,nz, *blkn ;
47 int *nsum ; float *sx,*sxx, *sy,*syy, *sxy,*sw , *pval ;
48 float psum , wsum ; int ninserted ;
49 } INCORR_BLOK_set ; /* 25 Jun 2014 */
50
51 typedef struct {
52 int meth ;
53 int iibot,iitop , jjbot,jjtop , kkbot,kktop ;
54 int nii , njj , nkk ;
55 INCORR_BLOK_set *ibs ;
56 } INCOR_localpearson ; /* 25 Jun 2014 */
57
58 #define ALLOW_DEBUG_LPC
59 #ifdef ALLOW_DEBUG_LPC
60 static int debug_lpc = 0 ;
61 #else
62 #define debug_lpc 0
63 #endif
64
65 #undef INCOR_methcode
66 #define INCOR_methcode(vp) ( ((vp) != NULL) ? ((INCOR_generic *)vp)->meth : 0 )
67
68 #undef MYatanh
69 #define MYatanh(x) ( ((x)<-0.9993293) ? -4.0 \
70 :((x)>+0.9993293) ? +4.0 : atanh(x) )
71
72 /****************************************************************************/
73 /*** Local Pearson correlation [25 Jun 2014 -- at long last!] ***/
74 /****************************************************************************/
75
76 #undef MINCOR
77 #define MINCOR 9 /* min number of points to correlate */
78
79 #undef CMAX
80 #define CMAX 0.99f
81
82 /* is abs(a) <= s ?? */
83
84 #undef FAS
85 #define FAS(a,s) ( (a) <= (s) && (a) >= -(s) )
86
87 /** define inside of a ball; is point (a,b,c) inside? **/
88 /** volume of ball = 4*PI/3 * siz**3 = 4.1888 * siz**3 **/
89
90 #define IB_BLOK_inside_ball(a,b,c,siz) \
91 ( ((a)*(a)+(b)*(b)+(c)*(c)) <= (siz) )
92
93 /** define inside of a cube **/
94 /** volume of cube = 8 * siz**3 **/
95 /** lattice vectors = [2*siz,0,0] [0,2*siz,0] [0,0,2*siz] **/
96
97 #define IB_BLOK_inside_cube(a,b,c,siz) \
98 ( FAS((a),(siz)) && FAS((b),(siz)) && FAS((c),(siz)) )
99
100 /** define inside of a rhombic dodecahedron (RHDD) **/
101 /** volume of RHDD = 2 * siz**3 **/
102 /** lattice vectors = [siz,siz,0] [0,siz,siz] [siz,0,siz] **/
103
104 #define IB_BLOK_inside_rhdd(a,b,c,siz) \
105 ( FAS((a)+(b),(siz)) && FAS((a)-(b),(siz)) && \
106 FAS((a)+(c),(siz)) && FAS((a)-(c),(siz)) && \
107 FAS((b)+(c),(siz)) && FAS((b)-(c),(siz)) )
108
109 /** define inside of a truncated octahedron (TOHD) **/
110 /** volume of TOHD = 4 * siz**3 **/
111 /** lattice vectors = [-siz,siz,siz] [siz,-siz,siz] [siz,siz,-siz] **/
112
113 #define IB_BLOK_inside_tohd(a,b,c,siz) \
114 ( FAS((a),(siz)) && FAS((b),(siz)) && FAS((c),(siz)) && \
115 FAS((a)+(b)+(c),1.5f*(siz)) && FAS((a)-(b)+(c),1.5f*(siz)) && \
116 FAS((a)+(b)-(c),1.5f*(siz)) && FAS((a)-(b)-(c),1.5f*(siz)) )
117
118 /** define inside of an arbitrary blok type **/
119
120 #define IB_BLOK_inside(bt,a,b,c,s) \
121 ( ((bt)==GA_BLOK_BALL) ? IB_BLOK_inside_ball((a),(b),(c),(s)) \
122 : ((bt)==GA_BLOK_CUBE) ? IB_BLOK_inside_cube((a),(b),(c),(s)) \
123 : ((bt)==GA_BLOK_RHDD) ? IB_BLOK_inside_rhdd((a),(b),(c),(s)) \
124 : ((bt)==GA_BLOK_TOHD) ? IB_BLOK_inside_tohd((a),(b),(c),(s)) \
125 : 0 )
126
127 /** add 1 value to a dynamically allocated integer array **/
128
129 #define IB_BLOK_ADDTO_intar(nar,nal,ar,val) \
130 do{ if( (nar) == (nal) ){ \
131 (nal) = 1.5*(nal)+16; (ar) = (int *)realloc((ar),sizeof(int)*(nal)); \
132 } \
133 (ar)[(nar)++] = (val); \
134 } while(0)
135
136 /** truncate dynamically allocated integer array down to size **/
137
138 #define IB_BLOK_CLIP_intar(nar,nal,ar) \
139 do{ if( (nar) < (nal) && (nar) > 0 ){ \
140 (nal) = (nar); (ar) = (int *)realloc((ar),sizeof(int)*(nal)); \
141 }} while(0)
142
143 /*----------------------------------------------------------------------------*/
144 /*! Fill a struct with list of points contained in sub-bloks of the base.
145
146 - nx,ny,nz = 3D grid dimensions
147 - dx,dy,dz = 3D grid spacings
148 - mask = byte mask of points to use (can be NULL)
149 - bloktype = one of GA_BLOK_BALL, GA_BLOK_CUBE, GA_BLOK_RHDD, GA_BLOK_TOHD
150 - blokrad = radius parameter for the bloks to be built
151 - minel = minimum number of points to put in a blok
152 (if 0, function will pick a value)
153 *//*--------------------------------------------------------------------------*/
154
create_INCORR_BLOK_set(int nx,int ny,int nz,float dx,float dy,float dz,byte * mask,int bloktype,float blokrad,int minel)155 INCORR_BLOK_set * create_INCORR_BLOK_set(
156 int nx , int ny , int nz ,
157 float dx , float dy , float dz ,
158 byte *mask ,
159 int bloktype , float blokrad , int minel )
160 {
161 INCORR_BLOK_set *ibs ;
162 float dxp,dyp,dzp , dxq,dyq,dzq , dxr,dyr,dzr , xt,yt,zt ;
163 float xx,yy,zz , uu,vv,ww , siz ;
164 THD_mat33 latmat , invlatmat ; THD_fvec3 pqr , xyz ;
165 int pb,pt , qb,qt , rb,rt , pp,qq,rr , tnblok,nball , ii , nxy,nxyz ;
166 int aa,bb,cc , dd,ss , np,nq,nr,npq , *nelm,*nalm,**elm ;
167 int ntot,nsav,nblok , sgood , maxel ;
168 int *ilist ;
169 const float shfac = 1.0f ;
170
171 ENTRY("create_INCORR_BLOK_set") ;
172
173 if( nx < 3 || ny < 3 || nz < 1 ) RETURN(NULL) ;
174 if( dx <= 0.0f ) dx = 1.0f ;
175 if( dy <= 0.0f ) dy = 1.0f ;
176 if( dz <= 0.0f ) dz = 1.0f ;
177
178 /* Create lattice vectors to generate translated bloks:
179 The (p,q,r)-th blok -- for integral p,q,r -- is at (x,y,z) offset
180 (dxp,dyp,dzp)*p + (dxq,dyq,dzq)*q + (dxr,dyr,dzr)*r
181 Also set the 'siz' parameter for the blok, to test for inclusion. */
182
183 switch( bloktype ){
184
185 /* balls go on a hexagonal close packed lattice,
186 but with lattice spacing reduced to avoid gaps
187 (of course, then the balls overlap -- c'est la geometrie) */
188
189 case GA_BLOK_BALL:{
190 float s3=1.73205f , /* sqrt(3) */
191 s6=2.44949f , /* sqrt(6) */
192 a =blokrad*0.866025f ; /* shrink spacing to avoid gaps */
193 siz = blokrad*blokrad ;
194 /* hexagonal close packing basis vectors for sphere of radius a */
195 a *= shfac ;
196 dxp = 2.0f * a ; dyp = 0.0f ; dzp = 0.0f ;
197 dxq = a ; dyq = a * s3; dzq = 0.0f ;
198 dxr = a ; dyr = a / s3; dzr = a * 0.666667f*s6 ;
199 }
200 break ;
201
202 /* cubes go on a simple cubical lattice, spaced so faces touch */
203
204 case GA_BLOK_CUBE:{
205 float a = blokrad ;
206 siz = a ; a *= shfac ;
207 dxp = 2*a ; dyp = 0.0f; dzp = 0.0f ;
208 dxq = 0.0f; dyq = 2*a ; dzq = 0.0f ;
209 dxr = 0.0f; dyr = 0.0f; dzr = 2*a ;
210 }
211 break ;
212
213 /* rhombic dodecahedra go on a FCC lattice,
214 spaced so that faces touch (i.e., no volumetric overlap) */
215
216 case GA_BLOK_RHDD:{
217 float a = blokrad ;
218 siz = a ; a *= shfac ;
219 dxp = a ; dyp = a ; dzp = 0.0f ;
220 dxq = 0.0f; dyq = a ; dzq = a ;
221 dxr = a ; dyr = 0.0f; dzr = a ;
222 }
223 break ;
224
225 /* truncated octahedra go on a BCC lattice,
226 spaced so that faces touch (i.e., no volumetric overlap) */
227
228 case GA_BLOK_TOHD:{
229 float a = blokrad ;
230 siz = a ; a *= shfac ;
231 dxp = -a ; dyp = a ; dzp = a ;
232 dxq = a ; dyq = -a ; dzq = a ;
233 dxr = a ; dyr = a ; dzr = -a ;
234 }
235 break ;
236
237 default: RETURN(NULL) ; /** should not happen! **/
238 }
239
240 /* find range of (p,q,r) indexes needed to cover volume,
241 by checking out all 7 corners besides (0,0,0) (where p=q=r=0) */
242
243 LOAD_MAT( latmat, dxp , dxq , dxr ,
244 dyp , dyq , dyr ,
245 dzp , dzq , dzr ) ; invlatmat = MAT_INV(latmat) ;
246
247 xt = (nx-1)*dx ; yt = (ny-1)*dy ; zt = (nz-1)*dz ;
248 pb = pt = qb = qt = rb = rt = 0 ; /* initialize (p,q,r) bot, top values */
249
250 LOAD_FVEC3(xyz , xt,0.0f,0.0f ); pqr = MATVEC( invlatmat , xyz ) ;
251 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
252 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
253 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
254
255 LOAD_FVEC3(xyz , xt,yt,0.0f ) ; pqr = MATVEC( invlatmat , xyz ) ;
256 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
257 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
258 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
259
260 LOAD_FVEC3(xyz , xt,0.0f,zt ) ; pqr = MATVEC( invlatmat , xyz ) ;
261 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
262 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
263 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
264
265 LOAD_FVEC3(xyz , xt,yt,zt ) ; pqr = MATVEC( invlatmat , xyz ) ;
266 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
267 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
268 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
269
270 LOAD_FVEC3(xyz , 0.0f,yt,0.0f ); pqr = MATVEC( invlatmat , xyz ) ;
271 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
272 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
273 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
274
275 LOAD_FVEC3(xyz , 0.0f,0.0f,zt ); pqr = MATVEC( invlatmat , xyz ) ;
276 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
277 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
278 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
279
280 LOAD_FVEC3(xyz , 0.0f,yt,zt ) ; pqr = MATVEC( invlatmat , xyz ) ;
281 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
282 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
283 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
284
285 /* Lattice index range is (p,q,r) = (pb..pt,qb..qt,rb..rt) inclusive */
286
287 np = pt-pb+1 ; /* number of p values to consider */
288 nq = qt-qb+1 ; npq = np*nq ;
289 nr = rt-rb+1 ;
290 tnblok = npq*nr ; /* total number of bloks to consider */
291
292 /* Now have list of bloks, so put points into each blok list */
293
294 nelm = (int *) calloc(sizeof(int) ,tnblok) ; /* # pts in each blok */
295 nalm = (int *) calloc(sizeof(int) ,tnblok) ; /* # malloc-ed in each blok */
296 elm = (int **)calloc(sizeof(int *),tnblok) ; /* list of pts in each blok */
297
298 nxy = nx*ny ; nxyz = nxy*nz ;
299
300 /* macro to test if (xx,yy,zz) goes into blok (ta,tb,tc), and put it there */
301
302 #undef TEST_BLOK_xyz
303 #define TEST_BLOK_xyz(ta,tb,tc) \
304 do{ LOAD_FVEC3( pqr , ta,tb,tc ) ; \
305 xyz = MATVEC( latmat , pqr ) ; \
306 uu = xx-xyz.xyz[0] ; vv = yy-xyz.xyz[1] ; ww = zz-xyz.xyz[2] ; \
307 if( IB_BLOK_inside( bloktype , uu,vv,ww , siz ) ){ \
308 dd = (ta-pb) + (tb-qb)*np + (tc-rb)*npq ; /* blok index */ \
309 IB_BLOK_ADDTO_intar( nelm[dd], nalm[dd], elm[dd], ss ) ; \
310 ntot++ ; sgood = 1 ; \
311 } \
312 } while(0) ;
313
314 /* loop over points to insert; each point goes in (at most) ONE blok */
315
316 for( ntot=ii=0 ; ii < nxyz ; ii++ ){
317 if( mask != NULL && !mask[ii] ) continue ; /* bad point */
318 pp = ii%nx ; rr = ii/nxy ; qq = (ii-rr*nxy)/nx ; /* 3D indexes (pp,qq,rr) */
319 sgood = 0 ; ss = ii ; /* index in 1D array */
320 xx = pp*dx ; yy = qq*dy ; zz = rr*dz ; /* xyz spatial coordinates */
321 LOAD_FVEC3( xyz , xx,yy,zz ) ;
322 pqr = MATVEC( invlatmat , xyz ) ; /* float lattice coords */
323 pp = (int)floorf(pqr.xyz[0]+.499f) ; /* integer lattice coords */
324 qq = (int)floorf(pqr.xyz[1]+.499f) ;
325 rr = (int)floorf(pqr.xyz[2]+.499f) ;
326 /* test if this blok holds the current point */
327 if( pp >= pb && pp <= pt && qq >= qb && qq <= qt && rr >= rb && rr <= rt ){
328 TEST_BLOK_xyz(pp,qq,rr) ; if( sgood ) goto TEST_BLOK_DONE ;
329 }
330 /* otherwise, search nearby bloks for inclusion of (xx,yy,zz) */
331 for( cc=rr-1 ; cc <= rr+1 ; cc++ ){
332 if( cc < rb || cc > rt ) continue ;
333 for( bb=qq-1 ; bb <= qq+1 ; bb++ ){
334 if( bb < qb || bb > qt ) continue ;
335 for( aa=pp-1 ; aa <= pp+1 ; aa++ ){
336 if( aa < pb || aa > pt ) continue ;
337 if( aa==pp && bb==qq && cc==rr ) continue ; /* already tested */
338 TEST_BLOK_xyz(aa,bb,cc) ; if( sgood ) goto TEST_BLOK_DONE ;
339 }
340 }
341 }
342 TEST_BLOK_DONE: /*nada*/ ;
343 }
344 #undef TEST_BLOK_xyz
345
346 /* compute the min number of points allowed per blok? */
347
348 if( minel < MINCOR ){
349 for( minel=dd=0 ; dd < tnblok ; dd++ ) minel = MAX(minel,nelm[dd]) ;
350 minel = (int)(0.321f*minel)+1 ; if( minel < MINCOR ) minel = MINCOR ;
351 }
352 for( maxel=dd=0 ; dd < tnblok ; dd++ ) maxel = MAX(maxel,nelm[dd]) ;
353
354 /* now cast out bloks that have too few points */
355
356 for( ntot=nblok=dd=0 ; dd < tnblok ; dd++ ){
357 if( nelm[dd] < minel ){
358 if( elm[dd] != NULL ){ free(elm[dd]); elm[dd] = NULL; }
359 nelm[dd] = 0 ; /* mark as an unsaved blok */
360 } else {
361 nblok++ ; /* count of saved bloks */
362 ntot += nelm[dd] ; /* count of saved points */
363 }
364 }
365 free(nalm) ;
366
367 if( nblok == 0 ){ /* didn't find any arrays to keep!? */
368 ERROR_message("create_INCORR_BLOK_set can't find bloks with at least %d nodes",minel);
369 free(nelm) ; free(elm) ; RETURN(NULL) ;
370 }
371
372 /* make a map from point index to (saved) blok index */
373
374 ilist = (int *)malloc(sizeof(int)*nxyz) ;
375 for( pp=0 ; pp < nxyz ; pp++ ) ilist[pp] = -666 ; /* mark as not in any blok */
376 for( rr=qq=dd=0 ; dd < tnblok ; dd++ ){
377 if( nelm[dd] == 0 ) continue ; /* not a saved blok */
378 for( pp=0 ; pp < nelm[dd] ; pp++ ) ilist[ elm[dd][pp] ] = rr ;
379 free(elm[dd]) ; /* no longer needed */
380 rr++ ; /* increment saved blok index */
381 }
382
383 /* create output struct */
384
385 ibs = (INCORR_BLOK_set *)malloc(sizeof(INCORR_BLOK_set)) ;
386 ibs->numblok = nblok ;
387 ibs->nelm = (int *) calloc(sizeof(int) ,nblok) ;
388 ibs->nsum = (int *) calloc(sizeof(int) ,nblok) ;
389 ibs->sx = (float *)calloc(sizeof(float),nblok) ;
390 ibs->sxx = (float *)calloc(sizeof(float),nblok) ;
391 ibs->sy = (float *)calloc(sizeof(float),nblok) ;
392 ibs->syy = (float *)calloc(sizeof(float),nblok) ;
393 ibs->sxy = (float *)calloc(sizeof(float),nblok) ;
394 ibs->sw = (float *)calloc(sizeof(float),nblok) ;
395 ibs->pval = (float *)calloc(sizeof(float),nblok) ;
396
397 ibs->nx = nx ; ibs->ny = ny ; ibs->nz = nz ; ibs->blkn = ilist ;
398 ibs->psum = ibs->wsum = 0.0f ; ibs->ninserted = 0 ;
399
400 for( rr=dd=0 ; dd < tnblok ; dd++ ){
401 if( nelm[dd] > 0 ) ibs->nelm[rr++] = nelm[dd] ;
402 }
403
404 free(nelm) ; free(elm) ;
405
406 RETURN(ibs) ;
407 }
408
409 /*---------------------------------------------------------------------------*/
410
destroy_INCORR_BLOK_set(INCORR_BLOK_set * ibs)411 void destroy_INCORR_BLOK_set( INCORR_BLOK_set *ibs )
412 {
413 if( ibs == NULL ) return ;
414
415 if( ibs->nelm != NULL ) free(ibs->nelm) ;
416 if( ibs->blkn != NULL ) free(ibs->blkn) ;
417 if( ibs->nsum != NULL ) free(ibs->nsum) ;
418 if( ibs->sx != NULL ) free(ibs->sx ) ;
419 if( ibs->sxx != NULL ) free(ibs->sxx ) ;
420 if( ibs->sy != NULL ) free(ibs->sy ) ;
421 if( ibs->syy != NULL ) free(ibs->syy ) ;
422 if( ibs->sxy != NULL ) free(ibs->sxy ) ;
423 if( ibs->sw != NULL ) free(ibs->sw ) ;
424 if( ibs->pval != NULL ) free(ibs->pval) ;
425 free(ibs) ; \
426 return ;
427 }
428
429 /*---------------------------------------------------------------------------*/
430
clear_INCORR_BLOK_set(INCORR_BLOK_set * ibs)431 void clear_INCORR_BLOK_set( INCORR_BLOK_set *ibs )
432 {
433 if( ibs == NULL ) return ;
434
435 AAmemset(ibs->nsum,0,sizeof(int) *ibs->numblok) ;
436 AAmemset(ibs->sx ,0,sizeof(float)*ibs->numblok) ;
437 AAmemset(ibs->sxx ,0,sizeof(float)*ibs->numblok) ;
438 AAmemset(ibs->sy ,0,sizeof(float)*ibs->numblok) ;
439 AAmemset(ibs->syy ,0,sizeof(float)*ibs->numblok) ;
440 AAmemset(ibs->sxy ,0,sizeof(float)*ibs->numblok) ;
441 AAmemset(ibs->sw ,0,sizeof(float)*ibs->numblok) ;
442 AAmemset(ibs->pval,0,sizeof(float)*ibs->numblok) ;
443 ibs->psum = ibs->wsum = 0.0f ; ibs->ninserted = 0 ;
444 return ;
445 }
446
447 /*---------------------------------------------------------------------------*/
448 /* Insert data into the struct. */
449
addto_INCORR_BLOK_set(INCORR_BLOK_set * ibs,int ibot,int itop,int jbot,int jtop,int kbot,int ktop,float * xval,float * yval,float * wt)450 void addto_INCORR_BLOK_set( INCORR_BLOK_set *ibs ,
451 int ibot , int itop , /* indexes are inclusive */
452 int jbot , int jtop ,
453 int kbot , int ktop ,
454 float *xval , float *yval , float *wt )
455 {
456 int ii,jj,kk , nx,ny,nz,nxy , pp,dd , nblok=ibs->numblok ;
457 int *nelm=ibs->nelm , *nsum=ibs->nsum , *blkn=ibs->blkn ;
458 float *sx = ibs->sx , *sxx = ibs->sxx , *sy = ibs->sy ,
459 *syy = ibs->syy , *sxy = ibs->sxy , *sw = ibs->sw ,
460 *pval = ibs->pval ;
461 double xx,yy,xy,ww ; float rval ;
462
463 nx = ibs->nx ; ny = ibs->ny ; nz = ibs->nz ; nxy = nx*ny ;
464
465 if( wt == NULL ){
466 for( pp=0,kk=kbot ; kk <= ktop ; kk++ ){
467 for( jj=jbot ; jj <= jtop ; jj++ ){
468 for( ii=ibot ; ii <= itop ; ii++,pp++ ){
469 dd = blkn[ii + jj*nx + kk*nxy] ;
470 if( dd >= 0 && nsum[dd] < nelm[dd] ){ /* in incomplete blok? */
471 xx = (double)xval[pp] ; yy = (double)yval[pp] ;
472 sx[dd] += xx ; sxx[dd] += xx*xx ;
473 sy[dd] += yy ; syy[dd] += yy*yy ; sxy[dd] += xx*yy ; sw[dd] += 1.0 ;
474 nsum[dd]++ ; ibs->ninserted++ ;
475 }
476 }}}
477 } else {
478 for( pp=0,kk=kbot ; kk <= ktop ; kk++ ){
479 for( jj=jbot ; jj <= jtop ; jj++ ){
480 for( ii=ibot ; ii <= itop ; ii++,pp++ ){
481 dd = blkn[ii + jj*nx + kk*nxy] ;
482 if( dd >= 0 && nsum[dd] < nelm[dd] ){ /* in incomplete blok? */
483 ww = wt[pp] ; if( ww <= 0.0 ) continue ; /* no weight? */
484 xx = (double)xval[pp] ; yy = (double)yval[pp] ;
485 sx[dd] += xx*ww; sxx[dd] += xx*xx*ww;
486 sy[dd] += yy*ww; syy[dd] += yy*yy*ww; sxy[dd] += xx*yy*ww; sw[dd] += ww;
487 nsum[dd]++ ; ibs->ninserted++ ;
488 }
489 }}}
490 }
491
492 /* find newly completed bloks and sum them up */
493
494 for( dd=0 ; dd < nblok ; dd++ ){
495 if( nsum[dd] == nelm[dd] && sw[dd] > 0.0f ){ /* now completed? */
496 ww = 1.0 / sw[dd] ;
497 xx = sxx[dd] - sx[dd] * sx[dd] * ww ;
498 yy = syy[dd] - sy[dd] * sy[dd] * ww ;
499 xy = sxy[dd] - sx[dd] * sy[dd] * ww ;
500 if( xx <= 0.0 || yy <= 0.0 ){
501 pval[dd] = 0.0f ;
502 } else {
503 rval = (float)(xy/sqrt(xx*yy)) ;
504 if( rval > CMAX ) rval = CMAX ;
505 else if( rval < -CMAX ) rval = -CMAX ;
506 pval[dd] = logf( (1.0f+rval)/(1.0f-rval) ) ;
507 ibs->psum += (float)sw[dd] * pval[dd] * fabsf(pval[dd]) ;
508 }
509 ibs->wsum += (float)sw[dd] ;
510 sw[dd] = 0.0f ; /* marked as completed and summed up */
511 }
512 }
513
514 return ;
515 }
516
517 /*---------------------------------------------------------------------------*/
518 /* Compute correlation, including the extra data, if any
519 (but this extra data is not stored into the struct).
520 The intended usage is
521 (1) setup struct with unchanging data with calls to addto_INCORR_BLOK_set
522 (2) finalize the correlation with calls to correlate_INCORR_BLOK_set
523 with data that changes -- from the warp patch currently being
524 optimized, that is
525 *//*-------------------------------------------------------------------------*/
526
correlate_INCORR_BLOK_set(INCORR_BLOK_set * ibs,int ibot,int itop,int jbot,int jtop,int kbot,int ktop,float * xval,float * yval,float * wt)527 float correlate_INCORR_BLOK_set( INCORR_BLOK_set *ibs ,
528 int ibot , int itop ,
529 int jbot , int jtop ,
530 int kbot , int ktop ,
531 float *xval , float *yval , float *wt )
532 {
533 int dd , nblok=ibs->numblok ;
534 int *nelm=ibs->nelm , *nsum=NULL ;
535 float *sx , *sxx , *sy , *syy , *sxy , *sw ;
536 double xx,yy,xy,ww ; float rval , wsum,psum ;
537
538 if( xval == NULL || yval == NULL ){ /* no extra data to add in */
539
540 sx = ibs->sx ; /* pointers to external struct data */
541 sxx = ibs->sxx ; /* this data will not be altered herein */
542 sy = ibs->sy ;
543 syy = ibs->syy ;
544 sxy = ibs->sxy ;
545 sw = ibs->sw ;
546
547 } else { /* add the extra data into a local copy of the struct's data */
548
549 int *blkn = ibs->blkn ; /* spatial index to blok index array */
550 int nx,ny,nz,nxy , ii,jj,kk , pp ;
551
552 nx = ibs->nx ; ny = ibs->ny ; nz = ibs->nz ; nxy = nx*ny ;
553
554 sx = (float *)calloc(sizeof(float),nblok) ; /* local copies */
555 sxx = (float *)calloc(sizeof(float),nblok) ; /* to be thrown */
556 sy = (float *)calloc(sizeof(float),nblok) ; /* away at end */
557 syy = (float *)calloc(sizeof(float),nblok) ;
558 sxy = (float *)calloc(sizeof(float),nblok) ;
559 sw = (float *)calloc(sizeof(float),nblok) ;
560 nsum = (int * )calloc(sizeof(int) ,nblok) ;
561 if( ibs->ninserted > 0 ){
562 AAmemcpy(sxx ,ibs->sxx ,sizeof(float)*nblok) ;
563 AAmemcpy(sy ,ibs->sy ,sizeof(float)*nblok) ;
564 AAmemcpy(syy ,ibs->syy ,sizeof(float)*nblok) ;
565 AAmemcpy(sxy ,ibs->sxy ,sizeof(float)*nblok) ;
566 AAmemcpy(sw ,ibs->sw ,sizeof(float)*nblok) ;
567 AAmemcpy(nsum,ibs->nsum,sizeof(int) *nblok) ;
568 }
569
570 if( debug_lpc )
571 fprintf(stderr,"++++++++++ Debug LPC output: %d bloks inserting x=%d..%d y=%d..%d z=%d..%d wt=%p\n",
572 nblok,ibot,itop,jbot,jtop,kbot,ktop,(void *)wt) ;
573
574 if( wt == NULL ){
575 for( pp=0,kk=kbot ; kk <= ktop ; kk++ ){
576 for( jj=jbot ; jj <= jtop ; jj++ ){
577 for( ii=ibot ; ii <= itop ; ii++,pp++ ){
578 dd = blkn[ii + jj*nx + kk*nxy] ;
579 if( dd >= 0 && nsum[dd] < nelm[dd] ){ /* in incomplete blok? */
580 xx = (double)xval[pp] ; yy = (double)yval[pp] ;
581 sx[dd] += xx ; sxx[dd] += xx*xx ;
582 sy[dd] += yy ; syy[dd] += yy*yy ; sxy[dd] += xx*yy ; sw[dd] += 1.0 ;
583 nsum[dd]++ ;
584 }
585 }}}
586 } else {
587 for( pp=0,kk=kbot ; kk <= ktop ; kk++ ){
588 for( jj=jbot ; jj <= jtop ; jj++ ){
589 for( ii=ibot ; ii <= itop ; ii++,pp++ ){
590 dd = blkn[ii + jj*nx + kk*nxy] ;
591 if( dd >= 0 && nsum[dd] < nelm[dd] ){ /* in incomplete blok? */
592 ww = wt[pp] ; if( ww <= 0.0 ) continue ; /* no weight? */
593 xx = (double)xval[pp] ; yy = (double)yval[pp] ;
594 sx[dd] += xx*ww; sxx[dd] += xx*xx*ww;
595 sy[dd] += yy*ww; syy[dd] += yy*yy*ww; sxy[dd] += xx*yy*ww; sw[dd] += ww;
596 nsum[dd]++ ;
597 }
598 }}}
599 }
600
601 if( debug_lpc )
602 fprintf(stderr," inserted %d points in temp copies\n",pp) ;
603
604 } /* end of putting data into local copy */
605
606 /*-- now finalize the summation --*/
607
608 psum = ibs->psum ; wsum = ibs->wsum ;
609
610 if( debug_lpc )
611 fprintf(stderr,"++++++++++ Debug LPC output: %d bloks init psum=%g wsum=%g\n pvals:",nblok,psum,wsum) ;
612
613 for( dd=0 ; dd < nblok ; dd++ ){ /* scan and finish incomplete bloks */
614 if( sw[dd] > 0.0f ){ /* incomplete but with data? */
615 ww = 1.0 / sw[dd] ;
616 xx = sxx[dd] - sx[dd] * sx[dd] * ww ;
617 yy = syy[dd] - sy[dd] * sy[dd] * ww ;
618 xy = sxy[dd] - sx[dd] * sy[dd] * ww ;
619 if( xx > 0.0 && yy > 0.0 ){
620 rval = (float)(xy/sqrt(xx*yy)) ;
621 if( rval > CMAX ) rval = CMAX ;
622 else if( rval < -CMAX ) rval = -CMAX ;
623 rval = logf( (1.0f+rval)/(1.0f-rval) ) ;
624 psum += (float)sw[dd] * rval * fabsf(rval) ;
625 if( debug_lpc ) fprintf(stderr," %g",rval) ;
626 } else if( debug_lpc ){
627 fprintf(stderr," {sx=%g sxx=%g sy=%g syy=%g sxy=%g sw=%g xx=%g yy=%g xy=%g}",
628 sx[dd],sxx[dd],sy[dd],syy[dd],sxy[dd],sw[dd],xx,yy,xy) ;
629 }
630 wsum += (float)sw[dd] ;
631 } else if( debug_lpc ){
632 fprintf(stderr," %g",ibs->pval[dd]) ;
633 }
634 }
635
636 if( sx != ibs->sx ){ /* free local copies */
637 free(sx); free(sxx); free(sy); free(syy); free(sxy); free(sw); free(nsum);
638 }
639
640 rval = (wsum <= 0.0f) ? 0.0f : psum/wsum ;
641
642 if( debug_lpc )
643 fprintf(stderr,"\n--------- final psum=%g wsum=%g rval=%g ----------------------------\n",psum,wsum,rval) ;
644
645 return rval ;
646 }
647
648 /****************************************************************************/
649 /*** Histogram-based measurements of dependence between two float arrays. ***/
650 /****************************************************************************/
651
652 #undef WW
653 #define WW(i) ((w==NULL) ? 1.0f : w[i]) /* weight function for i'th datum */
654
655 #undef XYC
656 #define XYC(p,q) xyc[(p)+(q)*nbp]
657
658 #ifndef WAY_BIG
659 # define WAY_BIG 1.e+10
660 #endif
661
662 #undef GOODVAL
663 #define GOODVAL(x) ((x) < WAY_BIG) /* x is not preposterous */
664
665 #undef RANGVAL
666 #define RANGVAL(x,b,t) ((x) >= (b) && (x) <= (t)) /* x between b and t */
667
668 /*--------------------------------------------------------------------------*/
669
INCOR_clipate(int nval,float * xar)670 float_pair INCOR_clipate( int nval , float *xar )
671 {
672 MRI_IMAGE *qim; float cbot,ctop, mmm , *qar; float_pair rr; int ii,nq;
673
674 ENTRY("INCOR_clipate") ;
675
676 qim = mri_new_vol( nval,1,1 , MRI_float ) ; qar = MRI_FLOAT_PTR(qim) ;
677 for( ii=nq=0 ; ii < nval ; ii++ ) if( GOODVAL(xar[ii]) ) qar[nq++] = xar[ii];
678 qim->nx = qim->nvox = nq ;
679 if( nq < 666 ){ rr.a = 1.0f; rr.b = 0.0f; mri_free(qim); RETURN(rr); }
680 mmm = mri_min( qim ) ;
681 if( mmm >= 0.0f ){ /* for positive images */
682 cbot = THD_cliplevel( qim , 0.321f ) ;
683 ctop = mri_quantile ( qim , 0.987f ) ;
684 if( ctop > 6.543f*cbot ) ctop = 6.543f*cbot ;
685 } else { /* for images including negative values: no go */
686 cbot = 1.0f; ctop = 0.0f;
687 }
688 mri_free(qim) ;
689 rr.a = cbot ; rr.b = ctop ; RETURN(rr) ;
690 }
691
692 /*--------------------------------------------------------------------------*/
693
INCOR_2Dhist_xyclip(int nval,float * xval,float * yval)694 float_quad INCOR_2Dhist_xyclip( int nval , float *xval , float *yval )
695 {
696 float_pair xcc , ycc ; float_quad xxyycc={0.0f,0.0f,0.0f,0.0f} ;
697
698 ENTRY("INCOR_2Dhist_xyclip") ;
699
700 if( nval < 666 || xval == NULL || yval == NULL ) RETURN(xxyycc) ;
701
702 xcc = INCOR_clipate( nval , xval ) ;
703 ycc = INCOR_clipate( nval , yval ) ;
704
705 if( xcc.a >= xcc.b || ycc.a >= ycc.b ) RETURN(xxyycc) ;
706
707 xxyycc.a = xcc.a ; xxyycc.b = xcc.b ;
708 xxyycc.c = ycc.a ; xxyycc.d = ycc.b ; RETURN(xxyycc) ;
709 }
710
711 /*--------------------------------------------------------------------------*/
712
INCOR_2Dhist_minmax(int nval,float * xval,float * yval)713 float_quad INCOR_2Dhist_minmax( int nval , float *xval , float *yval )
714 {
715 float_quad xxyy={0.0f,0.0f,0.0f,0.0f} ;
716 int ii ; float xb,xt,yb,yt ;
717
718 ENTRY("INCOR_2Dhist_minmax") ;
719
720 if( nval < 1 || xval == NULL || yval == NULL ) RETURN(xxyy) ;
721
722 xb = xt = xval[0] ; yb = yt = yval[0] ;
723 for( ii=1 ; ii < nval ; ii++ ){
724 if( xval[ii] < xb ) xb = xval[ii] ;
725 else if( xval[ii] > xt ) xt = xval[ii] ;
726 if( yval[ii] < yb ) yb = yval[ii] ;
727 else if( yval[ii] > yt ) yt = yval[ii] ;
728 }
729 xxyy.a = xb ; xxyy.b = xt ; xxyy.c = yb ; xxyy.d = yt ; RETURN(xxyy) ;
730 }
731
732 /*--------------------------------------------------------------------------*/
733
734 static byte *good=NULL ;
735 static int agood=0 ;
736
INCOR_setup_good(int ng)737 void INCOR_setup_good( int ng )
738 {
739 if( ng <= 0 ){
740 if( good != NULL ){ free(good); good = NULL; }
741 agood = 0 ;
742 } else if( ng > agood ){
743 good = realloc( good , sizeof(byte)*ng ) ; agood = ng ;
744 }
745 if( agood > 0 && good != NULL ) AAmemset(good,0,sizeof(byte)*agood) ;
746 return ;
747 }
748
749 /*--------------------------------------------------------------------------*/
750 /*! Load 2D histogram of x[0..n-1] and y[0..n-1], each point optionally
751 weighted by w[0..n-1] (weights are all 1 if w==NULL).
752 Used in the histogram-based measures of dependence between x[] and y[i].
753 If something is bad on input, nbin is set to 0. Otherwise, these global
754 variables are set:
755 - nbin = # of bins, nbp = nbin+1
756 - nww = sum of the weights used
757 - xc = marginal histogram of x[], for xc[0..nbin] (nbp points in)
758 - yc = marginal histogram of y[], for yc[0..nbin] (each direction)
759 - xyc = joint histogram of (x[],y[]), for XYC(0..nbin,0..nbin)
760 - The histograms can be later normalized (by 1/nww) to have sum==1
761 - Histogram can be retrieved by retrieve_2Dhist() and can be
762 erased by clear_2Dhist()
763 - Default number of equal-spaced bins in each direction is n^(1/3)
764 - the exponent can be changed with INCOR_set_2Dhist_hpower()
765 - you can set the number of bins with INCOR_set_2Dhist_hbin()
766 - x[] values outside the range xbot..xtop (inclusive) will not be
767 used in the histogram; mutatis mutandum for y[]
768 *//*------------------------------------------------------------------------*/
769
INCOR_addto_2Dhist(INCOR_2Dhist * tdh,int n,float * x,float * y,float * w)770 void INCOR_addto_2Dhist( INCOR_2Dhist *tdh , int n , float *x , float *y , float *w )
771 {
772 int ii , ngood , xyclip ;
773 float xxbot,xxtop , yybot,yytop ;
774 float xcbot,xctop , ycbot,yctop ;
775 int nbin,nbp,nbm ;
776 float *xc , *yc , *xyc ; float nww ;
777 #if 0
778 int use_omp , nthr ;
779 #else
780 # define use_omp 0
781 # define nthr 1
782 #endif
783
784 ENTRY("INCOR_addto_2Dhist") ;
785
786 if( tdh == NULL || tdh->nbin < 3 || n <= 0 || x == NULL || y == NULL ) EXRETURN ;
787
788 nbin = tdh->nbin ; nbp = nbin+1 ; nbm = nbin-1 ;
789 xc = tdh->xc ; yc = tdh->yc ; xyc = tdh->xyc ; nww = tdh->nww ;
790
791 /* get the min..max range for x and y data? */
792
793 INCOR_setup_good(n+4) ;
794 for( ngood=ii=0 ; ii < n ; ii++ ){
795 if( GOODVAL(x[ii]) && GOODVAL(y[ii]) && (WW(ii) > 0.0f) ){
796 good[ii] = 1 ; ngood++ ;
797 }
798 }
799 if( ngood == 0 ) EXRETURN ;
800
801 xxbot = tdh->xxbot ; xxtop = tdh->xxtop ;
802 yybot = tdh->yybot ; yytop = tdh->yytop ;
803 xcbot = tdh->xcbot ; xctop = tdh->xctop ;
804 ycbot = tdh->ycbot ; yctop = tdh->yctop ;
805
806 if( (xxbot >= xxtop) || (yybot >= yytop) ){ /* data ranges undefined */
807
808 xxbot = WAY_BIG ; xxtop = -WAY_BIG ;
809 for( ii=0 ; ii < n ; ii++ ){
810 if( good[ii] ){
811 if( x[ii] > xxtop ) xxtop = x[ii] ;
812 else if( x[ii] < xxbot ) xxbot = x[ii] ;
813 }
814 }
815 if( xxbot >= xxtop ) EXRETURN ;
816
817 yybot = WAY_BIG ; yytop = -WAY_BIG ;
818 for( ii=0 ; ii < n ; ii++ ){
819 if( good[ii] ){
820 if( y[ii] > yytop ) yytop = y[ii] ;
821 else if( y[ii] < yybot ) yybot = y[ii] ;
822 }
823 }
824 if( yybot >= yytop ) EXRETURN ;
825
826 tdh->xxbot = xxbot ; tdh->xxtop = xxtop ;
827 tdh->yybot = yybot ; tdh->yytop = yytop ;
828
829 xcbot = ycbot = tdh->xcbot = tdh->ycbot = WAY_BIG ; /* disable */
830 xctop = yctop = tdh->xctop = tdh->yctop = -WAY_BIG ; /* clipping */
831 }
832
833 /*-- count number of good values left in range (in both x and y) --*/
834
835 AAmemset(good,0,n) ;
836
837 for( ngood=ii=0 ; ii < n ; ii++ ){
838 if( RANGVAL(x[ii],xxbot,xxtop) && RANGVAL(y[ii],yybot,yytop) && (WW(ii) > 0.0f) ){
839 good[ii] = 1 ; ngood++ ;
840 }
841 }
842 if( ngood == 0 ) EXRETURN ;
843
844 /*--------------- add to the 2D and 1D histograms ---------------*/
845
846 xyclip = (xxbot < xcbot) && (xcbot < xctop) && (xctop < xxtop) &&
847 (yybot < ycbot) && (ycbot < yctop) && (yctop < yytop) ;
848
849 #ifndef nthr
850 nthr = omp_get_max_threads() ;
851 use_omp = (nthr > 1) ;
852 #endif
853
854 if( !use_omp ){ /*** serial code ***/
855
856 /* AFNI_do_nothing() ; fprintf(stderr,"h") ; */
857
858 if( !xyclip ){ /*------------ equal size bins ------------*/
859
860 float xb,xi , yb,yi , xx,yy , x1,y1 , ww ;
861 int jj,kk ;
862
863 xb = xxbot ; xi = nbm/(xxtop-xxbot) ;
864 yb = yybot ; yi = nbm/(yytop-yybot) ;
865 for( ii=0 ; ii < n ; ii++ ){
866 if( !good[ii] ) continue ;
867 xx = (x[ii]-xb)*xi ;
868 jj = (int)xx ; xx = xx - jj ; x1 = 1.0f-xx ;
869 yy = (y[ii]-yb)*yi ;
870 kk = (int)yy ; yy = yy - kk ; y1 = 1.0f-yy ;
871 ww = WW(ii) ; nww += ww ;
872
873 if( jj < 0 || kk < 0 || jj >= nbin || kk >= nbin ) continue ;
874
875 xc[jj] += (x1*ww); xc[jj+1] += (xx*ww);
876 yc[kk] += (y1*ww); yc[kk+1] += (yy*ww);
877
878 XYC(jj ,kk ) += x1*(y1*ww) ;
879 XYC(jj+1,kk ) += xx*(y1*ww) ;
880 XYC(jj ,kk+1) += x1*(yy*ww) ;
881 XYC(jj+1,kk+1) += xx*(yy*ww) ;
882 }
883
884 } else if( xyclip ){ /*------------ mostly equal bins ----------------*/
885
886 int jj,kk ;
887 float xbc=xcbot , xtc=xctop , ybc=ycbot , ytc=yctop ;
888 float xi,yi , xx,yy , x1,y1 , ww ;
889
890 AFNI_do_nothing() ; /* fprintf(stderr,"c") ; */
891
892 xi = (nbin-2.000001f)/(xtc-xbc) ;
893 yi = (nbin-2.000001f)/(ytc-ybc) ;
894 for( ii=0 ; ii < n ; ii++ ){
895 if( !good[ii] ) continue ;
896 xx = x[ii] ;
897 if( xx < xbc ){ jj = 0 ; xx = 0.0f ; }
898 else if( xx > xtc ){ jj = nbm ; xx = 1.0f ; }
899 else { xx = 1.0f+(xx-xbc)*xi; jj = (int)xx; xx = xx - jj; }
900 yy = y[ii] ;
901 if( yy < ybc ){ kk = 0 ; yy = 0.0f ; }
902 else if( yy > ytc ){ kk = nbm ; yy = 1.0f ; }
903 else { yy = 1.0f+(yy-ybc)*yi; kk = (int)yy; yy = yy - kk; }
904
905 if( jj < 0 || kk < 0 || jj >= nbin || kk >= nbin ) continue ;
906
907 x1 = 1.0f-xx ; y1 = 1.0f-yy ; ww = WW(ii) ; nww += ww ;
908
909 xc[jj] += (x1*ww); xc[jj+1] += (xx*ww);
910 yc[kk] += (y1*ww); yc[kk+1] += (yy*ww);
911
912 XYC(jj ,kk ) += x1*(y1*ww) ;
913 XYC(jj+1,kk ) += xx*(y1*ww) ;
914 XYC(jj ,kk+1) += x1*(yy*ww) ;
915 XYC(jj+1,kk+1) += xx*(yy*ww) ;
916 }
917
918 AFNI_do_nothing() ; /* fprintf(stderr,".") ; */
919
920 } /* end of clipped code */
921
922 } else { /*** parallelized using OpenMP ***/
923
924 float **xccar , **yccar , **xyccar , *nwwar ; int nbpq=nbp*nbp,itt ;
925
926 xccar = (float **)calloc(sizeof(float *),nthr) ; /* arrays for */
927 yccar = (float **)calloc(sizeof(float *),nthr) ; /* accumulation */
928 xyccar = (float **)calloc(sizeof(float *),nthr) ; /* in separate */
929 nwwar = (float * )calloc(sizeof(float) ,nthr) ; /* threads */
930
931 for( itt=0 ; itt < nthr ; itt++ ){
932 xccar [itt] = (float *)calloc(sizeof(float),nbp ) ;
933 yccar [itt] = (float *)calloc(sizeof(float),nbp ) ;
934 xyccar[itt] = (float *)calloc(sizeof(float),nbpq) ;
935 }
936
937 AFNI_do_nothing() ; fprintf(stderr,"H") ;
938
939 #undef XYCC
940 #define XYCC(p,q) xycc[(p)+(q)*nbp]
941
942 if( !xyclip ){ /*------------ equal size bins ------------*/
943
944 AFNI_do_nothing() ; fprintf(stderr,"y") ;
945
946 AFNI_OMP_START ;
947 #pragma omp parallel /*** start of parallel code ***/
948 {
949 float *xcc, *ycc , *xycc ;
950 float xb,xi , yb,yi , xx,yy , x1,y1 , ww ;
951 int ii,jj,kk , ithr ;
952
953 ithr = omp_get_thread_num() ;
954 #pragma omp barrier
955 fprintf(stderr,"%d",ithr) ;
956 xcc = xccar[ithr] ; ycc = yccar[ithr] ; xycc = xyccar[ithr] ;
957 xb = xxbot ; xi = nbm/(xxtop-xxbot) ;
958 yb = yybot ; yi = nbm/(yytop-yybot) ;
959 #pragma omp for
960 for( ii=0 ; ii < n ; ii++ ){
961 if( !good[ii] ) continue ;
962 xx = (x[ii]-xb)*xi ;
963 jj = (int)xx ; xx = xx-jj ; x1 = 1.0f-xx ;
964 yy = (y[ii]-yb)*yi ;
965 kk = (int)yy ; yy = yy-kk ; y1 = 1.0f-yy ;
966 ww = WW(ii) ; nwwar[ithr] += ww ;
967
968 if( jj < 0 || kk < 0 || jj >= nbin || kk >= nbin ) continue ;
969
970 xcc[jj] += (x1*ww); xcc[jj+1] += (xx*ww);
971 ycc[kk] += (y1*ww); ycc[kk+1] += (yy*ww);
972
973 XYCC(jj ,kk ) += x1*(y1*ww) ;
974 XYCC(jj+1,kk ) += xx*(y1*ww) ;
975 XYCC(jj ,kk+1) += x1*(yy*ww) ;
976 XYCC(jj+1,kk+1) += xx*(yy*ww) ;
977 }
978 } /*** end of parallel code ***/
979 AFNI_OMP_END ;
980
981 } else if( xyclip ){ /*------------ mostly equal bins ----------------*/
982
983 AFNI_do_nothing() ; fprintf(stderr,"x") ;
984
985 AFNI_OMP_START ;
986 #pragma omp parallel /*** start of parallel code ***/
987 {
988 float *xcc, *ycc , *xycc ;
989 int ii,jj,kk , ithr ;
990 float xbc=xcbot , xtc=xctop , ybc=ycbot , ytc=yctop ;
991 float xi,yi , xx,yy , x1,y1 , ww ;
992
993 ithr = omp_get_thread_num() ;
994 fprintf(stderr,"%d",ithr) ;
995 #pragma omp barrier
996 fprintf(stderr,":") ;
997 xcc = xccar[ithr] ; ycc = yccar[ithr] ; xycc = xyccar[ithr] ;
998 xi = (nbin-2.000001f)/(xtc-xbc) ;
999 yi = (nbin-2.000001f)/(ytc-ybc) ;
1000 #pragma omp for
1001 for( ii=0 ; ii < n ; ii++ ){
1002 if( !good[ii] ) continue ;
1003 xx = x[ii] ;
1004 if( xx < xbc ){ jj = 0 ; xx = 0.0f ; }
1005 else if( xx > xtc ){ jj = nbm ; xx = 1.0f ; }
1006 else { xx = 1.0f+(xx-xbc)*xi; jj = (int)xx; xx = xx-jj; }
1007 yy = y[ii] ;
1008 if( yy < ybc ){ kk = 0 ; yy = 0.0f ; }
1009 else if( yy > ytc ){ kk = nbm ; yy = 1.0f ; }
1010 else { yy = 1.0f+(yy-ybc)*yi; kk = (int)yy; yy = yy-kk; }
1011
1012 if( jj < 0 || kk < 0 || jj >= nbin || kk >= nbin ) continue ;
1013
1014 x1 = 1.0f-xx ; y1 = 1.0f-yy ; ww = WW(ii) ; nwwar[ithr] += ww ;
1015
1016 xcc[jj] += (x1*ww); xcc[jj+1] += (xx*ww);
1017 ycc[kk] += (y1*ww); ycc[kk+1] += (yy*ww);
1018
1019 XYCC(jj ,kk ) += x1*(y1*ww) ;
1020 XYCC(jj+1,kk ) += xx*(y1*ww) ;
1021 XYCC(jj ,kk+1) += x1*(yy*ww) ;
1022 XYCC(jj+1,kk+1) += xx*(yy*ww) ;
1023 }
1024 #pragma omp barrier
1025 fprintf(stderr,"%d",ithr) ;
1026 } /*** end of parallel code ***/
1027 AFNI_OMP_END ;
1028
1029 } /*-- end of mostly equal bins --*/
1030
1031 /* now merge the parallel thread results */
1032
1033 for( itt=0 ; itt < nthr ; itt++ ){
1034 if( nwwar[itt] > 0.0f ){
1035 nww += nwwar[itt] ;
1036 for( ii=0 ; ii < nbp ; ii++ ){ xc[ii] += xccar[itt][ii]; yc[ii] += yccar[itt][ii]; }
1037 for( ii=0 ; ii < nbpq ; ii++ ){ xyc[ii] += xyccar[itt][ii] ; }
1038 }
1039 free(xccar[itt]) ; free(yccar[itt]) ; free(xyccar[itt]) ;
1040 }
1041 free(xccar) ; free(yccar) ; free(xyccar) ; free(nwwar) ;
1042
1043 } /* end of using OpenMP */
1044
1045 /* AFNI_do_nothing() ; fprintf(stderr,".") ; */
1046
1047 tdh->nww = nww ;
1048 EXRETURN ;
1049 }
1050
1051 /*----------------------------------------------------------------------------*/
1052 /* scale histogram to have sum==1 */
1053
INCOR_normalize_2Dhist(INCOR_2Dhist * tdh)1054 void INCOR_normalize_2Dhist( INCOR_2Dhist *tdh )
1055 {
1056 float nww , *xc, *yc, *xyc ; int nbp ;
1057 if( tdh == NULL ) return ;
1058 nww = tdh->nww; xc = tdh->xc; yc = tdh->yc; xyc = tdh->xyc; nbp = tdh->nbin+1;
1059 if( nww > 0.0f && nww != 1.0f && xyc != NULL && xc != NULL && yc != NULL ){
1060 float ni ; int nbq , ii ;
1061 ni = 1.0f / nww ;
1062 for( ii=0 ; ii < nbp ; ii++ ){ xc[ii] *= ni; yc[ii] *= ni; }
1063 nbq = nbp*nbp ;
1064 for( ii=0 ; ii < nbq ; ii++ ){ xyc[ii] *= ni; }
1065 }
1066 return ;
1067 }
1068
1069 /*--------------------------------------------------------------------------*/
1070 /*! Compute the mutual info from a histogram (which also normalizes it).
1071 ----------------------------------------------------------------------------*/
1072
INCOR_mutual_info(INCOR_2Dhist * tdh)1073 double INCOR_mutual_info( INCOR_2Dhist *tdh )
1074 {
1075 int ii,jj ;
1076 double val , nww ;
1077 float *xc, *yc, *xyc ; int nbp ;
1078
1079 if( tdh == NULL ) return 0.0 ;
1080
1081 nww = (double)tdh->nww; if( nww <= 0.0 ) return 0.0 ;
1082 xc = tdh->xc; yc = tdh->yc; xyc = tdh->xyc; nbp = tdh->nbin+1;
1083
1084 INCOR_normalize_2Dhist(tdh) ;
1085
1086 /*-- compute MI from histogram --*/
1087
1088 val = 0.0 ;
1089 for( ii=0 ; ii < nbp ; ii++ ){
1090 for( jj=0 ; jj < nbp ; jj++ ){
1091 if( XYC(ii,jj) > 0.0f )
1092 val += XYC(ii,jj) * log( XYC(ii,jj)/(xc[ii]*yc[jj]) ) ;
1093 }}
1094 return (1.4427*val) ; /* units are bits, just for fun */
1095 }
1096
1097 /*--------------------------------------------------------------------------*/
1098 /*! Compute the normalized mutual info from a 2D histogram.
1099 Actually, returns H(x,y) / [ H(x)+H(y) ], which should be small if
1100 x and y are redundant and should be large if they are independent.
1101 ----------------------------------------------------------------------------*/
1102
INCOR_norm_mutinf(INCOR_2Dhist * tdh)1103 double INCOR_norm_mutinf( INCOR_2Dhist *tdh )
1104 {
1105 int ii,jj ;
1106 double numer , denom , nww ;
1107 float *xc, *yc, *xyc ; int nbp ;
1108
1109 if( tdh == NULL ) return 0.0 ;
1110
1111 nww = (double)tdh->nww; if( nww <= 0.0 ) return 0.0 ;
1112 xc = tdh->xc; yc = tdh->yc; xyc = tdh->xyc; nbp = tdh->nbin+1;
1113
1114 INCOR_normalize_2Dhist(tdh) ;
1115
1116 /*-- compute NMI from histogram --*/
1117
1118 denom = numer = 0.0 ;
1119 for( ii=0 ; ii < nbp ; ii++ ){
1120 if( xc[ii] > 0.0f ) denom += xc[ii] * log( xc[ii] ) ;
1121 if( yc[ii] > 0.0f ) denom += yc[ii] * log( yc[ii] ) ;
1122 for( jj=0 ; jj < nbp ; jj++ ){
1123 if( XYC(ii,jj) > 0.0f ) numer += XYC(ii,jj) * log( XYC(ii,jj) );
1124 }
1125 }
1126 if( denom != 0.0 ) denom = numer / denom ;
1127 return denom ;
1128 }
1129
1130 /*--------------------------------------------------------------------------*/
1131 /*! Compute the correlation ratio from a 2D histogram.
1132 ----------------------------------------------------------------------------*/
1133
INCOR_corr_ratio(INCOR_2Dhist * tdh,int crmode)1134 double INCOR_corr_ratio( INCOR_2Dhist *tdh , int crmode )
1135 {
1136 int ii,jj ;
1137 double vv,mm , val , cyvar , uyvar , yrat,xrat , nww ;
1138 float *xc, *yc, *xyc ; int nbp ;
1139
1140 if( tdh == NULL ) return 0.0 ;
1141
1142 nww = (double)tdh->nww; if( nww <= 0.0 ) return 0.0 ;
1143 xc = tdh->xc; yc = tdh->yc; xyc = tdh->xyc; nbp = tdh->nbin+1;
1144
1145 INCOR_normalize_2Dhist(tdh) ;
1146
1147 /*-- compute CR(y|x) from histogram --*/
1148
1149 cyvar = 0.0 ;
1150 for( ii=0 ; ii < nbp ; ii++ ){
1151 if( xc[ii] > 0.0f ){
1152 vv = mm = 0.0 ; /* mm=E(y|x) vv=E(y^2|x) */
1153 for( jj=1 ; jj < nbp ; jj++ ){
1154 mm += (jj * XYC(ii,jj)) ; vv += jj * (jj * XYC(ii,jj)) ;
1155 }
1156 cyvar += (vv - mm*mm/xc[ii] ) ; /* Var(y|x) */
1157 }
1158 }
1159 vv = mm = uyvar = 0.0 ;
1160 for( jj=1 ; jj < nbp ; jj++ ){ /* mm=E(y) vv=E(y^2) */
1161 mm += (jj * yc[jj]) ; vv += jj * (jj * yc[jj]) ;
1162 }
1163 uyvar = vv - mm*mm ; /* Var(y) */
1164 yrat = (uyvar > 0.0 ) ? cyvar/uyvar /* Var(y|x) / Var(y) */
1165 : 1.0 ;
1166
1167 if( crmode == 0 ) return (1.0-yrat) ; /** unsymmetric **/
1168
1169 /** compute CR(x|y) also, for symmetrization **/
1170
1171 cyvar = 0.0 ;
1172 for( jj=0 ; jj < nbp ; jj++ ){
1173 if( yc[jj] > 0.0 ){
1174 vv = mm = 0.0 ; /* mm=E(x|y) vv=E(x^2|y) */
1175 for( ii=1 ; ii < nbp ; ii++ ){
1176 mm += (ii * XYC(ii,jj)) ; vv += ii * (ii * XYC(ii,jj)) ;
1177 }
1178 cyvar += (vv - mm*mm/yc[jj] ) ; /* Var(x|y) */
1179 }
1180 }
1181 vv = mm = uyvar = 0.0 ;
1182 for( ii=1 ; ii < nbp ; ii++ ){ /* mm=E(x) vv=E(x^2) */
1183 mm += (ii * xc[ii]) ; vv += ii * (ii * xc[ii]) ;
1184 }
1185 uyvar = vv - mm*mm ; /* Var(x) */
1186 xrat = (uyvar > 0.0) ? cyvar/uyvar /* Var(x|y) / Var(x) */
1187 : 1.0 ;
1188
1189 if( crmode == 2 ) return (1.0 - 0.5*(xrat+yrat)) ; /** additive **/
1190
1191 return (1.0 - xrat*yrat) ; /** multiplicative **/
1192 }
1193
1194 /*--------------------------------------------------------------------------*/
1195 /*! Compute the Hellinger metric from a 2D histogram.
1196 ----------------------------------------------------------------------------*/
1197
INCOR_hellinger(INCOR_2Dhist * tdh)1198 double INCOR_hellinger( INCOR_2Dhist *tdh )
1199 {
1200 int ii,jj ;
1201 double val , pq , nww ;
1202 float *xc, *yc, *xyc ; int nbp ;
1203
1204 if( tdh == NULL ) return 0.0 ;
1205
1206 nww = (double)tdh->nww; if( nww <= 0.0 ) return 0.0 ;
1207 xc = tdh->xc; yc = tdh->yc; xyc = tdh->xyc; nbp = tdh->nbin+1;
1208
1209 INCOR_normalize_2Dhist(tdh) ;
1210
1211 /*-- compute Hell metric from histogram --*/
1212
1213 val = 0.0 ;
1214 for( ii=0 ; ii < nbp ; ii++ ){
1215 for( jj=0 ; jj < nbp ; jj++ ){
1216 pq = XYC(ii,jj) ;
1217 if( pq > 0.0 ) val += sqrt( pq * xc[ii] * yc[jj] ) ;
1218 }}
1219 return (1.0-val) ;
1220 }
1221
1222 /*--------------------------------------------------------------------------*/
1223 /*! Compute the Hellinger metric, mutual info, normalized MI, and
1224 (additively) symmetrized correlation ratio, and return all 4
1225 (in that order), from a 2D histogram.
1226
1227 The first 3 values all measure the closeness of the joint histogram to
1228 the product of the marginals:
1229 - Hellinger is smaller when the joint is closer to the marginals' product
1230 - MI is also smaller when the joint is closer to the marginal's product
1231 - NMI is larger when the joint is closer to the marginal's product
1232 Correlation ratio (symmetrized by addition == CRA) is larger when
1233 the two variables are nonlinearly correlated.
1234
1235 As measures of association (generalized correlation): more closely
1236 associated variables correspond to larger Hellinger and MI and CRA,
1237 and to smaller NMI.
1238 *//*------------------------------------------------------------------------*/
1239
INCOR_helmicra(INCOR_2Dhist * tdh)1240 double_quad INCOR_helmicra( INCOR_2Dhist *tdh )
1241 {
1242 int ii,jj ;
1243 double hel , pq , vv,uu , nww ;
1244 double val , cyvar , uyvar , yrat,xrat ;
1245 double_quad hmc = {0.0,0.0,0.0,0.0} ;
1246 float *xc, *yc, *xyc ; int nbp ;
1247
1248 if( tdh == NULL ) return hmc ;
1249
1250 nww = (double)tdh->nww; if( nww <= 0.0 ) return hmc ;
1251 xc = tdh->xc; yc = tdh->yc; xyc = tdh->xyc; nbp = tdh->nbin+1;
1252
1253 INCOR_normalize_2Dhist(tdh) ;
1254
1255 /*-- compute Hel, MI, NMI from histogram --*/
1256
1257 hel = vv = uu = 0.0 ;
1258 for( ii=0 ; ii < nbp ; ii++ ){
1259 if( xc[ii] > 0.0 ) vv += xc[ii] * log( xc[ii] ) ;
1260 if( yc[ii] > 0.0 ) vv += yc[ii] * log( yc[ii] ) ;
1261 for( jj=0 ; jj < nbp ; jj++ ){
1262 pq = XYC(ii,jj) ;
1263 if( pq > 0.0 ){
1264 hel += sqrt( pq * xc[ii] * yc[jj] ) ;
1265 uu += pq * log( pq );
1266 }
1267 }
1268 }
1269 hmc.a = 1.0 - hel ; /* Hellinger */
1270 hmc.b = uu - vv ; /* MI */
1271 hmc.c = (vv != 0.0) ? uu/vv : 0.0 ; /* NMI */
1272
1273 /*-- compute CR(y|x) from histogram --*/
1274
1275 cyvar = 0.0 ;
1276 for( ii=0 ; ii < nbp ; ii++ ){
1277 if( xc[ii] > 0.0 ){
1278 vv = uu = 0.0 ; /* uu=E(y|x) vv=E(y^2|x) */
1279 for( jj=1 ; jj < nbp ; jj++ ){
1280 uu += (jj * XYC(ii,jj)) ; vv += jj * (jj * XYC(ii,jj)) ;
1281 }
1282 cyvar += (vv - uu*uu/xc[ii] ) ; /* Var(y|x) */
1283 }
1284 }
1285 vv = uu = uyvar = 0.0 ;
1286 for( jj=1 ; jj < nbp ; jj++ ){ /* uu=E(y) vv=E(y^2) */
1287 uu += (jj * yc[jj]) ; vv += jj * (jj * yc[jj]) ;
1288 }
1289 uyvar = vv - uu*uu ; /* Var(y) */
1290 yrat = (uyvar > 0.0) ? cyvar/uyvar /* Var(y|x) / Var(y) */
1291 : 1.0 ;
1292
1293 /** compute CR(x|y) also, for symmetrization **/
1294
1295 cyvar = 0.0 ;
1296 for( jj=0 ; jj < nbp ; jj++ ){
1297 if( yc[jj] > 0.0 ){
1298 vv = uu = 0.0 ; /* uu=E(x|y) vv=E(x^2|y) */
1299 for( ii=1 ; ii < nbp ; ii++ ){
1300 uu += (ii * XYC(ii,jj)) ; vv += ii * (ii * XYC(ii,jj)) ;
1301 }
1302 cyvar += (vv - uu*uu/yc[jj] ) ; /* Var(x|y) */
1303 }
1304 }
1305 vv = uu = uyvar = 0.0 ;
1306 for( ii=1 ; ii < nbp ; ii++ ){ /* uu=E(x) vv=E(x^2) */
1307 uu += (ii * xc[ii]) ; vv += ii * (ii * xc[ii]) ;
1308 }
1309 uyvar = vv - uu*uu ; /* Var(x) */
1310 xrat = (uyvar > 0.0) ? cyvar/uyvar /* Var(x|y) / Var(x) */
1311 : 1.0 ;
1312
1313 hmc.d = 1.0 - 0.5*(xrat+yrat) ; /** additive symmetrization **/
1314 return hmc ;
1315 }
1316
1317 /*----------------------------------------------------------------------------*/
1318
1319 static double hpow = 0.3333333333321 ;
INCOR_set_2Dhist_hpower(double hh)1320 void INCOR_set_2Dhist_hpower( double hh )
1321 {
1322 hpow = (hh > 0.0 && hh < 1.0) ? hh : 0.3333333333321 ;
1323 }
1324
1325 static int nhbin = 0 ;
INCOR_set_2Dhist_hbin(int nn)1326 void INCOR_set_2Dhist_hbin( int nn ){ nhbin = nn; }
1327
INCOR_2Dhist_compute_nbin(int ndata)1328 int INCOR_2Dhist_compute_nbin( int ndata )
1329 {
1330 int nbin ;
1331
1332 nbin = (nhbin > 4) ? nhbin : (int)rint(pow((double)ndata,hpow)) ;
1333 if( nbin > 255 ) nbin = 255 ; else if( nbin < 5 ) nbin = 5 ;
1334 return nbin ;
1335 }
1336
1337 /*----------------------------------------------------------------------------*/
1338
INCOR_create_2Dhist(int nbin,float xbot,float xtop,float ybot,float ytop,float xcbot,float xctop,float ycbot,float yctop)1339 INCOR_2Dhist * INCOR_create_2Dhist( int nbin ,
1340 float xbot , float xtop ,
1341 float ybot , float ytop ,
1342 float xcbot, float xctop,
1343 float ycbot, float yctop )
1344 {
1345 INCOR_2Dhist *tdh ; int nbp ;
1346
1347 ENTRY("INCOR_create_2Dhist") ;
1348
1349 if( nbin < 3 ) nbin = 3 ;
1350
1351 tdh = (INCOR_2Dhist *)calloc(1,sizeof(INCOR_2Dhist)) ;
1352
1353 tdh->meth = 0 ; /* undefined as yet */
1354 tdh->nbin = nbin ;
1355 tdh->xxbot = xbot ; tdh->yybot = ybot ;
1356 tdh->xxtop = xtop ; tdh->yytop = ytop ;
1357 tdh->xcbot = xcbot ; tdh->ycbot = ycbot ;
1358 tdh->xctop = xctop ; tdh->yctop = yctop ;
1359
1360 nbp = nbin+1 ;
1361 tdh->xc = (float *)calloc(sizeof(float),nbp) ;
1362 tdh->yc = (float *)calloc(sizeof(float),nbp) ;
1363 tdh->xyc = (float *)calloc(sizeof(float),nbp*nbp) ;
1364 tdh->nww = 0.0f ;
1365
1366 RETURN(tdh) ;
1367 }
1368
1369 /*----------------------------------------------------------------------------*/
1370
INCOR_destroy_2Dhist(INCOR_2Dhist * tin)1371 void INCOR_destroy_2Dhist( INCOR_2Dhist *tin )
1372 {
1373 if( tin == NULL ) return ;
1374 if( tin->xc != NULL ) free(tin->xc) ;
1375 if( tin->yc != NULL ) free(tin->yc) ;
1376 if( tin->xyc != NULL ) free(tin->xyc) ;
1377 free(tin) ;
1378 return ;
1379 }
1380
1381 /*----------------------------------------------------------------------------*/
1382
INCOR_copyover_2Dhist(INCOR_2Dhist * tin,INCOR_2Dhist * tout)1383 void INCOR_copyover_2Dhist( INCOR_2Dhist *tin , INCOR_2Dhist *tout )
1384 {
1385 int nbp ;
1386
1387 ENTRY("INCOR_copyover_2Dhist") ;
1388
1389 if( tin == NULL || tout == NULL || tin == tout ) EXRETURN ;
1390
1391 if( tout->xc != NULL ) free(tout->xc) ;
1392 if( tout->yc != NULL ) free(tout->yc) ;
1393 if( tout->xyc != NULL ) free(tout->xyc) ;
1394
1395 tout->meth = tin->meth ;
1396 tout->nbin = tin->nbin ;
1397 tout->xxbot = tin->xxbot ; tout->yybot = tin->yybot ;
1398 tout->xxtop = tin->xxtop ; tout->yytop = tin->yytop ;
1399 tout->xcbot = tin->xcbot ; tout->ycbot = tin->ycbot ;
1400 tout->xctop = tin->xctop ; tout->yctop = tin->yctop ;
1401 tout->nww = tin->nww ;
1402
1403 nbp = tin->nbin + 1 ;
1404 tout->xc = (float *)malloc(sizeof(float)*nbp) ;
1405 tout->yc = (float *)malloc(sizeof(float)*nbp) ;
1406 tout->xyc = (float *)malloc(sizeof(float)*nbp*nbp) ;
1407
1408 AAmemcpy( tout->xc , tin->xc , sizeof(float)*nbp ) ;
1409 AAmemcpy( tout->yc , tin->yc , sizeof(float)*nbp ) ;
1410 AAmemcpy( tout->xyc, tin->xyc, sizeof(float)*nbp*nbp ) ;
1411
1412 EXRETURN ;
1413 }
1414
1415 /*----------------------------------------------------------------------------*/
1416 /*----------------------------------------------------------------------------*/
1417
INCOR_addto_incomplete_pearson(int n,float * x,float * y,float * w,INCOR_pearson * inpear)1418 void INCOR_addto_incomplete_pearson( int n, float *x, float *y,
1419 float *w, INCOR_pearson *inpear )
1420 {
1421 int ii ; double sx,sxx , sy,syy,sxy , sw ;
1422
1423 if( n <= 0 || x == NULL || y == NULL || inpear == NULL ) return ;
1424
1425 sx = inpear->sx ; sxx = inpear->sxx ;
1426 sy = inpear->sy ; syy = inpear->syy ; sxy = inpear->sxy ; sw = inpear->sw ;
1427
1428 if( w == NULL ){
1429 double xx , yy ;
1430 for( ii=0 ; ii < n ; ii++ ){
1431 xx = (double)x[ii] ; yy = (double)y[ii] ;
1432 sx += xx ; sxx += xx*xx ; sy += yy ; syy += yy*yy ; sxy += xx*yy ;
1433 }
1434 sw += (double)n ;
1435 } else {
1436 double xx , yy , ww ;
1437 for( ii=0 ; ii < n ; ii++ ){
1438 ww = (double)w[ii] ;
1439 if( ww > 0.0 ){
1440 xx = (double)x[ii] ; yy = (double)y[ii] ;
1441 sx += xx*ww ; sxx += xx*xx*ww ;
1442 sy += yy*ww ; syy += yy*yy*ww ; sxy += xx*yy*ww ; sw += ww ;
1443 }
1444 }
1445 }
1446
1447 inpear->npt += n ;
1448 inpear->sx = sx ; inpear->sxx = sxx ;
1449 inpear->sy = sy ; inpear->syy = syy ; inpear->sxy = sxy ; inpear->sw = sw ;
1450
1451 return ;
1452 }
1453
1454 /*----------------------------------------------------------------------------*/
1455
INCOR_destroy_incomplete_pearson(INCOR_pearson * inpear)1456 void INCOR_destroy_incomplete_pearson( INCOR_pearson *inpear )
1457 {
1458 if( inpear != NULL ) free((void *)inpear) ;
1459 }
1460
1461 /*----------------------------------------------------------------------------*/
1462
INCOR_create_incomplete_pearson(void)1463 INCOR_pearson * INCOR_create_incomplete_pearson(void)
1464 {
1465 INCOR_pearson *inpear ;
1466
1467 inpear = (INCOR_pearson *)calloc(1,sizeof(INCOR_pearson)) ;
1468 inpear->sx = 0.0 ; inpear->sxx = 0.0 ;
1469 inpear->sy = 0.0 ; inpear->syy = 0.0 ;
1470 inpear->sxy = 0.0 ; inpear->sw = 0.0 ; inpear->npt = 0 ;
1471
1472 inpear->meth = GA_MATCH_PEARSON_SCALAR ;
1473 return inpear ;
1474 }
1475
1476 /*----------------------------------------------------------------------------*/
1477
INCOR_incomplete_pearson(INCOR_pearson * inpear)1478 double INCOR_incomplete_pearson( INCOR_pearson *inpear )
1479 {
1480 double xv , yv , xy , swi , val ;
1481
1482 if( inpear->sw <= 0.0 ) return 0.0 ;
1483
1484 swi = 1.0 / inpear->sw ;
1485
1486 xv = inpear->sxx - inpear->sx * inpear->sx * swi ;
1487 yv = inpear->syy - inpear->sy * inpear->sy * swi ;
1488 xy = inpear->sxy - inpear->sx * inpear->sy * swi ;
1489
1490 if( xv <= 0.0 || yv <= 0.0 ) return 0.0 ;
1491 val = (xy/sqrt(xv*yv)) ; val = MYatanh(val) ; return val ;
1492 }
1493
1494 /*----------------------------------------------------------------------------*/
1495 /*----------------------------------------------------------------------------*/
1496
1497 #ifndef USE_OMP /*--- Serial code version ---*/
1498
INCOR_addto_incomplete_pearclp_SS(int n,float * x,float * y,float * w,INCOR_pearclp * inpear)1499 void INCOR_addto_incomplete_pearclp_SS( int n, float *x, float *y,
1500 float *w, INCOR_pearclp *inpear )
1501 {
1502 int ii ; double sx,sxx , sy,syy,sxy , sw ;
1503 double xcb,xct , ycb,yct , xmid,ymid ;
1504 double xdb,xdt , ydb,ydt ;
1505
1506 sx = inpear->sx ; sxx = inpear->sxx ;
1507 sy = inpear->sy ; syy = inpear->syy ; sxy = inpear->sxy ; sw = inpear->sw ;
1508
1509 xcb = inpear->xcbot ; xct = inpear->xctop ; xmid = 0.5f*(xcb+xct) ;
1510 ycb = inpear->ycbot ; yct = inpear->yctop ; ymid = 0.5f*(ycb+yct) ;
1511 xdb = inpear->xdbot ; xdt = inpear->xdtop ;
1512 ydb = inpear->ydbot ; ydt = inpear->ydtop ;
1513
1514 if( w == NULL ){
1515 double xx , yy , ww ; int cl ;
1516 for( ii=0 ; ii < n ; ii++ ){
1517 cl = 1 ;
1518 xx = (double)x[ii] ;
1519 if( xx <= xcb ){ xx = xdb; cl++; } else if( xx >= xct ){ xx = xdt; cl++; }
1520 yy = (double)y[ii] ;
1521 if( yy <= ycb ){ yy = ydb; cl++; } else if( yy >= yct ){ yy = ydt; cl++; }
1522 ww = 1.0 / cl ; xx -= xmid ; yy -= ymid ;
1523 sx += xx*ww ; sxx += xx*xx*ww ; sy += yy*ww ; syy += yy*yy*ww ; sxy += xx*yy*ww ; sw += ww ;
1524 }
1525 } else {
1526 double xx , yy , ww ; int cl ;
1527 for( ii=0 ; ii < n ; ii++ ){
1528 ww = (double)w[ii] ;
1529 if( ww > 0.0 ){
1530 cl = 1 ;
1531 xx = (double)x[ii] ;
1532 if( xx <= xcb ){ xx = xdb; cl++; } else if( xx >= xct ){ xx = xdt; cl++; }
1533 yy = (double)y[ii] ;
1534 if( yy <= ycb ){ yy = ydb; cl++; } else if( yy >= yct ){ yy = ydt; cl++; }
1535 ww /= cl ; xx -= xmid ; yy -= ymid ;
1536 sx += xx*ww ; sxx += xx*xx*ww ;
1537 sy += yy*ww ; syy += yy*yy*ww ; sxy += xx*yy*ww ; sw += ww ;
1538 }
1539 }
1540 }
1541
1542 inpear->npt += n ;
1543 inpear->sx = sx ; inpear->sxx = sxx ;
1544 inpear->sy = sy ; inpear->syy = syy ; inpear->sxy = sxy ; inpear->sw = sw ;
1545
1546 return ;
1547 }
1548
1549 /*----------------------------------------------------------------------------*/
1550 #else /*--- Parallel-ized verison of the above ----------------------*/
1551
INCOR_addto_incomplete_pearclp_PP(int n,float * x,float * y,float * w,INCOR_pearclp * inpear)1552 void INCOR_addto_incomplete_pearclp_PP( int n, float *x, float *y,
1553 float *w, INCOR_pearclp *inpear )
1554 {
1555 int jj ; double sx,sxx , sy,syy,sxy , sw ;
1556 double xcb,xct , ycb,yct , xmid,ymid ;
1557 double xdb,xdt , ydb,ydt ;
1558
1559 xcb = inpear->xcbot ; xct = inpear->xctop ; xmid = 0.5f*(xcb+xct) ;
1560 ycb = inpear->ycbot ; yct = inpear->yctop ; ymid = 0.5f*(ycb+yct) ;
1561 xdb = inpear->xdbot ; xdt = inpear->xdtop ;
1562 ydb = inpear->ydbot ; ydt = inpear->ydtop ;
1563
1564 for( jj=0 ; jj < nthmax ; jj++ ){
1565 dhaar[jj] = dhbbr[jj] = dhccr[jj] = dhddr[jj] = dheer[jj] = dhffr[jj] = 0.0 ;
1566 }
1567
1568 if( w == NULL ){
1569
1570 #pragma omp parallel
1571 { double xx , yy , ww ; int cl,ii ; int ith=omp_get_thread_num() ;
1572 double tx=0.0,txx=0.0,ty=0.0,tyy=0.0,txy=0.0,tw=0.0 ;
1573 #pragma omp for
1574 for( ii=0 ; ii < n ; ii++ ){
1575 cl = 1 ;
1576 xx = (double)x[ii] ;
1577 if( xx <= xcb ){ xx = xdb; cl++; } else if( xx >= xct ){ xx = xdt; cl++; }
1578 yy = (double)y[ii] ;
1579 if( yy <= ycb ){ yy = ydb; cl++; } else if( yy >= yct ){ yy = ydt; cl++; }
1580 ww = 1.0 / cl ; xx -= xmid ; yy -= ymid ;
1581 tx += ww*xx ; txx += ww*xx*xx ; ty += ww*yy ;
1582 tyy += ww*yy*yy ; txy += ww*xx*yy ; tw += ww ;
1583 }
1584 #pragma omp critical
1585 { dhaar[ith] = tx ; dhbbr[ith] = txx ; dhccr[ith] = ty ;
1586 dhddr[ith] = tyy ; dheer[ith] = txy ; dhffr[ith] = tw ; }
1587 } /* end parallel */
1588
1589 } else {
1590
1591 #pragma omp parallel
1592 { double xx , yy , ww ; int cl,ii ; int ith=omp_get_thread_num() ;
1593 double tx=0.0,txx=0.0,ty=0.0,tyy=0.0,txy=0.0,tw=0.0 ;
1594 #pragma omp for
1595 for( ii=0 ; ii < n ; ii++ ){
1596 ww = (double)w[ii] ;
1597 if( ww > 0.0 ){
1598 cl = 1 ;
1599 xx = (double)x[ii] ;
1600 if( xx <= xcb ){ xx = xdb; cl++; } else if( xx >= xct ){ xx = xdt; cl++; }
1601 yy = (double)y[ii] ;
1602 if( yy <= ycb ){ yy = ydb; cl++; } else if( yy >= yct ){ yy = ydt; cl++; }
1603 ww /= cl ; xx -= xmid ; yy -= ymid ;
1604 tx += ww*xx ; txx += ww*xx*xx ; ty += ww*yy ;
1605 tyy += ww*yy*yy ; txy += ww*xx*yy ; tw += ww ;
1606 }
1607 }
1608 #pragma omp critical
1609 { dhaar[ith] = tx ; dhbbr[ith] = txx ; dhccr[ith] = ty ;
1610 dhddr[ith] = tyy ; dheer[ith] = txy ; dhffr[ith] = tw ; }
1611 } /* end parallel */
1612 }
1613
1614 /*-- add partial sums from each thread to the results --*/
1615
1616 sx = inpear->sx ; sxx = inpear->sxx ;
1617 sy = inpear->sy ; syy = inpear->syy ;
1618 sxy = inpear->sxy ; sw = inpear->sw ;
1619 for( jj=0 ; jj < nthmax ; jj++ ){
1620 sx += dhaar[jj] ; sxx += dhbbr[jj] ;
1621 sy += dhccr[jj] ; syy += dhddr[jj] ;
1622 sxy += dheer[jj] ; sw += dhffr[jj] ;
1623 }
1624
1625 inpear->npt += n ;
1626 inpear->sx = sx ; inpear->sxx = sxx ;
1627 inpear->sy = sy ; inpear->syy = syy ; inpear->sxy = sxy ; inpear->sw = sw ;
1628
1629 return ;
1630 }
1631 #endif /* USE_OMP */
1632
1633 /*----------------------------------------------------------------------------*/
1634
INCOR_addto_incomplete_pearclp(int n,float * x,float * y,float * w,INCOR_pearclp * inpear)1635 void INCOR_addto_incomplete_pearclp( int n, float *x, float *y,
1636 float *w, INCOR_pearclp *inpear )
1637 {
1638 if( n <= 0 || x == NULL || y == NULL || inpear == NULL ) return ;
1639 #ifdef USE_OMP
1640 INCOR_addto_incomplete_pearclp_PP( n , x,y,w , inpear ) ;
1641 #else
1642 INCOR_addto_incomplete_pearclp_SS( n , x,y,w , inpear ) ;
1643 #endif
1644 return ;
1645 }
1646
1647 /*----------------------------------------------------------------------------*/
1648
INCOR_destroy_incomplete_pearclp(INCOR_pearclp * inpear)1649 void INCOR_destroy_incomplete_pearclp( INCOR_pearclp *inpear )
1650 {
1651 if( inpear != NULL ) free((void *)inpear) ;
1652 }
1653
1654 /*----------------------------------------------------------------------------*/
1655
INCOR_create_incomplete_pearclp(void)1656 INCOR_pearclp * INCOR_create_incomplete_pearclp(void)
1657 {
1658 INCOR_pearclp *inpear ;
1659
1660 inpear = (INCOR_pearclp *)calloc(1,sizeof(INCOR_pearclp)) ;
1661 inpear->sx = 0.0 ; inpear->sxx = 0.0 ;
1662 inpear->sy = 0.0 ; inpear->syy = 0.0 ;
1663 inpear->sxy = 0.0 ; inpear->sw = 0.0 ; inpear->npt = 0 ;
1664
1665 inpear->xcbot = inpear->xctop = inpear->ycbot = inpear->yctop = 0.0 ;
1666 inpear->xdbot = inpear->xdtop = inpear->ydbot = inpear->ydtop = 0.0 ;
1667
1668 inpear->meth = GA_MATCH_PEARCLP_SCALAR ;
1669 return inpear ;
1670 }
1671
1672 /*----------------------------------------------------------------------------*/
1673
INCOR_incomplete_pearclp(INCOR_pearclp * inpear)1674 double INCOR_incomplete_pearclp( INCOR_pearclp *inpear )
1675 {
1676 double xv , yv , xy , swi , val ;
1677
1678 if( inpear->sw <= 0.0 ) return 0.0 ;
1679
1680 swi = 1.0 / inpear->sw ;
1681
1682 xv = inpear->sxx - inpear->sx * inpear->sx * swi ;
1683 yv = inpear->syy - inpear->sy * inpear->sy * swi ;
1684 xy = inpear->sxy - inpear->sx * inpear->sy * swi ;
1685
1686 if( xv <= 0.0 || yv <= 0.0 ) return 0.0 ;
1687 val = (xy/sqrt(xv*yv)) ; val = MYatanh(val) ; return val ;
1688 }
1689
1690 /*============================================================================*/
1691 /* Generic INCOR functions */
1692 /*============================================================================*/
1693
1694 /*----------------------------------------------------------------------------*/
1695 /* Check if a method code is implemented. */
1696
INCOR_check_meth_code(int meth)1697 int INCOR_check_meth_code( int meth )
1698 {
1699 switch( meth ){
1700 case GA_MATCH_PEARSON_SCALAR: return 1 ;
1701
1702 case GA_MATCH_PEARCLP_SCALAR: return 3 ;
1703
1704 case GA_MATCH_MUTINFO_SCALAR:
1705 case GA_MATCH_CORRATIO_SCALAR:
1706 case GA_MATCH_NORMUTIN_SCALAR:
1707 case GA_MATCH_HELLINGER_SCALAR:
1708 case GA_MATCH_CRAT_SADD_SCALAR:
1709 case GA_MATCH_CRAT_USYM_SCALAR: return 2 ; /* uses 2Dhist */
1710
1711 case GA_MATCH_PEARSON_LOCALS:
1712 case GA_MATCH_PEARSON_LOCALA: return 4 ; /* 25 Jun 2014 */
1713 }
1714
1715 return 0 ;
1716 }
1717
1718 /*----------------------------------------------------------------------------*/
1719
1720 static byte *lpc_mask = NULL ;
1721
INCOR_set_lpc_mask(byte * mmm)1722 void INCOR_set_lpc_mask( byte *mmm ){ lpc_mask = mmm ; }
1723
1724 /*----------------------------------------------------------------------------*/
1725
1726 static INCORR_BLOK_set *lpc_ibs = NULL ;
1727
INCOR_reset_lpc_ibs(int nx,int ny,int nz)1728 void INCOR_reset_lpc_ibs( int nx , int ny , int nz )
1729 {
1730 int need_new ;
1731
1732 need_new = (lpc_ibs == NULL) ||
1733 (lpc_ibs->nx != nx ) ||
1734 (lpc_ibs->ny != ny ) ||
1735 (lpc_ibs->nz != nz ) ;
1736
1737 if( need_new ){
1738 if( lpc_ibs != NULL ) destroy_INCORR_BLOK_set(lpc_ibs) ;
1739 if( nx > 3 && ny > 3 && nz > 0 )
1740 lpc_ibs = create_INCORR_BLOK_set( nx,ny,nz, 1.0f,1.0f,1.0f,
1741 lpc_mask, GA_BLOK_TOHD, 6.54321f, 0 ) ;
1742 else
1743 lpc_ibs = NULL ;
1744 } else {
1745 clear_INCORR_BLOK_set( lpc_ibs ) ;
1746 }
1747
1748 return ;
1749 }
1750
1751 /*----------------------------------------------------------------------------*/
1752 /* Create an INCOR object, return a pointer to it.
1753 meth = method to use
1754 mpar = floatvec with parameters:
1755
1756 For meth == GA_MATCH_PEARSON_SCALAR, mpar is not used.
1757
1758 For meth == 2D histogram method, mpar is used as follows:
1759 mpar->ar[0] = nbin
1760 mpar->ar[1] = xbot mpar->ar[2] = xtop
1761 mpar->ar[3] = ybot mpar->ar[4] = ytop
1762 mpar->ar[5] = xcbot mpar->ar[6] = xctop
1763 mpar->ar[7] = ycbot mpar->ar[8] = yctop
1764 * If you have a good idea of how many data points are coming in,
1765 you can use function INCOR_2Dhist_compute_nbin() to compute
1766 a useful nbin from the number of data points.
1767 * If xbot >= xtop or if ybot >= ytop, then these values will be
1768 computed from the first set of data input (and then fixed).
1769 * If we have the relationships
1770 xbot < xcbot < xctop < xtop
1771 AND ybot < ycbot < yctop < ytop
1772 then the first and last bins are 'oddly' sized, and all the interior
1773 bins are evenly spaced. Otherwise, all the bins are evenly spaced.
1774 * If you have some sample data to supply, you can use function
1775 INCOR_2Dhist_minmax() to compute xbot, xtop, ybot, and ytop.
1776 In addition, you can use function INCOR_2Dhist_xyclip() to compute
1777 xcbot, xctop, ycbot, and yctop.
1778
1779 For meth == Local Pearson method, mpar is used as follows:
1780 mpar->ar[0] = nx ar[1] = ny ar[2] = nz (of total volume)
1781 ar[3] = ibot ar[4] = itop (of small volume)
1782 ar[5] = jbot ar[6] = jtop
1783 ar[7] = kbot ar[8] = ktop
1784 *//*--------------------------------------------------------------------------*/
1785
INCOR_create(int meth,floatvec * mpar)1786 void * INCOR_create( int meth , floatvec *mpar )
1787 {
1788 void *vinc = NULL;
1789
1790 ENTRY("INCOR_create") ;
1791
1792 switch( meth ){
1793
1794 case GA_MATCH_PEARSON_SCALAR:
1795 vinc = (void *)INCOR_create_incomplete_pearson() ;
1796 break ;
1797
1798 case GA_MATCH_PEARCLP_SCALAR:{
1799 INCOR_pearclp *pc ;
1800 pc = INCOR_create_incomplete_pearclp() ; vinc = (void *)pc ;
1801 if( mpar != NULL && mpar->nar > 8 ){
1802 pc->xdbot = mpar->ar[1] ; pc->xdtop = mpar->ar[2] ;
1803 pc->ydbot = mpar->ar[3] ; pc->ydtop = mpar->ar[4] ;
1804 pc->xcbot = mpar->ar[5] ; pc->xctop = mpar->ar[6] ;
1805 pc->ycbot = mpar->ar[7] ; pc->yctop = mpar->ar[8] ;
1806 }
1807 }
1808 break ;
1809
1810 case GA_MATCH_MUTINFO_SCALAR: /* methods that use 2D histograms */
1811 case GA_MATCH_CORRATIO_SCALAR:
1812 case GA_MATCH_NORMUTIN_SCALAR:
1813 case GA_MATCH_HELLINGER_SCALAR:
1814 case GA_MATCH_CRAT_SADD_SCALAR:
1815 case GA_MATCH_CRAT_USYM_SCALAR:{
1816 INCOR_2Dhist *tdh ;
1817 int nbin ;
1818 float xbot,xtop, ybot,ytop, xcbot,xctop, ycbot,yctop;
1819 nbin = (mpar == NULL) ? 0 : (int)mpar->ar[0] ;
1820 if( nbin < 0 ) nbin = INCOR_2Dhist_compute_nbin(-nbin) ;
1821 if( nbin > 0 && mpar != NULL && mpar->nar > 8 ){
1822 xbot = mpar->ar[1] ; xtop = mpar->ar[2] ;
1823 ybot = mpar->ar[3] ; ytop = mpar->ar[4] ;
1824 xcbot = mpar->ar[5] ; xctop = mpar->ar[6] ;
1825 ycbot = mpar->ar[7] ; yctop = mpar->ar[8] ;
1826 } else {
1827 xbot = xtop = ybot = ytop = xcbot = xctop = ycbot = yctop = 0.0f ;
1828 }
1829 tdh = INCOR_create_2Dhist( nbin , xbot ,xtop , ybot ,ytop ,
1830 xcbot,xctop, ycbot,yctop ) ;
1831 if( tdh != NULL ){
1832 tdh->meth = meth ; vinc = (void *)tdh ;
1833 }
1834 }
1835 break ;
1836
1837 case GA_MATCH_PEARSON_LOCALS: /* local Pearson stuff [25 Jun 2014] */
1838 case GA_MATCH_PEARSON_LOCALA:{
1839 INCOR_localpearson *ilp ;
1840 int nx,ny,nz ;
1841 if( mpar == NULL || mpar->nar < 9 ) RETURN(vinc) ; /* bad */
1842 ilp = (INCOR_localpearson *)calloc(1,sizeof(INCOR_localpearson)) ;
1843 ilp->meth = meth ;
1844 ilp->iibot = (int)mpar->ar[3] ; ilp->iitop = (int)mpar->ar[4] ;
1845 ilp->jjbot = (int)mpar->ar[5] ; ilp->jjtop = (int)mpar->ar[6] ;
1846 ilp->kkbot = (int)mpar->ar[7] ; ilp->kktop = (int)mpar->ar[8] ;
1847 ilp->nii = ilp->iitop - ilp->iibot + 1 ;
1848 ilp->njj = ilp->jjtop - ilp->jjbot + 1 ;
1849 ilp->nkk = ilp->kktop - ilp->kkbot + 1 ;
1850 nx = (int)mpar->ar[0] ; ny = (int)mpar->ar[1] ; nz = (int)mpar->ar[2] ;
1851 INCOR_reset_lpc_ibs(nx,ny,nz) ; ilp->ibs = lpc_ibs ;
1852 vinc = (void *)ilp ;
1853
1854 #ifdef ALLOW_DEBUG_LPC
1855 debug_lpc = AFNI_yesenv("AFNI_DEBUG_LPC") ;
1856 #endif
1857 }
1858 break ;
1859
1860 }
1861
1862 RETURN(vinc) ;
1863 }
1864
1865 /*----------------------------------------------------------------------------*/
1866 /* Erase an INCOR struct from the Macrocosmic All. */
1867
INCOR_destroy(void * vp)1868 void INCOR_destroy( void *vp )
1869 {
1870 ENTRY("INCOR_destroy") ;
1871
1872 if( vp == NULL ) EXRETURN ;
1873
1874 switch( INCOR_methcode(vp) ){
1875
1876 case GA_MATCH_PEARSON_SCALAR:
1877 INCOR_destroy_incomplete_pearson(vp) ;
1878 break ;
1879
1880 case GA_MATCH_PEARCLP_SCALAR:
1881 INCOR_destroy_incomplete_pearclp(vp) ;
1882 break ;
1883
1884 case GA_MATCH_MUTINFO_SCALAR:
1885 case GA_MATCH_CORRATIO_SCALAR:
1886 case GA_MATCH_NORMUTIN_SCALAR:
1887 case GA_MATCH_HELLINGER_SCALAR:
1888 case GA_MATCH_CRAT_SADD_SCALAR:
1889 case GA_MATCH_CRAT_USYM_SCALAR:
1890 INCOR_destroy_2Dhist(vp) ;
1891 break ;
1892
1893 case GA_MATCH_PEARSON_LOCALS: /* 25 Jun 2014 */
1894 case GA_MATCH_PEARSON_LOCALA:{
1895 if( lpc_mask == NULL ){
1896 INCOR_localpearson *ilp = (INCOR_localpearson *)vp ;
1897 destroy_INCORR_BLOK_set(ilp->ibs) ;
1898 free(ilp) ; lpc_ibs = NULL ;
1899 }
1900 }
1901 break ;
1902
1903 }
1904 EXRETURN ;
1905 }
1906
1907 /*----------------------------------------------------------------------------*/
1908 /* Copy the internal data of an INCOR struct 'vin' over that of 'vout'.
1909 This is used for adding data to vin while keeping vout pristine for re-use.
1910 If vin is NULL, then vout is cleared to the empty INCOR struct.
1911 *//*--------------------------------------------------------------------------*/
1912
INCOR_copyover(void * vin,void * vout)1913 void INCOR_copyover( void *vin , void *vout )
1914 {
1915 int meth ;
1916
1917 ENTRY("INCOR_copyover") ;
1918
1919 if( vout == NULL || vin == vout ) EXRETURN ;
1920
1921 if( vin != NULL ) meth = INCOR_methcode(vin) ;
1922 else meth = INCOR_methcode(vout) ;
1923
1924 switch( meth ){
1925
1926 case GA_MATCH_PEARSON_SCALAR:
1927 if( vin != NULL ){
1928 AAmemcpy( vout , vin , sizeof(INCOR_pearson) ) ;
1929 } else {
1930 INCOR_pearson *vp = (INCOR_pearson *)vout ;
1931 vp->sx = 0.0 ; vp->sxx = 0.0 ;
1932 vp->sy = 0.0 ; vp->syy = 0.0 ;
1933 vp->sxy = 0.0 ; vp->sw = 0.0 ; vp->npt = 0 ;
1934 }
1935 break ;
1936
1937 case GA_MATCH_PEARCLP_SCALAR:
1938 if( vin != NULL ){
1939 AAmemcpy( vout , vin , sizeof(INCOR_pearclp) ) ;
1940 } else {
1941 INCOR_pearclp *vp = (INCOR_pearclp *)vout ;
1942 vp->sx = 0.0 ; vp->sxx = 0.0 ;
1943 vp->sy = 0.0 ; vp->syy = 0.0 ;
1944 vp->sxy = 0.0 ; vp->sw = 0.0 ; vp->npt = 0 ;
1945 vp->xcbot = vp->xctop = vp->ycbot = vp->yctop = 0.0 ;
1946 vp->xdbot = vp->xdtop = vp->ydbot = vp->ydtop = 0.0 ;
1947 }
1948 break ;
1949
1950 case GA_MATCH_MUTINFO_SCALAR:
1951 case GA_MATCH_CORRATIO_SCALAR:
1952 case GA_MATCH_NORMUTIN_SCALAR:
1953 case GA_MATCH_HELLINGER_SCALAR:
1954 case GA_MATCH_CRAT_SADD_SCALAR:
1955 case GA_MATCH_CRAT_USYM_SCALAR:
1956 if( vin != NULL ){
1957 INCOR_copyover_2Dhist( vin , vout ) ;
1958 } else {
1959 INCOR_2Dhist *tdh=(INCOR_2Dhist *)vout ; int nbp=1+tdh->nbin ;
1960 AAmemset(tdh->xc ,0,sizeof(float)*nbp) ;
1961 AAmemset(tdh->yc ,0,sizeof(float)*nbp) ;
1962 AAmemset(tdh->xyc,0,sizeof(float)*nbp*nbp) ;
1963 tdh->nww = 0.0f ;
1964 }
1965 break ;
1966
1967 }
1968
1969 EXRETURN ;
1970 }
1971
1972 /*----------------------------------------------------------------------------*/
1973 /* Add data to an INCOR struct, for use in computing the 'correlation' later. */
1974
INCOR_addto(void * vin,int n,float * x,float * y,float * w)1975 void INCOR_addto( void *vin , int n , float *x , float *y , float *w )
1976 {
1977 ENTRY("INCOR_addto") ;
1978
1979 if( vin == NULL || n <= 0 || x == NULL || y == NULL ) EXRETURN ;
1980
1981 switch( INCOR_methcode(vin) ){
1982
1983 case GA_MATCH_PEARSON_SCALAR:
1984 INCOR_addto_incomplete_pearson( n , x , y , w , vin ) ;
1985 break ;
1986
1987 case GA_MATCH_PEARCLP_SCALAR:
1988 INCOR_addto_incomplete_pearclp( n , x , y , w , vin ) ;
1989 break ;
1990
1991 case GA_MATCH_MUTINFO_SCALAR:
1992 case GA_MATCH_CORRATIO_SCALAR:
1993 case GA_MATCH_NORMUTIN_SCALAR:
1994 case GA_MATCH_HELLINGER_SCALAR:
1995 case GA_MATCH_CRAT_SADD_SCALAR:
1996 case GA_MATCH_CRAT_USYM_SCALAR:
1997 INCOR_addto_2Dhist( vin , n , x , y, w ) ;
1998 break ;
1999
2000 case GA_MATCH_PEARSON_LOCALS: /* 25 Jun 2014 */
2001 case GA_MATCH_PEARSON_LOCALA:{
2002 INCOR_localpearson *ilp = (INCOR_localpearson *)vin ;
2003 if( n != ilp->ibs->nx * ilp->ibs->ny * ilp->ibs->nz )
2004 ERROR_exit("INCOR_addto mismatch for localpearson: n=%d nx=%d ny=%d nz=%d",
2005 n,ilp->ibs->nx,ilp->ibs->ny,ilp->ibs->nz) ;
2006 addto_INCORR_BLOK_set( ilp->ibs ,
2007 0 , ilp->ibs->nx-1 ,
2008 0 , ilp->ibs->ny-1 ,
2009 0 , ilp->ibs->nz-1 , x , y , w ) ;
2010 }
2011 break ;
2012
2013 }
2014
2015 EXRETURN ;
2016 }
2017
2018 /*----------------------------------------------------------------------------*/
2019 /* Evaluate the value of an INCOR struct, as currently constituted.
2020 * The data specified in arguments {n,x,y,w} is used (if n > 0), but is not
2021 added into the vin INCOR struct.
2022 * This feature lets you do correlations with a base set of unchanging
2023 data, plus a set of data that changes, and thence let you optimize the
2024 correlation as that varying data is munged around.
2025 *//*--------------------------------------------------------------------------*/
2026
INCOR_evaluate(void * vin,int n,float * x,float * y,float * w)2027 double INCOR_evaluate( void *vin , int n , float *x , float *y , float *w )
2028 {
2029 void *vtmp=NULL ; double val=0.0 ; int meth ;
2030
2031 ENTRY("INCOR_evaluate") ;
2032
2033 if( vin == NULL ) RETURN(val) ; /* should never transpire */
2034
2035 meth = INCOR_methcode(vin) ;
2036
2037 if( !(meth == GA_MATCH_PEARSON_LOCALA || meth == GA_MATCH_PEARSON_LOCALS) ){
2038 vtmp = INCOR_create( meth , NULL ) ;
2039 INCOR_copyover( vin , vtmp ) ;
2040 INCOR_addto( vtmp , n , x , y , w ) ;
2041 }
2042
2043 switch( meth ){
2044 case GA_MATCH_PEARSON_SCALAR: val = INCOR_incomplete_pearson(vtmp); break;
2045 case GA_MATCH_PEARCLP_SCALAR: val = INCOR_incomplete_pearclp(vtmp); break;
2046 case GA_MATCH_MUTINFO_SCALAR: val = INCOR_mutual_info(vtmp) ; break;
2047 case GA_MATCH_NORMUTIN_SCALAR: val = INCOR_norm_mutinf(vtmp) ; break;
2048 case GA_MATCH_HELLINGER_SCALAR: val = INCOR_hellinger(vtmp) ; break;
2049 case GA_MATCH_CORRATIO_SCALAR: val = INCOR_corr_ratio(vtmp,0) ; break;
2050 case GA_MATCH_CRAT_SADD_SCALAR: val = INCOR_corr_ratio(vtmp,2) ; break;
2051 case GA_MATCH_CRAT_USYM_SCALAR: val = INCOR_corr_ratio(vtmp,1) ; break;
2052
2053 case GA_MATCH_PEARSON_LOCALS: /* 25 Jun 2014 */
2054 case GA_MATCH_PEARSON_LOCALA:{
2055 INCOR_localpearson *ilp = (INCOR_localpearson *)vin ;
2056 if( n != ilp->nii * ilp->njj * ilp->nkk )
2057 ERROR_exit("INCOR_evaluate mismatch for localpearson: n=%d nii=%d njj=%d nkk=%d",
2058 n,ilp->nii,ilp->njj,ilp->nkk ) ;
2059 val = correlate_INCORR_BLOK_set( ilp->ibs ,
2060 ilp->iibot , ilp->iitop ,
2061 ilp->jjbot , ilp->jjtop ,
2062 ilp->kkbot , ilp->kktop , x,y,w ) ;
2063 }
2064 break ;
2065 }
2066
2067 if( vtmp != NULL ) INCOR_destroy(vtmp) ;
2068 RETURN(val) ;
2069 }
2070