Changeset 2197 for trunk/psLib/src/dataManip/psStats.c

Timestamp:

Oct 26, 2004, 11:25:59 AM (22 years ago)

Author:

gusciora

Message:

...

File:

• trunk/psLib/src/dataManip/psStats.c (modified) (29 diffs)

trunk/psLib/src/dataManip/psStats.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.67 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-13 02:40:13 $
+ *  @version $Revision: 1.68 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-26 21:24:43 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -40 +40 @@
 /*****************************************************************************/
 // will use robust statistical methods.
-#define GAUSS_WIDTH 20                     // The width of the Gaussian or boxcar smoothing.
+#define GAUSS_WIDTH 5       // The width of the Gaussian or boxcar smoothing.
 #define CLIPPED_NUM_ITER_LB 1
 #define CLIPPED_NUM_ITER_UB 10
@@ -57 +57 @@
 /* TYPE DEFINITIONS                                                          */
 /*****************************************************************************/
+psVector* p_psConvertToF32(psVector* in);
 
 /*****************************************************************************/
@@ -423 +424 @@
     float rangeMax = 0.0;       // Exclude date above this
 
-    // Determine if the number of data points exceed a threshold which will
-    // cause to generate robust stats, as opposed to exact stats.
-    // XXX: This code is no longer used.  We now calculate the median exactly
-    // regardless of the vector size.
-    /*
-     * // Calculate the robust quartiles. stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
-     * p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
-     * 
-     * // Store the robust quartiles into the sample quartile members. stats->sampleMedian =
-     * stats2->robustMedian;
-     * 
-     * // Free temporary data buffers. psFree(stats2);
-     * 
-     * // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. stats->options = stats->options |
-     * PS_STAT_ROBUST_FOR_SAMPLE;
-     * 
-     * return; } */
-
     // Determine how many data points fit inside this min/max range
     // and are not masked, IF the maskVector is not NULL>
