Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	src/ten/bimod.c	Lines:	0	268	0.0 %
Date:	2017-05-26	Branches:	0	102	0.0 %


/*
  Teem: Tools to process and visualize scientific data and images             .
  Copyright (C) 2013, 2012, 2011, 2010, 2009  University of Chicago
  Copyright (C) 2008, 2007, 2006, 2005  Gordon Kindlmann
  Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998  University of Utah

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public License
  (LGPL) as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  The terms of redistributing and/or modifying this software also
  include exceptions to the LGPL that facilitate static linking.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public License
  along with this library; if not, write to Free Software Foundation, Inc.,
  51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
*/

#include "ten.h"
#include "privateTen.h"

tenEMBimodalParm*
tenEMBimodalParmNew() {
  tenEMBimodalParm *biparm;

  biparm = (tenEMBimodalParm*)calloc(1, sizeof(tenEMBimodalParm));
  if (biparm) {
    biparm->minProb = 0.0001;
    biparm->minProb2 = 0.0001;
    biparm->minDelta = 0.00001;
    biparm->minFraction = 0.05;  /* 5% */
    biparm->minConfidence = 0.7;
    biparm->maxIteration = 200;
    biparm->verbose = AIR_FALSE;

    biparm->histo = NULL;
    biparm->pp1 = biparm->pp2 = NULL;
    biparm->vmin = biparm->vmax = AIR_NAN;
    biparm->N = 0;
  }
  return biparm;
}

tenEMBimodalParm*
tenEMBimodalParmNix(tenEMBimodalParm *biparm) {

  if (biparm) {
    biparm->histo = (double *)airFree(biparm->histo);
    biparm->pp1 = (double *)airFree(biparm->pp1);
    biparm->pp2 = (double *)airFree(biparm->pp2);
  }
  airFree(biparm);
  return NULL;
}

int
_tenEMBimodalInit(tenEMBimodalParm *biparm, const Nrrd *_nhisto) {
  static const char me[]="_tenEMBimodalInit";
  int i, median;
  Nrrd *nhisto;
  double medianD, sum;
  airArray *mop;

  if (!( biparm->maxIteration > 5 )) {
    biffAddf(TEN, "%s: biparm->maxIteration = %d too small", me,
             biparm->maxIteration);
    return 1;
  }

  mop = airMopNew();
  nhisto = nrrdNew();
  airMopAdd(mop, nhisto, (airMopper)nrrdNuke, airMopOnError);
  airMopAdd(mop, nhisto, (airMopper)nrrdNix, airMopOnOkay);
  if (nrrdConvert(nhisto, _nhisto, nrrdTypeDouble)) {
    biffMovef(TEN, NRRD, "%s: trouble converting histogram to double", me);
    airMopError(mop); return 1;
  }
  biparm->N = nhisto->axis[0].size;
  biparm->histo = (double*)(nhisto->data);
  biparm->vmin = (AIR_EXISTS(nhisto->axis[0].min)
                  ? nhisto->axis[0].min
                  : -0.5);
  biparm->vmax = (AIR_EXISTS(nhisto->axis[0].max)
                  ? nhisto->axis[0].max
                  : biparm->N - 0.5);

  (nrrdMeasureLine[nrrdMeasureHistoMedian])
    (&medianD, nrrdTypeDouble,
     biparm->histo, nrrdTypeDouble, biparm->N,
     AIR_NAN, AIR_NAN);
  (nrrdMeasureLine[nrrdMeasureSum])
    (&sum, nrrdTypeDouble,
     biparm->histo, nrrdTypeDouble, biparm->N,
     AIR_NAN, AIR_NAN);
  for (i=0; i<biparm->N; i++) {
    biparm->histo[i] /= sum;
  }
  if (!AIR_EXISTS(medianD)) {
    biffMovef(TEN, NRRD,
              "%s: got empty histogram? (median calculation failed)", me);
    airMopError(mop); return 1;
  }
  median = (int)medianD;

  biparm->pp1 = (double*)calloc(biparm->N, sizeof(double));
  biparm->pp2 = (double*)calloc(biparm->N, sizeof(double));
  if (!( biparm->pp1 && biparm->pp2 )) {
    biffAddf(TEN, "%s: couldn't allocate posterior prob. buffers", me);
    airMopError(mop); return 1;
  }

  /* get mean and stdv of bins below median */
  (nrrdMeasureLine[nrrdMeasureHistoMean])
    (&(biparm->mean1), nrrdTypeDouble,
     biparm->histo, nrrdTypeDouble, median,
     AIR_NAN, AIR_NAN);
  (nrrdMeasureLine[nrrdMeasureHistoSD])
    (&(biparm->stdv1), nrrdTypeDouble,
     biparm->histo, nrrdTypeDouble, median,
     AIR_NAN, AIR_NAN);

  /* get mean (shift upwards by median) and stdv of bins above median */
  (nrrdMeasureLine[nrrdMeasureHistoMean])
    (&(biparm->mean2), nrrdTypeDouble,
     biparm->histo + median, nrrdTypeDouble, biparm->N - median,
     AIR_NAN, AIR_NAN);
  (nrrdMeasureLine[nrrdMeasureHistoSD])
    (&(biparm->stdv2), nrrdTypeDouble,
     biparm->histo + median, nrrdTypeDouble, biparm->N - median,
     AIR_NAN, AIR_NAN);

  biparm->mean2 += median;
  biparm->fraction1 = 0.5;

  if (biparm->verbose) {
    fprintf(stderr, "%s: median = %d\n", me, median);
    fprintf(stderr, "%s: m1, s1 = %g, %g; m2, s2 = %g, %g\n", me,
            biparm->mean1, biparm->stdv1,
            biparm->mean2, biparm->stdv2);
  }

  airMopOkay(mop);
  return 0;
}

