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

Timestamp:

Dec 8, 2006, 1:38:54 AM (19 years ago)

Author:

magnier

Message:

• moved psHistogram functions to psHistogram.c,h (pslib_strict.h, MAkefile.am, psMemory.c)
• changed psVectorStats return value to bool (psImageGeomManip.c, psImagePixelExtract.c, psImageStats.c, psMinPoly.c,
• added FITTED_MEAN,STDEV_V2 (quadratic fit to log-peak)
• added stats->results
• set stats->results when stats are measured
• return stats calculated even if not requested
• moved vectorFittedStats call out of vectorRobustStats to psVectorStats
• internal calls to private stats function XXX now tests stats->results for the corresponding flag, NOT if value of stats->XXX == NAN.
• dropped unused function NonEmpty (commented out)
• rolled up vectorSampleMedian with if tests in the loop
• rolled up vectorSampleStdev with if tests in the loop
• call psError on failure for any internal stats call (the user should call psErrorClear if they ignore errors from psVectorStats)

File:

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

trunk/psLib/src/math/psStats.c

@@ -7 +7 @@
  *  on those data structures.
  *
- *  @author GLG, MHPCC
- *
- *  XXX: The following stats members are never used, or set in this code.
- *      stats->clippedNvalues
+ *  @author GLG (MHPCC), EAM (IfA)
  *
  * XXX: Must do
@@ -16 +13 @@
  * use ->min and ->max (PS_STAT_USE_RANGE)
  *
- *  @version $Revision: 1.193 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-12-05 20:05:57 $
+ *  @version $Revision: 1.194 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-12-08 11:38:54 $
  *
- *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *  Copyright 2006 IfA, University of Hawaii
  */
 
@@ -35 +32 @@
 #include "psMemory.h"
 #include "psAbort.h"
-#include "psImage.h"
+//#include "psImage.h"
 #include "psVector.h"
 #include "psTrace.h"
 #include "psLogMsg.h"
 #include "psError.h"
+#include "psHistogram.h"
 #include "psStats.h"
 #include "psMinimizeLMM.h"
@@ -64 +62 @@
 #define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
 (0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
+
 /*****************************************************************************/
 /* TYPE DEFINITIONS                                                          */
@@ -87 +86 @@
 NOTE: it is assumed that any call to these statistical functions will have
 been preceded by a call to the psVectorStats() function.  Various sanity tests
-will only be performed in psVectorStats().  Should we perform the sanity
-checks in each routine anyway?
+will only be performed in psVectorStats().
  
-XXX: For many of these private stats routines, what should be done if there
-are no acceptable elements in the input vector (if no elements lie within
-range, or there are no unmasked elements, or the input vector is NULL)?
-Currently we set the value to NAN.
- 
-XXX: Optimization: many routines have an "empty" boolean variable which keeps
-track of whether or not the vector has any valid elements.  This code can
-possibly be optimized away.
- *****************************************************************************/
-/******************************************************************************
-vectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
-mean of the input vector.  If there was a problem with the mean calculation,
-this routine sets stats->sampleMean to NAN.
- *****************************************************************************/
-# if (0)
-    static bool vectorSampleMean(const psVector* myVector,
-                                 const psVector* errors,
-                                 const psVector* maskVector,
-                                 psMaskType maskVal,
-                                 psStats* stats)
-{
-    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
-
-    psF32 mean = 0.0;                   // The mean
-    long count = 0;                     // Number of points contributing to this mean
-
-    psF32 *data = myVector->data.F32;   // Dereference
-    int numData = myVector->n;          // Number of data points
-
-    // If PS_STAT_USE_RANGE is requested, then we enter a slightly different loop.
-    if (!errors) {
-        if (stats->options & PS_STAT_USE_RANGE) {
-            if (maskVector) {
-                // No errors, check range, check mask
-                psU8 *maskData = maskVector->data.U8;
-                for (long i = 0; i < numData; i++) {
-                    // Check if the data is with the specified range
-                    if (!(maskVal & maskData[i]) && stats->min <= data[i] && data[i] <= stats->max) {
-                        mean += data[i];
-                        count++;
-                    }
-                }
-            } else {
-                // No errors, check range, no mask
-                for (long i = 0; i < numData; i++) {
-                    if (stats->min <= data[i] && data[i] <= stats->max) {
-                        mean += data[i];
-                        count++;
-                    }
-                }
-            }
-        } else {
-            if (maskVector) {
-                // No errors, no range check, check mask
-                psU8 *maskData = maskVector->data.U8;
-                for (long i = 0; i < numData; i++) {
-                    if (!(maskVal & maskData[i])) {
-                        mean += data[i];
-                        count++;
-                    }
-                }
-
-            } else {
-                // No errors, no range check, no mask
-                for (long i = 0; i < numData; i++) {
-                    mean += data[i];
-                }
-                count = numData;
-            }
-        }
-
-        if (count > 0) {
-            mean /= (psF32)count;
-        } else {
-            mean = NAN;
-        }
-
-    } else {
-        psF32 sumWeights = 0.0;         // The sum of the weights
-        psF32 *errorsData = errors->data.F32;
-        if (stats->options & PS_STAT_USE_RANGE) {
-            if (maskVector) {
-                psU8 *maskData = maskVector->data.U8;
-                for (long i = 0; i < numData; i++) {
-                    // Check if the data is with the specified range
-                    if (!(maskVal & maskData[i]) && stats->min <= data[i] && data[i] <= stats->max) {
-                        float weight = 1.0 / PS_SQR(errorsData[i]);
-                        mean += data[i] * weight;
-                        sumWeights += weight;
-                        count++;
-                    }
-                }
-            } else {
-                for (long i = 0; i < myVector->n; i++) {
-                    if (stats->min <= data[i] && data[i] <= stats->max) {
-                        float weight = 1.0 / PS_SQR(errorsData[i]);
-                        mean += data[i] * weight;
-                        sumWeights += weight;
-                        count++;
-                    }
-                }
-            }
-        } else {
-            if (maskVector) {
-                psU8 *maskData = maskVector->data.U8;
-                for (long i = 0; i < myVector->n; i++) {
-                    if (!(maskVal & maskData[i])) {
-                        float weight = 1.0 / PS_SQR(errorsData[i]);
-                        mean += data[i] * weight;
-                        sumWeights += weight;
-                        count++;
-                    }
-                }
-            } else {
-                for (long i = 0; i < myVector->n; i++) {
-                    float weight = 1.0 / PS_SQR(errorsData[i]);
-                    mean += data[i] * weight;
-                    sumWeights += weight;
-                }
-                count = myVector->n;
-            }
-        }
-
-        if (count > 0) {
-            mean /= sumWeights;
-        } else {
-            mean = NAN;
-        }
-    }
-
-    stats->sampleMean = mean;
-    if (isnan(mean)) {
-        psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
-        return false;
-    }
-
-    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
-    return true;
-}
-#endif
+For many of these private stats routines, it is possible that there are no acceptable elements
+in the input vector (if no elements lie within range, or there are no unmasked elements, or the
+input vector is NULL).  In such cases, we set the desired value to NAN and return an error.
+The calling function may clear the error.
+*****************************************************************************/
+
+// static prototypes:
+static psF32 minimizeLMChi2Gauss1D(psVector *deriv, const psVector *params, const psVector *coords);
+static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
 
 /******************************************************************************
@@ -236 +103 @@
 this routine sets stats->sampleMean to NAN.
  
-XXX : using the method below, with a single loop for various options
-      costs only a small amount and is much easier to debug. running 10000 tests
-      of 1000 point vectors for the two methods gives:
-                     separate    single
-(mask: 0, range: 0): 0.067 sec   0.073 sec
-(mask: 1, range: 0): 0.098 sec  0.102 sec
-(mask: 0, range: 0): 0.136 sec  0.162 sec
-(mask: 1, range: 0): 0.177 sec  0.181 sec
- *****************************************************************************/
+using the method below, with a single loop for various options costs only a small amount and is
+much easier to debug. running 10000 tests of 1000 point vectors for the two methods gives:
+ 
+                     separate    single          w/errors
+(mask: 0, range: 0): 0.067 sec  0.073 sec (10%)  0.072 sec
+(mask: 1, range: 0): 0.098 sec  0.102 sec ( 4%)  0.119 sec
+(mask: 0, range: 0): 0.136 sec  0.162 sec (20%)  0.170 sec
+(mask: 1, range: 0): 0.177 sec  0.181 sec ( 2%)  0.198 sec
+ 
+(effect of errors not tested)
+ 
+To optmize this, use a macro and ifdef in or out the three states (errors, mask, range)
+*****************************************************************************/
 static bool vectorSampleMean(const psVector* myVector,
                              const psVector* errors,
@@ -253 +124 @@
     psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
 
+    long count = 0;                     // Number of points contributing to this mean
     psF32 mean = 0.0;                   // The mean
-    long count = 0;                     // Number of points contributing to this mean
+    psF32 weight;
 
     psF32 *data = myVector->data.F32;   // Dereference
@@ -262 +134 @@
     bool useRange = stats->options & PS_STAT_USE_RANGE;
 
-    // If errors exist, we enter a slightly different loop.
-    // XXX collapse to a single loop?
-    if (!errors) {
-        for (long i = 0; i < numData; i++) {
-            // Check if the data is with the specified range
-            if (useRange && (data[i] < stats->min))
-                continue;
-            if (useRange && (data[i] > stats->max))
-                continue;
-            if (maskData && (maskData[i] & maskVal))
-                continue;
-            mean += data[i];
-            count++;
-        }
-        mean = (count > 0) ? mean / (psF32) count : NAN;
-    } else {
-        psF32 sumWeights = 0.0;         // The sum of the weights
-        psF32 *errorsData = errors->data.F32;
-        for (long i = 0; i < numData; i++) {
-            // Check if the data is with the specified range
-            if (useRange && (data[i] < stats->min))
-                continue;
-            if (useRange && (data[i] > stats->max))
-                continue;
-            if (maskData && (maskData[i] & maskVal))
-                continue;
-
-            float weight = 1.0 / PS_SQR(errorsData[i]);
+    psF32 sumWeights = 0.0;  // The sum of the weights
+    psF32 *errorsData = (errors == NULL) ? NULL : errors->data.F32;
+
+    for (long i = 0; i < numData; i++) {
+        // Check if the data is with the specified range
+        if (useRange && (data[i] < stats->min))
+            continue;
+        if (useRange && (data[i] > stats->max))
+            continue;
+        if (maskData && (maskData[i] & maskVal))
+            continue;
+        if (errors) {
+            weight = (errorsData[i] == 0) ? 0.0 : PS_SQR(errorsData[i]);
             mean += data[i] * weight;
             sumWeights += weight;
-            count++;
-        }
+        } else {
+            mean += data[i];
+        }
+        count++;
+    }
+    if (errors) {
         mean = (count > 0) ? mean / sumWeights : NAN;
-    }
-
+    } else {
+        mean = (count > 0) ? mean / count : NAN;
+    }
     stats->sampleMean = mean;
     if (isnan(mean)) {
+        // XXX raise an error here?
         psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
         return false;
     }
 
+    stats->results |= PS_STAT_SAMPLE_MEAN;
     psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
     return true;
@@ -319 +183 @@
 (mask: 0, range: 0):  0.179 sec  0.208 sec
 (mask: 1, range: 0):  0.200 sec  0.244 sec
- *****************************************************************************/
+*****************************************************************************/
 static long vectorMinMax(const psVector* myVector,
                          const psVector* maskVector,
@@ -350 +214 @@
     }
 
+    // XXX save numValid in psStats?
     if (numValid == 0) {
         stats->max = NAN;
@@ -356 +221 @@
         stats->max = max;
         stats->min = min;
+        stats->results |= PS_STAT_MIN;
+        stats->results |= PS_STAT_MAX;
     }
     psTrace("psLib.math", 4, "---- %s(%d) end ----\n", __func__, numValid);
     return numValid;
 }
-
-#if 0
-/******************************************************************************
-vectorCheckNonEmpty(myVector, maskVector, maskVal, stats): This routine returns
-"true" if the inputPsVector has 1 or more valid elements (a valid element is an
-unmasked element within the specified min/max range).  Otherwise, return
-"false".
- *****************************************************************************/
-static bool vectorCheckNonEmpty(const psVector* myVector,
-                                const psVector* maskVector,
-                                psMaskType maskVal,
-                                psStats* stats)
-{
-    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_VECTOR_NON_NULL(myVector, false);
-    PS_ASSERT_PTR_NON_NULL(stats, false);
-
-    if (stats->options & PS_STAT_USE_RANGE) {
-        if (maskVector) {
-            for (long i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i]) &&
-                        (stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
-                    return true;
-                }
-            }
-        } else {
-            for (long i = 0; i < myVector->n; i++) {
-                if ((stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
-                    return true;
-                }
-            }
-        }
-    } else {
-        if (maskVector != NULL) {
-            for (long i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
-                    return true;
-                }
-            }
-        } else {
-            if (myVector->n > 0) {
-                psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
-                return true;
-            }
-        }
-    }
-    psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
-    return(false);
-}
-#endif
 
 /******************************************************************************
 vectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the median
 and quartiles of the input vector.  Returns true on success (including if there
-were no valid input vector elements).
- *****************************************************************************/
-static bool vectorSampleMedian(const psVector* myVector,
+were no valid input vector elements). Expects F32 vector for input.
+*****************************************************************************/
+static bool vectorSampleMedian(const psVector* inVector,
                                const psVector* maskVector,
                                psMaskType maskVal,
@@ -428 +240 @@
     psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
 
+    bool useRange = stats->options & PS_STAT_USE_RANGE;
+    psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8; // Dereference the vector
+    psF32 *input = inVector->data.F32; // Dereference the vector
+
     // Allocate temporary vectors for the data.
-    psVector* vector = psVectorAllocEmpty(myVector->n, PS_TYPE_F32); // Vector containing the unmasked values
+    psVector *outVector = psVectorAllocEmpty(inVector->n, PS_TYPE_F32); // Vector containing the unmasked values
+    psF32 *output = outVector->data.F32; // Dereference the vector
+
     long count = 0;                     // Number of valid entries
 
-    // Determine if we must only use data points within a min/max range.
-    if (stats->options & PS_STAT_USE_RANGE) {
-        // Store all non-masked data points within the min/max range
-        // into the temporary vectors.
-        if (maskVector != NULL) {
-            for (long i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i]) &&
-                        (stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    vector->data.F32[count] = myVector->data.F32[i];
-                    count++;
-                }
-            }
-        } else {
-            for (long i = 0; i < myVector->n; i++) {
-                if ((stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    vector->data.F32[count] = myVector->data.F32[i];
-                    count++;
-                }
-            }
-        }
-    } else {
-        // Store all non-masked data points into the temporary vectors.
-        if (maskVector) {
-            for (long i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    vector->data.F32[count] = myVector->data.F32[i];
-                    count++;
-                }
-            }
-        } else {
-            vector = psVectorCopy(vector, myVector, PS_TYPE_F32);
-            count = myVector->n;
-        }
-    }
-    vector->n = count;
+    // Store all non-masked data points within the min/max range
+    // into the temporary vectors.
+    for (long i = 0; i < inVector->n; i++) {
+        if (useRange && (input[i] < stats->min))
+            continue;
+        if (useRange && (input[i] > stats->max))
+            continue;
+        if (maskData && (maskData[i] & maskVal))
+            continue;
+
+        output[count] = input[i];
+        count++;
+    }
+    outVector->n = count;
+
     if (count == 0) {
         psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No valid data in input vector.\n");
-        psFree(vector);
+        psFree(outVector);
+        stats->sampleUQ = NAN;
+        stats->sampleLQ = NAN;
         stats->sampleMedian = NAN;
         return false;
@@ -477 +275 @@
 
     // Sort the temporary vector.
-    vector = psVectorSort(vector, vector); // Sort in-place (since it's a copy, it's OK)
-    if (!vector) {
-        psError(PS_ERR_UNEXPECTED_NULL,
-                false,
-                _("Failed to sort input data."));
-        psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
-        psFree(vector);
+    // XXX this is messed up: should psVectorSort free on error or not?
+    if (!psVectorSort(outVector, outVector)) { // Sort in-place (since it's a copy, it's OK)
+        psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to sort input data."));
+        stats->sampleUQ = NAN;
+        stats->sampleLQ = NAN;
+        stats->sampleMedian = NAN;
         return false;
     }
 
-    // Calculate the median exactly.  Use the average if the number of samples
-    // is even.
+    // Calculate the median.  Use the average if the number of samples if even.
     if (count % 2 == 0) {
-        stats->sampleMedian = 0.5 * (vector->data.F32[(count / 2) - 1] +
-                                     vector->data.F32[count / 2]);
+        stats->sampleMedian = 0.5 * (output[(count/2) - 1] + output[count/2]);
     } else {
-        stats->sampleMedian = vector->data.F32[count / 2];
+        stats->sampleMedian = output[count/2];
     }
 
     // Calculate the quartile points exactly.
-    stats->sampleUQ = vector->data.F32[3 * (count / 4)];
-    stats->sampleLQ = vector->data.F32[count / 4];
+    stats->sampleUQ = output[(int)(0.75*count)];
+    stats->sampleLQ = output[(int)(0.25*count)];
+
+    stats->results |= PS_STAT_SAMPLE_MEDIAN;
+    stats->results |= PS_STAT_SAMPLE_QUARTILE;
 
     // Free the temporary data structures.
-    psFree(vector);
+    psFree(outVector);
 
     // Return "true" on success.
     psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
     return true;
-}
-
-/******************************************************************************
-vectorSmoothHistGaussian(): This routine smoothes the data in the input
-robustHistogram with a Gaussian of width sigma.  It returns a psVector of the
-smoothed data.
- *****************************************************************************/
-psVector *vectorSmoothHistGaussian(psHistogram *histogram,
-                                   psF32 sigma)
-{
-    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
-    psTrace("psLib.math", 5, "(histogram->nums->n, sigma) is (%d, %.2f\n", (int) histogram->nums->n, sigma);
-    PS_ASSERT_PTR_NON_NULL(histogram, NULL);
-    PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL);
-    PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL);
-    if (psTraceGetLevel("psLib.math") >= 8) {
-        PS_VECTOR_PRINT_F32(histogram->nums);
-    }
-
-    long numBins = histogram->nums->n;  // Number of histogram bins
-    psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); // Smoothed version of histogram bins
-    const psVector *bounds = histogram->bounds; // The bounds for the histogram bins
-
-    if (!histogram->uniform) {
-        //
-        // We get here if the histogram is non-uniform.  This code is not tested.
-        // However, it is also not used anywhere, yet.
-        //
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSmoothHistGaussian() on non-uniform histograms has not been tested or used.\n");
-
-        for (long i = 0; i < numBins; i++) {
-            // Determine the midpoint of bin i.
-            float iMid = PS_BIN_MIDPOINT(histogram, i);
-
-            //
-            // We determine the bin numbers (jMin:jMax) corresponding to a
-            // range of data values surrounding iMid.  The range is of size:
-            // 2*PS_GAUSS_WIDTH*sigma
-            long jMin, jMax;
-            psF32 x = iMid - (PS_GAUSS_WIDTH * sigma);
-            for (jMin = i; jMin > 0 && bounds->data.F32[jMin] > x; jMin--)
-                ;
-            x = iMid + (PS_GAUSS_WIDTH * sigma);
-            for (jMax = i; jMax < bounds->n - 1 && bounds->data.F32[jMax + 1] > x; jMax++)
-                ;
-
-            //
-            // Loop from jMin to jMax, computing the gaussian of data i.
-            //
-            smooth->data.F32[i] = 0.0;
-            for (long j = jMin ; j <= jMax ; j++) {
-                float jMid = PS_BIN_MIDPOINT(histogram, j);
-                smooth->data.F32[i] +=
-                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
-            }
-        }
-    } else {
-        //
-        // We get here if the histogram is uniform.
-        //
-        for (long i = 0; i < numBins; i++) {
-            psF32 binSize = bounds->data.F32[1] - bounds->data.F32[0];
-            psS32 gaussWidth = ((PS_GAUSS_WIDTH * sigma) / binSize);
-
-            //
-            // XXX: The following is wrong.  However, in practice, the sigma was too
-            // large compared to the binSize.  This meant that the smoothing was done
-            // over 500 bins in the robust stats algorithm.  This mean that the smoothing
-            // took way too long.
-            //
-            #if 0
-
-            gaussWidth = 10;
-            #endif
-            //
-            // We determine the bin numbers (jMin:jMax) corresponding to a
-            // range of data values surrounding iMid.  The range is of size:
-            // 2*PS_GAUSS_WIDTH*sigma
-            //
-            psS32 jMin = PS_MAX(i - gaussWidth, 0);
-            psS32 jMax = PS_MIN(i + gaussWidth, bounds->n - 1);
-
-            //
-            // Loop from jMin to jMax, computing the gaussian of data i.
-            //
-            smooth->data.F32[i] = 0.0;
-            float iMid = PS_BIN_MIDPOINT(histogram, i);
-            for (long j = jMin ; j <= jMax ; j++) {
-                float jMid = PS_BIN_MIDPOINT(histogram, j);
-                smooth->data.F32[i] +=
-                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
-            }
-        }
-    }
-
-    if (psTraceGetLevel("psLib.math") >= 8) {
-        PS_VECTOR_PRINT_F32(smooth);
-    }
-    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
-    return(smooth);
 }
 