@@ -507 +490 @@
 
 /******************************************************************************
-    This routine smoothes the data in the input robustHistogram with a
-    Gaussian of width sigma.
- 
-    XXX: Must consult with IfA on the proper values for the Gaussian
-    smoothing.  Currently, the Gaussian coefficients are such that only
-    the center coefficient is non-zero.  This is because "e" is being
-    raised to a very large negative number for all points except the
-    center point.
- 
-    XXX: use a static variable for gaussianCoefs[] and compute them once.
- *****************************************************************************/
-psVector* p_psVectorsmoothHistGaussian(psHistogram* robustHistogram,
+p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input
+robustHistogram with a Gaussian of width sigma.
+ *****************************************************************************/
+#define GAUSS_WIDTH_OTHER 5.0
+psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram,
                                        float sigma)
 {
     int i = 0;                  // Loop index variable
     int j = 0;                  // Loop index variable
-    float denom = 0.0;          // Temporary variable
-    float expo = 0.0;           // Temporary variable
-    float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients
-    psVector* smooth = psVectorAlloc(robustHistogram->nums->n, PS_TYPE_F32);
-
-    for (i = 0; i < (1 + (2 * GAUSS_WIDTH)); i++) {
-        if (fabs(sigma) >= FLT_EPSILON) {
-            // If sigma does not equal zero, then we use Gaussian smoothing.
-            #ifdef  DARWIN
-            denom = (float)sqrt(2.0 * M_PI * sigma * sigma);
-            #else
-
-            denom = sqrtf(2.0 * M_PI * sigma * sigma);
-            #endif
-
-            expo = -(float)((i - GAUSS_WIDTH) * (i - GAUSS_WIDTH));
-            expo /= (2.0 * sigma * sigma);
-            gaussianCoefs[i] = exp(expo / denom);
-
-            //printf("p_psVectorsmoothHistGaussian(): gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
-            // NOTE: Gaussian smoothing just isn't working with low sigma
-            // values.  The problem is that the Gaussian coefficients are
-            // all zero, except for the middle coefficient, which is exactly
-            // one.  Therefore, I'm using boxcar smoothing.
-            gaussianCoefs[i] = 1.0 / (1.0 + (2.0 * (float)GAUSS_WIDTH));
-            // printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
-        } else {
-            /* If sigma equals zero (all pixels have the same value) the above code will divide by zero.
-             * Therefore, we don't need to smooth the data. */
-            for (i = 0; i < robustHistogram->nums->n; i++) {
-                smooth->data.F32[i] = (float)robustHistogram->nums->data.S32[i];
-            }
-            return (smooth);
-        }
-    }
-
-    // Perform the actual smoothing.
-    for (i = 0; i < robustHistogram->nums->n; i++) {
+    float iMid;
+    float jMid;
+    int numBins = robustHistogram->nums->n;
+    int numBounds = robustHistogram->bounds->n;
+    psVector* smooth = psVectorAlloc(numBins, PS_TYPE_F32);
+    int jMin = 0;
+    int jMax = 0;
+    float firstBound = robustHistogram->bounds->data.F32[0];
+    float lastBound = robustHistogram->bounds->data.F32[numBounds-1];
+    psScalar x;
+
+    x.type.type = PS_TYPE_F32;
+    for (i = 0; i < numBins; i++) {
+        iMid = (robustHistogram->bounds->data.F32[i] +
+                robustHistogram->bounds->data.F32[i+1]) / 2.0;
+
+        // YYY: The p_psVectorBinDisect() routine does much of the work of
+        // the following conditionals, however, it also reports a warning
+        // message.  I don't want the warning message so I reproduce the
+        // conditionals here.  Maybe p_psVectorBinDisect() should not produce
+        // warnings?
+
+        x.data.F32 = iMid - (GAUSS_WIDTH * sigma);
+        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
+            jMin = p_psVectorBinDisect(robustHistogram->bounds, &x);
+        } else if (x.data.F32 <= firstBound) {
+            jMin = 0;
+        } else if (x.data.F32 >= lastBound) {
+            jMin = robustHistogram->bounds->n - 1;
+
+        }
+
+        x.data.F32 = iMid + (GAUSS_WIDTH * sigma);
+        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
+            jMax = p_psVectorBinDisect(robustHistogram->bounds, &x);
+        } else if (x.data.F32 <= firstBound) {
+            jMax = 0;
+        } else if (x.data.F32 >= lastBound) {
+            jMax = robustHistogram->bounds->n - 1;
+
+        }
+
         smooth->data.F32[i] = 0.0;
-        for (j = -GAUSS_WIDTH; j <= +GAUSS_WIDTH; j++) {
-            if (((j + i) >= 0) && ((j + i) < smooth->n)) {
-                smooth->data.F32[i] += (gaussianCoefs[j + GAUSS_WIDTH] *
-                                        (float)robustHistogram->nums->data.S32[j + i]);
-            }
-        }
-    }
-    return (smooth);
+        for (j = jMin ; j <= jMax ; j++) {
+            jMid = (robustHistogram->bounds->data.F32[j] +
+                    robustHistogram->bounds->data.F32[j+1]) / 2.0;
+            smooth->data.F32[i] +=
+                ((float) robustHistogram->nums->data.U32[j]) *
+                psGaussian(jMid, iMid, sigma, true);
+        }
+    }
+
+    return(smooth);
 }
 
@@ -913 +895 @@
 specified data point.
  *****************************************************************************/
