1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #undef MAIN
8
9 #include "afni.h"
10
11 /*-----------------------------------------------------------------------------
12 Register a function to be called from the "FIM+" menu.
13 menu_name = label for menu
14 nbrik = number of values returned by user_func
15 user_func = function to call for each data time series
16 user_data = extra data to go to user_func
17
18 void user_func( int n, float *ts, void *user_data, int nbrik, void *v )
19
20 n = length of time series
21 ts[] = time series data
22 user_data = whatever you want (same as above)
23 nbrik = same as above
24 val = (float *) v = output array of values
25 val[] = user_func should fill val[0..nbrik-1] with values to save
26 (presumably computed from ts[])
27
28 user_func should not destroy the data in ts[], since it will be re-used
29 if more than one fimfunc is used at a time.
30
31 Before the first call with time series data, user_func will be called like so:
32 user_func( n,NULL,user_data,nbrik, FIMdata *fd ) ;
33 where FIMdata is the type
34 typedef struct {
35 MRI_IMAGE *ref_ts , *ort_ts ;
36 int nvox , ignore , polort ;
37 } FIMdata ;
38 and fd is used to pass in the parameters set by the user:
39 fd->ref_ts = float image of reference ("ideal") time series
40 fd->ref_ts->nx = length of data
41 fd->ref_ts->ny = number of ref time series
42 MRI_FLOAT_PTR(fd->ref_ts) = pointer to ref timeseries data
43 fd->ort_ts = similarly for the orts (but this may be NULL)
44 fd->nvox = number of voxels that will be passed in
45 = number of "normal" calls to user_data
46 fd->ignore = ignore setting (but ignored data is included in ts[])
47 fd->polort = polort setting (but no detrending is done to ts[])
48 N.B.: FIMdata is typedef-ed in afni.h.
49
50 After the last (nvox-th) time series is processed, user_func will be called
51 one last time, with the form
52 user_func( -k,NULL,user_data,nbrik, dset ) ;
53 where dset is a pointer to the dataset (THD_3dim_dataset *) that has
54 just been constructed. The value k (negative of the first argument) will
55 be the index of the first sub-brick in the dataset - it contains the data
56 from val[0]. If desired, you can use this final call to set sub-brick
57 parameters for sub-bricks k through k+nbrik-1 (for val[0] .. val[nbrik-1]).
58
59 Both special calls have the second argument NULL, which is how user_func
60 can distinguish them from a normal call with timeseries data in ts[].
61
62 The initialization call will have the first argument be the number of
63 points in the timeseries to come ('n') - this will be positive. The
64 final call will have the first argument be the negative of the sub-brick
65 index of the first sub-brick created by this routine. This value will
66 be non-positive (will be 0 or less), which is how the initialization
67 and final calls can be distinguished. In each of these cases, the
68 last argument must be properly cast to the correct pointer type before
69 being used.
70
71 A sample fimfunc for the Spearman rank correlation, spearman_fimfunc(),
72 is given at the end of this file.
73 --------------------------------------------------------------------------------*/
74
AFNI_register_fimfunc(char * menu_name,int nbrik,generic_func * user_func,void * user_data)75 void AFNI_register_fimfunc( char *menu_name , int nbrik ,
76 generic_func *user_func , void *user_data )
77 {
78 MCW_function_list *rlist = &(GLOBAL_library.registered_fim) ;
79 int num = rlist->num ;
80
81 ENTRY("AFNI_register_fimfunc") ;
82
83 if( menu_name == NULL || menu_name[0] == '\0' ||
84 nbrik <= 0 || user_func == NULL ) EXRETURN ; /* bad inputs! */
85
86 if( num == sizeof(int) ) EXRETURN ; /* too many already! */
87
88 if( num == 0 ){
89 rlist->flags=NULL; rlist->labels=NULL; rlist->funcs=NULL;
90 rlist->func_data=NULL; rlist->func_code=NULL; rlist->func_init=NULL;
91 }
92
93 rlist->flags = (int *) XtRealloc( (char *)rlist->flags, sizeof(int)*(num+1) ) ;
94
95 rlist->labels = (char **) XtRealloc( (char *)rlist->labels ,
96 sizeof(char *)*(num+1) ) ;
97
98 rlist->funcs = (generic_func **) XtRealloc( (char *)rlist->funcs ,
99 sizeof(generic_func *)*(num+1) ) ;
100
101 rlist->func_data = (void **) XtRealloc( (char *)rlist->func_data ,
102 sizeof(void *)*(num+1) ) ;
103
104 rlist->func_code = (int *) XtRealloc( (char *)rlist->func_code, sizeof(int)*(num+1) ) ;
105
106 rlist->func_init = (generic_func **) XtRealloc( (char *)rlist->func_init ,
107 sizeof(generic_func *)*(num+1) ) ;
108
109 rlist->flags[num] = nbrik ;
110 rlist->labels[num] = XtNewString(menu_name) ;
111 rlist->funcs[num] = user_func ;
112 rlist->func_data[num] = user_data ;
113 rlist->func_code[num] = FUNC_FIM ;
114 rlist->func_init[num] = NULL ;
115
116 rlist->num = num+1 ;
117 EXRETURN ;
118 }
119
120 /*******************************************************************************
121 Code for the Spearman rank and quandrant correlation fimfuncs
122 ********************************************************************************/
123
spearman_rank_manycorr(int n,float * x,int nr,float * rv,float ** r)124 static float spearman_rank_manycorr( int n , float *x ,
125 int nr, float *rv, float **r )
126 {
127 register int ii ;
128 register float ss ;
129 int jj ;
130 float bb , xv ;
131
132 if( nr == 1 ) return spearman_rank_corr( n,x,rv[0],r[0] ) ;
133
134 xv = spearman_rank_prepare( n , x ) ; if( xv <= 0.0 ) return 0.0 ;
135
136 for( jj=0,bb=0.0 ; jj < nr ; jj++ ){
137 for( ii=0,ss=0.0 ; ii < n ; ii++ ) ss += x[ii] * r[jj][ii] ;
138 ss /= sqrt(rv[jj]*xv) ;
139 if( fabs(ss) > fabs(bb) ) bb = ss ;
140 }
141
142 return bb ;
143 }
144
quadrant_manycorr(int n,float * x,int nr,float * rv,float ** r)145 static float quadrant_manycorr( int n , float *x ,
146 int nr, float *rv, float **r )
147 {
148 register int ii ;
149 register float ss ;
150 int jj ;
151 float bb , xv ;
152
153 if( nr == 1 ) return quadrant_corr( n,x,rv[0],r[0] ) ;
154
155 xv = quadrant_corr_prepare( n , x ) ; if( xv <= 0.0 ) return 0.0 ;
156
157 for( jj=0,bb=0.0 ; jj < nr ; jj++ ){
158 for( ii=0,ss=0.0 ; ii < n ; ii++ ) ss += x[ii] * r[jj][ii] ;
159 ss /= sqrt(rv[jj]*xv) ;
160 if( fabs(ss) > fabs(bb) ) bb = ss ;
161 }
162
163 return bb ;
164 }
165
166 /*--------------------------------------------------------------------------
167 A sample fimfunc to compute the Spearman rank correlation
168 ----------------------------------------------------------------------------*/
169
spearman_fimfunc(int n,float * ts,void * ud,int nb,void * val)170 void spearman_fimfunc( int n, float *ts, void *ud, int nb, void *val )
171 {
172 static float *tsc=NULL , *refc=NULL , *refv=NULL , **rr ;
173 static int *keep=NULL ;
174 static int ntsc , nref , nkeep , polort,ignore , ncall;
175
176 int ii , kk ;
177 float *v ;
178
179 /*--- handle special cases ---*/
180
181 if( ts == NULL ){
182
183 /*--- the initialization call ---*/
184
185 if( n > 0 ){
186
187 FIMdata *fd = (FIMdata *) val ;
188
189 polort = fd->polort ; /* not used here */
190 ignore = fd->ignore ;
191 ncall = 0 ; /* how many times we have been called */
192
193 /* make workspace for copy of ts data when it arrives */
194
195 ntsc = n ;
196 if( tsc != NULL ) free(tsc) ;
197 tsc = (float *) malloc(sizeof(float)*ntsc) ;
198
199 /* make copy of ref data */
200
201 nref = fd->ref_ts->ny ;
202 if( refc != NULL ) free(refc) ;
203 if( refv != NULL ) free(refv) ;
204 if( keep != NULL ) free(keep) ;
205 if( rr != NULL ) free(rr) ;
206 refc = (float *) malloc(sizeof(float)*ntsc*nref) ; /* copy of ref */
207 refv = (float *) malloc(sizeof(float)*nref) ; /* rank variances */
208 keep = (int *) malloc(sizeof(int)*ntsc) ; /* keeper indices */
209 rr = (float **)malloc(sizeof(float *)*nref) ; /* convenience ptrs */
210
211 for( kk=0 ; kk < nref ; kk++ ){
212 rr[kk] = refc+kk*ntsc ; /* compute ptr */
213 memcpy( rr[kk] , /* copy data */
214 MRI_FLOAT_PTR(fd->ref_ts) + kk*fd->ref_ts->nx ,
215 sizeof(float)*ntsc ) ;
216 }
217
218 /* mark those places we will keep (where ref is OK) */
219
220 for( nkeep=0,ii=ignore ; ii < ntsc ; ii++ ){ /* for each time point */
221 for( kk=0 ; kk < nref ; kk++ ) /* check each ref */
222 if( rr[kk][ii] >= WAY_BIG ) break ;
223
224 if( kk == nref ) keep[nkeep++] = ii ; /* mark if all are OK */
225 }
226
227 /* compress ref, eliminating non-keepers */
228
229 if( nkeep < ntsc ){
230 for( ii=0 ; ii < nkeep ; ii++ ){
231 for( kk=0 ; kk < nref ; kk++ )
232 rr[kk][ii] = rr[kk][keep[ii]] ; /* copy backwards */
233 }
234 }
235
236 /* prepare each ref vector for rank processing */
237
238 for( kk=0 ; kk < nref ; kk++ )
239 refv[kk] = spearman_rank_prepare( nkeep , rr[kk] ) ;
240
241 #if 0
242 fprintf(stderr,"spearman_fimfunc: initialize ntsc=%d nkeep=%d nref=%d\n",
243 ntsc,nkeep,nref);
244 #endif
245
246 return ;
247
248 /*--- the ending call ---*/
249
250 } else {
251 THD_3dim_dataset *dset = (THD_3dim_dataset *) val ;
252 int kb = -n ;
253
254 free(tsc) ; tsc = NULL ; /* free workspaces */
255 free(refc) ; refc = NULL ;
256 free(refv) ; refv = NULL ;
257 free(keep) ; keep = NULL ;
258 free(rr) ; rr = NULL ;
259 rank_order_float(0,NULL) ;
260
261 /* edit sub-brick statistical parameters and name */
262
263 EDIT_BRICK_TO_FICO( dset , kb , nkeep , (nref==1)?1:2 , 1 ) ;
264 EDIT_BRICK_LABEL( dset , kb , "Spearman CC" ) ;
265
266 #if 0
267 fprintf(stderr,"spearman_fimfunc: finalize with kb=%d\n",kb);
268 #endif
269
270 return ;
271 }
272 }
273
274 /*--- the normal case [with data] ---*/
275
276 ncall++ ;
277 #if 0
278 if(ncall%1000==0) fprintf(stderr,"spearman_fimfunc: ncall=%d\n",ncall);
279 #endif
280
281 /* copy input timeseries data (since we will mangle it) */
282
283 if( nkeep < ntsc )
284 for( ii=0 ; ii < nkeep ; ii++ ) tsc[ii] = ts[keep[ii]]; /* the hard way */
285 else
286 memcpy(tsc,ts,sizeof(float)*ntsc) ; /* the easy way */
287
288 /* compute our single result */
289
290 v = (float *) val ;
291 v[0] = spearman_rank_manycorr( nkeep , tsc , nref , refv , rr ) ;
292
293 return ;
294 }
295
296 /*--------------------------------------------------------------------------
297 A sample fimfunc to compute the quadrant correlation
298 ----------------------------------------------------------------------------*/
299
quadrant_fimfunc(int n,float * ts,void * ud,int nb,void * val)300 void quadrant_fimfunc( int n, float *ts, void *ud, int nb, void *val )
301 {
302 static float *tsc=NULL , *refc=NULL , *refv=NULL , **rr ;
303 static int *keep=NULL ;
304 static int ntsc , nref , nkeep , polort,ignore , ncall;
305
306 int ii , kk ;
307 float *v ;
308
309 /*--- handle special cases ---*/
310
311 if( ts == NULL ){
312
313 /*--- the initialization call ---*/
314
315 if( n > 0 ){
316
317 FIMdata *fd = (FIMdata *) val ;
318
319 polort = fd->polort ; /* not used here */
320 ignore = fd->ignore ;
321 ncall = 0 ; /* how many times we have been called */
322
323 /* make workspace for copy of ts data when it arrives */
324
325 ntsc = n ;
326 if( tsc != NULL ) free(tsc) ;
327 tsc = (float *) malloc(sizeof(float)*ntsc) ;
328
329 /* make copy of ref data */
330
331 nref = fd->ref_ts->ny ;
332 if( refc != NULL ) free(refc) ;
333 if( refv != NULL ) free(refv) ;
334 if( keep != NULL ) free(keep) ;
335 if( rr != NULL ) free(rr) ;
336 refc = (float *) malloc(sizeof(float)*ntsc*nref) ; /* copy of ref */
337 refv = (float *) malloc(sizeof(float)*nref) ; /* rank variances */
338 keep = (int *) malloc(sizeof(int)*ntsc) ; /* keeper indices */
339 rr = (float **)malloc(sizeof(float *)*nref) ; /* convenience ptrs */
340
341 for( kk=0 ; kk < nref ; kk++ ){
342 rr[kk] = refc+kk*ntsc ; /* compute ptr */
343 memcpy( rr[kk] , /* copy data */
344 MRI_FLOAT_PTR(fd->ref_ts) + kk*fd->ref_ts->nx ,
345 sizeof(float)*ntsc ) ;
346 }
347
348 /* mark those places we will keep (where ref is OK) */
349
350 for( nkeep=0,ii=ignore ; ii < ntsc ; ii++ ){ /* for each time point */
351 for( kk=0 ; kk < nref ; kk++ ) /* check each ref */
352 if( rr[kk][ii] >= WAY_BIG ) break ;
353
354 if( kk == nref ) keep[nkeep++] = ii ; /* mark if all are OK */
355 }
356
357 /* compress ref, eliminating non-keepers */
358
359 if( nkeep < ntsc ){
360 for( ii=0 ; ii < nkeep ; ii++ ){
361 for( kk=0 ; kk < nref ; kk++ )
362 rr[kk][ii] = rr[kk][keep[ii]] ; /* copy backwards */
363 }
364 }
365
366 /* prepare each ref vector for rank processing */
367
368 for( kk=0 ; kk < nref ; kk++ )
369 refv[kk] = quadrant_corr_prepare( nkeep , rr[kk] ) ;
370
371 #if 0
372 fprintf(stderr,"quadrant_fimfunc: initialize ntsc=%d nkeep=%d nref=%d\n",
373 ntsc,nkeep,nref);
374 #endif
375
376 return ;
377
378 /*--- the ending call ---*/
379
380 } else {
381 THD_3dim_dataset *dset = (THD_3dim_dataset *) val ;
382 int kb = -n ;
383
384 free(tsc) ; tsc = NULL ; /* free workspaces */
385 free(refc) ; refc = NULL ;
386 free(refv) ; refv = NULL ;
387 free(keep) ; keep = NULL ;
388 free(rr) ; rr = NULL ;
389 rank_order_float(0,NULL) ;
390
391 /* edit sub-brick statistical parameters and name */
392
393 EDIT_BRICK_TO_FICO( dset , kb , nkeep , (nref==1)?1:2 , 1 ) ;
394 EDIT_BRICK_LABEL( dset , kb , "Quadrant CC" ) ;
395
396 #if 0
397 fprintf(stderr,"quadrant_fimfunc: finalize with kb=%d\n",kb);
398 #endif
399
400 return ;
401 }
402 }
403
404 /*--- the normal case [with data] ---*/
405
406 ncall++ ;
407 #if 0
408 if(ncall%1000==0) fprintf(stderr,"quadrant_fimfunc: ncall=%d\n",ncall);
409 #endif
410
411 /* copy input timeseries data (since we will mangle it) */
412
413 if( nkeep < ntsc )
414 for( ii=0 ; ii < nkeep ; ii++ ) tsc[ii] = ts[keep[ii]]; /* the hard way */
415 else
416 memcpy(tsc,ts,sizeof(float)*ntsc) ; /* the easy way */
417
418 /* compute our single result */
419
420 v = (float *) val ;
421 v[0] = quadrant_manycorr( nkeep , tsc , nref , refv , rr ) ;
422
423 return ;
424 }
425