/*****************************************************************************
Major portions of this software are copyrighted by the Medical College
of Wisconsin, 1994-2000, and are released under the Gnu General Public
License, Version 2. See the file README.Copyright for details.
******************************************************************************/
/*---------------------------------------------------------------------------*/
/*
This file contains routines for implementing the fim+ functions.
File: fim+.c
Author: B. Douglas Ward
Date: 28 April 2000
Mod: Use output_type array in regression_analysis routine to avoid
some unnecessary calculations.
Date: 18 May 2000
*/
/*---------------------------------------------------------------------------*/
/*
Include software for linear regression analysis and sorting numbers.
*/
#include "RegAna.c"
#include "ranks.c"
/*---------------------------------------------------------------------------*/
/*
FIM output parameter definitions
*/
#define MAX_OUTPUT_TYPE 12
static char * OUTPUT_TYPE_name[MAX_OUTPUT_TYPE] =
{ "Fit Coef", "Best Index", "% Change", "% From Ave", "Baseline", "Average",
"Correlation", "% From Top", "Topline", "Sigma Resid",
"Spearman CC", "Quadrant CC" } ;
#define FIM_FitCoef (0)
#define FIM_BestIndex (1)
#define FIM_PrcntChange (2)
#define FIM_Baseline (4)
#define FIM_Correlation (6)
#define FIM_PrcntFromAve (3)
#define FIM_Average (5)
#define FIM_PrcntFromTop (7)
#define FIM_Topline (8)
#define FIM_SigmaResid (9)
#define FIM_SpearmanCC (10)
#define FIM_QuadrantCC (11)
/*---------------------------------------------------------------------------*/
#define USE_LEGENDRE
#ifdef USE_LEGENDRE
double legendre( double x , int m ) /* Legendre polynomials over [-1,1] */
{
if( m < 0 ) return 1.0 ; /* bad input */
switch( m ){ /*** P_m(x) for m=0..20 ***/
case 0: return 1.0 ;
case 1: return x ;
case 2: return (3.0*x*x-1.0)/2.0 ;
case 3: return (5.0*x*x-3.0)*x/2.0 ;
case 4: return ((35.0*x*x-30.0)*x*x+3.0)/8.0 ;
case 5: return ((63.0*x*x-70.0)*x*x+15.0)*x/8.0 ;
case 6: return (((231.0*x*x-315.0)*x*x+105.0)*x*x-5.0)/16.0 ;
case 7: return (((429.0*x*x-693.0)*x*x+315.0)*x*x-35.0)*x/16.0 ;
case 8: return ((((6435.0*x*x-12012.0)*x*x+6930.0)*x*x-1260.0)*x*x+35.0)/128.0;
/** 07 Feb 2005: this part generated by Maple, then hand massaged **/
case 9:
return (0.24609375e1 + (-0.3609375e2 + (0.140765625e3 + (-0.20109375e3
+ 0.949609375e2 * x * x) * x * x) * x * x) * x * x) * x;
case 10:
return -0.24609375e0 + (0.1353515625e2 + (-0.1173046875e3 +
(0.3519140625e3 + (-0.42732421875e3 + 0.18042578125e3 * x * x)
* x * x) * x * x) * x * x) * x * x;
case 11:
return (-0.270703125e1 + (0.5865234375e2 + (-0.3519140625e3 +
(0.8546484375e3 + (-0.90212890625e3 + 0.34444921875e3 * x * x)
* x * x) * x * x) * x * x) * x * x) * x;
case 12:
return 0.2255859375e0 + (-0.17595703125e2 + (0.2199462890625e3 +
(-0.99708984375e3 + (0.20297900390625e4 + (-0.1894470703125e4
+ 0.6601943359375e3 * x * x) * x * x) * x * x) * x * x) * x * x)
* x * x;
case 13:
return (0.29326171875e1 + (-0.87978515625e2 + (0.7478173828125e3 +
(-0.270638671875e4 + (0.47361767578125e4 + (-0.3961166015625e4
+ 0.12696044921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
* x * x) * x;
case 14:
return -0.20947265625e0 + (0.2199462890625e2 + (-0.37390869140625e3 +
(0.236808837890625e4 + (-0.710426513671875e4 +
(0.1089320654296875e5 + (-0.825242919921875e4 +
0.244852294921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
* x * x) * x * x;
case 15:
return (-0.314208984375e1 + (0.12463623046875e3 + (-0.142085302734375e4
+ (0.710426513671875e4 + (-0.1815534423828125e5 +
(0.2475728759765625e5 + (-0.1713966064453125e5 +
0.473381103515625e4 * x * x) * x * x) * x * x) * x * x)
* x * x) * x * x) * x * x) * x;
case 16:
return 0.196380615234375e0 + (-0.26707763671875e2 + (0.5920220947265625e3
+ (-0.4972985595703125e4 + (0.2042476226806641e5 +
(-0.4538836059570312e5 + (0.5570389709472656e5 +
(-0.3550358276367188e5 + 0.9171758880615234e4 * x * x) * x * x)
* x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
case 17:
return (0.3338470458984375e1 + (-0.169149169921875e3 +
(0.2486492797851562e4 + (-0.1633980981445312e5 +
(0.5673545074462891e5 + (-0.1114077941894531e6 +
(0.1242625396728516e6 + (-0.7337407104492188e5 +
0.1780400253295898e5 * x * x) * x * x) * x * x) * x * x)
* x * x) * x * x) * x * x) * x * x) * x;
case 18:
return -0.1854705810546875e0 + (0.3171546936035156e2 +
(-0.8880331420898438e3 + (0.9531555725097656e4 +
(-0.5106190567016602e5 + (0.153185717010498e6 +
(-0.2692355026245117e6 + (0.275152766418457e6 +
(-0.1513340215301514e6 + 0.3461889381408691e5 * x * x) * x * x)
* x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
case 19:
return (-0.3523941040039062e1 + (0.2220082855224609e3 +
(-0.4084952453613281e4 + (0.3404127044677734e5 +
(-0.153185717010498e6 + (0.4038532539367676e6 +
(-0.6420231216430664e6 + (0.6053360861206055e6 +
(-0.3115700443267822e6 + 0.6741574058532715e5 * x * x) * x * x)
* x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x;
case 20:
return 0.1761970520019531e0 + (-0.3700138092041016e2 +
(0.127654764175415e4 + (-0.1702063522338867e5 +
(0.1148892877578735e6 + (-0.4442385793304443e6 +
(0.1043287572669983e7 + (-0.1513340215301514e7 +
(0.1324172688388824e7 + (-0.6404495355606079e6 +
0.1314606941413879e6 * x * x) * x * x) * x * x) * x * x) * x * x)
* x * x) * x * x) * x * x) * x * x) * x * x;
}
#if 0
/* order out of range: return Chebyshev instead (it's easy) */
if( x >= 1.0 ) x = 0.0 ;
else if ( x <= -1.0 ) x = 3.14159265358979323846 ;
else x = acos(x) ;
return cos(m*x) ;
#else
/** if here, m > 20 ==> use recurrence relation **/
{ double pk, pkm1, pkm2 ; int k ;
pkm2 = legendre( x , 19 ) ;
pkm1 = legendre( x , 20 ) ;
for( k=21 ; k <= m ; k++ , pkm2=pkm1 , pkm1=pk )
pk = ((2.0*k-1.0)*x*pkm1 - (k-1.0)*pkm2)/k ;
return pk ;
}
#endif
}
#endif /* USE_LEGENDRE */
/*---------------------------------------------------------------------------*/
/*
Initialize independent variable X matrix
*/
int init_indep_var_matrix
(
int q, /* number of parameters in the baseline model */
int p, /* number of parameters in the baseline model
plus number of ideals */
int NFirst, /* index of first image to use in the analysis */
int N, /* total number of images used in the analysis */
int polort, /* degree of polynomial ort */
int num_ort_files, /* number of ort time series files */
int num_ideal_files, /* number of ideal time series files */
MRI_IMAGE ** ort_array, /* ort time series arrays */
int ** ort_list, /* list of ort time series */
MRI_IMAGE ** ideal_array, /* ideal time series arrays */
int ** ideal_list, /* list of ideal time series */
float * x_bot, /* minimum of stimulus time series */
float * x_ave, /* average of stimulus time series */
float * x_top, /* maximum of stimulus time series */
int * good_list, /* list of good time points */
matrix * x /* independent variable matrix */
)
{
const int BIG_NUMBER = 33333;
int i; /* time index */
int m; /* X matrix column index */
int n; /* X matrix row index */
int is; /* input ideal index */
float * far = NULL;
int nx, ny, iq, nq;
int Ngood;
matrix xgood;
#ifdef USE_LEGENDRE
double nfac,nsub ;
#endif
/*----- Initialize X matrix -----*/
matrix_create (N, p, x);
matrix_initialize (&xgood);
#ifdef USE_LEGENDRE
nsub = 0.5*(N-1) ;
nfac = 1.0/nsub ;
#endif
/*----- Set up columns of X matrix corresponding to polynomial orts -----*/
for (m = 0; m <= polort; m++)
for (n = 0; n < N; n++)
{
#ifndef USE_LEGENDRE /** the old polort way: t^m **/
i = n + NFirst;
(*x).elts[n][m] = pow ((double)i, (double)m);
#else /** the new polort way: Legendre P_m(t) - 29 Mar 2005 **/
(*x).elts[n][m] = legendre( nfac*(n-nsub) , m ) ;
#endif
}
/*----- Set up columns of X matrix corresponding to ort time series -----*/
for (is = 0; is < num_ort_files; is++)
{
far = MRI_FLOAT_PTR (ort_array[is]);
nx = ort_array[is]->nx;
ny = ort_array[is]->ny;
if (ort_list[is] == NULL)
for (iq = 0; iq < ny; iq++)
{
for (n = 0; n < N; n++)
{
i = n + NFirst;
(*x).elts[n][m] = *(far + iq*nx + i);
}
m++;
}
else
{
nq = ort_list[is][0];
for (iq = 1; iq <= nq; iq++)
{
for (n = 0; n < N; n++)
{
i = n + NFirst;
(*x).elts[n][m] = *(far + ort_list[is][iq]*nx + i);
}
m++;
}
}
}
/*----- Set up columns of X matrix corresponding to ideal time series -----*/
for (is = 0; is < num_ideal_files; is++)
{
far = MRI_FLOAT_PTR (ideal_array[is]);
nx = ideal_array[is]->nx;
ny = ideal_array[is]->ny;
if (ideal_list[is] == NULL)
for (iq = 0; iq < ny; iq++)
{
for (n = 0; n < N; n++)
{
i = n + NFirst;
(*x).elts[n][m] = *(far + iq*nx + i);
}
m++;
}
else
{
nq = ideal_list[is][0];
for (iq = 1; iq <= nq; iq++)
{
for (n = 0; n < N; n++)
{
i = n + NFirst;
(*x).elts[n][m] = *(far + ideal_list[is][iq]*nx + i);
}
m++;
}
}
}
/*----- Remove row if ideal contains value over 33333 -----*/
Ngood = N;
m = 0;
for (n = 0; n < N; n++)
{
for (is = q; is < p; is++)
{
if ((*x).elts[n][is] >= BIG_NUMBER) break;
}
if (is < p)
{
Ngood--;
}
else
{
good_list[m] = n;
m++;
}
}
matrix_extract_rows ((*x), Ngood, good_list, &xgood);
matrix_equate (xgood, x);
/*----- Find min, max, and ave for each column of the X matrix -----*/
for (is = 0; is < p; is++)
{
x_bot[is] = x_top[is] = (*x).elts[0][is];
x_ave[is] = 0.0;
for (n = 0; n < Ngood; n++)
{
if ((*x).elts[n][is] < x_bot[is]) x_bot[is] = (*x).elts[n][is];
if ((*x).elts[n][is] > x_top[is]) x_top[is] = (*x).elts[n][is];
x_ave[is] += (*x).elts[n][is] / Ngood;
}
}
matrix_destroy (&xgood);
return (Ngood);
}
/*---------------------------------------------------------------------------*/
/*
Initialization for the regression analysis.
*/
int init_regression_analysis
(
int q, /* number of parameters in the baseline model */
int p, /* number of parameters in the baseline model
plus number of ideals */
int N, /* total number of images used in the analysis */
int num_idealts, /* number of ideal time series */
matrix xdata, /* independent variable matrix */
matrix * x_base, /* extracted X matrix for baseline model */
matrix * xtxinvxt_base, /* matrix: (1/(X'X))X' for baseline model */
matrix * x_ideal, /* extracted X matrices for ideal models */
matrix * xtxinvxt_ideal, /* matrix: (1/(X'X))X' for ideal models */
float ** rarray /* ranked arrays of ideal time series */
)
{
int * plist = NULL; /* list of model parameters */
int ip, it; /* parameter indices */
int is, js; /* ideal indices */
int jm; /* lag index */
int ok; /* flag for successful matrix calculation */
matrix xtxinv_temp; /* intermediate results */
vector ideal; /* ideal vector */
vector coef_temp; /* intermediate results */
vector xres; /* vector of residuals */
float sse_base; /* baseline model error sum of squares */
/*----- Initialize matrix -----*/
matrix_initialize (&xtxinv_temp);
vector_initialize (&ideal);
vector_initialize (&coef_temp);
vector_initialize (&xres);
/*----- Allocate memory -----*/
plist = (int *) malloc (sizeof(int) * p); MTEST (plist);
/*----- Initialize matrices for the baseline model -----*/
for (ip = 0; ip < q; ip++)
plist[ip] = ip;
ok = calc_matrices (xdata, q, plist, x_base, &xtxinv_temp, xtxinvxt_base);
if (!ok) { matrix_destroy (&xtxinv_temp); return (0); };
/*----- Initialize matrices for ideal functions -----*/
for (is = 0; is < num_idealts; is++)
{
for (ip = 0; ip < q; ip++)
{
plist[ip] = ip;
}
plist[q] = q + is;
ok = calc_matrices (xdata, q+1, plist,
&(x_ideal[is]), &xtxinv_temp, &(xtxinvxt_ideal[is]));
if (!ok) { matrix_destroy (&xtxinv_temp); return (0); };
}
/*----- Set up the ranked array for each ideal -----*/
for (is = 0; is < num_idealts; is++)
{
/*----- Convert column of matrix to vector -----*/
column_to_vector (xdata, q+is, &ideal);
/*----- Calculate regression coefficients for baseline model -----*/
calc_coef (*xtxinvxt_base, ideal, &coef_temp);
/*----- Calculate the error sum of squares for the baseline model -----*/
sse_base = calc_resids (*x_base, coef_temp, ideal, &xres);
/*----- Form rank array from residual array -----*/
rarray[is] = rank_darray (N, xres.elts);
}
/*----- Destroy matrix -----*/
matrix_destroy (&xtxinv_temp);
vector_destroy (&ideal);
vector_destroy (&coef_temp);
vector_destroy (&xres);
/*----- Deallocate memory -----*/
free (plist); plist = NULL;
return (1);
}
/*---------------------------------------------------------------------------*/
/*
Calculate a correlation coefficient.
*/
float calc_CC
(
int N,
float * x,
float * y
)
{
const float EPSILON = 1.0e-10; /* protect against divide by zero */
float cc;
int i;
float csum, xsum, ysum, prod;
csum = xsum = ysum = 0.0;
for (i = 0; i < N; i++)
{
csum += x[i] * y[i];
xsum += x[i] * x[i];
ysum += y[i] * y[i];
}
prod = xsum * ysum;
if (prod < EPSILON)
cc = 0.0;
else
cc = csum / sqrt(prod);
return (cc);
}
/*---------------------------------------------------------------------------*/
/*
Calculate the Spearman correlation coefficient.
*/
float calc_SpearmanCC
(
int N,
float * r,
float * s
)
{
float cc;
float rbar;
float * dr = NULL;
float * ds = NULL;
int i;
dr = (float *) malloc (sizeof(float) * N);
ds = (float *) malloc (sizeof(float) * N);
rbar = (N+1.0) / 2.0;
for (i = 0; i < N; i++)
{
dr[i] = r[i] - rbar;
ds[i] = s[i] - rbar;
}
cc = calc_CC(N, dr, ds);
free (dr); dr = NULL;
free (ds); ds = NULL;
return (cc);
}
/*---------------------------------------------------------------------------*/
/*
Calculate the sign function.
*/
float sign (float x)
{
if (x > 0.0) return (1.0);
if (x < 0.0) return (-1.0);
return (0.0);
}
/*---------------------------------------------------------------------------*/
/*
Calculate the Quadrant correlation coefficient.
*/
float calc_QuadrantCC
(
int N,
float * r,
float * s
)
{
float cc;
float rbar;
float * dr = NULL;
float * ds = NULL;
int i;
dr = (float *) malloc (sizeof(float) * N);
ds = (float *) malloc (sizeof(float) * N);
rbar = (N+1.0) / 2.0;
for (i = 0; i < N; i++)
{
dr[i] = sign(r[i] - rbar);
ds[i] = sign(s[i] - rbar);
}
cc = calc_CC(N, dr, ds);
free (dr); dr = NULL;
free (ds); ds = NULL;
return (cc);
}
/*---------------------------------------------------------------------------*/
/*
Calculate percentage change relative to baseline.
*/
float percent_change (float num, float den)
{
const float EPSILON = 1.0e-10; /* guard against divide by zero */
const float MAX_PERCENT = 1000.0; /* limit maximum percent change */
float PrcntChange;
if (fabs(den) < EPSILON)
PrcntChange = sign(num) * MAX_PERCENT;
else
PrcntChange = 100.0 * num / den;
if (PrcntChange > MAX_PERCENT) PrcntChange = MAX_PERCENT;
if (PrcntChange < -MAX_PERCENT) PrcntChange = -MAX_PERCENT;
return (PrcntChange);
}
/*---------------------------------------------------------------------------*/
/*
Calculate results for this voxel.
*/
void regression_analysis
(
int N, /* number of usable data points */
int q, /* number of parameters in baseline model */
int num_idealts, /* number of ideal time series */
matrix x_base, /* extracted X matrix for baseline model */
matrix xtxinvxt_base, /* matrix: (1/(X'X))X' for baseline model */
matrix * x_ideal, /* extracted X matrices for ideal models */
matrix * xtxinvxt_ideal, /* matrix: (1/(X'X))X' for ideal models */
vector y, /* vector of measured data */
float * x_bot, /* minimum of stimulus time series */
float * x_ave, /* average of stimulus time series */
float * x_top, /* maximum of stimulus time series */
float ** rarray, /* ranked arrays of ideal time series */
int * output_type, /* list of operator requested outputs */
float * FimParams /* output fim parameters */
)
{
const float EPSILON = 1.0e-05; /* protection against divide by zero */
int is; /* input ideal index */
float sse_base; /* error sum of squares, baseline model */
float sse_ideal; /* error sum of squares, ideal model */
vector coef_temp; /* intermediate results */
vector coef_best; /* best results */
vector yres; /* vector of residuals */
float rtemp, rbest; /* best correlation coefficient */
float stemp, sbest; /* best Spearman correlation coefficient */
float qtemp, qbest; /* best quadrant correlation coefficient */
float mse; /* mean square error (sample variance) */
float * sarray = NULL; /* ranked array of measurement data */
/* fim output parameters */
float FitCoef = 0.0;
int BestIndex = 0;
float PrcntChange = 0.0;
float Baseline = 0.0;
float Correlation = 0.0;
float PrcntFromAve = 0.0;
float Average = 0.0;
float PrcntFromTop = 0.0;
float Topline = 0.0;
float SigmaResid = 0.0;
float SpearmanCC = 0.0;
float QuadrantCC = 0.0;
/*----- Initialization -----*/
vector_initialize (&coef_temp);
vector_initialize (&coef_best);
vector_initialize (&yres);
/*----- Calculate regression coefficients for baseline model -----*/
calc_coef (xtxinvxt_base, y, &coef_temp);
/*----- Calculate the error sum of squares for the baseline model -----*/
sse_base = calc_resids (x_base, coef_temp, y, &yres);
/*----- Form rank array from y array -----*/
if (output_type[FIM_SpearmanCC] || output_type[FIM_QuadrantCC])
sarray = rank_darray (N, yres.elts);
/*----- Determine the best ideal reference for this voxel -----*/
rbest = 0.0; sbest = 0.0; qbest = 0.0; BestIndex = -1;
for (is = 0; is < num_idealts; is++)
{
/*----- Calculate regression coefficients for ideal model -----*/
calc_coef (xtxinvxt_ideal[is], y, &coef_temp);
/*----- Calculate the error sum of squares for the ideal model -----*/
sse_ideal = calc_sse (x_ideal[is], coef_temp, y);
/*----- Calculate partial R^2 for this ideal -----*/
rtemp = calc_rsqr (sse_ideal, sse_base);
if (rtemp >= rbest)
{
rbest = rtemp;
BestIndex = is;
vector_equate (coef_temp, &coef_best);
if (num_idealts == 1)
mse = sse_ideal / (N-q-1);
else
mse = sse_ideal / (N-q-2);
}
/*----- Calculate the Spearman rank correlation coefficient -----*/
if (output_type[FIM_SpearmanCC])
{
stemp = calc_SpearmanCC (N, rarray[is], sarray);
if (fabs(stemp) > fabs(sbest)) sbest = stemp;
}
/*----- Calculate the Quadrant correlation coefficient -----*/
if (output_type[FIM_QuadrantCC])
{
qtemp = calc_QuadrantCC (N, rarray[is], sarray);
if (fabs(qtemp) > fabs(qbest)) qbest = qtemp;
}
}
if ((0 <= BestIndex) && (BestIndex < num_idealts))
{
float Top, Ave, Base, Center;
int ip;
Top = x_top[q+BestIndex];
Ave = x_ave[q+BestIndex];
Base = x_bot[q+BestIndex];
FitCoef = coef_best.elts[q];
Correlation = sqrt(rbest);
if (FitCoef < 0.0) Correlation = -Correlation;
Center = 0.0;
for (ip = 0; ip < q; ip++)
Center += coef_best.elts[ip] * x_ave[ip];
Baseline = Center + FitCoef*Base;
Average = Center + FitCoef*Ave;
Topline = Center + FitCoef*Top;
PrcntChange = percent_change (FitCoef * (Top-Base), Baseline);
PrcntFromAve = percent_change (FitCoef * (Top-Base), Average);
PrcntFromTop = percent_change (FitCoef * (Top-Base), Topline);
SigmaResid = sqrt(mse);
SpearmanCC = sbest;
QuadrantCC = qbest;
}
/*----- Save output parameters -----*/
FimParams[FIM_FitCoef] = FitCoef;
FimParams[FIM_BestIndex] = BestIndex + 1.0;
FimParams[FIM_PrcntChange] = PrcntChange;
FimParams[FIM_Baseline] = Baseline;
FimParams[FIM_Correlation] = Correlation;
FimParams[FIM_PrcntFromAve] = PrcntFromAve;
FimParams[FIM_Average] = Average;
FimParams[FIM_PrcntFromTop] = PrcntFromTop;
FimParams[FIM_Topline] = Topline;
FimParams[FIM_SigmaResid] = SigmaResid;
FimParams[FIM_SpearmanCC] = SpearmanCC;
FimParams[FIM_QuadrantCC] = QuadrantCC;
/*----- Dispose of vectors -----*/
vector_destroy (&coef_temp);
vector_destroy (&coef_best);
vector_destroy (&yres);
/*----- Deallocate memory -----*/
if (sarray != NULL)
{ free (sarray); sarray = NULL; }
}
/*---------------------------------------------------------------------------*/
static char lbuf[4096]; /* character string containing statistical summary */
static char sbuf[256];
void report_results
(
int * output_type, /* output type options */
float * FimParams, /* output fim parameters */
char ** label /* statistical summary for ouput display */
)
{
int ip; /* parameter index */
if( label != NULL ){ /* assemble this 1 line at a time from sbuf */
lbuf[0] = '\0' ; /* make this a 0 length string to start */
/** for each reference, make a string into sbuf **/
for (ip = 0; ip < MAX_OUTPUT_TYPE; ip++)
if (output_type[ip])
{
sprintf (sbuf, "%12s = %10.4f \n",
OUTPUT_TYPE_name[ip], FimParams[ip]);
strcat (lbuf, sbuf);
}
*label = lbuf ; /* send address of lbuf back in what label points to */
}
}
/*---------------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1