@@ -619 +317 @@
     NULL
  
- *****************************************************************************/
+using the method below, with a single loop for various options costs only a small amount and is
+much easier to debug. running 10000 tests of 1000 point vectors for the two methods gives:
+ 
+                     single
+(mask: 0, range: 0): 0.193 sec
+(mask: 1, range: 0): 0.257 sec
+(mask: 0, range: 0): 0.349 sec
+(mask: 1, range: 0): 0.401 sec
+ 
+*****************************************************************************/
 static bool vectorSampleStdev(const psVector* myVector,
                               const psVector* errors,
@@ -630 +337 @@
     // This procedure requires the mean.  If it has not been already
     // calculated, then call vectorSampleMean()
-    if (isnan(stats->sampleMean)) {
+    if (!(stats->results & PS_STAT_SAMPLE_MEAN)) {
         vectorSampleMean(myVector, errors, maskVector, maskVal, stats);
     }
+
     // If the mean is NAN, then generate a warning and set the stdev to NAN.
     if (isnan(stats->sampleMean)) {
@@ -641 +349 @@
     }
 
+    psF32 *data = myVector->data.F32;   // Dereference
+    psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;
+    bool useRange = stats->options & PS_STAT_USE_RANGE;
+    psF32 *errorsData = (errors == NULL) ? NULL : errors->data.F32;
+
     // Accumulate the sums
     double mean = stats->sampleMean;    // The mean
@@ -647 +360 @@
     double errorDivisor = 0.0;          // Division for errors
     long count = 0;                     // Number of data points being used
-    if (stats->options & PS_STAT_USE_RANGE) {
-        if (maskVector) {
-            for (long i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i]) &&
-                        (stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    double diff = myVector->data.F32[i] - mean;
-                    sumSquares += PS_SQR(diff);
-                    sumDiffs += diff;
-                    count++;
-                    if (errors) {
-                        errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
-                    }
-                }
-            }
-        } else {
-            for (long i = 0; i < myVector->n; i++) {
-                if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) {
-                    double diff = myVector->data.F32[i] - mean;
-                    sumSquares += PS_SQR(diff);
-                    sumDiffs += diff;
-                    count++;
-                    if (errors) {
-                        errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
-                    }
-                }
-            }
-            count = myVector->n;
-        }
-    } else {
-        if (maskVector) {
-            for (long i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    double diff = myVector->data.F32[i] - mean;
-                    sumSquares += PS_SQR(diff);
-                    sumDiffs += diff;
-                    count++;
-                    if (errors) {
-                        errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
-                    }
-                }
-            }
-        } else {
-            for (long i = 0; i < myVector->n; i++) {
-                double diff = myVector->data.F32[i] - mean;
-                sumSquares += PS_SQR(diff);
-                sumDiffs += diff;
-                count++;
-                if (errors) {
-                    errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
-                }
-            }
-            count = myVector->n;
+    for (long i = 0; i < myVector->n; i++) {
+        // Check if the data is with the specified range
+        if (useRange && (data[i] < stats->min))
+            continue;
+        if (useRange && (data[i] > stats->max))
+            continue;
+        if (maskData && (maskData[i] & maskVal))
+            continue;
+
+        double diff = data[i] - mean;
+        sumSquares += PS_SQR(diff);
+        sumDiffs += diff;
+        count++;
+        if (errors) {
+            errorDivisor += 1.0 / PS_SQR(errorsData[i]);
         }
     }
@@ -717 +392 @@
         stats->sampleStdev = (1.0 / sqrtf(errorDivisor));
     } else {
-        stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) /
-                                  (float)(count - 1));
-    }
+        stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) / (float)(count - 1));
+    }
+    stats->results |= PS_STAT_SAMPLE_STDEV;
+
     psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
 