-float p_psGaussian(const psVector* restrict myData, const psVector* restrict myParams)
+float p_psGaussian(const psVector* restrict myData,
+                   const psVector* restrict myParams)
 {
     float x = myData->data.F32[0];
@@ -930 +913 @@
 calculates the specified partial derivative of the above Gaussian function.
  *****************************************************************************/
-float p_psGaussianDeriv(const psVector* restrict myData, const psVector* restrict myParams, int whichParam)
+float p_psGaussianDeriv(const psVector* restrict myData,
+                        const psVector* restrict myParams, int whichParam)
 {
     float x = myData->data.F32[0];
@@ -1036 +1020 @@
 
 /******************************************************************************
-p_psFitQuadratic(robustHistogram. cumulativeSums, binNum, fitFloat): given
-the specified histogram, with the specified pre-calculated cumulativeSums,
-this routine fits a quadratic f(x) to the 3 bins surrounding bin "binNum",
-and then finds the point x such that f(x) == fitFloat.
+fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): This routine
+takes psVectors of x/y pairs as input, and fits a quadratic to the 3 points
+surrounding element binNum in the vectors (the midpoint between element i
+and i+1 is used for x[i]).  It then determines for what value x does that
+quadratic f(x) = yVal (the input parameter).
  
-XXX: This function is currently not being used.
- *****************************************************************************/
-float p_psFitQuadratic(psHistogram* histogram,
-                       psVector* cumulativeSums,
-                       int binNum,
-                       float fitFloat)
-{
+// XXX: Use static variables.
+*****************************************************************************/
+float fitQuadraticSearchForYThenReturnX(psVector *xVec,
+                                        psVector *yVec,
+                                        int binNum,
+                                        float yVal)
+{
+    //    static psVector* x = NULL;
+    //    static psVector* y = NULL;
+    //    static psVector* yErr = NULL;
     psVector* x = psVectorAlloc(3, PS_TYPE_F64);
     psVector* y = psVectorAlloc(3, PS_TYPE_F64);
     psVector* yErr = psVectorAlloc(3, PS_TYPE_F64);
+
+    /*
+        if (x == NULL) {
+            x = psVectorAlloc(3, PS_TYPE_F64);
+            p_psMemSetPersistent(x, true);
+        }
+        if (y == NULL) {
+            y = psVectorAlloc(3, PS_TYPE_F64);
+            p_psMemSetPersistent(y, true);
+        }
+        if (yErr == NULL) {
+            yErr = psVectorAlloc(3, PS_TYPE_F64);
+            p_psMemSetPersistent(yErr, true);
+        }
+    */
     psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
-
-    if ((binNum > 0) && (binNum < (histogram->nums->n + 1))) {
-        x->data.F64[0] = (double)0.5 *
-                         (histogram->bounds->data.F32[binNum - 1] + histogram->bounds->data.F32[binNum]);
-        x->data.F64[1] = (double)0.5 *
-                         (histogram->bounds->data.F32[binNum] + histogram->bounds->data.F32[binNum + 1]);
-        x->data.F64[2] = (double)0.5 *
-                         (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum + 2]);
-
-        y->data.F64[0] = cumulativeSums->data.F32[binNum - 1];
-        y->data.F64[1] = cumulativeSums->data.F32[binNum];
-        y->data.F64[2] = cumulativeSums->data.F32[binNum + 1];
-
-        if (!((y->data.F64[0] <= fitFloat) && (fitFloat <= y->data.F64[2]))) {
-            psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n");
-        }
-
+    float tmpFloat;
+
+    if ((binNum > 0) && (binNum < (yVec->n - 2))) {
+        x->data.F64[0] = (double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
+        x->data.F64[1] = (double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
+        x->data.F64[2] = (double) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2]));
+        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];
+
+        // Ensure that yVal is within the range of the bins we are using.
+        if (!((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2]))) {
+            printf("((%f), %f, %f)\n", yVal, y->data.F64[0], y->data.F64[2]);
+            psError(__func__, "yVal not within y-range\n");
+        }
         yErr->data.F64[0] = 1.0;
         yErr->data.F64[1] = 1.0;
         yErr->data.F64[2] = 1.0;
 
+        // Determine the coefficients of the polynomial.
         myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
-        return (p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat));
+        // Call p_ps1DPolyMedian(), which does a binary search on the
+        // polynomial, looking for the value x such that f(x) = yVal
+        tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal);
     } else {
-        return (0.5 * (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum]));
+        if (binNum == 0) {
+            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
+                              xVec->data.F32[binNum + 1]);
+        } else if (binNum == (xVec->n - 1)) {
+            // XXX: Is this right?
+            tmpFloat = xVec->data.F32[binNum];
+        } else if (binNum == (xVec->n - 2)) {
+            // XXX: Is this right?
+            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
+                              xVec->data.F32[binNum + 1]);
+        }
     }
 
@@ -1083 +1097 @@
     psFree(yErr);
     psFree(myPoly);
-    return (0.0);
-}
+    return(tmpFloat);
+}
+
+#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
+(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
 
 /******************************************************************************
@@ -1111 +1128 @@
 *****************************************************************************/
 void p_psVectorRobustStats(const psVector* restrict myVector,
-                           const psVector* restrict maskVector, unsigned int maskVal, psStats* stats)
+                           const psVector* restrict maskVector,
+                           unsigned int maskVal,
+                           psStats* stats)
 {
     psHistogram* robustHistogram = NULL;
@@ -1118 +1137 @@
     int LQBinNum = -1;          // Bin num for lower quartile
     int UQBinNum = -1;          // Bin num for upper quartile
+    int medianBinNum = -1;
     int i = 0;                  // Loop index variable
-    int maxBinNum = 0;
-    float maxBinCount = 0.0;
+    int modeBinNum = 0;
+    float modeBinCount = 0.0;
     float dL = 0.0;
     int numBins = 0;
-    psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN);
-
-    // psImage* domain;
-    // psVector* errors;
-    // psVector* data;
-    // psVector* initialGuess;
-    // psVector* theParams;
-    // float chiSq=0.0;
-    // float max = -HUGE;
-    // float max_pos;
     float myMean = 0.0;
     float myStdev = 0.0;
