00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef USE_QUIET
00027 # define quiet 0
00028 #endif
00029
00030
00031
00032
00033
00034
00035
00036 #include "pdf.h"
00037
00038 float calc_error (float * vertex);
00039
00040
00041
00042
00043
00044
00045
00046 #define DIMENSION 9
00047
00048 pdf p;
00049
00050
00051
00052
00053
00054
00055
00056 #include "randgen.c"
00057 #include "pdf.c"
00058 #include "Simplexx.c"
00059
00060
00061
00062
00063
00064
00065
00066 void estpdf_short_initialize
00067 (
00068 int nxyz,
00069 short * sfim,
00070 float * gpeak,
00071 float * wpeak
00072 )
00073
00074 {
00075 pdf ps;
00076 int gmax, wmax;
00077 int kk;
00078 int ok = 1;
00079
00080
00081
00082 PDF_initialize (&p);
00083 PDF_initialize (&ps);
00084
00085
00086
00087 PDF_short_to_pdf (nxyz, sfim, &p);
00088 PDF_sprint ("\nOriginal PDF:", p);
00089
00090
00091
00092 PDF_trim (0.01, 0.99, &p);
00093 PDF_sprint ("\nTrimmed PDF:", p);
00094
00095
00096
00097 PDF_copy (p, &ps);
00098 PDF_smooth (&ps);
00099 PDF_sprint ("\nSmoothed PDF:", ps);
00100
00101
00102
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
00133
00134
00135 void estpdf_float_initialize
00136 (
00137 int nxyz,
00138 float * ffim,
00139 int nbin,
00140 float * gpeak,
00141 float * wpeak
00142 )
00143
00144 {
00145 pdf ps;
00146 int gmax, wmax;
00147 int kk;
00148 int ok = 1;
00149
00150
00151
00152 PDF_initialize (&p);
00153 PDF_initialize (&ps);
00154
00155
00156
00157 PDF_float_to_pdf (nxyz, ffim, nbin, &p);
00158 PDF_sprint ("\nOriginal PDF:", p);
00159
00160
00161
00162 PDF_trim (0.01, 0.99, &p);
00163 PDF_sprint ("\nTrimmed PDF:", p);
00164
00165
00166
00167 PDF_copy (p, &ps);
00168 PDF_smooth (&ps);
00169 PDF_sprint ("\nSmoothed PDF:", ps);
00170
00171
00172
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
00203
00204
00205 void generate_initial_guess (float gpeak, float wpeak, float * parameters)
00206 {
00207 float b;
00208 float bmean;
00209 float bsigma;
00210 float g;
00211 float gmean;
00212 float gsigma;
00213 float w;
00214 float wmean;
00215 float wsigma;
00216
00217
00218
00219 b = 0.75;
00220 g = 0.25;
00221 w = 0.25;
00222
00223
00224
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
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
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
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
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
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
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
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
00330
00331
00332 float calc_error (float * vertex)
00333
00334 {
00335 const float BIG_NUMBER = 1.0e+10;
00336
00337 float b, bmean, bsigma, g, gmean, gsigma, w, wmean, wsigma;
00338 float deltah, deltam;
00339
00340 int i;
00341 float t;
00342 float diff, sse;
00343
00344 count += 1;
00345
00346
00347
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
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
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
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
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
00437
00438
00439 void estpdf_short (int nxyz, short * sfim, float * parameters)
00440 {
00441 float gpeak;
00442 float wpeak;
00443 float sse;
00444
00445
00446
00447 if( !quiet )
00448 printf ("\nEstimating PDF of voxel intensities \n");
00449
00450
00451
00452 estpdf_short_initialize (nxyz, sfim, &gpeak, &wpeak);
00453
00454
00455 generate_initial_guess (gpeak, wpeak, parameters);
00456
00457
00458
00459 simplex_optimization (parameters, &sse);
00460
00461
00462
00463 output_pdf_results (parameters, sse);
00464
00465
00466
00467
00468
00469
00470
00471 return;
00472 }
00473
00474
00475
00476
00477
00478
00479
00480 void estpdf_float (int nxyz, float * ffim, int nbin, float * parameters)
00481 {
00482 float gpeak;
00483 float wpeak;
00484 float sse;
00485
00486
00487
00488 if( !quiet )
00489 printf ("\nEstimating PDF of voxel intensities \n");
00490
00491
00492
00493 estpdf_float_initialize (nxyz, ffim, nbin, &gpeak, &wpeak);
00494
00495
00496
00497 generate_initial_guess (gpeak, wpeak, parameters);
00498
00499
00500
00501 simplex_optimization (parameters, &sse);
00502
00503
00504
00505 output_pdf_results (parameters, sse);
00506
00507
00508
00509
00510
00511
00512
00513 return ;
00514 }
00515
00516
00517