@@ -736 +412 @@
 Returns
     true for success; false otherwise
- *****************************************************************************/
+*****************************************************************************/
 static bool vectorClippedStats(const psVector* myVector,
                                const psVector* errors,
@@ -757 +433 @@
                                  PS_CLIPPED_SIGMA_UB, -1);
 
-    // Allocate a psStats structure for calculating the mean, median, and
-    // stdev.
-    psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+    // unless we succeed, these will have NAN values
+    stats->clippedMean = NAN;
+    stats->clippedStdev = NAN;
+    stats->clippedNvalues = 0;
 
     // We copy the mask vector, to preserve the original
@@ -773 +450 @@
     }
 
-    // 1. Compute the sample median.
-    vectorSampleMedian(myVector, tmpMask, maskVal, statsTmp);
-    if (isnan(statsTmp->sampleMedian)) {
+    // 1. Compute the sample median, which we save for output
+    vectorSampleMedian(myVector, tmpMask, maskVal, stats);
+    if (isnan(stats->sampleMedian)) {
         psLogMsg(__func__, PS_LOG_WARN, "Call to vectorSampleMedian returned NAN\n");
-        stats->clippedMean = NAN;
-        stats->clippedStdev = NAN;
         psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
         psFree(tmpMask);
-        psFree(statsTmp);
         return false;
     }
-    psTrace("psLib.math", 6, "The initial sample median is %f\n", statsTmp->sampleMedian);
-
-    // 2. Compute the sample standard deviation.
-    vectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp);
-    if (isnan(statsTmp->sampleStdev)) {
+    psTrace("psLib.math", 6, "The initial sample median is %f\n", stats->sampleMedian);
+
+    // 2. Compute the sample standard deviation, which we save for output
+    vectorSampleStdev(myVector, errors, tmpMask, maskVal, stats);
+    if (isnan(stats->sampleStdev)) {
         psLogMsg(__func__, PS_LOG_WARN, "Call to vectorSampleStdev returned NAN\n");
-        stats->clippedMean = NAN;
-        stats->clippedStdev = NAN;
         psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
         psFree(tmpMask);
-        psFree(statsTmp);
         return false;
     }
-    psTrace("psLib.math", 6, "The initial sample stdev is %f\n", statsTmp->sampleStdev);
+    psTrace("psLib.math", 6, "The initial sample stdev is %f\n", stats->sampleStdev);
 
     // 3. Use the sample median as the first estimator of the mean X.
-    psF32 clippedMean = statsTmp->sampleMedian;
+    psF32 clippedMean = stats->sampleMedian;
 
     // 4. Use the sample stdev as the first estimator of the mean stdev.
-    psF32 clippedStdev = statsTmp->sampleStdev;
+    psF32 clippedStdev = stats->sampleStdev;
 
     // 5. Repeat N (stats->clipIter) times:
@@ -836 +507 @@
 
         // b) compute new mean and stdev
+        // Allocate a psStats structure for calculating the mean and stdev.
+        psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
         vectorSampleMean(myVector, errors, tmpMask, maskVal, statsTmp);
         vectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp);
@@ -853 +526 @@
             clippedStdev = statsTmp->sampleStdev;
         }
+        psFree(statsTmp);
     }
 
@@ -859 +533 @@
 
     // 7. The last calcuated value of x is the clipped mean.
-    if (stats->options & PS_STAT_CLIPPED_MEAN) {
-        stats->clippedMean = clippedMean;
-        psTrace("psLib.math", 6, "The final clipped mean is %f\n", clippedMean);
-    }
     // 8. The last calcuated value of stdev is the clipped stdev.
-    if (stats->options & PS_STAT_CLIPPED_STDEV) {
-        stats->clippedStdev = clippedStdev;
-        psTrace("psLib.math", 6, "The final clipped stdev is %f\n", clippedStdev);
-    }
+    // we always return both stats even if only one was requested
+    stats->clippedMean = clippedMean;
+    stats->clippedStdev = clippedStdev;
+
+    stats->results |= PS_STAT_CLIPPED_MEAN;
+    stats->results |= PS_STAT_CLIPPED_STDEV;
+
+    psTrace("psLib.math", 6, "The final clipped mean is %f\n", clippedMean);
+    psTrace("psLib.math", 6, "The final clipped stdev is %f\n", clippedStdev);
 
     psFree(tmpMask);
-    psFree(statsTmp);
     psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
     return true;
 }