@@ -1139 +1149 @@
     float sumSquares = 0.0;
     float sumDiffs = 0.0;
-    psVector* x = psVectorAlloc(3, PS_TYPE_F64);
-    psVector* y = psVectorAlloc(3, PS_TYPE_F64);
-    psVector* yErr = psVectorAlloc(3, PS_TYPE_F64);
-    psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
-    psVector* cumulativeRobustSumsFullRange = NULL;
-    psVector* cumulativeRobustSumsDlRange = NULL;
+    psVector* cumulativeRobustSums = NULL;
     float sumRobust = 0.0;
     float sumN50 = 0.0;
     float sumNfit = 0.0;
-    float cumulativeMedian = 0.0;
+    psScalar tmpScalar;
+    tmpScalar.type.type = PS_TYPE_F32;
+    psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN);
 
     // Compute the initial bin size of the robust histogram.  This is done
@@ -1188 +1195 @@
         p_psVectorMax(myVector, maskVector, maskVal, stats);
     }
+
     // Create the histogram structure.  NOTE: we can not specify the bin size
     // precisely since the argument to psHistogramAlloc() is the number of
@@ -1198 +1206 @@
     psVectorHistogram(robustHistogram, myVector, maskVector, maskVal);
 
-    // Smooth the histogram.
-    // XXX: is that the right stdev?
-    robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram,
+    // Smooth the histogram, Gaussian-style.
+    robustHistogramVector = p_psVectorSmoothHistGaussian(robustHistogram,
                             tmpStats->clippedStdev / 4.0f);
 
@@ -1208 +1215 @@
 
     /**************************************************************************
-    Determine the lower/upper quartiles.
+    Determine the median/lower/upper quartile bin numbers.
+
+    We define a vector called "cumulativeRobustSums" where the value at
+    index position i is equal to the sum of bins 0:i.  This will be used in
+    determining the median and lower/upper quartiles.
     **************************************************************************/
-    // We define a vector called "cumulativeRobustSums..." where the value at
-    // index position i is equal to the sum of bins 0:i.  This will be used
-    // now and later in determining the lower/upper quartiles.
-    cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
-    cumulativeRobustSumsFullRange->data.F32[0] = robustHistogramVector->data.F32[0];
+    cumulativeRobustSums = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
+    cumulativeRobustSums->data.F32[0] = robustHistogramVector->data.F32[0];
     for (i = 1; i < robustHistogramVector->n; i++) {
-        cumulativeRobustSumsFullRange->data.F32[i] =
-            cumulativeRobustSumsFullRange->data.F32[i - 1] + robustHistogramVector->data.F32[i];
-    }
-    sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1];
-
-    // Determine the bin number containing the lower quartile point.
-    LQBinNum = -1;
-    for (i = 0; i < cumulativeRobustSumsFullRange->n; i++) {
-        if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust / 4.0)) {
-            LQBinNum = i;
-            break;
-        }
-    }
-
-    // Determine the bin number containing the upper quartile point.
-    UQBinNum = -1;
-    for (i = cumulativeRobustSumsFullRange->n - 1; i >= 0; i--) {
-        if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0 * sumRobust / 4.0)) {
-            UQBinNum = i;
-            break;
-        }
-    }
-
-    if ((LQBinNum == -1) || (UQBinNum == -1)) {
-        psAbort(__func__, "Could not determine the robust lower/upper quartiles.");
-    }
-
+        cumulativeRobustSums->data.F32[i] = cumulativeRobustSums->data.F32[i - 1] +
+                                            robustHistogramVector->data.F32[i];
+    }
+    sumRobust = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1];
+
+    tmpScalar.data.F32 = sumRobust / 4.0;
+    LQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar);
+    tmpScalar.data.F32 = 3.0 * sumRobust / 4.0;
+    UQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar);
+    tmpScalar.data.F32 = sumRobust / 2.0;
+    medianBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar);
+
+    if ((LQBinNum < 0) || (UQBinNum < 0)) {
+        psError(__func__, "Could not determine the robust lower/upper quartile bin numbers.");
+    }
+    if (medianBinNum < 0) {
+        psError(__func__, "Could not determine the robust lower/upper quartile bin numbers.");
+    }
     /**************************************************************************
     Determine the mode in the range LQ:UQ.
     **************************************************************************/
     // Determine the bin with the peak value in the range LQ to UQ.