void
_tenEMBimodalBoost(double *pp1P, double *pp2P, double b) {
  double p1, p2, tmp;
  int sw=AIR_FALSE;

  if (*pp1P < *pp2P) {
    ELL_SWAP2(*pp1P, *pp2P, tmp);
    sw = AIR_TRUE;
  }
  p1 = 1 - pow(1 - *pp1P, b);
  p2 = 1 - p1;
  if (sw) {
    *pp1P = p2;
    *pp2P = p1;
  } else {
    *pp1P = p1;
    *pp2P = p2;
  }
}

/*
** what is posterior probability that measured value x comes from
** material 1 and 2, stored in pp1 and pp2
*/
void
_tenEMBimodalPP(tenEMBimodalParm *biparm) {
  int i;
  double g1, g2, pp1, pp2, f1, min;

  min = (1 == biparm->stage
         ? biparm->minProb
         : biparm->minProb2);
  f1 = biparm->fraction1;
  for (i=0; i<biparm->N; i++) {
    g1 = airGaussian(i, biparm->mean1, biparm->stdv1);
    g2 = airGaussian(i, biparm->mean2, biparm->stdv2);
    if (g1 <= min && g2 <= min) {
      pp1 = pp2 = 0;
    } else {
      pp1 = f1*g1 / (f1*g1 + (1-f1)*g2);
      pp2 = 1 - pp1;
    }
    biparm->pp1[i] = pp1;
    biparm->pp2[i] = pp2;
  }

  if (biparm->verbose > 1) {
    Nrrd *ntmp = nrrdNew();
    nrrdWrap_va(ntmp, biparm->pp1, nrrdTypeDouble, 1,
                AIR_CAST(size_t, biparm->N));
    nrrdSave("pp1.nrrd", ntmp, NULL);
    nrrdWrap_va(ntmp, biparm->pp2, nrrdTypeDouble, 1,
                AIR_CAST(size_t, biparm->N));
    nrrdSave("pp2.nrrd", ntmp, NULL);
    nrrdNix(ntmp);
  }

  return;
}

double
_tenEMBimodalNewFraction1(tenEMBimodalParm *biparm) {
  int i;
  double pp1, pp2, h, sum1, sum2;

  sum1 = sum2 = 0.0;
  for (i=0; i<biparm->N; i++) {
    pp1 = biparm->pp1[i];
    pp2 = biparm->pp2[i];
    h = biparm->histo[i];
    sum1 += pp1*h;
    sum2 += pp2*h;
  }
  return sum1/(sum1 + sum2);
}

void
_tenEMBimodalNewMean(double *m1P, double *m2P,
                     tenEMBimodalParm *biparm) {
  int i;
  double pp1, pp2, h, sum1, isum1, sum2, isum2;

  sum1 = isum1 = sum2 = isum2 = 0.0;
  for (i=0; i<biparm->N; i++) {
    pp1 = biparm->pp1[i];
    pp2 = biparm->pp2[i];
    h = biparm->histo[i];
    isum1 += i*pp1*h;
    isum2 += i*pp2*h;
    sum1 += pp1*h;
    sum2 += pp2*h;
  }
  *m1P = isum1/sum1;
  *m2P = isum2/sum2;
}

void
_tenEMBimodalNewSigma(double *s1P, double *s2P,
                      double m1, double m2,
                      tenEMBimodalParm *biparm) {
  int i;
  double pp1, pp2, h, sum1, isum1, sum2, isum2;

  sum1 = isum1 = sum2 = isum2 = 0.0;
  for (i=0; i<biparm->N; i++) {
    pp1 = biparm->pp1[i];
    pp2 = biparm->pp2[i];
    h = biparm->histo[i];
    isum1 += (i-m1)*(i-m1)*pp1*h;
    isum2 += (i-m2)*(i-m2)*pp2*h;
    sum1 += pp1*h;
    sum2 += pp2*h;
  }
  *s1P = sqrt(isum1/sum1);
  *s2P = sqrt(isum2/sum2);
}

