Doxygen Source Code Documentation
fim+.c File Reference
#include "RegAna.c"#include "ranks.c"Go to the source code of this file.
| Defines | |
| #define | MAX_OUTPUT_TYPE 12 | 
| #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 | 
| Functions | |
| double | legendre (double x, int m) | 
| int | init_indep_var_matrix (int q, int p, int NFirst, int N, int polort, int num_ort_files, int num_ideal_files, MRI_IMAGE **ort_array, int **ort_list, MRI_IMAGE **ideal_array, int **ideal_list, float *x_bot, float *x_ave, float *x_top, int *good_list, matrix *x) | 
| int | init_regression_analysis (int q, int p, int N, int num_idealts, matrix xdata, matrix *x_base, matrix *xtxinvxt_base, matrix *x_ideal, matrix *xtxinvxt_ideal, float **rarray) | 
| float | calc_CC (int N, float *x, float *y) | 
| float | calc_SpearmanCC (int N, float *r, float *s) | 
| float | sign (float x) | 
| float | calc_QuadrantCC (int N, float *r, float *s) | 
| float | percent_change (float num, float den) | 
| void | regression_analysis (int N, int q, int num_idealts, matrix x_base, matrix xtxinvxt_base, matrix *x_ideal, matrix *xtxinvxt_ideal, vector y, float *x_bot, float *x_ave, float *x_top, float **rarray, int *output_type, float *FimParams) | 
| void | report_results (int *output_type, float *FimParams, char **label) | 
| Variables | |
| char * | OUTPUT_TYPE_name [MAX_OUTPUT_TYPE] | 
| char | lbuf [4096] | 
| char | sbuf [256] | 
Define Documentation
| 
 | 
| 
 Definition at line 48 of file fim+.c. Referenced by display_help_menu(), and regression_analysis(). | 
| 
 | 
| 
 Definition at line 45 of file fim+.c. Referenced by display_help_menu(), and regression_analysis(). | 
| 
 | 
| 
 Definition at line 43 of file fim+.c. Referenced by check_for_valid_inputs(), display_help_menu(), and regression_analysis(). | 
| 
 | 
| 
 Definition at line 46 of file fim+.c. Referenced by calculate_results(), display_help_menu(), regression_analysis(), and write_bucket_data(). | 
| 
 | 
| 
 Definition at line 42 of file fim+.c. Referenced by display_help_menu(), and regression_analysis(). | 
| 
 | 
| 
 Definition at line 44 of file fim+.c. Referenced by display_help_menu(), and regression_analysis(). | 
| 
 | 
| 
 Definition at line 47 of file fim+.c. Referenced by display_help_menu(), and regression_analysis(). | 
| 
 | 
| 
 Definition at line 49 of file fim+.c. Referenced by display_help_menu(), and regression_analysis(). | 
| 
 | 
| 
 Definition at line 53 of file fim+.c. Referenced by display_help_menu(), regression_analysis(), and write_bucket_data(). | 
| 
 | 
| 
 Definition at line 51 of file fim+.c. Referenced by display_help_menu(), and regression_analysis(). | 
| 
 | 
| 
 Definition at line 52 of file fim+.c. Referenced by display_help_menu(), regression_analysis(), and write_bucket_data(). | 
| 
 | 
| 
 Definition at line 50 of file fim+.c. Referenced by display_help_menu(), and regression_analysis(). | 
| 
 | 
| 
 Definition at line 34 of file fim+.c. Referenced by allocate_memory(), calculate_results(), get_options(), initialize_options(), report_results(), save_voxel(), terminate_program(), and write_bucket_data(). | 
| 
 | 
| 
 | 
Function Documentation
| 
 | ||||||||||||||||
| 
 Definition at line 462 of file fim+.c. References i. Referenced by calc_QuadrantCC(), and calc_SpearmanCC(). 
 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 }
 | 
| 
 | ||||||||||||||||
| 
 Definition at line 558 of file fim+.c. References calc_CC(), free, i, malloc, r, and sign. Referenced by regression_analysis(). 
 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 }
 | 
| 
 | ||||||||||||||||
| 
 Definition at line 503 of file fim+.c. References calc_CC(), free, i, malloc, and r. Referenced by regression_analysis(). 
 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 }
 | 
| 
 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| if here, m > 20 ==> use recurrence relation * Definition at line 185 of file fim+.c. References BIG_NUMBER, far, i, legendre(), matrix_create(), matrix_destroy(), matrix_equate(), matrix_extract_rows(), matrix_initialize(), MRI_FLOAT_PTR, p, and q. Referenced by calculate_results(). 
 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 }
 | 
| 
 | ||||||||||||||||||||||||||||||||||||||||||||
| the new polort way: Legendre P_m(t) - 29 Mar 2005 * Definition at line 362 of file fim+.c. References calc_coef(), calc_matrices(), calc_resids(), column_to_vector(), vector::elts, free, malloc, matrix_destroy(), matrix_initialize(), MTEST, p, q, rank_darray(), vector_destroy(), and vector_initialize(). Referenced by calculate_results(). 
 00370 : (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 } | 
| 
 | ||||||||||||
| 
 Definition at line 60 of file fim+.c. Referenced by init_indep_var_matrix(), and legendre(). 
 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 }
 | 
| 
 | ||||||||||||
| 
 Definition at line 599 of file fim+.c. References sign. Referenced by regression_analysis(). 
 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 }
 | 
| 
 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 
 Definition at line 624 of file fim+.c. References calc_coef(), calc_QuadrantCC(), calc_resids(), calc_rsqr(), calc_SpearmanCC(), calc_sse(), vector::elts, FIM_Average, FIM_Baseline, FIM_BestIndex, FIM_Correlation, FIM_FitCoef, FIM_PrcntChange, FIM_PrcntFromAve, FIM_PrcntFromTop, FIM_QuadrantCC, FIM_SigmaResid, FIM_SpearmanCC, FIM_Topline, free, percent_change(), q, rank_darray(), vector_destroy(), vector_equate(), and vector_initialize(). Referenced by calculate_results(). 
 00629 : (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 } | 
| 
 | ||||||||||||||||
| 
 Definition at line 806 of file fim+.c. References lbuf, MAX_OUTPUT_TYPE, OUTPUT_TYPE_name, and sbuf. 
 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 }
 | 
| 
 | 
| 
 Definition at line 543 of file fim+.c. 
 00545 {
00546   if (x > 0.0)  return (1.0);
00547   if (x < 0.0)  return (-1.0);
00548   return (0.0);
00549 }
 | 
Variable Documentation
| 
 | 
| 
 Definition at line 801 of file fim+.c. Referenced by report_results(). | 
| 
 | 
| Initial value: 
  { "Fit Coef", "Best Index", "% Change", "% From Ave", "Baseline", "Average", 
    "Correlation", "% From Top", "Topline", "Sigma Resid",
    "Spearman CC", "Quadrant CC" }Definition at line 36 of file fim+.c. Referenced by report_results(). | 
| 
 | 
| 
 Definition at line 802 of file fim+.c. Referenced by report_results(). | 
 
                             
                             
                             
                             
                             
                             
                             
                             
                             
                             
                             
                             
 
 
 
 
       
	   
	   
	   
	  