-    maxBinNum = LQBinNum;
-    maxBinCount = robustHistogramVector->data.F32[LQBinNum];
-    sumN50 = (float)robustHistogram->nums->data.S32[LQBinNum];
+    modeBinNum = LQBinNum;
+    modeBinCount = robustHistogramVector->data.F32[LQBinNum];
+    sumN50 = (float)robustHistogram->nums->data.U32[LQBinNum];
     for (i = LQBinNum + 1; i <= UQBinNum; i++) {
-        if (robustHistogramVector->data.F32[i] > maxBinCount) {
-            maxBinNum = i;
-            maxBinCount = robustHistogramVector->data.F32[i];
-        }
-        sumN50 += (float)robustHistogram->nums->data.S32[i];
-    }
-
-    // XXX: is dL defined as the value at the LQ/UQ, or the bin number?
+        if (robustHistogramVector->data.F32[i] > modeBinCount) {
+            modeBinNum = i;
+            modeBinCount = robustHistogramVector->data.F32[i];
+        }
+        sumN50 += (float)robustHistogram->nums->data.U32[i];
+    }
+
     dL = (UQBinNum - LQBinNum) / 4;
-
-    printf("(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum);
-
     /**************************************************************************
     Determine the mean/stdev for the bins in the range mode-dL to mode+dL
     **************************************************************************/
-    cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
-    for (i = 0; i < robustHistogramVector->n; i++) {
-        cumulativeRobustSumsDlRange->data.F32[i] = 0.0;
-    }
-    sumNfit = 0.0;
-    cumulativeMedian = 0.0;
-    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
-        if ((0 <= i) && (i < robustHistogramVector->n)) {
-            cumulativeRobustSumsDlRange->data.F32[i] =
-                cumulativeRobustSumsDlRange->data.F32[i - 1] + robustHistogramVector->data.F32[i];
-            cumulativeMedian += robustHistogramVector->data.F32[i];
-            sumNfit += (float)robustHistogram->nums->data.S32[i];
-        }
-    }
-
     // Calculate the mean of the smoothed robust histogram in the range
     // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
     // that bin (this is a non-exact approximation).
+    sumNfit = 0.0;
     myMean = 0.0;
-    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
+    for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
         if ((0 <= i) && (i < robustHistogramVector->n)) {
-            myMean += (robustHistogramVector->data.F32[i]) * 0.5 *
-                      (robustHistogram->bounds->data.F32[i + 1] + robustHistogram->bounds->data.F32[i]);
+            myMean += (robustHistogramVector->data.F32[i]) * PS_BIN_MIDPOINT(robustHistogram, i);
             countFloat += robustHistogramVector->data.F32[i];
         }
+
+        sumNfit += (float)robustHistogram->nums->data.U32[i];
     }
     myMean /= countFloat;
@@ -1297 +1279 @@
     // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
     // that bin.
-    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
+    for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
         if ((0 <= i) && (i < robustHistogramVector->n)) {
-            diff = (0.5 * (robustHistogram->bounds->data.F32[i + 1] +
-                           robustHistogram->bounds->data.F32[i])) - myMean;
+            diff = PS_BIN_MIDPOINT(robustHistogram, i) - myMean;
             sumSquares += diff * diff * robustHistogramVector->data.F32[i];
             sumDiffs += diff * robustHistogramVector->data.F32[i];
@@ -1315 +1296 @@
 
     if (stats->options & PS_STAT_ROBUST_MODE) {
-        stats->robustMode = 0.5 *
-                            (robustHistogram->bounds->data.F32[maxBinNum] + robustHistogram->bounds->data.F32[maxBinNum + 1]);
+        stats->robustMode = PS_BIN_MIDPOINT(robustHistogram, modeBinNum);
     }
 
@@ -1322 +1302 @@
         stats->robustStdev = myStdev;
     }
