1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*---------------------------------------------------------------------------*/
8 /*
9 These routines estimate the probability density function (PDF)
10 corresponding to the distribution of gray-matter and white-matter
11 voxel intensities. The estimate is formed as the sum of three normal
12 distributions, using the Simplex algorithm for non-linear optimization
13 to estimate the unknown parameters.
14
15 File: estpdf.c
16 Author: B. Douglas Ward
17 Date: 27 May 1999
18
19 Mod: Extensive restructuring of the code.
20 Date: 28 January 2000
21
22 Mod: USE_QUIET (RWCox)
23 Date: 16 Apr 2003
24 */
25
26 #ifndef USE_QUIET /* 16 Apr 2003 */
27 # define quiet 0
28 #endif
29
30
31 /*---------------------------------------------------------------------------*/
32 /*
33 Include header files and forward declarations.
34 */
35
36 #include "pdf.h"
37
38 float calc_error (float * vertex);
39
40
41 /*---------------------------------------------------------------------------*/
42 /*
43 Global variables and constants.
44 */
45
46 #define DIMENSION 9 /* number of parameters to be estimated */
47
48 pdf p; /* empirical pdf */
49
50
51 /*---------------------------------------------------------------------------*/
52 /*
53 Include source code files.
54 */
55
56 #include "randgen.c"
57 #include "pdf.c"
58 #include "Simplexx.c"
59
60
61 /*---------------------------------------------------------------------------*/
62 /*
63 Perform initialization for estimation of PDF from short array.
64 */
65
estpdf_short_initialize(int nxyz,short * sfim,float * gpeak,float * wpeak)66 void estpdf_short_initialize
67 (
68 int nxyz,
69 short * sfim,
70 float * gpeak, /* estimated peak of gray-matter distribution */
71 float * wpeak /* estimated peak of white-matter distribution */
72 )
73
74 {
75 pdf ps;
76 int gmax, wmax;
77 int kk;
78 int ok = 1;
79
80
81 /*---- Initialize pdf's -----*/
82 PDF_initialize (&p);
83 PDF_initialize (&ps);
84
85
86 /*----- Convert short array to pdf estimate -----*/
87 PDF_short_to_pdf (nxyz, sfim, &p);
88 PDF_sprint ("\nOriginal PDF:", p);
89
90
91 /*----- Trim extreme values from pdf estimate -----*/
92 PDF_trim (0.01, 0.99, &p);
93 PDF_sprint ("\nTrimmed PDF:", p);
94
95
96 /*----- Smooth the pdf estimate -----*/
97 PDF_copy (p, &ps);
98 PDF_smooth (&ps);
99 PDF_sprint ("\nSmoothed PDF:", ps);
100
101
102 /*----- Try to locate bimodality of the pdf -----*/
103 ok = PDF_find_bimodal (ps, &gmax, &wmax);
104 if (ok)
105 {
106 *gpeak = PDF_ibin_to_xvalue (ps, gmax);
107 *wpeak = PDF_ibin_to_xvalue (ps, wmax);
108 }
109 else
110 {
111 printf ("Unable to find bimodal distribution \n");
112 *gpeak = (2.0/3.0)*p.lower_bnd + (1.0/3.0)*p.upper_bnd;
113 *wpeak = (1.0/3.0)*p.lower_bnd + (2.0/3.0)*p.upper_bnd;
114 }
115
116
117 if( !quiet ){
118 printf ("\nInitial PDF estimates: \n");
119 printf ("Lower Bnd = %8.3f Upper Bnd = %8.3f \n",
120 p.lower_bnd, p.upper_bnd);
121 printf ("Gray Peak = %8.3f White Peak = %8.3f \n", *gpeak, *wpeak);
122 }
123
124
125 PDF_destroy (&ps);
126
127 }
128
129
130 /*---------------------------------------------------------------------------*/
131 /*
132 Perform initialization for estimation of PDF from float array.
133 */
134
estpdf_float_initialize(int nxyz,float * ffim,int nbin,float * gpeak,float * wpeak)135 void estpdf_float_initialize
136 (
137 int nxyz,
138 float * ffim,
139 int nbin,
140 float * gpeak, /* estimated peak of gray-matter distribution */
141 float * wpeak /* estimated peak of white-matter distribution */
142 )
143
144 {
145 pdf ps;
146 int gmax, wmax;
147 int kk;
148 int ok = 1;
149
150
151 /*---- Initialize pdf's -----*/
152 PDF_initialize (&p);
153 PDF_initialize (&ps);
154
155
156 /*----- Convert float array to pdf estimate -----*/
157 PDF_float_to_pdf (nxyz, ffim, nbin, &p);
158 PDF_sprint ("\nOriginal PDF:", p);
159
160
161 /*----- Trim extreme values from pdf estimate -----*/
162 PDF_trim (0.01, 0.99, &p);
163 PDF_sprint ("\nTrimmed PDF:", p);
164
165
166 /*----- Smooth the pdf estimate -----*/
167 PDF_copy (p, &ps);
168 PDF_smooth (&ps);
169 PDF_sprint ("\nSmoothed PDF:", ps);
170
171
172 /*----- Try to locate bimodality of the pdf -----*/
173 ok = PDF_find_bimodal (ps, &gmax, &wmax);
174 if (ok)
175 {
176 *gpeak = PDF_ibin_to_xvalue (ps, gmax);
177 *wpeak = PDF_ibin_to_xvalue (ps, wmax);
178 }
179 else
180 {
181 printf ("Unable to find bimodal distribution \n");
182 *gpeak = (2.0/3.0)*p.lower_bnd + (1.0/3.0)*p.upper_bnd;
183 *wpeak = (1.0/3.0)*p.lower_bnd + (2.0/3.0)*p.upper_bnd;
184 }
185
186
187 if( !quiet ){
188 printf ("\nInitial PDF estimates: \n");
189 printf ("Lower Bnd = %8.3f Upper Bnd = %8.3f \n",
190 p.lower_bnd, p.upper_bnd);
191 printf ("Gray Peak = %8.3f White Peak = %8.3f \n", *gpeak, *wpeak);
192 }
193
194
195 PDF_destroy (&ps);
196
197 }
198
199
200 /*---------------------------------------------------------------------------*/
201 /*
202 Generate the initial guess for the parameter vector.
203 */
204
generate_initial_guess(float gpeak,float wpeak,float * parameters)205 void generate_initial_guess (float gpeak, float wpeak, float * parameters)
206 {
207 float b; /* coefficient for background distribution */
208 float bmean; /* mean for background distribution */
209 float bsigma; /* std. dev. for background distribution */
210 float g; /* coefficient for gray-matter distribution */
211 float gmean; /* mean for gray-matter distribution */
212 float gsigma; /* std. dev. for gray-matter distribution */
213 float w; /* coefficient for white-matter distribution */
214 float wmean; /* mean for white-matter distribution */
215 float wsigma; /* std. dev. for white-matter distribution */
216
217
218 /*----- Initialize distribution coefficients -----*/
219 b = 0.75;
220 g = 0.25;
221 w = 0.25;
222
223
224 /*----- Initialize distribution means -----*/
225 bmean = p.lower_bnd;
226
227 if ((gpeak > p.lower_bnd) && (gpeak < p.upper_bnd) && (gpeak < wpeak))
228 gmean = gpeak;
229 else
230 gmean = p.lower_bnd;
231
232 if ((wpeak > p.lower_bnd) && (wpeak < p.upper_bnd) && (wpeak > gpeak))
233 wmean = wpeak;
234 else
235 wmean = p.upper_bnd;
236
237 if ((gmean-bmean) < 0.25*(wmean-bmean)) gmean = bmean + 0.25*(wmean-bmean);
238 if ((wmean-gmean) < 0.25*(wmean-bmean)) gmean = wmean - 0.25*(wmean-bmean);
239
240
241 /*----- Initialize distribution standard deviations -----*/
242 bsigma = 0.25 * (p.upper_bnd - p.lower_bnd);
243 gsigma = 0.25 * (wmean - gmean);
244 wsigma = 0.25 * (wmean - gmean);
245
246
247 /*----- Set parameter vector -----*/
248 parameters[0] = b;
249 parameters[1] = bmean;
250 parameters[2] = bsigma;
251 parameters[3] = g;
252 parameters[4] = gmean;
253 parameters[5] = gsigma;
254 parameters[6] = w;
255 parameters[7] = wmean;
256 parameters[8] = wsigma;
257
258 }
259
260
261 /*---------------------------------------------------------------------------*/
262 /*
263 Write parameter vector.
264 */
265
write_parameter_vector(float * parameters)266 void write_parameter_vector (float * parameters)
267 {
268 int i;
269
270 printf ("Dimension = %d \n", DIMENSION);
271 for (i = 0; i < DIMENSION; i++)
272 printf ("parameter[%d] = %f \n", i, parameters[i]);
273 }
274
275
276 /*---------------------------------------------------------------------------*/
277 /*
278 Normal probability density function.
279 */
280
normal(float x,float mean,float sigma)281 float normal (float x, float mean, float sigma)
282
283 {
284 float z;
285
286 z = (x - mean) / sigma;
287
288 return ( (1.0/(sqrt(2.0*PI)*sigma)) * exp (-0.5 * z * z) );
289 }
290
291
292 /*---------------------------------------------------------------------------*/
293 /*
294 Estimate the voxel intensity distribution as the sum of three normal PDF's.
295 */
296
estimate(float * parameters,float x)297 float estimate (float * parameters, float x)
298
299 {
300 float b, bmean, bsigma, g, gmean, gsigma, w, wmean, wsigma;
301 float z, fval;
302
303
304 /*----- Initialize local variables -----*/
305 b = parameters[0];
306 bmean = parameters[1];
307 bsigma = parameters[2];
308 g = parameters[3];
309 gmean = parameters[4];
310 gsigma = parameters[5];
311 w = parameters[6];
312 wmean = parameters[7];
313 wsigma = parameters[8];
314
315
316 /*----- Calculate the sum of three normal PDF's -----*/
317 fval = b * normal (x, bmean, bsigma);
318 fval += g * normal (x, gmean, gsigma);
319 fval += w * normal (x, wmean, wsigma);
320
321
322 return (fval);
323
324 }
325
326
327 /*---------------------------------------------------------------------------*/
328 /*
329 Calculate the error sum of squares for the PDF estimate.
330 */
331
calc_error(float * vertex)332 float calc_error (float * vertex)
333
334 {
335 const float BIG_NUMBER = 1.0e+10; /* return when constraints are violated */
336
337 float b, bmean, bsigma, g, gmean, gsigma, w, wmean, wsigma; /* parameters */
338 float deltah, deltam; /* rough estimate of spread of distribution */
339
340 int i;
341 float t;
342 float diff, sse;
343
344 count += 1;
345
346
347 /*----- Assign local variables -----*/
348 b = vertex[0];
349 bmean = vertex[1];
350 bsigma = vertex[2];
351 g = vertex[3];
352 gmean = vertex[4];
353 gsigma = vertex[5];
354 w = vertex[6];
355 wmean = vertex[7];
356 wsigma = vertex[8];
357
358 deltah = p.upper_bnd - p.lower_bnd;
359 deltam = wmean - gmean;
360
361
362 /*----- Apply constraints? -----*/
363 if ((b < 0.05) || (b > 1.5)) return (BIG_NUMBER);
364 if ((g < 0.05) || (g > 1.0)) return (BIG_NUMBER);
365 if ((w < 0.05) || (w > 1.0)) return (BIG_NUMBER);
366 if ((b+g+w < 1.0) || (b+g+w > 2.0)) return (BIG_NUMBER);
367
368 if ((bmean < p.lower_bnd) || (bmean > p.upper_bnd)) return (BIG_NUMBER);
369 if ((gmean < p.lower_bnd) || (gmean > p.upper_bnd)) return (BIG_NUMBER);
370 if ((wmean < p.lower_bnd) || (wmean > p.upper_bnd)) return (BIG_NUMBER);
371 if ((gmean < bmean) || (gmean > wmean)) return (BIG_NUMBER);
372
373 if ((gmean-bmean) < 0.10*(wmean-bmean)) return (BIG_NUMBER);
374 if ((wmean-gmean) < 0.10*(wmean-bmean)) return (BIG_NUMBER);
375
376 if ((bsigma < 0.01*deltah) || (bsigma > 0.5*deltah)) return (BIG_NUMBER);
377 if ((gsigma < 0.01*deltam) || (gsigma > 0.5*deltam)) return (BIG_NUMBER);
378 if ((wsigma < 0.01*deltam) || (wsigma > 0.5*deltam)) return (BIG_NUMBER);
379
380
381 /*----- Not constrained, so calculate actual error sum of squares -----*/
382 sse = 0.0;
383
384 for (i = 0; i < p.nbin; i++)
385 {
386 t = PDF_ibin_to_xvalue (p, i);
387 diff = p.prob[i] - estimate (vertex, t)*p.width;
388 sse += diff * diff;
389 }
390
391
392 return (sse);
393 }
394
395
396 /*---------------------------------------------------------------------------*/
397 /*
398 Write the parameter estimates.
399 */
400
output_pdf_results(float * vertex,float sse)401 void output_pdf_results (float * vertex, float sse)
402 {
403 float b, bmean, bsigma, g, gmean, gsigma, w, wmean, wsigma;
404
405
406 /*----- Assign variables -----*/
407 b = vertex[0];
408 bmean = vertex[1];
409 bsigma = vertex[2];
410 g = vertex[3];
411 gmean = vertex[4];
412 gsigma = vertex[5];
413 w = vertex[6];
414 wmean = vertex[7];
415 wsigma = vertex[8];
416
417 if( !quiet ){
418 printf ("\nProbability Density Function Estimates: \n");
419 printf ("Background Coef = %f \n", b);
420 printf ("Background Mean = %f \n", bmean);
421 printf ("Background Std Dev = %f \n", bsigma);
422 printf ("Gray Matter Coef = %f \n", g);
423 printf ("Gray Matter Mean = %f \n", gmean);
424 printf ("Gray Matter Std Dev = %f \n", gsigma);
425 printf ("White Matter Coef = %f \n", w);
426 printf ("White Matter Mean = %f \n", wmean);
427 printf ("White Matter Std Dev = %f \n", wsigma);
428 printf ("\nrmse = %f \n", sqrt (sse / p.nbin ));
429 }
430
431 }
432
433
434 /*---------------------------------------------------------------------------*/
435 /*
436 Estimate the PDF of the voxel intensities.
437 */
438
estpdf_short(int nxyz,short * sfim,float * parameters)439 void estpdf_short (int nxyz, short * sfim, float * parameters)
440 {
441 float gpeak; /* estimated peak of gray-matter distribution */
442 float wpeak; /* estimated peak of white-matter distribution */
443 float sse;
444
445
446 /*----- Progress report -----*/
447 if( !quiet )
448 printf ("\nEstimating PDF of voxel intensities \n");
449
450
451 /*----- Initialization for PDF estimation -----*/
452 estpdf_short_initialize (nxyz, sfim, &gpeak, &wpeak);
453
454
455 generate_initial_guess (gpeak, wpeak, parameters);
456
457
458 /*----- Get least squares estimate for PDF parameters -----*/
459 simplex_optimization (parameters, &sse);
460
461
462 /*----- Report PDF parameters -----*/
463 output_pdf_results (parameters, sse);
464
465
466 /*----- Free memory -----*/
467 /*
468 PDF_destroy (&p);
469 */
470
471 return;
472 }
473
474
475 /*---------------------------------------------------------------------------*/
476 /*
477 Estimate the PDF of the voxel intensities.
478 */
479
estpdf_float(int nxyz,float * ffim,int nbin,float * parameters)480 void estpdf_float (int nxyz, float * ffim, int nbin, float * parameters)
481 {
482 float gpeak; /* estimated peak of gray-matter distribution */
483 float wpeak; /* estimated peak of white-matter distribution */
484 float sse;
485
486
487 /*----- Progress report -----*/
488 if( !quiet )
489 printf ("\nEstimating PDF of voxel intensities \n");
490
491
492 /*----- Initialization for PDF estimation -----*/
493 estpdf_float_initialize (nxyz, ffim, nbin, &gpeak, &wpeak);
494
495
496 /*----- Make initial estimate of the parameters from previous results -----*/
497 generate_initial_guess (gpeak, wpeak, parameters);
498
499
500 /*----- Get least squares estimate for PDF parameters -----*/
501 simplex_optimization (parameters, &sse);
502
503
504 /*----- Report PDF parameters -----*/
505 output_pdf_results (parameters, sse);
506
507
508 /*----- Free memory -----*/
509 /*
510 PDF_destroy (&p);
511 */
512
513 return ;
514 }
515
516
517 /*---------------------------------------------------------------------------*/
518