Changeset 6305 for trunk/psLib/src/math/psStats.c

Timestamp:

Feb 2, 2006, 11:09:08 AM (20 years ago)

Author:

gusciora

Message:

This file contains some of the generic math functions which
were previously in psStats and psSpline and psMinimize.

File:

• trunk/psLib/src/math/psStats.c (modified) (47 diffs)

trunk/psLib/src/math/psStats.c

@@ -3 +3 @@
  *  @ingroup Stat
  *
- *  This file will hold the definition of the histogram and stats data
+ *  This file will hold the definitions of the histogram and stats data
  *  structures.  It also contains prototypes for procedures which operate
  *  on those data structures.
@@ -10 +10 @@
  *
  *  XXX: The following stats members are never used, or set in this code.
- *      stats->robustN50
  *      stats->clippedNvalues
- *      stats->binsize
  *
  * XXX: Must do
  * nSubsample points
  * use ->min and ->max (PS_STAT_USE_RANGE)
- * use ->binsize (PS_STAT_USE_BINSIZE)
  *
- *
- *
- *
- *  @version $Revision: 1.164 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-01 02:07:24 $
+ *  @version $Revision: 1.165 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-02 21:09:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -49 +43 @@
 #include "psPolynomial.h"
 #include "psConstants.h"
+#include "psMathUtils.h"
 
 #include "psErrorText.h"
@@ -61 +56 @@
 #define PS_CLIPPED_SIGMA_LB 1.0
 #define PS_CLIPPED_SIGMA_UB 10.0
-#define PS_POLY_MEDIAN_MAX_ITERATIONS 10
+#define PS_POLY_MEDIAN_MAX_ITERATIONS 30
 
 #define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
@@ -847 +842 @@
 
     // Calculate the quartile points exactly.
-    // XXX: We should probably do a bin-midpoint if the quartile is not an
-    // integer.
     stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
     stats->sampleLQ = sortedVector->data.F32[nValues / 4];
@@ -859 +852 @@
 }
 