-    // To determine the median (and later, the lower/upper quartile), we fit
-    // a quadratic to the three bins surrounding the bin containing the median.
-    // The quadratic y=f(x) with x being the midpoint of each bin, and y being
-    // the cumulative number of data points in all bins up to, and including,
-    // this bin.  We then solve the quadratic for
-
-    // XXX: Create a "p_psQuadraticFitAndSolveForX() function.
-
+
+    // To determine the median, we fit a quadratic y=f(x) to the three bins
+    // surrounding the bin containing the median (x is the midpoint of each
+    // bin and y is the value of each bin).  Then we figure out what value
+    // of x corresponds to f(x) being the median (half of all points).
     if (stats->options & PS_STAT_ROBUST_MEDIAN) {
-        if ((maxBinNum > 0) && (maxBinNum < (robustHistogram->nums->n - 1))) {
-            x->data.F64[0] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[maxBinNum - 1] +
-                              robustHistogram->bounds->data.F32[maxBinNum]);
-            x->data.F64[1] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[maxBinNum] +
-                              robustHistogram->bounds->data.F32[maxBinNum + 1]);
-            x->data.F64[2] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[maxBinNum + 1] +
-                              robustHistogram->bounds->data.F32[maxBinNum + 2]);
-
-            y->data.F64[0] = cumulativeRobustSumsDlRange->data.F32[maxBinNum - 1];
-            y->data.F64[1] = cumulativeRobustSumsDlRange->data.F32[maxBinNum];
-            y->data.F64[2] = cumulativeRobustSumsDlRange->data.F32[maxBinNum + 1];
-
-            // Ensure that cumulativeMedian/2 is actually within the range of the bins
-            // we are using.
-            cumulativeMedian *= 0.5;
-            if (!((y->data.F64[0] <= cumulativeMedian) && (cumulativeMedian <= y->data.F64[2]))) {
-                printf("((%f), %f, %f)\n", cumulativeMedian, y->data.F64[0], y->data.F64[2]);
-                psAbort(__func__, "p_psVectorRobustStats(1): midpoint not within y-range\n");
-            }
-            // XXX: yErr is not currently used by psVectorFitPolynomial1D().  We
-            // may have to set this meaningfully later.
-            yErr->data.F64[0] = 1.0;
-            yErr->data.F64[1] = 1.0;
-            yErr->data.F64[2] = 1.0;
-
-            // Determine the coefficients of the polynomial.
-            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
-            // Call p_ps1DPolyMedian(), which does a binary search on the
-            // polynomial, looking for the value x such that
-            // f(x) = cumulativeMedian.
-            stats->robustMedian = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], cumulativeMedian);
-        } else {
-            // If the mode is the first/last histogram bin, then simply use
-            // the midpoint of that bin.
-            stats->robustMedian = 0.5 * (robustHistogram->bounds->data.F32[maxBinNum + 1] +
-                                         robustHistogram->bounds->data.F32[maxBinNum]);
-        }
-    }
-    // The lower/upper quartile calculations are very similar to the median
-    // calculations.  We fit a quadratic to the array containing the
-    // cumulative data points in each bin, then search for x such that
-    // f(x) equals the lower/upper quartile exactly.
-    //
+        // Take a psVector.  Fit a polynomial y = f(x) to the 3 data elements
+        // surrounding medianBinNum, then find the x-value corresponding y = sumRobust/2.0.
+        stats->robustMedian = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds,
+                              cumulativeRobustSums,
+                              medianBinNum,
+                              sumRobust/2.0);
+    }
+
+    // To determine the quartiles, we fit a quadratic y=f(x) to the three bins
+    // surrounding the bin containing LQ/UQ (x is the midpoint of each
+    // bin and y is the value of each bin).  Then we figure out what value
+    // of x corresponds to f(x) being the LQ/UQ.
     if (stats->options & PS_STAT_ROBUST_QUARTILE) {
-        countFloat = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1];
-
-        if ((LQBinNum > 0) && (LQBinNum < (robustHistogram->nums->n - 1))) {
-            x->data.F64[0] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[LQBinNum - 1] +
-                              robustHistogram->bounds->data.F32[LQBinNum]);
-            x->data.F64[1] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[LQBinNum] +
-                              robustHistogram->bounds->data.F32[LQBinNum + 1]);
-            x->data.F64[2] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[LQBinNum + 1] +
-                              robustHistogram->bounds->data.F32[LQBinNum + 2]);
-
-            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[LQBinNum - 1];
-            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[LQBinNum];
-            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[LQBinNum + 1];
-
-            if (!((y->data.F64[0] <= (countFloat / 4.0)) && ((countFloat / 4.0) <= y->data.F64[2]))) {
-                psAbort(__func__, "p_psVectorRobustStats(2): midpoint not within y-range\n");
-            }
-
-            yErr->data.F64[0] = 1.0;
-            yErr->data.F64[1] = 1.0;
-            yErr->data.F64[2] = 1.0;
-
-            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
-            stats->robustLQ = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], countFloat / 4.0);
-
-        } else {
-            // If the LQ is the first/last histogram bin, then simply use
-            // the midpoint of that bin.
-            stats->robustLQ = 0.5 * (robustHistogram->bounds->data.F32[LQBinNum + 1] +
-                                     robustHistogram->bounds->data.F32[LQBinNum]);
-        }
-
-        if ((UQBinNum > 0) && (UQBinNum < (robustHistogram->nums->n - 1))) {
-            x->data.F64[0] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[UQBinNum - 1] +
-                              robustHistogram->bounds->data.F32[UQBinNum]);
-            x->data.F64[1] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[UQBinNum] +
-                              robustHistogram->bounds->data.F32[UQBinNum + 1]);
-            x->data.F64[2] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[UQBinNum + 1] +
-                              robustHistogram->bounds->data.F32[UQBinNum + 2]);
-
-            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[UQBinNum - 1];
-            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[UQBinNum];
-            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[UQBinNum + 1];
-
-            if (!((y->data.F64[0] <= (3.0 * countFloat / 4.0)) &&
-                    ((3.0 * countFloat / 4.0) <= y->data.F64[2]))) {
-                psAbort(__func__, "p_psVectorRobustStats(3): midpoint not within y-range\n");
-            }
-
-            yErr->data.F64[0] = 1.0;
-            yErr->data.F64[1] = 1.0;
-            yErr->data.F64[2] = 1.0;
-
-            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
-            stats->robustUQ = p_ps1DPolyMedian(myPoly,
-                                               x->data.F64[0], x->data.F64[2], 3.0 * countFloat / 4.0);
-        } else {
-            // If the UQ is the first/last histogram bin, then simply use
-            // the midpoint of that bin.
-            stats->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum + 1] +
-                                     robustHistogram->bounds->data.F32[UQBinNum]);
-        }
-    }
+        countFloat = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1];
+
+        stats->robustLQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds,
+                          cumulativeRobustSums,
+                          LQBinNum,
+                          countFloat/4.0);
+        stats->robustUQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds,
+                          cumulativeRobustSums,
+                          UQBinNum,
+                          3.0 * countFloat/4.0);
+    }
+    // XXX: I think sumNfit == sumN50 here.
     stats->robustNfit = sumNfit;
     stats->robustN50 = sumN50;
 
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
     psFree(tmpStats);
     psFree(robustHistogram);