void
_tenEMBimodalSaveImage(tenEMBimodalParm *biparm) {
  char name[AIR_STRLEN_MED];
  Nrrd *nh, *nm, *nhi, *nmi, *ni;
  NrrdRange *range;
  const Nrrd *nhmhi[3];
  double *m, max;
  int i;

  nh = nrrdNew();
  nm = nrrdNew();
  nhi = nrrdNew();
  nmi = nrrdNew();
  ni = nrrdNew();
  nrrdWrap_va(nh, biparm->histo, nrrdTypeDouble, 1,
              AIR_CAST(size_t, biparm->N));
  range = nrrdRangeNewSet(nh, nrrdBlind8BitRangeFalse);
  max = range->max*1.1;
  nrrdRangeNix(range);
  nrrdCopy(nm, nh);
  m = (double*)(nm->data);
  for (i=0; i<biparm->N; i++) {
    m[i] = biparm->fraction1*airGaussian(i, biparm->mean1, biparm->stdv1);
    m[i] += (1-biparm->fraction1)*airGaussian(i, biparm->mean2, biparm->stdv2);
  }
  nrrdHistoDraw(nmi, nm, 400, AIR_FALSE, max);
  nrrdHistoDraw(nhi, nh, 400, AIR_FALSE, max);
  ELL_3V_SET(nhmhi, nhi, nmi, nhi);
  nrrdJoin(ni, nhmhi, 3, 0, AIR_TRUE);
  sprintf(name, "%04d-%d.png", biparm->iteration, biparm->stage);
  nrrdSave(name, ni, NULL);
  nh = nrrdNix(nh);
  nm = nrrdNuke(nm);
  nhi = nrrdNuke(nhi);
  nmi = nrrdNuke(nmi);
  ni = nrrdNuke(ni);
  return;
}


int
_tenEMBimodalIterate(tenEMBimodalParm *biparm) {
  static const char me[]="_tenEMBimodalIterate";
  double om1, os1, om2, os2, of1, m1, s1, m2, s2, f1;

  /* copy old values */
  om1 = biparm->mean1;
  os1 = biparm->stdv1;
  of1 = biparm->fraction1;
  om2 = biparm->mean2;
  os2 = biparm->stdv2;

  /* find new values, and calculate delta */
  _tenEMBimodalPP(biparm);
  f1 = _tenEMBimodalNewFraction1(biparm);
  /*   if (1 == biparm->stage) { */
    _tenEMBimodalNewMean(&m1, &m2, biparm);
    /*   } */
  _tenEMBimodalNewSigma(&s1, &s2, m1, m2, biparm);

  biparm->delta = ((fabs(m1 - om1) + fabs(m2 - om2)
                    + fabs(s1 - os1) + fabs(s2 - os2))/biparm->N
                   + fabs(f1 - of1));

  /* set new values */
  biparm->mean1 = m1;
  biparm->stdv1 = s1;
  biparm->fraction1 = f1;
  biparm->mean2 = m2;
  biparm->stdv2 = s2;

  if (biparm->verbose) {
    fprintf(stderr, "%s(%d:%d):\n", me, biparm->stage, biparm->iteration);
    fprintf(stderr, "  m1, s1 = %g, %g\n", m1, s1);
    fprintf(stderr, "  m2, s2 = %g, %g\n", m2, s2);
    fprintf(stderr, "  f1 = %g ; delta = %g\n", f1, biparm->delta);
  }
  if (biparm->verbose > 1) {
    _tenEMBimodalSaveImage(biparm);
  }
  return 0;
}

int
_tenEMBimodalConfThresh(tenEMBimodalParm *biparm) {
  static const char me[]="_tenEMBimodalConfThresh";
  double m1, s1, m2, s2, f1, f2, A, B, C, D, t1, t2;

  biparm->confidence = ((biparm->mean2 - biparm->mean1)
                        / (biparm->stdv1 + biparm->stdv2));
  m1 = biparm->mean1;
  s1 = biparm->stdv1;
  f1 = biparm->fraction1;
  m2 = biparm->mean2;
  s2 = biparm->stdv2;
  f2 = 1 - f1;
  A = s1*s1 - s2*s2;
  B = 2*(m1*s2*s2 - m2*s1*s1);
  C = s1*s1*m2*m2 - s2*s2*m1*m1 + 4*s1*s1*s2*s2*log(s2*f1/(s1*f2));
  D = B*B - 4*A*C;
  if (D < 0) {
    biffAddf(TEN, "%s: threshold descriminant went negative (%g)", me, D);
    return 1;
  }
  t1 = (-B + sqrt(D))/(2*A);
  if (AIR_IN_OP(m1, t1, m2)) {
    biparm->threshold = t1;
  } else {
    t2 = (-B - sqrt(D))/(2*A);
    if (AIR_IN_OP(m1, t2, m2)) {
      biparm->threshold = t2;
    } else {
      biffAddf(TEN,
               "%s: neither computed threshold %g,%g inside open interval "
               "between means (%g,%g)", me, t1, t2, m1, m2);
      return 1;
    }
  }

  if (biparm->verbose) {
    fprintf(stderr, "%s: conf = %g, thresh = %g\n", me,
            biparm->confidence, biparm->threshold);
  }
  return 0;
}

