00001 #include "mrilib.h"
00002 
00003 /*==============================================================================*/
00004 /*========== The following functions are lifted from afni_fimfunc.c ============*/
00005 /*------------------------------------------------------------------------------
00006    Rank-order a float array, with ties getting the average rank.
00007    The output overwrites the input.
00008 --------------------------------------------------------------------------------*/
00009 
00010 static void rank_order_float( int n , float * a )
00011 {
00012    register int ii , ns , n1 , ib ;
00013    static int    nb = 0 ;
00014    static int *   b = NULL ;  /* workspaces */
00015    static float * c = NULL ;
00016    float cs ;
00017 
00018    /*- handle special cases -*/
00019 
00020    if( a == NULL ){
00021       if( b != NULL ){ free(b); free(c); b=NULL ; c=NULL; nb=0; }  /* free workspaces */
00022       return ;
00023    }
00024 
00025    if( n < 1 ) return ;                     /* meaningless input */
00026    if( n == 1 ){ a[0] = 0.0 ; return ; }    /* only one point!? */
00027 
00028    /*- make workspaces, if needed -*/
00029 
00030    if( nb < n ){
00031       if( b != NULL ){ free(b); free(c); }
00032       b  = (int   *) malloc(sizeof(int  )*n) ;
00033       c  = (float *) malloc(sizeof(float)*n) ;
00034       nb = n ;
00035    }
00036 
00037    for( ii=0 ; ii < n ; ii++ ) c[ii] = b[ii] = ii ;
00038 
00039    /*- sort input, carrying b along -*/
00040 
00041    qsort_floatint( n , a , b ) ;  /* see cs_sort_fi.c */
00042 
00043    /* compute ranks into c[] */
00044 
00045    n1 = n-1 ;
00046    for( ii=0 ; ii < n1 ; ii++ ){
00047       if( a[ii] == a[ii+1] ){                  /* handle ties */
00048          cs = 2*ii+1 ; ns = 2 ; ib=ii ; ii++ ;
00049          while( ii < n1 && a[ii] == a[ii+1] ){ ii++ ; ns++ ; cs += ii ; }
00050          for( cs/=ns ; ib <= ii ; ib++ ) c[ib] = cs ;
00051       }
00052    }
00053 
00054    for( ii=0 ; ii < n ; ii++ ) a[b[ii]] = c[ii] ;
00055 
00056    return ;
00057 }
00058 
00059 /*---------------------------------------------------------------------------
00060    Rank orders a[], subtracts the mean rank, and returns the sum-of-squares
00061 -----------------------------------------------------------------------------*/
00062 
00063 static float spearman_rank_prepare( int n , float * a )
00064 {
00065    register int ii ;
00066    register float rb , rs ;
00067 
00068    rank_order_float( n , a ) ;
00069 
00070    rb = 0.5*(n-1) ; rs=0.0 ;
00071    for( ii=0 ; ii < n ; ii++ ){
00072       a[ii] -= rb ;
00073       rs    += a[ii]*a[ii] ;
00074    }
00075 
00076    return rs ;
00077 }
00078 
00079 /*------------------------------------------------------------------------------*/
00080 
00081 static float quadrant_corr_prepare( int n , float * a )
00082 {
00083    register int ii ;
00084    register float rb , rs ;
00085 
00086    rank_order_float( n , a ) ;
00087 
00088    rb = 0.5*(n-1) ; rs=0.0 ;
00089    for( ii=0 ; ii < n ; ii++ ){
00090       a[ii] = (a[ii] > rb) ? 1.0
00091                            : (a[ii] < rb) ? -1.0 : 0.0 ;
00092       rs    += a[ii]*a[ii] ;
00093    }
00094 
00095    return rs ;
00096 }
00097 
00098 /*-----------------------------------------------------------------------------
00099     To correlate x[] with r[], do
00100       rv = spearman_rank_prepare(n,r) ;
00101     then
00102       corr = spearman_rank_corr(n,x,rv,r) ;
00103     Note that these 2 routines are destructive (r and x are replaced by ranks)
00104 -------------------------------------------------------------------------------*/
00105 
00106 static float spearman_rank_corr( int n , float * x , float rv , float * r )
00107 {
00108    register int ii ;
00109    register float ss ; float xv ;
00110 
00111    xv = spearman_rank_prepare( n , x ) ; if( xv <= 0.0 ) return 0.0 ;
00112 
00113    for( ii=0,ss=0.0 ; ii < n ; ii++ ) ss += x[ii] * r[ii] ;
00114 
00115    return ( ss/sqrt(rv*xv) ) ;
00116 }
00117 
00118 /*------------------------------------------------------------------------------*/
00119 
00120 static float quadrant_corr( int n , float * x , float rv , float * r )
00121 {
00122    register int ii ;
00123    register float ss ; float xv ;
00124 
00125    xv = quadrant_corr_prepare( n , x ) ; if( xv <= 0.0 ) return 0.0 ;
00126 
00127    for( ii=0,ss=0.0 ; ii < n ; ii++ ) ss += x[ii] * r[ii] ;
00128 
00129    return ( ss/sqrt(rv*xv) ) ;
00130 }
00131 
00132 /*=============================================================================
00133   Compute correlations, destructively
00134 ===============================================================================*/
00135 
00136 float THD_spearman_corr( int n , float *x , float *y )
00137 {
00138    float xv = spearman_rank_prepare(n,x) ;
00139    if( xv <= 0.0 ) return 0.0 ;
00140    return spearman_rank_corr( n,y,xv,x ) ;
00141 }
00142 
00143 float THD_quadrant_corr( int n , float *x , float *y )
00144 {
00145    float xv = quadrant_corr_prepare(n,x) ;
00146    if( xv <= 0.0 ) return 0.0 ;
00147    return quadrant_corr( n,y,xv,x ) ;
00148 }
00149 
00150 float THD_pearson_corr( int n, float *x , float *y )
00151 {
00152    float xv=0 , yv=0 , xy=0 ;
00153    int ii ;
00154 
00155    for( ii=0 ; ii < n ; ii++ ){
00156       xv += x[ii]*x[ii] ; yv += y[ii]*y[ii] ; xy += x[ii]*y[ii] ;
00157    }
00158 
00159    if( xv <= 0.0 || yv <= 0.0 ) return 0.0 ;
00160    return xy/sqrt(xv*yv) ;
00161 }