-/******************************************************************************
-p_psVectorSampleStdevOLD(myVector, maskVector, maskVal, stats): calculates the
-stdev of the input vector.
-Inputs
-    myVector
-    maskVector
-    maskVal
-    stats
-Returns
-    NULL
- 
-XXX: remove this
- *****************************************************************************/
-void p_psVectorSampleStdevOLD(const psVector* myVector,
-                              const psVector* errors,
-                              const psVector* maskVector,
-                              psU32 maskVal,
-                              psStats* stats)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    psS32 i = 0;                  // Loop index variable
-    psS32 countInt = 0;           // # of data points being used
-    psF32 countFloat = 0.0;     // # of data points being used
-    psF32 mean = 0.0;           // The mean
-    psF32 diff = 0.0;           // Used in calculating stdev
-    psF32 sumSquares = 0.0;     // temporary variable
-    psF32 sumDiffs = 0.0;       // temporary variable
-
-    // This procedure requires the mean.  If it has not been already
-    // calculated, then call p_psVectorSampleMean()
-    if (0 != isnan(stats->sampleMean)) {
-        p_psVectorSampleMean(myVector, errors, maskVector, maskVal, stats);
-    }
-    // If the mean is NAN, then generate a warning and set the stdev to NAN.
-    if (0 != isnan(stats->sampleMean)) {
-        stats->sampleStdev = NAN;
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): p_psVectorSampleMean() reported a NAN mean.\n");
-        psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-        return;
-    }
-
-    mean = stats->sampleMean;
-    if (stats->options & PS_STAT_USE_RANGE) {
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i]) &&
-                        (stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    diff = myVector->data.F32[i] - mean;
-                    sumSquares += (diff * diff);
-                    sumDiffs += diff;
-                    countInt++;
-                }
-            }
-        } else {
-            for (i = 0; i < myVector->n; i++) {
-                if ((stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    diff = myVector->data.F32[i] - mean;
-                    sumSquares += (diff * diff);
-                    sumDiffs += diff;
-                    countInt++;
-                }
-            }
-            countInt = myVector->n;
-        }
-    } else {
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    diff = myVector->data.F32[i] - mean;
-                    sumSquares += (diff * diff);
-                    sumDiffs += diff;
-                    countInt++;
-                }
-            }
-        } else {
-            for (i = 0; i < myVector->n; i++) {
-                diff = myVector->data.F32[i] - mean;
-                sumSquares += (diff * diff);
-                sumDiffs += diff;
-                countInt++;
-            }
-            countInt = myVector->n;
-        }
-    }
-    if (countInt == 0) {
-        stats->sampleStdev = NAN;
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): no valid psVector elements (%d).  Setting stats->sampleStdev = NAN.\n", countInt);
-    } else if (countInt == 1) {
-        stats->sampleStdev = 0.0;
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d).  Setting stats->sampleStdev = 0.0.\n", countInt);
-    } else {
-        countFloat = (psF32)countInt;
-        stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
-    }
-    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-}
 /******************************************************************************
 p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
@@ -1094 +989 @@
     -2: warning
  
-XXX: Do we really need to calculate median and mean?
- 
 XXX: Use static vectors for tmpMask.
- 
-XXX: This has not been tested.
  *****************************************************************************/
 psS32 p_psVectorClippedStats(const psVector* myVector,
@@ -1223 +1114 @@
         // Otherwise, use the new results and continue.
         if (isnan(statsTmp->sampleMean) || isnan(statsTmp->sampleStdev)) {
-            // Exit loop.  XXX: throw an error/warning here?
+            // Exit loop.
             iter = stats->clipIter;
-            rc = -1;
+            psError(PS_ERR_UNKNOWN, true, "p_psVectorSampleMean() or p_psVectorSampleStdev() returned a NAN.\n");
+            psFree(tmpMask);
+            return(-1);
         } else {
             clippedMean = statsTmp->sampleMean;
@@ -1246 +1139 @@
 }
 
-/*****************************************************************************
-These macros and functions define the following functions:
- 
-    p_psNormalizeVectorRange(myData, low, high)
- 
-That assumes that the low/high arguments are PS_TYPE_F64; the vector myData
-can be of any type.  Arguments low/high will be converted to the appropriate
-type and one of the type-specific functions below will be called:
- 
-    p_psNormalizeVectorRangeU8(myData, low, high)
-    p_psNormalizeVectorRangeU16(myData, low, high)
-    p_psNormalizeVectorRangeU32(myData, low, high)
-    p_psNormalizeVectorRangeU64(myData, low, high)
-    p_psNormalizeVectorRangeS8(myData, low, high)
-    p_psNormalizeVectorRangeS16(myData, low, high)
-    p_psNormalizeVectorRangeS32(myData, low, high)
-    p_psNormalizeVectorRangeS64(myData, low, high)
-    p_psNormalizeVectorRangeF32(myData, low, high)
-    p_psNormalizeVectorRangeF64(myData, low, high)
- *****************************************************************************/
-#define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \
-void p_psNormalizeVectorRange##TYPE(psVector* myData, \
-                                    ps##TYPE outLow, \
-                                    ps##TYPE outHigh) \
-{ \
-    ps##TYPE min = (ps##TYPE) PS_MAX_##TYPE; \
-    ps##TYPE max = (ps##TYPE) -PS_MAX_##TYPE; \
-    psS32 i = 0; \
-    \
-    for (i = 0; i < myData->n; i++) { \
-        if (myData->data.TYPE[i] < min) { \
-            min = myData->data.TYPE[i]; \
-        } \
-        if (myData->data.TYPE[i] > max) { \
-            max = myData->data.TYPE[i]; \
-        } \
-    } \
-    \
-    /* Ensure that max!=min before we divide by (max-min) */ \
-    if (max != min) { \
-        for (i = 0; i < myData->n; i++) { \
-            myData->data.TYPE[i] = (outLow + (myData->data.TYPE[i] - min) * \
-                                    (outHigh - outLow) / (max - min)); \
-        } \
-    } else { \
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: (max==min).  Setting all elements to min.\n"); \
-        for (i = 0; i < myData->n; i++) \
-        { \
-            \
-            myData->data.TYPE[i] = outLow; \
-            \
-        } \
-    } \
-} \
-
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U8)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U16)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U32)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U64)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S8)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S16)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S32)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S64)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F32)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64)
-
-void p_psNormalizeVectorRange(psVector* myData,
-                              psF64 outLow,
-                              psF64 outHigh)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    switch (myData->type.type) {
-    case PS_TYPE_U8:
-        p_psNormalizeVectorRangeU8(myData, (psU8) outLow, (psU8) outHigh);
-        break;
-    case PS_TYPE_U16:
-        p_psNormalizeVectorRangeU16(myData, (psU16) outLow, (psU16) outHigh);
-        break;
-    case PS_TYPE_U32:
-        p_psNormalizeVectorRangeU32(myData, (psU32) outLow, (psU32) outHigh);
-        break;
-    case PS_TYPE_U64:
-        p_psNormalizeVectorRangeU64(myData, (psU64) outLow, (psU64) outHigh);
-        break;
-    case PS_TYPE_S8:
-        p_psNormalizeVectorRangeS8(myData, (psS8) outLow, (psS8) outHigh);
-        break;
-    case PS_TYPE_S16:
-        p_psNormalizeVectorRangeS16(myData, (psS16) outLow, (psS16) outHigh);
-        break;
-    case PS_TYPE_S32:
-        p_psNormalizeVectorRangeS32(myData, (psS32) outLow, (psS32) outHigh);
-        break;
-    case PS_TYPE_S64:
-        p_psNormalizeVectorRangeS64(myData, (psS64) outLow, (psS64) outHigh);
-        break;
-    case PS_TYPE_F32:
-        p_psNormalizeVectorRangeF32(myData, (psF32) outLow, (psF32) outHigh);
-        break;
-    case PS_TYPE_F64:
-        p_psNormalizeVectorRangeF64(myData, (psF64) outLow, (psF64) outHigh);
-        break;
-    case PS_TYPE_C32:
-    case PS_TYPE_C64:
-    default:
-        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
-                "Unallowable operation: %s has incorrect type.",
-                myData);
-        break;
-    }
-    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-}
+
+
 
 /******************************************************************************
@@ -1371 +1154 @@
 XXX: Create a 2nd-order polynomial version and solve for X analytically.
  *****************************************************************************/