-
-static psF32 QuadraticInverse(psF32 a,
-                              psF32 b,
-                              psF32 c,
-                              psF32 y,
-                              psF32 xLo,
-                              psF32 xHi
-                             )
-{
-    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;
-
-    #if 0
-
-    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);
-    #endif
-
-    if (xLo <= x1 && x1 <= xHi) {
-        return x1;
-    }
-    if (xLo <= x2 && x2 <= xHi) {
-        return x2;
-    }
-    return 0.5 * (xLo + xHi);
-}
-
-
-/******************************************************************************
-fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): A general
-routine which fits a quadratic to three points and returns the x-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 (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).
- 
-*****************************************************************************/
-static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec,
-        psVector *yVec,
-        psS32 binNum,
-        psF32 yVal
-                                              )
-{
-    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
-    psTrace("psLib.math", 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_VECTORS_SIZE_EQUAL(xVec, yVec, 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] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
-        x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
-        x->data.F64[2] = (psF64) (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];
-        psTrace("psLib.math", 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("psLib.math", 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
-        psTrace("psLib.math", 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]);
-        }
-
-        //
-        // 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);
-            psTrace("psLib.math", 5, "---- %s() end ----\n", __func__);
-            return NAN;
-        }
-
-        // Determine the coefficients of the polynomial.
-        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
-        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x);
-        if (myPoly == NULL) {
-            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);
-            psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__);
-            return NAN;
-        }
-        psTrace("psLib.math", 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
-        psTrace("psLib.math", 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
-        psTrace("psLib.math", 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
-        psTrace("psLib.math", 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("psLib.math", 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
-        tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal,
-                                    x->data.F64[0], x->data.F64[2]);
-        psFree(myPoly);
-
-        if (isnan(tmpFloat)) {
-            psError(PS_ERR_UNEXPECTED_NULL,
-                    false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
-            psFree(x);
-            psFree(y);
-            psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__);
-            return(NAN);
-        }
-    } 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.
-            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
-                              xVec->data.F32[binNum + 1]);
-        } else if (binNum == (xVec->n - 1)) {
-            // The special case where we have two points only at the end of
-            // the vectors x and y.
-            // 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]);
-        }
-    }
-
-    psTrace("psLib.math", 6, "FIT: return %f\n", tmpFloat);
-    psFree(x);
-    psFree(y);
-
-    psTrace("psLib.math", 5, "---- %s(%f) end ----\n", __func__, tmpFloat);
-    return tmpFloat;
-}
-
-/******************************************************************************
-NOTE: We assume unnormalized gaussians.
- *****************************************************************************/
-static psF32 minimizeLMChi2Gauss1D(psVector *deriv,
-                                   const psVector *params,
-                                   const psVector *coords
-                                  )
-{
-    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_SIZE(params, (long)2, NAN);
-    PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(coords, NAN);
-    PS_ASSERT_VECTOR_SIZE(coords, (long)1, NAN);
-    PS_ASSERT_VECTOR_TYPE(coords, PS_TYPE_F32, NAN);
-
-    psF32 x = coords->data.F32[0];
-    psF32 mean = params->data.F32[0];
-    psF32 var = params->data.F32[1];
-    psF32 dx = (x - mean);
-
-    psF32 gauss = exp (-0.5*PS_SQR(dx)/var);
-    if (deriv) {
-        PS_ASSERT_VECTOR_SIZE(deriv, (long)2, NAN);
-        PS_ASSERT_VECTOR_TYPE(deriv, PS_TYPE_F32, NAN);
-        psF32 tmp = dx * gauss;
-        deriv->data.F32[0] = tmp / var;
-        deriv->data.F32[1] = tmp * dx / (var * var);
-    }
-
-
-    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
-    return gauss;
-}
-
-/*
- * vectorFittedStats requires guess for fittedMean and fittedStdev
- * robustN50 should also be set
- */
-static bool vectorFittedStats (const psVector* myVector,
-                               const psVector* errors,
-                               psVector* mask,
-                               psMaskType maskVal,
-                               psStats* stats)
-{
-
-    psTrace("psLib.math", 6, "The guess mean  is %f.\n", stats->fittedMean);
-    psTrace("psLib.math", 6, "The guess stdev is %f.\n", stats->fittedStdev);
-
-    psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
-
-    psScalar tmpScalar;                 // Temporary scalar variable, for p_psVectorBinDisect
-    tmpScalar.type.type = PS_TYPE_F32;
-
-    psF32 binSize = 1;
-    if (stats->options & PS_STAT_USE_BINSIZE) {
-        // Set initial bin size to the specified value.
-        binSize = stats->binsize;
-        psTrace("psLib.math", 6, "Setting initial robust bin size to %.2f\n", binSize);
-    } else {
-        // construct a histogram with (sigma/10 < binsize < sigma)
-        // set roughly so that the lowest bins have about 2 cnts
-        // Nsmallest ~ N50 / (4*dN))
-        psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8));
-        binSize = stats->fittedStdev / dN;
-    }
-
-    // Determine the min/max of the vector (which prior outliers masked out)
-    int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values
-    float min = statsMinMax->min;
-    float max = statsMinMax->max;
-    if (numValid == 0 || isnan(min) || isnan(max)) {
-        psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
-        psFree(statsMinMax);
-        psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
-        return false;
-    }
-
-    // Calculate the number of bins.
-    long numBins = (max - min) / binSize;
-    psTrace("psLib.math", 6, "The new min/max values are (%f, %f).\n", min, max);
-    psTrace("psLib.math", 6, "The new bin size is %f.\n", binSize);
-    psTrace("psLib.math", 6, "The numBins is %ld\n", numBins);
-
-    psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
-    histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal);
-    if (psTraceGetLevel("psLib.math") >= 8) {
-        PS_VECTOR_PRINT_F32(histogram->nums);
-    }
-
-    long binMin = 0;                // Low bin to check
-    long binMax = 0;                // High bin to check
-
-    // Fit a Gaussian to the bins in the requested range about the guess mean
-    // min,max overrides clipSigma
-    psF32 maxFitSigma = 2.0;
-    if (isfinite(stats->clipSigma)) {
-        maxFitSigma = fabs(stats->clipSigma);
-    }
-    if (isfinite(stats->max)) {
-        maxFitSigma = fabs(stats->max);
-    }
-
-    psF32 minFitSigma = 2.0;
-    if (isfinite(stats->clipSigma)) {
-        minFitSigma = fabs(stats->clipSigma);
-    }
-    if (isfinite(stats->min)) {
-        minFitSigma = fabs(stats->min);
-    }
-
-    tmpScalar.data.F32 = stats->fittedMean - minFitSigma*stats->fittedStdev;
-    if (tmpScalar.data.F32 < min) {
-        binMin = 0;
-    } else {
-        binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
-    }
-    tmpScalar.data.F32 = stats->fittedMean + maxFitSigma*stats->fittedStdev;
-    if (tmpScalar.data.F32 > max) {
-        binMax = histogram->bounds->n - 1;
-    } else {
-        binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
-    }
-    if (binMin < 0 || binMax < 0) {
-        psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma);
-        psFree(statsMinMax);
-        psFree(histogram);
-        psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
-        return false;
-    }
-
-    // Generate the variables that will be used in the Gaussian fitting
-    // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins
-    psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates
-    psArray *x = psArrayAlloc((1 + (binMax - binMin))); // Array of ordinates
-    for (long i = binMin, j = 0; i <= binMax ; i++, j++) {
-        y->data.F32[j] = histogram->nums->data.F32[i];
-        psVector *ordinate = psVectorAlloc(1, PS_TYPE_F32); // The ordinate value
-        ordinate->data.F32[0] = PS_BIN_MIDPOINT(histogram, i);
-        x->data[j] = ordinate;
-    }
-    if (psTraceGetLevel("psLib.math") >= 8) {
-        // XXX: Print the x array somehow.
-        PS_VECTOR_PRINT_F32(y);
-    }
-    psTrace("psLib.math", 6, "The clipped numBins is %ld\n", y->n);
-
-    // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
-    // XXX EAM : why not just fit the normalization as well??
-    p_psNormalizeVectorRange(y, 0.0, 1.0);
-
-    // Fit a Gaussian to the data.
-    psMinimization *minimizer = psMinimizationAlloc(100, 0.01); // The minimizer information
-    psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian
-    // Initial guess for the mean (index 0) and var (index 1).
-    params->data.F32[0] = stats->fittedMean;
-    params->data.F32[1] = PS_SQR(stats->fittedStdev);
-    if (psTraceGetLevel("psLib.math") >= 8) {
-        PS_VECTOR_PRINT_F32(params);
-        PS_VECTOR_PRINT_F32(y);
-    }
-    if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) {
-        psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
-        psFree(statsMinMax);
-        psFree(histogram);
-        psFree(x);
-        psFree(y);
-        psFree(minimizer);
-        psFree(params);
-        psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
-        return false;
-    }
-    if (psTraceGetLevel("psLib.math") >= 8) {
-        PS_VECTOR_PRINT_F32(params);
-    }
-
-    // The fitted mean is the Gaussian mean.
-    stats->fittedMean = params->data.F32[0];
-    psTrace("psLib.math", 6, "The fitted mean is %f.\n", stats->fittedMean);
-
-    // The fitted standard deviation
-    stats->fittedStdev = sqrt(params->data.F32[1]);
-    psTrace("psLib.math", 6, "The fitted stdev is %f.\n", stats->fittedStdev);
-
-    // Clean up after fitting
-    psFree (histogram);
-    psFree (x);
-    psFree (y);
-    psFree (minimizer);
-    psFree (params);
-    psFree (statsMinMax);
-    return true;
-}
-
 
 /******************************************************************************
@@ -1278 +596 @@
     }
 
+    // statsMinMax is only applied to a subset of the data points
     psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
     psHistogram *histogram = NULL;      // Histogram of the data
@@ -1287 +606 @@
     // Iterate to get the best bin size
     for (int iterate = 1; iterate > 0; iterate++) {
-        psTrace("psLib.math", 6,
-                "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n",
-                iterate);
+        psTrace("psLib.math", 6, "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n", iterate);
 
         // Get the minimum and maximum values
@@ -1308 +625 @@
         if (fabs(max - min) <= FLT_EPSILON) {
             psTrace("psLib.math", 5, "All data points have the same value: %f.\n", min);
-            if (stats->options & PS_STAT_ROBUST_MEDIAN) {
-                stats->robustMedian = min;
-            }
-            if (stats->options & PS_STAT_ROBUST_QUARTILE) {
-                stats->robustUQ = min;
-                stats->robustLQ = min;
-            }
+            stats->robustMedian = min;
+            stats->robustUQ = min;
+            stats->robustLQ = min;
             stats->robustN50 = numValid;
             psFree(statsMinMax);
@@ -1433 +746 @@
             return false;
         }
-        // XXXX we can probably just use the bin values to get sigma without interpolating...
-        // ADD step 4b: Interpolate Sigma to find these two positions exactly: these are the 1sigma positions.
-        #if 0
-        // XXX: I am overriding the ADD for now.  The following code implements the ADD exactly and is having
-        // problems fitting a 2nd-order polynomial to data y-values suchs as (1, 1, 100).  Therefore, I
-        // commented out the code and am doing simply linear interpolation.
-
-        // Value for the 15.8655% mark
-        float binLoF32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, binLo,
-                         totalDataPoints * 0.158655f);
-        // Value for the 84.1345% mark
-        float binHiF32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, binHi,
-                         totalDataPoints * 0.841345f);
-        psTrace("psLib.math", 6, "The exact 15.8655%% and 84.1345%% data point positions are: (%f, %f)\n",
-                binLoF32, binHiF32);
-        #else
-        // This code basically interpolates to find the positions exactly.
+
+        // ADD step 4b: Interpolate Sigma (linearly) to find these two positions exactly: these are the 1sigma positions.
         psTrace("psLib.math", 6, "binLo is %ld.  Nums at that bin and the next are (%.2f, %.2f)\n",
                 binLo, cumulative->nums->data.F32[binLo], cumulative->nums->data.F32[binLo+1]);
@@ -1455 +753 @@
                 binHi, cumulative->nums->data.F32[binHi], cumulative->nums->data.F32[binHi+1]);
 
-        float deltaBounds = cumulative->bounds->data.F32[binLo+1] - cumulative->bounds->data.F32[binLo];
-        float deltaNums;
-        float prevPixels;
+        float deltaBounds, deltaNums, prevPixels;
+        float percentNums, base, binLoF32;
+
+        // find the -1 sigma points with linear interpolation
+        deltaBounds = cumulative->bounds->data.F32[binLo+1] - cumulative->bounds->data.F32[binLo];
         if (binLo == 0) {
             deltaNums = cumulative->nums->data.F32[0];
@@ -1465 +765 @@
             prevPixels = cumulative->nums->data.F32[binLo-1];
         }
-        float percentNums = (totalDataPoints * 0.158655f) - prevPixels;
-        float base = cumulative->bounds->data.F32[binLo];
-        float binLoF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 15.8655% mark
+        percentNums = (totalDataPoints * 0.158655f) - prevPixels;
+        base = cumulative->bounds->data.F32[binLo];
+        binLoF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 15.8655% mark
         psTrace("psLib.math", 6,
                 "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
                 base, deltaBounds, deltaNums, prevPixels, percentNums);
 
+        // find the +1 sigma points with linear interpolation
         deltaBounds = cumulative->bounds->data.F32[binHi+1] - cumulative->bounds->data.F32[binHi];
         if (binHi == 0) {
@@ -1486 +787 @@
                 "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
                 base, deltaBounds, deltaNums, prevPixels, percentNums);
-        psTrace("psLib.math", 6,
+
+        // report +/- 1 sigma points
+        psTrace("psLib.math", 5,
                 "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n",
                 binLoF32, binHiF32);
-        #endif
 
         // ADD step 5: Determine SIGMA as 1/2 of the distance between these positions.
@@ -1585 +887 @@
     psTrace("psLib.math", 6, "The robustN50 is %ld.\n", N50);
 
-    // Fitted statistics
-    if (stats->options & PS_STAT_FITTED_MEAN || stats->options & PS_STAT_FITTED_STDEV) {
-        stats->fittedStdev = stats->robustStdev;  // pass the guess sigma
-        stats->fittedMean = stats->robustMedian;  // pass the guess mean
-        if (!vectorFittedStats(myVector, errors, mask, maskVal, stats)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to estimate statistics (pass 1)");
-            psFree(statsMinMax);
-            psFree(mask);
-
-            return false;
-        }
-        // if there is a large swing in sigma, try a second time
-        if ((stats->fittedStdev/stats->robustStdev) < 0.75) {
-            if (!vectorFittedStats (myVector, errors, mask, maskVal, stats)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to estimate statistics (pass 2)");
-                psFree(statsMinMax);
-                psFree(mask);
-
-                return false;
-            }
-        }
-    }
-
     // Clean up
     psFree(statsMinMax);
     psFree(mask);
 
+    stats->results |= PS_STAT_ROBUST_MEDIAN;
+    stats->results |= PS_STAT_ROBUST_STDEV;
+
     psTrace("psLib.math", 4, "---- %s(0) end  ----\n", __func__);
     return true;
 }
 
+/*
+ * vectorFittedStats requires guess for fittedMean and fittedStdev
+ * robustN50 should also be set
+ */
+static bool vectorFittedStats (const psVector* myVector,
+                               const psVector* errors,
+                               psVector* mask,
+                               psMaskType maskVal,
+                               psStats* stats)
+{
+
+    // This procedure requires the mean.  If it has not been already
+    // calculated, then call vectorSampleMean()
+    if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
+        vectorRobustStats(myVector, errors, mask, maskVal, stats);
+    }
+
+    // If the mean is NAN, then generate a warning and set the stdev to NAN.
+    if (isnan(stats->robustMedian)) {
+        stats->fittedStdev = NAN;
+        stats->fittedStdev = NAN;
+        psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
+        return false;
+    }
+
+    float guessStdev = stats->robustStdev;  // pass the guess sigma
+    float guessMean = stats->robustMedian;  // pass the guess mean
+
+    psTrace("psLib.math", 6, "The guess mean  is %f.\n", guessMean);
+    psTrace("psLib.math", 6, "The guess stdev is %f.\n", guessStdev);
+
+    bool done = false;
+    for (int iteration = 0; !done && (iteration < 2); iteration ++) {
+        psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
+
+        psScalar tmpScalar;                 // Temporary scalar variable, for p_psVectorBinDisect
+        tmpScalar.type.type = PS_TYPE_F32;
+
+        psF32 binSize = 1;
+        if (stats->options & PS_STAT_USE_BINSIZE) {
+            // Set initial bin size to the specified value.
+            binSize = stats->binsize;
+            psTrace("psLib.math", 6, "Setting initial robust bin size to %.2f\n", binSize);
+        } else {
+            // construct a histogram with (sigma/10 < binsize < sigma)
+            // set roughly so that the lowest bins have about 2 cnts
+            // Nsmallest ~ N50 / (4*dN))
+            psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8));
+            binSize = guessStdev / dN;
+        }
+
+        // Determine the min/max of the vector (which prior outliers masked out)
+        int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values
+        float min = statsMinMax->min;
+        float max = statsMinMax->max;
+        if (numValid == 0 || isnan(min) || isnan(max)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
+            psFree(statsMinMax);
+            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        // Calculate the number of bins.
+        // XXX can we calculate the binMin, binMax **before** building this histogram?
+        long numBins = (max - min) / binSize;
+        psTrace("psLib.math", 6, "The new min/max values are (%f, %f).\n", min, max);
+        psTrace("psLib.math", 6, "The new bin size is %f.\n", binSize);
+        psTrace("psLib.math", 6, "The numBins is %ld\n", numBins);
+
+        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
+        histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal);
+        if (psTraceGetLevel("psLib.math") >= 8) {
+            PS_VECTOR_PRINT_F32(histogram->nums);
+        }
+
+        long binMin = 0;                // Low bin to check
+        long binMax = 0;                // High bin to check
+
+        // Fit a Gaussian to the bins in the requested range about the guess mean
+        // min,max overrides clipSigma
+        psF32 maxFitSigma = 2.0;
+        if (isfinite(stats->clipSigma)) {
+            maxFitSigma = fabs(stats->clipSigma);
+        }
+        if (isfinite(stats->max)) {
+            maxFitSigma = fabs(stats->max);
+        }
+
+        psF32 minFitSigma = 2.0;
+        if (isfinite(stats->clipSigma)) {
+            minFitSigma = fabs(stats->clipSigma);
+        }
+        if (isfinite(stats->min)) {
+            minFitSigma = fabs(stats->min);
+        }
+
+        tmpScalar.data.F32 = guessMean - minFitSigma*guessStdev;
+        if (tmpScalar.data.F32 < min) {
+            binMin = 0;
+        } else {
+            binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
+        }
+        tmpScalar.data.F32 = guessMean + maxFitSigma*guessStdev;
+        if (tmpScalar.data.F32 > max) {
+            binMax = histogram->bounds->n - 1;
+        } else {
+            binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
+        }
+        if (binMin < 0 || binMax < 0) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma);
+            psFree(statsMinMax);
+            psFree(histogram);
+            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        // Generate the variables that will be used in the Gaussian fitting
+        // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins
+        psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates
+        psArray *x = psArrayAlloc((1 + (binMax - binMin))); // Array of ordinates
+        for (long i = binMin, j = 0; i <= binMax ; i++, j++) {
+            y->data.F32[j] = histogram->nums->data.F32[i];
+            psVector *ordinate = psVectorAlloc(1, PS_TYPE_F32); // The ordinate value
+            ordinate->data.F32[0] = PS_BIN_MIDPOINT(histogram, i);
+            x->data[j] = ordinate;
+        }
+        if (psTraceGetLevel("psLib.math") >= 8) {
+            // XXX: Print the x array somehow.
+            PS_VECTOR_PRINT_F32(y);
+        }
+        psTrace("psLib.math", 6, "The clipped numBins is %ld\n", y->n);
+        psTrace("psLib.math", 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin);
+        psTrace("psLib.math", 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax);
+
+        // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
+        // XXX EAM : why not just fit the normalization as well??
+        p_psNormalizeVectorRange(y, 0.0, 1.0);
+
+        // Fit a Gaussian to the data.
+        psMinimization *minimizer = psMinimizationAlloc(100, 0.01); // The minimizer information
+        psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian
+        // Initial guess for the mean (index 0) and var (index 1).
+        params->data.F32[0] = guessMean;
+        params->data.F32[1] = PS_SQR(guessStdev);
+        if (psTraceGetLevel("psLib.math") >= 8) {
+            PS_VECTOR_PRINT_F32(params);
+            PS_VECTOR_PRINT_F32(y);
+        }
+        if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
+            psFree(params);
+            psFree(minimizer);
+            psFree(x);
+            psFree(y);
+            psFree(histogram);
+            psFree(statsMinMax);
+            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+        if (psTraceGetLevel("psLib.math") >= 8) {
+            PS_VECTOR_PRINT_F32(params);
+        }
+
+        guessMean = params->data.F32[0];
+        guessStdev = sqrt(params->data.F32[1]);
+        if (guessStdev > 0.75*stats->robustStdev) {
+            done = true;
+        }
+
+        // Clean up after fitting
+        psFree (params);
+        psFree (minimizer);
+        psFree (x);
+        psFree (y);
+        psFree (histogram);
+        psFree (statsMinMax);
+    }
+
+    // The fitted mean is the Gaussian mean.
+    stats->fittedMean = guessMean;
+    psTrace("psLib.math", 6, "The fitted mean is %f.\n", stats->fittedMean);
+
+    // The fitted standard deviation
+    stats->fittedStdev = guessStdev;
+    psTrace("psLib.math", 6, "The fitted stdev is %f.\n", stats->fittedStdev);
+
+    stats->results |= PS_STAT_FITTED_MEAN;
+    stats->results |= PS_STAT_FITTED_STDEV;
+
+    return true;
+}
+
+
+/********************
+ * vectorFittedStats_v2 requires guess for fittedMean and fittedStdev
+ * robustN50 should also be set
+ * gaussian fit is performed using 2D polynomial to ln(y)
+ ********************/
+static bool vectorFittedStats_v2 (const psVector* myVector,
+                                  const psVector* errors,
+                                  psVector* mask,
+                                  psMaskType maskVal,
+                                  psStats* stats)
+{
+
+    // This procedure requires the mean.  If it has not been already
+    // calculated, then call vectorSampleMean()
+    if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
+        vectorRobustStats(myVector, errors, mask, maskVal, stats);
+    }
+
+    // If the mean is NAN, then generate a warning and set the stdev to NAN.
+    if (isnan(stats->robustMedian)) {
+        stats->fittedStdev = NAN;
+        stats->fittedStdev = NAN;
+        psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
+        return false;
+    }
+
+    float guessStdev = stats->robustStdev;  // pass the guess sigma
+    float guessMean = stats->robustMedian;  // pass the guess mean
+
+    psTrace("psLib.math", 6, "The guess mean  is %f.\n", guessMean);
+    psTrace("psLib.math", 6, "The guess stdev is %f.\n", guessStdev);
+
+    bool done = false;
+    for (int iteration = 0; !done && (iteration < 2); iteration ++) {
+        psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
+
+        psScalar tmpScalar;                 // Temporary scalar variable, for p_psVectorBinDisect
+        tmpScalar.type.type = PS_TYPE_F32;
+
+        psF32 binSize = 1;
+        if (stats->options & PS_STAT_USE_BINSIZE) {
+            // Set initial bin size to the specified value.
+            binSize = stats->binsize;
+            psTrace("psLib.math", 6, "Setting initial robust bin size to %.2f\n", binSize);
+        } else {
+            // construct a histogram with (sigma/10 < binsize < sigma)
+            // set roughly so that the lowest bins have about 2 cnts
+            // Nsmallest ~ N50 / (4*dN))
+            psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8));
+            binSize = guessStdev / dN;
+        }
+
+        // Determine the min/max of the vector (which prior outliers masked out)
+        int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values
+        float min = statsMinMax->min;
+        float max = statsMinMax->max;
+        if (numValid == 0 || isnan(min) || isnan(max)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
+            psFree(statsMinMax);
+            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        // Calculate the number of bins.
+        // XXX can we calculate the binMin, binMax **before** building this histogram?
+        long numBins = (max - min) / binSize;
+        psTrace("psLib.math", 6, "The new min/max values are (%f, %f).\n", min, max);
+        psTrace("psLib.math", 6, "The new bin size is %f.\n", binSize);
+        psTrace("psLib.math", 6, "The numBins is %ld\n", numBins);
+
+        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
+        histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal);
+        if (psTraceGetLevel("psLib.math") >= 8) {
+            PS_VECTOR_PRINT_F32(histogram->nums);
+        }
+
+        long binMin = 0;                // Low bin to check
+        long binMax = 0;                // High bin to check
+
+        // Fit a Gaussian to the bins in the requested range about the guess mean
+        // min,max overrides clipSigma
+        psF32 maxFitSigma = 2.0;
+        if (isfinite(stats->clipSigma)) {
+            maxFitSigma = fabs(stats->clipSigma);
+        }
+        if (isfinite(stats->max)) {
+            maxFitSigma = fabs(stats->max);
+        }
+
+        psF32 minFitSigma = 2.0;
+        if (isfinite(stats->clipSigma)) {
+            minFitSigma = fabs(stats->clipSigma);
+        }
+        if (isfinite(stats->min)) {
+            minFitSigma = fabs(stats->min);
+        }
+
+        tmpScalar.data.F32 = guessMean - minFitSigma*guessStdev;
+        if (tmpScalar.data.F32 < min) {
+            binMin = 0;
+        } else {
+            binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
+        }
+        tmpScalar.data.F32 = guessMean + maxFitSigma*guessStdev;
+        if (tmpScalar.data.F32 > max) {
+            binMax = histogram->bounds->n - 1;
+        } else {
+            binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
+        }
+        if (binMin < 0 || binMax < 0) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma);
+            psFree(statsMinMax);
+            psFree(histogram);
+            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        // Generate the variables that will be used in the Gaussian fitting
+        // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins
+        psVector *y = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates
+        psVector *x = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of ordinates
+        long j = 0;
+        for (long i = binMin; i <= binMax ; i++) {
+            if (histogram->nums->data.F32[i] <= 0.0)
+                continue;
+            x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i);
+            y->data.F32[j] = log(histogram->nums->data.F32[i]);
+            // XXX note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2)
+            j++;
+        }
+        y->n = x->n = j;
+        if (psTraceGetLevel("psLib.math") >= 8) {
+            // XXX: Print the x array somehow.
+            PS_VECTOR_PRINT_F32(y);
+        }
+        psTrace("psLib.math", 6, "The clipped numBins is %ld\n", y->n);
+        psTrace("psLib.math", 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin);
+        psTrace("psLib.math", 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax);
+
+        // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
+        psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
+
+        // XXX how can we protect against some extreme outliers?  the robust histogram
+        // should have already dealt with those, but we are again sensitive to them...
+        // psStats *fitStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+        // fitStats->clipIter = 3.0;
+        // fitStats->clipSigma = 3.0;
+        // psVector *fitMask = psVectorAlloc(y->n, PS_TYPE_U8);
+        // psVectorInit (fitMask, 0);
+
+        if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
+            psFree(x);
+            psFree(y);
+            psFree(poly);
+            // psFree(fitStats);
+            // psFree(fitMask);
+            psFree(histogram);
+            psFree(statsMinMax);
+            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        guessStdev = sqrt(-0.5/poly->coeff[2]);
+        guessMean = poly->coeff[1]*PS_SQR(guessStdev);
+        if (guessStdev > 0.75*stats->robustStdev) {
+            done = true;
+        } else {
+            psTrace("psLib.math", 6, "The new guess mean  is %f.\n", guessMean);
+            psTrace("psLib.math", 6, "The new guess stdev is %f.\n", guessStdev);
+        }
+
+        // Clean up after fitting
+        psFree (x);
+        psFree (y);
+        psFree (poly);
+        // psFree (fitStats);
+        // psFree (fitMask);
+        psFree (histogram);
+        psFree (statsMinMax);
+    }
+
+    // The fitted mean is the Gaussian mean.
+    stats->fittedMean = guessMean;
+    psTrace("psLib.math", 6, "The fitted mean is %f.\n", stats->fittedMean);
+
+    // The fitted standard deviation
+    stats->fittedStdev = guessStdev;
+    psTrace("psLib.math", 6, "The fitted stdev is %f.\n", stats->fittedStdev);
+
+    stats->results |= PS_STAT_FITTED_MEAN_V2;
+    stats->results |= PS_STAT_FITTED_STDEV_V2;
+
+    return true;
+}
+
+
+/******************************************************************************
+vectorSmoothHistGaussian(): This routine smoothes the data in the input
+robustHistogram with a Gaussian of width sigma.  It returns a psVector of the
+smoothed data.
+ 
+XXX this function is unused
+*****************************************************************************/
+psVector *vectorSmoothHistGaussian(psHistogram *histogram,
+                                   psF32 sigma)
+{
+    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
+    psTrace("psLib.math", 5, "(histogram->nums->n, sigma) is (%d, %.2f\n", (int) histogram->nums->n, sigma);
+    PS_ASSERT_PTR_NON_NULL(histogram, NULL);
+    PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL);
+    PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL);
+    if (psTraceGetLevel("psLib.math") >= 8) {
+        PS_VECTOR_PRINT_F32(histogram->nums);
+    }
+
+    long numBins = histogram->nums->n;  // Number of histogram bins
+    psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); // Smoothed version of histogram bins
+    const psVector *bounds = histogram->bounds; // The bounds for the histogram bins
+
+    if (!histogram->uniform) {
+        //
+        // We get here if the histogram is non-uniform.  This code is not tested.
+        // However, it is also not used anywhere, yet.
+        //
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSmoothHistGaussian() on non-uniform histograms has not been tested or used.\n");
+
+        for (long i = 0; i < numBins; i++) {
+            // Determine the midpoint of bin i.
+            float iMid = PS_BIN_MIDPOINT(histogram, i);
+
+            //
+            // We determine the bin numbers (jMin:jMax) corresponding to a
+            // range of data values surrounding iMid.  The range is of size:
+            // 2*PS_GAUSS_WIDTH*sigma
+            long jMin, jMax;
+            psF32 x = iMid - (PS_GAUSS_WIDTH * sigma);
+            for (jMin = i; jMin > 0 && bounds->data.F32[jMin] > x; jMin--)
+                ;
+            x = iMid + (PS_GAUSS_WIDTH * sigma);
+            for (jMax = i; jMax < bounds->n - 1 && bounds->data.F32[jMax + 1] > x; jMax++)
+                ;
+
+            //
+            // Loop from jMin to jMax, computing the gaussian of data i.
+            //
+            smooth->data.F32[i] = 0.0;
+            for (long j = jMin ; j <= jMax ; j++) {
+                float jMid = PS_BIN_MIDPOINT(histogram, j);
+                smooth->data.F32[i] +=
+                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
+            }
+        }
+    } else {
+        //
+        // We get here if the histogram is uniform.
+        //
+        for (long i = 0; i < numBins; i++) {
+            psF32 binSize = bounds->data.F32[1] - bounds->data.F32[0];
+            psS32 gaussWidth = ((PS_GAUSS_WIDTH * sigma) / binSize);
+
+            //
+            // XXX: The following is wrong.  However, in practice, the sigma was too
+            // large compared to the binSize.  This meant that the smoothing was done
+            // over 500 bins in the robust stats algorithm.  This mean that the smoothing
+            // took way too long.
+            //
+            #if 0
+
+            gaussWidth = 10;
+            #endif
+            //
+            // We determine the bin numbers (jMin:jMax) corresponding to a
+            // range of data values surrounding iMid.  The range is of size:
+            // 2*PS_GAUSS_WIDTH*sigma
+            //
+            psS32 jMin = PS_MAX(i - gaussWidth, 0);
+            psS32 jMax = PS_MIN(i + gaussWidth, bounds->n - 1);
+
+            //
+            // Loop from jMin to jMax, computing the gaussian of data i.
+            //
+            smooth->data.F32[i] = 0.0;
+            float iMid = PS_BIN_MIDPOINT(histogram, i);
+            for (long j = jMin ; j <= jMax ; j++) {
+                float jMid = PS_BIN_MIDPOINT(histogram, j);
+                smooth->data.F32[i] +=
+                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
+            }
+        }
+    }
+
+    if (psTraceGetLevel("psLib.math") >= 8) {
+        PS_VECTOR_PRINT_F32(smooth);
+    }
+    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
+    return(smooth);
+}
+
 /*****************************************************************************/
 
