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 This program estimates 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 The Simplex algorithm is adapted from: A Jump Start Course in C++ Programming
16
17 File: estpdf.c
18 Author: B. Douglas Ward
19 Date: 27 May 1999
20 */
21
22
23 /*---------------------------------------------------------------------------*/
24 /*
25 Include source code.
26 */
27
28 #include "randgen.c"
29
30 /*---------------------------------------------------------------------------*/
31 /*
32 Global variables and constants.
33 */
34
35 #define DIMENSION 9 /* number of parameters to be estimated */
36
37 int * hist_data; /* histogram of voxel gray-scale intensities */
38 int hist_min; /* minimum histogram bin value; this is set to
39 exclude the large population near zero */
40 int hist_max; /* maximum histogram bin value; this is set to
41 exclude voxels with very high intensities */
42 int gpeak; /* estimated peak of gray-matter distribution */
43 int wpeak; /* estimated peak of white-matter distribution */
44 int numpts; /* number of voxels within the histogram limits */
45
46 int count = 0; /* count of sse evaluations */
47 int number_restarts = 0; /* count of simplex algorithm restarts */
48
49
50 /*---------------------------------------------------------------------------*/
51 /*
52 Set up the histogram and generate some of the initial parameter estimates.
53
54 This routine was adapted from: program 3dstrip.c
55 */
56
57 #undef MAX
58 #define MAX(a,b) (((a)<(b)) ? (b) : (a))
59
60 #undef MIN
61 #define MIN(a,b) (((a)>(b)) ? (b) : (a))
62
63 #define NPMAX 128
64
65
find_base_value(int nxyz,short * sfim)66 void find_base_value( int nxyz , short * sfim )
67 {
68 float bper = 50.0 , bmin = 1 ;
69 float tper = 98.0;
70
71 int ii , kk , nbin , sval , sum , nbot , a,b,c , npeak,ntop , nvox ;
72 int * fbin ;
73 int kmin[NPMAX] , kmax[NPMAX] ;
74 int nmin , nmax ;
75
76 /*-- make histogram of shorts --*/
77
78 fbin = (int *) malloc( sizeof(int) * 32768 ) ;
79
80 for( kk=0 ; kk < 32768 ; kk++ ) fbin[kk] = 0 ;
81
82 nvox = 0 ;
83
84 for( ii=0 ; ii < nxyz ; ii++ ){
85 kk = sfim[ii] ; if( kk >= 0 ){ fbin[kk]++ ; nvox++ ; }
86 }
87
88 /*-- find largest value --*/
89
90 for( kk=32767 ; kk > 0 ; kk-- ) if( fbin[kk] > 0 ) break ;
91 if( kk == 0 ){fprintf(stderr,"** find_base_value: All voxels are zero!\n");exit(1);}
92 nbin = kk+1 ;
93
94 /*-- find bper point in cumulative distribution --*/
95
96 sval = 0.01 * bper * nvox ;
97 sum = 0 ;
98 for( kk=0 ; kk < nbin ; kk++ ){
99 sum += fbin[kk] ; if( sum >= sval ) break ;
100 }
101 nbot = kk ; if( nbot == 0 ) nbot = 1 ; if( bmin > nbot ) nbot = bmin ;
102 if( nbot >= nbin-9 ){
103 fprintf(stderr,"** find_base_value: Base point on histogram too high\n");
104 exit(1);
105 }
106
107 hist_min = nbot;
108
109
110 /*-- find tper point in cumulative distribution --*/
111
112 sval = 0.01 * tper * nvox ;
113 sum = 0 ;
114 for( kk=0 ; kk < nbin ; kk++ ){
115 sum += fbin[kk] ; if( sum >= sval ) break ;
116 }
117
118 hist_max = kk ;
119
120
121 /*----- Save original histogram data -----*/
122 hist_data = (int *) malloc( sizeof(int) * nbin ) ;
123 for (kk = 0; kk < nbin; kk++)
124 hist_data[kk] = fbin[kk];
125
126
127 /*-- smooth histogram --*/
128
129 b = fbin[nbot-1] ; c = fbin[nbot] ;
130 for( kk=nbot ; kk < nbin ; kk++ ){
131 a = b ; b = c ; c = fbin[kk+1] ; fbin[kk] = 0.25*(a+c+2*b) ;
132 }
133
134 /*-- find minima and maxima above bper point --*/
135
136 nmin = nmax = 0 ;
137 for( kk=nbot+1 ; kk < nbin ; kk++ ){
138 if( fbin[kk] < fbin[kk-1] && fbin[kk] < fbin[kk+1] && nmin < NPMAX ){
139 kmin[nmin++] = kk ;
140 } else if( fbin[kk] > fbin[kk-1] && fbin[kk] > fbin[kk+1] && nmax < NPMAX ){
141 kmax[nmax++] = kk ;
142 }
143 }
144
145
146 /*-- find the largest two maxima --*/
147
148 if( nmax == 0 ){
149 fprintf(stderr,"** find_base_value: No histogram maxima above base point\n");
150 exit(1);
151 }
152
153 if( nmax == 1 ){
154 npeak = kmax[0] ; ntop = 0 ;
155 } else {
156 int f1,f2 , k1,k2 , fk , klow,kup ;
157
158 k1 = 0 ; f1 = fbin[kmax[0]] ;
159 k2 = 1 ; f2 = fbin[kmax[1]] ;
160 if( f1 < f2 ){
161 k1 = 1 ; f1 = fbin[kmax[1]] ;
162 k2 = 0 ; f2 = fbin[kmax[0]] ;
163 }
164
165 for( kk=2 ; kk < nmax ; kk++ ){
166 fk = fbin[kmax[kk]] ;
167 if( fk > f1 ){
168 f2 = f1 ; k2 = k1 ;
169 f1 = fk ; k1 = kk ;
170 } else if( fk > f2 ){
171 f2 = fk ; k2 = kk ;
172 }
173 }
174 npeak = MIN( kmax[k1] , kmax[k2] ) ; /* smaller bin of the 2 top peaks */
175
176 /* find valley between 2 peaks */
177
178 ntop = MAX( kmax[k1] , kmax[k2] ) ;
179
180 gpeak = npeak;
181 wpeak = ntop;
182
183
184 fk = fbin[ntop] ; klow = ntop ;
185 for( kk=ntop-1 ; kk >= npeak ; kk-- ){
186 if( fbin[kk] < fk ){ fk = fbin[kk] ; klow = kk ; }
187 }
188 fk = MAX( 0.10*fk , 0.05*fbin[ntop] ) ;
189 kup = MIN( nbin-1 , ntop+3*(ntop-klow+2) ) ;
190 for( kk=ntop+1 ; kk <= kup ; kk++ ) if( fbin[kk] < fk ) break ;
191
192 ntop = kk ;
193 }
194
195 for( kk=npeak-1 ; kk > 0 ; kk-- )
196 if( fbin[kk] < fbin[kk-1] && fbin[kk] < fbin[kk+1] ) break ;
197
198
199 if (kk > hist_min) hist_min = kk ;
200
201 if( ntop == 0 )
202 {
203 ntop = npeak + (npeak-kk) ;
204 gpeak = npeak;
205 wpeak = ntop;
206 }
207
208 free (fbin);
209
210 return ;
211 }
212
213
214 /*---------------------------------------------------------------------------*/
215 /*
216 Perform initialization for estimation of PDF.
217 */
218
219
estpdf_initialize()220 RwcBoolean estpdf_initialize ()
221
222 {
223 int nx, ny, nz, nxy, nxyz, ixyz; /* voxel counters */
224 int n; /* histogram bin index */
225 short * sfim; /* pointer to anat data */
226
227
228 /*----- Initialize local variables -----*/
229 if (anat == NULL) return (FALSE);
230 nx = DSET_NX(anat); ny = DSET_NY(anat); nz = DSET_NZ(anat);
231 nxy = nx*ny; nxyz = nxy*nz;
232 sfim = (short *) DSET_BRICK_ARRAY(anat,0) ;
233 if (sfim == NULL) return (FALSE);
234
235
236 /*----- Set up histogram, and get initial parameter estimates -----*/
237 find_base_value (nxyz, sfim);
238
239
240 /*----- Count number of voxels inside histogram intensity limits -----*/
241 numpts = 0;
242 for (n = hist_min; n < hist_max; n++)
243 numpts += hist_data[n];
244
245
246 /*----- Initialize random number generator -----*/
247 srand48 (1234567);
248
249
250 return (TRUE);
251 }
252
253
254 /*---------------------------------------------------------------------------*/
255 /*
256 Generate the initial guess for the parameter vector.
257 */
258
generate_initial_guess(float * parameters)259 void generate_initial_guess (float * parameters)
260 {
261 float b; /* coefficient for background distribution */
262 float bmean; /* mean for background distribution */
263 float bsigma; /* std. dev. for background distribution */
264 float g; /* coefficient for gray-matter distribution */
265 float gmean; /* mean for gray-matter distribution */
266 float gsigma; /* std. dev. for gray-matter distribution */
267 float w; /* coefficient for white-matter distribution */
268 float wmean; /* mean for white-matter distribution */
269 float wsigma; /* std. dev. for white-matter distribution */
270
271
272 /*----- Initialize distribution coefficients -----*/
273 b = 0.75;
274 g = 0.25;
275 w = 0.25;
276
277
278 /*----- Initialize distribution means -----*/
279 bmean = hist_min;
280
281 if ((gpeak > hist_min) && (gpeak < hist_max) && (gpeak < wpeak))
282 gmean = gpeak;
283 else
284 gmean = hist_min;
285
286 if ((wpeak > hist_min) && (wpeak < hist_max) && (wpeak > gpeak))
287 wmean = wpeak;
288 else
289 wmean = hist_max;
290
291 if ((gmean-bmean) < 0.25*(wmean-bmean)) gmean = bmean + 0.25*(wmean-bmean);
292 if ((wmean-gmean) < 0.25*(wmean-bmean)) gmean = wmean - 0.25*(wmean-bmean);
293
294
295 /*----- Initialize distribution standard deviations -----*/
296 bsigma = 0.25 * (hist_max - hist_min);
297 gsigma = 0.25 * (wmean - gmean);
298 wsigma = 0.25 * (wmean - gmean);
299
300
301 /*----- Set parameter vector -----*/
302 parameters[0] = b;
303 parameters[1] = bmean;
304 parameters[2] = bsigma;
305 parameters[3] = g;
306 parameters[4] = gmean;
307 parameters[5] = gsigma;
308 parameters[6] = w;
309 parameters[7] = wmean;
310 parameters[8] = wsigma;
311
312 }
313
314
315 /*---------------------------------------------------------------------------*/
316 /*
317 Write parameter vector.
318 */
319
write_parameter_vector(float * parameters)320 void write_parameter_vector (float * parameters)
321 {
322 int i;
323
324 printf ("Dimension = %d \n", DIMENSION);
325 for (i = 0; i < DIMENSION; i++)
326 printf ("parameter[%d] = %f \n", i, parameters[i]);
327 }
328
329
330 /*---------------------------------------------------------------------------*/
331 /*
332 Normal probability density function.
333 */
334
normal(float x,float mean,float sigma)335 float normal (float x, float mean, float sigma)
336
337 {
338 float z;
339
340 z = (x - mean) / sigma;
341
342 return ( (1.0/(sqrt(2.0*PI)*sigma)) * exp (-0.5 * z * z) );
343 }
344
345
346 /*---------------------------------------------------------------------------*/
347 /*
348 Estimate the voxel intensity distribution as the sum of three normal PDF's.
349 */
350
estimate(float * parameters,float x)351 float estimate (float * parameters, float x)
352
353 {
354 float b, bmean, bsigma, g, gmean, gsigma, w, wmean, wsigma;
355 float z, fval;
356
357
358 /*----- Initialize local variables -----*/
359 b = parameters[0];
360 bmean = parameters[1];
361 bsigma = parameters[2];
362 g = parameters[3];
363 gmean = parameters[4];
364 gsigma = parameters[5];
365 w = parameters[6];
366 wmean = parameters[7];
367 wsigma = parameters[8];
368
369
370 /*----- Calculate the sum of three normal PDF's -----*/
371 fval = b * normal (x, bmean, bsigma);
372 fval += g * normal (x, gmean, gsigma);
373 fval += w * normal (x, wmean, wsigma);
374
375
376 return (fval);
377
378 }
379
380
381 /*---------------------------------------------------------------------------*/
382 /*
383 Calculate the error sum of squares for the PDF estimate.
384 */
385
calc_sse(float * vertex)386 float calc_sse (float * vertex)
387
388 {
389 const float delt = 1.0; /* histogram bin size */
390 const float BIG_NUMBER = 1.0e+10; /* return when constraints are violated */
391
392 float b, bmean, bsigma, g, gmean, gsigma, w, wmean, wsigma; /* parameters */
393 float deltah, deltam; /* rough estimate of spread of distribution */
394
395 int i;
396 float t;
397 float diff, sse;
398
399 count += 1;
400
401
402 /*----- Assign local variables -----*/
403 b = vertex[0];
404 bmean = vertex[1];
405 bsigma = vertex[2];
406 g = vertex[3];
407 gmean = vertex[4];
408 gsigma = vertex[5];
409 w = vertex[6];
410 wmean = vertex[7];
411 wsigma = vertex[8];
412
413 deltah = hist_max - hist_min;
414 deltam = wmean - gmean;
415
416
417 /*----- Apply constraints? -----*/
418 if ((b < 0.05) || (b > 1.5)) return (BIG_NUMBER);
419 if ((g < 0.05) || (g > 1.0)) return (BIG_NUMBER);
420 if ((w < 0.05) || (w > 1.0)) return (BIG_NUMBER);
421 if ((b+g+w < 1.0) || (b+g+w > 2.0)) return (BIG_NUMBER);
422
423 if ((bmean < hist_min) || (bmean > hist_max)) return (BIG_NUMBER);
424 if ((gmean < hist_min) || (gmean > hist_max)) return (BIG_NUMBER);
425 if ((wmean < hist_min) || (wmean > hist_max)) return (BIG_NUMBER);
426 if ((gmean < bmean) || (gmean > wmean)) return (BIG_NUMBER);
427
428 if ((gmean-bmean) < 0.10*(wmean-bmean)) return (BIG_NUMBER);
429 if ((wmean-gmean) < 0.10*(wmean-bmean)) return (BIG_NUMBER);
430
431 if ((bsigma < 0.01*deltah) || (bsigma > 0.5*deltah)) return (BIG_NUMBER);
432 if ((gsigma < 0.01*deltam) || (gsigma > 0.5*deltam)) return (BIG_NUMBER);
433 if ((wsigma < 0.01*deltam) || (wsigma > 0.5*deltam)) return (BIG_NUMBER);
434
435
436 /*----- Not constrained, so calculate actual error sum of squares -----*/
437 sse = 0.0;
438
439 for (i = hist_min; i < hist_max; i++)
440 {
441 t = i * delt;
442 diff = ((float) hist_data[i])/numpts - estimate (vertex, t);
443 sse += diff * diff;
444 }
445
446
447 return (sse);
448 }
449
450
451 /*---------------------------------------------------------------------------*/
452
allocate_arrays(float *** simplex,float ** centroid,float ** response,float ** step_size,float ** test1,float ** test2)453 void allocate_arrays (float *** simplex, float ** centroid,
454 float ** response, float ** step_size,
455 float ** test1, float ** test2)
456 {
457 int i;
458
459 *centroid = (float *) malloc (sizeof(float) * DIMENSION);
460 *response = (float *) malloc (sizeof(float) * (DIMENSION+1));
461 *step_size = (float *) malloc (sizeof(float) * DIMENSION);
462 *test1 = (float *) malloc (sizeof(float) * DIMENSION);
463 *test2 = (float *) malloc (sizeof(float) * DIMENSION);
464
465 *simplex = (float **) malloc (sizeof(float *) * (DIMENSION+1));
466
467 for (i = 0; i < DIMENSION+1; i++)
468 (*simplex)[i] = (float *) malloc (sizeof(float) * DIMENSION);
469
470 }
471
472
473 /*---------------------------------------------------------------------------*/
474
deallocate_arrays(float *** simplex,float ** centroid,float ** response,float ** step_size,float ** test1,float ** test2)475 void deallocate_arrays (float *** simplex, float ** centroid,
476 float ** response, float ** step_size,
477 float ** test1, float ** test2)
478 {
479 int i;
480
481 free (*centroid); *centroid = NULL;
482 free (*response); *response = NULL;
483 free (*step_size); *step_size = NULL;
484 free (*test1); *test1 = NULL;
485 free (*test2); *test2 = NULL;
486
487 for (i = 0; i < DIMENSION+1; i++)
488 {
489 free ((*simplex)[i]);
490 (*simplex)[i] = NULL;
491 }
492
493 free (*simplex); *simplex = NULL;
494
495 }
496
497
498 /*---------------------------------------------------------------------------*/
499
eval_vertices(float * response,int * worst,int * next,int * best)500 void eval_vertices (float * response, int * worst, int * next, int * best)
501 {
502 int i;
503
504 /* initialize values */
505 *worst = 0;
506 *best = 0;
507
508 /* find the best and worst */
509 for (i = 1; i < DIMENSION+1; i++)
510 {
511 if (response[i] > response[*worst])
512 *worst = i;
513 if (response[i] < response[*best])
514 *best = i;
515 }
516
517 /* find the next worst index */
518 if (*worst == 0)
519 *next = 1;
520 else
521 *next = 0;
522
523 for (i = 0; i < DIMENSION+1; i++)
524 if ((i != *worst) && (response[i] > response[*next]))
525 *next = i;
526 }
527
528
529 /*---------------------------------------------------------------------------*/
530
restart(float ** simplex,float * response,float * step_size)531 void restart (float ** simplex, float * response,
532 float * step_size)
533 {
534 const float STEP_FACTOR = 0.9;
535 int i, j;
536 int worst, next, best;
537 float minval, maxval;
538
539
540 /* find the current best vertex */
541 eval_vertices (response, &worst, &next, &best);
542
543 /* set the first vertex to the current best */
544 for (i = 0; i < DIMENSION; i++)
545 simplex[0][i] = simplex[best][i];
546
547 /* decrease step size */
548 for (i = 0; i < DIMENSION; i++)
549 step_size[i] *= STEP_FACTOR;
550
551 /* set up remaining vertices of simplex using new step size */
552 for (i = 1; i < DIMENSION+1; i++)
553 for (j = 0; j < DIMENSION; j++)
554 {
555 minval = simplex[0][j] - step_size[j];
556 maxval = simplex[0][j] + step_size[j];
557 simplex[i][j] = rand_uniform (minval, maxval);
558 }
559
560 /* initialize response for each vector */
561 for (i = 0; i < DIMENSION+1; i++)
562 response[i] = calc_sse (simplex[i]);
563 }
564
565
566 /*---------------------------------------------------------------------------*/
567 /*
568 Calculate the centroid of the simplex, ignoring the worst vertex.
569 */
570
calc_centroid(float ** simplex,int worst,float * centroid)571 void calc_centroid (float ** simplex, int worst, float * centroid)
572 {
573 int i, j;
574
575 for (i = 0; i < DIMENSION; i++)
576 {
577 centroid[i] = 0.0;
578
579 /* add each vertex, except the worst */
580 for (j = 0; j < DIMENSION+1; j++)
581 if (j != worst)
582 centroid[i] += simplex[j][i];
583 }
584
585 /* divide by the number of vertices */
586 for (i = 0; i < DIMENSION; i++)
587 centroid[i] /= DIMENSION;
588 }
589
590
591 /*---------------------------------------------------------------------------*/
592 /*
593 Calculate the reflection of the worst vertex about the centroid.
594 */
595
calc_reflection(float ** simplex,float * centroid,int worst,float coef,float * vertex)596 void calc_reflection (float ** simplex, float * centroid,
597 int worst, float coef, float * vertex)
598 {
599 int i;
600
601 for (i = 0; i < DIMENSION; i++)
602 vertex[i] = centroid[i] + coef*(centroid[i] - simplex[worst][i]);
603 }
604
605
606 /*---------------------------------------------------------------------------*/
607 /*
608 Replace a vertex of the simplex.
609 */
610
replace(float ** simplex,float * response,int index,float * vertex,float resp)611 void replace (float ** simplex, float * response,
612 int index, float * vertex, float resp)
613 {
614 int i;
615
616 for (i = 0; i < DIMENSION; i++)
617 simplex[index][i] = vertex[i];
618
619 response[index] = resp;
620 }
621
622
623 /*---------------------------------------------------------------------------*/
624 /*
625 Calculate goodness of fit.
626 */
627
calc_good_fit(float * response)628 float calc_good_fit (float * response)
629 {
630 int i;
631
632 float avg, sd, tmp;
633
634 /* average the response values */
635 avg = 0.0;
636 for (i = 0; i < DIMENSION+1; i++)
637 avg += response[i];
638 avg /= DIMENSION+1;
639
640 /* compute standard deviation of response */
641 sd = 0.0;
642 for (i = 0; i < DIMENSION+1; i++)
643 {
644 tmp = response[i] - avg;
645 sd += tmp*tmp;
646 }
647 sd /= DIMENSION;
648
649 return (sqrt(sd));
650 }
651
652
653 /*---------------------------------------------------------------------------*/
654 /*
655 Perform initialization for the Simplex algorithm.
656 */
657
simplex_initialize(float * parameters,float ** simplex,float * response,float * step_size)658 void simplex_initialize (float * parameters, float ** simplex,
659 float * response, float * step_size)
660 {
661 int i, j;
662 int worst, next, best;
663 float resp;
664 float minval, maxval;
665
666
667 for (j = 0; j < DIMENSION; j++)
668 {
669 simplex[0][j] = parameters[j];
670 step_size[j] = 0.5 * parameters[j];
671 }
672
673 for (i = 1; i < DIMENSION+1; i++)
674 for (j = 0; j < DIMENSION; j++)
675 {
676 minval = simplex[0][j] - step_size[j];
677 maxval = simplex[0][j] + step_size[j];
678 simplex[i][j] = rand_uniform (minval, maxval);
679 }
680
681 for (i = 0; i < DIMENSION+1; i++)
682 response[i] = calc_sse (simplex[i]);
683
684 for (i = 1; i < 500; i++)
685 {
686 for (j = 0; j < DIMENSION; j++)
687 {
688 minval = simplex[0][j] - step_size[j];
689 maxval = simplex[0][j] + step_size[j];
690 parameters[j] = rand_uniform (minval, maxval);
691 }
692
693 resp = calc_sse (parameters);
694 eval_vertices (response, &worst, &next, &best);
695 if (resp < response[worst])
696 replace (simplex, response, worst, parameters, resp);
697 }
698 }
699
700
701 /*---------------------------------------------------------------------------*/
702 /*
703 Use Simplex algorithm to estimate the voxel intensity distribution.
704
705 The Simplex algorithm is adapted from: A Jump Start Course in C++ Programming
706 */
707
simplex_optimization(float * parameters,float * sse)708 void simplex_optimization (float * parameters, float * sse)
709 {
710 const int MAX_ITERATIONS = 100;
711 const int MAX_RESTARTS = 25;
712 const float EXPANSION_COEF = 2.0;
713 const float REFLECTION_COEF = 1.0;
714 const float CONTRACTION_COEF = 0.5;
715 const float TOLERANCE = 1.0e-10;
716
717 float ** simplex = NULL;
718 float * centroid = NULL;
719 float * response = NULL;
720 float * step_size = NULL;
721 float * test1 = NULL;
722 float * test2 = NULL;
723 float resp1, resp2;
724 int i, worst, best, next;
725 int num_iter, num_restarts;
726 int done;
727 float fit;
728
729
730 allocate_arrays (&simplex, ¢roid, &response, &step_size, &test1, &test2);
731
732 simplex_initialize (parameters, simplex, response, step_size);
733
734 /* start loop to do simplex optimization */
735 num_iter = 0;
736 num_restarts = 0;
737 done = 0;
738
739 while (!done)
740 {
741 /* find the worst vertex and compute centroid of remaining simplex,
742 discarding the worst vertex */
743 eval_vertices (response, &worst, &next, &best);
744 calc_centroid (simplex, worst, centroid);
745
746 /* reflect the worst point through the centroid */
747 calc_reflection (simplex, centroid, worst,
748 REFLECTION_COEF, test1);
749 resp1 = calc_sse (test1);
750
751 /* test the reflection against the best vertex and expand it if the
752 reflection is better. if not, keep the reflection */
753 if (resp1 < response[best])
754 {
755 /* try expanding */
756 calc_reflection (simplex, centroid, worst, EXPANSION_COEF, test2);
757 resp2 = calc_sse (test2);
758 if (resp2 <= resp1) /* keep expansion */
759 replace (simplex, response, worst, test2, resp2);
760 else /* keep reflection */
761 replace (simplex, response, worst, test1, resp1);
762 }
763 else if (resp1 < response[next])
764 {
765 /* new response is between the best and next worst
766 so keep reflection */
767 replace (simplex, response, worst, test1, resp1);
768 }
769 else
770 {
771 /* try contraction */
772 if (resp1 >= response[worst])
773 calc_reflection (simplex, centroid, worst,
774 -CONTRACTION_COEF, test2);
775 else
776 calc_reflection (simplex, centroid, worst,
777 CONTRACTION_COEF, test2);
778 resp2 = calc_sse (test2);
779
780 /* test the contracted response against the worst response */
781 if (resp2 > response[worst])
782 {
783 /* new contracted response is worse, so decrease step size
784 and restart */
785 num_iter = 0;
786 num_restarts += 1;
787 restart (simplex, response, step_size);
788 }
789 else /* keep contraction */
790 replace (simplex, response, worst, test2, resp2);
791 }
792
793 /* test to determine when to stop.
794 first, check the number of iterations */
795 num_iter += 1; /* increment iteration counter */
796 if (num_iter >= MAX_ITERATIONS)
797 {
798 /* restart with smaller steps */
799 num_iter = 0;
800 num_restarts += 1;
801 restart (simplex, response, step_size);
802 }
803
804 /* limit the number of restarts */
805 if (num_restarts == MAX_RESTARTS) done = 1;
806
807 /* compare relative standard deviation of vertex responses
808 against a defined tolerance limit */
809 fit = calc_good_fit (response);
810 if (fit <= TOLERANCE) done = 1;
811
812 /* if done, copy the best solution to the output array */
813 if (done)
814 {
815 eval_vertices (response, &worst, &next, &best);
816 for (i = 0; i < DIMENSION; i++)
817 parameters[i] = simplex[best][i];
818 *sse = response[best];
819 }
820
821 } /* while (!done) */
822
823 number_restarts = num_restarts;
824 deallocate_arrays (&simplex, ¢roid, &response, &step_size,
825 &test1, &test2);
826
827 }
828
829
830 /*---------------------------------------------------------------------------*/
831 /*
832 Write the parameter estimates.
833 */
834
output_pdf_results(float * vertex,float sse)835 void output_pdf_results (float * vertex, float sse)
836 {
837 float b, bmean, bsigma, g, gmean, gsigma, w, wmean, wsigma;
838
839
840 /*----- Assign variables -----*/
841 b = vertex[0];
842 bmean = vertex[1];
843 bsigma = vertex[2];
844 g = vertex[3];
845 gmean = vertex[4];
846 gsigma = vertex[5];
847 w = vertex[6];
848 wmean = vertex[7];
849 wsigma = vertex[8];
850
851 printf ("hist_min = %d hist_max = %d \n", hist_min, hist_max);
852 printf ("gpeak = %d wpeak = %d \n", gpeak, wpeak);
853 printf ("rmse = %f \n", sqrt (sse / (hist_max - hist_min) ) );
854 printf ("\nProbability Density Function Estimates: \n");
855 printf ("Background Coef = %f \n", b);
856 printf ("Background Mean = %f \n", bmean);
857 printf ("Background Std Dev = %f \n", bsigma);
858 printf ("Gray Matter Coef = %f \n", g);
859 printf ("Gray Matter Mean = %f \n", gmean);
860 printf ("Gray Matter Std Dev = %f \n", gsigma);
861 printf ("White Matter Coef = %f \n", w);
862 printf ("White Matter Mean = %f \n", wmean);
863 printf ("White Matter Std Dev = %f \n", wsigma);
864 }
865
866
867 /*---------------------------------------------------------------------------*/
868 /*
869 Estimate the PDF of the voxel intensities.
870 */
871
estpdf(float * parameters)872 RwcBoolean estpdf (float * parameters)
873 {
874 float sse;
875 RwcBoolean ok;
876
877
878 /*----- Progress report -----*/
879 if (! quiet) printf ("\nEstimating PDF of voxel intensities \n");
880
881
882 /*----- Initialization for PDF estimation -----*/
883 ok = estpdf_initialize ();
884 if (! ok) return (FALSE);
885
886 generate_initial_guess (parameters);
887
888
889 /*----- Get least squares estimate for PDF parameters -----*/
890 simplex_optimization (parameters, &sse);
891
892
893 /*----- Report PDF parameters -----*/
894 if (! quiet) output_pdf_results (parameters, sse);
895
896
897 /*----- Free memory -----*/
898 free (hist_data); hist_data = NULL;
899
900
901 return (TRUE);
902 }
903
904
905 /*---------------------------------------------------------------------------*/
906