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 program estimates the probability density function (PDF) 
00010   corresponding to the distribution of gray-matter and white-matter
00011   voxel intensities.  The estimate is formed as the sum of three normal
00012   distributions, using the Simplex algorithm for non-linear optimization
00013   to estimate the unknown parameters.
00014 
00015   The Simplex algorithm is adapted from: A Jump Start Course in C++ Programming
00016   
00017   File:    estpdf.c
00018   Author:  B. Douglas Ward
00019   Date:    27 May 1999
00020 */
00021 
00022 
00023 /*---------------------------------------------------------------------------*/
00024 /*
00025   Include source code.
00026 */
00027 
00028 #include "randgen.c"
00029 
00030 /*---------------------------------------------------------------------------*/
00031 /*
00032   Global variables and constants.
00033 */
00034 
00035 #define DIMENSION 9      /* number of parameters to be estimated */
00036 
00037 int * hist_data;         /* histogram of voxel gray-scale intensities */
00038 int hist_min;            /* minimum histogram bin value;  this is set to 
00039                             exclude the large population near zero */
00040 int hist_max;            /* maximum histogram bin value;  this is set to
00041                             exclude voxels with very high intensities */
00042 int gpeak;               /* estimated peak of gray-matter distribution */
00043 int wpeak;               /* estimated peak of white-matter distribution */
00044 int numpts;              /* number of voxels within the histogram limits */
00045 
00046 int count = 0;               /* count of sse evaluations */
00047 int number_restarts = 0;     /* count of simplex algorithm restarts */
00048 
00049 
00050 /*---------------------------------------------------------------------------*/
00051 /*
00052   Set up the histogram and generate some of the initial parameter estimates. 
00053 
00054   This routine was adapted from:  program 3dstrip.c
00055 */
00056 
00057 #undef MAX
00058 #define MAX(a,b) (((a)<(b)) ? (b) : (a))
00059 
00060 #undef MIN
00061 #define MIN(a,b) (((a)>(b)) ? (b) : (a))
00062 
00063 #define NPMAX 128
00064 
00065 
00066 void find_base_value( int nxyz , short * sfim )
00067 {
00068    float bper = 50.0 , bmin = 1 ;
00069    float tper = 98.0;
00070 
00071    int ii , kk , nbin , sval , sum , nbot , a,b,c , npeak,ntop , nvox ;
00072    int * fbin ;
00073    int   kmin[NPMAX] , kmax[NPMAX] ;
00074    int   nmin        , nmax        ;
00075 
00076    /*-- make histogram of shorts --*/
00077 
00078    fbin = (int *) malloc( sizeof(int) * 32768 ) ;
00079 
00080    for( kk=0 ; kk < 32768 ; kk++ ) fbin[kk] = 0 ;
00081 
00082    nvox = 0 ;
00083 
00084    for( ii=0 ; ii < nxyz ; ii++ ){
00085       kk = sfim[ii] ; if( kk >= 0 ){ fbin[kk]++ ; nvox++ ; }
00086    }
00087 
00088    /*-- find largest value --*/
00089 
00090    for( kk=32767 ; kk > 0 ; kk-- ) if( fbin[kk] > 0 ) break ;
00091    if( kk == 0 ){fprintf(stderr,"** find_base_value: All voxels are zero!\n");exit(1);}
00092    nbin = kk+1 ;
00093 
00094    /*-- find bper point in cumulative distribution --*/
00095 
00096    sval = 0.01 * bper * nvox ;
00097    sum  = 0 ;
00098    for( kk=0 ; kk < nbin ; kk++ ){
00099       sum += fbin[kk] ; if( sum >= sval ) break ;
00100    }
00101    nbot = kk ; if( nbot == 0 ) nbot = 1 ; if( bmin > nbot ) nbot = bmin ;
00102    if( nbot >= nbin-9 ){
00103       fprintf(stderr,"** find_base_value: Base point on histogram too high\n");
00104       exit(1);
00105    }
00106 
00107    hist_min = nbot;
00108 
00109 
00110    /*-- find tper point in cumulative distribution --*/
00111 
00112    sval = 0.01 * tper * nvox ;
00113    sum  = 0 ;
00114    for( kk=0 ; kk < nbin ; kk++ ){
00115       sum += fbin[kk] ; if( sum >= sval ) break ;
00116    }
00117 
00118    hist_max = kk ; 
00119 
00120 
00121    /*----- Save original histogram data -----*/
00122    hist_data = (int *) malloc( sizeof(int) * nbin ) ;
00123    for (kk = 0;  kk < nbin;  kk++)
00124      hist_data[kk] = fbin[kk];
00125 
00126 
00127    /*-- smooth histogram --*/
00128 
00129    b = fbin[nbot-1] ; c = fbin[nbot] ;
00130    for( kk=nbot ; kk < nbin ; kk++ ){
00131       a = b ; b = c ; c = fbin[kk+1] ; fbin[kk] = 0.25*(a+c+2*b) ;
00132    }
00133 
00134    /*-- find minima and maxima above bper point --*/
00135 
00136    nmin = nmax = 0 ;
00137    for( kk=nbot+1 ; kk < nbin ; kk++ ){
00138       if( fbin[kk] < fbin[kk-1] && fbin[kk] < fbin[kk+1] && nmin < NPMAX ){
00139          kmin[nmin++] = kk ;
00140       } else if( fbin[kk] > fbin[kk-1] && fbin[kk] > fbin[kk+1] && nmax < NPMAX ){
00141          kmax[nmax++] = kk ;
00142       }
00143    }
00144    
00145 
00146    /*-- find the largest two maxima --*/
00147 
00148    if( nmax == 0 ){
00149       fprintf(stderr,"** find_base_value: No histogram maxima above base point\n");
00150       exit(1);
00151    }
00152 
00153    if( nmax == 1 ){
00154       npeak = kmax[0] ; ntop = 0 ;
00155    } else {
00156       int f1,f2 , k1,k2 , fk , klow,kup ;
00157 
00158       k1 = 0 ; f1 = fbin[kmax[0]] ;
00159       k2 = 1 ; f2 = fbin[kmax[1]] ;
00160       if( f1 < f2 ){
00161          k1 = 1 ; f1 = fbin[kmax[1]] ;
00162          k2 = 0 ; f2 = fbin[kmax[0]] ;
00163       }
00164 
00165       for( kk=2 ; kk < nmax ; kk++ ){
00166          fk = fbin[kmax[kk]] ;
00167          if( fk > f1 ){
00168             f2 = f1 ; k2 = k1 ;
00169             f1 = fk ; k1 = kk ;
00170          } else if( fk > f2 ){
00171             f2 = fk ; k2 = kk ;
00172          }
00173       }
00174       npeak = MIN( kmax[k1] , kmax[k2] ) ;  /* smaller bin of the 2 top peaks */
00175 
00176       /* find valley between 2 peaks */
00177 
00178       ntop  = MAX( kmax[k1] , kmax[k2] ) ;
00179 
00180       gpeak = npeak;
00181       wpeak = ntop;
00182 
00183       
00184       fk = fbin[ntop] ; klow = ntop ;
00185       for( kk=ntop-1 ; kk >= npeak ; kk-- ){
00186          if( fbin[kk] < fk ){ fk = fbin[kk] ; klow = kk ; }
00187       }
00188       fk  = MAX( 0.10*fk , 0.05*fbin[ntop] ) ;
00189       kup = MIN( nbin-1 , ntop+3*(ntop-klow+2) ) ;
00190       for( kk=ntop+1 ; kk <= kup ; kk++ ) if( fbin[kk] < fk ) break ;
00191 
00192       ntop = kk ;
00193    }
00194 
00195    for( kk=npeak-1 ; kk > 0 ; kk-- )
00196       if( fbin[kk] < fbin[kk-1] && fbin[kk] < fbin[kk+1] ) break ;
00197 
00198 
00199    if (kk > hist_min)  hist_min = kk ;
00200 
00201    if( ntop == 0 )
00202      {
00203        ntop = npeak + (npeak-kk) ;
00204        gpeak = npeak;
00205        wpeak = ntop;
00206      }
00207 
00208    free (fbin);
00209 
00210    return ;
00211 }
00212 
00213 
00214 /*---------------------------------------------------------------------------*/
00215 /*
00216   Perform initialization for estimation of PDF. 
00217 */
00218 
00219 
00220 Boolean estpdf_initialize ()
00221 
00222 {
00223   int nx, ny, nz, nxy, nxyz, ixyz;       /* voxel counters */
00224   int n;                                 /* histogram bin index */
00225   short * sfim;                          /* pointer to anat data */
00226 
00227 
00228   /*----- Initialize local variables -----*/
00229   if (anat == NULL)  return (FALSE);
00230   nx = DSET_NX(anat);   ny = DSET_NY(anat);   nz = DSET_NZ(anat);
00231   nxy = nx*ny;   nxyz = nxy*nz;
00232   sfim  = (short *) DSET_BRICK_ARRAY(anat,0) ;
00233   if (sfim == NULL)  return (FALSE);
00234 
00235 
00236   /*----- Set up histogram, and get initial parameter estimates -----*/
00237   find_base_value (nxyz, sfim);
00238 
00239 
00240   /*----- Count number of voxels inside histogram intensity limits -----*/
00241   numpts = 0;
00242   for (n = hist_min;  n < hist_max;  n++)
00243     numpts += hist_data[n];
00244 
00245 
00246   /*----- Initialize random number generator -----*/
00247   srand48 (1234567);
00248 
00249 
00250   return (TRUE);
00251 }
00252 
00253 
00254 /*---------------------------------------------------------------------------*/
00255 /*
00256   Generate the initial guess for the parameter vector.
00257 */
00258 
00259 void generate_initial_guess (float * parameters)
00260 {
00261   float b;                   /* coefficient for background distribution */
00262   float bmean;               /* mean for background distribution */ 
00263   float bsigma;              /* std. dev. for background distribution */
00264   float g;                   /* coefficient for gray-matter distribution */
00265   float gmean;               /* mean for gray-matter distribution */
00266   float gsigma;              /* std. dev. for gray-matter distribution */
00267   float w;                   /* coefficient for white-matter distribution */
00268   float wmean;               /* mean for white-matter distribution */
00269   float wsigma;              /* std. dev. for white-matter distribution */
00270 
00271 
00272   /*----- Initialize distribution coefficients -----*/
00273   b = 0.75;
00274   g = 0.25;
00275   w = 0.25;
00276 
00277 
00278   /*----- Initialize distribution means -----*/
00279   bmean = hist_min;
00280 
00281   if ((gpeak > hist_min) && (gpeak < hist_max) && (gpeak < wpeak)) 
00282     gmean = gpeak;
00283   else
00284     gmean = hist_min;
00285 
00286   if ((wpeak > hist_min) && (wpeak < hist_max) && (wpeak > gpeak))
00287     wmean = wpeak;
00288   else
00289     wmean = hist_max;
00290 
00291   if ((gmean-bmean) < 0.25*(wmean-bmean))  gmean = bmean + 0.25*(wmean-bmean);
00292   if ((wmean-gmean) < 0.25*(wmean-bmean))  gmean = wmean - 0.25*(wmean-bmean);
00293 
00294 
00295   /*----- Initialize distribution standard deviations -----*/
00296   bsigma = 0.25 * (hist_max - hist_min);
00297   gsigma = 0.25 * (wmean - gmean);
00298   wsigma = 0.25 * (wmean - gmean);
00299 
00300 
00301   /*----- Set parameter vector -----*/
00302   parameters[0] = b;
00303   parameters[1] = bmean;
00304   parameters[2] = bsigma;
00305   parameters[3] = g;
00306   parameters[4] = gmean;
00307   parameters[5] = gsigma;
00308   parameters[6] = w;
00309   parameters[7] = wmean;
00310   parameters[8] = wsigma;
00311 
00312 }
00313 
00314 
00315 /*---------------------------------------------------------------------------*/
00316 /*
00317   Write parameter vector.
00318 */
00319 
00320 void write_parameter_vector (float * parameters)
00321 {
00322   int i;
00323 
00324   printf ("Dimension = %d \n", DIMENSION);
00325   for (i = 0;  i < DIMENSION;  i++)
00326     printf ("parameter[%d] = %f \n", i, parameters[i]);
00327 }
00328 
00329 
00330 /*---------------------------------------------------------------------------*/
00331 /*
00332   Normal probability density function.
00333 */
00334 
00335 float normal (float x, float mean, float sigma)
00336 
00337 {
00338   float z;
00339 
00340   z = (x - mean) / sigma;
00341 
00342   return ( (1.0/(sqrt(2.0*PI)*sigma)) * exp (-0.5 * z * z) );
00343 }
00344 
00345 
00346 /*---------------------------------------------------------------------------*/
00347 /*
00348   Estimate the voxel intensity distribution as the sum of three normal PDF's.
00349 */
00350 
00351 float estimate (float * parameters, float x)
00352 
00353 {
00354   float b, bmean, bsigma, g, gmean, gsigma, w, wmean, wsigma;
00355   float z, fval;
00356 
00357 
00358   /*----- Initialize local variables -----*/
00359   b      = parameters[0];
00360   bmean  = parameters[1];
00361   bsigma = parameters[2];
00362   g      = parameters[3];
00363   gmean  = parameters[4];
00364   gsigma = parameters[5];
00365   w      = parameters[6];
00366   wmean  = parameters[7];
00367   wsigma = parameters[8];
00368 
00369 
00370   /*----- Calculate the sum of three normal PDF's -----*/ 
00371   fval  = b * normal (x, bmean, bsigma);
00372   fval += g * normal (x, gmean, gsigma);
00373   fval += w * normal (x, wmean, wsigma);
00374 
00375 
00376   return (fval);
00377   
00378 }
00379 
00380 
00381 /*---------------------------------------------------------------------------*/
00382 /*
00383   Calculate the error sum of squares for the PDF estimate.
00384 */
00385 
00386 float calc_sse (float * vertex)
00387 
00388 {
00389   const float delt = 1.0;            /* histogram bin size */
00390   const float BIG_NUMBER = 1.0e+10;  /* return when constraints are violated */
00391 
00392   float b, bmean, bsigma, g, gmean, gsigma, w, wmean, wsigma;  /* parameters */
00393   float deltah, deltam;          /* rough estimate of spread of distribution */
00394 
00395   int i;
00396   float t;
00397   float diff, sse;
00398 
00399   count += 1;
00400 
00401 
00402   /*----- Assign local variables -----*/
00403   b      = vertex[0];
00404   bmean  = vertex[1];
00405   bsigma = vertex[2];
00406   g      = vertex[3];
00407   gmean  = vertex[4];
00408   gsigma = vertex[5];
00409   w      = vertex[6];
00410   wmean  = vertex[7];
00411   wsigma = vertex[8];
00412 
00413   deltah = hist_max - hist_min;
00414   deltam = wmean - gmean;
00415 
00416 
00417   /*----- Apply constraints? -----*/
00418   if ((b < 0.05) || (b > 1.5))          return (BIG_NUMBER);
00419   if ((g < 0.05) || (g > 1.0))          return (BIG_NUMBER);
00420   if ((w < 0.05) || (w > 1.0))          return (BIG_NUMBER);
00421   if ((b+g+w < 1.0) || (b+g+w > 2.0))   return (BIG_NUMBER);
00422 
00423   if ((bmean < hist_min) || (bmean > hist_max))  return (BIG_NUMBER);
00424   if ((gmean < hist_min) || (gmean > hist_max))  return (BIG_NUMBER);
00425   if ((wmean < hist_min) || (wmean > hist_max))  return (BIG_NUMBER);
00426   if ((gmean < bmean)    || (gmean > wmean))     return (BIG_NUMBER);
00427 
00428   if ((gmean-bmean) < 0.10*(wmean-bmean))        return (BIG_NUMBER);
00429   if ((wmean-gmean) < 0.10*(wmean-bmean))        return (BIG_NUMBER);
00430 
00431   if ((bsigma < 0.01*deltah) || (bsigma > 0.5*deltah))  return (BIG_NUMBER);
00432   if ((gsigma < 0.01*deltam) || (gsigma > 0.5*deltam))  return (BIG_NUMBER);
00433   if ((wsigma < 0.01*deltam) || (wsigma > 0.5*deltam))  return (BIG_NUMBER);
00434 
00435 
00436   /*----- Not constrained, so calculate actual error sum of squares -----*/
00437   sse = 0.0;
00438 
00439   for (i = hist_min;  i < hist_max;  i++)
00440     {
00441       t = i * delt;
00442       diff = ((float) hist_data[i])/numpts - estimate (vertex, t);
00443       sse += diff * diff;
00444     }
00445 
00446   
00447   return (sse);
00448 }
00449 
00450 
00451 /*---------------------------------------------------------------------------*/
00452 
00453 void allocate_arrays (float *** simplex, float ** centroid,
00454                       float ** response, float ** step_size, 
00455                       float ** test1, float ** test2)
00456 {
00457   int i;
00458 
00459   *centroid = (float *) malloc (sizeof(float) * DIMENSION);
00460   *response = (float *) malloc (sizeof(float) * (DIMENSION+1));
00461   *step_size = (float *) malloc (sizeof(float) * DIMENSION);
00462   *test1 = (float *) malloc (sizeof(float) * DIMENSION);
00463   *test2 = (float *) malloc (sizeof(float) * DIMENSION);
00464    
00465   *simplex = (float **) malloc (sizeof(float *) * (DIMENSION+1));
00466 
00467   for (i = 0;  i < DIMENSION+1;  i++)
00468     (*simplex)[i] = (float *) malloc (sizeof(float) * DIMENSION);
00469        
00470 }
00471 
00472 
00473 /*---------------------------------------------------------------------------*/
00474 
00475 void deallocate_arrays (float *** simplex, float ** centroid,
00476                         float ** response, float ** step_size, 
00477                         float ** test1, float ** test2)
00478 {
00479   int i;
00480 
00481   free (*centroid);    *centroid = NULL;
00482   free (*response);    *response = NULL;
00483   free (*step_size);   *step_size = NULL;
00484   free (*test1);       *test1 = NULL;
00485   free (*test2);       *test2 = NULL;
00486 
00487   for (i = 0;  i < DIMENSION+1;  i++)
00488     {
00489       free ((*simplex)[i]);
00490       (*simplex)[i] = NULL;
00491     }
00492 
00493   free (*simplex);     *simplex = NULL;
00494        
00495 }
00496 
00497 
00498 /*---------------------------------------------------------------------------*/
00499 
00500 void eval_vertices (float * response, int * worst, int * next, int * best)
00501 {
00502   int i;
00503 
00504   /* initialize values */
00505   *worst = 0;
00506   *best = 0;
00507 
00508   /* find the best and worst */
00509   for (i = 1;  i < DIMENSION+1;  i++)
00510     {
00511       if (response[i] > response[*worst])
00512         *worst = i;
00513       if (response[i] < response[*best])
00514         *best = i;
00515     }
00516 
00517   /* find the next worst index */
00518   if (*worst == 0)
00519     *next = 1;
00520   else
00521     *next = 0;
00522   
00523   for (i = 0;  i < DIMENSION+1;  i++)
00524     if ((i != *worst) && (response[i] > response[*next]))
00525       *next = i;
00526 }
00527 
00528 
00529 /*---------------------------------------------------------------------------*/
00530 
00531 void restart (float ** simplex, float * response, 
00532               float * step_size)
00533 {
00534   const float STEP_FACTOR = 0.9;
00535   int i, j;
00536   int worst, next, best;
00537   float minval, maxval;
00538 
00539 
00540   /* find the current best vertex */
00541   eval_vertices (response, &worst, &next, &best); 
00542 
00543   /* set the first vertex to the current best */
00544   for (i = 0; i < DIMENSION;  i++)
00545     simplex[0][i] = simplex[best][i];
00546 
00547   /* decrease step size */
00548   for (i = 0;  i < DIMENSION;  i++)
00549     step_size[i] *= STEP_FACTOR;
00550 
00551   /* set up remaining vertices of simplex using new step size */
00552   for (i = 1;  i < DIMENSION+1;  i++)
00553     for (j = 0;  j < DIMENSION;  j++)
00554       {
00555         minval = simplex[0][j] - step_size[j];
00556         maxval = simplex[0][j] + step_size[j];
00557       simplex[i][j] = rand_uniform (minval, maxval);
00558       }
00559 
00560   /* initialize response for each vector */
00561   for (i = 0;  i < DIMENSION+1;  i++)
00562     response[i] = calc_sse (simplex[i]);
00563 }
00564 
00565 
00566 /*---------------------------------------------------------------------------*/
00567 /*
00568   Calculate the centroid of the simplex, ignoring the worst vertex.
00569 */
00570 
00571 void calc_centroid (float ** simplex, int worst, float * centroid)
00572 {
00573   int i, j;
00574 
00575   for (i = 0;  i < DIMENSION;  i++)
00576     {
00577       centroid[i] = 0.0;
00578 
00579       /* add each vertex, except the worst */
00580       for (j = 0;  j < DIMENSION+1;  j++)
00581         if (j != worst)
00582           centroid[i] += simplex[j][i];
00583     }
00584 
00585   /* divide by the number of vertices */
00586   for (i = 0;  i < DIMENSION;  i++)
00587     centroid[i] /= DIMENSION;
00588 }
00589 
00590 
00591 /*---------------------------------------------------------------------------*/
00592 /*
00593   Calculate the reflection of the worst vertex about the centroid.
00594 */
00595 
00596 void calc_reflection (float ** simplex, float * centroid, 
00597                       int worst, float coef, float * vertex)
00598 {
00599   int i;
00600 
00601   for (i = 0;  i < DIMENSION;  i++)
00602     vertex[i] = centroid[i] + coef*(centroid[i] - simplex[worst][i]);
00603 }
00604 
00605 
00606 /*---------------------------------------------------------------------------*/
00607 /*
00608   Replace a vertex of the simplex.
00609 */
00610 
00611 void replace (float ** simplex, float * response, 
00612               int index, float * vertex, float resp)
00613 {
00614   int i;
00615 
00616   for (i = 0;  i < DIMENSION;  i++)
00617     simplex[index][i] = vertex[i];
00618 
00619   response[index] = resp;
00620 }
00621 
00622 
00623 /*---------------------------------------------------------------------------*/
00624 /*
00625   Calculate goodness of fit.
00626 */
00627 
00628 float calc_good_fit (float * response)
00629 {
00630   int i;
00631 
00632   float avg, sd, tmp;
00633 
00634   /* average the response values */
00635   avg = 0.0;
00636   for (i = 0;  i < DIMENSION+1;  i++)
00637     avg += response[i];
00638   avg /= DIMENSION+1;
00639 
00640   /* compute standard deviation of response */
00641   sd = 0.0;
00642   for (i = 0;  i < DIMENSION+1;  i++)
00643     {
00644       tmp = response[i] - avg;
00645       sd += tmp*tmp;
00646     }
00647   sd /= DIMENSION;
00648 
00649   return (sqrt(sd));
00650 }
00651 
00652 
00653 /*---------------------------------------------------------------------------*/
00654 /*
00655   Perform initialization for the Simplex algorithm.
00656 */
00657 
00658 void simplex_initialize (float * parameters, float ** simplex, 
00659                          float * response, float * step_size)
00660 {
00661   int i, j;
00662   int worst, next, best;
00663   float resp;
00664   float minval, maxval;
00665 
00666 
00667   for (j = 0;  j < DIMENSION;  j++)
00668     {
00669       simplex[0][j] = parameters[j];
00670       step_size[j] = 0.5 * parameters[j];
00671     }
00672 
00673   for (i = 1;  i < DIMENSION+1;  i++)
00674     for (j = 0;  j < DIMENSION;  j++)
00675       {
00676         minval = simplex[0][j] - step_size[j];
00677         maxval = simplex[0][j] + step_size[j];
00678         simplex[i][j] = rand_uniform (minval, maxval);
00679       }
00680 
00681   for (i = 0;  i < DIMENSION+1;  i++)
00682       response[i] = calc_sse (simplex[i]);
00683 
00684   for (i = 1;  i < 500;  i++)
00685     {
00686       for (j = 0;  j < DIMENSION;  j++)
00687         {
00688           minval = simplex[0][j] - step_size[j];
00689           maxval = simplex[0][j] + step_size[j];
00690           parameters[j] = rand_uniform (minval, maxval);
00691         }
00692 
00693       resp = calc_sse (parameters);
00694       eval_vertices (response, &worst, &next, &best);
00695       if (resp < response[worst])
00696         replace (simplex, response, worst, parameters, resp);
00697     }
00698 }
00699 
00700 
00701 /*---------------------------------------------------------------------------*/
00702 /*
00703   Use Simplex algorithm to estimate the voxel intensity distribution.
00704 
00705   The Simplex algorithm is adapted from: A Jump Start Course in C++ Programming
00706 */
00707 
00708 void simplex_optimization (float * parameters, float * sse)
00709 {
00710   const int MAX_ITERATIONS = 100;
00711   const int MAX_RESTARTS = 25;
00712   const float EXPANSION_COEF = 2.0;
00713   const float REFLECTION_COEF = 1.0;
00714   const float CONTRACTION_COEF = 0.5;
00715   const float TOLERANCE = 1.0e-10;
00716 
00717   float ** simplex   = NULL;
00718   float * centroid   = NULL;
00719   float * response   = NULL;
00720   float * step_size  = NULL;
00721   float * test1      = NULL;
00722   float * test2      = NULL;
00723   float resp1, resp2;
00724   int i, worst, best, next;
00725   int num_iter, num_restarts;
00726   int done;
00727   float fit;
00728 
00729 
00730   allocate_arrays (&simplex, &centroid, &response, &step_size, &test1, &test2);
00731   
00732   simplex_initialize (parameters, simplex, response, step_size);
00733 
00734   /* start loop to do simplex optimization */
00735   num_iter = 0;
00736   num_restarts = 0;
00737   done = 0;
00738   
00739   while (!done)
00740     {
00741       /* find the worst vertex and compute centroid of remaining simplex, 
00742          discarding the worst vertex */
00743       eval_vertices (response, &worst, &next, &best);
00744       calc_centroid (simplex, worst, centroid);
00745       
00746       /* reflect the worst point through the centroid */
00747       calc_reflection (simplex, centroid, worst, 
00748                        REFLECTION_COEF, test1);
00749       resp1 = calc_sse (test1);
00750 
00751       /* test the reflection against the best vertex and expand it if the
00752          reflection is better.  if not, keep the reflection */
00753       if (resp1 < response[best])
00754         {
00755           /* try expanding */
00756           calc_reflection (simplex, centroid, worst, EXPANSION_COEF, test2);
00757           resp2 = calc_sse (test2);
00758           if (resp2 <= resp1)    /* keep expansion */     
00759             replace (simplex, response, worst, test2, resp2);
00760           else                   /* keep reflection */
00761             replace (simplex, response, worst, test1, resp1);
00762         }
00763       else if (resp1 < response[next])
00764         {
00765           /* new response is between the best and next worst 
00766              so keep reflection */
00767           replace (simplex, response, worst, test1, resp1); 
00768         }
00769           else
00770         {
00771           /* try contraction */
00772           if (resp1 >= response[worst])
00773             calc_reflection (simplex, centroid, worst, 
00774                              -CONTRACTION_COEF, test2);
00775           else
00776             calc_reflection (simplex, centroid, worst, 
00777                              CONTRACTION_COEF, test2);
00778           resp2 =  calc_sse (test2);
00779           
00780           /* test the contracted response against the worst response */
00781           if (resp2 > response[worst])
00782             {
00783               /* new contracted response is worse, so decrease step size
00784                  and restart */
00785               num_iter = 0;
00786               num_restarts += 1;
00787               restart (simplex, response, step_size);
00788             }
00789           else       /* keep contraction */
00790             replace (simplex, response, worst, test2, resp2);
00791         }
00792 
00793       /* test to determine when to stop.  
00794          first, check the number of iterations */
00795       num_iter += 1;    /* increment iteration counter */
00796       if (num_iter >= MAX_ITERATIONS)
00797         {
00798           /* restart with smaller steps */
00799           num_iter = 0;
00800           num_restarts += 1;
00801           restart (simplex, response, step_size);
00802         }
00803 
00804       /* limit the number of restarts */
00805       if (num_restarts == MAX_RESTARTS)  done = 1;
00806 
00807       /* compare relative standard deviation of vertex responses 
00808          against a defined tolerance limit */
00809       fit = calc_good_fit (response);
00810       if (fit <= TOLERANCE)  done = 1;
00811 
00812       /* if done, copy the best solution to the output array */
00813       if (done) 
00814         {
00815           eval_vertices (response, &worst, &next, &best);
00816           for (i = 0;  i < DIMENSION;  i++)
00817             parameters[i] = simplex[best][i];
00818           *sse = response[best];
00819         }
00820 
00821     }  /* while (!done) */
00822  
00823   number_restarts = num_restarts;
00824   deallocate_arrays (&simplex, &centroid, &response, &step_size,
00825                      &test1, &test2);
00826 
00827 }
00828 
00829 
00830 /*---------------------------------------------------------------------------*/
00831 /*
00832   Write the parameter estimates.
00833 */
00834 
00835 void output_pdf_results (float * vertex, float sse)
00836 {
00837   float b, bmean, bsigma, g, gmean, gsigma, w, wmean, wsigma;
00838 
00839 
00840   /*----- Assign variables -----*/
00841   b      = vertex[0];
00842   bmean  = vertex[1];
00843   bsigma = vertex[2];
00844   g      = vertex[3];
00845   gmean  = vertex[4];
00846   gsigma = vertex[5];
00847   w      = vertex[6];
00848   wmean  = vertex[7];
00849   wsigma = vertex[8];
00850 
00851   printf ("hist_min = %d   hist_max = %d \n", hist_min, hist_max);
00852   printf ("gpeak    = %d   wpeak    = %d \n", gpeak, wpeak);
00853   printf ("rmse = %f \n", sqrt (sse / (hist_max - hist_min) ) );
00854   printf ("\nProbability Density Function Estimates: \n");
00855   printf ("Background Coef      = %f \n", b);
00856   printf ("Background Mean      = %f \n", bmean);
00857   printf ("Background Std Dev   = %f \n", bsigma);
00858   printf ("Gray Matter Coef     = %f \n", g);
00859   printf ("Gray Matter Mean     = %f \n", gmean);
00860   printf ("Gray Matter Std Dev  = %f \n", gsigma);
00861   printf ("White Matter Coef    = %f \n", w);
00862   printf ("White Matter Mean    = %f \n", wmean);
00863   printf ("White Matter Std Dev = %f \n", wsigma);
00864 }
00865 
00866 
00867 /*---------------------------------------------------------------------------*/
00868 /*
00869   Estimate the PDF of the voxel intensities.
00870 */
00871 
00872 Boolean estpdf (float * parameters)
00873 {
00874   float sse;
00875   Boolean ok;
00876 
00877 
00878   /*----- Progress report -----*/
00879   if (! quiet)  printf ("\nEstimating PDF of voxel intensities \n");
00880 
00881   
00882   /*----- Initialization for PDF estimation -----*/
00883   ok = estpdf_initialize ();
00884   if (! ok)  return (FALSE);
00885 
00886   generate_initial_guess (parameters);
00887  
00888 
00889   /*----- Get least squares estimate for PDF parameters -----*/
00890   simplex_optimization (parameters, &sse);
00891 
00892 
00893   /*----- Report PDF parameters -----*/
00894   if (! quiet)  output_pdf_results (parameters, sse);
00895 
00896 
00897   /*----- Free memory -----*/
00898   free (hist_data);    hist_data  = NULL;
00899   
00900 
00901   return (TRUE);
00902 }
00903 
00904 
00905 /*---------------------------------------------------------------------------*/