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 Simplex algorithm for 
00010   non-linear optimization to estimate the unknown parameters.
00011 
00012   The Simplex algorithm is adapted from: A Jump Start Course in C++ Programming
00013   
00014   File:    Simplex.c
00015   Author:  B. Douglas Ward
00016   Date:    28 January 2000
00017 
00018 
00019   Note:  The calling program should include the following statements
00020          prior to "#include "Simplex.c"
00021 
00022   #define DIMENSION p    where p = number of parameters to be estimated
00023   #include "randgen.c"
00024   float calc_error (float * vertex);
00025 */
00026 
00027 
00028 /*---------------------------------------------------------------------------*/
00029 /*
00030   Global variables and constants.
00031 */
00032 
00033 int count = 0;               /* count of error evaluations */
00034 int number_restarts = 0;     /* count of simplex algorithm restarts */
00035 
00036 
00037 /*---------------------------------------------------------------------------*/
00038 
00039 void allocate_arrays (float *** simplex, float ** centroid,
00040                       float ** response, float ** step_size, 
00041                       float ** test1, float ** test2)
00042 {
00043   int i;
00044 
00045   *centroid = (float *) malloc (sizeof(float) * DIMENSION);
00046   *response = (float *) malloc (sizeof(float) * (DIMENSION+1));
00047   *step_size = (float *) malloc (sizeof(float) * DIMENSION);
00048   *test1 = (float *) malloc (sizeof(float) * DIMENSION);
00049   *test2 = (float *) malloc (sizeof(float) * DIMENSION);
00050    
00051   *simplex = (float **) malloc (sizeof(float *) * (DIMENSION+1));
00052 
00053   for (i = 0;  i < DIMENSION+1;  i++)
00054     (*simplex)[i] = (float *) malloc (sizeof(float) * DIMENSION);
00055        
00056 }
00057 
00058 
00059 /*---------------------------------------------------------------------------*/
00060 
00061 void deallocate_arrays (float *** simplex, float ** centroid,
00062                         float ** response, float ** step_size, 
00063                         float ** test1, float ** test2)
00064 {
00065   int i;
00066 
00067   free (*centroid);    *centroid = NULL;
00068   free (*response);    *response = NULL;
00069   free (*step_size);   *step_size = NULL;
00070   free (*test1);       *test1 = NULL;
00071   free (*test2);       *test2 = NULL;
00072 
00073   for (i = 0;  i < DIMENSION+1;  i++)
00074     {
00075       free ((*simplex)[i]);
00076       (*simplex)[i] = NULL;
00077     }
00078 
00079   free (*simplex);     *simplex = NULL;
00080        
00081 }
00082 
00083 
00084 /*---------------------------------------------------------------------------*/
00085 
00086 void eval_vertices (float * response, int * worst, int * next, int * best)
00087 {
00088   int i;
00089 
00090   /* initialize values */
00091   *worst = 0;
00092   *best = 0;
00093 
00094   /* find the best and worst */
00095   for (i = 1;  i < DIMENSION+1;  i++)
00096     {
00097       if (response[i] > response[*worst])
00098         *worst = i;
00099       if (response[i] < response[*best])
00100         *best = i;
00101     }
00102 
00103   /* find the next worst index */
00104   if (*worst == 0)
00105     *next = 1;
00106   else
00107     *next = 0;
00108   
00109   for (i = 0;  i < DIMENSION+1;  i++)
00110     if ((i != *worst) && (response[i] > response[*next]))
00111       *next = i;
00112 }
00113 
00114 
00115 /*---------------------------------------------------------------------------*/
00116 
00117 void restart (float ** simplex, float * response, 
00118               float * step_size)
00119 {
00120   const float STEP_FACTOR = 0.9;
00121   int i, j;
00122   int worst, next, best;
00123   float minval, maxval;
00124 
00125 
00126   /* find the current best vertex */
00127   eval_vertices (response, &worst, &next, &best); 
00128 
00129   /* set the first vertex to the current best */
00130   for (i = 0; i < DIMENSION;  i++)
00131     simplex[0][i] = simplex[best][i];
00132 
00133   /* decrease step size */
00134   for (i = 0;  i < DIMENSION;  i++)
00135     step_size[i] *= STEP_FACTOR;
00136 
00137   /* set up remaining vertices of simplex using new step size */
00138   for (i = 1;  i < DIMENSION+1;  i++)
00139     for (j = 0;  j < DIMENSION;  j++)
00140       {
00141         minval = simplex[0][j] - step_size[j];
00142         maxval = simplex[0][j] + step_size[j];
00143       simplex[i][j] = rand_uniform (minval, maxval);
00144       }
00145 
00146   /* initialize response for each vector */
00147   for (i = 0;  i < DIMENSION+1;  i++)
00148     response[i] = calc_error (simplex[i]);
00149 }
00150 
00151 
00152 /*---------------------------------------------------------------------------*/
00153 /*
00154   Calculate the centroid of the simplex, ignoring the worst vertex.
00155 */
00156 
00157 void calc_centroid (float ** simplex, int worst, float * centroid)
00158 {
00159   int i, j;
00160 
00161   for (i = 0;  i < DIMENSION;  i++)
00162     {
00163       centroid[i] = 0.0;
00164 
00165       /* add each vertex, except the worst */
00166       for (j = 0;  j < DIMENSION+1;  j++)
00167         if (j != worst)
00168           centroid[i] += simplex[j][i];
00169     }
00170 
00171   /* divide by the number of vertices */
00172   for (i = 0;  i < DIMENSION;  i++)
00173     centroid[i] /= DIMENSION;
00174 }
00175 
00176 
00177 /*---------------------------------------------------------------------------*/
00178 /*
00179   Calculate the reflection of the worst vertex about the centroid.
00180 */
00181 
00182 void calc_reflection (float ** simplex, float * centroid, 
00183                       int worst, float coef, float * vertex)
00184 {
00185   int i;
00186 
00187   for (i = 0;  i < DIMENSION;  i++)
00188     vertex[i] = centroid[i] + coef*(centroid[i] - simplex[worst][i]);
00189 }
00190 
00191 
00192 /*---------------------------------------------------------------------------*/
00193 /*
00194   Replace a vertex of the simplex.
00195 */
00196 
00197 void replace (float ** simplex, float * response, 
00198               int index, float * vertex, float resp)
00199 {
00200   int i;
00201 
00202   for (i = 0;  i < DIMENSION;  i++)
00203     simplex[index][i] = vertex[i];
00204 
00205   response[index] = resp;
00206 }
00207 
00208 
00209 /*---------------------------------------------------------------------------*/
00210 /*
00211   Calculate goodness of fit.
00212 */
00213 
00214 float calc_good_fit (float * response)
00215 {
00216   int i;
00217 
00218   float avg, sd, tmp;
00219 
00220   /* average the response values */
00221   avg = 0.0;
00222   for (i = 0;  i < DIMENSION+1;  i++)
00223     avg += response[i];
00224   avg /= DIMENSION+1;
00225 
00226   /* compute standard deviation of response */
00227   sd = 0.0;
00228   for (i = 0;  i < DIMENSION+1;  i++)
00229     {
00230       tmp = response[i] - avg;
00231       sd += tmp*tmp;
00232     }
00233   sd /= DIMENSION;
00234 
00235   return (sqrt(sd));
00236 }
00237 
00238 
00239 /*---------------------------------------------------------------------------*/
00240 /*
00241   Perform initialization for the Simplex algorithm.
00242 */
00243 
00244 void simplex_initialize (float * parameters, float ** simplex, 
00245                          float * response, float * step_size)
00246 {
00247   int i, j;
00248   int worst, next, best;
00249   float resp;
00250   float minval, maxval;
00251 
00252 
00253   for (j = 0;  j < DIMENSION;  j++)
00254     {
00255       simplex[0][j] = parameters[j];
00256       step_size[j] = 0.5 * parameters[j];
00257     }
00258 
00259   for (i = 1;  i < DIMENSION+1;  i++)
00260     for (j = 0;  j < DIMENSION;  j++)
00261       {
00262         minval = simplex[0][j] - step_size[j];
00263         maxval = simplex[0][j] + step_size[j];
00264         simplex[i][j] = rand_uniform (minval, maxval);
00265       }
00266 
00267   for (i = 0;  i < DIMENSION+1;  i++)
00268       response[i] = calc_error (simplex[i]);
00269 
00270   for (i = 1;  i < 500;  i++)
00271     {
00272       for (j = 0;  j < DIMENSION;  j++)
00273         {
00274           minval = simplex[0][j] - step_size[j];
00275           maxval = simplex[0][j] + step_size[j];
00276           parameters[j] = rand_uniform (minval, maxval);
00277         }
00278 
00279       resp = calc_error (parameters);
00280       eval_vertices (response, &worst, &next, &best);
00281       if (resp < response[worst])
00282         replace (simplex, response, worst, parameters, resp);
00283     }
00284 }
00285 
00286 
00287 /*---------------------------------------------------------------------------*/
00288 /*
00289   Use Simplex algorithm to estimate the voxel intensity distribution.
00290 
00291   The Simplex algorithm is adapted from: A Jump Start Course in C++ Programming
00292 */
00293 
00294 void simplex_optimization (float * parameters, float * sse)
00295 {
00296   const int MAX_ITERATIONS = 100;
00297   const int MAX_RESTARTS = 25;
00298   const float EXPANSION_COEF = 2.0;
00299   const float REFLECTION_COEF = 1.0;
00300   const float CONTRACTION_COEF = 0.5;
00301   const float TOLERANCE = 1.0e-10;
00302 
00303   float ** simplex   = NULL;
00304   float * centroid   = NULL;
00305   float * response   = NULL;
00306   float * step_size  = NULL;
00307   float * test1      = NULL;
00308   float * test2      = NULL;
00309   float resp1, resp2;
00310   int i, worst, best, next;
00311   int num_iter, num_restarts;
00312   int done;
00313   float fit;
00314 
00315 
00316   allocate_arrays (&simplex, &centroid, &response, &step_size, &test1, &test2);
00317   
00318   simplex_initialize (parameters, simplex, response, step_size);
00319 
00320   /* start loop to do simplex optimization */
00321   num_iter = 0;
00322   num_restarts = 0;
00323   done = 0;
00324   
00325   while (!done)
00326     {
00327       /* find the worst vertex and compute centroid of remaining simplex, 
00328          discarding the worst vertex */
00329       eval_vertices (response, &worst, &next, &best);
00330       calc_centroid (simplex, worst, centroid);
00331       
00332       /* reflect the worst point through the centroid */
00333       calc_reflection (simplex, centroid, worst, 
00334                        REFLECTION_COEF, test1);
00335       resp1 = calc_error (test1);
00336 
00337       /* test the reflection against the best vertex and expand it if the
00338          reflection is better.  if not, keep the reflection */
00339       if (resp1 < response[best])
00340         {
00341           /* try expanding */
00342           calc_reflection (simplex, centroid, worst, EXPANSION_COEF, test2);
00343           resp2 = calc_error (test2);
00344           if (resp2 <= resp1)    /* keep expansion */     
00345             replace (simplex, response, worst, test2, resp2);
00346           else                   /* keep reflection */
00347             replace (simplex, response, worst, test1, resp1);
00348         }
00349       else if (resp1 < response[next])
00350         {
00351           /* new response is between the best and next worst 
00352              so keep reflection */
00353           replace (simplex, response, worst, test1, resp1); 
00354         }
00355           else
00356         {
00357           /* try contraction */
00358           if (resp1 >= response[worst])
00359             calc_reflection (simplex, centroid, worst, 
00360                              -CONTRACTION_COEF, test2);
00361           else
00362             calc_reflection (simplex, centroid, worst, 
00363                              CONTRACTION_COEF, test2);
00364           resp2 =  calc_error (test2);
00365           
00366           /* test the contracted response against the worst response */
00367           if (resp2 > response[worst])
00368             {
00369               /* new contracted response is worse, so decrease step size
00370                  and restart */
00371               num_iter = 0;
00372               num_restarts += 1;
00373               restart (simplex, response, step_size);
00374             }
00375           else       /* keep contraction */
00376             replace (simplex, response, worst, test2, resp2);
00377         }
00378 
00379       /* test to determine when to stop.  
00380          first, check the number of iterations */
00381       num_iter += 1;    /* increment iteration counter */
00382       if (num_iter >= MAX_ITERATIONS)
00383         {
00384           /* restart with smaller steps */
00385           num_iter = 0;
00386           num_restarts += 1;
00387           restart (simplex, response, step_size);
00388         }
00389 
00390       /* limit the number of restarts */
00391       if (num_restarts == MAX_RESTARTS)  done = 1;
00392 
00393       /* compare relative standard deviation of vertex responses 
00394          against a defined tolerance limit */
00395       fit = calc_good_fit (response);
00396       if (fit <= TOLERANCE)  done = 1;
00397 
00398       /* if done, copy the best solution to the output array */
00399       if (done) 
00400         {
00401           eval_vertices (response, &worst, &next, &best);
00402           for (i = 0;  i < DIMENSION;  i++)
00403             parameters[i] = simplex[best][i];
00404           *sse = response[best];
00405         }
00406 
00407     }  /* while (!done) */
00408  
00409   number_restarts = num_restarts;
00410   deallocate_arrays (&simplex, &centroid, &response, &step_size,
00411                      &test1, &test2);
00412 
00413 }
00414 
00415 
00416 /*---------------------------------------------------------------------------*/