-    psFree(myPoly);
-    psFree(cumulativeRobustSumsFullRange);
-    psFree(cumulativeRobustSumsDlRange);
+    psFree(robustHistogramVector);
+    psFree(cumulativeRobustSums);
 }
 
@@ -1468 +1350 @@
     // I am commenting this code out, and replacing it with code which
     // calculates the mean directly on the robustHistogram.
-    //
-    // XXX: Replace this with the LaGrange interpolation stuff.
  
     domain = psImageAlloc(1, robustHistogramVector->n, PS_TYPE_F32);
@@ -1479 +1359 @@
     max_pos = -1;
     for (i=0;i<robustHistogramVector->n;i++) {
-        domain->data.F32[i][0] = 0.5 * (robustHistogram->bounds->data.F32[i+1] +
-                                        robustHistogram->bounds->data.F32[i]);
+        domain->data.F32[i][0] =  PS_BIN_MIDPOINT(robustHistogram, i);bounds->data.F32[modeBinNum + 1]);
         data->data.F32[i] = (float) robustHistogramVector->data.F32[i];
         errors->data.F32[i] = 1.0;
@@ -1695 +1574 @@
     psScalar tmpScalar;
     tmpScalar.type.type = PS_TYPE_F32;
+    psVector* inF32;
+    int mustFreeTmp = 1;
 
     // NOTE: Verify that this is the correct action.
@@ -1706 +1587 @@
 
     if (out->nums->type.type != PS_TYPE_U32) {
-        psAbort(__func__, "Only data type PS_TYPE_U32 for output.nums member.");
+        psAbort(__func__, "Only data type PS_TYPE_U32 for out->nums member.");
     }
     // NOTE: Verify that this is the correct action.
@@ -1713 +1594 @@
     }
 
-    if (in->type.type != PS_TYPE_F32) {
-        psAbort(__func__, "Only data type PS_TYPE_F32 is currently supported (0x%x).", in->type.type);
-    }
-
     if (mask != NULL) {
         if (in->n != mask->n) {
-            psAbort(__func__, "Vector data and vector mask are of different sizes.");
+            psError(__func__, "Vector data and vector mask are of different sizes.");
         }
         if (mask->type.type != PS_TYPE_U8) {
-            psAbort(__func__, "Vector mask must be type PS_TYPE_U8");
-        }
+            psError(__func__, "Vector mask must be type PS_TYPE_U8");
+        }
+    }
+
+    inF32 = p_psConvertToF32((psVector *) in);
+    if (inF32 == NULL) {
+        inF32 = (psVector *) in;
+        mustFreeTmp = 0;
     }
 
