Changeset 42719 for branches/2dbias/psLib/src/math/psStats.c

Timestamp:

Sep 13, 2024, 10:19:51 PM (22 months ago)

Author:

hgao

Message:

finetune overscan region and stats for reducing and smoothing overscan array/vector

File:

• branches/2dbias/psLib/src/math/psStats.c (modified) (140 diffs)

branches/2dbias/psLib/src/math/psStats.c

@@ -19 +19 @@
  */
 
-
 // Note: choice of return values in many places is quite grey.  For example, calculating the sample standard
 // deviation: here it's an error to get an array of size zero, but it's not an error to get an array of size
@@ -63 +62 @@
 #include "psString.h"
 
-
-
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
 /*****************************************************************************/
-#define PS_GAUSS_WIDTH 3                // The width of the Gaussian smoothing.
+#define PS_GAUSS_WIDTH 3         // The width of the Gaussian smoothing.
 #define PS_CLIPPED_NUM_ITER_LB 1 // This corresponds to N in the ADD.
 #define PS_CLIPPED_NUM_ITER_UB 10
@@ -77 +74 @@
 #define TRACE "psLib.math"
 
-#define MASK_MARK 0x80   // XXX : can we change this? bit to use internally to mark data as bad
-#define PS_ROBUST_MAX_ITERATIONS 20     // Maximum number of iterations for robust statistics
+#define MASK_MARK 0x80              // XXX : can we change this? bit to use internally to mark data as bad
+#define PS_ROBUST_MAX_ITERATIONS 20 // Maximum number of iterations for robust statistics
 
 #define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
-(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
+    (0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM) + 1]))
 
 // set the bin closest to the corresponding value.  if USE_END is +/- 1,
 // out-of-range saturates on lower/upper bin REGARDLESS of actual value
-#define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE, USE_END) { \
-        psVectorBinaryDisectResult result; \
-        psScalar tmpScalar; \
-        tmpScalar.type.type = PS_TYPE_F32; \
-        tmpScalar.data.F32 = (VALUE); \
-        RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar); \
-        switch (result) { \
-          case PS_BINARY_DISECT_PASS: \
-            break; \
-          case PS_BINARY_DISECT_OUTSIDE_RANGE: \
-            psTrace(TRACE, 6, "selected bin outside range"); \
-            if (USE_END == -1) { RESULT = 0; } \
-            if (USE_END == +1) { RESULT = VECTOR->n - 1; } \
-            break; \
-          case PS_BINARY_DISECT_INVALID_INPUT: \
-          case PS_BINARY_DISECT_INVALID_TYPE: \
-            psAbort ("programming error"); \
-            break; \
-        } }
-
-# define PS_BIN_INTERPOLATE(RESULT, VECTOR, BOUNDS, BIN, VALUE) { \
-        float dX, dY, Xo, Yo, Xt; \
-        if (BIN == BOUNDS->n - 1) { \
-            dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]); \
-            dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1]; \
-            Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \
-            Yo = VECTOR->data.F32[BIN]; \
-        } else { \
-            dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]); \
-            dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN]; \
-            Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \
-            Yo = VECTOR->data.F32[BIN]; \
-        } \
-        if (dY != 0) { \
-            Xt = (VALUE - Yo)*dX/dY + Xo; \
-        } else { \
-            Xt = Xo; \
-        } \
-        Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \
+#define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE, USE_END)            \
+    {                                                               \
+        psVectorBinaryDisectResult result;                          \
+        psScalar tmpScalar;                                         \
+        tmpScalar.type.type = PS_TYPE_F32;                          \
+        tmpScalar.data.F32 = (VALUE);                               \
+        RESULT = psVectorBinaryDisect(&result, VECTOR, &tmpScalar); \
+        switch (result)                                             \
+        {                                                           \
+        case PS_BINARY_DISECT_PASS:                                 \
+            break;                                                  \
+        case PS_BINARY_DISECT_OUTSIDE_RANGE:                        \
+            psTrace(TRACE, 6, "selected bin outside range");        \
+            if (USE_END == -1)                                      \
+            {                                                       \
+                RESULT = 0;                                         \
+            }                                                       \
+            if (USE_END == +1)                                      \
+            {                                                       \
+                RESULT = VECTOR->n - 1;                             \
+            }                                                       \
+            break;                                                  \
+        case PS_BINARY_DISECT_INVALID_INPUT:                        \
+        case PS_BINARY_DISECT_INVALID_TYPE:                         \
+            psAbort("programming error");                           \
+            break;                                                  \
+        }                                                           \
+    }
+
+#define PS_BIN_INTERPOLATE(RESULT, VECTOR, BOUNDS, BIN, VALUE)                             \
+    {                                                                                      \
+        float dX, dY, Xo, Yo, Xt;                                                          \
+        if (BIN == BOUNDS->n - 1)                                                          \
+        {                                                                                  \
+            dX = 0.5 * (BOUNDS->data.F32[BIN + 1] - BOUNDS->data.F32[BIN - 1]);            \
+            dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN - 1];                        \
+            Xo = 0.5 * (BOUNDS->data.F32[BIN + 1] + BOUNDS->data.F32[BIN]);                \
+            Yo = VECTOR->data.F32[BIN];                                                    \
+        }                                                                                  \
+        else                                                                               \
+        {                                                                                  \
+            dX = 0.5 * (BOUNDS->data.F32[BIN + 2] - BOUNDS->data.F32[BIN]);                \
+            dY = VECTOR->data.F32[BIN + 1] - VECTOR->data.F32[BIN];                        \
+            Xo = 0.5 * (BOUNDS->data.F32[BIN + 1] + BOUNDS->data.F32[BIN]);                \
+            Yo = VECTOR->data.F32[BIN];                                                    \
+        }                                                                                  \
+        if (dY != 0)                                                                       \
+        {                                                                                  \
+            Xt = (VALUE - Yo) * dX / dY + Xo;                                              \
+        }                                                                                  \
+        else                                                                               \
+        {                                                                                  \
+            Xt = Xo;                                                                       \
+        }                                                                                  \
+        Xt = PS_MIN(BOUNDS->data.F32[BIN + 1], PS_MAX(BOUNDS->data.F32[BIN], Xt));         \
         psTrace(TRACE, 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n", \
-                Xo, Yo, dX, dY, Xt, VALUE); \
-        RESULT = Xt; }
-
-# define COUNT_WARNING(LIMIT, INTERVAL, ...) { \
-        static int nCalls = 1; \
-        if (nCalls < LIMIT) { \
-            psWarning(__VA_ARGS__); \
-        } \
-        if (!(nCalls % INTERVAL)) { \
-            psWarning(__VA_ARGS__); \
+                Xo, Yo, dX, dY, Xt, VALUE);                                                \
+        RESULT = Xt;                                                                       \
+    }
+
+#define COUNT_WARNING(LIMIT, INTERVAL, ...)                 \
+    {                                                       \
+        static int nCalls = 1;                              \
+        if (nCalls < LIMIT)                                 \
+        {                                                   \
+            psWarning(__VA_ARGS__);                         \
+        }                                                   \
+        if (!(nCalls % INTERVAL))                           \
+        {                                                   \
+            psWarning(__VA_ARGS__);                         \
             psWarning("(warning raised %d times)", nCalls); \
-        } \
-        nCalls ++; \
-}
+        }                                                   \
+        nCalls++;                                           \
+    }
 
 // Debug information
@@ -195 +212 @@
 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,