@@ -1621 +1395 @@
 
 /*****************************************************************************/
-
-static void histogramFree(psHistogram* myHist);
 
 // We keep statsFree so that we can identify statistics pointers from the memblock
@@ -1632 +1404 @@
 /******************************************************************************
     psStatsAlloc(): This routine must create a new psStats data structure.
- *****************************************************************************/
+*****************************************************************************/
 psStats* psStatsAlloc(psStatsOptions options)
 {
@@ -1675 +1447 @@
 }
 
-
-
-/******************************************************************************
-psHistogramAlloc(lower, upper, n): allocate a uniform histogram structure
-with the specifed upper and lower limits, and the specifed number of bins.
-This routine will also set the bounds for each of the bins.
- 
-Input:
-    lower
-    upper
-    n
-Returns:
-    The histogram structure
- *****************************************************************************/
-psHistogram* psHistogramAlloc(float lower, float upper, int n)
-{
-    psTrace("psLib.math", 3, "---- %s() begin  ----\n", __func__);
-    psTrace("psLib.math", 5, "(lower, upper, n) is (%f, %f, %d)\n", lower, upper, n);
-    PS_ASSERT_INT_POSITIVE(n, NULL);
-    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(upper, lower, NULL);
-
-    // Allocate memory for the new histogram structure.  If there are N bins, then there are N+1 bounds to
-    // those bins.
-    psHistogram *newHist = (psHistogram* ) psAlloc(sizeof(psHistogram)); // The new histogram structure
-    psMemSetDeallocator(newHist, (psFreeFunc) histogramFree);
-    psVector* newBounds = psVectorAlloc(n + 1, PS_TYPE_F32);
-    newHist->bounds = newBounds;
-
-    // Calculate the bounds for each bin.
-    psF32 binSize = (upper - lower) / (psF32)n; // The histogram bin size
-    // XXX: Is the following necessary? It prevents the max data point from being in a non-existant bin.
-    binSize += FLT_EPSILON;
-    for (long i = 0; i < n + 1; i++) {
-        newBounds->data.F32[i] = lower + (binSize * (psF32)i);
-    }
-
-    // Allocate the bins, and initialize them to zero.
-    newHist->nums = psVectorAlloc(n, PS_TYPE_F32);
-    psVectorInit(newHist->nums, 0.0);
-
-    // Initialize the other members.
-    newHist->minNum = 0;
-    newHist->maxNum = 0;
-    newHist->uniform = true;
-
-    psTrace("psLib.math", 3, "---- %s() end  ----\n", __func__);
-    return newHist;
-}
-
-/******************************************************************************
-psHistogramAllocGeneric(bounds): allocate a non-uniform histogram structure
-with the specifed bounds.
- 
-Input:
-    bounds
-Returns:
-    The histogram structure
- *****************************************************************************/
-psHistogram* psHistogramAllocGeneric(const psVector* bounds)
-{
-    psTrace("psLib.math", 3, "---- %s() begin  ----\n", __func__);
-    PS_ASSERT_VECTOR_NON_NULL(bounds, NULL);
-    PS_ASSERT_VECTOR_TYPE(bounds, PS_TYPE_F32, NULL);
-    PS_ASSERT_LONG_LARGER_THAN_OR_EQUAL(bounds->n, (long)2, NULL);
-
-    // Allocate memory for the new histogram structure.
-    psHistogram *newHist = (psHistogram* ) psAlloc(sizeof(psHistogram)); // The new histogram structure
-    psMemSetDeallocator(newHist, (psFreeFunc) histogramFree);
-    psVector* newBounds = psVectorCopy(NULL, bounds, PS_TYPE_F32);
-    newHist->bounds = newBounds;
-
-    // Allocate the bins, and initialize them to zero.  If there are N bounds,
-    // then there are N-1 bins.
-    newHist->nums = psVectorAlloc((bounds->n) - 1, PS_TYPE_F32);
-    psVectorInit(newHist->nums, 0.0);
-
-    // Initialize the other members.
-    newHist->minNum = 0;
-    newHist->maxNum = 0;
-    newHist->uniform = false;
-
-    psTrace("psLib.math", 3, "---- %s() end  ----\n", __func__);
-    return (newHist);
-}
-
-static void histogramFree(psHistogram* myHist)
-{
-    psFree(myHist->bounds);
-    psFree(myHist->nums);
-}
-
-
-bool psMemCheckHistogram(psPtr ptr)
-{
-    PS_ASSERT_PTR(ptr, false);
-    return ( psMemGetDeallocator(ptr) == (psFreeFunc)histogramFree );
-}
-
-
-
-/*****************************************************************************
-UpdateHistogramBins(binNum, out, data, error): This routine is to be used when
-updating the histogram in the presence of errors in the input data.  We treat
-the data point as a boxcar PDF and update a range of points surrounding the
-histogram bin which contains the point.  The width of that boxcar is defined
-as 2.35 * error.
- 
-XXX: Must test this.
- *****************************************************************************/
-static bool UpdateHistogramBins(long binNum, // Bin number of the data point
-                                psHistogram* out, // The histogram to be updated
-                                psF32 data, // The data point value
-                                psF32 error // The error in the data point
-                               )
-{
-    psTrace("psLib.math", 3, "---- %s() begin  ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(out, false);
-    PS_ASSERT_PTR_NON_NULL(out->bounds, false);
-    PS_ASSERT_PTR_NON_NULL(out->nums, false);
-    PS_ASSERT_LONG_WITHIN_RANGE(binNum, (long)0, (long)((out->nums->n)-1), false);
-    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(error, 0.0, false);
-    PS_ASSERT_FLOAT_WITHIN_RANGE(data, out->bounds->data.F32[0],
-                                 out->bounds->data.F32[(out->bounds->n)-1], false);
-
-    psF32 boxcarWidth = 2.35 * error;   // Width of the boxcar
-    psF32 boxcarCenter = (out->bounds->data.F32[binNum] +
-                          out->bounds->data.F32[binNum+1]) / 2.0; // Centre of the boxcar
-    psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0); // Left endpoint of the boxcar for the PDF
-    psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0); // Right endpoint of the boxcar for the PDF
-    psS32 boxcarLeftBinNum = 0;         // Bin number for left endpoint
-    psS32 boxcarRightBinNum = 0;        // Bin number for right endpoint
-
-    // Determine the left endpoint of the boxcar for the PDF.
-    for (long bin = binNum; bin >= 0; bin--) {
-        if (out->nums->data.F32[bin] <= boxcarLeft) {
-            boxcarLeftBinNum = bin;
-            break;
-        }
-    }
-
-    // Determine the right endpoint of the boxcar for the PDF.
-    for (long bin = binNum; bin < out->nums->n; bin++) {
-        if (out->nums->data.F32[bin] >= boxcarRight) {
-            boxcarRightBinNum = bin;
-            break;
-        }
-    }
-
-    // If the boxcar fits entirely inside this bin, then simply add 1.0 to the
-    // bin and return.
-    if (boxcarLeftBinNum == boxcarRightBinNum) {
-        out->nums->data.F32[binNum] += 1.0;
-        psTrace("psLib.math", 3, "---- %s(true) end  ----\n", __func__);
-        return true;
-    }
-
-    // If we get here, multiple bins must be updated.  We handle the left-most endpoint, and right-most
-    // endpoints differently.
-    out->nums->data.F32[boxcarLeftBinNum] +=
-        (out->bounds->data.F32[boxcarLeftBinNum+1] - boxcarLeft) / boxcarWidth;
-
-    // Loop through the center bins, if any.
-    for (long bin = boxcarLeftBinNum + 1; bin < (boxcarRightBinNum - 1); bin++) {
-        out->nums->data.F32[bin] +=
-            (out->bounds->data.F32[bin+1] - out->bounds->data.F32[bin]) / boxcarWidth;
-    }
-
-    // Handle the right endpoint differently.
-    out->nums->data.F32[boxcarRightBinNum]+=
-        (boxcarRight - out->bounds->data.F32[boxcarRightBinNum]) / boxcarWidth;
-
-    psTrace("psLib.math", 3, "---- %s(true) end  ----\n", __func__);
-    return true;
-}
-
-
-/*****************************************************************************
-psVectorHistogram(out, in, errors, mask, maskVal): this procedure takes as
-input a preallocated and initialized histogram structure.  It fills the bins
-in that histogram structure in accordance with the input data "in" and the,
-possibly NULL, mask vector.
- 
-Inputs:
-    out
-    in
-    mask
-    maskVal
-Returns:
-    The histogram structure "out".
- *****************************************************************************/
-psHistogram* psVectorHistogram(psHistogram* out,
-                               const psVector* values,
-                               const psVector* errors,
-                               const psVector* mask,
-                               psMaskType maskVal)
-{
-    psTrace("psLib.math", 3, "---- %s() begin  ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(out, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(out->bounds, NULL);
-    PS_ASSERT_VECTOR_TYPE(out->bounds, PS_TYPE_F32, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(out->bounds->n, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(out->nums, NULL);
-    PS_ASSERT_VECTOR_TYPE(out->nums, PS_TYPE_F32, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(out->nums->n, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(values, out);
-    if (mask) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(values, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    if (errors) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(values, errors, NULL);
-        PS_ASSERT_VECTOR_TYPE(errors, values->type.type, NULL);
-    }
-
-    long binNum = 0;                    // A temporary bin number
-    long numBins = out->nums->n;        // The total number of bins
-    psScalar tmpScalar;
-    tmpScalar.type.type = PS_TYPE_F32;
-
-    // Convert input and errors vectors to F32 if necessary.
-    psVector* inF32 = NULL;             // F32 version of input vector
-    if (values->type.type == PS_TYPE_F32) {
-        inF32 = psMemIncrRefCounter((psPtr)values);
-    } else {
-        inF32 = psVectorCopy(NULL, values, PS_TYPE_F32);
-    }
-    psVector* errorsF32 = NULL;         // F32 version of errors vector
-    if (errors) {
-        if (errors->type.type == PS_TYPE_F32) {
-            errorsF32 = psMemIncrRefCounter((psPtr)errors);
-        } else {
-            errorsF32 = psVectorCopy(NULL, errors, PS_TYPE_F32);
-        }
-    }
-
-    for (long i = 0; i < inF32->n; i++) {
-        // Check if this pixel is masked, and if so, skip it.
-        if (!mask || (mask && (!(mask->data.U8[i] & maskVal)))) {
-            if (inF32->data.F32[i] < out->bounds->data.F32[0]) {
-                // If this pixel is below minimum value, count it, then skip.
-                out->minNum++;
-            } else if (inF32->data.F32[i] > out->bounds->data.F32[numBins]) {
-                // If this pixel is above maximum value, count it, then skip.
-                out->maxNum++;
-            } else {
-                // If this is a uniform histogram, determining the correct
-                // number is trivial.
-                if (out->uniform == true) {
-                    float binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; // Histogram bin size
-                    binNum = (inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize;
-                    if (errorsF32) {
-                        if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errorsF32->data.F32[i])) {
-                            psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram "
-                                     "bins with the errors vector.\n");
-                        }
-                    } else {
-                        // This if-statement really shouldn't be necessary.
-                        // However, due to numerical lack of precision, we
-                        // occasionally produce a binNum outside the range.
-                        if (binNum >= out->nums->n) {
-                            binNum = out->nums->n - 1;
-                        }
-                        (out->nums->data.F32[binNum])+= 1.0;
-                    }
-
-                } else {
-                    // If this is a non-uniform histogram, determining the
-                    // correct bin number requires a bit more work.
-                    tmpScalar.data.F32 = inF32->data.F32[i];
-                    binNum = p_psVectorBinDisect( *(psVector* *)&out->bounds, &tmpScalar);
-                    if (binNum < 0) {
-                        psLogMsg(__func__, PS_LOG_WARN,
-                                 "WARNING: psVectorHistogram(): element outside histogram bounds.\n");
-                    } else {
-                        if (errorsF32 != NULL) {
-                            if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errors->data.F32[i])) {
-                                psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram "
-                                         "bins with the errors vector.\n");
-                            }
-                        } else {
-                            out->nums->data.F32[binNum] += 1.0;
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-    psFree(inF32);
-    psFree(errorsF32);
-
-    psTrace("psLib.math", 3, "---- %s() end  ----\n", __func__);
-    return (out);
-}
-
 /******************************************************************************
 psVectorStats(in, mask, maskVal, stats): this is the public API
@@ -1984 +1461 @@
  
 XXX: Should we free stats if the asserts fail? NO; we don't own it (RHL).
- *****************************************************************************/
-psStats* psVectorStats(
-    psStats* stats,
-    const psVector* in,
-    const psVector* errors,
-    const psVector* mask,
-    psMaskType maskVal)
+*****************************************************************************/
+bool psVectorStats(psStats* stats,
+                   const psVector* in,
+                   const psVector* errors,
+                   const psVector* mask,
+                   psMaskType maskVal)
 {
     psTrace("psLib.math", 3,"---- %s() begin  ----\n", __func__);
@@ -2036 +1512 @@
     }
 
+    // init the value of stats->results: this is used internally to check if
+    // prior functions have been called
+    stats->results = PS_STAT_NONE;
+    bool status = true;
+
     // ************************************************************************
     if (stats->options & PS_STAT_SAMPLE_MEAN) {
         if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMean() returned an error.\n");
-            stats->sampleMean = NAN;
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector sample mean");
+            status &= false;
         }
     }
@@ -2047 +1528 @@
     if (stats->options & (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE)) {
         if (!vectorSampleMedian(inF32, maskU8, maskVal, stats)) {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMedian() returned an error.\n");
-            stats->sampleMedian = NAN;
-            stats->sampleUQ = NAN;
-            stats->sampleLQ = NAN;
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample median");
+            status &= false;
         }
     }
@@ -2056 +1535 @@
     // ************************************************************************
     if (stats->options & PS_STAT_SAMPLE_STDEV) {
-        if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMean() returned an error.\n");
-            stats->sampleMean = NAN;
-        } else {
-            vectorSampleStdev(inF32, errorsF32, maskU8, maskVal, stats);
-        }
-    }
-
-    // ************************************************************************
-    if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE |
-                          PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
-        if (!vectorRobustStats(inF32, errorsF32, maskU8, maskVal, stats)) {
-            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate the specified statistic."));
-            psFree(inF32);
-            psFree(errorsF32);
-            psTrace("psLib.math", 3,"---- %s(NULL) end  ----\n", __func__);
-            return(NULL);
-        }
-    }
-
-    // ************************************************************************
-    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
-        if (!vectorClippedStats(inF32, errorsF32, maskU8, maskVal, stats)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics for input psVector.\n");
-            stats->clippedMean = NAN;
-            stats->clippedStdev = NAN;
+        if (!vectorSampleStdev(inF32, errorsF32, maskU8, maskVal, stats)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample stdev");
+            status &= false;
         }
     }
@@ -2089 +1545 @@
         if (vectorMinMax(inF32, maskU8, maskVal, stats) == 0) {
             psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum and maximum");
+            status &= false;
+        }
+    }
+
+    // ************************************************************************
+    if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE)) {
+        if (!vectorRobustStats(inF32, errorsF32, maskU8, maskVal, stats)) {
+            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate robust statistics"));
+            status &= false;
+        }
+    }
+
+    // ************************************************************************
+    if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
+        if (!vectorFittedStats(inF32, errorsF32, maskU8, maskVal, stats)) {
+            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
+            status &= false;
+        }
+    }
+
+    // ************************************************************************
+    if (stats->options & (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2)) {
+        if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
+            psAbort ("stats", "you may not specify both FITTED_MEAN and FITTED_MEAN_V2");
+        }
+        if (!vectorFittedStats_v2(inF32, errorsF32, maskU8, maskVal, stats)) {
+            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
+            status &= false;
+        }
+    }
+
+    // ************************************************************************
+    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
+        if (!vectorClippedStats(inF32, errorsF32, maskU8, maskVal, stats)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics\n");
+            status &= false;
         }
     }
