Index: trunk/psLib/src/math/psStats.c
===================================================================
--- trunk/psLib/src/math/psStats.c	(revision 23988)
+++ trunk/psLib/src/math/psStats.c	(revision 24092)
@@ -162,5 +162,6 @@
 static psF32 minimizeLMChi2Gauss1D(psVector *deriv, const psVector *params, const psVector *coords);
 // static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
-static psF32 fitQuadraticSearchForYThenReturnXusingValues(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
+// static psF32 fitQuadraticSearchForYThenReturnXusingValues(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
+static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
 
 /******************************************************************************
@@ -822,9 +823,12 @@
         PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints/2.0, 0);
 
-        psTrace(TRACE, 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian,
-                cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]);
-
-        // ADD step 3: Interpolate to the exact 50% position: this is the robust histogram median.
-        stats->robustMedian = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
+        psTrace(TRACE, 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian, cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]);
+
+        // ADD step 3: Interpolate to the exact 50% position in bin units
+	stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
+	// float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
+	// fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian);
+
+	// convert bin to bin value: this is the robust histogram median.
         if (isnan(stats->robustMedian)) {
             psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to fit a quadratic and calculate the 50-percent position.\n");
@@ -971,8 +975,8 @@
     // ADD step 8: Interpolate to find these two positions exactly: these are the upper and lower quartile
     // positions.
-    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f);
-    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f);
+    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f);
+    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f);
     if (isnan(binLo25F32) || isnan(binHi25F32)) {
-        psLogMsg(TRACE, PS_LOG_DETAIL, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n");
+        psLogMsg(TRACE, PS_LOG_DETAIL, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnBin() returned a NAN.\n");
         goto escape;
     }
@@ -2737,5 +2741,4 @@
     return tmpFloat;
 }
-# endif
 
 /******************************************************************************
@@ -2865,4 +2868,135 @@
     return tmpFloat;
 }
+# endif
+
+/******************************************************************************
+fitQuadraticSearchForYThenReturnXusingValues(*xVec, *yVec, binNum, yVal): A general routine
+which fits a quadratic to three points and returns the x bin value corresponding to the input
+y-value.  This routine takes psVectors of x/y pairs as input, and fits a quadratic to the 3
+points surrounding element binNum in the vectors.  This version uses the values of x[i] for the
+x coordinates (not the midpoints).  This is appropriate for a cumulative histogram.  It then
+determines for what value x does that quadratic f(x) = yVal (the input parameter).
+
+XXX this function is used a fair amount in an inner loop: the polynomial fitting and evaluation
+could easily be done with statically allocated doubles, skipping the psLib versions of
+polynomial fitting, etc.
+
+*****************************************************************************/
+static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec,
+						 psVector *yVec,
+						 psS32 binNum,
+						 psF32 yVal
+    )
+{
+    psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
+    if (psTraceGetLevel("psLib.math") >= 8) {
+        PS_VECTOR_PRINT_F32(xVec);
+        PS_VECTOR_PRINT_F32(yVec);
+    }
+
+    PS_ASSERT_VECTOR_NON_NULL(xVec, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(yVec, NAN);
+    PS_ASSERT_VECTOR_TYPE(xVec, PS_TYPE_F32, NAN);
+    PS_ASSERT_VECTOR_TYPE(yVec, PS_TYPE_F32, NAN);
+    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(xVec->n - 1), NAN);
+    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN);
+
+    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
+    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
+    psF32 tmpFloat = 0.0f;
+
+    if ((binNum >= 1) && (binNum <= (yVec->n - 2)) && (binNum <= (xVec->n - 2))) {
+        // The general case.  We have all three points.
+        x->data.F64[0] = binNum - 1;
+        x->data.F64[1] = binNum;
+        x->data.F64[2] = binNum + 1;
+        y->data.F64[0] = yVec->data.F32[binNum - 1];
+        y->data.F64[1] = yVec->data.F32[binNum];
+        y->data.F64[2] = yVec->data.F32[binNum + 1];
+        psTrace(TRACE, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
+        psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
+        psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
+
+        // Ensure that the y value lies within range of the y values.
+        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
+               ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    _("Specified yVal, %g, is not within y-range, %g to %g."),
+                    (psF64)yVal, y->data.F64[0], y->data.F64[2]);
+            return NAN;
+        }
+
+        // Ensure that the y values are monotonic.
+        if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) ||
+            ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
+            psError(PS_ERR_UNKNOWN, true,
+                    "This routine must be called with monotonically increasing or decreasing data points.\n");
+            psFree(x);
+            psFree(y);
+            return NAN;
+        }
+
+        // Determine the coefficients of the polynomial.
+        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
+        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) {
+            psError(PS_ERR_UNEXPECTED_NULL, false,
+                    _("Failed to fit a 1-dimensional polynomial to the three specified data points.  "
+                      "Returning NAN."));
+            psFree(x);
+            psFree(y);
+            return NAN;
+        }
+
+        psTrace(TRACE, 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
+        psTrace(TRACE, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
+        psTrace(TRACE, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
+        psTrace(TRACE, 6, "Fitted y vec is (%f %f %f)\n",
+                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
+                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
+                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
+
+        psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
+        float binValue = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[2]);
+        psFree(myPoly);
+
+        if (isnan(binValue)) {
+            psError(PS_ERR_UNEXPECTED_NULL,
+                    false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
+            psFree(x);
+            psFree(y);
+            return(NAN);
+        }
+
+	// I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1
+	assert (binValue >= binNum - 1);
+	assert (binValue <= binNum + 1);
+
+	int fitBin = binValue;
+	float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
+	float dY = binValue - fitBin;
+	tmpFloat = xVec->data.F32[fitBin] + dY * dX;
+    } else {
+        // These are special cases where the bin is at the beginning or end of the vector.
+        if (binNum == 0) {
+            // We have two points only at the beginning of the vectors x and y.
+	    // X = (dX/dY)(Y - Yo) + Xo
+	    float dX = xVec->data.F32[1] - xVec->data.F32[0];
+	    float dY = yVec->data.F32[1] - yVec->data.F32[0];	    
+	    tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0];
+        } else if (binNum == (xVec->n - 1)) {
+            // We have two points only at the end of the vectors x and y.
+	    // X = (dX/dY)(Y - Yo) + Xo
+	    float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1];
+	    float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1];	    
+	    tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1];
+        } 
+    }
+
+    psTrace(TRACE, 6, "FIT: return %f\n", tmpFloat);
+    psFree(x);
+    psFree(y);
+
+    return tmpFloat;
+}
 
 /******************************************************************************