/***************************************************************************** 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 */ } } /*---------------------------------------------------------------------------*/