@@ -2096 +1588 @@
     psFree(maskU8);
     psTrace("psLib.math", 3,"---- %s() end  ----\n", __func__);
-    return (stats);
+    return (status);
 }
 
@@ -2106 +1598 @@
     }
 
-    READ_STAT("MEAN",     PS_STAT_SAMPLE_MEAN);
-    READ_STAT("STDEV",    PS_STAT_SAMPLE_STDEV);
-    READ_STAT("MEDIAN",   PS_STAT_SAMPLE_MEDIAN);
-    READ_STAT("QUARTILE", PS_STAT_SAMPLE_QUARTILE);
+    READ_STAT("MEAN",       PS_STAT_SAMPLE_MEAN);
+    READ_STAT("STDEV",      PS_STAT_SAMPLE_STDEV);
+    READ_STAT("MEDIAN",     PS_STAT_SAMPLE_MEDIAN);
+    READ_STAT("QUARTILE",   PS_STAT_SAMPLE_QUARTILE);
     READ_STAT("SAMPLE_MEAN",     PS_STAT_SAMPLE_MEAN);
     READ_STAT("SAMPLE_STDEV",    PS_STAT_SAMPLE_STDEV);
@@ -2118 +1610 @@
     READ_STAT("ROBUST_STDEV",    PS_STAT_ROBUST_STDEV);
     READ_STAT("ROBUST_QUARTILE", PS_STAT_ROBUST_QUARTILE);
