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   These routines estimate 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   File:    estpdf.c
00016   Author:  B. Douglas Ward
00017   Date:    27 May 1999
00018 
00019   Mod:     Extensive restructuring of the code.
00020   Date:    28 January 2000
00021 
00022   Mod:     USE_QUIET (RWCox)
00023   Date:    16 Apr 2003
00024 */
00025 
00026 #ifndef USE_QUIET  /* 16 Apr 2003 */
00027 # define quiet 0
00028 #endif
00029 
00030 
00031 /*---------------------------------------------------------------------------*/
00032 /*
00033   Include header files and forward declarations.
00034 */
00035 
00036 #include "pdf.h"
00037 
00038 float calc_error (float * vertex);
00039 
00040 
00041 /*---------------------------------------------------------------------------*/
00042 /*
00043   Global variables and constants.
00044 */
00045 
00046 #define DIMENSION 9      /* number of parameters to be estimated */
00047 
00048 pdf p;                   /* empirical pdf */
00049 
00050 
00051 /*---------------------------------------------------------------------------*/
00052 /*
00053   Include source code files.
00054 */
00055 
00056 #include "randgen.c"
00057 #include "pdf.c"
00058 #include "Simplexx.c"
00059 
00060 
00061 /*---------------------------------------------------------------------------*/
00062 /*
00063   Perform initialization for estimation of PDF from short array. 
00064 */
00065 
00066 void estpdf_short_initialize 
00067 (
00068   int nxyz, 
00069   short * sfim,
00070   float * gpeak,             /* estimated peak of gray-matter distribution */
00071   float * wpeak              /* estimated peak of white-matter distribution */
00072 )
00073 
00074 {
00075   pdf ps;
00076   int gmax, wmax;
00077   int kk;
00078   int ok = 1;
00079 
00080 
00081   /*---- Initialize pdf's -----*/
00082   PDF_initialize (&p);
00083   PDF_initialize (&ps);
00084 
00085 
00086   /*----- Convert short array to pdf estimate -----*/
00087   PDF_short_to_pdf (nxyz, sfim, &p);
00088   PDF_sprint ("\nOriginal PDF:", p);
00089 
00090 
00091   /*----- Trim extreme values from pdf estimate -----*/
00092   PDF_trim (0.01, 0.99, &p);
00093   PDF_sprint ("\nTrimmed PDF:", p);
00094 
00095 
00096   /*----- Smooth the pdf estimate -----*/
00097   PDF_copy (p, &ps);
00098   PDF_smooth (&ps);
00099   PDF_sprint ("\nSmoothed PDF:", ps);
00100 
00101 
00102   /*----- Try to locate bimodality of the pdf -----*/
00103   ok = PDF_find_bimodal (ps, &gmax, &wmax);
00104   if (ok)
00105     {
00106       *gpeak = PDF_ibin_to_xvalue (ps, gmax);
00107       *wpeak = PDF_ibin_to_xvalue (ps, wmax);
00108     }
00109   else
00110     {
00111       printf ("Unable to find bimodal distribution \n");
00112       *gpeak = (2.0/3.0)*p.lower_bnd + (1.0/3.0)*p.upper_bnd;
00113       *wpeak = (1.0/3.0)*p.lower_bnd + (2.0/3.0)*p.upper_bnd;
00114     }
00115 
00116 
00117   if( !quiet ){
00118    printf ("\nInitial PDF estimates: \n");
00119    printf ("Lower Bnd = %8.3f   Upper Bnd  = %8.3f \n", 
00120            p.lower_bnd, p.upper_bnd);
00121    printf ("Gray Peak = %8.3f   White Peak = %8.3f \n", *gpeak, *wpeak);
00122   }
00123 
00124 
00125   PDF_destroy (&ps);
00126 
00127 }
00128 
00129 
00130 /*---------------------------------------------------------------------------*/
00131 /*
00132   Perform initialization for estimation of PDF from float array. 
00133 */
00134 
00135 void estpdf_float_initialize 
00136 (
00137   int nxyz, 
00138   float * ffim,
00139   int nbin,
00140   float * gpeak,             /* estimated peak of gray-matter distribution */
00141   float * wpeak              /* estimated peak of white-matter distribution */
00142 )
00143 
00144 {
00145   pdf ps;
00146   int gmax, wmax;
00147   int kk;
00148   int ok = 1;
00149 
00150 
00151   /*---- Initialize pdf's -----*/
00152   PDF_initialize (&p);
00153   PDF_initialize (&ps);
00154 
00155 
00156   /*----- Convert float array to pdf estimate -----*/
00157   PDF_float_to_pdf (nxyz, ffim, nbin, &p);
00158   PDF_sprint ("\nOriginal PDF:", p);
00159 
00160 
00161   /*----- Trim extreme values from pdf estimate -----*/
00162   PDF_trim (0.01, 0.99, &p);
00163   PDF_sprint ("\nTrimmed PDF:", p);
00164 
00165 
00166   /*----- Smooth the pdf estimate -----*/
00167   PDF_copy (p, &ps);
00168   PDF_smooth (&ps);
00169   PDF_sprint ("\nSmoothed PDF:", ps);
00170 
00171 
00172   /*----- Try to locate bimodality of the pdf -----*/
00173   ok = PDF_find_bimodal (ps, &gmax, &wmax);
00174   if (ok)
00175     {
00176       *gpeak = PDF_ibin_to_xvalue (ps, gmax);
00177       *wpeak = PDF_ibin_to_xvalue (ps, wmax);
00178     }
00179   else
00180     {
00181       printf ("Unable to find bimodal distribution \n");
00182       *gpeak = (2.0/3.0)*p.lower_bnd + (1.0/3.0)*p.upper_bnd;
00183       *wpeak = (1.0/3.0)*p.lower_bnd + (2.0/3.0)*p.upper_bnd;
00184     }
00185 
00186 
00187   if( !quiet ){
00188     printf ("\nInitial PDF estimates: \n");
00189     printf ("Lower Bnd = %8.3f   Upper Bnd  = %8.3f \n", 
00190             p.lower_bnd, p.upper_bnd);
00191     printf ("Gray Peak = %8.3f   White Peak = %8.3f \n", *gpeak, *wpeak);
00192   }
00193 
00194 
00195   PDF_destroy (&ps);
00196 
00197 }
00198 
00199 
00200 /*---------------------------------------------------------------------------*/
00201 /*
00202   Generate the initial guess for the parameter vector.
00203 */
00204 
00205 void generate_initial_guess (float gpeak, float wpeak, float * parameters)
00206 {
00207   float b;                   /* coefficient for background distribution */
00208   float bmean;               /* mean for background distribution */ 
00209   float bsigma;              /* std. dev. for background distribution */
00210   float g;                   /* coefficient for gray-matter distribution */
00211   float gmean;               /* mean for gray-matter distribution */
00212   float gsigma;              /* std. dev. for gray-matter distribution */
00213   float w;                   /* coefficient for white-matter distribution */
00214   float wmean;               /* mean for white-matter distribution */
00215   float wsigma;              /* std. dev. for white-matter distribution */
00216 
00217 
00218   /*----- Initialize distribution coefficients -----*/
00219   b = 0.75;
00220   g = 0.25;
00221   w = 0.25;
00222 
00223 
00224   /*----- Initialize distribution means -----*/
00225   bmean = p.lower_bnd;
00226 
00227   if ((gpeak > p.lower_bnd) && (gpeak < p.upper_bnd) && (gpeak < wpeak)) 
00228     gmean = gpeak;
00229   else
00230     gmean = p.lower_bnd;
00231 
00232   if ((wpeak > p.lower_bnd) && (wpeak < p.upper_bnd) && (wpeak > gpeak))
00233     wmean = wpeak;
00234   else
00235     wmean = p.upper_bnd;
00236 
00237   if ((gmean-bmean) < 0.25*(wmean-bmean))  gmean = bmean + 0.25*(wmean-bmean);
00238   if ((wmean-gmean) < 0.25*(wmean-bmean))  gmean = wmean - 0.25*(wmean-bmean);
00239 
00240 
00241   /*----- Initialize distribution standard deviations -----*/
00242   bsigma = 0.25 * (p.upper_bnd - p.lower_bnd);
00243   gsigma = 0.25 * (wmean - gmean);
00244   wsigma = 0.25 * (wmean - gmean);
00245 
00246 
00247   /*----- Set parameter vector -----*/
00248   parameters[0] = b;
00249   parameters[1] = bmean;
00250   parameters[2] = bsigma;
00251   parameters[3] = g;
00252   parameters[4] = gmean;
00253   parameters[5] = gsigma;
00254   parameters[6] = w;
00255   parameters[7] = wmean;
00256   parameters[8] = wsigma;
00257 
00258 }
00259 
00260 
00261 /*---------------------------------------------------------------------------*/
00262 /*
00263   Write parameter vector.
00264 */
00265 
00266 void write_parameter_vector (float * parameters)
00267 {
00268   int i;
00269 
00270   printf ("Dimension = %d \n", DIMENSION);
00271   for (i = 0;  i < DIMENSION;  i++)
00272     printf ("parameter[%d] = %f \n", i, parameters[i]);
00273 }
00274 
00275 
00276 /*---------------------------------------------------------------------------*/
00277 /*
00278   Normal probability density function.
00279 */
00280 
00281 float normal (float x, float mean, float sigma)
00282 
00283 {
00284   float z;
00285 
00286   z = (x - mean) / sigma;
00287 
00288   return ( (1.0/(sqrt(2.0*PI)*sigma)) * exp (-0.5 * z * z) );
00289 }
00290 
00291 
00292 /*---------------------------------------------------------------------------*/
00293 /*
00294   Estimate the voxel intensity distribution as the sum of three normal PDF's.
00295 */
00296 
00297 float estimate (float * parameters, float x)
00298 
00299 {
00300   float b, bmean, bsigma, g, gmean, gsigma, w, wmean, wsigma;
00301   float z, fval;
00302 
00303 
00304   /*----- Initialize local variables -----*/
00305   b      = parameters[0];
00306   bmean  = parameters[1];
00307   bsigma = parameters[2];
00308   g      = parameters[3];
00309   gmean  = parameters[4];
00310   gsigma = parameters[5];
00311   w      = parameters[6];
00312   wmean  = parameters[7];
00313   wsigma = parameters[8];
00314 
00315 
00316   /*----- Calculate the sum of three normal PDF's -----*/ 
00317   fval  = b * normal (x, bmean, bsigma);
00318   fval += g * normal (x, gmean, gsigma);
00319   fval += w * normal (x, wmean, wsigma);
00320 
00321 
00322   return (fval);
00323   
00324 }
00325 
00326 
00327 /*---------------------------------------------------------------------------*/
00328 /*
00329   Calculate the error sum of squares for the PDF estimate.
00330 */
00331 
00332 float calc_error (float * vertex)
00333 
00334 {
00335   const float BIG_NUMBER = 1.0e+10;  /* return when constraints are violated */
00336 
00337   float b, bmean, bsigma, g, gmean, gsigma, w, wmean, wsigma;  /* parameters */
00338   float deltah, deltam;          /* rough estimate of spread of distribution */
00339 
00340   int i;
00341   float t;
00342   float diff, sse;
00343 
00344   count += 1;
00345 
00346 
00347   /*----- Assign local variables -----*/
00348   b      = vertex[0];
00349   bmean  = vertex[1];
00350   bsigma = vertex[2];
00351   g      = vertex[3];
00352   gmean  = vertex[4];
00353   gsigma = vertex[5];
00354   w      = vertex[6];
00355   wmean  = vertex[7];
00356   wsigma = vertex[8];
00357 
00358   deltah = p.upper_bnd - p.lower_bnd;
00359   deltam = wmean - gmean;
00360 
00361 
00362   /*----- Apply constraints? -----*/
00363   if ((b < 0.05) || (b > 1.5))          return (BIG_NUMBER);
00364   if ((g < 0.05) || (g > 1.0))          return (BIG_NUMBER);
00365   if ((w < 0.05) || (w > 1.0))          return (BIG_NUMBER);
00366   if ((b+g+w < 1.0) || (b+g+w > 2.0))   return (BIG_NUMBER);
00367 
00368   if ((bmean < p.lower_bnd) || (bmean > p.upper_bnd))  return (BIG_NUMBER);
00369   if ((gmean < p.lower_bnd) || (gmean > p.upper_bnd))  return (BIG_NUMBER);
00370   if ((wmean < p.lower_bnd) || (wmean > p.upper_bnd))  return (BIG_NUMBER);
00371   if ((gmean < bmean)    || (gmean > wmean))     return (BIG_NUMBER);
00372 
00373   if ((gmean-bmean) < 0.10*(wmean-bmean))        return (BIG_NUMBER);
00374   if ((wmean-gmean) < 0.10*(wmean-bmean))        return (BIG_NUMBER);
00375 
00376   if ((bsigma < 0.01*deltah) || (bsigma > 0.5*deltah))  return (BIG_NUMBER);
00377   if ((gsigma < 0.01*deltam) || (gsigma > 0.5*deltam))  return (BIG_NUMBER);
00378   if ((wsigma < 0.01*deltam) || (wsigma > 0.5*deltam))  return (BIG_NUMBER);
00379 
00380 
00381   /*----- Not constrained, so calculate actual error sum of squares -----*/
00382   sse = 0.0;
00383 
00384   for (i = 0;  i < p.nbin;  i++)
00385     {
00386       t = PDF_ibin_to_xvalue (p, i);
00387       diff = p.prob[i] - estimate (vertex, t)*p.width;
00388       sse += diff * diff;
00389     }
00390 
00391   
00392   return (sse);
00393 }
00394 
00395 
00396 /*---------------------------------------------------------------------------*/
00397 /*
00398   Write the parameter estimates.
00399 */
00400 
00401 void output_pdf_results (float * vertex, float sse)
00402 {
00403   float b, bmean, bsigma, g, gmean, gsigma, w, wmean, wsigma;
00404 
00405 
00406   /*----- Assign variables -----*/
00407   b      = vertex[0];
00408   bmean  = vertex[1];
00409   bsigma = vertex[2];
00410   g      = vertex[3];
00411   gmean  = vertex[4];
00412   gsigma = vertex[5];
00413   w      = vertex[6];
00414   wmean  = vertex[7];
00415   wsigma = vertex[8];
00416 
00417   if( !quiet ){
00418    printf ("\nProbability Density Function Estimates: \n");
00419    printf ("Background Coef      = %f \n", b);
00420    printf ("Background Mean      = %f \n", bmean);
00421    printf ("Background Std Dev   = %f \n", bsigma);
00422    printf ("Gray Matter Coef     = %f \n", g);
00423    printf ("Gray Matter Mean     = %f \n", gmean);
00424    printf ("Gray Matter Std Dev  = %f \n", gsigma);
00425    printf ("White Matter Coef    = %f \n", w);
00426    printf ("White Matter Mean    = %f \n", wmean);
00427    printf ("White Matter Std Dev = %f \n", wsigma);
00428    printf ("\nrmse = %f \n", sqrt (sse / p.nbin ));
00429   }
00430 
00431 }
00432 
00433 
00434 /*---------------------------------------------------------------------------*/
00435 /*
00436   Estimate the PDF of the voxel intensities.
00437 */
00438 
00439 void estpdf_short (int nxyz, short * sfim, float * parameters)
00440 {
00441   float gpeak;               /* estimated peak of gray-matter distribution */
00442   float wpeak;               /* estimated peak of white-matter distribution */
00443   float sse;
00444 
00445 
00446   /*----- Progress report -----*/
00447   if( !quiet )
00448    printf ("\nEstimating PDF of voxel intensities \n");
00449 
00450   
00451   /*----- Initialization for PDF estimation -----*/
00452   estpdf_short_initialize (nxyz, sfim, &gpeak, &wpeak);
00453 
00454 
00455   generate_initial_guess (gpeak, wpeak, parameters);
00456  
00457 
00458   /*----- Get least squares estimate for PDF parameters -----*/
00459   simplex_optimization (parameters, &sse);
00460 
00461 
00462   /*----- Report PDF parameters -----*/
00463   output_pdf_results (parameters, sse);
00464 
00465 
00466   /*----- Free memory -----*/
00467   /*
00468   PDF_destroy (&p);
00469   */
00470 
00471   return;
00472 }
00473 
00474 
00475 /*---------------------------------------------------------------------------*/
00476 /*
00477   Estimate the PDF of the voxel intensities.
00478 */
00479 
00480 void estpdf_float (int nxyz, float * ffim, int nbin, float * parameters)
00481 {
00482   float gpeak;               /* estimated peak of gray-matter distribution */
00483   float wpeak;               /* estimated peak of white-matter distribution */
00484   float sse;
00485 
00486 
00487   /*----- Progress report -----*/
00488   if( !quiet )
00489    printf ("\nEstimating PDF of voxel intensities \n");
00490 
00491   
00492   /*----- Initialization for PDF estimation -----*/
00493   estpdf_float_initialize (nxyz, ffim, nbin, &gpeak, &wpeak);
00494 
00495 
00496   /*----- Make initial estimate of the parameters from previous results -----*/
00497   generate_initial_guess (gpeak, wpeak, parameters);
00498  
00499 
00500   /*----- Get least squares estimate for PDF parameters -----*/
00501   simplex_optimization (parameters, &sse);
00502 
00503 
00504   /*----- Report PDF parameters -----*/
00505   output_pdf_results (parameters, sse);
00506 
00507 
00508   /*----- Free memory -----*/
00509   /*
00510   PDF_destroy (&p);
00511   */
00512  
00513   return ;
00514 }
00515 
00516 
00517 /*---------------------------------------------------------------------------*/