-psF32 p_ps1DPolyMedian(psPolynomial1D* myPoly,
-                       psF32 rangeLow,
-                       psF32 rangeHigh,
-                       psF32 getThisValue)
+psF32 p_ps1DPolyMedian(
+    psPolynomial1D* myPoly,
+    psF32 rangeLow,
+    psF32 rangeHigh,
+    psF32 getThisValue)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
@@ -1405 +1189 @@
 
         f = psPolynomial1DEval(myPoly, midpoint);
+        printf("p_ps1DPolyMedian(): f(%f) is %f\n", midpoint, f);
         if (fabs(f - getThisValue) <= FLT_EPSILON) {
             return (midpoint);
@@ -1419 +1204 @@
     return (midpoint);
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#define PS_SQRT(x) sqrt(x)
+static psF32 QuadraticInverse(
+    psF32 a,
+    psF32 b,
+    psF32 c,
+    psF32 y)
+{
+    psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a));
+
+    psF64 x1 = -b/(2.0*a) + tmp;
+    psF64 x2 = -b/(2.0*a) - tmp;
+
+    psF64 y1 = (a * x1 * x1) + (b * x1) + c;
+    psF64 y2 = (a * x2 * x2) + (b * x2) + c;
+
+    printf("QuadraticInverse: %fx^2 + %fx + %f\n", a, b, c);
+    printf("QuadraticInverse: y is %f\n", y);
+    printf("QuadraticInverse: (x1, x2) is (%f %f)\n", x1, x2);
+    printf("QuadraticInverse: (y1, y2) is (%f %f)\n", y1, y2);
+
+    return(x1);
+}
+
 
 /******************************************************************************
@@ -1438 +1284 @@
 arbitrary vectors.  We should probably test that condition.
 *****************************************************************************/
-psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec,
-                                        psVector *yVec,
-                                        psS32 binNum,
-                                        psF32 yVal)
+psF32 fitQuadraticSearchForYThenReturnX(
+    psVector *xVec,
+    psVector *yVec,
+    psS32 binNum,
+    psF32 yVal)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
@@ -1461 +1308 @@
     psVector *y = psVectorAlloc(3, PS_TYPE_F64);
     psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
-    psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
 
     psF32 tmpFloat = 0.0f;
@@ -1489 +1335 @@
             if (!(y->data.F64[1] <= y->data.F64[2])) {
                 psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n");
-                psFree(myPoly);
                 psFree(x);
                 psFree(y);
@@ -1505 +1350 @@
             if (!(y->data.F64[1] >= y->data.F64[2])) {
                 psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n");
-                psFree(myPoly);
                 psFree(x);
                 psFree(y);
@@ -1525 +1369 @@
 
         // Determine the coefficients of the polynomial.
-        //        myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
+        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
         myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x);
         if (myPoly == NULL) {
@@ -1531 +1375 @@
                     false,
                     PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLYNOMIAL_1D_FIT);
-            psFree(myPoly);
             psFree(x);
             psFree(y);
@@ -1541 +1384 @@
         psTrace(__func__, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
         psTrace(__func__, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
-        psTrace(__func__, 6, "Fitted y vec is (%f %f %f)\n", (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
+        psTrace(__func__, 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]));
@@ -1548 +1392 @@
         // polynomial, looking for the value x such that f(x) = yVal
         psTrace(__func__, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
+        printf("Calling p_ps1DPolyMedian() for the y-value (%.2f) that has an x in the range (%.2f - %.2f)\n", yVal, x->data.F64[0], x->data.F64[2]);
         tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal);
+        printf("Cool, tmpFloat is %.2f\n", tmpFloat);
+
+        QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal);
+        psFree(myPoly);
+
 
         if (isnan(tmpFloat)) {
@@ -1554 +1404 @@
                     false,
                     PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN);
-            psFree(myPoly);
             psFree(x);
             psFree(y);
@@ -1581 +1430 @@
 
     psTrace(__func__, 6, "FIT: return %f\n", tmpFloat);
-    psFree(myPoly);
     psFree(x);
     psFree(y);
@@ -1591 +1439 @@
 
 /******************************************************************************
-XXX: We assume unnormalized gaussians.
+NOTE: We assume unnormalized gaussians.
  *****************************************************************************/
-psF32 psMinimizeLMChi2Gauss1D(psVector *deriv,
-                              const psVector *params,
-                              const psVector *coords)
+psF32 psMinimizeLMChi2Gauss1D(
+    psVector *deriv,
+    const psVector *params,
+    const psVector *coords)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
@@ -1627 +1476 @@
 }
 
+//XXX: Use the psLib function?
 psF32 LinInterpolate(psF32 x0,
                      psF32 x1,
@@ -1643 +1493 @@
     return(x0 + y * (x1 - x0) / (y1 - y0));
 }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
 
 
@@ -1754 +1532 @@
 XXX: Review and ensure that all memory is free'ed at premature exits.
 *****************************************************************************/
-#define INITIAL_NUM_BINS 1000.0
+#define INITIAL_NUM_BINS 500.0
 psS32 p_psVectorRobustStats(const psVector* myVector,
                             const psVector* errors,
@@ -1770 +1548 @@
     psHistogram *cumulativeRobustHistogram = NULL;
     psS32 numBins = 0;
-    psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32);
-    tmpScalar->type.type = PS_TYPE_F32;
+    psScalar tmpScalar;
+    tmpScalar.type.type = PS_TYPE_F32;
     psS32 totalDataPoints = 0;
     psS32 rc = 0;
@@ -1806 +1584 @@
                 psFree(tmpStatsMinMax);
                 psFree(tmpMaskVec);
-                psFree(tmpScalar);
                 psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
                 return(1);
@@ -1833 +1610 @@
             psFree(tmpStatsMinMax);
             psFree(tmpMaskVec);
-            psFree(tmpScalar);
 
             psTrace(__func__, 4, "---- %s(0) end  ----\n", __func__);
@@ -1885 +1661 @@
             binMedian = 0;
         } else {
-            tmpScalar->data.F32 = totalDataPoints/2.0;
-            binMedian = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+            tmpScalar.data.F32 = totalDataPoints/2.0;
+            binMedian = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
             if (binMedian < 0) {
                 psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50 precent data point (%d).\n", binMedian);
@@ -1892 +1668 @@
                 psFree(robustHistogram);
                 psFree(cumulativeRobustHistogram);
-                psFree(tmpScalar);
                 psFree(tmpMaskVec);
                 psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
@@ -1905 +1680 @@
         // XXX: Check return codes.
         //
+        printf("ADD: step 3: Interpolate to the exact 50-percent position: this is the robust histogram median.\n");
         stats->robustMedian = fitQuadraticSearchForYThenReturnX(
                                   *(psVector* *)&cumulativeRobustHistogram->bounds,
@@ -1918 +1694 @@
         psS32 binLo = 0;
         psS32 binHi = 0;
-        tmpScalar->data.F32 = totalDataPoints * 0.158655f;
-        if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
+        tmpScalar.data.F32 = totalDataPoints * 0.158655f;
+        if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
             binLo = 0;
         } else {
             // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
-            binLo = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+            binLo = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
             if (binLo > cumulativeRobustHistogram->nums->n-1) {
                 binLo = cumulativeRobustHistogram->nums->n-1;
             }
         }
-        tmpScalar->data.F32 = totalDataPoints * 0.841345f;
-        if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
+        tmpScalar.data.F32 = totalDataPoints * 0.841345f;
+        if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
             binHi = 0;
         } else {
             // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
-            binHi = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+            binHi = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
             if (binHi > cumulativeRobustHistogram->nums->n-1) {
                 binHi = cumulativeRobustHistogram->nums->n-1;
@@ -1947 +1723 @@
             psFree(robustHistogram);
             psFree(cumulativeRobustHistogram);
-            psFree(tmpScalar);
             psFree(tmpMaskVec);
             psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
@@ -2014 +1789 @@
         if (sigma < (2 * binSize)) {
             psTrace(__func__, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
-            psF32 medianLo;
-            psF32 medianHi;
-            if ((binMedian - 25) > 0) {
-                medianLo = robustHistogram->bounds->data.F32[binMedian - 25];
-            } else {
-                medianLo = robustHistogram->bounds->data.F32[0];
-            }
-            if ((binMedian + 25) < robustHistogram->bounds->n) {
-                medianHi = robustHistogram->bounds->data.F32[binMedian + 25];
-            } else {
-                medianHi = robustHistogram->bounds->data.F32[robustHistogram->bounds->n - 1];
-            }
-
+            psF32 medianLo = robustHistogram->bounds->data.F32[PS_MAX(0, (binMedian - 25))];
+            psF32 medianHi = robustHistogram->bounds->data.F32[PS_MIN(robustHistogram->bounds->n - 1, (binMedian + 25))];
             psTrace(__func__, 6, "Masking data more than 25 bins from the median (%.2f)\n", stats->robustMedian);
             psTrace(__func__, 6, "The median is at bin number %d.  We mask bins outside the bin range (%d:%d)\n", binMedian, binMedian - 25, binMedian + 25);
@@ -2052 +1816 @@
     psS32 binHi25 = 0;
 
-    tmpScalar->data.F32 = totalDataPoints * 0.25f;
-    if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
+    tmpScalar.data.F32 = totalDataPoints * 0.25f;
+    if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
         // XXX: Special case where its in first bin.  Must code last bin.
         binLo25 = 0;
     } else {
-        binLo25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
-    }
-    tmpScalar->data.F32 = totalDataPoints * 0.75f;
-    if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
+        binLo25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
+    }
+    tmpScalar.data.F32 = totalDataPoints * 0.75f;
+    if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
         // XXX: Special case where its in first bin.  Must code last bin.
         binHi25 = 0;
     } else {
-        binHi25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+        binHi25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
     }
     if ((binLo25 < 0) || (binHi25 < 0)) {
@@ -2071 +1835 @@
         psFree(robustHistogram);
         psFree(cumulativeRobustHistogram);
-        psFree(tmpScalar);
         psFree(tmpMaskVec);
         psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
@@ -2084 +1847 @@
     // XXX: Check for errors.
     //
+    printf("ADD: step 8: Interpolate for the lower quartile positions.\n");
     psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(
                            *(psVector* *)&cumulativeRobustHistogram->bounds,
@@ -2089 +1853 @@
                            binLo25,
                            totalDataPoints * 0.25f);
+    printf("ADD: step 8: Interpolate for the upper quartile positions.\n");
     psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(
                            *(psVector* *)&cumulativeRobustHistogram->bounds,
@@ -2133 +1898 @@
             psFree(robustHistogram);
             psFree(cumulativeRobustHistogram);
-            psFree(tmpScalar);
             psFree(tmpMaskVec);
             psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
@@ -2170 +1934 @@
         psS32 binMin = 0;
         psS32 binMax = 0;
-        tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma);
-        if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) {
+        tmpScalar.data.F32 = stats->robustMedian - (2.0 * sigma);
+        if (tmpScalar.data.F32 <= newHistogram->bounds->data.F32[0]) {
             binMin = 0;
         } else {
-            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
-        }
-
-        tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma);
-        if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) {
+            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
+        }
+
+        tmpScalar.data.F32 = stats->robustMedian + (2.0 + sigma);
+        if (tmpScalar.data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) {
             binMax = newHistogram->bounds->n-1;
         } else {
-            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
+            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
         }
         if ((binMin < 0) || (binMax < 0)) {
@@ -2207 +1971 @@
         // histogram median.
         //
-        tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma);
-        if (tmpScalar->data.F32 < tmpStatsMinMax->min) {
+        tmpScalar.data.F32 = stats->robustMedian - (20.0 * sigma);
+        if (tmpScalar.data.F32 < tmpStatsMinMax->min) {
             binMin = 0;
         } else {
-            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
+            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
             // XXX: check for errors here.
         }
-        tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma);
-        if (tmpScalar->data.F32 > tmpStatsMinMax->max) {
+        tmpScalar.data.F32 = stats->robustMedian + (20.0 * sigma);
+        if (tmpScalar.data.F32 > tmpStatsMinMax->max) {
             binMax = newHistogramSmoothed->n - 1;
         } else {
-            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
+            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
             // XXX: check for errors here.
         }
@@ -2275 +2039 @@
             psFree(robustHistogram);
             psFree(cumulativeRobustHistogram);
-            psFree(tmpScalar);
             psFree(newHistogram);
             psFree(x);
@@ -2314 +2077 @@
     psFree(tmpStatsMinMax);
     psFree(cumulativeRobustHistogram);
-    psFree(tmpScalar);
     psFree(tmpMaskVec);
     psFree(robustHistogram);
@@ -2321 +2083 @@
     return(0);
 }
-
-
-
-
-
-
-
-
 
 /*****************************************************************************/
@@ -2362 +2116 @@
     newStruct->robustUQ = NAN;
     newStruct->robustLQ = NAN;
-    newStruct->robustN50 = -1;            // XXX: This is never used
+    newStruct->robustN50 = -1;
     newStruct->fittedMean = NAN;
     newStruct->fittedStdev = NAN;
@@ -2373 +2127 @@
     newStruct->min = NAN;
     newStruct->max = NAN;
-    newStruct->binsize = NAN;          // XXX: This is never used
+    newStruct->binsize = NAN;
     newStruct->nSubsample = 100000;
     newStruct->options = options;