-    READ_STAT("FITTED",       PS_STAT_FITTED_MEAN);
-    READ_STAT("FITTED_MEAN",  PS_STAT_FITTED_MEAN);
-    READ_STAT("FITTED_STDEV", PS_STAT_ROBUST_STDEV);
-    READ_STAT("CLIPPED",       PS_STAT_CLIPPED_MEAN);
-    READ_STAT("CLIPPED_MEAN",  PS_STAT_CLIPPED_MEAN);
-    READ_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV);
+    READ_STAT("FITTED",         PS_STAT_FITTED_MEAN);
+    READ_STAT("FITTED_MEAN",    PS_STAT_FITTED_MEAN);
+    READ_STAT("FITTED_STDEV",   PS_STAT_FITTED_STDEV);
+    READ_STAT("FITTED_V2",       PS_STAT_FITTED_MEAN_V2);
+    READ_STAT("FITTED_MEAN_V2",  PS_STAT_FITTED_MEAN_V2);
+    READ_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2);
+    READ_STAT("CLIPPED",         PS_STAT_CLIPPED_MEAN);
+    READ_STAT("CLIPPED_MEAN",    PS_STAT_CLIPPED_MEAN);
+    READ_STAT("CLIPPED_STDEV",   PS_STAT_CLIPPED_STDEV);
 
     psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to parse statistic: %s\n", string);