-                                 const psVector* maskVector,
-                                 psVectorMaskType maskVal,
-                                 psStats* stats)
+static bool vectorSampleMean(const psVector *myVector,
+                             const psVector *errors,
+                             const psVector *maskVector,
+                             psVectorMaskType maskVal,
+                             psStats *stats)
 {
-    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 mean = 0.0; // The mean
     psF32 weight;
 
-    psF32 *data = myVector->data.F32;   // Dereference
-    int numData = myVector->n;          // Number of data points
+    psF32 *data = myVector->data.F32; // Dereference
+    int numData = myVector->n;        // Number of data points
 
     psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA;
     bool useRange = stats->options & PS_STAT_USE_RANGE;
 
-    psF32 sumWeights = 0.0;  // The sum of the weights
+    psF32 sumWeights = 0.0; // The sum of the weights
     psF32 *errorsData = (errors == NULL) ? NULL : errors->data.F32;
 
-    for (long i = 0; i < numData; i++) {
+    for (long i = 0; i < numData; i++)
+    {
         // Check if the data is with the specified range
         if (!isfinite(data[i]))
@@ -224 +242 @@
         if (maskData && (maskData[i] & maskVal))
             continue;
-        if (errors) {
+        if (errors)
+        {
             weight = (errorsData[i] == 0) ? 0.0 : PS_SQR(errorsData[i]);
             mean += data[i] * weight;
             sumWeights += weight;
-        } else {
+        }
+        else
+        {
             mean += data[i];
         }
         count++;
-
-    }
-    if (errors) {
+    }
+    if (errors)
+    {
         mean = (count > 0) ? mean / sumWeights : NAN;
-    } else {
+    }
+    else
+    {
         mean = (count > 0) ? mean / count : NAN;
     }
     stats->sampleMean = mean;
 
-    if (!isnan(mean)) {
+    if (!isnan(mean))
+    {
         stats->results |= PS_STAT_SAMPLE_MEAN;
     }
@@ -260 +284 @@
 (mask: 1, range: 0):  0.200 sec  0.244 sec
 *****************************************************************************/
-    static long vectorMinMax(const psVector* myVector,
-                             const psVector* maskVector,
-                             psVectorMaskType maskVal,
-                             psStats* stats
-        )
+static long vectorMinMax(const psVector *myVector,
+                         const psVector *maskVector,
+                         psVectorMaskType maskVal,
+                         psStats *stats)
 {
     psF32 max = -PS_MAX_F32;            // The calculated maximum
@@ -270 +293 @@
     psF32 *vector = myVector->data.F32; // Dereference the vector
 
-    int num = myVector->n;              // Number of values
-    int numValid = 0;                   // Number of valid values
+    int num = myVector->n; // Number of values
+    int numValid = 0;      // Number of valid values
 
     psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA;
     bool useRange = stats->options & PS_STAT_USE_RANGE;
 
-    for (long i = 0; i < num; i++) {
+    for (long i = 0; i < num; i++)
+    {
         // Check if the data is with the specified range
         if (!isfinite(vector[i]))
@@ -288 +312 @@
 
         numValid++;
-        max = PS_MAX (vector[i], max);
-        min = PS_MIN (vector[i], min);
+        max = PS_MAX(vector[i], max);
+        min = PS_MIN(vector[i], min);
     }
 
     // XXX save numValid in psStats?
-    if (numValid == 0) {
+    if (numValid == 0)
+    {
         stats->max = NAN;
         stats->min = NAN;
-    } else {
+    }
+    else
+    {
         stats->max = max;
         stats->min = min;
@@ -310 +337 @@
 were no valid input vector elements). Expects F32 vector for input.
 *****************************************************************************/
-static bool vectorSampleMedian(const psVector* inVector,
-                               const psVector* maskVector,
+static bool vectorSampleMedian(const psVector *inVector,
+                               const psVector *maskVector,
                                psVectorMaskType maskVal,
-                               psStats* stats)
+                               psStats *stats)
 {
     bool useRange = stats->options & PS_STAT_USE_RANGE;
     psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; // Dereference the vector
-    psF32 *input = inVector->data.F32; // Dereference the vector
+    psF32 *input = inVector->data.F32;                                                                    // Dereference the vector
 
     // use the temporary vector for the sorted output
-    stats->tmpData = psVectorRecycle (stats->tmpData, inVector->n, PS_TYPE_F32);
+    stats->tmpData = psVectorRecycle(stats->tmpData, inVector->n, PS_TYPE_F32);
     psVector *outVector = stats->tmpData;
     psF32 *output = outVector->data.F32; // Dereference the vector
 
-    if (maskData) psAssert (maskVector->n == inVector->n, "oops");
-
-    long count = 0;                     // Number of valid entries
+    if (maskData)
+        psAssert(maskVector->n == inVector->n, "oops");
+
+    long count = 0; // Number of valid entries
 
     // Store all non-masked data points within the min/max range
     // into the temporary vectors.
-    for (long i = 0; i < inVector->n; i++) {
-        psAssert (count >= 0, "oops");
-        psAssert (count < outVector->n, "oops");
-        psAssert (i >= 0, "oops");
-        psAssert (i < inVector->n, "oops");
-
-        if (!isfinite(input[i])) continue;
-        if (useRange && (input[i] < stats->min)) continue;
-        if (useRange && (input[i] > stats->max)) continue;
-        if (maskData && (maskData[i] & maskVal)) continue;
+    for (long i = 0; i < inVector->n; i++)
+    {
+        psAssert(count >= 0, "oops");
+        psAssert(count < outVector->n, "oops");
+        psAssert(i >= 0, "oops");
+        psAssert(i < inVector->n, "oops");
+
+        if (!isfinite(input[i]))
+            continue;
+        if (useRange && (input[i] < stats->min))
+            continue;
+        if (useRange && (input[i] > stats->max))
+            continue;
+        if (maskData && (maskData[i] & maskVal))
+            continue;
 
         output[count] = input[i];
@@ -346 +379 @@
     outVector->n = count;
 
-    if (count == 0) {
+    if (count == 0)
+    {
         COUNT_WARNING(10, 100, "No valid data in input vector.\n");
         stats->sampleUQ = NAN;
@@ -355 +389 @@
 
     // Sort the temporary vector.
-    if (!psVectorSort(outVector, outVector)) { // Sort in-place (since it's a copy, it's OK)
+    if (!psVectorSort(outVector, outVector))
+    { // Sort in-place (since it's a copy, it's OK)
         // an error in psVectorSort is a serious error:
-        // NULL input vector, psVectorCopy failure, invalid vector type
+        // NULL input vector, psVectorCopy failure, invalid vector type
         psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to sort input data."));
         stats->sampleUQ = NAN;
@@ -366 +401 @@
 
     // Calculate the median.  Use the average if the number of samples if even.
-    int midPt = (count/2);
-    psAssert (midPt >=           0, "oops");
-    psAssert (midPt < outVector->n, "oops");
-    if (count % 2 == 0) {
-        psAssert ((midPt - 1) >=           0, "oops");
-        psAssert ((midPt - 1) < outVector->n, "oops");
+    int midPt = (count / 2);
+    psAssert(midPt >= 0, "oops");
+    psAssert(midPt < outVector->n, "oops");
+    if (count % 2 == 0)
+    {
+        psAssert((midPt - 1) >= 0, "oops");
+        psAssert((midPt - 1) < outVector->n, "oops");
         stats->sampleMedian = 0.5 * (output[midPt - 1] + output[midPt]);
-    } else {
+    }
+    else
+    {
         stats->sampleMedian = output[midPt];
     }
 
-    int Qmin = (int)(0.25*count);
-    int Qmax = (int)(0.75*count);
-    psAssert (Qmin >=           0, "oops");
-    psAssert (Qmin < outVector->n, "oops");
-    psAssert (Qmax >=           0, "oops");
-    psAssert (Qmax < outVector->n, "oops");
+    int Qmin = (int)(0.25 * count);
+    int Qmax = (int)(0.75 * count);
+    psAssert(Qmin >= 0, "oops");
+    psAssert(Qmin < outVector->n, "oops");
+    psAssert(Qmax >= 0, "oops");
+    psAssert(Qmax < outVector->n, "oops");
 
     // Calculate the quartile points exactly.
     stats->sampleUQ = output[Qmax];
     stats->sampleLQ = output[Qmin];
-     
+
     stats->results |= PS_STAT_SAMPLE_MEDIAN;
     stats->results |= PS_STAT_SAMPLE_QUARTILE;
@@ -409 +447 @@
 *****************************************************************************/
 
-static bool vectorSampleStdev(const psVector* myVector,
-                              const psVector* errors,
-                              const psVector* maskVector,
+static bool vectorSampleStdev(const psVector *myVector,
+                              const psVector *errors,
+                              const psVector *maskVector,
                               psVectorMaskType maskVal,
-                              psStats* stats)
+                              psStats *stats)
 {
     // This procedure requires the mean.  If it has not been already
     // calculated, then call vectorSampleMean()
-    if (!(stats->results & PS_STAT_SAMPLE_MEAN)) {
+    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)) {
+    if (isnan(stats->sampleMean))
+    {
         COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN.");
         stats->sampleStdev = NAN;
@@ -428 +468 @@
     }
 
-    psF32 *data = myVector->data.F32;   // Dereference
+    psF32 *data = myVector->data.F32; // Dereference
     psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA;
     bool useRange = stats->options & PS_STAT_USE_RANGE;
@@ -434 +474 @@
 
     // Accumulate the sums
-    double mean = stats->sampleMean;    // The mean
-    double sumSquares = 0.0;            // Sum of the squares
-    double sumDiffs = 0.0;              // Sum of the differences
-    double errorDivisor = 0.0;          // Division for errors
-    long count = 0;                     // Number of data points being used
-    for (long i = 0; i < myVector->n; i++) {
+    double mean = stats->sampleMean; // The mean
+    double sumSquares = 0.0;         // Sum of the squares
+    double sumDiffs = 0.0;           // Sum of the differences
+    double errorDivisor = 0.0;       // Division for errors
+    long count = 0;                  // Number of data points being used
+    for (long i = 0; i < myVector->n; i++)
+    {
         // Check if the data is with the specified range
         if (!isfinite(data[i]))
             continue;
-        if (useRange && (data[i] < stats->min)) {
+        if (useRange && (data[i] < stats->min))
+        {
             continue;
         }
-        if (useRange && (data[i] > stats->max)) {
+        if (useRange && (data[i] > stats->max))
+        {
             continue;
         }
-        if (maskData && (maskData[i] & maskVal)) {
+        if (maskData && (maskData[i] & maskVal))
+        {
             continue;
         }
@@ -457 +501 @@
         sumDiffs += diff;
         count++;
-        if (errors) {
+        if (errors)
+        {
             errorDivisor += 1.0 / PS_SQR(errorsData[i]);
         }
     }
 
-    if (count == 0) {
+    if (count == 0)
+    {
         // This is an ambiguous case: error or not?
         // It's not an empty array (that's been asserted on previously), but everything's been masked out.
@@ -470 +516 @@
         return true;
     }
-    if (count == 1) {
+    if (count == 1)
+    {
         stats->sampleStdev = 0.0;
         COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): only one valid psVector elements (%ld). Setting stats->sampleStdev = 0.0.\n", count);
@@ -476 +523 @@
     }
 
-    if (errors) {
+    if (errors)
+    {
         stats->sampleStdev = (1.0 / sqrtf(errorDivisor));
-    } else {
+    }
+    else
+    {
         stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) / (float)(count - 1));
     }
@@ -486 +536 @@
 }
 
-static bool vectorSampleMoments(const psVector* myVector,
-                                const psVector* maskVector,
+static bool vectorSampleMoments(const psVector *myVector,
+                                const psVector *maskVector,
                                 psVectorMaskType maskVal,
-                                psStats* stats)
+                                psStats *stats)
 {
     // This procedure requires the mean and standard deviation
-    if (!(stats->results & PS_STAT_SAMPLE_MEAN)) {
+    if (!(stats->results & PS_STAT_SAMPLE_MEAN))
+    {
         vectorSampleMean(myVector, NULL, maskVector, maskVal, stats);
     }
-    if (isnan(stats->sampleMean)) {
+    if (isnan(stats->sampleMean))
+    {
         COUNT_WARNING(10, 100, "WARNING: vectorSampleMoments(): sample mean is NAN.\n");
         goto SAMPLE_MOMENTS_BAD;
     }
-    if (!(stats->results & PS_STAT_SAMPLE_STDEV)) {
+    if (!(stats->results & PS_STAT_SAMPLE_STDEV))
+    {
         vectorSampleStdev(myVector, NULL, maskVector, maskVal, stats);
     }
-    if (isnan(stats->sampleStdev) || stats->sampleStdev == 0.0) {
+    if (isnan(stats->sampleStdev) || stats->sampleStdev == 0.0)
+    {
         COUNT_WARNING(10, 100, "WARNING: vectorSampleMoments(): sample stdev is NAN or 0.\n");
         goto SAMPLE_MOMENTS_BAD;
     }
 
-    psF32 *data = myVector->data.F32;   // Dereference
+    psF32 *data = myVector->data.F32; // Dereference
     psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA;
     bool useRange = stats->options & PS_STAT_USE_RANGE;
 
     // Accumulate the sums
-    double mean = stats->sampleMean;    // The mean
-    double sum3 = 0.0;                  // Sum of the cubes of the differences
-    double sum4 = 0.0;                  // Sum of the fourth powers of the differences
-    long count = 0;                     // Number of data points being used
-    for (long i = 0; i < myVector->n; i++) {
+    double mean = stats->sampleMean; // The mean
+    double sum3 = 0.0;               // Sum of the cubes of the differences
+    double sum4 = 0.0;               // Sum of the fourth powers of the differences
+    long count = 0;                  // Number of data points being used
+    for (long i = 0; i < myVector->n; i++)
+    {
         // Check if the data is with the specified range
         if (!isfinite(data[i]))
             continue;
-        if (useRange && (data[i] < stats->min)) {
+        if (useRange && (data[i] < stats->min))
+        {
             continue;
         }
-        if (useRange && (data[i] > stats->max)) {
+        if (useRange && (data[i] > stats->max))
+        {
             continue;
         }
-        if (maskData && (maskData[i] & maskVal)) {
+        if (maskData && (maskData[i] & maskVal))
+        {
             continue;
         }
 
-        double diff = data[i] - mean;   // Difference from the mean
-        double temp;                    // Temporary variable for accumulating
+        double diff = data[i] - mean; // Difference from the mean
+        double temp;                  // Temporary variable for accumulating
 
         sum3 += temp = PS_SQR(diff);
@@ -539 +597 @@
     }
 
-    psAssert(count > 1, "impossible");                  // It should be, because we have a mean and standard deviation
-
-    double stdev = stats->sampleStdev;  // Standard deviation
-    double variance = PS_SQR(stdev);    // Variance
+    psAssert(count > 1, "impossible"); // It should be, because we have a mean and standard deviation
+
+    double stdev = stats->sampleStdev; // Standard deviation
+    double variance = PS_SQR(stdev);   // Variance
 
     // Formula for skewness and kurtosis from Numerical Recipes in C, p 613.
@@ -552 +610 @@
     return true;
 
- SAMPLE_MOMENTS_BAD:
+SAMPLE_MOMENTS_BAD:
     // stats->sampleStdev has already been set
     stats->sampleSkewness = NAN;
@@ -571 +629 @@
     true for success; false otherwise
 *****************************************************************************/
-static bool vectorClippedStats(const psVector* myVector,
-                               const psVector* errors,
-                               psVector* maskInput,
+static bool vectorClippedStats(const psVector *myVector,
+                               const psVector *errors,
+                               psVector *maskInput,
                                psVectorMaskType maskValInput,
-                               psStats* stats
-    )
+                               psStats *stats)
 {
     // Ensure that stats->clipIter is within the proper range.
@@ -590 +647 @@
     // unless we succeed, these will have NAN values
     stats->clippedMean = NAN;
+    stats->clippedMedian = NAN;
     stats->clippedStdev = NAN;
     stats->clippedNvalues = 0;
@@ -597 +655 @@
 
     // use the temporary vector for local temporary mask
-    stats->tmpMask = psVectorRecycle (stats->tmpMask, myVector->n, PS_TYPE_VECTOR_MASK);
+    stats->tmpMask = psVectorRecycle(stats->tmpMask, myVector->n, PS_TYPE_VECTOR_MASK);
     psVector *tmpMask = stats->tmpMask;
     psVectorInit(tmpMask, 0);
-    if (maskInput) {
-        for (long i = 0; i < myVector->n; i++) {
-            if (maskInput->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValInput) {
+    if (maskInput)
+    {
+        for (long i = 0; i < myVector->n; i++)
+        {
+            if (maskInput->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValInput)
+            {
                 tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal;
             }
@@ -610 +671 @@
     // 1. Compute the sample median, which we save for output
     vectorSampleMedian(myVector, tmpMask, maskVal, stats);
-    if (isnan(stats->sampleMedian)) {
+    if (isnan(stats->sampleMedian))
+    {
         COUNT_WARNING(10, 100, "Call to vectorSampleMedian returned NAN\n");
         return true;
@@ -618 +680 @@
     // 2. Compute the sample standard deviation, which we save for output
     vectorSampleStdev(myVector, errors, tmpMask, maskVal, stats);
-    if (isnan(stats->sampleStdev)) {
+    if (isnan(stats->sampleStdev))
+    {
         COUNT_WARNING(10, 100, "Call to vectorSampleStdev returned NAN\n");
         return true;
@@ -625 +688 @@
 
     // 3. Use the sample median as the first estimator of the mean X.
-    psF32 clippedMean = stats->sampleMedian;
-
+    psF32 clippedMedian = stats->sampleMedian;
+    
     // 4. Use the sample stdev as the first estimator of the mean stdev.
     psF32 clippedStdev = stats->sampleStdev;
 
     // 5. Repeat N (stats->clipIter) times:
-    long numClipped = 0;                // Number of values clipped
-    bool clipped = true;                // Have we clipped anything in this iteration
-    for (int iter = 0; iter < stats->clipIter && clipped; iter++) {
+    long numClipped = 0; // Number of values clipped
+    bool clipped = true; // Have we clipped anything in this iteration
+    for (int iter = 0; iter < stats->clipIter && clipped; iter++)
+    {
         clipped = false;
         psTrace(TRACE, 6, "------------ Iteration %d ------------\n", iter);
         // a) Exclude all values x_i for which |x_i - x| > K * stdev
-        if (errors) {
+        if (errors)
+        {
             // XXXX if we convert errors to variance, this should square the other terms (A*A faster than
             // sqrt(A))
-            for (long j = 0; j < myVector->n; j++) {
+            for (long j = 0; j < myVector->n; j++)
+            {
                 if (!tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] &&
-                    fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) {
+                    fabsf(myVector->data.F32[j] - clippedMedian) > stats->clipSigma * errors->data.F32[j])
+                {
                     tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff;
                     psTrace(TRACE, 10, "Clipped %ld: %f +/- %f\n", j,
@@ -650 +717 @@
                 }
             }
-        } else {
-            for (long j = 0; j < myVector->n; j++) {
+        }
+        else
+        {
+            for (long j = 0; j < myVector->n; j++)
+            {
                 if (!tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] &&
-                    fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) {
+                    fabsf(myVector->data.F32[j] - clippedMedian) > (stats->clipSigma * clippedStdev))
+                {
                     tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff;
                     psTrace(TRACE, 10, "Clipped %ld: %f\n", j, myVector->data.F32[j]);
@@ -662 +733 @@
         }
 
-        // b) compute new mean and stdev
+        // b) compute new median and stdev
         // Allocate a psStats structure for calculating the mean and stdev.
-        // XXX Can we just use this psStats structure to calculate the SAMPLE MEAN and STDEV?
-        // psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
-        // vectorSampleMean(myVector, errors, tmpMask, maskVal, statsTmp);
+        // XXX Can we just use this psStats structure to calculate the SAMPLE MEDIAN and STDEV?
+        // psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+        // vectorSampleMedian(myVector, errors, tmpMask, maskVal, statsTmp);
         // vectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp);
-        vectorSampleMean(myVector, errors, tmpMask, maskVal, stats);
+        vectorSampleMedian(myVector, tmpMask, maskVal, stats);
         vectorSampleStdev(myVector, errors, tmpMask, maskVal, stats);
-        psTrace(TRACE, 6, "The new sample mean is %f\n", stats->sampleMean);
+        psTrace(TRACE, 6, "The new sample mean is %f\n", stats->sampleMedian);
         psTrace(TRACE, 6, "The new sample stdev is %f\n", stats->sampleStdev);
 
-        // If the new mean and stdev are NAN, we must exit the loop.
+        // If the new median and stdev are NAN, we must exit the loop.
         // Otherwise, use the new results and continue.
-        if (isnan(stats->sampleMean) || isnan(stats->sampleStdev)) {
+        if (isnan(stats->sampleMedian) || isnan(stats->sampleStdev))
+        {
             iter = stats->clipIter;
-            COUNT_WARNING(10, 100, "vectorSampleMean() or vectorSampleStdev() returned a NAN.\n");
-            clippedMean = NAN;
+            COUNT_WARNING(10, 100, "vectorSampleMedian() or vectorSampleStdev() returned a NAN.\n");
+            clippedMedian = NAN;
             clippedStdev = NAN;
             return true;
-        } else {
-            clippedMean = stats->sampleMean;
+        }
+        else
+        {
+            clippedMedian = stats->sampleMedian;
             clippedStdev = stats->sampleStdev;
         }
@@ -690 +764 @@
     // Number of values used in calculation is the total number of data values, minus those we clipped
     stats->clippedNvalues = myVector->n - numClipped;
-
-    // 7. The last calcuated value of x is the clipped mean.
+    // calculate the clipped mean
+    psF32 clippedMean = stats->sampleMean;
+    vectorSampleMean(myVector, errors, tmpMask, maskVal, stats);
+    if (isnan(stats->sampleMean))
+    {
+        COUNT_WARNING(10, 100, "vectorSampleMean() returned a NAN.\n");
+        clippedMean = NAN;
+        return true;
+    }
+    else
+    {
+        clippedMean = stats->sampleMean;
+    }
+
+    // 7. The last calcuated value of x is the clipped median.
     // 8. The last calcuated value of stdev is the clipped stdev.
-    // we always return both stats even if only one was requested
+    // we always return all stats even if only one was requested
     stats->clippedMean = clippedMean;
+    stats->clippedMedian = clippedMedian;
     stats->clippedStdev = clippedStdev;
 
     stats->results |= PS_STAT_CLIPPED_MEAN;
+    stats->results |= PS_STAT_CLIPPED_MEDIAN;
     stats->results |= PS_STAT_CLIPPED_STDEV;
 
     psTrace(TRACE, 6, "The final clipped mean is %f\n", clippedMean);
+    psTrace(TRACE, 6, "The final clipped median is %f\n", clippedMedian);
     psTrace(TRACE, 6, "The final clipped stdev is %f\n", clippedStdev);
 
@@ -725 +815 @@
 *****************************************************************************/
 #define INITIAL_NUM_BINS 1000.0
-static bool vectorRobustStats(const psVector* myVector,
-                              const psVector* errors,
-                              psVector* maskInput,
+static bool vectorRobustStats(const psVector *myVector,
+                              const psVector *errors,
+                              psVector *maskInput,
                               psVectorMaskType maskValInput,
-                              psStats* stats)
+                              psStats *stats)
 {
-    if (psTraceGetLevel("psLib.math") >= 8) {
+    if (psTraceGetLevel("psLib.math") >= 8)
+    {
         PS_VECTOR_PRINT_F32(myVector);
     }
@@ -741 +832 @@
     psVector *mask = psVectorAlloc(myVector->n, PS_TYPE_VECTOR_MASK); // The actual mask we will use
     psVectorInit(mask, 0);
-    if (maskInput) {
-        for (long i = 0; i < myVector->n; i++) {
-            if (maskInput->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValInput) {
+    if (maskInput)
+    {
+        for (long i = 0; i < myVector->n; i++)
+        {
+            if (maskInput->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValInput)
+            {
                 mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal;
             }
@@ -751 +845 @@
     // 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
-    psHistogram *cumulative = NULL;     // Cumulative histogram of the data
-    float min = NAN, max = NAN;         // Mimimum and maximum values
-    float sigma = NAN;                  // The robust standard deviation
-    long totalDataPoints = 0;           // Total number of (unmasked) data points
-
-    float binSize = 0.0;            // Size of bins for histogram
+    psHistogram *histogram = NULL;                                  // Histogram of the data
+    psHistogram *cumulative = NULL;                                 // Cumulative histogram of the data
+    float min = NAN, max = NAN;                                     // Mimimum and maximum values
+    float sigma = NAN;                                              // The robust standard deviation
+    long totalDataPoints = 0;                                       // Total number of (unmasked) data points
+
+    float binSize = 0.0; // Size of bins for histogram
     long binLo, binHi;
     long binL2, binH2;
@@ -764 +858 @@
 
     // Iterate to get the best bin size; an iteration limit is enforced at the bottom of the loop.
-    for (int iterate = 1; iterate > 0; iterate++) {
+    for (int iterate = 1; iterate > 0; iterate++)
+    {
         psTrace(TRACE, 6, "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n", iterate);
 
-        if (iterate >= PS_ROBUST_MAX_ITERATIONS) {
-          // This occurs when a large number of the values are identical --- a bin size cannot be found
-          // that will spread out the distribution.  Therefore, set what we can, and fall over
-          // gracefully.
-          COUNT_WARNING(10, 100, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS);
-          goto escape;
+        if (iterate >= PS_ROBUST_MAX_ITERATIONS)
+        {
+            // This occurs when a large number of the values are identical --- a bin size cannot be found
+            // that will spread out the distribution.  Therefore, set what we can, and fall over
+            // gracefully.
+            COUNT_WARNING(10, 100, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS);
+            goto escape;
         }
 
@@ -779 +875 @@
         min = statsMinMax->min;
         max = statsMinMax->max;
-        if (numValid == 0 || isnan(min) || isnan(max)) {
+        if (numValid == 0 || isnan(min) || isnan(max))
+        {
             // Data range calculation failed
             COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
             goto escape;
         }
-        if (!isfinite(max - min)) {
-            COUNT_WARNING(10, 100, "Range of of the input vector is too large: %lf.\n", (double)max - (double) min);
+        if (!isfinite(max - min))
+        {
+            COUNT_WARNING(10, 100, "Range of of the input vector is too large: %lf.\n", (double)max - (double)min);
             goto escape;
         }
@@ -791 +889 @@
 
         // If all data points have the same value, then we set the appropriate members of stats and return.
-        if (fabs(max - min) <= FLT_EPSILON) {
+        if (fabs(max - min) <= FLT_EPSILON)
+        {
             stats->robustMedian = min;
             stats->robustStdev = 0.0;
@@ -807 +906 @@
         }
 
-        if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE)) {
+        if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE))
+        {
             // Set initial bin size to the specified value.
             binSize = stats->binsize;
             psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize);
-        } else {
+        }
+        else
+        {
             // Determine the bin size of the robust histogram, using the pre-defined number of bins
-            binSize = (max - min) / INITIAL_NUM_BINS;
+            binSize = (max - min) / INITIAL_NUM_BINS;
         }
         psTrace(TRACE, 6, "Initial robust bin size is %.2f\n", binSize);
@@ -830 +932 @@
         // Assert here so we can get more information about what is going wrong.
         psAssert(numBins > 0, "Invalid numBins: %ld max: %f min: %f binSize: %f", numBins, max, min, binSize);
-        
+
         // Generate the histogram
-        histogram = psHistogramAlloc(min - 2.0*binSize, max + 2.0*binSize, numBins);
+        histogram = psHistogramAlloc(min - 2.0 * binSize, max + 2.0 * binSize, numBins);
         // XXXXX we need to consider this step if errors -> variance
-        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
+        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal))
+        {
             // if psVectorHistogram returns false, we have a programming error
             psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for robust statistics.\n");
@@ -843 +946 @@
             return false;
         }
-        if (psTraceGetLevel("psLib.math") >= 8) {
+        if (psTraceGetLevel("psLib.math") >= 8)
+        {
             PS_VECTOR_PRINT_F32(histogram->bounds);
             PS_VECTOR_PRINT_F32(histogram->nums);
@@ -854 +958 @@
         int nMaxBin = histogram->nums->data.F32[0];
         int iMaxBin = 0;
-        for (long i = 1; i < histogram->nums->n; i++) {
-            if (histogram->nums->data.F32[i] > nMaxBin) {
+        for (long i = 1; i < histogram->nums->n; i++)
+        {
+            if (histogram->nums->data.F32[i] > nMaxBin)
+            {
                 nMaxBin = histogram->nums->data.F32[i];
                 iMaxBin = i;
             }
         }
-        if (nMaxBin > numValid / 2) {
-            float minKeep = histogram->bounds->data.F32[iMaxBin] - 10*binSize;
-            float maxKeep = histogram->bounds->data.F32[iMaxBin + 1] + 10*binSize;
+        if (nMaxBin > numValid / 2)
+        {
+            float minKeep = histogram->bounds->data.F32[iMaxBin] - 10 * binSize;
+            float maxKeep = histogram->bounds->data.F32[iMaxBin + 1] + 10 * binSize;
             int nInvalid = 0;
-            for (long i = 0; i < myVector->n; i++) {
+            for (long i = 0; i < myVector->n; i++)
+            {
                 // skip the already-masked values
-                if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskVal) continue;
+                if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskVal)
+                    continue;
                 bool invalid = false;
                 invalid |= (myVector->data.F32[i] < minKeep);
                 invalid |= (myVector->data.F32[i] > maxKeep);
                 invalid |= (!isfinite(myVector->data.F32[i]));
-                if (!invalid) continue;
+                if (!invalid)
+                    continue;
                 mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal;
-                nInvalid ++;
-            }
-
-            if (nInvalid) {
-              psTrace(TRACE, 6, "data is concentrated in a single bin (%d = %f - %f), masking %d extreme outliers and retrying\n", 
-                      iMaxBin, histogram->bounds->data.F32[iMaxBin], histogram->bounds->data.F32[iMaxBin+1], nInvalid);
-              psFree(histogram);
-              psFree(cumulative);
-              histogram = NULL;
-              cumulative = NULL;
-              continue;
+                nInvalid++;
+            }
+
+            if (nInvalid)
+            {
+                psTrace(TRACE, 6, "data is concentrated in a single bin (%d = %f - %f), masking %d extreme outliers and retrying\n",
+                        iMaxBin, histogram->bounds->data.F32[iMaxBin], histogram->bounds->data.F32[iMaxBin + 1], nInvalid);
+                psFree(histogram);
+                psFree(cumulative);
+                histogram = NULL;
+                cumulative = NULL;
+                continue;
             }
             // if we did not mask anything, give up.
         }
 
-        // We were causing psHistogramAlloc to assert. 
+        // We were causing psHistogramAlloc to assert.
         // Assert here so we can get more information about what is going wrong.
         psAssert(numBins > 0, "Invalid numBins %ld max: %f min: %f binSize: %f", numBins, max, min, binSize);
@@ -898 +1009 @@
         cumulative->bounds->data.F32[0] = histogram->bounds->data.F32[1];
 
-        // Correctly fill the cumulative distribution with monotonically increasing values (skip zero valued bins).
-        long Nc = 1;  // track the current bin of cumulative
-        // the boundaries for the current cumulative bin are from upper end of the last valid histogram bin to the 
-        // upper end of the current histogram bin
-        for (long i = 1; i < histogram->nums->n - 1; i++) {
-            if (histogram->nums->data.F32[i] == 0.0) continue;
-            cumulative->nums->data.F32[Nc] = cumulative->nums->data.F32[Nc - 1] + histogram->nums->data.F32[i];
-            cumulative->bounds->data.F32[Nc] = histogram->bounds->data.F32[i+1];
-            Nc ++;
-        }
-        long Nlast = Nc - 1;  // last valid cumulative bin 
-        for (long i = Nc; i < histogram->nums->n; i++) { // Ensure the unused entries are filled.
-            cumulative->nums->data.F32[i] = cumulative->nums->data.F32[Nlast];
-            cumulative->bounds->data.F32[i] = cumulative->bounds->data.F32[i-1] + 1.0;
-        }
-        
-        if (psTraceGetLevel("psLib.math") >= 8) {
+        // Correctly fill the cumulative distribution with monotonically increasing values (skip zero valued bins).
+        long Nc = 1; // track the current bin of cumulative
+        // the boundaries for the current cumulative bin are from upper end of the last valid histogram bin to the
+        // upper end of the current histogram bin
+        for (long i = 1; i < histogram->nums->n - 1; i++)
+        {
+            if (histogram->nums->data.F32[i] == 0.0)
+                continue;
+            cumulative->nums->data.F32[Nc] = cumulative->nums->data.F32[Nc - 1] + histogram->nums->data.F32[i];
+            cumulative->bounds->data.F32[Nc] = histogram->bounds->data.F32[i + 1];
+            Nc++;
+        }
+        long Nlast = Nc - 1; // last valid cumulative bin
+        for (long i = Nc; i < histogram->nums->n; i++)
+        { // Ensure the unused entries are filled.
+            cumulative->nums->data.F32[i] = cumulative->nums->data.F32[Nlast];
+            cumulative->bounds->data.F32[i] = cumulative->bounds->data.F32[i - 1] + 1.0;
+        }
+
+        if (psTraceGetLevel("psLib.math") >= 8)
+        {
             PS_VECTOR_PRINT_F32(cumulative->bounds);
             PS_VECTOR_PRINT_F32(cumulative->nums);
@@ -924 +1039 @@
 
         // find bin which is the lower bound of median (value[binMedian] < median < value[binMedian+1]
-        PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints/2.0, 0);
-
-        psTrace(TRACE, 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian, cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]);
+        PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints / 2.0, 0);
+
+        psTrace(TRACE, 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian, cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian + 1]);
 
         // ADD step 3: Interpolate to the exact 50% position in bin units
         // stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
-        // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
+        // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
         // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian);
-        // There's no reason to do a quadratic fit near the 50% bin, as it's approximately linear there.
-        // Instead, do a 5-point linear fit.
-        stats->robustMedian = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
+        // There's no reason to do a quadratic fit near the 50% bin, as it's approximately linear there.
+        // Instead, do a 5-point linear fit.
+        stats->robustMedian = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints / 2.0);
 
         // convert bin to bin value: this is the robust histogram median.
-        if (isnan(stats->robustMedian)) {
+        if (isnan(stats->robustMedian))
+        {
             COUNT_WARNING(10, 100, "Failed to fit a quadratic and calculate the 50-percent position.\n");
             goto escape;
@@ -949 +1065 @@
         PS_BIN_FOR_VALUE(binL2, cumulative->nums, totalDataPoints * 0.308538f, 0);
         PS_BIN_FOR_VALUE(binH2, cumulative->nums, totalDataPoints * 0.691462f, 0);
-        PS_BIN_FOR_VALUE(binL4, cumulative->nums, totalDataPoints * 0.022750f, 0);
-        PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977250f, 0);
-        
-        
+        PS_BIN_FOR_VALUE(binL4, cumulative->nums, totalDataPoints * 0.022750f, 0);
+        PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977250f, 0);
+
         psTrace(TRACE, 6, "The 15.8655%% and 84.1345%% data point bins are (%ld, %ld).\n",
                 binLo, binHi);
@@ -960 +1075 @@
         psTrace(TRACE, 6, "binH2 midpoint is %f\n", PS_BIN_MIDPOINT(cumulative, binH2));
 
-        if ((binLo < 0) || (binHi < 0)) {
+        if ((binLo < 0) || (binHi < 0))
+        {
             COUNT_WARNING(10, 100, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n");
             goto escape;
         }
-    
+
         // ADD step 4b: Interpolate Sigma (linearly) to find these two positions exactly: these are the 1sigma
         // positions.
         psTrace(TRACE, 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]);
+                binLo, cumulative->nums->data.F32[binLo], cumulative->nums->data.F32[binLo + 1]);
         psTrace(TRACE, 6, "binHi is %ld.  Nums at that bin and the next are (%.2f, %.2f)\n",
-                binHi, cumulative->nums->data.F32[binHi], cumulative->nums->data.F32[binHi+1]);
+                binHi, cumulative->nums->data.F32[binHi], cumulative->nums->data.F32[binHi + 1]);
 
         // find the +0.5 and -0.5 sigma points with linear interpolation.  binLo and binHi are the bins
@@ -977 +1093 @@
         float binLoF32, binHiF32, binL2F32, binH2F32, binL4F32, binH4F32;
 #if (0)
-        PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binLo,
-                            totalDataPoints * 0.158655f);
-        PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binHi,
-                            totalDataPoints * 0.841345f);
-        PS_BIN_INTERPOLATE (binL2F32, cumulative->nums, cumulative->bounds, binL2,
-                            totalDataPoints * 0.308538f);
-        PS_BIN_INTERPOLATE (binH2F32, cumulative->nums, cumulative->bounds, binH2,
-                            totalDataPoints * 0.691462f);
-        PS_BIN_INTERPOLATE (binL4F32, cumulative->nums, cumulative->bounds, binL4,
-                            totalDataPoints * 0.022750f);
-        PS_BIN_INTERPOLATE (binH4F32, cumulative->nums, cumulative->bounds, binH4,
-                            totalDataPoints * 0.977250f);
+        PS_BIN_INTERPOLATE(binLoF32, cumulative->nums, cumulative->bounds, binLo,
+                           totalDataPoints * 0.158655f);
+        PS_BIN_INTERPOLATE(binHiF32, cumulative->nums, cumulative->bounds, binHi,
+                           totalDataPoints * 0.841345f);
+        PS_BIN_INTERPOLATE(binL2F32, cumulative->nums, cumulative->bounds, binL2,
+                           totalDataPoints * 0.308538f);
+        PS_BIN_INTERPOLATE(binH2F32, cumulative->nums, cumulative->bounds, binH2,
+                           totalDataPoints * 0.691462f);
+        PS_BIN_INTERPOLATE(binL4F32, cumulative->nums, cumulative->bounds, binL4,
+                           totalDataPoints * 0.022750f);
+        PS_BIN_INTERPOLATE(binH4F32, cumulative->nums, cumulative->bounds, binH4,
+                           totalDataPoints * 0.977250f);
 #else
         binLoF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo, totalDataPoints * 0.158655);
-        binHiF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi, totalDataPoints * 0.841345);          
+        binHiF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi, totalDataPoints * 0.841345);
         binL2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binL2, totalDataPoints * 0.308538);
-        binH2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH2, totalDataPoints * 0.691462);          
+        binH2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH2, totalDataPoints * 0.691462);
         binL4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binL4, totalDataPoints * 0.022750);
-        binH4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH4, totalDataPoints * 0.977250);
-#endif  
+        binH4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH4, totalDataPoints * 0.977250);
+#endif
         // report +/- 1 sigma points
         psTrace(TRACE, 5,
@@ -1008 +1124 @@
                 binL4F32, binH4F32);
 
-        // If some of the fits failed, attempt to fix this
-        if (!isfinite(binLoF32) && isfinite(binHiF32)) { binLoF32 = -1.0 * binHiF32; }
-        if (!isfinite(binHiF32) && isfinite(binLoF32)) { binHiF32 = -1.0 * binLoF32; }
-        if (!isfinite(binL2F32) && isfinite(binH2F32)) { binL2F32 = -1.0 * binH2F32; }
-        if (!isfinite(binH2F32) && isfinite(binL2F32)) { binH2F32 = -1.0 * binL2F32; }
-        if (!isfinite(binL4F32) && isfinite(binH4F32)) { binL4F32 = -1.0 * binH4F32; }
-        if (!isfinite(binH4F32) && isfinite(binL4F32)) { binH4F32 = -1.0 * binL4F32; }
-        
+        // If some of the fits failed, attempt to fix this
+        if (!isfinite(binLoF32) && isfinite(binHiF32))
+        {
+            binLoF32 = -1.0 * binHiF32;
+        }
+        if (!isfinite(binHiF32) && isfinite(binLoF32))
+        {
+            binHiF32 = -1.0 * binLoF32;
+        }
+        if (!isfinite(binL2F32) && isfinite(binH2F32))
+        {
+            binL2F32 = -1.0 * binH2F32;
+        }
+        if (!isfinite(binH2F32) && isfinite(binL2F32))
+        {
+            binH2F32 = -1.0 * binL2F32;
+        }
+        if (!isfinite(binL4F32) && isfinite(binH4F32))
+        {
+            binL4F32 = -1.0 * binH4F32;
+        }
+        if (!isfinite(binH4F32) && isfinite(binL4F32))
+        {
+            binH4F32 = -1.0 * binL4F32;
+        }
+
         // ADD step 5: Determine SIGMA as the distance between binL2 and binH2 (+/- 0.5 sigma)
-
 
         float sigma1 = (binH2F32 - binL2F32);
@@ -1023 +1156 @@
         float sigma4 = (binH4F32 - binL4F32) / 4.0;
 
-        // Fix again?
-        if (!isfinite(sigma1) && isfinite(sigma2) && isfinite(sigma4)) { sigma1 = (sigma2 + sigma4) / 2.0; }
-        if (!isfinite(sigma2) && isfinite(sigma1) && isfinite(sigma4)) { sigma2 = (sigma1 + sigma4) / 2.0; }
-        if (!isfinite(sigma4) && isfinite(sigma2) && isfinite(sigma1)) { sigma4 = (sigma2 + sigma1) / 2.0; }
-        
+        // Fix again?
+        if (!isfinite(sigma1) && isfinite(sigma2) && isfinite(sigma4))
+        {
+            sigma1 = (sigma2 + sigma4) / 2.0;
+        }
+        if (!isfinite(sigma2) && isfinite(sigma1) && isfinite(sigma4))
+        {
+            sigma2 = (sigma1 + sigma4) / 2.0;
+        }
+        if (!isfinite(sigma4) && isfinite(sigma2) && isfinite(sigma1))
+        {
+            sigma4 = (sigma2 + sigma1) / 2.0;
+        }
+
         // take the smallest of the three: if we have a clump with wide outliers, sigma2 and
         // sigma4 will be biased high; if we have a bi-modal distribution, sigma1 and sigma2
         // will be biased high.
-        //        sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4));
-        // CZW: Instead, take the median.  Taking the MIN forces a bias on unbiased data.
-        //      It seems like occasionally getting the wrong answer on a complex distribution
-        //      is more acceptable than always getting the wrong answer for simple ones.
-
-        
-        sigma = PS_MAX( PS_MIN(sigma1,sigma2),
-                        PS_MIN( PS_MAX(sigma1,sigma2),
-                                sigma4));
+        //        sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4));
+        // CZW: Instead, take the median.  Taking the MIN forces a bias on unbiased data.
+        //      It seems like occasionally getting the wrong answer on a complex distribution
+        //      is more acceptable than always getting the wrong answer for simple ones.
+
+        sigma = PS_MAX(PS_MIN(sigma1, sigma2),
+                       PS_MIN(PS_MAX(sigma1, sigma2),
+                              sigma4));
 
         psTrace(TRACE, 6, "The 1x sigma is %f.\n", sigma1);
@@ -1046 +1187 @@
 
         psTrace(TRACE, 6, "The current sigma is %f.\n", sigma);
-        //        stats->robustStdev = sigma;
-        stats->robustStdev = sigma;
-
-#if (CZW && 0) 
-        // Skewness check: Find least biased sample for each pair.
-        sigma1 = 2.0 * PS_MIN(binH2F32 - stats->robustMedian,
-                              stats->robustMedian - binL2F32);
-        sigma2 = 1.0 * PS_MIN(binHiF32 - stats->robustMedian,
-                              stats->robustMedian - binLoF32);
-        sigma4 = 0.5 * PS_MIN(binH4F32 - stats->robustMedian,
-                              stats->robustMedian - binL4F32);
-        // Kurtosis check: Take median sample as the solution.
-        stats->robustStdev = PS_MAX( PS_MIN(sigma1,sigma2),
-                                     PS_MIN( PS_MAX(sigma1,sigma2),
-                                             sigma4));
+        //        stats->robustStdev = sigma;
+        stats->robustStdev = sigma;
+
+#if (CZW && 0)
+        // Skewness check: Find least biased sample for each pair.
+        sigma1 = 2.0 * PS_MIN(binH2F32 - stats->robustMedian,
+                              stats->robustMedian - binL2F32);
+        sigma2 = 1.0 * PS_MIN(binHiF32 - stats->robustMedian,
+                              stats->robustMedian - binLoF32);
+        sigma4 = 0.5 * PS_MIN(binH4F32 - stats->robustMedian,
+                              stats->robustMedian - binL4F32);
+        // Kurtosis check: Take median sample as the solution.
+        stats->robustStdev = PS_MAX(PS_MIN(sigma1, sigma2),
+                                    PS_MIN(PS_MAX(sigma1, sigma2),
+                                           sigma4));
 #endif
 
-        
 #if (CZW)
-        //      printf("CZW: bad sigma?: %f %f  %f %f  %f %f  %f %f %f  %f\n",
-        //             binH2F32,binL2F32,binHiF32,binLoF32,binH4F32,binL4F32,
-        //             sigma1,sigma2,sigma4,sigma);
-        
-        printf("CZW Robust (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n",
-               iterate,
-               stats->robustMedian,stats->robustStdev,
-               fabs(cumulative->bounds->data.F32[binMedian] - cumulative->bounds->data.F32[binMedian + 1]),
-               
-               cumulative->bounds->data.F32[binMedian-3],cumulative->bounds->data.F32[binMedian-2],
-               cumulative->bounds->data.F32[binMedian-1],
-               cumulative->bounds->data.F32[binMedian],
-               cumulative->bounds->data.F32[binMedian+1],
-               cumulative->bounds->data.F32[binMedian+2],cumulative->bounds->data.F32[binMedian+3],
-               
-               cumulative->nums->data.F32[binMedian-3],cumulative->nums->data.F32[binMedian-2],
-               cumulative->nums->data.F32[binMedian-1],
-               cumulative->nums->data.F32[binMedian],
-               cumulative->nums->data.F32[binMedian+1],
-               cumulative->nums->data.F32[binMedian+2],cumulative->nums->data.F32[binMedian+3]);
-        //      PS_VECTOR_PRINT_F32(histogram->bounds);
-        //      PS_VECTOR_PRINT_F32(histogram->nums);
-        //      PS_VECTOR_PRINT_F32(cumulative->bounds);
-        //      PS_VECTOR_PRINT_F32(cumulative->nums);
+        //      printf("CZW: bad sigma?: %f %f  %f %f  %f %f  %f %f %f  %f\n",
+        //             binH2F32,binL2F32,binHiF32,binLoF32,binH4F32,binL4F32,
+        //             sigma1,sigma2,sigma4,sigma);
+
+        printf("CZW Robust (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n",
+               iterate,
+               stats->robustMedian, stats->robustStdev,
+               fabs(cumulative->bounds->data.F32[binMedian] - cumulative->bounds->data.F32[binMedian + 1]),
+
+               cumulative->bounds->data.F32[binMedian - 3], cumulative->bounds->data.F32[binMedian - 2],
+               cumulative->bounds->data.F32[binMedian - 1],
+               cumulative->bounds->data.F32[binMedian],
+               cumulative->bounds->data.F32[binMedian + 1],
+               cumulative->bounds->data.F32[binMedian + 2], cumulative->bounds->data.F32[binMedian + 3],
+
+               cumulative->nums->data.F32[binMedian - 3], cumulative->nums->data.F32[binMedian - 2],
+               cumulative->nums->data.F32[binMedian - 1],
+               cumulative->nums->data.F32[binMedian],
+               cumulative->nums->data.F32[binMedian + 1],
+               cumulative->nums->data.F32[binMedian + 2], cumulative->nums->data.F32[binMedian + 3]);
+        //      PS_VECTOR_PRINT_F32(histogram->bounds);
+        //      PS_VECTOR_PRINT_F32(histogram->nums);
+        //      PS_VECTOR_PRINT_F32(cumulative->bounds);
+        //      PS_VECTOR_PRINT_F32(cumulative->nums);
 #endif
 
         // ADD step 6: If the measured SIGMA is less than 2 times the bin size, exclude points which are more
         // than 25 bins from the median, recalculate the bin size, and perform the algorithm again.
-        if (sigma < (3.0 * binSize)) {
+        if (sigma < (3.0 * binSize))
+        {
             psTrace(TRACE, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
 
-            // these limits are supposed to be 25 x the raw bin size, NOT 25 of the cumulative histogram bins
-            psF32 medianLo = stats->robustMedian - 25*binSize;
-            psF32 medianHi = stats->robustMedian + 25*binSize;
+            // these limits are supposed to be 25 x the raw bin size, NOT 25 of the cumulative histogram bins
+            psF32 medianLo = stats->robustMedian - 25 * binSize;
+            psF32 medianHi = stats->robustMedian + 25 * binSize;
 
             // long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region
@@ -1107 +1248 @@
             // psTrace(TRACE, 6, "The median is at bin number %ld.  We mask bins outside the bin range (%ld:%ld)\n", binMedian, maskLo, maskHi);
             psTrace(TRACE, 6, "Masking data outside (%f %f)\n", medianLo, medianHi);
-            int Nmasked = 0;
-            for (long i = 0 ; i < myVector->n ; i++) {
-                if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) {
-                    if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & MASK_MARK) continue;
+            int Nmasked = 0;
+            for (long i = 0; i < myVector->n; i++)
+            {
+                if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi))
+                {
+                    if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & MASK_MARK)
+                        continue;
                     mask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= MASK_MARK;
                     psTrace(TRACE, 6, "Masking element %ld is %f\n", i, myVector->data.F32[i]);
-                    Nmasked ++;
+                    Nmasked++;
                 }
             }
 
-            if (Nmasked == 0) {
-                // no significant change to the sigma & binsize -- we are done here 
-                iterate = -1;
-                continue;
-            }
+            if (Nmasked == 0)
+            {
+                // no significant change to the sigma & binsize -- we are done here
+                iterate = -1;
+                continue;
+            }
 
             // Free the histograms; they will be recreated on the next iteration, with new bounds
@@ -1130 +1275 @@
             cumulative = NULL;
 
-            if (iterate >= PS_ROBUST_MAX_ITERATIONS) {
+            if (iterate >= PS_ROBUST_MAX_ITERATIONS)
+            {
                 // This occurs when a large number of the values are identical --- a bin size cannot be found
                 // that will spread out the distribution.  Therefore, set what we can, and fall over
@@ -1147 +1293 @@
                 return true;
             }
-        } else {
+        }
+        else
+        {
             // We've got the bin size correct now
             psTrace(TRACE, 6, "*************: No more iteration.  sigma is %f\n", sigma);
@@ -1153 +1301 @@
         }
     }
-    
+
     // XXX test lines while studying algorithm errors
     // fprintf (stderr, "robust stats test %7.1f +/- %7.1f : %4ld %4ld %4ld %4ld %4ld  : %f %f %f\n",
@@ -1161 +1309 @@
     // ADD step 7: Find the bins which contains the 25% and 75% data points.
     long binLo25, binHi25;
-    PS_BIN_FOR_VALUE (binLo25, cumulative->nums, totalDataPoints * 0.25f, 0);
-    PS_BIN_FOR_VALUE (binHi25, cumulative->nums, totalDataPoints * 0.75f, 0);
+    PS_BIN_FOR_VALUE(binLo25, cumulative->nums, totalDataPoints * 0.25f, 0);
+    PS_BIN_FOR_VALUE(binHi25, cumulative->nums, totalDataPoints * 0.75f, 0);
     psTrace(TRACE, 6, "The 25-percent and 75-percent data point bins are (%ld, %ld).\n", binLo25, binHi25);
 
@@ -1169 +1317 @@
     psF32 binLo25F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f);
     psF32 binHi25F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f);
-    if (isnan(binLo25F32) || isnan(binHi25F32)) {
+    if (isnan(binLo25F32) || isnan(binHi25F32))
+    {
         COUNT_WARNING(10, 100, "could not determine the robustUQ or LQ: fitLinearSearchForYThenReturnBin() returned a NAN.\n");
         goto escape;
@@ -1179 +1328 @@
             binLo25F32, binHi25F32);
     long N50 = 0;
-    for (long i = 0 ; i < myVector->n ; i++) {
+    for (long i = 0; i < myVector->n; i++)
+    {
         if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[i] &&
-            (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) {
+            (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i]))
+        {
             N50++;
         }
@@ -1226 +1377 @@
  * version follows the upper portion of the distribution until it passes 0.5*peak
  ********************/
-static bool vectorFittedStats (const psVector* myVector,
-                                  const psVector* errors,
-                                  psVector* mask,
-                                  psVectorMaskType maskVal,
-                                  psStats* stats)
+static bool vectorFittedStats(const psVector *myVector,
+                              const psVector *errors,
+                              psVector *mask,
+                              psVectorMaskType 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)) {
-        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
+    if (!(stats->results & PS_STAT_ROBUST_MEDIAN))
+    {
+        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats))
+        {
             psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
             return false;
@@ -1243 +1396 @@
 
     // If the mean is NAN, then generate a warning and set the stdev to NAN.
-    if (isnan(stats->robustMedian)) {
-        stats->fittedMean = NAN;
-        stats->fittedStdev = NAN;
-        stats->results |= PS_STAT_FITTED_MEAN;
-        stats->results |= PS_STAT_FITTED_STDEV;
+    if (isnan(stats->robustMedian))
+    {
+        stats->fittedMean = NAN;
+        stats->fittedStdev = NAN;
+        stats->results |= PS_STAT_FITTED_MEAN;
+        stats->results |= PS_STAT_FITTED_STDEV;
         return true;
     }
 
-    if (stats->robustStdev <= FLT_EPSILON) {
-        stats->fittedMean = stats->robustMedian;
-        stats->fittedStdev = stats->robustStdev;
-        stats->results |= PS_STAT_FITTED_MEAN;
-        stats->results |= PS_STAT_FITTED_STDEV;
+    if (stats->robustStdev <= FLT_EPSILON)
+    {
+        stats->fittedMean = stats->robustMedian;
+        stats->fittedStdev = stats->robustStdev;
+        stats->results |= PS_STAT_FITTED_MEAN;
+        stats->results |= PS_STAT_FITTED_STDEV;
         return true;
     }
-    if (myVector->n < 1) { printf("There are no elements in this vector.\n"); abort(); }
-    float guessStdev = stats->robustStdev;  // pass the guess sigma
-    float guessMean = stats->robustMedian;  // pass the guess mean
+    if (myVector->n < 1)
+    {
+        printf("There are no elements in this vector.\n");
+        abort();
+    }
+    float guessStdev = stats->robustStdev; // pass the guess sigma
+    float guessMean = stats->robustMedian; // pass the guess mean
 
     psTrace(TRACE, 6, "The ** starting ** guess mean  is %f.\n", guessMean);
@@ -1266 +1425 @@
 
     bool done = false;
-    for (int iteration = 0; !done && (iteration < 2); iteration ++) {
+    for (int iteration = 0; !done && (iteration < 2); iteration++)
+    {
         psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
 
         psF32 binSize = 1;
-        if (stats->options & PS_STAT_USE_BINSIZE) {
+        if (stats->options & PS_STAT_USE_BINSIZE)
+        {
             // Set initial bin size to the specified value.
             binSize = stats->binsize;
             psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize);
-        } else {
+        }
+        else
+        {
             // construct a histogram with (sigma/2 < binsize < sigma)
             // set roughly so that the lowest bins have about 2 cnts
             // Nsmallest ~ N50 / (4*dN))
-          //            psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8));
-
-          // CZW 2013-11-20: We know that the histogram is going to be basically Gaussian.
-          // Furthermore, we only use the inner +/- 2 sigma parts.  Therefore, define the
-          // binsize such that the bin at 2 sigma contains ~50 points (S/N ~ 7).  robustN50
-          // contains half the total points, so 2 * robustN50 / 50 is the fraction of all
-          // points in the 2 sigma bin.  Dance the erf() relations around, and it looks like
-          // there's a factor of about 1/20 to include.  Keep the PS_MAX to ensure we never bin
-          // wider than 1 sigma when the number of points is small.
-          psF32 dN = PS_MAX(1, (stats->robustN50 / 500.0));
-          binSize = guessStdev / dN;
+            //            psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8));
+
+            // CZW 2013-11-20: We know that the histogram is going to be basically Gaussian.
+            // Furthermore, we only use the inner +/- 2 sigma parts.  Therefore, define the
+            // binsize such that the bin at 2 sigma contains ~50 points (S/N ~ 7).  robustN50
+            // contains half the total points, so 2 * robustN50 / 50 is the fraction of all
+            // points in the 2 sigma bin.  Dance the erf() relations around, and it looks like
+            // there's a factor of about 1/20 to include.  Keep the PS_MAX to ensure we never bin
+            // wider than 1 sigma when the number of points is small.
+            psF32 dN = PS_MAX(1, (stats->robustN50 / 500.0));
+            binSize = guessStdev / dN;
         }
 
@@ -1295 +1458 @@
         float min = statsMinMax->min;
         float max = statsMinMax->max;
-        if (numValid == 0 || isnan(min) || isnan(max)) {
+        if (numValid == 0 || isnan(min) || isnan(max))
+        {
             COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
             psFree(statsMinMax);
@@ -1302 +1466 @@
 
         // If all data points have the same value, then we set the appropriate members of stats and return.
-        if (fabs(max - min) <= FLT_EPSILON) {
+        if (fabs(max - min) <= FLT_EPSILON)
+        {
             COUNT_WARNING(10, 100, "All data points have the same value: %f.\n", min);
             stats->fittedMean = min;
@@ -1314 +1479 @@
         // XXX can we calculate the binMin, binMax **before** building this histogram?
         // if the range is too absurd, adjust numBins & binSize
-        // We no longer want to reset the binSize here, as it can cause odd things.  Better to select
-        // a number of bins, and then set the min/max values to put those bins sanely around the mean.
-        //        long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize));
-        //        binSize = (max - min) / (float) numBins;
+        // We no longer want to reset the binSize here, as it can cause odd things.  Better to select
+        // a number of bins, and then set the min/max values to put those bins sanely around the mean.
+        //        long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize));
+        //        binSize = (max - min) / (float) numBins;
         psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max);
         psTrace(TRACE, 6, "The new bin size is %f.\n", binSize);
-        //        psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
-
+        //        psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
 
 #define FITTED_CLIPPING_NUM 5.0
-        if (min < guessMean - FITTED_CLIPPING_NUM * guessStdev) {
-          min = guessMean - FITTED_CLIPPING_NUM * guessStdev;
-        }
-        if (max > guessMean + FITTED_CLIPPING_NUM * guessStdev) {
-          max = guessMean + FITTED_CLIPPING_NUM * guessStdev;
-        }
-        long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize));
-        if (CZW) { printf("I've clipped: %f %f => %f %f ; %f %ld\n",guessMean,guessStdev,min,max,binSize,stats->robustN50); }
+        if (min < guessMean - FITTED_CLIPPING_NUM * guessStdev)
+        {
+            min = guessMean - FITTED_CLIPPING_NUM * guessStdev;
+        }
+        if (max > guessMean + FITTED_CLIPPING_NUM * guessStdev)
+        {
+            max = guessMean + FITTED_CLIPPING_NUM * guessStdev;
+        }
+        long numBins = PS_MAX(50, PS_MIN(10000, (max - min) / binSize));
+        if (CZW)
+        {
+            printf("I've clipped: %f %f => %f %f ; %f %ld\n", guessMean, guessStdev, min, max, binSize, stats->robustN50);
+        }
         psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
-        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
+        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal))
+        {
             COUNT_WARNING(10, 100, "Unable to generate histogram for fitted statistics v4.\n");
             psFree(histogram);
@@ -1339 +1509 @@
             goto escape;
         }
-        if (psTraceGetLevel("psLib.math") >= 8) {
+        if (psTraceGetLevel("psLib.math") >= 8)
+        {
             PS_VECTOR_PRINT_F32(histogram->nums);
         }
@@ -1347 +1518 @@
         // set the full-range upper and lower limits
         psF32 maxFitSigma = 2.0;
-        if (isfinite(stats->clipSigma)) {
+        if (isfinite(stats->clipSigma))
+        {
             maxFitSigma = fabs(stats->clipSigma);
         }
-        if (isfinite(stats->max)) {
+        if (isfinite(stats->max))
+        {
             maxFitSigma = fabs(stats->max);
         }
 
         psF32 minFitSigma = 2.0;
-        if (isfinite(stats->clipSigma)) {
+        if (isfinite(stats->clipSigma))
+        {
             minFitSigma = fabs(stats->clipSigma);
         }
-        if (isfinite(stats->min)) {
+        if (isfinite(stats->min))
+        {
             minFitSigma = fabs(stats->min);
         }
@@ -1364 +1539 @@
         // select the min and max bins, saturating on the lower and upper end-points
         long binMin, binMax;
-        PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, 0);
-        PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0);
-
-        if (binMin == binMax) {
+        PS_BIN_FOR_VALUE(binMin, histogram->bounds, guessMean - minFitSigma * guessStdev, 0);
+        PS_BIN_FOR_VALUE(binMax, histogram->bounds, guessMean + maxFitSigma * guessStdev, 0);
+
+        if (binMin == binMax)
+        {
             COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
             psFree(statsMinMax);
@@ -1375 +1551 @@
 
         // search for mode (peak of histogram within range mean-2sigma - mean+2sigma
-        long  binPeak = binMin;
+        long binPeak = binMin;
         float valPeak = histogram->nums->data.F32[binPeak];
-        for (int i = binMin; i < binMax; i++) {
-            if (histogram->nums->data.F32[i] > valPeak) {
+        for (int i = binMin; i < binMax; i++)
+        {
+            if (histogram->nums->data.F32[i] > valPeak)
+            {
                 binPeak = i;
                 valPeak = histogram->nums->data.F32[binPeak];
             }
-            psTrace (TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]);
-            if (CZW) { printf("CENTERED_HISTOGRAM: %f %f\n",
-                              PS_BIN_MIDPOINT(histogram,i),
-                              histogram->nums->data.F32[i]); }
-        }
-        psTrace (TRACE, 6, "\n");
-
-        if (CZW) { printf("Bin selection done: %ld %f %f %ld %f %f %ld %f %f\n",
-                          binMin,PS_BIN_MIDPOINT(histogram,binMin),histogram->nums->data.F32[binMin],
-                          binMax,PS_BIN_MIDPOINT(histogram,binMax),histogram->nums->data.F32[binMax],
-                          binPeak,PS_BIN_MIDPOINT(histogram,binPeak),histogram->nums->data.F32[binPeak]);
-        }
-        
+            psTrace(TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]);
+            if (CZW)
+            {
+                printf("CENTERED_HISTOGRAM: %f %f\n",
+                       PS_BIN_MIDPOINT(histogram, i),
+                       histogram->nums->data.F32[i]);
+            }
+        }
+        psTrace(TRACE, 6, "\n");
+
+        if (CZW)
+        {
+            printf("Bin selection done: %ld %f %f %ld %f %f %ld %f %f\n",
+                   binMin, PS_BIN_MIDPOINT(histogram, binMin), histogram->nums->data.F32[binMin],
+                   binMax, PS_BIN_MIDPOINT(histogram, binMax), histogram->nums->data.F32[binMax],
+                   binPeak, PS_BIN_MIDPOINT(histogram, binPeak), histogram->nums->data.F32[binPeak]);
+        }
+
         // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak
         psTrace(TRACE, 6, "The clipped numBins is %ld\n", binMax - binMin);
@@ -1402 +1585 @@
         psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]);
 
-        
-        float lowfitMean = NAN;
-        float lowfitStdev = NAN;
+        float lowfitMean = NAN;
+        float lowfitStdev = NAN;
         {
             // fit the lower half of the distribution
@@ -1410 +1592 @@
             // run up until we drop below 0.50*valPeak
             long binS = binMin;
-            long binE = PS_MIN (binPeak + 3, binMax);
-            for (int i = binPeak - 3; i >= binMin; i--) {
-                if (histogram->nums->data.F32[i] < 0.25*valPeak) {
+            long binE = PS_MIN(binPeak + 3, binMax);
+            for (int i = binPeak - 3; i >= binMin; i--)
+            {
+                if (histogram->nums->data.F32[i] < 0.25 * valPeak)
+                {
                     binS = i;
                     break;
                 }
             }
-            for (int i = binPeak + 3; i < binMax; i++) {
-                if (histogram->nums->data.F32[i] < 0.50*valPeak) {
+            for (int i = binPeak + 3; i < binMax; i++)
+            {
+                if (histogram->nums->data.F32[i] < 0.50 * valPeak)
+                {
                     binE = i;
                     break;
@@ -1431 +1617 @@
             psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates
             long j = 0;
-            for (long i = binS; i < binE; i++) {
+            for (long i = binS; i < binE; i++)
+            {
                 if (histogram->nums->data.F32[i] <= 0.0)
                     continue;
@@ -1440 +1627 @@
             }
             y->n = x->n = j;
-            
+
             // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
             // XXX this fit may fail with an error for an ill-conditioned matrix (bad data)
             // we probably should be able to handle the data errors gracefully
             psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
-            bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
+            bool status = psVectorFitPolynomial1D(poly, NULL, 0, y, NULL, x);
 #if (CZW && 1)
-            printf("CZW: LowfitPoly: %f %f %f\n",poly->coeff[0],poly->coeff[1],poly->coeff[2]);
-            for (long i = 0; i < x->n; i++) {
-              printf("CZW: Lowfit: %d %ld %f %f %f\n",
-                     status,i,x->data.F32[i],y->data.F32[i],
-                     poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
-                     poly->coeff[2] * pow(x->data.F32[i],2));
-            }
+            printf("CZW: LowfitPoly: %f %f %f\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
+            for (long i = 0; i < x->n; i++)
+            {
+                printf("CZW: Lowfit: %d %ld %f %f %f\n",
+                       status, i, x->data.F32[i], y->data.F32[i],
+                       poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
+                           poly->coeff[2] * pow(x->data.F32[i], 2));
+            }
 #endif
             psFree(x);
             psFree(y);
 
-            if (!status) {
+            if (!status)
+            {
                 psErrorClear();
                 COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
@@ -1467 +1656 @@
             }
 
-            if (poly->coeff[2] >= 0.0) {
+            if (poly->coeff[2] >= 0.0)
+            {
                 COUNT_WARNING(10, 100, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
 
@@ -1477 +1667 @@
                 // tends to be found in a single bin.  make one attempt to recover by dropping the guessStdev
                 // down by a jump and trying again
-                if (iteration == 0) {
-                    guessStdev = 0.25*guessStdev;
+                if (iteration == 0)
+                {
+                    guessStdev = 0.25 * guessStdev;
                     psTrace(TRACE, 6, "*** retry, new stdev is %f.\n", guessStdev);
                     continue;
@@ -1488 +1679 @@
 
             // calculate lower mean & stdev from parabolic fit -- use this as the result
-            lowfitStdev = sqrt(-0.5/poly->coeff[2]);
-            lowfitMean  = poly->coeff[1]*PS_SQR(lowfitStdev);
+            lowfitStdev = sqrt(-0.5 / poly->coeff[2]);
+            lowfitMean = poly->coeff[1] * PS_SQR(lowfitStdev);
 
             psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
@@ -1498 +1689 @@
         }
 
-        float fullfitMean  = NAN;
-        float fullfitStdev = NAN;
-        float minValueSym  = NAN;
-        float maxValueSym  = NAN;
-
-        // try the full fit as well:
-        {
+        float fullfitMean = NAN;
+        float fullfitStdev = NAN;
+        float minValueSym = NAN;
+        float maxValueSym = NAN;
+
+        // try the full fit as well:
+        {
             // fit a symmetric distribution
             // run up until we drop below 0.15*valPeak
@@ -1510 +1701 @@
             long binS = binMin;
             long binE = binMax;
-            for (int i = binPeak - 3; i >= binMin; i--) {
-                if (histogram->nums->data.F32[i] < 0.15*valPeak) {
+            for (int i = binPeak - 3; i >= binMin; i--)
+            {
+                if (histogram->nums->data.F32[i] < 0.15 * valPeak)
+                {
                     binS = i;
                     break;
                 }
             }
-            for (int i = binPeak + 3; i < binMax; i++) {
-                if (histogram->nums->data.F32[i] < 0.15*valPeak) {
+            for (int i = binPeak + 3; i < binMax; i++)
+            {
+                if (histogram->nums->data.F32[i] < 0.15 * valPeak)
+                {
                     binE = i;
                     break;
@@ -1531 +1726 @@
             psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates
             long j = 0;
-            for (long i = binS; i < binE; i++) {
+            for (long i = binS; i < binE; i++)
+            {
                 if (histogram->nums->data.F32[i] <= 0.0)
                     continue;
@@ -1543 +1739 @@
             // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
             psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
-            bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
+            bool status = psVectorFitPolynomial1D(poly, NULL, 0, y, NULL, x);
 #if (CZW && 1)
-            printf("CZW: FullfitPoly: %f %f %f\n",poly->coeff[0],poly->coeff[1],poly->coeff[2]);
-            for (long i = 0; i < x->n; i++) {
-              printf("CZW: Fullfit: %d %ld %f %f %f\n",
-                     status,i,x->data.F32[i],y->data.F32[i],
-                     poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
-                     poly->coeff[2] * pow(x->data.F32[i],2));
-            }
+            printf("CZW: FullfitPoly: %f %f %f\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
+            for (long i = 0; i < x->n; i++)
+            {
+                printf("CZW: Fullfit: %d %ld %f %f %f\n",
+                       status, i, x->data.F32[i], y->data.F32[i],
+                       poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
+                           poly->coeff[2] * pow(x->data.F32[i], 2));
+            }
 #endif
             psFree(x);
             psFree(y);
 
-            if (!status) {
+            if (!status)
+            {
                 psErrorClear();
                 COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
@@ -1566 +1764 @@
 
             // calculate upper mean & stdev from parabolic fit -- ignore this value
-            fullfitStdev = sqrt(-0.5/poly->coeff[2]);
-            fullfitMean = poly->coeff[1]*PS_SQR(fullfitStdev);
+            fullfitStdev = sqrt(-0.5 / poly->coeff[2]);
+            fullfitMean = poly->coeff[1] * PS_SQR(fullfitStdev);
 
 #ifndef PS_NO_TRACE
@@ -1580 +1778 @@
 
             // saturate on min or max value
-            if (fullfitMean < minValueSym) {
+            if (fullfitMean < minValueSym)
+            {
                 fullfitMean = minValueSym;
                 psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean);
@@ -1586 +1785 @@
 
             // saturate on min or max value
-            if (fullfitMean > maxValueSym) {
+            if (fullfitMean > maxValueSym)
+            {
                 fullfitMean = maxValueSym;
                 psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean);
             }
 
-            
-            psFree (poly);
-        }
-
-        // we now have the fullfit and the lowfit mean and stdev values
-        // accept the fullfit unless minValueSym < lowfitMean < fullfitMean
-
-        if (isfinite(lowfitMean) && isfinite(lowfitStdev) && (lowfitMean < fullfitMean) && (lowfitMean > minValueSym)) {
-            guessMean  = lowfitMean;
-            guessStdev = lowfitStdev;
-        } else {
-            guessMean  = fullfitMean;
-            guessStdev = fullfitStdev;
-        }
-
-        if (!isfinite(guessMean) || !isfinite(guessStdev)) {
-            guessMean  = stats->robustMedian;
-            guessStdev = stats->robustStdev;
-        }
-
-        if (guessStdev > 0.75*stats->robustStdev) {
-            done = true;
-        }
-
-        
+            psFree(poly);
+        }
+
+        // we now have the fullfit and the lowfit mean and stdev values
+        // accept the fullfit unless minValueSym < lowfitMean < fullfitMean
+
+        if (isfinite(lowfitMean) && isfinite(lowfitStdev) && (lowfitMean < fullfitMean) && (lowfitMean > minValueSym))
+        {
+            guessMean = lowfitMean;
+            guessStdev = lowfitStdev;
+        }
+        else
+        {
+            guessMean = fullfitMean;
+            guessStdev = fullfitStdev;
+        }
+
+        if (!isfinite(guessMean) || !isfinite(guessStdev))
+        {
+            guessMean = stats->robustMedian;
+            guessStdev = stats->robustStdev;
+        }
+
+        if (guessStdev > 0.75 * stats->robustStdev)
+        {
+            done = true;
+        }
+
 #if (CZW && 1)
-        printf("CZW IN FITTED: iter   %d %f \n"
-               "               low    %f %f \n"
-               "               full   %f %f \n"
-               "               robust %f %f \n"
-               "               final  %f %f\n",
-               iteration,minValueSym,
-               lowfitMean,lowfitStdev,
-               fullfitMean,fullfitStdev,
-               stats->robustMedian,stats->robustStdev,
-               guessMean,guessStdev);
+        printf("CZW IN FITTED: iter   %d %f \n"
+               "               low    %f %f \n"
+               "               full   %f %f \n"
+               "               robust %f %f \n"
+               "               final  %f %f\n",
+               iteration, minValueSym,
+               lowfitMean, lowfitStdev,
+               fullfitMean, fullfitStdev,
+               stats->robustMedian, stats->robustStdev,
+               guessMean, guessStdev);
 #endif
 
         // Clean up after fitting
-        psFree (histogram);
-        psFree (statsMinMax);
+        psFree(histogram);
+        psFree(statsMinMax);
     }
 
@@ -1655 +1858 @@
     return true;
 }
-
 
 /******************************************************************************
@@ -1668 +1870 @@
 {
     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
-    psTrace(TRACE, 5, "(histogram->nums->n, sigma) is (%d, %.2f\n", (int) histogram->nums->n, sigma);
+    psTrace(TRACE, 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) {
+    if (psTraceGetLevel("psLib.math") >= 8)
+    {
         PS_VECTOR_PRINT_F32(histogram->nums);
     }
 
-    long numBins = histogram->nums->n;  // Number of histogram bins
+    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) {
+    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.
@@ -1688 +1892 @@
                   "histograms has not been tested or used.\n");
 
-        for (long i = 0; i < numBins; i++) {
+        for (long i = 0; i < numBins; i++)
+        {
             // Determine the midpoint of bin i.
             float iMid = PS_BIN_MIDPOINT(histogram, i);
@@ -1708 +1913 @@
             //
             smooth->data.F32[i] = 0.0;
-            for (long j = jMin ; j <= jMax ; j++) {
+            for (long j = jMin; j <= jMax; j++)
+            {
                 float jMid = PS_BIN_MIDPOINT(histogram, j);
                 smooth->data.F32[i] +=
@@ -1714 +1920 @@
             }
         }
-    } else {
+    }
+    else
+    {
         //
         // We get here if the histogram is uniform.
         //
-        for (long i = 0; i < numBins; i++) {
+        for (long i = 0; i < numBins; i++)
+        {
             psF32 binSize = bounds->data.F32[1] - bounds->data.F32[0];
             psS32 gaussWidth = ((PS_GAUSS_WIDTH * sigma) / binSize);
@@ -1745 +1954 @@
             smooth->data.F32[i] = 0.0;
             float iMid = PS_BIN_MIDPOINT(histogram, i);
-            for (long j = jMin ; j <= jMax ; j++) {
+            for (long j = jMin; j <= jMax; j++)
+            {
                 float jMid = PS_BIN_MIDPOINT(histogram, j);
                 smooth->data.F32[i] +=
@@ -1753 +1963 @@
     }
 
-    if (psTraceGetLevel("psLib.math") >= 8) {
+    if (psTraceGetLevel("psLib.math") >= 8)
+    {
         PS_VECTOR_PRINT_F32(smooth);
     }
     psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
-    return(smooth);
+    return (smooth);
 }
 
@@ -1769 +1980 @@
 static void statsFree(psStats *stats)
 {
-    if (!stats) return;
-    psFree (stats->tmpData);
-    psFree (stats->tmpMask);
+    if (!stats)
+        return;
+    psFree(stats->tmpData);
+    psFree(stats->tmpMask);
     return;
 }
@@ -1778 +1990 @@
     psStatsAlloc(): This routine must create a new psStats data structure.
 *****************************************************************************/
-psStats* p_psStatsAlloc(const char *file, unsigned int lineno, const char *func, psStatsOptions options)
+psStats *p_psStatsAlloc(const char *file, unsigned int lineno, const char *func, psStatsOptions options)
 {
     psStats *stats = p_psAlloc(file, lineno, func, sizeof(psStats));
@@ -1784 +1996 @@
 
     // set initial, default values
-    psStatsInit (stats);
+    psStatsInit(stats);
 
     // these values are can be set as desired by the user.  they are not affected by
@@ -1818 +2030 @@
     stats->clippedMean = NAN;
     stats->clippedStdev = NAN;
-    stats->clippedNvalues = -1;     // XXX: This is never used
+    stats->clippedNvalues = -1; // XXX: This is never used
     stats->min = NAN;
     stats->max = NAN;
@@ -1825 +2037 @@
 }
 
-
 bool psMemCheckStats(psPtr ptr)
 {
     PS_ASSERT_PTR(ptr, false);
-    return ( psMemGetDeallocator(ptr) == (psFreeFunc)statsFree );
+    return (psMemGetDeallocator(ptr) == (psFreeFunc)statsFree);
 }
 
@@ -1847 +2058 @@
 XXX: Should we free stats if the asserts fail? NO; we don't own it (RHL).
 *****************************************************************************/
-bool psVectorStats(psStats* stats,
-                   const psVector* in,
-                   const psVector* errors,
-                   const psVector* mask,
+bool psVectorStats(psStats *stats,
+                   const psVector *in,
+                   const psVector *errors,
+                   const psVector *mask,
                    psVectorMaskType maskVal)
 {
@@ -1856 +2067 @@
     PS_ASSERT_VECTOR_NON_NULL(in, false);
     PS_ASSERT_VECTOR_NON_EMPTY(in, false);
-    if (mask) {
+    if (mask)
+    {
         PS_ASSERT_VECTORS_SIZE_EQUAL(mask, in, false);
         PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_VECTOR_MASK, false);
     }
-    if (errors) {
+    if (errors)
+    {
         PS_ASSERT_VECTORS_SIZE_EQUAL(errors, in, false);
         PS_ASSERT_VECTOR_TYPE(errors, in->type.type, false);
@@ -1866 +2079 @@
 
     // Convert types, as necessary
-    psVector *inF32 = NULL;             // Input vector of values, F32 version
-    if (in->type.type == PS_TYPE_F32) {
+    psVector *inF32 = NULL; // Input vector of values, F32 version
+    if (in->type.type == PS_TYPE_F32)
+    {
         inF32 = psMemIncrRefCounter((psPtr)in);
-    } else {
+    }
+    else
+    {
         inF32 = psVectorCopy(NULL, in, PS_TYPE_F32);
     }
-    psVector *errorsF32 = NULL;         // Input vector of errors, F32 version
-    if (errors) {
-        if (errors->type.type == PS_TYPE_F32) {
+    psVector *errorsF32 = NULL; // Input vector of errors, F32 version
+    if (errors)
+    {
+        if (errors->type.type == PS_TYPE_F32)
+        {
             errorsF32 = psMemIncrRefCounter((psPtr)errors);
-        } else {
+        }
+        else
+        {
             errorsF32 = psVectorCopy(NULL, errors, PS_TYPE_F32);
         }
     }
-    psVector *maskVector = NULL;            // Input mask vector, U8 version
-    if (mask) {
-        if (mask->type.type == PS_TYPE_VECTOR_MASK) {
+    psVector *maskVector = NULL; // Input mask vector, U8 version
+    if (mask)
+    {
+        if (mask->type.type == PS_TYPE_VECTOR_MASK)
+        {
             maskVector = psMemIncrRefCounter((psPtr)mask);
-        } else {
+        }
+        else
+        {
             maskVector = psVectorCopy(NULL, mask, PS_TYPE_VECTOR_MASK);
         }
     }
 
-    if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) {
+    if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max))
+    {
         PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->max, stats->min, false);
     }
 
-    if ((stats->options & PS_STAT_USE_BINSIZE) && (stats->min >= stats->max)) {
+    if ((stats->options & PS_STAT_USE_BINSIZE) && (stats->min >= stats->max))
+    {
         PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->binsize, 0.0, false);
     }
@@ -1903 +2129 @@
 
     // ************************************************************************
-    if (stats->options & PS_STAT_SAMPLE_MEAN) {
-        // NOTE: vectorSampleMean cannot return 'false'
-        if (!vectorSampleMean(inF32, errorsF32, maskVector, maskVal, stats)) {
+    if (stats->options & PS_STAT_SAMPLE_MEAN)
+    {
+        // NOTE: vectorSampleMean cannot return 'false'
+        if (!vectorSampleMean(inF32, errorsF32, maskVector, maskVal, stats))
+        {
             psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector sample mean");
             status &= false;
@@ -1912 +2140 @@
 
     // ************************************************************************
-    if (stats->options & (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE)) {
-        // NOTE: vectorSampleMedian only returns 'false' for very bad cases:
-        // NULL input vector, psVectorCopy failure, invalid vector type (likely programming errors)
-        if (!vectorSampleMedian(inF32, maskVector, maskVal, stats)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample median");
-            status &= false;
-        }
+    if (stats->options & (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE))
+    {
+        // NOTE: vectorSampleMedian only returns 'false' for very bad cases:
+        // NULL input vector, psVectorCopy failure, invalid vector type (likely programming errors)
+        if (!vectorSampleMedian(inF32, maskVector, maskVal, stats))
+        {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample median");
+            status &= false;
+        }
     }
 
     // ************************************************************************
-    if (stats->options & PS_STAT_SAMPLE_STDEV) {
-        // NOTE: vectorSampleStdev cannot return 'false'
-        if (!vectorSampleStdev(inF32, errorsF32, maskVector, maskVal, stats)) {
+    if (stats->options & PS_STAT_SAMPLE_STDEV)
+    {
+        // NOTE: vectorSampleStdev cannot return 'false'
+        if (!vectorSampleStdev(inF32, errorsF32, maskVector, maskVal, stats))
+        {
             psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample stdev");
             status &= false;
@@ -1930 +2162 @@
     }
 
-    if (stats->options & (PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS)) {
-        // NOTE: vectorSampleMoments cannot return 'false'
-        if (!vectorSampleMoments(inF32, maskVector, maskVal, stats)) {
+    if (stats->options & (PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS))
+    {
+        // NOTE: vectorSampleMoments cannot return 'false'
+        if (!vectorSampleMoments(inF32, maskVector, maskVal, stats))
+        {
             psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample moments");
             status &= false;
@@ -1939 +2173 @@
 
     // ************************************************************************
-    if (stats->options & (PS_STAT_MAX | PS_STAT_MIN)) {
-        // NOTE: vectorMinMax returns 0 if there are no valid values, 
-        // but this is not an error condition.  stats.min,max are set to NAN.
-        // vectorMinMax cannot raise an error
+    if (stats->options & (PS_STAT_MAX | PS_STAT_MIN))
+    {
+        // NOTE: vectorMinMax returns 0 if there are no valid values,
+        // but this is not an error condition.  stats.min,max are set to NAN.
+        // vectorMinMax cannot raise an error
         vectorMinMax(inF32, maskVector, maskVal, stats);
     }
 
     // ************************************************************************
-    if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE)) {
-        if (!vectorRobustStats(inF32, errorsF32, maskVector, maskVal, stats)) {
+    if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE))
+    {
+        if (!vectorRobustStats(inF32, errorsF32, maskVector, maskVal, stats))
+        {
             psError(PS_ERR_UNKNOWN, false, _("Failed to calculate robust statistics"));
             status &= false;
@@ -1955 +2192 @@
 
     // ************************************************************************
-    if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
-        if (!vectorFittedStats(inF32, errorsF32, maskVector, maskVal, stats)) {
+    if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV))
+    {
+        if (!vectorFittedStats(inF32, errorsF32, maskVector, maskVal, stats))
+        {
             psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
             status &= false;
@@ -1963 +2202 @@
 
     // ************************************************************************
-    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
-        if (!vectorClippedStats(inF32, errorsF32, maskVector, maskVal, stats)) {
+    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_MEDIAN) || (stats->options & PS_STAT_CLIPPED_STDEV))
+    {
+        if (!vectorClippedStats(inF32, errorsF32, maskVector, maskVal, stats))
+        {
             psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics\n");
             status &= false;
@@ -1980 +2221 @@
     PS_ASSERT_STRING_NON_EMPTY(string, PS_STAT_NONE);
 
-#define READ_STAT(NAME, SYMBOL) \
-    if (strcasecmp(string, NAME) == 0) { \
-        return SYMBOL; \
-    }
-
-    READ_STAT("MEAN",       PS_STAT_SAMPLE_MEAN);
-    READ_STAT("STDEV",      PS_STAT_SAMPLE_STDEV);
-    READ_STAT("SKEWNESS",   PS_STAT_SAMPLE_SKEWNESS);
-    READ_STAT("KURTOSIS",   PS_STAT_SAMPLE_KURTOSIS);
-    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);
-    READ_STAT("SAMPLE_MEDIAN",   PS_STAT_SAMPLE_MEDIAN);
+#define READ_STAT(NAME, SYMBOL)        \
+    if (strcasecmp(string, NAME) == 0) \
+    {                                  \
+        return SYMBOL;                 \
+    }
+
+    READ_STAT("MEAN", PS_STAT_SAMPLE_MEAN);
+    READ_STAT("STDEV", PS_STAT_SAMPLE_STDEV);
+    READ_STAT("SKEWNESS", PS_STAT_SAMPLE_SKEWNESS);
+    READ_STAT("KURTOSIS", PS_STAT_SAMPLE_KURTOSIS);
+    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);
+    READ_STAT("SAMPLE_MEDIAN", PS_STAT_SAMPLE_MEDIAN);
     READ_STAT("SAMPLE_QUARTILE", PS_STAT_SAMPLE_QUARTILE);
     READ_STAT("SAMPLE_SKEWNESS", PS_STAT_SAMPLE_SKEWNESS);
     READ_STAT("SAMPLE_KURTOSIS", PS_STAT_SAMPLE_KURTOSIS);
-    READ_STAT("ROBUST",          PS_STAT_ROBUST_MEDIAN);
-    READ_STAT("ROBUST_MEDIAN",   PS_STAT_ROBUST_MEDIAN);
-    READ_STAT("ROBUST_STDEV",    PS_STAT_ROBUST_STDEV);
+    READ_STAT("ROBUST", PS_STAT_ROBUST_MEDIAN);
+    READ_STAT("ROBUST_MEDIAN", PS_STAT_ROBUST_MEDIAN);
+    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_FITTED_STDEV);
-    READ_STAT("FITTED_V2",       PS_STAT_FITTED_MEAN);
-    READ_STAT("FITTED_MEAN_V2",  PS_STAT_FITTED_MEAN);
+    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);
+    READ_STAT("FITTED_MEAN_V2", PS_STAT_FITTED_MEAN);
     READ_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV);
-    READ_STAT("FITTED_V3",       PS_STAT_FITTED_MEAN);
-    READ_STAT("FITTED_MEAN_V3",  PS_STAT_FITTED_MEAN);
+    READ_STAT("FITTED_V3", PS_STAT_FITTED_MEAN);
+    READ_STAT("FITTED_MEAN_V3", PS_STAT_FITTED_MEAN);
     READ_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV);
-    READ_STAT("FITTED_V4",       PS_STAT_FITTED_MEAN);
-    READ_STAT("FITTED_MEAN_V4",  PS_STAT_FITTED_MEAN);
+    READ_STAT("FITTED_V4", PS_STAT_FITTED_MEAN);
+    READ_STAT("FITTED_MEAN_V4", PS_STAT_FITTED_MEAN);
     READ_STAT("FITTED_STDEV_V4", PS_STAT_FITTED_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("CLIPPED", PS_STAT_CLIPPED_MEAN);
+    READ_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN);
+    READ_STAT("CLIPPED_MEDIAN", PS_STAT_CLIPPED_MEDIAN);
+    READ_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV);
 
     psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to parse statistic: %s\n", string);
@@ -2023 +2266 @@
 psString psStatsOptionToString(psStatsOptions option)
 {
-    psString string = NULL;             // String to return
-
-#define WRITE_STAT(NAME, SYMBOL) \
-    if (option & SYMBOL) { \
+    psString string = NULL; // String to return
+
+#define WRITE_STAT(NAME, SYMBOL)              \
+    if (option & SYMBOL)                      \
+    {                                         \
         psStringAppend(&string, "%s ", NAME); \
     }
 
     // Same list as above (for psStatsFromString), but with repeat symbols removed
-    WRITE_STAT("SAMPLE_MEAN",     PS_STAT_SAMPLE_MEAN);
-    WRITE_STAT("SAMPLE_STDEV",    PS_STAT_SAMPLE_STDEV);
-    WRITE_STAT("SAMPLE_MEDIAN",   PS_STAT_SAMPLE_MEDIAN);
+    WRITE_STAT("SAMPLE_MEAN", PS_STAT_SAMPLE_MEAN);
+    WRITE_STAT("SAMPLE_STDEV", PS_STAT_SAMPLE_STDEV);
+    WRITE_STAT("SAMPLE_MEDIAN", PS_STAT_SAMPLE_MEDIAN);
     WRITE_STAT("SAMPLE_QUARTILE", PS_STAT_SAMPLE_QUARTILE);
     WRITE_STAT("SAMPLE_SKEWNESS", PS_STAT_SAMPLE_SKEWNESS);
     WRITE_STAT("SAMPLE_KURTOSIS", PS_STAT_SAMPLE_KURTOSIS);
-    WRITE_STAT("ROBUST_MEDIAN",   PS_STAT_ROBUST_MEDIAN);
-    WRITE_STAT("ROBUST_STDEV",    PS_STAT_ROBUST_STDEV);
+    WRITE_STAT("ROBUST_MEDIAN", PS_STAT_ROBUST_MEDIAN);
+    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_FITTED_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("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN);
+    WRITE_STAT("CLIPPED_MEDIAN", PS_STAT_CLIPPED_MEDIAN);
+    WRITE_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV);
 
     return string;
@@ -2051 +2296 @@
 {
     psList *subStrings = psStringSplit(string, " ,;", false); // List of sub-strings
-    if (!subStrings || psListLength(subStrings) == 0) {
+    if (!subStrings || psListLength(subStrings) == 0)
+    {
         // Nothing here
         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "No string to parse for statistics: %s\n", string);
@@ -2057 +2303 @@
         return NULL;
     }
-    psStats *stats = psStatsAlloc(0);   // Generate empty stats structure
+    psStats *stats = psStatsAlloc(0);                                                // Generate empty stats structure
     psListIterator *iterator = psListIteratorAlloc(subStrings, PS_LIST_HEAD, false); // Iterator
-    psString statString;                // Statistic string, from iteration
-    while ((statString = psListGetAndIncrement(iterator))) {
+    psString statString;                                                             // Statistic string, from iteration
+    while ((statString = psListGetAndIncrement(iterator)))
+    {
         psStatsOptions option = psStatsOptionFromString(statString);
-        if (option == 0) {
+        if (option == 0)
+        {
             psWarning("Unable to interpret statistic option: %s --- ignored.\n", statString);
             continue;
@@ -2080 +2328 @@
 psStatsOptions psStatsSingleOption(psStatsOptions option)
 {
-    switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) {
-      case PS_STAT_SAMPLE_MEAN:
-      case PS_STAT_SAMPLE_MEDIAN:
-      case PS_STAT_SAMPLE_STDEV:
-      case PS_STAT_SAMPLE_QUARTILE:
-      case PS_STAT_SAMPLE_SKEWNESS:
-      case PS_STAT_SAMPLE_KURTOSIS:
-      case PS_STAT_ROBUST_MEDIAN:
-      case PS_STAT_ROBUST_STDEV:
-      case PS_STAT_ROBUST_QUARTILE:
-      case PS_STAT_FITTED_MEAN:
-      case PS_STAT_FITTED_STDEV:
-      case PS_STAT_CLIPPED_MEAN:
-      case PS_STAT_CLIPPED_STDEV:
-      case PS_STAT_MAX:
-      case PS_STAT_MIN:
+    switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE))
+    {
+    case PS_STAT_SAMPLE_MEAN:
+    case PS_STAT_SAMPLE_MEDIAN:
+    case PS_STAT_SAMPLE_STDEV:
+    case PS_STAT_SAMPLE_QUARTILE:
+    case PS_STAT_SAMPLE_SKEWNESS:
+    case PS_STAT_SAMPLE_KURTOSIS:
+    case PS_STAT_ROBUST_MEDIAN:
+    case PS_STAT_ROBUST_STDEV:
+    case PS_STAT_ROBUST_QUARTILE:
+    case PS_STAT_FITTED_MEAN:
+    case PS_STAT_FITTED_STDEV:
+    case PS_STAT_CLIPPED_MEAN:
+    case PS_STAT_CLIPPED_MEDIAN:
+    case PS_STAT_CLIPPED_STDEV:
+    case PS_STAT_MAX:
+    case PS_STAT_MIN:
         return option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE);
-      default:
+    default:
         return 0;
     }
@@ -2107 +2357 @@
 {
     return options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN |
-                      PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN);
+                      PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_MEDIAN | PS_STAT_FITTED_MEAN);
 }
 
@@ -2116 +2366 @@
 }
 
-
 double psStatsGetValue(const psStats *stats, psStatsOptions option)
 {
     // We could call psStatsSingle to check, but it would be a waste since we effectively do it anyway
-    switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) {
-      case PS_STAT_SAMPLE_MEAN:
+    switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE))
+    {
+    case PS_STAT_SAMPLE_MEAN:
         return stats->sampleMean;
-      case PS_STAT_SAMPLE_MEDIAN:
+    case PS_STAT_SAMPLE_MEDIAN:
         return stats->sampleMedian;
-      case PS_STAT_SAMPLE_STDEV:
+    case PS_STAT_SAMPLE_STDEV:
         return stats->sampleStdev;
-      case PS_STAT_SAMPLE_SKEWNESS:
+    case PS_STAT_SAMPLE_SKEWNESS:
         return stats->sampleSkewness;
-      case PS_STAT_SAMPLE_KURTOSIS:
+    case PS_STAT_SAMPLE_KURTOSIS:
         return stats->sampleKurtosis;
-      case PS_STAT_ROBUST_MEDIAN:
+    case PS_STAT_ROBUST_MEDIAN:
         return stats->robustMedian;
-      case PS_STAT_ROBUST_STDEV:
+    case PS_STAT_ROBUST_STDEV:
         return stats->robustStdev;
-      case PS_STAT_FITTED_MEAN:
+    case PS_STAT_FITTED_MEAN:
         return stats->fittedMean;
-      case PS_STAT_FITTED_STDEV:
+    case PS_STAT_FITTED_STDEV:
         return stats->fittedStdev;
-      case PS_STAT_CLIPPED_MEAN:
+    case PS_STAT_CLIPPED_MEAN:
         return stats->clippedMean;
-      case PS_STAT_CLIPPED_STDEV:
+    case PS_STAT_CLIPPED_MEDIAN:
+        return stats->clippedMedian;
+    case PS_STAT_CLIPPED_STDEV:
         return stats->clippedStdev;
-      case PS_STAT_MAX:
+    case PS_STAT_MAX:
         return stats->max;
-      case PS_STAT_MIN:
+    case PS_STAT_MIN:
         return stats->min;
-      case PS_STAT_SAMPLE_QUARTILE:
-      case PS_STAT_ROBUST_QUARTILE:
+    case PS_STAT_SAMPLE_QUARTILE:
+    case PS_STAT_ROBUST_QUARTILE:
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot return a single quartile value; "
-                "get them yourself.\n");
+                                                  "get them yourself.\n");
         return NAN;
-      default:
+    default:
         return NAN;
     }
@@ -2161 +2413 @@
 // other private functions used above
 
-# if (0)
+#if (0)
 static psF32 QuadraticInverse(psF32 a,
                               psF32 b,
@@ -2167 +2419 @@
                               psF32 y,
                               psF32 xLo,
-                              psF32 xHi
-    )
+                              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) {
+    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) {
+    if (xLo <= x2 && x2 <= xHi)
+    {
         return x2;
     }
@@ -2185 +2438 @@
 
 static psF32 LinearInverse(psF32 a,
-                           psF32 b,
-                           psF32 y,
-                           psF32 xLo,
-                           psF32 xHi
-    )
+                           psF32 b,
+                           psF32 y,
+                           psF32 xLo,
+                           psF32 xHi)
 {
     psF64 x = (y - b) / a;
 
-    if (xLo <= x && x <= xHi) {
+    if (xLo <= x && x <= xHi)
+    {
         return x;
     }
     return 0.5 * (xLo + xHi);
 }
-# endif
-
-# if (0)
+#endif
+
+#if (0)
 /******************************************************************************
 fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): A general
@@ -2214 +2467 @@
                                                psVector *yVec,
                                                psS32 binNum,
-                                               psF32 yVal
-    )
+                                               psF32 yVal)
 {
     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
     psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
-    if (psTraceGetLevel("psLib.math") >= 8) {
+    if (psTraceGetLevel("psLib.math") >= 8)
+    {
         PS_VECTOR_PRINT_F32(xVec);
         PS_VECTOR_PRINT_F32(yVec);
@@ -2236 +2489 @@
     psF32 tmpFloat = 0.0f;
 
-    if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->n - 2))) {
+    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]));
+        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(TRACE, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1],
-                xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
+                xVec->data.F32[binNum], xVec->data.F32[binNum + 1], xVec->data.F32[binNum + 2]);
         psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
         psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
@@ -2252 +2506 @@
         // 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]))) ) {
+        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."),
@@ -2263 +2518 @@
         //
         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]))) {
+            ((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");
@@ -2274 +2530 @@
         // Determine the coefficients of the polynomial.
         psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
-        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) {
+        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x))
+        {
             psError(PS_ERR_UNEXPECTED_NULL, false,
                     _("Failed to fit a 1-dimensional polynomial to the three specified data points.  "
@@ -2287 +2544 @@
         psTrace(TRACE, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
         psTrace(TRACE, 6, "Fitted y vec is (%f %f %f)\n",
-                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
-                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
-                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
+                (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[0]),
+                (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[1]),
+                (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[2]));
 
         psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
@@ -2296 +2553 @@
         psFree(myPoly);
 
-        if (isnan(tmpFloat)) {
+        if (isnan(tmpFloat))
+        {
             psError(PS_ERR_UNEXPECTED_NULL,
                     false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
@@ -2302 +2560 @@
             psFree(y);
             psTrace(TRACE, 5, "---- %s(NAN) end ----\n", __func__);
-            return(NAN);
-        }
-    } else {
+            return (NAN);
+        }
+    }
+    else
+    {
         // These are special cases where the bin is at the beginning or end of the vector.
-        if (binNum == 0) {
+        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)) {
+        }
+        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)) {
+        }
+        else if (binNum == (xVec->n - 2))
+        {
             // XXX: Is this right?
             tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]);
@@ -2341 +2606 @@
                                                           psVector *yVec,
                                                           psS32 binNum,
-                                                          psF32 yVal
-    )
+                                                          psF32 yVal)
 {
     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
     psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
-    if (psTraceGetLevel("psLib.math") >= 8) {
+    if (psTraceGetLevel("psLib.math") >= 8)
+    {
         PS_VECTOR_PRINT_F32(xVec);
         PS_VECTOR_PRINT_F32(yVec);
@@ -2362 +2627 @@
     psF32 tmpFloat = 0.0f;
 
-    if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->n - 2))) {
+    if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->n - 2)))
+    {
         // The general case.  We have all three points.
         x->data.F64[0] = xVec->data.F32[binNum - 1];
         x->data.F64[1] = xVec->data.F32[binNum];
-        x->data.F64[2] = xVec->data.F32[binNum+1];
+        x->data.F64[2] = xVec->data.F32[binNum + 1];
         y->data.F64[0] = yVec->data.F32[binNum - 1];
         y->data.F64[1] = yVec->data.F32[binNum];
         y->data.F64[2] = yVec->data.F32[binNum + 1];
         psTrace(TRACE, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1],
-                xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
+                xVec->data.F32[binNum], xVec->data.F32[binNum + 1], xVec->data.F32[binNum + 2]);
         psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
         psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
@@ -2378 +2644 @@
         // 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]))) ) {
+        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."),
@@ -2389 +2656 @@
         //
         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]))) {
+            ((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");
@@ -2400 +2668 @@
         // Determine the coefficients of the polynomial.
         psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
-        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) {
+        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x))
+        {
             psError(PS_ERR_UNEXPECTED_NULL, false,
                     _("Failed to fit a 1-dimensional polynomial to the three specified data points.  "
@@ -2413 +2682 @@
         psTrace(TRACE, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
         psTrace(TRACE, 6, "Fitted y vec is (%f %f %f)\n",
-                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
-                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
-                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
+                (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[0]),
+                (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[1]),
+                (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[2]));
 
         psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
@@ -2422 +2691 @@
         psFree(myPoly);
 
-        if (isnan(tmpFloat)) {
+        if (isnan(tmpFloat))
+        {
             psError(PS_ERR_UNEXPECTED_NULL,
                     false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
@@ -2428 +2698 @@
             psFree(y);
             psTrace(TRACE, 5, "---- %s(NAN) end ----\n", __func__);
-            return(NAN);
-        }
-    } else {
+            return (NAN);
+        }
+    }
+    else
+    {
         // These are special cases where the bin is at the beginning or end of the vector.
-        if (binNum == 0) {
+        if (binNum == 0)
+        {
             // We have two points only at the beginning of the vectors x and y.
             // XXX this does not seem to be doing the linear interpolation / extrapolation
             tmpFloat = 0.5 * (xVec->data.F32[binNum] +
                               xVec->data.F32[binNum + 1]);
-        } else if (binNum == (xVec->n - 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)) {
+        }
+        else if (binNum == (xVec->n - 2))
+        {
             // XXX: Is this right?
             tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]);
@@ -2472 +2749 @@
                                                  psVector *yVec,
                                                  psS32 binNum,
-                                                 psF32 yVal
-    )
+                                                 psF32 yVal)
 {
     psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
-    if (psTraceGetLevel("psLib.math") >= 8) {
+    if (psTraceGetLevel("psLib.math") >= 8)
+    {
         PS_VECTOR_PRINT_F32(xVec);
         PS_VECTOR_PRINT_F32(yVec);
@@ -2495 +2772 @@
 
     //    if ((binNum >= 1) && (binNum <= (yVec->n - 2)) && (binNum <= (xVec->n - 2))) {
-    if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3))) {
+    if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3)))
+    {
         // The general case.  We have all three points.
-      //        x->data.F64[0] = binNum - 1;
-      //        x->data.F64[1] = binNum;
-      //        x->data.F64[2] = binNum + 1;
-      x->data.F64[0] = xVec->data.F32[binNum - 2];
-      x->data.F64[1] = xVec->data.F32[binNum - 1];
-      x->data.F64[2] = xVec->data.F32[binNum + 0];
-      x->data.F64[3] = xVec->data.F32[binNum + 1];
-      x->data.F64[4] = xVec->data.F32[binNum + 2];
+        //        x->data.F64[0] = binNum - 1;
+        //        x->data.F64[1] = binNum;
+        //        x->data.F64[2] = binNum + 1;
+        x->data.F64[0] = xVec->data.F32[binNum - 2];
+        x->data.F64[1] = xVec->data.F32[binNum - 1];
+        x->data.F64[2] = xVec->data.F32[binNum + 0];
+        x->data.F64[3] = xVec->data.F32[binNum + 1];
+        x->data.F64[4] = xVec->data.F32[binNum + 2];
         y->data.F64[0] = yVec->data.F32[binNum - 2];
         y->data.F64[1] = yVec->data.F32[binNum - 1];
         y->data.F64[2] = yVec->data.F32[binNum + 0];
-        y->data.F64[3] = yVec->data.F32[binNum + 1];
-        y->data.F64[4] = yVec->data.F32[binNum + 2];
-        psTrace(TRACE, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
+        y->data.F64[3] = yVec->data.F32[binNum + 1];
+        y->data.F64[4] = yVec->data.F32[binNum + 2];
+        psTrace(TRACE, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum + 1], xVec->data.F32[binNum + 2]);
         psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
         psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
 
-
         // Ensure that the y value lies within range of the y values.
-        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) ||
-               ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0]))) ) {
+        if (!(((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) ||
+              ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0]))))
+        {
             psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                     _("Specified yVal, %g, is not within y-range, %g to %g."),
@@ -2526 +2804 @@
         // 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]))) {
+            ((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");
@@ -2536 +2815 @@
         // Determine the coefficients of the polynomial.
         psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
-        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) {
+        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x))
+        {
             psError(PS_ERR_UNEXPECTED_NULL, false,
                     _("Failed to fit a 1-dimensional polynomial to the three specified data points.  "
@@ -2549 +2829 @@
         psTrace(TRACE, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
         psTrace(TRACE, 6, "Fitted y vec is (%f %f %f)\n",
-                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
-                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
-                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
+                (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[0]),
+                (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[1]),
+                (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[2]));
 
         psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
@@ -2557 +2837 @@
         psFree(myPoly);
 
-        if (isnan(binValue)) {
+        if (isnan(binValue))
+        {
             psError(PS_ERR_UNEXPECTED_NULL,
                     false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
             psFree(x);
             psFree(y);
-            return(NAN);
-        }
-        
+            return (NAN);
+        }
+
         // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1
-        //      assert (binValue >= binNum - 1);
-        //      assert (binValue <= binNum + 1);
-
-        //      int fitBin = binValue;
-        //        float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
-        //        float dY = binValue - fitBin;
-        //        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
-        tmpFloat = binValue;
-        
-    } else {
+        //      assert (binValue >= binNum - 1);
+        //      assert (binValue <= binNum + 1);
+
+        //      int fitBin = binValue;
+        //        float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
+        //        float dY = binValue - fitBin;
+        //        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
+        tmpFloat = binValue;
+    }
+    else
+    {
         // These are special cases where the bin is at the beginning or end of the vector.
-        if (binNum == 0) {
+        if (binNum == 0)
+        {
             // We have two points only at the beginning of the vectors x and y.
             // X = (dX/dY)(Y - Yo) + Xo
             float dX = xVec->data.F32[1] - xVec->data.F32[0];
             float dY = yVec->data.F32[1] - yVec->data.F32[0];
-            if (dY == 0.0) {
+            if (dY == 0.0)
+            {
                 tmpFloat = xVec->data.F32[0];
-            } else {
+            }
+            else
+            {
                 tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0];
             }
-        } else if (binNum == (xVec->n - 1)) {
+        }
+        else if (binNum == (xVec->n - 1))
+        {
             // We have two points only at the end of the vectors x and y.
             // X = (dX/dY)(Y - Yo) + Xo
-            float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1];
-            float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1];
-            if (dY == 0.0) {
-                tmpFloat = xVec->data.F32[binNum-1];
-            } else {
-                tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1];
+            float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum - 1];
+            float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum - 1];
+            if (dY == 0.0)
+            {
+                tmpFloat = xVec->data.F32[binNum - 1];
+            }
+            else
+            {
+                tmpFloat = (yVal - yVec->data.F32[binNum - 1]) * (dX / dY) + xVec->data.F32[binNum - 1];
             }
         }
@@ -2606 +2897 @@
     return tmpFloat;
 }
-# endif
-
+#endif
 
 /******************************************************************************
@@ -2623 +2913 @@
 *****************************************************************************/
 static psF32 fitLinearSearchForYThenReturnBin(const psVector *xVec,
-                                              psVector *yVec,
-                                              psS32 binNum,
-                                              psF32 yVal
-    )
+                                              psVector *yVec,
+                                              psS32 binNum,
+                                              psF32 yVal)
 {
 
-# if (1)
-# define HALF_SIZE 2
-  double Sx = 0.0;
-
-  double Sy = 0.0;
-  double Sxx = 0.0;
-  double Sxy = 0.0;
-  double deltaY = 0.0;
-  int N = 0;
-
-  for (int u = binNum - HALF_SIZE; u <= binNum + HALF_SIZE; u++) {
-    if ((u >= 0)&&(u < yVec->n)) {
-      if (u+1 < xVec->n) {
-        Sx += yVec->data.F32[u];
-        Sxx += PS_SQR(yVec->data.F32[u]);
-
-        deltaY = xVec->data.F32[u];
-        //deltaY = 0.5 * (xVec->data.F32[u] + xVec->data.F32[u+1]);
-        Sy += deltaY;
-        Sxy += yVec->data.F32[u] * deltaY;
-        N += 1;
-      }
-    }
-  }
-  double Det = N * Sxx - Sx * Sx;
-  if (Det == 0.0) return NAN;
-  if (N == 0) return NAN;
-
-  double C0 = (Sy*Sxx - Sx*Sxy) / Det;
-  double C1 = (Sxy*N - Sx*Sy) / Det;
-  
-  double value = C0 + yVal*C1;
-  return value;
-  
-  
-# else
+#if (1)
+#define HALF_SIZE 2
+    double Sx = 0.0;
+
+    double Sy = 0.0;
+    double Sxx = 0.0;
+    double Sxy = 0.0;
+    double deltaY = 0.0;
+    int N = 0;
+
+    for (int u = binNum - HALF_SIZE; u <= binNum + HALF_SIZE; u++)
+    {
+        if ((u >= 0) && (u < yVec->n))
+        {
+            if (u + 1 < xVec->n)
+            {
+                Sx += yVec->data.F32[u];
+                Sxx += PS_SQR(yVec->data.F32[u]);
+
+                deltaY = xVec->data.F32[u];
+                // deltaY = 0.5 * (xVec->data.F32[u] + xVec->data.F32[u+1]);
+                Sy += deltaY;
+                Sxy += yVec->data.F32[u] * deltaY;
+                N += 1;
+            }
+        }
+    }
+    double Det = N * Sxx - Sx * Sx;
+    if (Det == 0.0)
+        return NAN;
+    if (N == 0)
+        return NAN;
+
+    double C0 = (Sy * Sxx - Sx * Sxy) / Det;
+    double C1 = (Sxy * N - Sx * Sy) / Det;
+
+    double value = C0 + yVal * C1;
+    return value;
+
+#else
     psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
-    if (psTraceGetLevel("psLib.math") >= 8) {
+    if (psTraceGetLevel("psLib.math") >= 8)
+    {
         PS_VECTOR_PRINT_F32(xVec);
         PS_VECTOR_PRINT_F32(yVec);
@@ -2684 +2978 @@
     psF32 tmpFloat = 0.0f;
 
-    if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3))) {
-        x->data.F64[0] = xVec->data.F32[binNum - 2];
-        x->data.F64[1] = xVec->data.F32[binNum - 1];
-        x->data.F64[2] = xVec->data.F32[binNum + 0];
-        x->data.F64[3] = xVec->data.F32[binNum + 1];
-        x->data.F64[4] = xVec->data.F32[binNum + 2];
-
-        y->data.F64[0] = yVec->data.F32[binNum - 2];
-        y->data.F64[1] = yVec->data.F32[binNum - 1];
-        y->data.F64[2] = yVec->data.F32[binNum + 0];
-        y->data.F64[3] = yVec->data.F32[binNum + 1];
-        y->data.F64[4] = yVec->data.F32[binNum + 2];
-        psTrace(TRACE, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
-        psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
-        psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
-
-        // Ensure that the y value lies within range of the y values.
-        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) ||
-               ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0]))) ) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                    _("Specified yVal, %g, is not within y-range, %g to %g."),
+    if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3)))
+    {
+        x->data.F64[0] = xVec->data.F32[binNum - 2];
+        x->data.F64[1] = xVec->data.F32[binNum - 1];
+        x->data.F64[2] = xVec->data.F32[binNum + 0];
+        x->data.F64[3] = xVec->data.F32[binNum + 1];
+        x->data.F64[4] = xVec->data.F32[binNum + 2];
+
+        y->data.F64[0] = yVec->data.F32[binNum - 2];
+        y->data.F64[1] = yVec->data.F32[binNum - 1];
+        y->data.F64[2] = yVec->data.F32[binNum + 0];
+        y->data.F64[3] = yVec->data.F32[binNum + 1];
+        y->data.F64[4] = yVec->data.F32[binNum + 2];
+        psTrace(TRACE, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum + 1], xVec->data.F32[binNum + 2]);
+        psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
+        psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
+
+        // Ensure that the y value lies within range of the y values.
+        if (!(((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) ||
+              ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0]))))
+        {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    _("Specified yVal, %g, is not within y-range, %g to %g."),
                     (psF64)yVal, y->data.F64[0], y->data.F64[2]);
             return NAN;
@@ -2711 +3007 @@
         // 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]))) {
+            ((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");
@@ -2721 +3018 @@
         // Determine the coefficients of the polynomial.
         psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1);
-        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) {
+        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x))
+        {
             psError(PS_ERR_UNEXPECTED_NULL, false,
                     _("Failed to fit a 1-dimensional polynomial to the three specified data points.  "
@@ -2733 +3031 @@
         psTrace(TRACE, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
         psTrace(TRACE, 6, "Fitted y vec is (%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[0]),
+                (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[1]));
 
         psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
@@ -2740 +3038 @@
         psFree(myPoly);
 
-        if (isnan(binValue)) {
+        if (isnan(binValue))
+        {
             psError(PS_ERR_UNEXPECTED_NULL,
                     false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
             psFree(x);
             psFree(y);
-            return(NAN);
-        }
-        
+            return (NAN);
+        }
+
         // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1
-        //      assert (binValue >= binNum - 1);
-        //      assert (binValue <= binNum + 1);
-
-        //      int fitBin = binValue;
-        //        float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
-        //        float dY = binValue - fitBin;
-        //        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
-        tmpFloat = binValue;
-                
-        
-    } else {
+        //      assert (binValue >= binNum - 1);
+        //      assert (binValue <= binNum + 1);
+
+        //      int fitBin = binValue;
+        //        float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
+        //        float dY = binValue - fitBin;
+        //        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
+        tmpFloat = binValue;
+    }
+    else
+    {
         // These are special cases where the bin is at the beginning or end of the vector.
-        if (binNum == 0) {
+        if (binNum == 0)
+        {
             // We have two points only at the beginning of the vectors x and y.
             // X = (dX/dY)(Y - Yo) + Xo
             float dX = xVec->data.F32[1] - xVec->data.F32[0];
             float dY = yVec->data.F32[1] - yVec->data.F32[0];
-            if (dY == 0.0) {
+            if (dY == 0.0)
+            {
                 tmpFloat = xVec->data.F32[0];
-            } else {
+            }
+            else
+            {
                 tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0];
             }
-        } else if (binNum == (xVec->n - 1)) {
+        }
+        else if (binNum == (xVec->n - 1))
+        {
             // We have two points only at the end of the vectors x and y.
             // X = (dX/dY)(Y - Yo) + Xo
-            float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1];
-            float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1];
-            if (dY == 0.0) {
-                tmpFloat = xVec->data.F32[binNum-1];
-            } else {
-                tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1];
+            float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum - 1];
+            float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum - 1];
+            if (dY == 0.0)
+            {
+                tmpFloat = xVec->data.F32[binNum - 1];
+            }
+            else
+            {
+                tmpFloat = (yVal - yVec->data.F32[binNum - 1]) * (dX / dY) + xVec->data.F32[binNum - 1];
             }
         }
@@ -2789 +3097 @@
 
     return tmpFloat;
-# endif
+#endif
 }