@@ -1730 +1613 @@
 
     numBins = out->nums->n;
-    for (i = 0; i < in->n; i++) {
+    for (i = 0; i < inF32->n; i++) {
         // Check if this pixel is masked, and if so, skip it.
         if ((mask == NULL) || ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) {
             // Check if this pixel is below the minimum value, and if so
             // count it, then skip it.
-            if (in->data.F32[i] < out->bounds->data.F32[0]) {
+            if (inF32->data.F32[i] < out->bounds->data.F32[0]) {
                 out->minNum++;
                 // Check if this pixel is above the maximum value, and if so
                 // count it, then skip it.
-            } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
+            } else if (inF32->data.F32[i] > out->bounds->data.F32[numBins]) {
                 out->maxNum++;
             } else {
@@ -1746 +1629 @@
                 if (out->uniform == true) {
                     binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
-                    binNum = (int)((in->data.F32[i] - out->bounds->data.F32[0]) / binSize);
+                    binNum = (int)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize);
 
                     // NOTE: This next if-statement really shouldn't be necessary.
@@ -1755 +1638 @@
                     }
 
-                    (out->nums->data.S32[binNum])++;
+                    (out->nums->data.U32[binNum])++;
 
                     // If this is a non-uniform histogram, determining the correct
                     // bin number requires a bit more work.
                 } else {
-                    tmpScalar.data.F32 = in->data.F32[i];
-                    tmp = p_psVectorBinDisect(out->bounds,
-                                              &tmpScalar);
-                    if (tmp == -1) {
-                        //                      XXX: Generate warning message
+                    tmpScalar.data.F32 = inF32->data.F32[i];
+                    tmp = p_psVectorBinDisect(out->bounds, &tmpScalar);
+                    if (tmp < 0) {
+                        // XXX: Generate warning message
                     } else {
-                        (out->nums->data.S32[tmp])++;
+                        (out->nums->data.U32[tmp])++;
                     }
                 }
@@ -1772 +1654 @@
         }
     }
+
+    if (mustFreeTmp == 1) {
+        psFree(inF32);
+    }
     return (out);
 }
 
 /******************************************************************************
-p_psConvertToF32(myVector, maskVector, maskVal, stats): this is the cheap
-way to support a variety of vector data types: we simply convert the input
-vector to F32 at the beginning, and write all of our functions in F32.  If
-the vast majority of all vector stat operations are F32 (or any other single
-type), then this is probably the best way to go.  Otherwise, when the
-algorithms stablize, we will then macro everything and put type support in
-the various stat functions.
+p_psConvertToF32(in): this is the cheap way to support a variety of vector
+data types: we simply convert the input vector to F32 at the beginning, and
+write all of our functions in F32.  If the vast majority of all vector stat
+operations are F32 (or any other single type), then this is probably the
+best way to go.  Otherwise, when the algorithms stablize, we will then macro
+everything and put type support in the various stat functions.
  
 XXX: Should the default data type be F64?  Since we are buying Athlons...
  *****************************************************************************/
-psVector* p_psConvertToF32(psStats* stats, psVector* in, psVector* mask, unsigned int maskVal)
+psVector* p_psConvertToF32(psVector* in)
 {
     int i = 0;
@@ -1800 +1685 @@
         for (i = 0; i < in->n; i++) {
             tmp->data.F32[i] = (float)in->data.U16[i];
+        }
+    } else if (in->type.type == PS_TYPE_U8) {
+        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
+        for (i = 0; i < in->n; i++) {
+            tmp->data.F32[i] = (float)in->data.U8[i];
         }
     } else if (in->type.type == PS_TYPE_F64) {
@@ -1839 +1729 @@
     int mustFreeTmp = 1;
 
-    // NOTE: Verify that this is the correct action.
     if (in == NULL) {
-        return (stats);
-    }
-    if (stats == NULL) {
-        return (NULL);
-    }
-
-    inF32 = p_psConvertToF32(stats, in, mask, maskVal);
+        psError(__func__, "in is NULL\n");
+        return(stats);
+    }
+    PS_CHECK_NULL_PTR_RETURN_NULL(stats);
+
+    inF32 = p_psConvertToF32((psVector *) in);
     if (inF32 == NULL) {
         inF32 = in;
         mustFreeTmp = 0;
     }
-    // XXX: Should we abort if (stats->min == stats->max)?
     if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) {
-        psAbort(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max);
+        psError(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max);
+        return(stats);
     }