@@ -2146 +1641 @@
     WRITE_STAT("ROBUST_STDEV",    PS_STAT_ROBUST_STDEV);
     WRITE_STAT("ROBUST_QUARTILE", PS_STAT_ROBUST_QUARTILE);
-    WRITE_STAT("FITTED_MEAN",  PS_STAT_FITTED_MEAN);
-    WRITE_STAT("FITTED_STDEV", PS_STAT_ROBUST_STDEV);
-    WRITE_STAT("CLIPPED_MEAN",  PS_STAT_CLIPPED_MEAN);
-    WRITE_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV);
+    WRITE_STAT("FITTED_MEAN",     PS_STAT_FITTED_MEAN);
+    WRITE_STAT("FITTED_STDEV",    PS_STAT_FITTED_STDEV);
+    WRITE_STAT("FITTED_MEAN_V2",  PS_STAT_FITTED_MEAN_V2);
+    WRITE_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2);
+    WRITE_STAT("CLIPPED_MEAN",    PS_STAT_CLIPPED_MEAN);
+    WRITE_STAT("CLIPPED_STDEV",   PS_STAT_CLIPPED_STDEV);
 
     return string;
@@ -2197 +1694 @@
     case PS_STAT_FITTED_MEAN:
     case PS_STAT_FITTED_STDEV:
+    case PS_STAT_FITTED_MEAN_V2:
+    case PS_STAT_FITTED_STDEV_V2:
     case PS_STAT_CLIPPED_MEAN:
     case PS_STAT_CLIPPED_STDEV:
@@ -2227 +1726 @@
     case PS_STAT_FITTED_STDEV:
         return stats->fittedStdev;
+    case PS_STAT_FITTED_MEAN_V2:
+        return stats->fittedMean;
+    case PS_STAT_FITTED_STDEV_V2:
+        return stats->fittedStdev;
     case PS_STAT_CLIPPED_MEAN:
         return stats->clippedMean;
@@ -2246 +1749 @@
     return NAN;
 }
+
+// other private functions used above
+
+static psF32 QuadraticInverse(psF32 a,
+                              psF32 b,
+                              psF32 c,
+                              psF32 y,
+                              psF32 xLo,
+                              psF32 xHi
+                             )
+{
+    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;
+
+    if (xLo <= x1 && x1 <= xHi) {
+        return x1;
+    }
+    if (xLo <= x2 && x2 <= xHi) {
+        return x2;
+    }
+    return 0.5 * (xLo + xHi);
+}
+
+
+/******************************************************************************
+fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): A general
+routine which fits a quadratic to three points and returns the x-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 (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).
+ 
+*****************************************************************************/
+static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec,
+        psVector *yVec,
+        psS32 binNum,
+        psF32 yVal
+                                              )
+{
+    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
+    psTrace("psLib.math", 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_VECTORS_SIZE_EQUAL(xVec, yVec, 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] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
+        x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
+        x->data.F64[2] = (psF64) (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];
+        psTrace("psLib.math", 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("psLib.math", 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
+        psTrace("psLib.math", 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]);
+        }
+
+        //
+        // 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);
+            psTrace("psLib.math", 5, "---- %s() end ----\n", __func__);
+            return NAN;
+        }
+
+        // Determine the coefficients of the polynomial.
+        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
+        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x);
+        if (myPoly == NULL) {
+            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);
+            psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__);
+            return NAN;
+        }
+        psTrace("psLib.math", 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
+        psTrace("psLib.math", 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
+        psTrace("psLib.math", 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
+        psTrace("psLib.math", 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("psLib.math", 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
+        tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal,
+                                    x->data.F64[0], x->data.F64[2]);
+        psFree(myPoly);
+
+        if (isnan(tmpFloat)) {
+            psError(PS_ERR_UNEXPECTED_NULL,
+                    false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
+            psFree(x);
+            psFree(y);
+            psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__);
+            return(NAN);
+        }
+    } 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.
+            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
+                              xVec->data.F32[binNum + 1]);
+        } else if (binNum == (xVec->n - 1)) {
+            // The special case where we have two points only at the end of
+            // the vectors x and y.
+            // 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]);
+        }
+    }
+
+    psTrace("psLib.math", 6, "FIT: return %f\n", tmpFloat);
+    psFree(x);
+    psFree(y);
+
+    psTrace("psLib.math", 5, "---- %s(%f) end ----\n", __func__, tmpFloat);
+    return tmpFloat;
+}
+
+/******************************************************************************
+NOTE: We assume unnormalized gaussians.
+*****************************************************************************/
+static psF32 minimizeLMChi2Gauss1D(psVector *deriv,
+                                   const psVector *params,
+                                   const psVector *coords
+                                  )
+{
+    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
+    PS_ASSERT_VECTOR_SIZE(params, (long)2, NAN);
+    PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(coords, NAN);
+    PS_ASSERT_VECTOR_SIZE(coords, (long)1, NAN);
+    PS_ASSERT_VECTOR_TYPE(coords, PS_TYPE_F32, NAN);
+
+    psF32 x = coords->data.F32[0];
+    psF32 mean = params->data.F32[0];
+    psF32 var = params->data.F32[1];
+    psF32 dx = (x - mean);
+
+    psF32 gauss = exp (-0.5*PS_SQR(dx)/var);
+    if (deriv) {
+        PS_ASSERT_VECTOR_SIZE(deriv, (long)2, NAN);
+        PS_ASSERT_VECTOR_TYPE(deriv, PS_TYPE_F32, NAN);
+        psF32 tmp = dx * gauss;
+        deriv->data.F32[0] = tmp / var;
+        deriv->data.F32[1] = tmp * dx / (var * var);
+    }
+
+
+    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
+    return gauss;
+}