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 /*---------------------------------------------------------------------------*/
00008 /*
00009   This file contains routines for implementing the fim+ functions.
00010 
00011   File:    fim+.c
00012   Author:  B. Douglas Ward
00013   Date:    28 April 2000
00014 
00015 
00016   Mod:     Use output_type array in regression_analysis routine to avoid
00017            some unnecessary calculations.
00018   Date:    18 May 2000
00019 */
00020 
00021 /*---------------------------------------------------------------------------*/
00022 /*
00023   Include software for linear regression analysis and sorting numbers.
00024 */
00025 
00026 #include "RegAna.c"
00027 #include "ranks.c"
00028 
00029 /*---------------------------------------------------------------------------*/
00030 /*
00031   FIM output parameter definitions 
00032 */
00033 
00034 #define MAX_OUTPUT_TYPE 12
00035 
00036 static char * OUTPUT_TYPE_name[MAX_OUTPUT_TYPE] =
00037   { "Fit Coef", "Best Index", "% Change", "% From Ave", "Baseline", "Average", 
00038     "Correlation", "% From Top", "Topline", "Sigma Resid",
00039     "Spearman CC", "Quadrant CC" } ;
00040 
00041 
00042 #define FIM_FitCoef      (0)
00043 #define FIM_BestIndex    (1)
00044 #define FIM_PrcntChange  (2)
00045 #define FIM_Baseline     (4)
00046 #define FIM_Correlation  (6)
00047 #define FIM_PrcntFromAve (3)
00048 #define FIM_Average      (5)
00049 #define FIM_PrcntFromTop (7)
00050 #define FIM_Topline      (8)
00051 #define FIM_SigmaResid   (9)
00052 #define FIM_SpearmanCC   (10)
00053 #define FIM_QuadrantCC   (11)
00054 
00055 /*---------------------------------------------------------------------------*/
00056 
00057 #define USE_LEGENDRE
00058 
00059 #ifdef USE_LEGENDRE
00060 double legendre( double x , int m )   /* Legendre polynomials over [-1,1] */
00061 {
00062    if( m < 0 ) return 1.0 ;    /* bad input */
00063 
00064    switch( m ){                /*** P_m(x) for m=0..20 ***/
00065     case 0: return 1.0 ;
00066     case 1: return x ;
00067     case 2: return (3.0*x*x-1.0)/2.0 ;
00068     case 3: return (5.0*x*x-3.0)*x/2.0 ;
00069     case 4: return ((35.0*x*x-30.0)*x*x+3.0)/8.0 ;
00070     case 5: return ((63.0*x*x-70.0)*x*x+15.0)*x/8.0 ;
00071     case 6: return (((231.0*x*x-315.0)*x*x+105.0)*x*x-5.0)/16.0 ;
00072     case 7: return (((429.0*x*x-693.0)*x*x+315.0)*x*x-35.0)*x/16.0 ;
00073     case 8: return ((((6435.0*x*x-12012.0)*x*x+6930.0)*x*x-1260.0)*x*x+35.0)/128.0;
00074 
00075            /** 07 Feb 2005: this part generated by Maple, then hand massaged **/
00076 
00077     case 9:
00078       return (0.24609375e1 + (-0.3609375e2 + (0.140765625e3 + (-0.20109375e3
00079               + 0.949609375e2 * x * x) * x * x) * x * x) * x * x) * x;
00080 
00081     case 10:
00082       return -0.24609375e0 + (0.1353515625e2 + (-0.1173046875e3 +
00083               (0.3519140625e3 + (-0.42732421875e3 + 0.18042578125e3 * x * x)
00084              * x * x) * x * x) * x * x) * x * x;
00085 
00086     case 11:
00087       return (-0.270703125e1 + (0.5865234375e2 + (-0.3519140625e3 +
00088              (0.8546484375e3 + (-0.90212890625e3 + 0.34444921875e3 * x * x)
00089              * x * x) * x * x) * x * x) * x * x) * x;
00090 
00091     case 12:
00092       return 0.2255859375e0 + (-0.17595703125e2 + (0.2199462890625e3 +
00093              (-0.99708984375e3 + (0.20297900390625e4 + (-0.1894470703125e4
00094              + 0.6601943359375e3 * x * x) * x * x) * x * x) * x * x) * x * x)
00095              * x * x;
00096 
00097     case 13:
00098       return (0.29326171875e1 + (-0.87978515625e2 + (0.7478173828125e3 +
00099              (-0.270638671875e4 + (0.47361767578125e4 + (-0.3961166015625e4
00100              + 0.12696044921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
00101             * x * x) * x;
00102 
00103     case 14:
00104       return -0.20947265625e0 + (0.2199462890625e2 + (-0.37390869140625e3 +
00105              (0.236808837890625e4 + (-0.710426513671875e4 +
00106              (0.1089320654296875e5 + (-0.825242919921875e4 +
00107             0.244852294921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
00108            * x * x) * x * x;
00109 
00110     case 15:
00111       return (-0.314208984375e1 + (0.12463623046875e3 + (-0.142085302734375e4
00112             + (0.710426513671875e4 + (-0.1815534423828125e5 +
00113               (0.2475728759765625e5 + (-0.1713966064453125e5 +
00114                0.473381103515625e4 * x * x) * x * x) * x * x) * x * x)
00115              * x * x) * x * x) * x * x) * x;
00116 
00117     case 16:
00118       return 0.196380615234375e0 + (-0.26707763671875e2 + (0.5920220947265625e3
00119             + (-0.4972985595703125e4 + (0.2042476226806641e5 +
00120               (-0.4538836059570312e5 + (0.5570389709472656e5 +
00121                (-0.3550358276367188e5 + 0.9171758880615234e4 * x * x) * x * x)
00122             * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
00123 
00124     case 17:
00125       return (0.3338470458984375e1 + (-0.169149169921875e3 +
00126              (0.2486492797851562e4 + (-0.1633980981445312e5 +
00127              (0.5673545074462891e5 + (-0.1114077941894531e6 +
00128              (0.1242625396728516e6 + (-0.7337407104492188e5 +
00129               0.1780400253295898e5 * x * x) * x * x) * x * x) * x * x)
00130            * x * x) * x * x) * x * x) * x * x) * x;
00131 
00132     case 18:
00133       return -0.1854705810546875e0 + (0.3171546936035156e2 +
00134             (-0.8880331420898438e3 + (0.9531555725097656e4 +
00135             (-0.5106190567016602e5 + (0.153185717010498e6 +
00136             (-0.2692355026245117e6 + (0.275152766418457e6 +
00137             (-0.1513340215301514e6 + 0.3461889381408691e5 * x * x) * x * x)
00138            * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
00139 
00140     case 19:
00141       return (-0.3523941040039062e1 + (0.2220082855224609e3 +
00142              (-0.4084952453613281e4 + (0.3404127044677734e5 +
00143              (-0.153185717010498e6 + (0.4038532539367676e6 +
00144              (-0.6420231216430664e6 + (0.6053360861206055e6 +
00145              (-0.3115700443267822e6 + 0.6741574058532715e5 * x * x) * x * x)
00146           * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x;
00147 
00148     case 20:
00149       return 0.1761970520019531e0 + (-0.3700138092041016e2 +
00150             (0.127654764175415e4 + (-0.1702063522338867e5 +
00151             (0.1148892877578735e6 + (-0.4442385793304443e6 +
00152             (0.1043287572669983e7 + (-0.1513340215301514e7 +
00153             (0.1324172688388824e7 + (-0.6404495355606079e6 +
00154              0.1314606941413879e6 * x * x) * x * x) * x * x) * x * x) * x * x)
00155             * x * x) * x * x) * x * x) * x * x) * x * x;
00156    }
00157 
00158 #if 0
00159    /* order out of range: return Chebyshev instead (it's easy) */
00160 
00161         if(  x >=  1.0 ) x = 0.0 ;
00162    else if ( x <= -1.0 ) x = 3.14159265358979323846 ;
00163    else                  x = acos(x) ;
00164    return cos(m*x) ;
00165 #else
00166    /** if here, m > 20 ==> use recurrence relation **/
00167 
00168    { double pk, pkm1, pkm2 ; int k ;
00169      pkm2 = legendre( x , 19 ) ;
00170      pkm1 = legendre( x , 20 ) ;
00171      for( k=21 ; k <= m ; k++ , pkm2=pkm1 , pkm1=pk )
00172        pk = ((2.0*k-1.0)*x*pkm1 - (k-1.0)*pkm2)/k ;
00173      return pk ;
00174    }
00175 #endif
00176 }
00177 #endif /* USE_LEGENDRE */
00178 
00179 /*---------------------------------------------------------------------------*/
00180 /* 
00181    Initialize independent variable X matrix 
00182 */
00183 
00184 int init_indep_var_matrix 
00185 (
00186   int q,                      /* number of parameters in the baseline model */
00187   int p,                      /* number of parameters in the baseline model 
00188                                  plus number of ideals */
00189   int NFirst,                 /* index of first image to use in the analysis */
00190   int N,                      /* total number of images used in the analysis */
00191   int polort,                 /* degree of polynomial ort */
00192   int num_ort_files,          /* number of ort time series files */
00193   int num_ideal_files,        /* number of ideal time series files */
00194   MRI_IMAGE ** ort_array,     /* ort time series arrays */
00195   int ** ort_list,            /* list of ort time series */
00196   MRI_IMAGE ** ideal_array,   /* ideal time series arrays */
00197   int ** ideal_list,          /* list of ideal time series */
00198   float * x_bot,              /* minimum of stimulus time series */
00199   float * x_ave,              /* average of stimulus time series */
00200   float * x_top,              /* maximum of stimulus time series */
00201   int * good_list,            /* list of good time points */
00202   matrix * x                  /* independent variable matrix */
00203 )
00204 
00205 {
00206   const int BIG_NUMBER = 33333;
00207   int i;                    /* time index */
00208   int m;                    /* X matrix column index */
00209   int n;                    /* X matrix row index */
00210   int is;                   /* input ideal index */
00211   float * far = NULL;
00212   int nx, ny, iq, nq;
00213   int Ngood;
00214   matrix xgood;
00215 #ifdef USE_LEGENDRE
00216   double nfac,nsub ;
00217 #endif
00218 
00219 
00220   /*----- Initialize X matrix -----*/
00221   matrix_create (N, p, x);
00222   matrix_initialize (&xgood);
00223 
00224 #ifdef USE_LEGENDRE
00225   nsub = 0.5*(N-1) ;
00226   nfac = 1.0/nsub  ;
00227 #endif
00228 
00229   /*----- Set up columns of X matrix corresponding to polynomial orts -----*/
00230   for (m = 0;  m <= polort;  m++)
00231     for (n = 0;  n < N;  n++)
00232       {
00233 #ifndef USE_LEGENDRE   /** the old polort way: t^m **/
00234              i = n + NFirst;
00235              (*x).elts[n][m] = pow ((double)i, (double)m);
00236 
00237 #else                  /** the new polort way: Legendre P_m(t) - 29 Mar 2005 **/
00238 
00239              (*x).elts[n][m] = legendre( nfac*(n-nsub) , m ) ;
00240 #endif
00241       }
00242 
00243 
00244   /*----- Set up columns of X matrix corresponding to ort time series -----*/
00245   for (is = 0;  is < num_ort_files;  is++)
00246     {
00247       far = MRI_FLOAT_PTR (ort_array[is]);
00248       nx = ort_array[is]->nx;
00249       ny = ort_array[is]->ny;
00250 
00251       if (ort_list[is] == NULL)
00252         for (iq = 0;  iq < ny;  iq++)
00253           {
00254             for (n = 0;  n < N;  n++)
00255               {
00256                 i = n + NFirst;
00257                 (*x).elts[n][m] = *(far + iq*nx + i);
00258               }
00259             m++;
00260           }
00261       else
00262         {
00263           nq = ort_list[is][0];
00264           for (iq = 1;  iq <= nq;  iq++)
00265             {
00266               for (n = 0;  n < N;  n++)
00267                 {
00268                   i = n + NFirst;
00269                   (*x).elts[n][m] = *(far + ort_list[is][iq]*nx + i);
00270                 }
00271               m++;
00272             }
00273         }
00274     }
00275 
00276 
00277   /*----- Set up columns of X matrix corresponding to ideal time series -----*/
00278   for (is = 0;  is < num_ideal_files;  is++)
00279     {
00280       far = MRI_FLOAT_PTR (ideal_array[is]);
00281       nx = ideal_array[is]->nx;
00282       ny = ideal_array[is]->ny;
00283 
00284       if (ideal_list[is] == NULL)
00285         for (iq = 0;  iq < ny;  iq++)
00286           {
00287             for (n = 0;  n < N;  n++)
00288               {
00289                 i = n + NFirst;
00290                 (*x).elts[n][m] = *(far + iq*nx + i);
00291               }
00292             
00293             m++;
00294           }
00295       else
00296         {
00297           nq = ideal_list[is][0];
00298           for (iq = 1;  iq <= nq;  iq++)
00299             {
00300               for (n = 0;  n < N;  n++)
00301                 {
00302                   i = n + NFirst;
00303                   (*x).elts[n][m] = *(far + ideal_list[is][iq]*nx + i);
00304                 }
00305 
00306               m++;
00307             }
00308         }
00309     }
00310 
00311 
00312   /*----- Remove row if ideal contains value over 33333 -----*/
00313   Ngood = N;
00314   m = 0;
00315   for (n = 0;  n < N;  n++)
00316     {
00317       for (is = q;  is < p;  is++)
00318         {
00319           if ((*x).elts[n][is] >= BIG_NUMBER)  break;
00320         }
00321       if (is < p)
00322         {
00323           Ngood--;
00324         }
00325       else
00326         {
00327           good_list[m] = n;
00328           m++;
00329         }
00330     }
00331   matrix_extract_rows ((*x), Ngood, good_list, &xgood);
00332   matrix_equate (xgood, x);
00333 
00334 
00335   /*----- Find min, max, and ave for each column of the X matrix -----*/
00336   for (is = 0;  is < p;  is++)
00337     {      
00338       x_bot[is] = x_top[is] = (*x).elts[0][is];
00339       x_ave[is] = 0.0;
00340       for (n = 0;  n < Ngood;  n++)
00341         {
00342           if ((*x).elts[n][is] < x_bot[is])  x_bot[is] = (*x).elts[n][is];  
00343           if ((*x).elts[n][is] > x_top[is])  x_top[is] = (*x).elts[n][is];
00344           x_ave[is] += (*x).elts[n][is] / Ngood;
00345         }
00346     }
00347   
00348   
00349   matrix_destroy (&xgood);
00350 
00351   return (Ngood);
00352 
00353 }
00354 
00355 
00356 /*---------------------------------------------------------------------------*/
00357 /*
00358   Initialization for the regression analysis.
00359 */
00360 
00361 int init_regression_analysis 
00362 (
00363   int q,                      /* number of parameters in the baseline model */
00364   int p,                      /* number of parameters in the baseline model 
00365                                  plus number of ideals */
00366   int N,                      /* total number of images used in the analysis */
00367   int num_idealts,            /* number of ideal time series */
00368   matrix xdata,               /* independent variable matrix */
00369   matrix * x_base,            /* extracted X matrix    for baseline model */
00370   matrix * xtxinvxt_base,     /* matrix:  (1/(X'X))X'  for baseline model */
00371   matrix * x_ideal,           /* extracted X matrices  for ideal models */
00372   matrix * xtxinvxt_ideal,    /* matrix:  (1/(X'X))X'  for ideal models */
00373   float ** rarray             /* ranked arrays of ideal time series */
00374 )
00375 
00376 {
00377   int * plist = NULL;         /* list of model parameters */
00378   int ip, it;                 /* parameter indices */
00379   int is, js;                 /* ideal indices */ 
00380   int jm;                     /* lag index */
00381   int ok;                     /* flag for successful matrix calculation */
00382   matrix xtxinv_temp;         /* intermediate results */
00383   vector ideal;               /* ideal vector */
00384   vector coef_temp;           /* intermediate results */
00385   vector xres;                /* vector of residuals */
00386   float sse_base;             /* baseline model error sum of squares */ 
00387         
00388 
00389   /*----- Initialize matrix -----*/
00390   matrix_initialize (&xtxinv_temp);
00391   vector_initialize (&ideal);
00392   vector_initialize (&coef_temp);
00393   vector_initialize (&xres);
00394 
00395 
00396   /*----- Allocate memory -----*/
00397   plist = (int *) malloc (sizeof(int) * p);   MTEST (plist);
00398 
00399 
00400   /*----- Initialize matrices for the baseline model -----*/
00401   for (ip = 0;  ip < q;  ip++)
00402     plist[ip] = ip;
00403   ok = calc_matrices (xdata, q, plist, x_base, &xtxinv_temp, xtxinvxt_base);
00404   if (!ok)  { matrix_destroy (&xtxinv_temp);  return (0); };
00405 
00406 
00407   /*----- Initialize matrices for ideal functions -----*/
00408   for (is = 0;  is < num_idealts;  is++)
00409     {
00410       for (ip = 0;  ip < q;  ip++)
00411         {
00412           plist[ip] = ip;
00413         }
00414 
00415       plist[q] = q + is;
00416 
00417       ok = calc_matrices (xdata, q+1, plist, 
00418                           &(x_ideal[is]), &xtxinv_temp, &(xtxinvxt_ideal[is]));
00419       if (!ok)  { matrix_destroy (&xtxinv_temp);  return (0); };
00420     }
00421 
00422 
00423   /*----- Set up the ranked array for each ideal -----*/
00424   for (is = 0;  is < num_idealts;  is++)
00425     {
00426       /*----- Convert column of matrix to vector -----*/
00427       column_to_vector (xdata, q+is, &ideal);
00428 
00429       /*----- Calculate regression coefficients for baseline model -----*/
00430       calc_coef (*xtxinvxt_base, ideal, &coef_temp);
00431 
00432       /*----- Calculate the error sum of squares for the baseline model -----*/
00433       sse_base = calc_resids (*x_base, coef_temp, ideal, &xres);
00434     
00435       /*----- Form rank array from residual array -----*/
00436       rarray[is] = rank_darray (N, xres.elts);
00437 
00438     }
00439 
00440 
00441   /*----- Destroy matrix -----*/
00442   matrix_destroy (&xtxinv_temp);
00443   vector_destroy (&ideal);
00444   vector_destroy (&coef_temp);
00445   vector_destroy (&xres);
00446 
00447 
00448   /*----- Deallocate memory -----*/
00449   free (plist);   plist = NULL;
00450 
00451 
00452   return (1);
00453 }
00454 
00455 
00456 /*---------------------------------------------------------------------------*/
00457 /* 
00458    Calculate a correlation coefficient. 
00459 */
00460 
00461 float calc_CC
00462 (
00463   int N,
00464   float * x,
00465   float * y
00466 )
00467 
00468 {
00469   const float EPSILON = 1.0e-10;   /* protect against divide by zero */
00470   float cc;
00471   int i;
00472   float csum, xsum, ysum, prod;
00473 
00474 
00475   csum = xsum = ysum = 0.0;
00476 
00477   for (i = 0;  i < N;  i++)
00478     {
00479       csum += x[i] * y[i];
00480       xsum += x[i] * x[i];
00481       ysum += y[i] * y[i];
00482     }
00483 
00484 
00485   prod = xsum * ysum;
00486   if (prod < EPSILON)
00487     cc = 0.0;
00488   else
00489     cc = csum / sqrt(prod);
00490 
00491 
00492   return (cc);
00493  
00494 }
00495 
00496 
00497 /*---------------------------------------------------------------------------*/
00498 /* 
00499    Calculate the Spearman correlation coefficient. 
00500 */
00501 
00502 float calc_SpearmanCC
00503 (
00504   int N,
00505   float * r,
00506   float * s
00507 )
00508 
00509 {
00510   float cc;
00511   float rbar;
00512   float * dr = NULL;
00513   float * ds = NULL;
00514   int i;
00515 
00516   
00517   dr = (float *) malloc (sizeof(float) * N);
00518   ds = (float *) malloc (sizeof(float) * N);
00519 
00520 
00521   rbar = (N+1.0) / 2.0;
00522   for (i = 0;  i < N;  i++)
00523     {
00524       dr[i] = r[i] - rbar;
00525       ds[i] = s[i] - rbar;
00526     }
00527 
00528 
00529   cc = calc_CC(N, dr, ds);
00530 
00531   free (dr);   dr = NULL;
00532   free (ds);   ds = NULL;
00533  
00534   return (cc);
00535 }
00536 
00537 
00538 /*---------------------------------------------------------------------------*/
00539 /* 
00540    Calculate the sign function.
00541 */
00542 
00543 float sign (float x)
00544 
00545 {
00546   if (x > 0.0)  return (1.0);
00547   if (x < 0.0)  return (-1.0);
00548   return (0.0);
00549 }
00550 
00551 
00552 /*---------------------------------------------------------------------------*/
00553 /* 
00554    Calculate the Quadrant correlation coefficient. 
00555 */
00556 
00557 float calc_QuadrantCC
00558 (
00559   int N,
00560   float * r,
00561   float * s
00562 )
00563 
00564 {
00565   float cc;
00566   float rbar;
00567   float * dr = NULL;
00568   float * ds = NULL;
00569   int i;
00570 
00571   
00572   dr = (float *) malloc (sizeof(float) * N);
00573   ds = (float *) malloc (sizeof(float) * N);
00574 
00575 
00576   rbar = (N+1.0) / 2.0;
00577   for (i = 0;  i < N;  i++)
00578     {
00579       dr[i] = sign(r[i] - rbar);
00580       ds[i] = sign(s[i] - rbar);
00581     }
00582 
00583 
00584   cc = calc_CC(N, dr, ds);
00585 
00586   free (dr);   dr = NULL;
00587   free (ds);   ds = NULL;
00588  
00589   return (cc);
00590 }
00591 
00592 
00593 
00594 /*---------------------------------------------------------------------------*/
00595 /* 
00596    Calculate percentage change relative to baseline. 
00597 */
00598 
00599 float percent_change (float num, float den)
00600 {
00601   const float EPSILON = 1.0e-10;      /* guard against divide by zero */
00602   const float MAX_PERCENT = 1000.0;   /* limit maximum percent change */ 
00603   float PrcntChange;
00604 
00605 
00606   if (fabs(den) < EPSILON)
00607     PrcntChange = sign(num) * MAX_PERCENT;
00608   else
00609     PrcntChange  = 100.0 * num / den;
00610  
00611   if (PrcntChange >  MAX_PERCENT)   PrcntChange =  MAX_PERCENT;
00612   if (PrcntChange < -MAX_PERCENT)   PrcntChange = -MAX_PERCENT;
00613 
00614   return (PrcntChange);
00615 }
00616 
00617 
00618 /*---------------------------------------------------------------------------*/
00619 /* 
00620    Calculate results for this voxel. 
00621 */
00622 
00623 void regression_analysis 
00624 (
00625   int N,                    /* number of usable data points */
00626   int q,                    /* number of parameters in baseline model */
00627   int num_idealts,          /* number of ideal time series */
00628   matrix x_base,            /* extracted X matrix    for baseline model */
00629   matrix xtxinvxt_base,     /* matrix:  (1/(X'X))X'  for baseline model */
00630   matrix * x_ideal,         /* extracted X matrices  for ideal models */
00631   matrix * xtxinvxt_ideal,  /* matrix:  (1/(X'X))X'  for ideal models */
00632   vector y,                 /* vector of measured data */
00633   float * x_bot,            /* minimum of stimulus time series */
00634   float * x_ave,            /* average of stimulus time series */
00635   float * x_top,            /* maximum of stimulus time series */
00636   float ** rarray,          /* ranked arrays of ideal time series */
00637   int * output_type,        /* list of operator requested outputs */
00638   float * FimParams         /* output fim parameters */
00639 )
00640 
00641 {
00642   const float EPSILON = 1.0e-05;   /* protection against divide by zero */
00643   int is;                   /* input ideal index */
00644   float sse_base;           /* error sum of squares, baseline model */
00645   float sse_ideal;          /* error sum of squares, ideal model */
00646   vector coef_temp;         /* intermediate results */
00647   vector coef_best;         /* best results */
00648   vector yres;              /* vector of residuals */
00649   float rtemp, rbest;       /* best correlation coefficient */  
00650   float stemp, sbest;       /* best Spearman correlation coefficient */
00651   float qtemp, qbest;       /* best quadrant correlation coefficient */
00652   float mse;                /* mean square error (sample variance) */
00653 
00654   float * sarray = NULL;    /* ranked array of measurement data */
00655 
00656                             /* fim output parameters */
00657   float FitCoef      = 0.0; 
00658   int   BestIndex    = 0;
00659   float PrcntChange  = 0.0; 
00660   float Baseline     = 0.0; 
00661   float Correlation  = 0.0; 
00662   float PrcntFromAve = 0.0;
00663   float Average      = 0.0; 
00664   float PrcntFromTop = 0.0; 
00665   float Topline      = 0.0; 
00666   float SigmaResid   = 0.0; 
00667   float SpearmanCC   = 0.0; 
00668   float QuadrantCC   = 0.0;
00669 
00670 
00671   /*----- Initialization -----*/
00672   vector_initialize (&coef_temp);
00673   vector_initialize (&coef_best);
00674   vector_initialize (&yres);
00675 
00676 
00677   /*----- Calculate regression coefficients for baseline model -----*/
00678   calc_coef (xtxinvxt_base, y, &coef_temp);
00679 
00680 
00681   /*----- Calculate the error sum of squares for the baseline model -----*/ 
00682   sse_base = calc_resids (x_base, coef_temp, y, &yres);
00683 
00684     
00685   /*----- Form rank array from y array -----*/
00686   if (output_type[FIM_SpearmanCC] || output_type[FIM_QuadrantCC]) 
00687     sarray = rank_darray (N, yres.elts);
00688 
00689 
00690   /*----- Determine the best ideal reference for this voxel -----*/
00691   rbest = 0.0;  sbest = 0.0;  qbest = 0.0;  BestIndex = -1;
00692   for (is = 0;  is < num_idealts;  is++)
00693     {
00694 
00695       /*----- Calculate regression coefficients for ideal model -----*/
00696       calc_coef (xtxinvxt_ideal[is], y, &coef_temp);
00697 
00698 
00699       /*----- Calculate the error sum of squares for the ideal model -----*/ 
00700       sse_ideal = calc_sse (x_ideal[is], coef_temp, y);
00701 
00702 
00703       /*----- Calculate partial R^2 for this ideal -----*/
00704       rtemp = calc_rsqr (sse_ideal, sse_base);
00705 
00706 
00707       if (rtemp >= rbest)
00708         {
00709           rbest = rtemp;
00710           BestIndex = is;
00711           vector_equate (coef_temp, &coef_best);
00712           if (num_idealts == 1)
00713             mse = sse_ideal / (N-q-1);
00714           else
00715             mse = sse_ideal / (N-q-2);
00716         }
00717 
00718 
00719       /*----- Calculate the Spearman rank correlation coefficient -----*/
00720       if (output_type[FIM_SpearmanCC])
00721         { 
00722           stemp = calc_SpearmanCC (N, rarray[is], sarray);
00723           if (fabs(stemp) > fabs(sbest))  sbest = stemp;
00724         }
00725 
00726 
00727       /*----- Calculate the Quadrant correlation coefficient -----*/
00728       if (output_type[FIM_QuadrantCC])
00729         { 
00730           qtemp = calc_QuadrantCC (N, rarray[is], sarray);
00731           if (fabs(qtemp) > fabs(qbest))  qbest = qtemp;
00732         }
00733     }
00734   
00735   if ((0 <= BestIndex) && (BestIndex < num_idealts))
00736     {
00737       float Top, Ave, Base, Center;
00738       int ip;
00739       
00740       Top  = x_top[q+BestIndex];  
00741       Ave  = x_ave[q+BestIndex];  
00742       Base = x_bot[q+BestIndex];
00743       
00744       
00745       FitCoef = coef_best.elts[q];
00746       Correlation = sqrt(rbest);
00747       if (FitCoef < 0.0)  Correlation = -Correlation;
00748       
00749       Center = 0.0;
00750       for (ip = 0;  ip < q;  ip++)
00751         Center += coef_best.elts[ip] * x_ave[ip];
00752       
00753       Baseline = Center + FitCoef*Base;
00754       Average  = Center + FitCoef*Ave;
00755       Topline  = Center + FitCoef*Top;
00756 
00757       
00758       
00759       PrcntChange  = percent_change (FitCoef * (Top-Base), Baseline);
00760       PrcntFromAve = percent_change (FitCoef * (Top-Base), Average);
00761       PrcntFromTop = percent_change (FitCoef * (Top-Base), Topline);
00762       
00763       SigmaResid = sqrt(mse);
00764       
00765       SpearmanCC = sbest;
00766       
00767       QuadrantCC = qbest;
00768     }
00769   
00770 
00771   /*----- Save output parameters -----*/
00772   FimParams[FIM_FitCoef]      = FitCoef;
00773   FimParams[FIM_BestIndex]    = BestIndex + 1.0;
00774   FimParams[FIM_PrcntChange]  = PrcntChange;
00775   FimParams[FIM_Baseline]     = Baseline;
00776   FimParams[FIM_Correlation]  = Correlation;
00777   FimParams[FIM_PrcntFromAve] = PrcntFromAve;
00778   FimParams[FIM_Average]      = Average;
00779   FimParams[FIM_PrcntFromTop] = PrcntFromTop;
00780   FimParams[FIM_Topline]      = Topline;
00781   FimParams[FIM_SigmaResid]   = SigmaResid;
00782   FimParams[FIM_SpearmanCC]   = SpearmanCC;
00783   FimParams[FIM_QuadrantCC]   = QuadrantCC;
00784 
00785 
00786   /*----- Dispose of vectors -----*/
00787   vector_destroy (&coef_temp);
00788   vector_destroy (&coef_best);
00789   vector_destroy (&yres);
00790 
00791 
00792   /*----- Deallocate memory -----*/
00793   if (sarray != NULL)
00794     { free (sarray);  sarray = NULL; }
00795 
00796 }
00797   
00798 
00799 /*---------------------------------------------------------------------------*/
00800 
00801 static char lbuf[4096];   /* character string containing statistical summary */
00802 static char sbuf[256];
00803 
00804 
00805 void report_results 
00806 (
00807   int * output_type,          /* output type options */
00808   float * FimParams,          /* output fim parameters */
00809   char ** label               /* statistical summary for ouput display */
00810 )
00811 
00812 {
00813   int ip;                     /* parameter index */
00814 
00815   
00816   if( label != NULL ){  /* assemble this 1 line at a time from sbuf */
00817     
00818     lbuf[0] = '\0' ;   /* make this a 0 length string to start */
00819     
00820     /** for each reference, make a string into sbuf **/
00821 
00822     for (ip = 0;  ip < MAX_OUTPUT_TYPE;  ip++)
00823       if (output_type[ip])
00824           {
00825             sprintf (sbuf, "%12s = %10.4f \n", 
00826                      OUTPUT_TYPE_name[ip], FimParams[ip]);
00827             strcat (lbuf, sbuf);
00828           }    
00829 
00830     
00831     *label = lbuf ;  /* send address of lbuf back in what label points to */
00832   }
00833   
00834 }
00835 
00836 
00837 /*---------------------------------------------------------------------------*/
00838 
00839 
00840 
00841