1 /*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
4 Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5 Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6
7 This library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public License
9 (LGPL) as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
11 The terms of redistributing and/or modifying this software also
12 include exceptions to the LGPL that facilitate static linking.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with this library; if not, write to Free Software Foundation, Inc.,
21 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22 */
23
24 #include "ten.h"
25 #include "privateTen.h"
26
27 tenEMBimodalParm*
tenEMBimodalParmNew()28 tenEMBimodalParmNew() {
29 tenEMBimodalParm *biparm;
30
31 biparm = (tenEMBimodalParm*)calloc(1, sizeof(tenEMBimodalParm));
32 if (biparm) {
33 biparm->minProb = 0.0001;
34 biparm->minProb2 = 0.0001;
35 biparm->minDelta = 0.00001;
36 biparm->minFraction = 0.05; /* 5% */
37 biparm->minConfidence = 0.7;
38 biparm->maxIteration = 200;
39 biparm->verbose = AIR_FALSE;
40
41 biparm->histo = NULL;
42 biparm->pp1 = biparm->pp2 = NULL;
43 biparm->vmin = biparm->vmax = AIR_NAN;
44 biparm->N = 0;
45 }
46 return biparm;
47 }
48
49 tenEMBimodalParm*
tenEMBimodalParmNix(tenEMBimodalParm * biparm)50 tenEMBimodalParmNix(tenEMBimodalParm *biparm) {
51
52 if (biparm) {
53 biparm->histo = (double *)airFree(biparm->histo);
54 biparm->pp1 = (double *)airFree(biparm->pp1);
55 biparm->pp2 = (double *)airFree(biparm->pp2);
56 }
57 airFree(biparm);
58 return NULL;
59 }
60
61 int
_tenEMBimodalInit(tenEMBimodalParm * biparm,const Nrrd * _nhisto)62 _tenEMBimodalInit(tenEMBimodalParm *biparm, const Nrrd *_nhisto) {
63 static const char me[]="_tenEMBimodalInit";
64 int i, median;
65 Nrrd *nhisto;
66 double medianD, sum;
67 airArray *mop;
68
69 if (!( biparm->maxIteration > 5 )) {
70 biffAddf(TEN, "%s: biparm->maxIteration = %d too small", me,
71 biparm->maxIteration);
72 return 1;
73 }
74
75 mop = airMopNew();
76 nhisto = nrrdNew();
77 airMopAdd(mop, nhisto, (airMopper)nrrdNuke, airMopOnError);
78 airMopAdd(mop, nhisto, (airMopper)nrrdNix, airMopOnOkay);
79 if (nrrdConvert(nhisto, _nhisto, nrrdTypeDouble)) {
80 biffMovef(TEN, NRRD, "%s: trouble converting histogram to double", me);
81 airMopError(mop); return 1;
82 }
83 biparm->N = nhisto->axis[0].size;
84 biparm->histo = (double*)(nhisto->data);
85 biparm->vmin = (AIR_EXISTS(nhisto->axis[0].min)
86 ? nhisto->axis[0].min
87 : -0.5);
88 biparm->vmax = (AIR_EXISTS(nhisto->axis[0].max)
89 ? nhisto->axis[0].max
90 : biparm->N - 0.5);
91
92 (nrrdMeasureLine[nrrdMeasureHistoMedian])
93 (&medianD, nrrdTypeDouble,
94 biparm->histo, nrrdTypeDouble, biparm->N,
95 AIR_NAN, AIR_NAN);
96 (nrrdMeasureLine[nrrdMeasureSum])
97 (&sum, nrrdTypeDouble,
98 biparm->histo, nrrdTypeDouble, biparm->N,
99 AIR_NAN, AIR_NAN);
100 for (i=0; i<biparm->N; i++) {
101 biparm->histo[i] /= sum;
102 }
103 if (!AIR_EXISTS(medianD)) {
104 biffMovef(TEN, NRRD,
105 "%s: got empty histogram? (median calculation failed)", me);
106 airMopError(mop); return 1;
107 }
108 median = (int)medianD;
109
110 biparm->pp1 = (double*)calloc(biparm->N, sizeof(double));
111 biparm->pp2 = (double*)calloc(biparm->N, sizeof(double));
112 if (!( biparm->pp1 && biparm->pp2 )) {
113 biffAddf(TEN, "%s: couldn't allocate posterior prob. buffers", me);
114 airMopError(mop); return 1;
115 }
116
117 /* get mean and stdv of bins below median */
118 (nrrdMeasureLine[nrrdMeasureHistoMean])
119 (&(biparm->mean1), nrrdTypeDouble,
120 biparm->histo, nrrdTypeDouble, median,
121 AIR_NAN, AIR_NAN);
122 (nrrdMeasureLine[nrrdMeasureHistoSD])
123 (&(biparm->stdv1), nrrdTypeDouble,
124 biparm->histo, nrrdTypeDouble, median,
125 AIR_NAN, AIR_NAN);
126
127 /* get mean (shift upwards by median) and stdv of bins above median */
128 (nrrdMeasureLine[nrrdMeasureHistoMean])
129 (&(biparm->mean2), nrrdTypeDouble,
130 biparm->histo + median, nrrdTypeDouble, biparm->N - median,
131 AIR_NAN, AIR_NAN);
132 (nrrdMeasureLine[nrrdMeasureHistoSD])
133 (&(biparm->stdv2), nrrdTypeDouble,
134 biparm->histo + median, nrrdTypeDouble, biparm->N - median,
135 AIR_NAN, AIR_NAN);
136
137 biparm->mean2 += median;
138 biparm->fraction1 = 0.5;
139
140 if (biparm->verbose) {
141 fprintf(stderr, "%s: median = %d\n", me, median);
142 fprintf(stderr, "%s: m1, s1 = %g, %g; m2, s2 = %g, %g\n", me,
143 biparm->mean1, biparm->stdv1,
144 biparm->mean2, biparm->stdv2);
145 }
146
147 airMopOkay(mop);
148 return 0;
149 }
150
151 void
_tenEMBimodalBoost(double * pp1P,double * pp2P,double b)152 _tenEMBimodalBoost(double *pp1P, double *pp2P, double b) {
153 double p1, p2, tmp;
154 int sw=AIR_FALSE;
155
156 if (*pp1P < *pp2P) {
157 ELL_SWAP2(*pp1P, *pp2P, tmp);
158 sw = AIR_TRUE;
159 }
160 p1 = 1 - pow(1 - *pp1P, b);
161 p2 = 1 - p1;
162 if (sw) {
163 *pp1P = p2;
164 *pp2P = p1;
165 } else {
166 *pp1P = p1;
167 *pp2P = p2;
168 }
169 }
170
171 /*
172 ** what is posterior probability that measured value x comes from
173 ** material 1 and 2, stored in pp1 and pp2
174 */
175 void
_tenEMBimodalPP(tenEMBimodalParm * biparm)176 _tenEMBimodalPP(tenEMBimodalParm *biparm) {
177 int i;
178 double g1, g2, pp1, pp2, f1, min;
179
180 min = (1 == biparm->stage
181 ? biparm->minProb
182 : biparm->minProb2);
183 f1 = biparm->fraction1;
184 for (i=0; i<biparm->N; i++) {
185 g1 = airGaussian(i, biparm->mean1, biparm->stdv1);
186 g2 = airGaussian(i, biparm->mean2, biparm->stdv2);
187 if (g1 <= min && g2 <= min) {
188 pp1 = pp2 = 0;
189 } else {
190 pp1 = f1*g1 / (f1*g1 + (1-f1)*g2);
191 pp2 = 1 - pp1;
192 }
193 biparm->pp1[i] = pp1;
194 biparm->pp2[i] = pp2;
195 }
196
197 if (biparm->verbose > 1) {
198 Nrrd *ntmp = nrrdNew();
199 nrrdWrap_va(ntmp, biparm->pp1, nrrdTypeDouble, 1,
200 AIR_CAST(size_t, biparm->N));
201 nrrdSave("pp1.nrrd", ntmp, NULL);
202 nrrdWrap_va(ntmp, biparm->pp2, nrrdTypeDouble, 1,
203 AIR_CAST(size_t, biparm->N));
204 nrrdSave("pp2.nrrd", ntmp, NULL);
205 nrrdNix(ntmp);
206 }
207
208 return;
209 }
210
211 double
_tenEMBimodalNewFraction1(tenEMBimodalParm * biparm)212 _tenEMBimodalNewFraction1(tenEMBimodalParm *biparm) {
213 int i;
214 double pp1, pp2, h, sum1, sum2;
215
216 sum1 = sum2 = 0.0;
217 for (i=0; i<biparm->N; i++) {
218 pp1 = biparm->pp1[i];
219 pp2 = biparm->pp2[i];
220 h = biparm->histo[i];
221 sum1 += pp1*h;
222 sum2 += pp2*h;
223 }
224 return sum1/(sum1 + sum2);
225 }
226
227 void
_tenEMBimodalNewMean(double * m1P,double * m2P,tenEMBimodalParm * biparm)228 _tenEMBimodalNewMean(double *m1P, double *m2P,
229 tenEMBimodalParm *biparm) {
230 int i;
231 double pp1, pp2, h, sum1, isum1, sum2, isum2;
232
233 sum1 = isum1 = sum2 = isum2 = 0.0;
234 for (i=0; i<biparm->N; i++) {
235 pp1 = biparm->pp1[i];
236 pp2 = biparm->pp2[i];
237 h = biparm->histo[i];
238 isum1 += i*pp1*h;
239 isum2 += i*pp2*h;
240 sum1 += pp1*h;
241 sum2 += pp2*h;
242 }
243 *m1P = isum1/sum1;
244 *m2P = isum2/sum2;
245 }
246
247 void
_tenEMBimodalNewSigma(double * s1P,double * s2P,double m1,double m2,tenEMBimodalParm * biparm)248 _tenEMBimodalNewSigma(double *s1P, double *s2P,
249 double m1, double m2,
250 tenEMBimodalParm *biparm) {
251 int i;
252 double pp1, pp2, h, sum1, isum1, sum2, isum2;
253
254 sum1 = isum1 = sum2 = isum2 = 0.0;
255 for (i=0; i<biparm->N; i++) {
256 pp1 = biparm->pp1[i];
257 pp2 = biparm->pp2[i];
258 h = biparm->histo[i];
259 isum1 += (i-m1)*(i-m1)*pp1*h;
260 isum2 += (i-m2)*(i-m2)*pp2*h;
261 sum1 += pp1*h;
262 sum2 += pp2*h;
263 }
264 *s1P = sqrt(isum1/sum1);
265 *s2P = sqrt(isum2/sum2);
266 }
267
268 void
_tenEMBimodalSaveImage(tenEMBimodalParm * biparm)269 _tenEMBimodalSaveImage(tenEMBimodalParm *biparm) {
270 char name[AIR_STRLEN_MED];
271 Nrrd *nh, *nm, *nhi, *nmi, *ni;
272 NrrdRange *range;
273 const Nrrd *nhmhi[3];
274 double *m, max;
275 int i;
276
277 nh = nrrdNew();
278 nm = nrrdNew();
279 nhi = nrrdNew();
280 nmi = nrrdNew();
281 ni = nrrdNew();
282 nrrdWrap_va(nh, biparm->histo, nrrdTypeDouble, 1,
283 AIR_CAST(size_t, biparm->N));
284 range = nrrdRangeNewSet(nh, nrrdBlind8BitRangeFalse);
285 max = range->max*1.1;
286 nrrdRangeNix(range);
287 nrrdCopy(nm, nh);
288 m = (double*)(nm->data);
289 for (i=0; i<biparm->N; i++) {
290 m[i] = biparm->fraction1*airGaussian(i, biparm->mean1, biparm->stdv1);
291 m[i] += (1-biparm->fraction1)*airGaussian(i, biparm->mean2, biparm->stdv2);
292 }
293 nrrdHistoDraw(nmi, nm, 400, AIR_FALSE, max);
294 nrrdHistoDraw(nhi, nh, 400, AIR_FALSE, max);
295 ELL_3V_SET(nhmhi, nhi, nmi, nhi);
296 nrrdJoin(ni, nhmhi, 3, 0, AIR_TRUE);
297 sprintf(name, "%04d-%d.png", biparm->iteration, biparm->stage);
298 nrrdSave(name, ni, NULL);
299 nh = nrrdNix(nh);
300 nm = nrrdNuke(nm);
301 nhi = nrrdNuke(nhi);
302 nmi = nrrdNuke(nmi);
303 ni = nrrdNuke(ni);
304 return;
305 }
306
307
308 int
_tenEMBimodalIterate(tenEMBimodalParm * biparm)309 _tenEMBimodalIterate(tenEMBimodalParm *biparm) {
310 static const char me[]="_tenEMBimodalIterate";
311 double om1, os1, om2, os2, of1, m1, s1, m2, s2, f1;
312
313 /* copy old values */
314 om1 = biparm->mean1;
315 os1 = biparm->stdv1;
316 of1 = biparm->fraction1;
317 om2 = biparm->mean2;
318 os2 = biparm->stdv2;
319
320 /* find new values, and calculate delta */
321 _tenEMBimodalPP(biparm);
322 f1 = _tenEMBimodalNewFraction1(biparm);
323 /* if (1 == biparm->stage) { */
324 _tenEMBimodalNewMean(&m1, &m2, biparm);
325 /* } */
326 _tenEMBimodalNewSigma(&s1, &s2, m1, m2, biparm);
327
328 biparm->delta = ((fabs(m1 - om1) + fabs(m2 - om2)
329 + fabs(s1 - os1) + fabs(s2 - os2))/biparm->N
330 + fabs(f1 - of1));
331
332 /* set new values */
333 biparm->mean1 = m1;
334 biparm->stdv1 = s1;
335 biparm->fraction1 = f1;
336 biparm->mean2 = m2;
337 biparm->stdv2 = s2;
338
339 if (biparm->verbose) {
340 fprintf(stderr, "%s(%d:%d):\n", me, biparm->stage, biparm->iteration);
341 fprintf(stderr, " m1, s1 = %g, %g\n", m1, s1);
342 fprintf(stderr, " m2, s2 = %g, %g\n", m2, s2);
343 fprintf(stderr, " f1 = %g ; delta = %g\n", f1, biparm->delta);
344 }
345 if (biparm->verbose > 1) {
346 _tenEMBimodalSaveImage(biparm);
347 }
348 return 0;
349 }
350
351 int
_tenEMBimodalConfThresh(tenEMBimodalParm * biparm)352 _tenEMBimodalConfThresh(tenEMBimodalParm *biparm) {
353 static const char me[]="_tenEMBimodalConfThresh";
354 double m1, s1, m2, s2, f1, f2, A, B, C, D, t1, t2;
355
356 biparm->confidence = ((biparm->mean2 - biparm->mean1)
357 / (biparm->stdv1 + biparm->stdv2));
358 m1 = biparm->mean1;
359 s1 = biparm->stdv1;
360 f1 = biparm->fraction1;
361 m2 = biparm->mean2;
362 s2 = biparm->stdv2;
363 f2 = 1 - f1;
364 A = s1*s1 - s2*s2;
365 B = 2*(m1*s2*s2 - m2*s1*s1);
366 C = s1*s1*m2*m2 - s2*s2*m1*m1 + 4*s1*s1*s2*s2*log(s2*f1/(s1*f2));
367 D = B*B - 4*A*C;
368 if (D < 0) {
369 biffAddf(TEN, "%s: threshold descriminant went negative (%g)", me, D);
370 return 1;
371 }
372 t1 = (-B + sqrt(D))/(2*A);
373 if (AIR_IN_OP(m1, t1, m2)) {
374 biparm->threshold = t1;
375 } else {
376 t2 = (-B - sqrt(D))/(2*A);
377 if (AIR_IN_OP(m1, t2, m2)) {
378 biparm->threshold = t2;
379 } else {
380 biffAddf(TEN,
381 "%s: neither computed threshold %g,%g inside open interval "
382 "between means (%g,%g)", me, t1, t2, m1, m2);
383 return 1;
384 }
385 }
386
387 if (biparm->verbose) {
388 fprintf(stderr, "%s: conf = %g, thresh = %g\n", me,
389 biparm->confidence, biparm->threshold);
390 }
391 return 0;
392 }
393
394 int
_tenEMBimodalCheck(tenEMBimodalParm * biparm)395 _tenEMBimodalCheck(tenEMBimodalParm *biparm) {
396 static const char me[]="_tenEMBimodalCheck";
397
398 if (!( biparm->confidence > biparm->minConfidence )) {
399 biffAddf(TEN, "%s: confidence %g went below threshold %g", me,
400 biparm->confidence, biparm->minConfidence);
401 return 1;
402 }
403 if (!( biparm->stdv1 > 0 && biparm->stdv2 > 0 )) {
404 biffAddf(TEN, "%s: stdv of material 1 (%g) or 2 (%g) went negative", me,
405 biparm->stdv1, biparm->stdv2);
406 return 1;
407 }
408 if (!( biparm->mean1 > 0 && biparm->mean1 < biparm->N-1
409 && biparm->mean2 > 0 && biparm->mean2 < biparm->N-1 )) {
410 biffAddf(TEN, "%s: mean of material 1 (%g) or 2 (%g) went outside "
411 "given histogram range [0 .. %d]", me,
412 biparm->mean1, biparm->mean2, biparm->N-1);
413 return 1;
414 }
415 if (biparm->fraction1 < biparm->minFraction) {
416 biffAddf(TEN, "%s: material 1 fraction (%g) fell below threshold %g", me,
417 biparm->fraction1, biparm->minFraction);
418 return 1;
419 }
420 if (1 - biparm->fraction1 < biparm->minFraction) {
421 biffAddf(TEN, "%s: material 2 fraction (%g) fell below threshold %g", me,
422 1 - biparm->fraction1, biparm->minFraction);
423 return 1;
424 }
425 return 0;
426 }
427
428 int
tenEMBimodal(tenEMBimodalParm * biparm,const Nrrd * _nhisto)429 tenEMBimodal(tenEMBimodalParm *biparm, const Nrrd *_nhisto) {
430 static const char me[]="tenEMBimodal";
431 int done, _iter;
432
433 if (!(biparm && _nhisto)) {
434 biffAddf(TEN, "%s: got NULL pointer", me);
435 return 1;
436 }
437 if (!( 1 == _nhisto->dim )) {
438 biffAddf(TEN, "%s: histogram must be 1-D, not %d-D", me, _nhisto->dim);
439 return 1;
440 }
441
442 if (_tenEMBimodalInit(biparm, _nhisto)) {
443 biffAddf(TEN, "%s: trouble initializing parameters", me);
444 return 1;
445 }
446
447 done = AIR_FALSE;
448 biparm->iteration = 0;
449 for (biparm->stage = 1;
450 biparm->stage <= (biparm->twoStage ? 2 : 1);
451 biparm->stage++) {
452 for (_iter=0;
453 biparm->iteration <= biparm->maxIteration;
454 biparm->iteration++, _iter++) {
455 if (_tenEMBimodalIterate(biparm) /* sets delta */
456 || _tenEMBimodalConfThresh(biparm)
457 || _tenEMBimodalCheck(biparm)) {
458 biffAddf(TEN, "%s: problem with fitting (iter=%d)", me,
459 biparm->iteration);
460 return 1;
461 }
462 if (biparm->delta < biparm->minDelta
463 && (!biparm->twoStage || 1 == biparm->stage || _iter > 10) ) {
464 done = AIR_TRUE;
465 break;
466 }
467 }
468 }
469 if (!done) {
470 biffAddf(TEN, "%s: didn't converge after %d iterations", me,
471 biparm->maxIteration);
472 return 1;
473 }
474
475 return 0;
476 }
477