int
_tenEMBimodalCheck(tenEMBimodalParm *biparm) {
  static const char me[]="_tenEMBimodalCheck";

  if (!( biparm->confidence > biparm->minConfidence )) {
    biffAddf(TEN, "%s: confidence %g went below threshold %g", me,
             biparm->confidence, biparm->minConfidence);
    return 1;
  }
  if (!( biparm->stdv1 > 0 && biparm->stdv2 > 0 )) {
    biffAddf(TEN, "%s: stdv of material 1 (%g) or 2 (%g) went negative", me,
             biparm->stdv1, biparm->stdv2);
    return 1;
  }
  if (!( biparm->mean1 > 0 && biparm->mean1 < biparm->N-1
         && biparm->mean2 > 0 && biparm->mean2 < biparm->N-1 )) {
    biffAddf(TEN, "%s: mean of material 1 (%g) or 2 (%g) went outside "
             "given histogram range [0 .. %d]", me,
             biparm->mean1, biparm->mean2, biparm->N-1);
    return 1;
  }
  if (biparm->fraction1 < biparm->minFraction) {
    biffAddf(TEN, "%s: material 1 fraction (%g) fell below threshold %g", me,
             biparm->fraction1, biparm->minFraction);
    return 1;
  }
  if (1 - biparm->fraction1 < biparm->minFraction) {
    biffAddf(TEN, "%s: material 2 fraction (%g) fell below threshold %g", me,
             1 - biparm->fraction1, biparm->minFraction);
    return 1;
  }
  return 0;
}

int
tenEMBimodal(tenEMBimodalParm *biparm, const Nrrd *_nhisto) {
  static const char me[]="tenEMBimodal";
  int done, _iter;

  if (!(biparm && _nhisto)) {
    biffAddf(TEN, "%s: got NULL pointer", me);
    return 1;
  }
  if (!( 1 == _nhisto->dim )) {
    biffAddf(TEN, "%s: histogram must be 1-D, not %d-D", me, _nhisto->dim);
    return 1;
  }

  if (_tenEMBimodalInit(biparm, _nhisto)) {
    biffAddf(TEN, "%s: trouble initializing parameters", me);
    return 1;
  }

  done = AIR_FALSE;
  biparm->iteration = 0;
  for (biparm->stage = 1;
       biparm->stage <= (biparm->twoStage ? 2 : 1);
       biparm->stage++) {
    for (_iter=0;
         biparm->iteration <= biparm->maxIteration;
         biparm->iteration++, _iter++) {
      if (_tenEMBimodalIterate(biparm)    /* sets delta */
          || _tenEMBimodalConfThresh(biparm)
          || _tenEMBimodalCheck(biparm)) {
        biffAddf(TEN, "%s: problem with fitting (iter=%d)", me,
                 biparm->iteration);
        return 1;
      }
      if (biparm->delta < biparm->minDelta
          && (!biparm->twoStage || 1 == biparm->stage || _iter > 10) ) {
        done = AIR_TRUE;
        break;
      }
    }
  }
  if (!done) {
    biffAddf(TEN, "%s: didn't converge after %d iterations", me,
             biparm->maxIteration);
    return 1;
  }

  return 0;
}


Generated by: GCOVR (Version 3.3)

Line	Branch	Exec	Source
1			/*
2			Teem: Tools to process and visualize scientific data and images .
3			Copyright (C) 2013, 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*
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*
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
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
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
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 = f1g1 / (f1g1 + (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
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
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 += ipp1h;
239			isum2 += ipp2h;
240			sum1 += pp1*h;
241			sum2 += pp2*h;
242			}
243			*m1P = isum1/sum1;
244			*m2P = isum2/sum2;
245			}
246
247			void
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
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
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
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 = s1s1 - s2s2;
365			B = 2(m1s2s2 - m2s1*s1);
366			C = s1s1m2m2 - s2s2m1m1 + 4s1s1s2s2log(s2f1/(s1*f2));
367			D = BB - 4A*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
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
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			}