00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 #undef MAIN
00008 
00009 #include "afni.h"
00010 
00011 /*-----------------------------------------------------------------------------
00012    Register a function to be called from the "FIM+" menu.
00013      menu_name = label for menu
00014      nbrik     = number of values returned by user_func
00015      user_func = function to call for each data time series
00016      user_data = extra data to go to user_func
00017 
00018    void user_func( int n, float *ts, void *user_data, int nbrik, void *v )
00019 
00020      n         = length of time series
00021      ts[]      = time series data
00022      user_data = whatever you want (same as above)
00023      nbrik     = same as above
00024      val       = (float *) v = output array of values
00025      val[]     = user_func should fill val[0..nbrik-1] with values to save
00026                  (presumably computed from ts[])
00027 
00028    user_func should not destroy the data in ts[], since it will be re-used
00029    if more than one fimfunc is used at a time.
00030 
00031    Before the first call with time series data, user_func will be called like so:
00032      user_func( n,NULL,user_data,nbrik, FIMdata *fd ) ;
00033    where FIMdata is the type
00034      typedef struct {
00035         MRI_IMAGE *ref_ts , *ort_ts ;
00036         int nvox , ignore , polort ;
00037      } FIMdata ;
00038    and fd is used to pass in the parameters set by the user:
00039       fd->ref_ts = float image of reference ("ideal") time series
00040                    fd->ref_ts->nx = length of data
00041                    fd->ref_ts->ny = number of ref time series
00042                    MRI_FLOAT_PTR(fd->ref_ts) = pointer to ref timeseries data
00043       fd->ort_ts = similarly for the orts (but this may be NULL)
00044       fd->nvox   = number of voxels that will be passed in
00045                  = number of "normal" calls to user_data
00046       fd->ignore = ignore setting (but ignored data is included in ts[])
00047       fd->polort = polort setting (but no detrending is done to ts[])
00048    N.B.: FIMdata is typedef-ed in afni.h.
00049 
00050    After the last (nvox-th) time series is processed, user_func will be called
00051    one last time, with the form
00052      user_func( -k,NULL,user_data,nbrik, dset ) ;
00053    where dset is a pointer to the dataset (THD_3dim_dataset *) that has
00054    just been constructed.  The value k (negative of the first argument) will
00055    be the index of the first sub-brick in the dataset - it contains the data
00056    from val[0].  If desired, you can use this final call to set sub-brick
00057    parameters for sub-bricks k through k+nbrik-1 (for val[0] .. val[nbrik-1]).
00058 
00059    Both special calls have the second argument NULL, which is how user_func
00060    can distinguish them from a normal call with timeseries data in ts[].
00061 
00062    The initialization call will have the first argument be the number of
00063    points in the timeseries to come ('n') - this will be positive.  The
00064    final call will have the first argument be the negative of the sub-brick
00065    index of the first sub-brick created by this routine.  This value will
00066    be non-positive (will be 0 or less), which is how the initialization
00067    and final calls can be distinguished.  In each of these cases, the
00068    last argument must be properly cast to the correct pointer type before
00069    being used.
00070 
00071    A sample fimfunc for the Spearman rank correlation, spearman_fimfunc(),
00072    is given at the end of this file.
00073 --------------------------------------------------------------------------------*/
00074 
00075 void AFNI_register_fimfunc( char *menu_name , int nbrik ,
00076                             generic_func *user_func , void *user_data )
00077 {
00078    MCW_function_list *rlist = &(GLOBAL_library.registered_fim) ;
00079    int num = rlist->num ;
00080 
00081 ENTRY("AFNI_register_fimfunc") ;
00082 
00083    if( menu_name == NULL || menu_name[0] == '\0' ||
00084        nbrik <= 0        || user_func == NULL      ) EXRETURN ; /* bad inputs! */
00085 
00086    if( num == sizeof(int) ) EXRETURN ;  /* too many already! */
00087 
00088    if( num == 0 ){
00089      rlist->flags=NULL; rlist->labels=NULL; rlist->funcs=NULL;
00090      rlist->func_data=NULL; rlist->func_code=NULL; rlist->func_init=NULL;
00091    }
00092 
00093    rlist->flags = (int *) XtRealloc( (char *)rlist->flags, sizeof(int)*(num+1) ) ;
00094 
00095    rlist->labels = (char **) XtRealloc( (char *)rlist->labels ,
00096                                         sizeof(char *)*(num+1) ) ;
00097 
00098    rlist->funcs = (generic_func **) XtRealloc( (char *)rlist->funcs ,
00099                                                sizeof(generic_func *)*(num+1) ) ;
00100 
00101    rlist->func_data = (void **) XtRealloc( (char *)rlist->func_data ,
00102                                            sizeof(void *)*(num+1) ) ;
00103 
00104    rlist->func_code = (int *) XtRealloc( (char *)rlist->func_code, sizeof(int)*(num+1) ) ;
00105 
00106    rlist->func_init = (generic_func **) XtRealloc( (char *)rlist->func_init ,
00107                                                    sizeof(generic_func *)*(num+1) ) ;
00108 
00109    rlist->flags[num]     = nbrik ;
00110    rlist->labels[num]    = XtNewString(menu_name) ;
00111    rlist->funcs[num]     = user_func ;
00112    rlist->func_data[num] = user_data ;
00113    rlist->func_code[num] = FUNC_FIM  ;
00114    rlist->func_init[num] = NULL ;
00115 
00116    rlist->num = num+1 ;
00117    EXRETURN ;
00118 }
00119 
00120 /*******************************************************************************
00121   Code for the Spearman rank and quandrant correlation fimfuncs
00122 ********************************************************************************/
00123 
00124 /*------------------------------------------------------------------------------
00125    Rank-order a float array, with ties getting the average rank.
00126    The output overwrites the input.
00127 --------------------------------------------------------------------------------*/
00128 
00129 static void rank_order_float( int n , float *a )
00130 {
00131    register int ii , ns , n1 , ib ;
00132    static int    nb = 0 ;
00133    static int   *b  = NULL ;  /* workspaces */
00134    static float *c  = NULL ;
00135    float cs ;
00136 
00137    /*- handle special cases -*/
00138 
00139    if( a == NULL ){
00140       if( b != NULL ){ free(b); free(c); b=NULL ; c=NULL; nb=0; }  /* free workspaces */
00141       return ;
00142    }
00143 
00144    if( n < 1 ) return ;                     /* meaningless input */
00145    if( n == 1 ){ a[0] = 0.0 ; return ; }    /* only one point!? */
00146 
00147    /*- make workspaces, if needed -*/
00148 
00149    if( nb < n ){
00150       if( b != NULL ){ free(b); free(c); }
00151       b  = (int   *) malloc(sizeof(int  )*n) ;
00152       c  = (float *) malloc(sizeof(float)*n) ;
00153       nb = n ;
00154    }
00155 
00156    for( ii=0 ; ii < n ; ii++ ) c[ii] = b[ii] = ii ;
00157 
00158    /*- sort input, carrying b along -*/
00159 
00160    qsort_floatint( n , a , b ) ;  /* see cs_sort_fi.c */
00161 
00162    /* compute ranks into c[] */
00163 
00164    n1 = n-1 ;
00165    for( ii=0 ; ii < n1 ; ii++ ){
00166       if( a[ii] == a[ii+1] ){                  /* handle ties */
00167          cs = 2*ii+1 ; ns = 2 ; ib=ii ; ii++ ;
00168          while( ii < n1 && a[ii] == a[ii+1] ){ ii++ ; ns++ ; cs += ii ; }
00169          for( cs/=ns ; ib <= ii ; ib++ ) c[ib] = cs ;
00170       }
00171    }
00172 
00173    for( ii=0 ; ii < n ; ii++ ) a[b[ii]] = c[ii] ;
00174 
00175    return ;
00176 }
00177 
00178 /*---------------------------------------------------------------------------
00179    Rank orders a[], subtracts the mean rank, and returns the sum-of-squares
00180 -----------------------------------------------------------------------------*/
00181 
00182 static float spearman_rank_prepare( int n , float *a )
00183 {
00184    register int ii ;
00185    register float rb , rs ;
00186 
00187    rank_order_float( n , a ) ;
00188 
00189    rb = 0.5*(n-1) ; rs=0.0 ;
00190    for( ii=0 ; ii < n ; ii++ ){
00191       a[ii] -= rb ;
00192       rs    += a[ii]*a[ii] ;
00193    }
00194 
00195    return rs ;
00196 }
00197 
00198 static float quadrant_corr_prepare( int n , float *a )
00199 {
00200    register int ii ;
00201    register float rb , rs ;
00202 
00203    rank_order_float( n , a ) ;
00204 
00205    rb = 0.5*(n-1) ; rs=0.0 ;
00206    for( ii=0 ; ii < n ; ii++ ){
00207       a[ii] = (a[ii] > rb) ? 1.0
00208                            : (a[ii] < rb) ? -1.0 : 0.0 ;
00209       rs    += a[ii]*a[ii] ;
00210    }
00211 
00212    return rs ;
00213 }
00214 
00215 /*-----------------------------------------------------------------------------
00216     To correlate x[] with r[], do
00217       rv = spearman_rank_prepare(n,r) ;
00218     then
00219       corr = spearman_rank_corr(n,x,rv,r) ;
00220     Note that these 2 routines are destructive (r and x are replaced by ranks)
00221 -------------------------------------------------------------------------------*/
00222 
00223 static float spearman_rank_corr( int n , float *x , float rv , float *r )
00224 {
00225    register int ii ;
00226    register float ss ; float xv ;
00227 
00228    xv = spearman_rank_prepare( n , x ) ; if( xv <= 0.0 ) return 0.0 ;
00229 
00230    for( ii=0,ss=0.0 ; ii < n ; ii++ ) ss += x[ii] * r[ii] ;
00231 
00232    return ( ss/sqrt(rv*xv) ) ;
00233 }
00234 
00235 static float spearman_rank_manycorr( int n , float *x ,
00236                                      int nr, float *rv, float **r )
00237 {
00238    register int ii ;
00239    register float ss ;
00240    int jj ;
00241    float bb , xv ;
00242 
00243    if( nr == 1 ) return spearman_rank_corr( n,x,rv[0],r[0] ) ;
00244 
00245    xv = spearman_rank_prepare( n , x ) ; if( xv <= 0.0 ) return 0.0 ;
00246 
00247    for( jj=0,bb=0.0 ; jj < nr ; jj++ ){
00248       for( ii=0,ss=0.0 ; ii < n ; ii++ ) ss += x[ii] * r[jj][ii] ;
00249       ss /= sqrt(rv[jj]*xv) ;
00250       if( fabs(ss) > fabs(bb) ) bb = ss ;
00251    }
00252 
00253    return bb ;
00254 }
00255 
00256 static float quadrant_corr( int n , float *x , float rv , float *r )
00257 {
00258    register int ii ;
00259    register float ss ; float xv ;
00260 
00261    xv = quadrant_corr_prepare( n , x ) ; if( xv <= 0.0 ) return 0.0 ;
00262 
00263    for( ii=0,ss=0.0 ; ii < n ; ii++ ) ss += x[ii] * r[ii] ;
00264 
00265    return ( ss/sqrt(rv*xv) ) ;
00266 }
00267 
00268 static float quadrant_manycorr( int n , float *x ,
00269                                 int nr, float *rv, float **r )
00270 {
00271    register int ii ;
00272    register float ss ;
00273    int jj ;
00274    float bb , xv ;
00275 
00276    if( nr == 1 ) return quadrant_corr( n,x,rv[0],r[0] ) ;
00277 
00278    xv = quadrant_corr_prepare( n , x ) ; if( xv <= 0.0 ) return 0.0 ;
00279 
00280    for( jj=0,bb=0.0 ; jj < nr ; jj++ ){
00281       for( ii=0,ss=0.0 ; ii < n ; ii++ ) ss += x[ii] * r[jj][ii] ;
00282       ss /= sqrt(rv[jj]*xv) ;
00283       if( fabs(ss) > fabs(bb) ) bb = ss ;
00284    }
00285 
00286    return bb ;
00287 }
00288 
00289 /*--------------------------------------------------------------------------
00290   A sample fimfunc to compute the Spearman rank correlation
00291 ----------------------------------------------------------------------------*/
00292 
00293 void spearman_fimfunc( int n, float *ts, void *ud, int nb, void *val )
00294 {
00295    static float *tsc=NULL , *refc=NULL , *refv=NULL , **rr ;
00296    static int   *keep=NULL ;
00297    static int    ntsc , nref , nkeep , polort,ignore , ncall;
00298 
00299    int ii , kk ;
00300    float *v ;
00301 
00302    /*--- handle special cases ---*/
00303 
00304    if( ts == NULL ){
00305 
00306       /*--- the initialization call ---*/
00307 
00308       if( n > 0 ){
00309 
00310          FIMdata *fd = (FIMdata *) val ;
00311 
00312          polort = fd->polort ;  /* not used here */
00313          ignore = fd->ignore ;
00314          ncall  = 0 ;           /* how many times we have been called */
00315 
00316          /* make workspace for copy of ts data when it arrives */
00317 
00318          ntsc = n ;
00319          if( tsc != NULL ) free(tsc) ;
00320          tsc = (float *) malloc(sizeof(float)*ntsc) ;
00321 
00322          /* make copy of ref data */
00323 
00324          nref = fd->ref_ts->ny ;
00325          if( refc != NULL ) free(refc) ;
00326          if( refv != NULL ) free(refv) ;
00327          if( keep != NULL ) free(keep) ;
00328          if( rr   != NULL ) free(rr) ;
00329          refc = (float *) malloc(sizeof(float)*ntsc*nref) ;  /* copy of ref   */
00330          refv = (float *) malloc(sizeof(float)*nref) ;      /* rank variances */
00331          keep = (int *)   malloc(sizeof(int)*ntsc) ;       /* keeper indices  */
00332          rr   = (float **)malloc(sizeof(float *)*nref) ;  /* convenience ptrs */
00333 
00334          for( kk=0 ; kk < nref ; kk++ ){
00335             rr[kk] = refc+kk*ntsc ;                       /* compute ptr */
00336             memcpy( rr[kk] ,                              /* copy data  */
00337                     MRI_FLOAT_PTR(fd->ref_ts) + kk*fd->ref_ts->nx ,
00338                     sizeof(float)*ntsc  ) ;
00339          }
00340 
00341          /* mark those places we will keep (where ref is OK) */
00342 
00343          for( nkeep=0,ii=ignore ; ii < ntsc ; ii++ ){ /* for each time point */
00344             for( kk=0 ; kk < nref ; kk++ )            /* check each ref     */
00345                if( rr[kk][ii] >= WAY_BIG ) break ;
00346 
00347             if( kk == nref ) keep[nkeep++] = ii ;     /* mark if all are OK */
00348          }
00349 
00350          /* compress ref, eliminating non-keepers */
00351 
00352          if( nkeep < ntsc ){
00353             for( ii=0 ; ii < nkeep ; ii++ ){
00354                for( kk=0 ; kk < nref ; kk++ )
00355                   rr[kk][ii] = rr[kk][keep[ii]] ;  /* copy backwards */
00356             }
00357          }
00358 
00359          /* prepare each ref vector for rank processing */
00360 
00361          for( kk=0 ; kk < nref ; kk++ )
00362             refv[kk] = spearman_rank_prepare( nkeep , rr[kk] ) ;
00363 
00364 #if 0
00365 fprintf(stderr,"spearman_fimfunc: initialize ntsc=%d nkeep=%d nref=%d\n",
00366         ntsc,nkeep,nref);
00367 #endif
00368 
00369          return ;
00370 
00371       /*--- the ending call ---*/
00372 
00373       } else {
00374          THD_3dim_dataset *dset = (THD_3dim_dataset *) val ;
00375          int kb = -n ;
00376 
00377          free(tsc)  ; tsc  = NULL ;  /* free workspaces */
00378          free(refc) ; refc = NULL ;
00379          free(refv) ; refv = NULL ;
00380          free(keep) ; keep = NULL ;
00381          free(rr)   ; rr   = NULL ;
00382          rank_order_float(0,NULL) ;
00383 
00384          /* edit sub-brick statistical parameters and name */
00385 
00386          EDIT_BRICK_TO_FICO( dset , kb , nkeep , (nref==1)?1:2 , 1 ) ;
00387          EDIT_BRICK_LABEL( dset , kb , "Spearman CC" ) ;
00388 
00389 #if 0
00390 fprintf(stderr,"spearman_fimfunc: finalize with kb=%d\n",kb);
00391 #endif
00392 
00393          return ;
00394       }
00395    }
00396 
00397    /*--- the normal case [with data] ---*/
00398 
00399    ncall++ ;
00400 #if 0
00401 if(ncall%1000==0) fprintf(stderr,"spearman_fimfunc: ncall=%d\n",ncall);
00402 #endif
00403 
00404    /* copy input timeseries data (since we will mangle it) */
00405 
00406    if( nkeep < ntsc )
00407       for( ii=0 ; ii < nkeep ; ii++ ) tsc[ii] = ts[keep[ii]]; /* the hard way */
00408    else
00409       memcpy(tsc,ts,sizeof(float)*ntsc) ;                     /* the easy way */
00410 
00411    /* compute our single result */
00412 
00413    v    = (float *) val ;
00414    v[0] = spearman_rank_manycorr( nkeep , tsc , nref , refv , rr ) ;
00415 
00416    return ;
00417 }
00418 
00419 /*--------------------------------------------------------------------------
00420   A sample fimfunc to compute the quadrant correlation
00421 ----------------------------------------------------------------------------*/
00422 
00423 void quadrant_fimfunc( int n, float *ts, void *ud, int nb, void *val )
00424 {
00425    static float *tsc=NULL , *refc=NULL , *refv=NULL , **rr ;
00426    static int   *keep=NULL ;
00427    static int    ntsc , nref , nkeep , polort,ignore , ncall;
00428 
00429    int ii , kk ;
00430    float *v ;
00431 
00432    /*--- handle special cases ---*/
00433 
00434    if( ts == NULL ){
00435 
00436       /*--- the initialization call ---*/
00437 
00438       if( n > 0 ){
00439 
00440          FIMdata *fd = (FIMdata *) val ;
00441 
00442          polort = fd->polort ;  /* not used here */
00443          ignore = fd->ignore ;
00444          ncall  = 0 ;           /* how many times we have been called */
00445 
00446          /* make workspace for copy of ts data when it arrives */
00447 
00448          ntsc = n ;
00449          if( tsc != NULL ) free(tsc) ;
00450          tsc = (float *) malloc(sizeof(float)*ntsc) ;
00451 
00452          /* make copy of ref data */
00453 
00454          nref = fd->ref_ts->ny ;
00455          if( refc != NULL ) free(refc) ;
00456          if( refv != NULL ) free(refv) ;
00457          if( keep != NULL ) free(keep) ;
00458          if( rr   != NULL ) free(rr) ;
00459          refc = (float *) malloc(sizeof(float)*ntsc*nref) ;  /* copy of ref   */
00460          refv = (float *) malloc(sizeof(float)*nref) ;      /* rank variances */
00461          keep = (int *)   malloc(sizeof(int)*ntsc) ;       /* keeper indices  */
00462          rr   = (float **)malloc(sizeof(float *)*nref) ;  /* convenience ptrs */
00463 
00464          for( kk=0 ; kk < nref ; kk++ ){
00465             rr[kk] = refc+kk*ntsc ;                       /* compute ptr */
00466             memcpy( rr[kk] ,                              /* copy data  */
00467                     MRI_FLOAT_PTR(fd->ref_ts) + kk*fd->ref_ts->nx ,
00468                     sizeof(float)*ntsc  ) ;
00469          }
00470 
00471          /* mark those places we will keep (where ref is OK) */
00472 
00473          for( nkeep=0,ii=ignore ; ii < ntsc ; ii++ ){ /* for each time point */
00474             for( kk=0 ; kk < nref ; kk++ )            /* check each ref     */
00475                if( rr[kk][ii] >= WAY_BIG ) break ;
00476 
00477             if( kk == nref ) keep[nkeep++] = ii ;     /* mark if all are OK */
00478          }
00479 
00480          /* compress ref, eliminating non-keepers */
00481 
00482          if( nkeep < ntsc ){
00483             for( ii=0 ; ii < nkeep ; ii++ ){
00484                for( kk=0 ; kk < nref ; kk++ )
00485                   rr[kk][ii] = rr[kk][keep[ii]] ;  /* copy backwards */
00486             }
00487          }
00488 
00489          /* prepare each ref vector for rank processing */
00490 
00491          for( kk=0 ; kk < nref ; kk++ )
00492             refv[kk] = quadrant_corr_prepare( nkeep , rr[kk] ) ;
00493 
00494 #if 0
00495 fprintf(stderr,"quadrant_fimfunc: initialize ntsc=%d nkeep=%d nref=%d\n",
00496         ntsc,nkeep,nref);
00497 #endif
00498 
00499          return ;
00500 
00501       /*--- the ending call ---*/
00502 
00503       } else {
00504          THD_3dim_dataset *dset = (THD_3dim_dataset *) val ;
00505          int kb = -n ;
00506 
00507          free(tsc)  ; tsc  = NULL ;  /* free workspaces */
00508          free(refc) ; refc = NULL ;
00509          free(refv) ; refv = NULL ;
00510          free(keep) ; keep = NULL ;
00511          free(rr)   ; rr   = NULL ;
00512          rank_order_float(0,NULL) ;
00513 
00514          /* edit sub-brick statistical parameters and name */
00515 
00516          EDIT_BRICK_TO_FICO( dset , kb , nkeep , (nref==1)?1:2 , 1 ) ;
00517          EDIT_BRICK_LABEL( dset , kb , "Quadrant CC" ) ;
00518 
00519 #if 0
00520 fprintf(stderr,"quadrant_fimfunc: finalize with kb=%d\n",kb);
00521 #endif
00522 
00523          return ;
00524       }
00525    }
00526 
00527    /*--- the normal case [with data] ---*/
00528 
00529    ncall++ ;
00530 #if 0
00531 if(ncall%1000==0) fprintf(stderr,"quadrant_fimfunc: ncall=%d\n",ncall);
00532 #endif
00533 
00534    /* copy input timeseries data (since we will mangle it) */
00535 
00536    if( nkeep < ntsc )
00537       for( ii=0 ; ii < nkeep ; ii++ ) tsc[ii] = ts[keep[ii]]; /* the hard way */
00538    else
00539       memcpy(tsc,ts,sizeof(float)*ntsc) ;                     /* the easy way */
00540 
00541    /* compute our single result */
00542 
00543    v    = (float *) val ;
00544    v[0] = quadrant_manycorr( nkeep , tsc , nref , refv , rr ) ;
00545 
00546    return ;
00547 }