Changeset 7132

Timestamp:

May 17, 2006, 3:20:51 PM (20 years ago)

Author:

Paul Price

Message:

Cleaning up code in psStats.c and psMathUtils.[ch]

Location:

trunk/psLib/src/math

Files:

• psMathUtils.c (modified) (11 diffs)
• psMathUtils.h (modified) (3 diffs)
• psStats.c (modified) (66 diffs)

trunk/psLib/src/math/psMathUtils.c

@@ -1 +1 @@
 /** @file psMathUtils.c
  *
- *  This file contains standard math rotines.
+ *  This file contains standard math routines.
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-02 21:09:07 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-05-18 01:20:51 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -48 +48 @@
 
 /*****************************************************************************
-vectorBinDisectTYPE(): This is a macro for a private function which takes as
-input an array of data as well as a single value for that data.  The input
-vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all i).
-This routine does a binary disection of the vector and returns "i" such that
-(v[i] <= x <= v[i+1).  If x lies outside the range of v[], then this routine
-prints a warning message and returns (-2 or -1).
+This is a macro covering the various types for the below function.
  *****************************************************************************/
-#define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \
-static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \
-                                   psS32 numBins, \
-                                   ps##TYPE x) \
-{ \
-    psS32 min; \
-    psS32 max; \
-    psS32 mid; \
+#define VECTOR_BIN_DISECT_CASE(TYPE) \
+case PS_TYPE_##TYPE: { \
+    ps##TYPE *bounds = bins->data.TYPE; \
+    ps##TYPE value = x->data.TYPE; \
+    long min; \
+    long max; \
+    long mid; \
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__); \
     /* psTrace(__func__, 6, "Determining the bin for: %f\n", x); */\
-    if (x < bins[0]) { \
+    if (value < bounds[0]) { \
         psLogMsg(__func__, PS_LOG_WARN, \
                  "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \
-                 #TYPE, x, bins[0], bins[numBins-1]); \
+                 #TYPE, value, bounds[0], bounds[numBins-1]); \
         return(-2); \
     } \
-    if (x > bins[numBins-1]) { \
+    if (value > bounds[numBins-1]) { \
         psLogMsg(__func__, PS_LOG_WARN, \
                  "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \
-                 #TYPE, x, bins[0], bins[numBins-1]); \
+                 #TYPE, value, bounds[0], bounds[numBins-1]); \
         return(-1); \
     } \
@@ -82 +76 @@
     mid = ((max+1)-min)/2; \
     while (min != max) { \
-        /*psTrace(__func__, 6, "(min, mid, max) is (%u, %u, %u): (x, bins[mid]) is (%f, %f)\n", min, mid, max, x, bins[mid]); */\
         \
-        if (x == bins[mid]) { \
-            /*psTrace(__func__, 6, "%.2f should be in this range (%.2f - %.2f)\n", x, bins[mid], bins[mid+1]); */\
+        if (value == bounds[mid]) { \
             psTrace(__func__, 4, "---- %s(%d) end (1) ----\n", __func__, mid); \
             return(mid); \
-        } else if (x < bins[mid]) { \
+        } else if (value < bounds[mid]) { \
             max = mid-1; \
         } else { \
@@ -95 +87 @@
         mid = ((max+1)+min)/2; \
     } \
-    /*psTrace(__func__, 6, "%.2f should be in this range (%.2f - %.2f)\n", x, bins[min], bins[min+1]);*/ \
     psTrace(__func__, 4, "---- %s(%d) end (2) ----\n", __func__, min); \
     return(min); \
-} \
-
-FUNC_MACRO_VECTOR_BIN_DISECT(S8)
-FUNC_MACRO_VECTOR_BIN_DISECT(S16)
-FUNC_MACRO_VECTOR_BIN_DISECT(S32)
-FUNC_MACRO_VECTOR_BIN_DISECT(S64)
-FUNC_MACRO_VECTOR_BIN_DISECT(U8)
-FUNC_MACRO_VECTOR_BIN_DISECT(U16)
-FUNC_MACRO_VECTOR_BIN_DISECT(U32)
-FUNC_MACRO_VECTOR_BIN_DISECT(U64)
-FUNC_MACRO_VECTOR_BIN_DISECT(F32)
-FUNC_MACRO_VECTOR_BIN_DISECT(F64)
+}
 
 /*****************************************************************************
-p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect().
+p_psVectorBinDisect(): This function takes as input an array of data as well as a single value for that data.
+The input vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all i).  This routine does a
+binary disection of the vector and returns "i" such that (v[i] <= x <= v[i+1).  If x lies outside the range of
+v[], then this routine prints a warning message and returns (-2 or -1).
   *****************************************************************************/
 psS32 p_psVectorBinDisect(
-    psVector *bins,
-    psScalar *x)
+    const psVector *bins,
+    const psScalar *x)
 {
     PS_ASSERT_VECTOR_NON_NULL(bins, -4);
@@ -122 +105 @@
     PS_ASSERT_PTR_NON_NULL(x, -6);
     PS_ASSERT_PTR_TYPE_EQUAL(x, bins, -3);
-    char* strType;
+    long numBins = bins->n;             // Number of bins
 
     switch (x->type.type) {
-    case PS_TYPE_U8:
-        return(vectorBinDisectU8(bins->data.U8, bins->n, x->data.U8));
-    case PS_TYPE_U16:
-        return(vectorBinDisectU16(bins->data.U16, bins->n, x->data.U16));
-    case PS_TYPE_U32:
-        return(vectorBinDisectU32(bins->data.U32, bins->n, x->data.U32));
-    case PS_TYPE_U64:
-        return(vectorBinDisectU64(bins->data.U64, bins->n, x->data.U64));
-    case PS_TYPE_S8:
-        return(vectorBinDisectS8(bins->data.S8, bins->n, x->data.S8));
-    case PS_TYPE_S16:
-        return(vectorBinDisectS16(bins->data.S16, bins->n, x->data.S16));
-    case PS_TYPE_S32:
-        return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
-    case PS_TYPE_S64:
-        return(vectorBinDisectS64(bins->data.S64, bins->n, x->data.S64));
-    case PS_TYPE_F32:
-        return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
-    case PS_TYPE_F64:
-        return(vectorBinDisectF64(bins->data.F64, bins->n, x->data.F64));
-    default:
-        PS_TYPE_NAME(strType, x->type.type);
-        psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
-        return -1;
+        VECTOR_BIN_DISECT_CASE(S8);
+        VECTOR_BIN_DISECT_CASE(S16);
+        VECTOR_BIN_DISECT_CASE(S32);
+        VECTOR_BIN_DISECT_CASE(S64);
+        VECTOR_BIN_DISECT_CASE(U8);
+        VECTOR_BIN_DISECT_CASE(U16);
+        VECTOR_BIN_DISECT_CASE(U32);
+        VECTOR_BIN_DISECT_CASE(U64);
+        VECTOR_BIN_DISECT_CASE(F32);
+        VECTOR_BIN_DISECT_CASE(F64);
+    default: {
+            char* strType;
+            PS_TYPE_NAME(strType, x->type.type);
+            psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
+            return -1;
+        }
     }
     return(-3);
@@ -154 +129 @@
 
 
-
-/*****************************************************************************
-fullInterpolateTYPE(): This routine will take as input n-element psVectors
-domain and range, and the x value, assumed to lie with the domain vector.  It
-produces as output the (order-1)-order LaGrange interpolated value of x.
- 
-XXX: do we error check for non-distinct domain values?
- 
-This routine will only be called inside this file.  So, we do not have to
-do PS_ASSERTS.
- *****************************************************************************/
-#define FUNC_MACRO_FULL_INTERPOLATE(TYPE) \
-static psScalar *vectorInterpolate##TYPE( \
-        psScalar *out, \
-        psVector *domain, \
-        psVector *range, \
-        psS32 order, \
-        psScalar *x) \
-{ \
+/*************************************************************************************************************
+Helper macro for p_psVectorInterpolate: handles each of the types
+*************************************************************************************************************/
+#define VECTOR_INTERPOLATE_CASE(TYPE) \
+case PS_TYPE_##TYPE: { \
     psTrace(__func__, 4, "---- %s() begin %u-order.) (%d data points) ----\n", __func__, order, order+1); \
     if (x->data.TYPE < domain->data.TYPE[0]) { \
@@ -213 +174 @@
     psTrace(__func__, 4, "---- %s(....) end ----\n", __func__); \
     return(out); \
-} \
-
-FUNC_MACRO_FULL_INTERPOLATE(U8)
-FUNC_MACRO_FULL_INTERPOLATE(U16)
-FUNC_MACRO_FULL_INTERPOLATE(U32)
-FUNC_MACRO_FULL_INTERPOLATE(U64)
-FUNC_MACRO_FULL_INTERPOLATE(S8)
-FUNC_MACRO_FULL_INTERPOLATE(S16)
-FUNC_MACRO_FULL_INTERPOLATE(S32)
-FUNC_MACRO_FULL_INTERPOLATE(S64)
-FUNC_MACRO_FULL_INTERPOLATE(F32)
-FUNC_MACRO_FULL_INTERPOLATE(F64)
+}
 
 /*****************************************************************************
@@ -236 +186 @@
 psScalar *p_psVectorInterpolate(
     psScalar *out,
-    psVector *domain,
-    psVector *range,
+    const psVector *domain,
+    const psVector *range,
     psS32 order,
-    psScalar *x)
+    const psScalar *x)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
@@ -249 +199 @@
     PS_ASSERT_PTR_TYPE_EQUAL(domain, range, NULL);
     PS_ASSERT_PTR_TYPE_EQUAL(domain, x, NULL);
-    if (out == NULL) {
+    if (!out) {
         out = psScalarAlloc(0, x->type.type);
     } else {
         PS_ASSERT_PTR_TYPE_EQUAL(domain, out, NULL);
     }
-    char* strType;
 
     switch (x->type.type) {
-    case PS_TYPE_U8:
-        return(vectorInterpolateU8(out, domain, range, order, x));
-    case PS_TYPE_U16:
-        return(vectorInterpolateU16(out, domain, range, order, x));
-    case PS_TYPE_U32:
-        return(vectorInterpolateU32(out, domain, range, order, x));
-    case PS_TYPE_U64:
-        return(vectorInterpolateU64(out, domain, range, order, x));
-    case PS_TYPE_S8:
-        return(vectorInterpolateS8(out, domain, range, order, x));
-    case PS_TYPE_S16:
-        return(vectorInterpolateS16(out, domain, range, order, x));
-    case PS_TYPE_S32:
-        return(vectorInterpolateS32(out, domain, range, order, x));
-    case PS_TYPE_S64:
-        return(vectorInterpolateS64(out, domain, range, order, x));
-    case PS_TYPE_F32:
-        return(vectorInterpolateF32(out, domain, range, order, x));
-    case PS_TYPE_F64:
-        return(vectorInterpolateF64(out, domain, range, order, x));
-    default:
-        PS_TYPE_NAME(strType, x->type.type);
-        psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
-        return NULL;
+        VECTOR_INTERPOLATE_CASE(U8);
+        VECTOR_INTERPOLATE_CASE(U16);
+        VECTOR_INTERPOLATE_CASE(U32);
+        VECTOR_INTERPOLATE_CASE(U64);
+        VECTOR_INTERPOLATE_CASE(S8);
+        VECTOR_INTERPOLATE_CASE(S16);
+        VECTOR_INTERPOLATE_CASE(S32);
+        VECTOR_INTERPOLATE_CASE(S64);
+        VECTOR_INTERPOLATE_CASE(F32);
+        VECTOR_INTERPOLATE_CASE(F64);
+    default: {
+            char* strType;
+            PS_TYPE_NAME(strType, x->type.type);
+            psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
+            return NULL;
+        }
     }
 
@@ -287 +228 @@
 
 /*****************************************************************************
-These macros and functions define the following functions:
- 
-    p_psNormalizeVectorRange(myData, low, high)
- 
-That assumes that the low/high arguments are PS_TYPE_F64; the vector myData
-can be of any type.  Arguments low/high will be converted to the appropriate
-type and one of the type-specific functions below will be called:
- 
-    p_psNormalizeVectorRangeU8(myData, low, high)
-    p_psNormalizeVectorRangeU16(myData, low, high)
-    p_psNormalizeVectorRangeU32(myData, low, high)
-    p_psNormalizeVectorRangeU64(myData, low, high)
-    p_psNormalizeVectorRangeS8(myData, low, high)
-    p_psNormalizeVectorRangeS16(myData, low, high)
-    p_psNormalizeVectorRangeS32(myData, low, high)
-    p_psNormalizeVectorRangeS64(myData, low, high)
-    p_psNormalizeVectorRangeF32(myData, low, high)
-    p_psNormalizeVectorRangeF64(myData, low, high)
+Helper macro for p_psNormalizeVectorRange: handles each of the types
  *****************************************************************************/
-#define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \
-static void p_psNormalizeVectorRange##TYPE( \
-        psVector* myData, \
-        ps##TYPE outLow, \
-        ps##TYPE outHigh) \
-{ \
-    ps##TYPE min = (ps##TYPE) PS_MAX_##TYPE; \
-    ps##TYPE max = (ps##TYPE) -PS_MAX_##TYPE; \
-    psS32 i = 0; \
-    \
-    for (i = 0; i < myData->n; i++) { \
+#define NORMALIZE_VECTOR_RANGE_CASE(TYPE) \
+case PS_TYPE_##TYPE: { \
+    ps##TYPE low = outLow; \
+    ps##TYPE high = outHigh; \
+    ps##TYPE min = (ps##TYPE)PS_MAX_##TYPE; \
+    ps##TYPE max = (ps##TYPE)-PS_MAX_##TYPE; \
+    \
+    for (long i = 0; i < myData->n; i++) { \
         if (myData->data.TYPE[i] < min) { \
             min = myData->data.TYPE[i]; \
@@ -327 +248 @@
     /* Ensure that max!=min before we divide by (max-min) */ \
     if (max != min) { \
-        for (i = 0; i < myData->n; i++) { \
-            myData->data.TYPE[i] = (outLow + (myData->data.TYPE[i] - min) * \
-                                    (outHigh - outLow) / (max - min)); \
+        for (long i = 0; i < myData->n; i++) { \
+            myData->data.TYPE[i] = (low + (myData->data.TYPE[i] - min) * \
+                                    (high - low) / (max - min)); \
         } \
     } else { \
         psLogMsg(__func__, PS_LOG_WARN, "WARNING: (max==min).  Setting all elements to min.\n"); \
-        for (i = 0; i < myData->n; i++) \
+        for (long i = 0; i < myData->n; i++) \
         { \
             \
-            myData->data.TYPE[i] = outLow; \
+            myData->data.TYPE[i] = low; \
             \
         } \
     } \
+    break; \
 } \
 
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U8)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U16)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U32)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U64)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S8)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S16)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S32)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S64)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F32)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64)
-
-psS32 p_psNormalizeVectorRange(psVector* myData,
-                               psF64 outLow,
-                               psF64 outHigh)
+/*************************************************************************************************************
+p_psNormalizeVectorRange(myData, low, high): this function normalises the vector (myData) to the range
+low:
+high.
+*************************************************************************************************************/
+bool p_psNormalizeVectorRange(psVector* myData,
+                              psF64 outLow,
+                              psF64 outHigh)
 {
-    PS_ASSERT_VECTOR_NON_NULL(myData, -1);
-    char* strType;
-
+    PS_ASSERT_VECTOR_NON_NULL(myData, false);
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
+
     switch (myData->type.type) {
-    case PS_TYPE_U8:
-        p_psNormalizeVectorRangeU8(myData, (psU8) outLow, (psU8) outHigh);
-        break;
-    case PS_TYPE_U16:
-        p_psNormalizeVectorRangeU16(myData, (psU16) outLow, (psU16) outHigh);
-        break;
-    case PS_TYPE_U32:
-        p_psNormalizeVectorRangeU32(myData, (psU32) outLow, (psU32) outHigh);
-        break;
-    case PS_TYPE_U64:
-        p_psNormalizeVectorRangeU64(myData, (psU64) outLow, (psU64) outHigh);
-        break;
-    case PS_TYPE_S8:
-        p_psNormalizeVectorRangeS8(myData, (psS8) outLow, (psS8) outHigh);
-        break;
-    case PS_TYPE_S16:
-        p_psNormalizeVectorRangeS16(myData, (psS16) outLow, (psS16) outHigh);
-        break;
-    case PS_TYPE_S32:
-        p_psNormalizeVectorRangeS32(myData, (psS32) outLow, (psS32) outHigh);
-        break;
-    case PS_TYPE_S64:
-        p_psNormalizeVectorRangeS64(myData, (psS64) outLow, (psS64) outHigh);
-        break;
-    case PS_TYPE_F32:
-        p_psNormalizeVectorRangeF32(myData, (psF32) outLow, (psF32) outHigh);
-        break;
-    case PS_TYPE_F64:
-        p_psNormalizeVectorRangeF64(myData, (psF64) outLow, (psF64) outHigh);
-        break;
-    case PS_TYPE_C32:
-    case PS_TYPE_C64:
-    default:
-        PS_TYPE_NAME(strType, myData->type.type);
-        psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
-        break;
+        NORMALIZE_VECTOR_RANGE_CASE(U8);
+        NORMALIZE_VECTOR_RANGE_CASE(U16);
+        NORMALIZE_VECTOR_RANGE_CASE(U32);
+        NORMALIZE_VECTOR_RANGE_CASE(U64);
+        NORMALIZE_VECTOR_RANGE_CASE(S8);
+        NORMALIZE_VECTOR_RANGE_CASE(S16);
+        NORMALIZE_VECTOR_RANGE_CASE(S32);
+        NORMALIZE_VECTOR_RANGE_CASE(S64);
+        NORMALIZE_VECTOR_RANGE_CASE(F32);
+        NORMALIZE_VECTOR_RANGE_CASE(F64);
+    default: {
+            char* strType;
+            PS_TYPE_NAME(strType, myData->type.type);
+            psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
+            break;
+        }
     }
     psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-    return(0);
-}
-
-
+    return true;
+}
+
+

trunk/psLib/src/math/psMathUtils.h

@@ -7 +7 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-17 00:56:48 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-05-18 01:20:51 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -34 +34 @@
  */
 psS32 p_psVectorBinDisect(
-    psVector *bins,                    ///< Array of non-decreasing values
-    psScalar *x                        ///< Target value to find
+    const psVector *bins,               ///< Array of non-decreasing values
+    const psScalar *x                   ///< Target value to find
 );
 
@@ -44 +44 @@
  */
 psScalar *p_psVectorInterpolate(
-    psScalar *out,
-    psVector *domain,                  ///< Domain (x coords) for interpolation
-    psVector *range,                   ///< Range (y coords) for interpolation
-    psS32 order,                       ///< Order of interpolation function
-    psScalar *x                        ///< Location at which to evaluate
+    psScalar *out,                      ///< Output scalar, or NULL
+    const psVector *domain,             ///< Domain (x coords) for interpolation
+    const psVector *range,              ///< Range (y coords) for interpolation
+    psS32 order,                        ///< Order of interpolation function
+    const psScalar *x                   ///< Location at which to evaluate
 );
 
-psS32 p_psNormalizeVectorRange(psVector* myData,
-                               psF64 outLow,
-                               psF64 outHigh);
+bool p_psNormalizeVectorRange(psVector* myData,
+                              psF64 outLow,
+                              psF64 outHigh);
 
 /** \} */ // End of MathGroup Functions

trunk/psLib/src/math/psStats.c

@@ -16 +16 @@
  * use ->min and ->max (PS_STAT_USE_RANGE)
  *
- *  @version $Revision: 1.171 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-05-03 23:16:16 $
+ *  @version $Revision: 1.172 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-05-18 01:20:50 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -28 +28 @@
 #include <float.h>
 #include <math.h>
+#include <limits.h>
 
 /*****************************************************************************/
@@ -50 +51 @@
 /* DEFINE STATEMENTS                                                         */
 /*****************************************************************************/
-#define PS_GAUSS_WIDTH 5                // The width of the Gaussian smoothing.
+#define PS_GAUSS_WIDTH 3                // The width of the Gaussian smoothing.
 // This corresponds to N in the ADD.
 #define PS_CLIPPED_NUM_ITER_LB 1
@@ -63 +64 @@
 /* TYPE DEFINITIONS                                                          */
 /*****************************************************************************/
-psVector* p_psConvertToF32(psVector* in);
+
+// None
+
 /*****************************************************************************/
 /* GLOBAL VARIABLES                                                          */
@@ -154 +157 @@
  *****************************************************************************/
 /******************************************************************************
-p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
+vectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
 mean of the input vector.  If there was a problem with the mean calculation,
 this routine sets stats->sampleMean to NAN.
  *****************************************************************************/
-psS32 p_psVectorSampleMean(const psVector* myVector,
-                           const psVector* errors,
-                           const psVector* maskVector,
-                           psU32 maskVal,
-                           psStats* stats)
+static bool vectorSampleMean(const psVector* myVector,
+                             const psVector* errors,
+                             const psVector* maskVector,
+                             psMaskType maskVal,
+                             psStats* stats)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
 
-    psS32 i = 0;                // Loop index variable
-    psF32 mean = 0.0;           // The mean
-    psS32 count = 0;            // # of points in this mean
-
-    // If PS_STAT_USE_RANGE is requested, then we enter a slightly different
-    // loop.
-    if (errors == NULL) {
+    psF32 mean = 0.0;                   // The mean
+    long count = 0;                     // Number of points contributing to this mean
+
+    // If PS_STAT_USE_RANGE is requested, then we enter a slightly different loop.
+    if (!errors) {
         if (stats->options & PS_STAT_USE_RANGE) {
-            if (maskVector != NULL) {
-                for (i = 0; i < myVector->n; i++) {
+            if (maskVector) {
+                // No errors, check range, check mask
+                for (long i = 0; i < myVector->n; i++) {
                     // Check if the data is with the specified range
                     if (!(maskVal & maskVector->data.U8[i]) &&
@@ -184 +186 @@
                     }
                 }
-                if (count != 0) {
-                    mean /= (psF32)count;
-                } else {
-                    mean = NAN;
-                }
 
             } else {
-                for (i = 0; i < myVector->n; i++) {
+                // No errors, check range, no mask
+                for (long i = 0; i < myVector->n; i++) {
                     if ((stats->min <= myVector->data.F32[i]) &&
                             (myVector->data.F32[i] <= stats->max)) {
@@ -198 +196 @@
                     }
                 }
-                if (count != 0) {
-                    mean /= (psF32)count;
-                } else {
-                    mean = NAN;
-                }
-            }
-        } else {
-            if (maskVector != NULL) {
-                for (i = 0; i < myVector->n; i++) {
+            }
+        } else {
+            if (maskVector) {
+                // No errors, no range check, check mask
+                for (long i = 0; i < myVector->n; i++) {
                     if (!(maskVal & maskVector->data.U8[i])) {
                         mean += myVector->data.F32[i];
@@ -212 +206 @@
                     }
                 }
-                if (count != 0) {
-                    mean /= (psF32)count;
-                } else {
-                    mean = NAN;
-                }
+
             } else {
-                for (i = 0; i < myVector->n; i++) {
+                // No errors, no range check, no mask
+                for (long i = 0; i < myVector->n; i++) {
                     mean += myVector->data.F32[i];
                 }
-                mean /= (psF32)myVector->n;
-            }
-        }
+                count = myVector->n;
+            }
+        }
+
+        if (count > 0) {
+            mean /= (psF32)count;
+        } else {
+            mean = NAN;
+        }
+
     } else {
-        psF32 errorDivisor = 0.0;
-        psF32 errorSqr = 0.0;
+        psF32 sumWeights = 0.0;         // The sum of the weights
         if (stats->options & PS_STAT_USE_RANGE) {
             if (maskVector != NULL) {
-                for (i = 0; i < myVector->n; i++) {
+                for (long i = 0; i < myVector->n; i++) {
                     // Check if the data is with the specified range
                     if (!(maskVal & maskVector->data.U8[i]) &&
                             (stats->min <= myVector->data.F32[i]) &&
                             (myVector->data.F32[i] <= stats->max)) {
-                        errorSqr = errors->data.F32[i]*errors->data.F32[i];
-                        mean += myVector->data.F32[i]/errorSqr;
-                        errorDivisor+= (1.0 / errorSqr);
+                        float weight = 1.0 / PS_SQR(errors->data.F32[i]);
+                        mean += myVector->data.F32[i] * weight;
+                        sumWeights += weight;
                         count++;
                     }
                 }
-                if (count != 0) {
-                    mean /= errorDivisor;
-                } else {
-                    mean = NAN;
-                }
-
             } else {
-                for (i = 0; i < myVector->n; i++) {
+                for (long i = 0; i < myVector->n; i++) {
                     if ((stats->min <= myVector->data.F32[i]) &&
                             (myVector->data.F32[i] <= stats->max)) {
-                        errorSqr = errors->data.F32[i]*errors->data.F32[i];
-                        mean += myVector->data.F32[i]/errorSqr;
-                        errorDivisor+= (1.0 / errorSqr);
+                        float weight = 1.0 / PS_SQR(errors->data.F32[i]);
+                        mean += myVector->data.F32[i] * weight;
+                        sumWeights += weight;
                         count++;
                     }
                 }
-                if (count != 0) {
-                    mean /= errorDivisor;
-                } else {
-                    mean = NAN;
-                }
             }
         } else {
             if (maskVector != NULL) {
-                for (i = 0; i < myVector->n; i++) {
+                for (long i = 0; i < myVector->n; i++) {
                     if (!(maskVal & maskVector->data.U8[i])) {
-                        errorSqr = errors->data.F32[i]*errors->data.F32[i];
-                        mean += myVector->data.F32[i]/errorSqr;
-                        errorDivisor+= (1.0 / errorSqr);
+                        float weight = 1.0 / PS_SQR(errors->data.F32[i]);
+                        mean += myVector->data.F32[i] * weight;
+                        sumWeights += weight;
                         count++;
                     }
                 }
-                if (count != 0) {
-                    mean /= errorDivisor;
-                } else {
-                    mean = NAN;
-                }
             } else {
-                for (i = 0; i < myVector->n; i++) {
-                    errorSqr = errors->data.F32[i]*errors->data.F32[i];
-                    mean += myVector->data.F32[i]/errorSqr;
-                    errorDivisor+= (1.0 / errorSqr);
-                }
-                mean /= errorDivisor;
-            }
+                for (long i = 0; i < myVector->n; i++) {
+                    float weight = 1.0 / PS_SQR(errors->data.F32[i]);
+                    mean += myVector->data.F32[i] * weight;
+                    sumWeights += weight;
+                }
+                count = myVector->n;
+            }
+        }
+
+        if (count > 0) {
+            mean /= sumWeights;
+        } else {
+            mean = NAN;
         }
     }
@@ -290 +277 @@
     stats->sampleMean = mean;
     if (isnan(mean)) {
-        psTrace(__func__, 4, "---- %s(-1) end ----\n", __func__);
-        return(-1);
-    } else {
-        psTrace(__func__, 4, "---- %s(0) end ----\n", __func__);
-        return(0);
-    }
-
+        psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
+        return false;
+    }
+
+    psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
+    return true;
 }
 
 /******************************************************************************
-p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the
-max of the input vector.  If there was a problem with the max calculation,
-this routine sets stats->max to NAN.
+vectorMinMax(myVector, maskVector, maskVal, stats): calculates the minimum and maximum of the input vector.
+If there was a problem with the calculation, this routine sets stats->max and stats->min to NAN.  Return the
+number of valid values in the vector (those not masked or outside the specified range.
  *****************************************************************************/
-psS32 p_psVectorMax(const psVector* myVector,
-                    const psVector* maskVector,
-                    psU32 maskVal,
-                    psStats* stats)
+static long vectorMinMax(const psVector* myVector,
+                         const psVector* maskVector,
+                         psMaskType maskVal,
+                         psStats* stats
+                        )
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    psS32 i = 0;                // Loop index variable
-    psF32 max = -PS_MAX_F32;    // The calculated maximum
-    psS32 empty = true;         // Does this vector have valid elements?
+    psF32 max = -PS_MAX_F32;            // The calculated maximum
+    psF32 min = PS_MAX_F32;             // The calculated minimum
+    psF32 *vector = myVector->data.F32; // Dereference the vector
+    psU8 *mask = NULL;                  // Dereference the mask
+    if (maskVector) {
+        mask = maskVector->data.U8;
+    }
+    int num = myVector->n;              // Number of values
+    int numValid = 0;                   // Number of valid values
 
     // If PS_STAT_USE_RANGE is requested, then we enter a different loop.
     if (stats->options & PS_STAT_USE_RANGE) {
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    if ((myVector->data.F32[i] > max) &&
-                            (stats->min <= myVector->data.F32[i]) &&
-                            (myVector->data.F32[i] <= stats->max)) {
-                        max = myVector->data.F32[i];
-                        empty = false;
+        if (mask) {
+            // Need to check range and mask
+            for (long i = 0; i < num; i++) {
+                if (!(maskVal & mask[i]) && vector[i] >= stats->min && vector[i] <= stats->max) {
+                    numValid++;
+                    if (vector[i] > max) {
+                        max = vector[i];
                     }
-                }
-            }
-        } else {
-            for (i = 0; i < myVector->n; i++) {
-                if ((myVector->data.F32[i] > max) &&
-                        (stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    max = myVector->data.F32[i];
-                    empty = false;
+                    if (vector[i] < min) {
+                        min = vector[i];
+                    }
+                }
+            }
+        } else {
+            // Need to check range
+            for (long i = 0; i < num; i++) {
+                if (vector[i] >= stats->min && vector[i] <= stats->max) {
+                    numValid++;
+                    if (vector[i] > max) {
+                        max = vector[i];
+                    }
+                    if (vector[i] < min) {
+                        min = vector[i];
+                    }
                 }
             }
         }
     } else {
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    if (myVector->data.F32[i] > max) {
-                        max = myVector->data.F32[i];
-                        empty = false;
+        if (mask) {
+            // Need to check mask
+            for (long i = 0; i < num; i++) {
+                if (!(maskVal & mask[i])) {
+                    numValid++;
+                    if (vector[i] > max) {
+                        max = vector[i];
                     }
-                }
-            }
-        } else {
-            for (i = 0; i < myVector->n; i++) {
-                if (myVector->data.F32[i] > max) {
-                    max = myVector->data.F32[i];
-                    empty = false;
-                }
-            }
-        }
-    }
-    if (empty == false) {
+                    if (vector[i] < min) {
+                        min = vector[i];
+                    }
+                }
+            }
+        } else {
+            // No checks
+            for (long i = 0; i < num; i++) {
+                if (vector[i] > max) {
+                    max = vector[i];
+                }
+                if (vector[i] < min) {
+                    min = vector[i];
+                }
+            }
+            numValid = num;
+        }
+    }
+
+
+    if (numValid == 0) {
+        stats->max = NAN;
+        stats->min = NAN;
+    } else {
         stats->max = max;
-    } else {
-        stats->max = NAN;
-        psTrace(__func__, 4, "---- %s(1) end ----\n", __func__);
-        return(1);
-    }
-    psTrace(__func__, 4, "---- %s(0) end ----\n", __func__);
-    return(0);
+        stats->min = min;
+    }
+    psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, numValid);
+    return numValid;
 }
 
+#if 0
 /******************************************************************************
-p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the
-minimum of the input vector.  If there was a problem with the min calculation,
-this routine sets stats->min to NAN.
- *****************************************************************************/
-psS32 p_psVectorMin(const psVector* myVector,
-                    const psVector* maskVector,
-                    psU32 maskVal,
-                    psStats* stats)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    psS32 i = 0;                // Loop index variable
-    psF32 min = PS_MAX_F32;   // The calculated maximum
-    psS32 empty = true;         // Does this vector have valid elements?
-
-    if (stats->options & PS_STAT_USE_RANGE) {
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    if ((myVector->data.F32[i] < min) &&
-                            (stats->min <= myVector->data.F32[i]) &&
-                            (myVector->data.F32[i] <= stats->max)) {
-                        min = myVector->data.F32[i];
-                        empty = false;
-                    }
-                }
-            }
-        } else {
-            for (i = 0; i < myVector->n; i++) {
-                if ((myVector->data.F32[i] < min) &&
-                        (stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    min = myVector->data.F32[i];
-                    empty = false;
-                }
-            }
-        }
-    } else {
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    if (myVector->data.F32[i] < min) {
-                        min = myVector->data.F32[i];
-                        empty = false;
-                    }
-                }
-            }
-        } else {
-            for (i = 0; i < myVector->n; i++) {
-                if (myVector->data.F32[i] < min) {
-                    min = myVector->data.F32[i];
-                    empty = false;
-                }
-            }
-        }
-    }
-
-    if (empty == false) {
-        stats->min = min;
-    } else {
-        stats->min = NAN;
-        psTrace(__func__, 4, "---- %s(1) end ----\n", __func__);
-        return(1);
-    }
-    psTrace(__func__, 4, "---- %s(0) end ----\n", __func__);
-    return(0);
-}
-
-/******************************************************************************
-p_psVectorNValues(myVector, maskVector, maskVal, stats): This routine returns
+vectorCheckNonEmpty(myVector, maskVector, maskVal, stats): This routine returns
 "true" if the inputPsVector has 1 or more valid elements (a valid element is an
 unmasked element within the specified min/max range).  Otherwise, return
 "false".
  *****************************************************************************/
-bool p_psVectorCheckNonEmpty(
-    const psVector* myVector,
-    const psVector* maskVector,
-    psU32 maskVal,
-    psStats* stats)
+static bool vectorCheckNonEmpty(const psVector* myVector,
+                                const psVector* maskVector,
+                                psMaskType maskVal,
+                                psStats* stats)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_VECTOR_NON_NULL(myVector, -1);
-    PS_ASSERT_PTR_NON_NULL(stats, -1);
-    psS32 i = 0;                // Loop index variable
+    PS_ASSERT_VECTOR_NON_NULL(myVector, false);
+    PS_ASSERT_PTR_NON_NULL(stats, false);
 
     if (stats->options & PS_STAT_USE_RANGE) {
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
+        if (maskVector) {
+            for (long i = 0; i < myVector->n; i++) {
                 if (!(maskVal & maskVector->data.U8[i]) &&
                         (stats->min <= myVector->data.F32[i]) &&
                         (myVector->data.F32[i] <= stats->max)) {
                     psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
-                    return(true);
-                }
-            }
-        } else {
-            for (i = 0; i < myVector->n; i++) {
+                    return true;
+                }
+            }
+        } else {
+            for (long i = 0; i < myVector->n; i++) {
                 if ((stats->min <= myVector->data.F32[i]) &&
                         (myVector->data.F32[i] <= stats->max)) {
                     psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
-                    return(true);
+                    return true;
                 }
             }
@@ -473 +413 @@
     } else {
         if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
+            for (long i = 0; i < myVector->n; i++) {
                 if (!(maskVal & maskVector->data.U8[i])) {
                     psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
-                    return(true);
+                    return true;
                 }
             }
@@ -482 +422 @@
             if (myVector->n > 0) {
                 psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
-                return(true);
+                return true;
             }
         }
@@ -489 +429 @@
     return(false);
 }
-
+#endif
 
 /******************************************************************************
-p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the
-number of non-masked pixels in the vector that fall within the min/max
-range, if specified.
+vectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the median
+and quartiles of the input vector.  Returns true on success (including if there
+were no valid input vector elements).
  *****************************************************************************/
-psS32 p_psVectorNValues(const psVector* myVector,
-                        const psVector* maskVector,
-                        psU32 maskVal,
-                        psStats* stats)
+static bool vectorSampleMedian(const psVector* myVector,
+                               const psVector* maskVector,
+                               psMaskType maskVal,
+                               psStats* stats)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_VECTOR_NON_NULL(myVector, -1);
-    PS_ASSERT_PTR_NON_NULL(stats, -1);
-    psS32 i = 0;                // Loop index variable
-    psS32 numData = 0;          // The number of data points
-
-    if (stats->options & PS_STAT_USE_RANGE) {
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i]) &&
-                        (stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    numData++;
-                }
-            }
-        } else {
-            for (i = 0; i < myVector->n; i++) {
-                if ((stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    numData++;
-                }
-            }
-        }
-    } else {
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    numData++;
-                }
-            }
-        } else {
-            numData = myVector->n;
-        }
-    }
-
-    psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, numData);
-    return (numData);
-}
-
-/******************************************************************************
-p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the
-median of the input vector.  Returns true on success (including if there were
-no valid input vector elements).
- *****************************************************************************/
-bool p_psVectorSampleMedian(const psVector* myVector,
-                            const psVector* maskVector,
-                            psU32 maskVal,
-                            psStats* stats)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    psVector* unsortedVector = NULL;    // Temporary vector
-    psVector* sortedVector = NULL;      // Temporary vector
-    psS32 i = 0;                        // Loop index variable
-    psS32 count = 0;
-    psS32 nValues = 0;
-
-    // Determine how many data points fit inside this min/max range
-    // and are not masked, if the maskVector is not NULL.
-    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
-
-    if (nValues <= 0) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian(): no valid elements in input vector.\n");
-        return(true);
-        stats->sampleMedian = NAN;
-    }
 
     // Allocate temporary vectors for the data.
-    unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
-    unsortedVector->n = nValues;
+    psVector* vector = psVectorAlloc(myVector->n, PS_TYPE_F32); // Vector containing the unmasked values
+    long count = 0;                     // Number of valid entries
 
     // Determine if we must only use data points within a min/max range.
@@ -575 +451 @@
         // Store all non-masked data points within the min/max range
         // into the temporary vectors.
-        count = 0;
         if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
+            for (long i = 0; i < myVector->n; i++) {
                 if (!(maskVal & maskVector->data.U8[i]) &&
                         (stats->min <= myVector->data.F32[i]) &&
                         (myVector->data.F32[i] <= stats->max)) {
-                    unsortedVector->data.F32[count] = myVector->data.F32[i];
+                    vector->data.F32[count] = myVector->data.F32[i];
                     count++;
                 }
             }
         } else {
-            for (i = 0; i < myVector->n; i++) {
+            for (long i = 0; i < myVector->n; i++) {
                 if ((stats->min <= myVector->data.F32[i]) &&
                         (myVector->data.F32[i] <= stats->max)) {
-                    unsortedVector->data.F32[count] = myVector->data.F32[i];
+                    vector->data.F32[count] = myVector->data.F32[i];
                     count++;
                 }
@@ -596 +471 @@
     } else {
         // Store all non-masked data points into the temporary vectors.
-        count = 0;
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
+        if (maskVector) {
+            for (long i = 0; i < myVector->n; i++) {
                 if (!(maskVal & maskVector->data.U8[i])) {
-                    unsortedVector->data.F32[count] = myVector->data.F32[i];
+                    vector->data.F32[count] = myVector->data.F32[i];
                     count++;
                 }
             }
         } else {
-            for (i = 0; i < myVector->n; i++) {
-                unsortedVector->data.F32[i] = myVector->data.F32[i];
-            }
-        }
-    }
-    // Sort the temporary vectors.
-    sortedVector = psVectorSort(sortedVector, unsortedVector);
-    if (sortedVector == NULL) {
+            vector = psVectorCopy(vector, myVector, PS_TYPE_F32);
+            count = myVector->n;
+        }
+    }
+    vector->n = count;
+    if (count == 0) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No valid data in input vector.\n");
+        psFree(vector);
+        stats->sampleMedian = NAN;
+        return false;
+    }
+
+    // Sort the temporary vector.
+    vector = psVectorSort(vector, vector); // Sort in-place (since it's a copy, it's OK)
+    if (!vector) {
         psError(PS_ERR_UNEXPECTED_NULL,
                 false,
                 PS_ERRORTEXT_psStats_STATS_SAMPLE_MEDIAN_SORT_PROBLEM);
         psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
-        psFree(unsortedVector);
-        return(false);
+        psFree(vector);
+        return false;
     }
 
     // Calculate the median exactly.  Use the average if the number of samples
     // is even.
-    if (0 == (nValues % 2)) {
-        stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues / 2) - 1] +
-                                     sortedVector->data.F32[nValues / 2]);
+    if (count % 2 == 0) {
+        stats->sampleMedian = 0.5 * (vector->data.F32[(count / 2) - 1] +
+                                     vector->data.F32[count / 2]);
     } else {
-        stats->sampleMedian = sortedVector->data.F32[nValues / 2];
-    }
+        stats->sampleMedian = vector->data.F32[count / 2];
+    }
+
+    // Calculate the quartile points exactly.
+    stats->sampleUQ = vector->data.F32[3 * (count / 4)];
+    stats->sampleLQ = vector->data.F32[count / 4];
 
     // Free the temporary data structures.
-    psFree(unsortedVector);
-    psFree(sortedVector);
+    psFree(vector);
 
     // Return "true" on success.
     psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
-    return(true);
+    return true;
 }
 
 /******************************************************************************
-p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input
+vectorSmoothHistGaussian(): This routine smoothes the data in the input
 robustHistogram with a Gaussian of width sigma.  It returns a psVector of the
 smoothed data.
  *****************************************************************************/
-psVector *p_psVectorSmoothHistGaussian(
-    psHistogram *histogram,
-    psF32 sigma)
+psVector *vectorSmoothHistGaussian(psHistogram *histogram,
+                                   psF32 sigma)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
@@ -657 +540 @@
     }
 
-    psS32 numBins = histogram->nums->n;
-    psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32);
+    long numBins = histogram->nums->n;  // Number of histogram bins
+    psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); // Smoothed version of histogram bins
     smooth->n = smooth->nalloc;
-    psS32 jMin = 0;
-    psS32 jMax = 0;
-
-    if (histogram->uniform == false) {
+    const psVector *bounds = histogram->bounds; // The bounds for the histogram bins
+
+    if (!histogram->uniform) {
         //
         // We get here if the histogram is non-uniform.  This code is not tested.
         // However, it is also not used anywhere, yet.
         //
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSmoothHistGaussian() on non-uniform histograms has not been tested or used.\n");
-
-        for (psS32 i = 0; i < numBins; i++) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSmoothHistGaussian() on non-uniform histograms has not been tested or used.\n");
+
+        for (long i = 0; i < numBins; i++) {
             // Determine the midpoint of bin i.
-            psS32 iMid = PS_BIN_MIDPOINT(histogram, i);
+            float iMid = PS_BIN_MIDPOINT(histogram, i);
 
             //
@@ -678 +560 @@
             // range of data values surrounding iMid.  The range is of size:
             // 2*PS_GAUSS_WIDTH*sigma
+            long jMin, jMax;
             psF32 x = iMid - (PS_GAUSS_WIDTH * sigma);
-            jMin = i;
-            while (jMin > 0) {
-                if ((histogram->bounds->data.F32[jMin] <= x) &&
-                        (histogram->bounds->data.F32[jMin+1] >= x)) {
-                    break;
-                }
-                jMin--;
-            }
-
+            for (jMin = i; jMin > 0 && bounds->data.F32[jMin] > x; jMin--)
+                ;
             x = iMid + (PS_GAUSS_WIDTH * sigma);
-            jMax = i;
-            while (jMax < (histogram->bounds->n - 1)) {
-                if ((histogram->bounds->data.F32[jMax] <= x) &&
-                        (histogram->bounds->data.F32[jMax+1] >= x)) {
-                    break;
-                }
-                jMax++;
-            }
-
+            for (jMax = i; jMax < bounds->n - 1 && bounds->data.F32[jMax + 1] > x; jMax++)
+                ;
 
             //
@@ -703 +572 @@
             //
             smooth->data.F32[i] = 0.0;
-            for (psS32 j = jMin ; j <= jMax ; j++) {
-                psS32 jMid = PS_BIN_MIDPOINT(histogram, j);
+            for (long j = jMin ; j <= jMax ; j++) {
+                float jMid = PS_BIN_MIDPOINT(histogram, j);
                 smooth->data.F32[i] +=
                     histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
@@ -713 +582 @@
         // We get here if the histogram is uniform.
         //
-        for (psS32 i = 0; i < numBins; i++) {
-            psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0];
-            psS32 gaussWidth = (psS32) ((PS_GAUSS_WIDTH * sigma) / binSize);
+        for (long i = 0; i < numBins; i++) {
+            psF32 binSize = bounds->data.F32[1] - bounds->data.F32[0];
+            psS32 gaussWidth = ((PS_GAUSS_WIDTH * sigma) / binSize);
 
             //
@@ -723 +592 @@
             // took way too long.
             //
+            #if 0
+
             gaussWidth = 10;
+            #endif
             //
             // We determine the bin numbers (jMin:jMax) corresponding to a
@@ -729 +601 @@
             // 2*PS_GAUSS_WIDTH*sigma
             //
-            psS32 jMin = i - gaussWidth;
-            if (jMin < 0 ) {
-                jMin = 0;
-            }
-            psS32 jMax = i + gaussWidth;
-            if (jMax > (histogram->bounds->n - 1)) {
-                jMax = (histogram->bounds->n - 1);
-            }
+            psS32 jMin = PS_MAX(i - gaussWidth, 0);
+            psS32 jMax = PS_MIN(i + gaussWidth, bounds->n - 1);
 
             //
@@ -742 +608 @@
             //
             smooth->data.F32[i] = 0.0;
-            psS32 iMid = PS_BIN_MIDPOINT(histogram, i);
-            for (psS32 j = jMin ; j <= jMax ; j++) {
-                psS32 jMid = PS_BIN_MIDPOINT(histogram, j);
+            float iMid = PS_BIN_MIDPOINT(histogram, i);
+            for (long j = jMin ; j <= jMax ; j++) {
+                float jMid = PS_BIN_MIDPOINT(histogram, j);
                 smooth->data.F32[i] +=
                     histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
@@ -757 +623 @@
     return(smooth);
 }
+
 /******************************************************************************
-p_psVectorSampleQuartiles(myVector, maskVector, maskVal, statsalculates
-the upper and/or lower quartiles of the input vector.
-Inputs
-    myVector
-    maskVector
-    maskVal
-    stats
-Returns
-    NULL
- *****************************************************************************/
-bool p_psVectorSampleQuartiles(const psVector* myVector,
-                               const psVector* maskVector,
-                               psU32 maskVal,
-                               psStats* stats)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    psVector* unsortedVector = NULL;
-    psVector* sortedVector = NULL;
-    psS32 i = 0;
-    psS32 count = 0;
-    psS32 nValues = 0;
-
-    // Determine how many data points fit inside this min/max range
-    // and are not masked, if the maskVector is not NULL.
-    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
-
-    // Allocate temporary vectors for the data.
-    unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
-    unsortedVector->n = unsortedVector->nalloc;
-
-    if (stats->options & PS_STAT_USE_RANGE) {
-        // Store all non-masked data points within the min/max range
-        // into the temporary vectors.
-        count = 0;
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i]) &&
-                        (stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) {
-                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
-                }
-            }
-        } else {
-            for (i = 0; i < myVector->n; i++) {
-                if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) {
-                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
-                }
-            }
-        }
-    } else {
-        // Store all non-masked data points into the temporary vectors.
-        count = 0;
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
-                }
-            }
-        } else {
-            for (i = 0; i < myVector->n; i++) {
-                unsortedVector->data.F32[i] = myVector->data.F32[i];
-            }
-        }
-    }
-
-    // Sort the temporary vectors.
-    sortedVector = psVectorSort(sortedVector, unsortedVector);
-    if (sortedVector == NULL) {
-        psError(PS_ERR_UNEXPECTED_NULL,
-                false,
-                PS_ERRORTEXT_psStats_STATS_SAMPLE_MEDIAN_SORT_PROBLEM);
-        psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
-        return(false);
-    }
-
-    // Calculate the quartile points exactly.
-    stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
-    stats->sampleLQ = sortedVector->data.F32[nValues / 4];
-
-    // Free the temporary data structures.
-    psFree(unsortedVector);
-    psFree(sortedVector);
-    psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
-    return(true);
-}
-
-/******************************************************************************
-p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
+vectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
 stdev of the input vector.
 Inputs
@@ -855 +636 @@
  
  *****************************************************************************/
-void p_psVectorSampleStdev(const psVector* myVector,
-                           const psVector* errors,
-                           const psVector* maskVector,
-                           psU32 maskVal,
-                           psStats* stats)
+static bool vectorSampleStdev(const psVector* myVector,
+                              const psVector* errors,
+                              const psVector* maskVector,
+                              psMaskType maskVal,
+                              psStats* stats)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    psS32 i = 0;                  // Loop index variable
-    psS32 countInt = 0;           // # of data points being used
-    psF32 countFloat = 0.0;     // # of data points being used
-    psF32 mean = 0.0;           // The mean
-    psF32 diff = 0.0;           // Used in calculating stdev
-    psF32 sumSquares = 0.0;     // temporary variable
-    psF32 sumDiffs = 0.0;       // temporary variable
-    psF32 errorDivisor = 0.0f;
 
     // This procedure requires the mean.  If it has not been already
-    // calculated, then call p_psVectorSampleMean()
-    if (0 != isnan(stats->sampleMean)) {
-        p_psVectorSampleMean(myVector, errors, maskVector, maskVal, stats);
+    // calculated, then call vectorSampleMean()
+    if (isnan(stats->sampleMean)) {
+        vectorSampleMean(myVector, errors, maskVector, maskVal, stats);
     }
     // If the mean is NAN, then generate a warning and set the stdev to NAN.
-    if (0 != isnan(stats->sampleMean)) {
+    if (isnan(stats->sampleMean)) {
         stats->sampleStdev = NAN;
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): p_psVectorSampleMean() reported a NAN mean.\n");
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleStdev(): vectorSampleMean() reported a NAN mean.\n");
         psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-        return;
-    }
-
-    mean = stats->sampleMean;
+        return false;
+    }
+
+    // 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
     if (stats->options & PS_STAT_USE_RANGE) {
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
+        if (maskVector) {
+            for (long i = 0; i < myVector->n; i++) {
                 if (!(maskVal & maskVector->data.U8[i]) &&
                         (stats->min <= myVector->data.F32[i]) &&
                         (myVector->data.F32[i] <= stats->max)) {
-                    diff = myVector->data.F32[i] - mean;
-                    sumSquares += (diff * diff);
+                    double diff = myVector->data.F32[i] - mean;
+                    sumSquares += PS_SQR(diff);
                     sumDiffs += diff;
-                    countInt++;
-                    if (errors != NULL) {
-                        errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i]));
+                    count++;
+                    if (errors) {
+                        errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
                     }
                 }
             }
         } else {
-            for (i = 0; i < myVector->n; i++) {
+            for (long i = 0; i < myVector->n; i++) {
                 if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) {
-                    diff = myVector->data.F32[i] - mean;
-                    sumSquares += (diff * diff);
+                    double diff = myVector->data.F32[i] - mean;
+                    sumSquares += PS_SQR(diff);
                     sumDiffs += diff;
-                    countInt++;
-                    if (errors != NULL) {
-                        errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i]));
+                    count++;
+                    if (errors) {
+                        errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
                     }
                 }
             }
-            countInt = myVector->n;
+            count = myVector->n;
         }
     } else {
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
+        if (maskVector) {
+            for (long i = 0; i < myVector->n; i++) {
                 if (!(maskVal & maskVector->data.U8[i])) {
-                    diff = myVector->data.F32[i] - mean;
-                    sumSquares += (diff * diff);
+                    double diff = myVector->data.F32[i] - mean;
+                    sumSquares += PS_SQR(diff);
                     sumDiffs += diff;
-                    countInt++;
-                    if (errors != NULL) {
-                        errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i]));
+                    count++;
+                    if (errors) {
+                        errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
                     }
                 }
             }
         } else {
-            for (i = 0; i < myVector->n; i++) {
-                diff = myVector->data.F32[i] - mean;
-                sumSquares += (diff * diff);
+            for (long i = 0; i < myVector->n; i++) {
+                double diff = myVector->data.F32[i] - mean;
+                sumSquares += PS_SQR(diff);
                 sumDiffs += diff;
-                countInt++;
-            }
-            countInt = myVector->n;
-            if (errors != NULL) {
-                errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i]));
-            }
-        }
-    }
-
-    if (countInt == 0) {
+                count++;
+                if (errors) {
+                    errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
+                }
+            }
+            count = myVector->n;
+        }
+    }
+
+    if (count == 0) {
         stats->sampleStdev = NAN;
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): no valid psVector elements (%d).  Setting stats->sampleStdev = NAN.\n", countInt);
-    } else if (countInt == 1) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleStdev(): no valid psVector elements (%d).  Setting stats->sampleStdev = NAN.\n", count);
+        return false;
+    }
+    if (count == 1) {
         stats->sampleStdev = 0.0;
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d).  Setting stats->sampleStdev = 0.0.\n", countInt);
-    } else {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleStdev(): only one valid psVector elements (%d).  Setting stats->sampleStdev = 0.0.\n", count);
+        return false;
+    }
+
+    if (errors) {
+        #if 1
         // XXX: The ADD specifies this as the definition of the standard
         // deviation if the errors are known.  Verify this with IfA: none of
         // the data points in the vector are used.  Verify that the masks and
         // data ranges are used correctly.
-        if (errors != NULL) {
-            if (0) {
-                stats->sampleStdev = (1.0 / sqrtf(errorDivisor));
-            } else {
-                countFloat = (psF32)countInt;
-                stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
-            }
-        } else {
-            countFloat = (psF32)countInt;
-            stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
-
-        }
+        stats->sampleStdev = (1.0 / sqrtf(errorDivisor));
+        #else
+
+        stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) /
+                                  (float)(count - 1));
+        #endif
+
+    } else {
+        stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) /
+                                  (float)(count - 1));
+
     }
     psTrace(__func__, 4, "---- %s() end ----\n", __func__);
+
+    return true;
 }
 
 /******************************************************************************
-p_psVectorClippedStats(myVector, errors, maskVector, maskVal, stats): calculates the
+vectorClippedStats(myVector, errors, maskVector, maskVal, stats): calculates the
 clipped stats (mean or stdev) of the input vector.
  
@@ -978 +763 @@
     stats
 Returns
-    0 for success.
-    -1: error
-    -2: warning
+    true for success; false otherwise
  *****************************************************************************/
-psS32 p_psVectorClippedStats(const psVector* myVector,
-                             const psVector* errors,
-                             const psVector* maskVector,
-                             psU32 maskVal,
-                             psStats* stats)
+static bool vectorClippedStats(const psVector* myVector,
+                               const psVector* errors,
+                               psVector* maskVector,
+                               psMaskType maskVal,
+                               psStats* stats
+                              )
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     psTrace(__func__, 4, "Trace level is %d\n", psTraceGetLevel(__func__));
-    psF32 clippedMean = 0.0;
-    psF32 clippedStdev = 0.0;
-    psVector* tmpMask = NULL;
-    static psStats *statsTmp = NULL;
-    psS32 rc = 0;
 
     // Ensure that stats->clipIter is within the proper range.
@@ -1008 +787 @@
     // Allocate a psStats structure for calculating the mean, median, and
     // stdev.
-    if (statsTmp == NULL) {
-        statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
-        p_psMemSetPersistent(statsTmp, true);
+    psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+
+    // We've already copied the mask vector, so it's safe to simply increment the reference counter.
+    psVector *tmpMask = NULL;
+    if (maskVector) {
+        tmpMask = psMemIncrRefCounter(maskVector);
     } else {
-        // Initialize structure if already allocated
-        statsTmp->sampleMean = NAN;
-        statsTmp->sampleMedian = NAN;
-        statsTmp->sampleStdev = NAN;
-        statsTmp->sampleUQ = NAN;
-        statsTmp->sampleLQ = NAN;
-        statsTmp->robustMedian = NAN;
-        statsTmp->robustStdev = NAN;
-        statsTmp->robustUQ = NAN;
-        statsTmp->robustLQ = NAN;
-        statsTmp->robustN50 = -1;
-        statsTmp->fittedMean = NAN;
-        statsTmp->fittedStdev = NAN;
-        statsTmp->fittedNfit = -1;
-        statsTmp->clippedMean = NAN;
-        statsTmp->clippedStdev = NAN;
-        statsTmp->clippedNvalues = -1;
-        statsTmp->clipSigma = 3.0;
-        statsTmp->clipIter = 3;
-        statsTmp->min = NAN;
-        statsTmp->max = NAN;
-        statsTmp->binsize = NAN;
-        statsTmp->nSubsample = 100000;
-    }
-
-    // We allocate a temporary mask vector since during the iterative
-    // steps that follow, we will be masking off additional data points.
-    // However, we do no want to modify the original mask vector.
-    tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
-    tmpMask->n = tmpMask->nalloc;
-
-    // If we were called with a mask vector, then initialize the temporary
-    // mask vector with those values.  Otherwise, initialize to zero.
-    if (maskVector != NULL) {
-        for (psS32 i = 0; i < tmpMask->n; i++) {
-            tmpMask->data.U8[i] = maskVector->data.U8[i];
-        }
-    } else {
-        for (psS32 i = 0; i < tmpMask->n; i++) {
-            tmpMask->data.U8[i] = 0;
-        }
+        tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
+        tmpMask->n = tmpMask->nalloc;
+        psVectorInit(tmpMask, 0);
     }
 
     // 1. Compute the sample median.
-    p_psVectorSampleMedian(myVector, maskVector, maskVal, statsTmp);
+    vectorSampleMedian(myVector, tmpMask, maskVal, statsTmp);
     if (isnan(statsTmp->sampleMedian)) {
-        psLogMsg(__func__, PS_LOG_WARN, "Call to p_psVectorSampleMedian returned NAN\n");
+        psLogMsg(__func__, PS_LOG_WARN, "Call to vectorSampleMedian returned NAN\n");
         stats->clippedMean = NAN;
         stats->clippedStdev = NAN;
-        psTrace(__func__, 4, "---- %s(-2) end ----\n", __func__);
-        return(-2);
+        psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
+        psFree(tmpMask);
+        psFree(statsTmp);
+        return false;
     }
     psTrace(__func__, 6, "The initial sample median is %f\n", statsTmp->sampleMedian);
 
     // 2. Compute the sample standard deviation.
-    p_psVectorSampleStdev(myVector, errors, maskVector, maskVal, statsTmp);
+    vectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp);
     if (isnan(statsTmp->sampleStdev)) {
-        psLogMsg(__func__, PS_LOG_WARN, "Call to p_psVectorSampleStdev returned NAN\n");
+        psLogMsg(__func__, PS_LOG_WARN, "Call to vectorSampleStdev returned NAN\n");
         stats->clippedMean = NAN;
         stats->clippedStdev = NAN;
-        psTrace(__func__, 4, "---- %s(-2) end ----\n", __func__);
-        return(-2);
+        psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
+        psFree(tmpMask);
+        psFree(statsTmp);
+        return false;
     }
     psTrace(__func__, 6, "The initial sample stdev is %f\n", statsTmp->sampleStdev);
 
     // 3. Use the sample median as the first estimator of the mean X.
-    clippedMean = statsTmp->sampleMedian;
+    psF32 clippedMean = statsTmp->sampleMedian;
 
     // 4. Use the sample stdev as the first estimator of the mean stdev.
-    clippedStdev = statsTmp->sampleStdev;
+    psF32 clippedStdev = statsTmp->sampleStdev;
 
     // 5. Repeat N (stats->clipIter) times:
-    for (psS32 iter = 0; iter < stats->clipIter; iter++) {
+    long numClipped = LONG_MAX;         // Number of values clipped in this iteration
+    for (int iter = 0; iter < stats->clipIter && numClipped > 0; iter++) {
+        numClipped = 0;
         psTrace(__func__, 6, "------------ Iteration %d ------------\n", iter);
         // a) Exclude all values x_i for which |x_i - x| > K * stdev
-        if (errors != NULL) {
-            for (psS32 j = 0; j < myVector->n; j++) {
-                if (fabs(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * errors->data.F32[j])) {
+        if (errors) {
+            for (long j = 0; j < myVector->n; j++) {
+                if (fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * errors->data.F32[j])) {
                     tmpMask->data.U8[j] = 0xff;
-                }
-            }
-        } else {
-            for (psS32 j = 0; j < myVector->n; j++) {
-                if (fabs(myVector->data.F32[j] - clippedMean) >
-                        (stats->clipSigma * clippedStdev)) {
+                    psTrace(__func__, 10, "Clipped %d: %f +/- %f\n", j,
+                            myVector->data.F32[j], errors->data.F32[j]);
+                    numClipped++;
+                }
+            }
+        } else {
+            for (long j = 0; j < myVector->n; j++) {
+                if (fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) {
                     tmpMask->data.U8[j] = 0xff;
+                    psTrace(__func__, 10, "Clipped %d: %f\n", j, myVector->data.F32[j]);
+                    numClipped++;
                 }
             }
@@ -1103 +857 @@
 
         // b) compute new mean and stdev
-        p_psVectorSampleMean(myVector, errors, tmpMask, 0xff, statsTmp);
-        p_psVectorSampleStdev(myVector, errors, tmpMask, 0xff, statsTmp);
+        vectorSampleMean(myVector, errors, tmpMask, maskVal, statsTmp);
+        vectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp);
         psTrace(__func__, 6, "The new sample mean is %f\n", statsTmp->sampleMean);
         psTrace(__func__, 6, "The new sample stdev is %f\n", statsTmp->sampleStdev);
@@ -1113 +867 @@
             // Exit loop.
             iter = stats->clipIter;
-            psError(PS_ERR_UNKNOWN, true, "p_psVectorSampleMean() or p_psVectorSampleStdev() returned a NAN.\n");
+            psError(PS_ERR_UNKNOWN, true, "vectorSampleMean() or vectorSampleStdev() returned a NAN.\n");
             psFree(tmpMask);
-            return(-1);
+            return false;
         } else {
             clippedMean = statsTmp->sampleMean;
@@ -1122 +876 @@
     }
 
-    // 7. The last calcuated value of x is the cliped mean.
+    // 7. The last calcuated value of x is the clipped mean.
     if (stats->options & PS_STAT_CLIPPED_MEAN) {
         stats->clippedMean = clippedMean;
         psTrace(__func__, 6, "The final clipped mean is %f\n", clippedMean);
     }
-    // 8. The last calcuated value of stdev is the cliped stdev.
+    // 8. The last calcuated value of stdev is the clipped stdev.
     if (stats->options & PS_STAT_CLIPPED_STDEV) {
         stats->clippedStdev = clippedStdev;
@@ -1134 +888 @@
 
     psFree(tmpMask);
-    psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, rc);
-    return(rc);
+    psFree(statsTmp);
+    psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
+    return true;
 }
 
-static psF32 QuadraticInverse(
-    psF32 a,
-    psF32 b,
-    psF32 c,
-    psF32 y,
-    psF32 xLo,
-    psF32 xHi)
+static psF32 QuadraticInverse(psF32 a,
+                              psF32 b,
+                              psF32 c,
+                              psF32 y,
+                              psF32 xLo,
+                              psF32 xHi
+                             )
 {
     psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a));
@@ -1151 +906 @@
     psF64 x2 = -b/(2.0*a) - tmp;
 
-    if (0) {
-        psF64 y1 = (a * x1 * x1) + (b * x1) + c;
-        psF64 y2 = (a * x2 * x2) + (b * x2) + c;
-        printf("QuadraticInverse: %fx^2 + %fx + %f\n", a, b, c);
-        printf("QuadraticInverse: y is %f\n", y);
-        printf("QuadraticInverse: (x1, x2) is (%f %f)\n", x1, x2);
-        printf("QuadraticInverse: (y1, y2) is (%f %f)\n", y1, y2);
-    }
-
-    if ((xLo <= x1) && (x1 <= xHi)) {
-        return(x1);
-    } else if ((xLo <= x2) && (x2 <= xHi)) {
-        return(x2);
-    } else {
-        return(0.5 * (xLo + xHi));
-    }
+    #if 0
+
+    psF64 y1 = (a * x1 * x1) + (b * x1) + c;
+    psF64 y2 = (a * x2 * x2) + (b * x2) + c;
+    printf("QuadraticInverse: %fx^2 + %fx + %f\n", a, b, c);
+    printf("QuadraticInverse: y is %f\n", y);
+    printf("QuadraticInverse: (x1, x2) is (%f %f)\n", x1, x2);
+    printf("QuadraticInverse: (y1, y2) is (%f %f)\n", y1, y2);
+    #endif
+
+    if (xLo <= x1 && x1 <= xHi) {
+        return x1;
+    }
+    if (xLo <= x2 && x2 <= xHi) {
+        return x2;
+    }
+    return 0.5 * (xLo + xHi);
 }
 
@@ -1180 +936 @@
  
 *****************************************************************************/
-psF32 fitQuadraticSearchForYThenReturnX(
-    psVector *xVec,
-    psVector *yVec,
-    psS32 binNum,
-    psF32 yVal)
+static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec,
+        psVector *yVec,
+        psS32 binNum,
+        psF32 yVal
+                                              )
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
@@ -1235 +991 @@
         if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) ||
                 ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
-            psError(PS_ERR_UNKNOWN, true, "This routine must be called with monotonically increasing or decreasing data points.\n");
+            psError(PS_ERR_UNKNOWN, true,
+                    "This routine must be called with monotonically increasing or decreasing data points.\n");
             psFree(x);
             psFree(y);
             psTrace(__func__, 5, "---- %s() end ----\n", __func__);
-            return(NAN);
+            return NAN;
         }
 
@@ -1251 +1008 @@
             psFree(y);
             psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__);
-            return(NAN);
+            return NAN;
         }
         psTrace(__func__, 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
@@ -1262 +1019 @@
 
         psTrace(__func__, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
-        tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[2]);
+        tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal,
+                                    x->data.F64[0], x->data.F64[2]);
         psFree(myPoly);
 
@@ -1295 +1053 @@
 
     psTrace(__func__, 5, "---- %s(%f) end ----\n", __func__, tmpFloat);
-    return(tmpFloat);
+    return tmpFloat;
 }
 
@@ -1301 +1059 @@
 NOTE: We assume unnormalized gaussians.
  *****************************************************************************/
-psF32 psMinimizeLMChi2Gauss1D(
-    psVector *deriv,
-    const psVector *params,
-    const psVector *coords)
+static psF32 minimizeLMChi2Gauss1D(psVector *deriv,
+                                   const psVector *params,
+                                   const psVector *coords
+                                  )
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
@@ -1318 +1076 @@
     psF32 stdev = params->data.F32[1];
 
-    if (deriv == NULL) {
-        deriv = psVectorAlloc(2, PS_TYPE_F32);
-        deriv->n = 2;
-    } else {
+    psF32 gauss = psGaussian(x, mean, stdev, false);
+    if (deriv) {
         PS_ASSERT_VECTOR_SIZE(deriv, 2, NAN);
         PS_ASSERT_VECTOR_TYPE(deriv, PS_TYPE_F32, NAN);
-    }
-
-    psF32 tmp = (x - mean) * psGaussian(x, mean, stdev, false);
-    deriv->data.F32[0] = tmp / (stdev * stdev);
-    tmp = (x - mean) * (x - mean) * psGaussian(x, mean, stdev, 0);
-    deriv->data.F32[1] = tmp / (stdev * stdev * stdev);
+        psF32 tmp = (x - mean) * gauss;
+        deriv->data.F32[0] = tmp / (stdev * stdev);
+        deriv->data.F32[1] = (x - mean) * tmp / (stdev * stdev * stdev);
+    }
+
 
     psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-    return(psGaussian(x, mean, stdev, false));
+    return gauss;
 }
 
 /******************************************************************************
-p_psVectorRobustStats(myVector, maskVector, maskVal, stats): This is the new
+vectorRobustStats(myVector, maskVector, maskVal, stats): This is the new
 version of the robust stats routine.
  
@@ -1354 +1109 @@
 *****************************************************************************/
 #define INITIAL_NUM_BINS 1000.0
-psS32 p_psVectorRobustStats(const psVector* myVector,
-                            const psVector* errors,
-                            const psVector* maskVector,
-                            psU32 maskVal,
-                            psStats* stats)
+static bool vectorRobustStats(const psVector* myVector,
+                              const psVector* errors,
+                              psVector* maskVector,
+                              psMaskType maskVal,
+                              psStats* stats)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
@@ -1364 +1119 @@
         PS_VECTOR_PRINT_F32(myVector);
     }
-    psHistogram *robustHistogram = NULL;
-    psHistogram *cumulativeRobustHistogram = NULL;
-    psS32 numBins = 0;
-    psScalar tmpScalar;
+
+    psScalar tmpScalar;                 // Temporary scalar variable, for p_psVectorBinDisect
     tmpScalar.type.type = PS_TYPE_F32;
-    psS32 totalDataPoints = 0;
-    psS32 rc = 0;
-    psS32 rcBool = false;
-    psVector *tmpMaskVec = psVectorAlloc(myVector->n, PS_TYPE_U8);
-    tmpMaskVec->n = tmpMaskVec->nalloc;
-    if (maskVector != NULL) {
-        for (psS32 i = 0 ; i < myVector->n ; i++) {
-            tmpMaskVec->data.U8[i] = maskVector->data.U8[i];
-        }
+
+    // We need a mask for these statistics, so if the one we're given is NULL, make our own.
+    psVector *mask = NULL;              // The actual mask we will use
+    if (maskVector) {
+        mask = psMemIncrRefCounter(maskVector);
     } else {
-        for (psS32 i = 0 ; i < myVector->n ; i++) {
-            tmpMaskVec->data.U8[i] = 0;
-        }
-    }
-    psS32 iterate = 1;
-    psF32 sigma;
-    psStats *tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX);
-
-    while (iterate > 0) {
-        psTrace(__func__, 6, "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n", iterate);
-        psF32 binSize = 0.0;
+        mask = psVectorAlloc(myVector->n, PS_TYPE_U8);
+        mask->n = myVector->n;
+        psVectorInit(mask, 0);
+    }
+
+    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
+
+    // Iterate to get the best bin size
+    for (int iterate = 1; iterate > 0; iterate++) {
+        psTrace(__func__, 6,
+                "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n",
+                iterate);
+
+        // Get the minimum and maximum values
+        int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of valid values
+        min = statsMinMax->min;
+        max = statsMinMax->max;
+        if (numValid == 0 || isnan(min) || isnan(max)) {
+            // Data range calculation failed
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
+            psFree(statsMinMax);
+            psFree(mask);
+            psTrace(__func__, 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+        psTrace(__func__, 6, "Data min/max is (%.2f, %.2f)\n", min, max);
+
+        // If all data points have the same value, then we set the appropiate members of stats and return.
+        if (fabs(max - min) <= FLT_EPSILON) {
+            psTrace(__func__, 5, "All data points have the same value: %f.\n", min);
+            if (stats->options & PS_STAT_ROBUST_MEDIAN) {
+                stats->robustMedian = min;
+            }
+            if (stats->options & PS_STAT_ROBUST_QUARTILE) {
+                stats->robustUQ = min;
+                stats->robustLQ = min;
+            }
+            stats->robustN50 = numValid;
+            psFree(statsMinMax);
+            psFree(mask);
+
+            psTrace(__func__, 4, "---- %s(0) end  ----\n", __func__);
+            return false;
+        }
+
+        float binSize = 0.0;            // Size of bins for histogram
         if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE)) {
             // Set initial bin size to the specified value.
@@ -1395 +1184 @@
             psTrace(__func__, 6, "Setting initial robust bin size to %.2f\n", binSize);
         } else {
-            // Determine the bin size of the robust histogram.  This is done
-            // by computing the total range of data values and dividing by 1000.0.
-
-            rc = p_psVectorMin(myVector, tmpMaskVec, 1, tmpStatsMinMax);
-            rc|= p_psVectorMax(myVector, tmpMaskVec, 1, tmpStatsMinMax);
-            if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
-                psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
-                psFree(tmpStatsMinMax);
-                psFree(tmpMaskVec);
-                psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-                return(1);
-            }
-            psTrace(__func__, 6, "Data min/max is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
-            binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / INITIAL_NUM_BINS;
+            // Determine the bin size of the robust histogram, using the pre-defined number of bins
+            binSize = (max - min) / INITIAL_NUM_BINS;
         }
         psTrace(__func__, 6, "Initial robust bin size is %.2f\n", binSize);
 
-        //
-        // If all data points have the same value, then we set the appropiate
-        // members of stats and return.
-        //
-        if (fabs(tmpStatsMinMax->max - tmpStatsMinMax->min) <= FLT_EPSILON) {
-            psTrace(__func__, 5, "All data points have the same value.\n", binSize);
-            if (stats->options & PS_STAT_ROBUST_MEDIAN) {
-                stats->robustMedian = tmpStatsMinMax->min;
-            }
-            if (stats->options & PS_STAT_ROBUST_QUARTILE) {
-                stats->robustUQ = tmpStatsMinMax->min;
-                stats->robustLQ = tmpStatsMinMax->min;
-            }
-            // XXX: Set these to the number of unmasked data points?
-            stats->robustN50 = 0;
-            psFree(tmpStatsMinMax);
-            psFree(tmpMaskVec);
-
-            psTrace(__func__, 4, "---- %s(0) end  ----\n", __func__);
-            return(0);
-        }
-
-        //
-        // ADD: Step 0.
-        // Construct the histogram with the specified bin size.
-        //
-        // NOTE: we can not specify the bin size precisely since the argument
-        // to psHistogramAlloc() is the number of bins, not the binSize.
-        // If we get here, we know that binSize != 0.0.
-        //
-        numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / binSize);
+        // ADD step 0: Construct the histogram with the specified bin size.  NOTE: we can not specify the bin
+        // size precisely since the argument to psHistogramAlloc() is the number of bins, not the binSize.  If
+        // we get here, we know that binSize != 0.0.
+        long numBins = (max - min) / binSize; // Number of bins
         psTrace(__func__, 6, "Numbins is %d\n", numBins);
-        psTrace(__func__, 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
-        robustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
-        cumulativeRobustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
-
-        // Populate the histogram array.
-        robustHistogram = psVectorHistogram(robustHistogram, myVector, errors, tmpMaskVec, maskVal);
+        psTrace(__func__, 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
+        // Generate the histogram
+        histogram = psHistogramAlloc(min, max, numBins);
+        histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal);
         if (psTraceGetLevel(__func__) >= 8) {
-            PS_VECTOR_PRINT_F32(robustHistogram->bounds);
-            PS_VECTOR_PRINT_F32(robustHistogram->nums);
-        }
-        //
-        // ADD: Step 1.
-        // Construct the cumulative histogram from the specific histogram
-        //
-        cumulativeRobustHistogram->nums->data.F32[0] = robustHistogram->nums->data.F32[0];
-        for (psS32 i = 1 ; i < robustHistogram->nums->n ; i++) {
-            cumulativeRobustHistogram->nums->data.F32[i] = cumulativeRobustHistogram->nums->data.F32[i-1] +
-                    robustHistogram->nums->data.F32[i];
+            PS_VECTOR_PRINT_F32(histogram->bounds);
+            PS_VECTOR_PRINT_F32(histogram->nums);
+        }
+
+        // ADD step 1: Convert the specific histogram to a cumulative histogram
+        cumulative = psHistogramAlloc(min, max, numBins);
+        cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
+        for (long i = 1; i < histogram->nums->n; i++) {
+            cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
         }
         if (psTraceGetLevel(__func__) >= 8) {
-            PS_VECTOR_PRINT_F32(cumulativeRobustHistogram->bounds);
-            PS_VECTOR_PRINT_F32(cumulativeRobustHistogram->nums);
-        }
-
-        //
-        // ADD: Step 2.
-        // Find the bin which contains the 50% data point.
-        //
-        totalDataPoints = cumulativeRobustHistogram->nums->data.F32[numBins - 1];
+            PS_VECTOR_PRINT_F32(cumulative->bounds);
+            PS_VECTOR_PRINT_F32(cumulative->nums);
+        }
+
+        // ADD step 2: Find the bin which contains the 50% data point.
+        totalDataPoints = cumulative->nums->data.F32[numBins - 1];
         psTrace(__func__, 6, "Total data points is %d\n", totalDataPoints);
-        psS32 binMedian;
-        if (totalDataPoints/2.0 < cumulativeRobustHistogram->nums->data.F32[0]) {
+        long binMedian;
+        if (totalDataPoints/2.0 < cumulative->nums->data.F32[0]) {
             // Special case: the median is in the first bin.  Perhaps we should recode this.
             // XXX: Must try a special case where the median in in the last bin.
@@ -1481 +1224 @@
         } else {
             tmpScalar.data.F32 = totalDataPoints/2.0;
-            binMedian = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
+            binMedian = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
             if (binMedian < 0) {
-                psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50 precent data point (%d).\n", binMedian);
-                psFree(tmpStatsMinMax);
-                psFree(robustHistogram);
-                psFree(cumulativeRobustHistogram);
-                psFree(tmpMaskVec);
-                psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-                return(1);
-            }
-        }
-        psTrace(__func__, 6, "The median bin is %d (%.2f to %.2f)\n", binMedian, cumulativeRobustHistogram->bounds->data.F32[binMedian], cumulativeRobustHistogram->bounds->data.F32[binMedian+1]);
-
-        //
-        // ADD: Step 3.
-        // Interpolate to the exact 50% position: this is the robust histogram median.
-        //
-        stats->robustMedian = fitQuadraticSearchForYThenReturnX(
-                                  *(psVector* *)&cumulativeRobustHistogram->bounds,
-                                  *(psVector* *)&cumulativeRobustHistogram->nums,
-                                  binMedian,
-                                  totalDataPoints/2.0);
+                psError(PS_ERR_UNKNOWN, false,
+                        "Failed to calculate the 50 precent data point (%d).\n", binMedian);
+                psFree(statsMinMax);
+                psFree(histogram);
+                psFree(cumulative);
+                psFree(mask);
+                psTrace(__func__, 4, "---- %s(false) end  ----\n", __func__);
+                return false;
+            }
+        }
+        psTrace(__func__, 6, "The median bin is %d (%.2f to %.2f)\n", binMedian,
+                cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]);
+
+        // ADD step 3: Interpolate to the exact 50% position: this is the robust histogram median.
+        stats->robustMedian = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums,
+                              binMedian, totalDataPoints/2.0);
         if (isnan(stats->robustMedian)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to fit a quadratic and calculate the 50-percent position.\n");
-            psFree(tmpStatsMinMax);
-            psFree(robustHistogram);
-            psFree(cumulativeRobustHistogram);
-            psFree(tmpMaskVec);
-            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-            return(1);
+            psError(PS_ERR_UNKNOWN, false,
+                    "Failed to fit a quadratic and calculate the 50-percent position.\n");
+            psFree(statsMinMax);
+            psFree(histogram);
+            psFree(cumulative);
+            psFree(mask);
+            psTrace(__func__, 4, "---- %s(false) end  ----\n", __func__);
+            return false;
         }
         psTrace(__func__, 6, "Current robust median is %f\n", stats->robustMedian);
 
-        //
-        // ADD: Step 4.
-        // Find the bins which contains the 15.8655% and 84.1345% data points.
-        //
-        psS32 binLo = 0;
-        psS32 binHi = 0;
+        // ADD step 4: Find the bins which contains the 15.8655% and 84.1345% data points.
+        long binLo = 0;                 // Bin containing the 15.8655% mark
         tmpScalar.data.F32 = totalDataPoints * 0.158655f;
-        if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
-            binLo = 0;
-        } else {
+        if (tmpScalar.data.F32 > cumulative->nums->data.F32[0]) {
             // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
-            binLo = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
-            if (binLo > cumulativeRobustHistogram->nums->n-1) {
-                binLo = cumulativeRobustHistogram->nums->n-1;
-            }
-        }
+            binLo = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
+            if (binLo > cumulative->nums->n-1) {
+                binLo = cumulative->nums->n-1;
+            }
+        }
+
+        long binHi = 0;                 // Bin containing the 84.1345% mark
         tmpScalar.data.F32 = totalDataPoints * 0.841345f;
-        if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
-            binHi = 0;
-        } else {
+        if (tmpScalar.data.F32 > cumulative->nums->data.F32[0]) {
             // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
-            binHi = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
-            if (binHi > cumulativeRobustHistogram->nums->n-1) {
-                binHi = cumulativeRobustHistogram->nums->n-1;
-            }
-        }
-        psTrace(__func__, 6, "The 15.8655-percent and 84.1345-percent data point bins are (%d, %d).\n",
+            binHi = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
+            if (binHi > cumulative->nums->n-1) {
+                binHi = cumulative->nums->n-1;
+            }
+        }
+        psTrace(__func__, 6, "The 15.8655%% and 84.1345%% data point bins are (%d, %d).\n",
                 binLo, binHi);
-        psTrace(__func__, 6, "binLo midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binLo));
-        psTrace(__func__, 6, "binHi midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binHi));
+        psTrace(__func__, 6, "binLo midpoint is %f\n", PS_BIN_MIDPOINT(cumulative, binLo));
+        psTrace(__func__, 6, "binHi midpoint is %f\n", PS_BIN_MIDPOINT(cumulative, binHi));
 
         if ((binLo < 0) || (binHi < 0)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 15.8655 and 84.1345 percent data point\n");
-            psFree(tmpStatsMinMax);
-            psFree(robustHistogram);
-            psFree(cumulativeRobustHistogram);
-            psFree(tmpMaskVec);
-            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-            return(1);
-        }
-
-        //
-        // ADD: Step 4b.
-        // Interpolate Sigma to find these two positions exactly: these are the 1sigma positions.
-        //
-        // XXX: I am overriding the ADD for now.  The following code implements
-        // the ADD exactly and is having problems fitting a 2nd-order polynomial
-        // to data y-values suchs as (1, 1, 100).  Therefore, I commented out the
-        // code and am doing simply linear interpolation.
-        //
-        psF32 binLoF32;
-        psF32 binHiF32;
-        if (0) {
-            binLoF32 = fitQuadraticSearchForYThenReturnX(
-                           *(psVector* *)&cumulativeRobustHistogram->bounds,
-                           *(psVector* *)&cumulativeRobustHistogram->nums,
-                           binLo,
-                           totalDataPoints * 0.158655f);
-            binHiF32 = fitQuadraticSearchForYThenReturnX(
-                           *(psVector* *)&cumulativeRobustHistogram->bounds,
-                           *(psVector* *)&cumulativeRobustHistogram->nums,
-                           binHi,
-                           totalDataPoints * 0.841345f);
-            psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n",
-                    binLoF32, binHiF32);
-        }
-
-        //
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n");
+            psFree(statsMinMax);
+            psFree(histogram);
+            psFree(cumulative);
+            psFree(mask);
+            psTrace(__func__, 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        // ADD step 4b: Interpolate Sigma to find these two positions exactly: these are the 1sigma positions.
+        #if 0
+        // XXX: I am overriding the ADD for now.  The following code implements the ADD exactly and is having
+        // problems fitting a 2nd-order polynomial to data y-values suchs as (1, 1, 100).  Therefore, I
+        // commented out the code and am doing simply linear interpolation.
+
+        // Value for the 15.8655% mark
+        float binLoF32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, binLo,
+                         totalDataPoints * 0.158655f);
+        // Value for the 84.1345% mark
+        float binHiF32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, binHi,
+                         totalDataPoints * 0.841345f);
+        psTrace(__func__, 6, "The exact 15.8655%% and 84.1345%% data point positions are: (%f, %f)\n",
+                binLoF32, binHiF32);
+        #else
         // This code basically interpolates to find the positions exactly.
-        //
-        if (1) {
-            psTrace(__func__, 6, "binLo is %d.  Nums at that bin and the next are (%.2f, %.2f)\n",
-                    binLo, cumulativeRobustHistogram->nums->data.F32[binLo],
-                    cumulativeRobustHistogram->nums->data.F32[binLo+1]);
-            psTrace(__func__, 6, "binHi is %d.  Nums at that bin and the next are (%.2f, %.2f)\n",
-                    binHi, cumulativeRobustHistogram->nums->data.F32[binHi],
-                    cumulativeRobustHistogram->nums->data.F32[binHi+1]);
-
-            psF32 deltaNums;
-            psF32 deltaBounds;
-            psF32 prevPixels;
-            psF32 percentNums;
-            psF32 base;
-            deltaBounds = cumulativeRobustHistogram->bounds->data.F32[binLo+1] - cumulativeRobustHistogram->bounds->data.F32[binLo];
-            if (binLo == 0) {
-                deltaNums = cumulativeRobustHistogram->nums->data.F32[0];
-                prevPixels = 0;
-            } else {
-                deltaNums = cumulativeRobustHistogram->nums->data.F32[binLo] - cumulativeRobustHistogram->nums->data.F32[binLo-1];
-                prevPixels = cumulativeRobustHistogram->nums->data.F32[binLo-1];
-            }
-            percentNums = (totalDataPoints * 0.158655f) - prevPixels;
-            base = cumulativeRobustHistogram->bounds->data.F32[binLo];
-            binLoF32 = base + (deltaBounds / deltaNums) * percentNums;
-            psTrace(__func__, 6, "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
-                    base, deltaBounds, deltaNums, prevPixels, percentNums);
-
-            deltaBounds = cumulativeRobustHistogram->bounds->data.F32[binHi+1] - cumulativeRobustHistogram->bounds->data.F32[binHi];
-            if (binHi == 0) {
-                deltaNums = cumulativeRobustHistogram->nums->data.F32[0];
-                prevPixels = 0;
-            } else {
-                deltaNums = cumulativeRobustHistogram->nums->data.F32[binHi] - cumulativeRobustHistogram->nums->data.F32[binHi-1];
-                prevPixels = cumulativeRobustHistogram->nums->data.F32[binHi-1];
-            }
-            percentNums = (totalDataPoints * 0.841345f) - prevPixels;
-            base = cumulativeRobustHistogram->bounds->data.F32[binHi];
-            binHiF32 = base + (deltaBounds / deltaNums) * percentNums;
-            psTrace(__func__, 6, "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
-                    base, deltaBounds, deltaNums, prevPixels, percentNums);
-            psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", binLoF32, binHiF32);
-        }
-
-        //
-        // ADD: Step 5.
-        // Determine SIGMA as 1/2 of the distance between these positions.
-        //
+        psTrace(__func__, 6, "binLo is %d.  Nums at that bin and the next are (%.2f, %.2f)\n",
+                binLo, cumulative->nums->data.F32[binLo], cumulative->nums->data.F32[binLo+1]);
+        psTrace(__func__, 6, "binHi is %d.  Nums at that bin and the next are (%.2f, %.2f)\n",
+                binHi, cumulative->nums->data.F32[binHi], cumulative->nums->data.F32[binHi+1]);
+
+        float deltaBounds = cumulative->bounds->data.F32[binLo+1] - cumulative->bounds->data.F32[binLo];
+        float deltaNums;
+        float prevPixels;
+        if (binLo == 0) {
+            deltaNums = cumulative->nums->data.F32[0];
+            prevPixels = 0;
+        } else {
+            deltaNums = cumulative->nums->data.F32[binLo] - cumulative->nums->data.F32[binLo-1];
+            prevPixels = cumulative->nums->data.F32[binLo-1];
+        }
+        float percentNums = (totalDataPoints * 0.158655f) - prevPixels;
+        float base = cumulative->bounds->data.F32[binLo];
+        float binLoF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 15.8655% mark
+        psTrace(__func__, 6,
+                "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
+                base, deltaBounds, deltaNums, prevPixels, percentNums);
+
+        deltaBounds = cumulative->bounds->data.F32[binHi+1] - cumulative->bounds->data.F32[binHi];
+        if (binHi == 0) {
+            deltaNums = cumulative->nums->data.F32[0];
+            prevPixels = 0;
+        } else {
+            deltaNums = cumulative->nums->data.F32[binHi] - cumulative->nums->data.F32[binHi-1];
+            prevPixels = cumulative->nums->data.F32[binHi-1];
+        }
+        percentNums = (totalDataPoints * 0.841345f) - prevPixels;
+        base = cumulative->bounds->data.F32[binHi];
+        float binHiF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 84.1345% mark
+        psTrace(__func__, 6,
+                "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
+                base, deltaBounds, deltaNums, prevPixels, percentNums);
+        psTrace(__func__, 6,
+                "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n",
+                binLoF32, binHiF32);
+        #endif
+
+        // ADD step 5: Determine SIGMA as 1/2 of the distance between these positions.
         sigma = (binHiF32 - binLoF32) / 2.0;
         psTrace(__func__, 6, "The current sigma is %f.\n", sigma);
         stats->robustStdev = sigma;
 
-        //
-        // 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.
-        //
+        // 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 < (2 * binSize)) {
             psTrace(__func__, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
-            psS32 maskLo = PS_MAX(0, (binMedian - 25));
-            psS32 maskHi = PS_MIN(robustHistogram->bounds->n - 1, (binMedian + 25));
-            psF32 medianLo = robustHistogram->bounds->data.F32[maskLo];
-            psF32 medianHi = robustHistogram->bounds->data.F32[maskHi];
-            psTrace(__func__, 6, "Masking data more than 25 bins from the median (%.2f)\n", stats->robustMedian);
-            psTrace(__func__, 6, "The median is at bin number %d.  We mask bins outside the bin range (%d:%d)\n",
+            long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region
+            long maskHi = PS_MIN(histogram->bounds->n - 1, (binMedian + 25)); // High index for masking
+            psF32 medianLo = histogram->bounds->data.F32[maskLo]; // Value at low index
+            psF32 medianHi = histogram->bounds->data.F32[maskHi]; // Value at high index
+            psTrace(__func__, 6, "Masking data more than 25 bins from the median\n");
+            psTrace(__func__, 6,
+                    "The median is at bin number %d.  We mask bins outside the bin range (%d:%d)\n",
                     binMedian, maskLo, maskHi);
             psTrace(__func__, 6, "Masking data outside (%f %f)\n", medianLo, medianHi);
-            for (psS32 i = 0 ; i < myVector->n ; i++) {
+            for (long i = 0 ; i < myVector->n ; i++) {
                 if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) {
-                    tmpMaskVec->data.U8[i] = 1;
+                    mask->data.U8[i] = 0xff;
                     psTrace(__func__, 6, "Masking element %d is %f\n", i, myVector->data.F32[i]);
                 }
             }
-            psFree(robustHistogram);
-            psFree(cumulativeRobustHistogram);
-            iterate++;
-        } else {
+            // Free the histograms; they will be recreated on the next iteration, with new bounds
+            psFree(histogram);
+            psFree(cumulative);
+        } else {
+            // We've got the bin size correct now
             psTrace(__func__, 6, "*************: No more iteration.  sigma is %f\n", sigma);
-            iterate = 0;
-        }
-    }
-
-    //
-    // ADD: Step 7.
-    // Find the bins which contains the 25% and 75% data points.
-    //
-    psS32 binLo25 = 0;
-    psS32 binHi25 = 0;
+            iterate = -1;
+        }
+    }
+    psFree(histogram);
+
+    // ADD step 7: Find the bins which contains the 25% and 75% data points.
+    long binLo25 = 0;                   // Bin containing the 25% value
+    long binHi25 = 0;                   // Bin containing the 75% value
 
     tmpScalar.data.F32 = totalDataPoints * 0.25f;
-    if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
+    if (tmpScalar.data.F32 <= cumulative->nums->data.F32[0]) {
         // Special case where its in first bin.  The last bin should be handled automatically.
         binLo25 = 0;
     } else {
-        binLo25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
+        binLo25 = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
     }
     tmpScalar.data.F32 = totalDataPoints * 0.75f;
-    if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
+    if (tmpScalar.data.F32 <= cumulative->nums->data.F32[0]) {
         // Special case where its in first bin.  The last bin should be handled automatically.
         binHi25 = 0;
     } else {
-        binHi25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
+        binHi25 = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
     }
     if ((binLo25 < 0) || (binHi25 < 0)) {
         psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 25 and 75 percent data points\n");
-        psFree(tmpStatsMinMax);
-        psFree(robustHistogram);
-        psFree(cumulativeRobustHistogram);
-        psFree(tmpMaskVec);
+        psFree(statsMinMax);
+        psFree(cumulative);
+        psFree(mask);
+        psTrace(__func__, 4, "---- %s(false) end  ----\n", __func__);
+        return false;
+    }
+    psTrace(__func__, 6, "The 25-percent and 75-precent data point bins are (%d, %d).\n", binLo25, binHi25);
+
+    // ADD step 8: Interpolate to find these two positions exactly: these are the upper and lower quartile
+    // positions.
+    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums,
+                       binLo25, totalDataPoints * 0.25f);
+    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums,
+                       binHi25, totalDataPoints * 0.75f);
+    psFree(cumulative);
+
+    if (isnan(binLo25F32) || isnan(binHi25F32)) {
+        psError(PS_ERR_UNKNOWN, false,
+                "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n");
+        psFree(statsMinMax);
+        psFree(mask);
         psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-        return(1);
-    }
-    psTrace(__func__, 6, "The 25-percent and 75-precent data point bins are (%d, %d).\n", binLo25, binHi25);
-
-    //
-    // ADD: Step 8.
-    // Interpolate to find these two positions exactly: these are the upper
-    // and lower quartile positions.
-    //
-    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(
-                           *(psVector* *)&cumulativeRobustHistogram->bounds,
-                           *(psVector* *)&cumulativeRobustHistogram->nums,
-                           binLo25,
-                           totalDataPoints * 0.25f);
-
-    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(
-                           *(psVector* *)&cumulativeRobustHistogram->bounds,
-                           *(psVector* *)&cumulativeRobustHistogram->nums,
-                           binHi25,
-                           totalDataPoints * 0.75f);
-
-    if (isnan(binLo25F32) || isnan(binHi25F32)) {
-        psError(PS_ERR_UNKNOWN, false, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n");
-        psFree(tmpStatsMinMax);
-        psFree(robustHistogram);
-        psFree(cumulativeRobustHistogram);
-        psFree(tmpMaskVec);
-        psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-        return(1);
+        return false;
     }
 
     stats->robustLQ = binLo25F32;
     stats->robustUQ = binHi25F32;
-    psTrace(__func__, 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n", binLo25F32, binHi25F32);
-    psS32 N50 = 0;
-    for (psS32 i = 0 ; i < myVector->n ; i++) {
-        // XXX: use maskVal here?
-        if ((0 == tmpMaskVec->data.U8[i]) &&
-                (binLo25F32 <= myVector->data.F32[i]) &&
-                (binHi25F32 >= myVector->data.F32[i])) {
+    psTrace(__func__, 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n",
+            binLo25F32, binHi25F32);
+    long N50 = 0;
+    for (long i = 0 ; i < myVector->n ; i++) {
+        if (!(mask->data.U8[i] & maskVal) &&
+                (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) {
             N50++;
         }
@@ -1740 +1440 @@
     psTrace(__func__, 6, "The robustN50 is %d.\n", N50);
 
-    // ************************************************************************
-    if ((stats->options & PS_STAT_FITTED_MEAN) ||
-            (stats->options & PS_STAT_FITTED_STDEV)) {
-        //
+
+    // Fitted statistics
+    if (stats->options & PS_STAT_FITTED_MEAN || stats->options & PS_STAT_FITTED_STDEV) {
         // New algorithm for calculated mean and stdev.
-        //
 
         // XXX: Where is this documented?
-        psF32 dN = (psF32) (0.17 * N50);
+        psF32 dN = 0.17 * N50;
         if (dN < 1.0) {
             dN = 1.0;
@@ -1756 +1454 @@
         psF32 newBinSize = sigma / dN;
 
-        //
         // Determine the min/max of the vector (which prior outliers masked out)
-        //
-        rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
-        rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
-        if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
+        int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values
+        min = statsMinMax->min;
+        max = statsMinMax->max;
+        if (numValid == 0 || isnan(min) || isnan(max)) {
             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
-            psFree(tmpStatsMinMax);
-            psFree(robustHistogram);
-            psFree(cumulativeRobustHistogram);
-            psFree(tmpMaskVec);
-            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-            return(1);
-        }
-
-        //
+            psFree(statsMinMax);
+            psFree(mask);
+            psTrace(__func__, 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
         // Calculate the number of bins.
-        //
-        numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize);
-        psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
+        long numBins = (max - min) / newBinSize;
+        psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", min, max);
         psTrace(__func__, 6, "The new bin size is %f.\n", newBinSize);
         psTrace(__func__, 6, "The numBins is %d\n", numBins);
 
-        psHistogram *newHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
-        newHistogram = psVectorHistogram(newHistogram, myVector, errors, tmpMaskVec, maskVal|1);
+        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
+        histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal);
         if (psTraceGetLevel(__func__) >= 8) {
-            PS_VECTOR_PRINT_F32(newHistogram->nums);
-        }
-
-        //
-        // FITTED STATISTICS HERE
-        //
-        // Smooth the resulting histogram with a Gaussian with sigma_s = 1
-        // bin.
-        //
-        psF32 sigma_s = newBinSize;
-
-        psVector *newHistogramSmoothed = p_psVectorSmoothHistGaussian(newHistogram, sigma_s);
+            PS_VECTOR_PRINT_F32(histogram->nums);
+        }
+
+        // Smooth the resulting histogram with a Gaussian with sigma_s = 1 bin.
+        psVector *smoothed = vectorSmoothHistGaussian(histogram, newBinSize); // Smoothed histogram
         if (psTraceGetLevel(__func__) >= 8) {
-            PS_VECTOR_PRINT_F32(newHistogramSmoothed);
-        }
-
-        //
-        // Find the bin with the peak value in the range 2 sigma of the
-        // robust histogram median.
-        //
-
-        psS32 binMin = 0;
-        psS32 binMax = 0;
+            PS_VECTOR_PRINT_F32(smoothed);
+        }
+
+        // Find the bin with the peak value in the range 2 sigma of the robust histogram median.
+        long binMin = 0;                // Low bin to check
+        long binMax = 0;                // High bin to check
         tmpScalar.data.F32 = stats->robustMedian - (2.0 * sigma);
-        if (tmpScalar.data.F32 <= newHistogram->bounds->data.F32[0]) {
+        if (tmpScalar.data.F32 <= histogram->bounds->data.F32[0]) {
             binMin = 0;
         } else {
-            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
+            binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
         }
 
         tmpScalar.data.F32 = stats->robustMedian + (2.0 + sigma);
-        if (tmpScalar.data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) {
-            binMax = newHistogram->bounds->n-1;
-        } else {
-            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
+        if (tmpScalar.data.F32 >= histogram->bounds->data.F32[histogram->bounds->n-1]) {
+            binMax = histogram->bounds->n-1;
+        } else {
+            binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
         }
         if ((binMin < 0) || (binMax < 0)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +- 2.0 * sigma bins\n");
-            psFree(tmpMaskVec);
-            psFree(robustHistogram);
-            psFree(cumulativeRobustHistogram);
-            psFree(tmpStatsMinMax);
-            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-            return(1);
-        }
-
-        psS32 binNum = binMin;
-        psF32 binMaxNums = newHistogramSmoothed->data.F32[binNum];
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +/- 2.0 * sigma bins\n");
+            psFree(statsMinMax);
+            psFree(histogram);
+            psFree(mask);
+            psTrace(__func__, 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        long binNum = binMin;           // Bin number containing the peak
+        psF32 binMaxNums = smoothed->data.F32[binNum]; // The peak value
         for (psS32 i = binMin+1 ; i <= binMax ; i++) {
-            if (newHistogramSmoothed->data.F32[i] > binMaxNums) {
+            if (smoothed->data.F32[i] > binMaxNums) {
                 binNum = i;
-                binMaxNums = newHistogramSmoothed->data.F32[i];
+                binMaxNums = smoothed->data.F32[i];
             }
         }
         psTrace(__func__, 6, "The peak bin is %d, with %f data.n", binNum, binMaxNums);
 
-        //
-        // Fit a Gaussian to the bins in the range 20 sigma of the robust
-        // histogram median.
-        //
+        // Fit a Gaussian to the bins in the range 20 sigma of the robust histogram median.
         tmpScalar.data.F32 = stats->robustMedian - (20.0 * sigma);
-        if (tmpScalar.data.F32 < tmpStatsMinMax->min) {
+        if (tmpScalar.data.F32 < min) {
             binMin = 0;
         } else {
-            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
-            // XXX: check for errors here.
+            binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
         }
         tmpScalar.data.F32 = stats->robustMedian + (20.0 * sigma);
-        if (tmpScalar.data.F32 > tmpStatsMinMax->max) {
-            binMax = newHistogramSmoothed->n - 1;
-        } else {
-            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
-            // XXX: check for errors here.
-        }
-        psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
-        psArray *x = psArrayAlloc((1 + (binMax - binMin)));
-        psS32 j = 0;
+        if (tmpScalar.data.F32 > max) {
+            binMax = smoothed->n - 1;
+        } else {
+            binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
+        }
+        if (binMin < 0 || binMax < 0) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +/- 20.0 * sigma bins\n");
+            psFree(statsMinMax);
+            psFree(histogram);
+            psFree(mask);
+            psTrace(__func__, 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        // Generate the variables that will be used in the Gaussian fitting
+        psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates
+        psArray *x = psArrayAlloc((1 + (binMax - binMin))); // Array of ordinates
         y->n = y->nalloc;
-
-        for (psS32 i = binMin ; i <= binMax ; i++) {
-            y->data.F32[j] = newHistogramSmoothed->data.F32[i];
-            x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32);
-            x->n++;
-            ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i);
-            ((psVector *) x->data[j])->n++;
-            j++;
+        x->n = x->nalloc;
+        for (long i = binMin, j = 0; i <= binMax ; i++, j++) {
+            y->data.F32[j] = smoothed->data.F32[i];
+            psVector *ordinate = psVectorAlloc(1, PS_TYPE_F32); // The ordinate value
+            ordinate->n = 1;
+            ordinate->data.F32[0] = PS_BIN_MIDPOINT(histogram, i);
+            x->data[j] = ordinate;
         }
         if (psTraceGetLevel(__func__) >= 8) {
@@ -1873 +1557 @@
             PS_VECTOR_PRINT_F32(y);
         }
-
-        rc = p_psVectorMin(y, NULL, 0, tmpStatsMinMax);
-        rc|= p_psVectorMax(y, NULL, 0, tmpStatsMinMax);
-        if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
+        psFree(smoothed);
+        psFree(histogram);
+
+        numValid = vectorMinMax(y, NULL, 0, statsMinMax); // Number of valid values
+        min = statsMinMax->min;
+        max = statsMinMax->max;
+        if (numValid == 0 || isnan(min) || isnan(max)) {
             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
-            psFree(tmpMaskVec);
-            psFree(robustHistogram);
-            psFree(cumulativeRobustHistogram);
-            psFree(tmpStatsMinMax);
-            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-            return(1);
-        }
-
-        //
+            psFree(mask);
+            psFree(histogram);
+            psFree(statsMinMax);
+            psTrace(__func__, 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
         // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
-        // XXX: Use the normalize routines for this.
-        //
-        for (psS32 i = 0 ; i < y->n ; i++) {
-            y->data.F32[i]= (y->data.F32[i] - tmpStatsMinMax->min) / (tmpStatsMinMax->max - tmpStatsMinMax->min);
-        }
-
-        //
-        psMinimization *min = psMinimizationAlloc(100, 0.01);
-        psVector *params = psVectorAlloc(2, PS_TYPE_F32);
+        p_psNormalizeVectorRange(y, 0.0, 1.0);
+
+        // Fit a Gaussian to the data.
+        psMinimization *minimizer = psMinimizationAlloc(100, 0.01); // The minimizer information
+        psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian
         params->n = params->nalloc;
-        // Initial guess for the mean ([0]) and standard dev.
+        // Initial guess for the mean (index 0) and standard dev (index 1).
         params->data.F32[0] = stats->robustMedian;
         params->data.F32[1] = sigma;
@@ -1905 +1586 @@
             PS_VECTOR_PRINT_F32(y);
         }
-
-        //
-        // Fit a Gaussian to the data.
-        //
-        rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D);
-        if (rcBool != true) {
+        if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) {
             psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
-            psFree(tmpStatsMinMax);
-            psFree(robustHistogram);
-            psFree(cumulativeRobustHistogram);
-            psFree(newHistogram);
+            psFree(statsMinMax);
             psFree(x);
             psFree(y);
-            psFree(min);
+            psFree(minimizer);
             psFree(params);
-            psFree(tmpMaskVec);
-            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-            return(1);
+            psFree(mask);
+            psTrace(__func__, 4, "---- %s(false) end  ----\n", __func__);
+            return false;
         }
         if (psTraceGetLevel(__func__) >= 8) {
@@ -1928 +1601 @@
         }
 
-        //
-        // The fitted mean mean_r is derived directly from the fitted
-        // Gaussian mean.
-        //
+        // The fitted mean is the Gaussian mean.
         stats->fittedMean = params->data.F32[0];
         psTrace(__func__, 6, "The fitted mean is %f.\n", params->data.F32[0]);
 
-        //
-        // The fitted standard deviation, SIGMA_r is determined by
-        // subtracting the smoothing scale in quadrature:
-        //     SIGMA_r^2 = SIGMA^2 - sigma_s^2
-        //
-        stats->fittedStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(sigma_s));
+        // The fitted standard deviation, SIGMA_r is determined by subtracting the smoothing scale in
+        // quadrature: SIGMA_r^2 = SIGMA^2 - sigma_s^2
+        stats->fittedStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(newBinSize));
         psTrace(__func__, 6, "The fitted stdev is %f.\n", stats->fittedStdev);
 
-        psFree(newHistogramSmoothed);
-        psFree(newHistogram);
+        // Clean up after fitting
         psFree(x);
         psFree(y);
-        psFree(min);
+        psFree(minimizer);
         psFree(params);
     }
 
-    psFree(tmpStatsMinMax);
-    psFree(cumulativeRobustHistogram);
-    psFree(tmpMaskVec);
-    psFree(robustHistogram);
+    // Clean up
+    psFree(statsMinMax);
+    psFree(mask);
 
     psTrace(__func__, 4, "---- %s(0) end  ----\n", __func__);
-    return(0);
+    return true;
 }
 
@@ -1968 +1633 @@
 static void histogramFree(psHistogram* myHist);
 
+// We keep statsFree so that we can identify statistics pointers from the memblock
 static void statsFree(psStats *newStruct)
 {
@@ -2038 +1704 @@
     PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(upper, lower, NULL);
 
-    psS32 i = 0;                  // Loop index variable
-    psHistogram* newHist = NULL;        // The new histogram structure
-    psF32 binSize = 0.0;        // The histogram bin size
-
-    // Allocate memory for the new histogram structure.  If there are N
-    // bins, then there are N+1 bounds to those bins.
-    newHist = (psHistogram* ) psAlloc(sizeof(psHistogram));
+    // Allocate memory for the new histogram structure.  If there are N bins, then there are N+1 bounds to
+    // those bins.
+    psHistogram *newHist = (psHistogram* ) psAlloc(sizeof(psHistogram)); // The new histogram structure
     psMemSetDeallocator(newHist, (psFreeFunc) histogramFree);
     psVector* newBounds = psVectorAlloc(n + 1, PS_TYPE_F32);
@@ -2051 +1713 @@
 
     // Calculate the bounds for each bin.
-    binSize = (upper - lower) / (psF32)n;
-    // XXX: Is the following necessary? It prevents the max data point
-    // from being in a non-existant bin.
+    psF32 binSize = (upper - lower) / (psF32)n; // The histogram bin size
+    // XXX: Is the following necessary? It prevents the max data point from being in a non-existant bin.
     binSize += FLT_EPSILON;
-    for (i = 0; i < n + 1; i++) {
+    for (long i = 0; i < n + 1; i++) {
         newBounds->data.F32[i] = lower + (binSize * (psF32)i);
     }
@@ -2061 +1722 @@
     // Allocate the bins, and initialize them to zero.
     newHist->nums = psVectorAlloc(n, PS_TYPE_F32);
-    for (i = 0; i < newHist->nums->nalloc; i++) {
+    for (long i = 0; i < newHist->nums->nalloc; i++) {
         newHist->nums->data.F32[i] = 0.0;
         newHist->nums->n++;
@@ -2072 +1733 @@
 
     psTrace(__func__, 3, "---- %s() end  ----\n", __func__);
-    return (newHist);
+    return newHist;
 }
 
@@ -2091 +1752 @@
     PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(bounds->n, 2, NULL);
 
-    psHistogram* newHist = NULL;        // The new histogram structure
-    psS32 i;                      // Loop index variable
-
     // Allocate memory for the new histogram structure.
-    newHist = (psHistogram* ) psAlloc(sizeof(psHistogram));
+    psHistogram *newHist = (psHistogram* ) psAlloc(sizeof(psHistogram)); // The new histogram structure
     psMemSetDeallocator(newHist, (psFreeFunc) histogramFree);
     psVector* newBounds = psVectorAlloc(bounds->n, PS_TYPE_F32);
     newHist->bounds = newBounds;
     newBounds->n = newHist->bounds->nalloc;
-    for (i = 0; i < bounds->n; i++) {
+    for (long i = 0; i < bounds->n; i++) {
         newBounds->data.F32[i] = bounds->data.F32[i];
     }
@@ -2107 +1765 @@
     // then there are N-1 bins.
     newHist->nums = psVectorAlloc((bounds->n) - 1, PS_TYPE_F32);
-    for (i = 0; i < newHist->nums->nalloc; i++) {
+    for (long i = 0; i < newHist->nums->nalloc; i++) {
         newHist->nums->data.F32[i] = 0.0;
         newHist->nums->n++;
@@ -2140 +1798 @@
 the data point as a boxcar PDF and update a range of points surrounding the
 histogram bin which contains the point.  The width of that boxcar is defined
-as 2.35 * error.  Inputs:
-    binNum: the bin number of the data point in the histogram
-    out: the histogram structure
-    data: the data point value
-    error: the error in that data point
+as 2.35 * error.
  
 XXX: Must test this.
  *****************************************************************************/
-psS32 UpdateHistogramBins(psS32 binNum,
-                          psHistogram* out,
-                          psF32 data,
-                          psF32 error)
+static bool UpdateHistogramBins(long binNum, // Bin number of the data point
+                                psHistogram* out, // The histogram to be updated
+                                psF32 data, // The data point value
+                                psF32 error // The error in the data point
+                               )
 {
     psTrace(__func__, 3, "---- %s() begin  ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(out, -1);
-    PS_ASSERT_PTR_NON_NULL(out->bounds, -1);
-    PS_ASSERT_PTR_NON_NULL(out->nums, -1);
-    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, ((out->nums->n)-1), -2);
-    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(error, 0.0, -3);
-
-    PS_ASSERT_FLOAT_WITHIN_RANGE(data, out->bounds->data.F32[0], out->bounds->data.F32[(out->bounds->n)-1], -4);
-
-    psF32 boxcarWidth = 2.35 * error;
+    PS_ASSERT_PTR_NON_NULL(out, false);
+    PS_ASSERT_PTR_NON_NULL(out->bounds, false);
+    PS_ASSERT_PTR_NON_NULL(out->nums, false);
+    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, ((out->nums->n)-1), false);
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(error, 0.0, false);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(data, out->bounds->data.F32[0],
+                                 out->bounds->data.F32[(out->bounds->n)-1], false);
+
+    psF32 boxcarWidth = 2.35 * error;   // Width of the boxcar
     psF32 boxcarCenter = (out->bounds->data.F32[binNum] +
-                          out->bounds->data.F32[binNum+1]) / 2.0;
-    psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0);
-    psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0);
-    psS32 bin;
-    psS32 boxcarLeftBinNum = 0;
-    psS32 boxcarRightBinNum = 0;
+                          out->bounds->data.F32[binNum+1]) / 2.0; // Centre of the boxcar
+    psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0); // Left endpoint of the boxcar for the PDF
+    psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0); // Right endpoint of the boxcar for the PDF
+    psS32 boxcarLeftBinNum = 0;         // Bin number for left endpoint
+    psS32 boxcarRightBinNum = 0;        // Bin number for right endpoint
 
     // Determine the left endpoint of the boxcar for the PDF.
-    for (bin=binNum ; bin >= 0 ; bin--) {
+    for (long bin = binNum; bin >= 0; bin--) {
         if (out->nums->data.F32[bin] <= boxcarLeft) {
             boxcarLeftBinNum = bin;
@@ -2180 +1834 @@
 
     // Determine the right endpoint of the boxcar for the PDF.
-    for (bin=binNum ; bin < out->nums->n ; bin++) {
+    for (long bin = binNum; bin < out->nums->n; bin++) {
         if (out->nums->data.F32[bin] >= boxcarRight) {
             boxcarRightBinNum = bin;
@@ -2187 +1841 @@
     }
 
-    //
     // If the boxcar fits entirely inside this bin, then simply add 1.0 to the
     // bin and return.
-    //
     if (boxcarLeftBinNum == boxcarRightBinNum) {
-        out->nums->data.F32[binNum]+= 1.0;
-        psTrace(__func__, 3, "---- %s(0) end  ----\n", __func__);
-        return(0);
-    }
-
-    //
-    // If we get here, multiple bins must be updated.  We handle the left
-    // endpoint, and right endpoint differently.
-    //
-    out->nums->data.F32[boxcarLeftBinNum]+=
+        out->nums->data.F32[binNum] += 1.0;
+        psTrace(__func__, 3, "---- %s(true) end  ----\n", __func__);
+        return true;
+    }
+
+    // If we get here, multiple bins must be updated.  We handle the left-most endpoint, and right-most
+    // endpoints differently.
+    out->nums->data.F32[boxcarLeftBinNum] +=
         (out->bounds->data.F32[boxcarLeftBinNum+1] - boxcarLeft) / boxcarWidth;
 
-    //
     // Loop through the center bins, if any.
-    //
-    for (bin = boxcarLeftBinNum + 1 ; bin < (boxcarRightBinNum - 1) ; bin++) {
-        out->nums->data.F32[bin]+=
+    for (long bin = boxcarLeftBinNum + 1; bin < (boxcarRightBinNum - 1); bin++) {
+        out->nums->data.F32[bin] +=
             (out->bounds->data.F32[bin+1] - out->bounds->data.F32[bin]) / boxcarWidth;
     }
 
-    //
     // Handle the right endpoint differently.
-    //
     out->nums->data.F32[boxcarRightBinNum]+=
         (boxcarRight - out->bounds->data.F32[boxcarRightBinNum]) / boxcarWidth;
 
-    //
-    // Return 0 on success.
-    //
-    psTrace(__func__, 3, "---- %s(0) end  ----\n", __func__);
-    return(0);
+    psTrace(__func__, 3, "---- %s(true) end  ----\n", __func__);
+    return true;
 }
 
@@ -2255 +1898 @@
     PS_ASSERT_INT_NONNEGATIVE(out->nums->n, NULL);
     PS_ASSERT_VECTOR_NON_NULL(values, out);
-    if (mask != NULL) {
+    if (mask) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(values, mask, NULL);
         PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
     }
-    if (errors != NULL) {
+    if (errors) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(values, errors, NULL);
         PS_ASSERT_VECTOR_TYPE(errors, values->type.type, NULL);
     }
 
-    psS32 i = 0;                  // Loop index variable
-    psF32 binSize = 0.0;          // Histogram bin size
-    psS32 binNum = 0;             // A temporary bin number
-    psS32 numBins = 0;            // The total number of bins
+    long binNum = 0;                    // A temporary bin number
+    long numBins = out->nums->n;        // The total number of bins
     psScalar tmpScalar;
     tmpScalar.type.type = PS_TYPE_F32;
-    psVector* inF32 = NULL;
-    psVector* errorsF32 = NULL;
-    psS32 mustFreeVectorIn = 1;
-    psS32 mustFreeVectorErrors = 1;
 
     // Convert input and errors vectors to F32 if necessary.
-    inF32 = p_psConvertToF32((psVector *) values);
-    if (inF32 == NULL) {
-        inF32 = (psVector *) values;
-        mustFreeVectorIn = 0;
-    }
-    errorsF32 = p_psConvertToF32((psVector *) errors);
-    if (errorsF32 == NULL) {
-        errorsF32 = (psVector *) errors;
-        mustFreeVectorErrors = 0;
-    }
-
-    numBins = out->nums->n;
-    for (i = 0; i < inF32->n; i++) {
+    psVector* inF32 = NULL;             // F32 version of input vector
+    if (values->type.type == PS_TYPE_F32) {
+        inF32 = psMemIncrRefCounter((psPtr)values);
+    } else {
+        values = psVectorCopy(NULL, values, PS_TYPE_F32);
+    }
+    psVector* errorsF32 = NULL;         // F32 version of errors vector
+    if (errors) {
+        if (errors->type.type == PS_TYPE_F32) {
+            errorsF32 = psMemIncrRefCounter((psPtr)errors);
+        } else {
+            errorsF32 = psVectorCopy(NULL, errors, PS_TYPE_F32);
+        }
+    }
+
+    for (long i = 0; i < inF32->n; i++) {
         // Check if this pixel is masked, and if so, skip it.
-        if ((mask == NULL) || ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) {
+        if (!mask || (mask && (!(mask->data.U8[i] & maskVal)))) {
             if (inF32->data.F32[i] < out->bounds->data.F32[0]) {
                 // If this pixel is below minimum value, count it, then skip.
@@ -2301 +1941 @@
                 // number is trivial.
                 if (out->uniform == true) {
-                    binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
-                    binNum = (psS32)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize);
-                    if (errorsF32 != NULL) {
-                        psS32 rc = UpdateHistogramBins(binNum, out, inF32->data.F32[i],
-                                                       errorsF32->data.F32[i]);
-                        if (rc < 0) {
-                            psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram bins with the errors vector.\n");
+                    float binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; // Histogram bin size
+                    binNum = (inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize;
+                    if (errorsF32) {
+                        if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errorsF32->data.F32[i])) {
+                            psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram "
+                                     "bins with the errors vector.\n");
                         }
                     } else {
@@ -2329 +1968 @@
                     } else {
                         if (errorsF32 != NULL) {
-                            psS32 rc = UpdateHistogramBins(binNum, out,
-                                                           inF32->data.F32[i],
-                                                           errors->data.F32[i]);
-                            if (rc < 0) {
-                                psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram bins with the errors vector.\n");
+                            if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errors->data.F32[i])) {
+                                psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram "
+                                         "bins with the errors vector.\n");
                             }
                         } else {
-                            (out->nums->data.F32[binNum])+= 1.0;
+                            out->nums->data.F32[binNum] += 1.0;
                         }
                     }
@@ -2344 +1981 @@
     }
 
-    if (mustFreeVectorIn == 1) {
-        psFree(inF32);
-    }
-    if (mustFreeVectorErrors == 1) {
-        psFree(errorsF32);
-    }
+    psFree(inF32);
+    psFree(errorsF32);
 
     psTrace(__func__, 3, "---- %s() end  ----\n", __func__);
     return (out);
-}
-
-/******************************************************************************
-p_psConvertToF32(in): this is the cheap way to support a variety of vector
-data types: we simply convert the input vector to F32 at the beginning, and
-write all of our functions in F32.  If the vast majority of all vector stat
-operations are F32 (or any other single type), then this is probably the
-best way to go.  Otherwise, when the algorithms stablize, we will then macro
-everything and put type support in the various stat functions.
- 
-XXX: Should the default data type be F64?  Since we are buying Opterons...
- 
-XXX: Use psLib functions instead.
- *****************************************************************************/
-psVector* p_psConvertToF32(psVector* in)
-{
-    psTrace(__func__, 4,"---- %s() begin  ----\n", __func__);
-    if (in == NULL) {
-        return(NULL);
-    }
-    psS32 i = 0;
-    psVector* tmp = NULL;
-
-    if (in->type.type == PS_TYPE_S8) {
-        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
-        for (i = 0; i < in->n; i++) {
-            tmp->data.F32[i] = (psF32)in->data.S8[i];
-            tmp->n++;
-        }
-    } else if (in->type.type == PS_TYPE_S16) {
-        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
-        for (i = 0; i < in->n; i++) {
-            tmp->data.F32[i] = (psF32) in->data.S16[i];
-            tmp->n++;
-        }
-    } else if (in->type.type == PS_TYPE_S32) {
-        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
-        for (i = 0; i < in->n; i++) {
-            tmp->data.F32[i] = (psF32)in->data.S32[i];
-            tmp->n++;
-        }
-    } else if (in->type.type == PS_TYPE_S64) {
-        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
-        for (i = 0; i < in->n; i++) {
-            tmp->data.F32[i] = (psF32)in->data.S64[i];
-            tmp->n++;
-        }
-    } else if (in->type.type == PS_TYPE_U8) {
-        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
-        for (i = 0; i < in->n; i++) {
-            tmp->data.F32[i] = (psF32)in->data.U8[i];
-            tmp->n++;
-        }
-    } else if (in->type.type == PS_TYPE_U16) {
-        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
-        for (i = 0; i < in->n; i++) {
-            tmp->data.F32[i] = (psF32)in->data.U16[i];
-            tmp->n++;
-        }
-    } else if (in->type.type == PS_TYPE_U32) {
-        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
-        for (i = 0; i < in->n; i++) {
-            tmp->data.F32[i] = (psF32)in->data.U32[i];
-            tmp->n++;
-        }
-    } else if (in->type.type == PS_TYPE_U64) {
-        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
-        for (i = 0; i < in->n; i++) {
-            tmp->data.F32[i] = (psF32)in->data.U64[i];
-            tmp->n++;
-        }
-    } else if (in->type.type == PS_TYPE_F64) {
-        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
-        for (i = 0; i < in->n; i++) {
-            tmp->data.F32[i] = (psF32)in->data.F64[i];
-            tmp->n++;
-        }
-    } else if (in->type.type == PS_TYPE_F32) {
-        // do nothing
-    } else {
-        char* strType;
-        PS_TYPE_NAME(strType, in->type.type);
-        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
-                PS_ERRORTEXT_psStats_VECTOR_TYPE_UNSUPPORTED,
-                strType);
-    }
-
-    psTrace(__func__, 4,"---- %s() end  ----\n", __func__);
-    return (tmp);
 }
 
@@ -2469 +2013 @@
     PS_ASSERT_PTR_NON_NULL(stats, NULL);
     PS_ASSERT_VECTOR_NON_NULL(in, NULL);
-    if (mask != NULL) {
+    if (mask) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(mask, in, stats);
         PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, stats);
     }
-    if (errors != NULL) {
+    if (errors) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(errors, in, stats);
         PS_ASSERT_VECTOR_TYPE(errors, in->type.type, stats);
     }
-    // XXX: Assert that "in" is F64, F32, U16, or S8
-
-    psVector* inF32 = NULL;
-    psVector* errorsF32 = NULL;
-    psS32 mustFreeVectorIn = 1;
-    psS32 mustFreeVectorErrors = 1;
-
-    inF32 = p_psConvertToF32((psVector *) in);
-    if (inF32 == NULL) {
-        //        printf("\n\ninF32 is NULL\n\n");
-        //        inF32 = (psVector *) in;
-        //        mustFreeVectorIn = 0;
-        inF32 = psVectorCopy(inF32, in, PS_TYPE_F32);
-        inF32->n = inF32->nalloc;
-        mustFreeVectorIn = 1;
-    }
-    //    else printf("\ninF32 has n=%ld, nalloc=%ld\n", inF32->n, inF32->nalloc);
-    errorsF32 = p_psConvertToF32((psVector *) errors);
-    if (errorsF32 == NULL) {
-        errorsF32 = (psVector *) errors;
-        mustFreeVectorErrors = 0;
-    }
-    //    inF32->n = inF32->nalloc;
-    //    errorsF32->n = errorsF32->nalloc;
+
+    // Convert types, as necessary
+    psVector *inF32 = NULL;             // Input vector of values, F32 version
+    if (in->type.type == PS_TYPE_F32) {
+        inF32 = psMemIncrRefCounter((psPtr)in);
+    } 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) {
+            errorsF32 = psMemIncrRefCounter((psPtr)errors);
+        } else {
+            errorsF32 = psVectorCopy(NULL, errors, PS_TYPE_F32);
+        }
+    }
+    psVector *maskU8 = NULL;            // Input mask vector, U8 version
+    if (mask) {
+        // We want a copy of the mask, since we may change values
+        maskU8 = psVectorCopy(NULL, mask, PS_TYPE_U8);
+    }
+
+
     if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) {
         PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->max, stats->min, stats);
@@ -2511 +2054 @@
     // ************************************************************************
     if (stats->options & PS_STAT_SAMPLE_MEAN) {
-        if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMean() returned an error.\n");
+        if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) {
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMean() returned an error.\n");
             stats->sampleMean = NAN;
         }
@@ -2518 +2061 @@
 
     // ************************************************************************
-    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
-        if (false == p_psVectorSampleMedian(inF32, mask, maskVal, stats)) {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian() returned an error.\n");
+    if (stats->options & (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE)) {
+        if (!vectorSampleMedian(inF32, maskU8, maskVal, stats)) {
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMedian() returned an error.\n");
             stats->sampleMedian = NAN;
+            stats->sampleUQ = NAN;
+            stats->sampleLQ = NAN;
         }
     }
@@ -2527 +2072 @@
     // ************************************************************************
     if (stats->options & PS_STAT_SAMPLE_STDEV) {
-        if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMean() returned an error.\n");
+        if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) {
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMean() returned an error.\n");
             stats->sampleMean = NAN;
         } else {
-            p_psVectorSampleStdev(inF32, errorsF32, mask, maskVal, stats);
+            vectorSampleStdev(inF32, errorsF32, maskU8, maskVal, stats);
         }
     }
 
     // ************************************************************************
-    if (stats->options & PS_STAT_SAMPLE_QUARTILE) {
-        if (false == p_psVectorSampleQuartiles(inF32, mask, maskVal, stats)) {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleQuartiles() returned an error.\n");
-            stats->sampleLQ = NAN;
-            stats->sampleUQ = NAN;
-        }
-    }
-
-    // ************************************************************************
-    if ((stats->options & PS_STAT_ROBUST_MEDIAN) ||
-            (stats->options & PS_STAT_ROBUST_STDEV) ||
-            (stats->options & PS_STAT_ROBUST_QUARTILE) ||
-            (stats->options & PS_STAT_FITTED_MEAN) ||
-            (stats->options & PS_STAT_FITTED_STDEV)) {
-        if (0 != p_psVectorRobustStats(inF32, errorsF32, mask, maskVal, stats)) {
+    if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE |
+                          PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
+        if (!vectorRobustStats(inF32, errorsF32, maskU8, maskVal, stats)) {
             psError(PS_ERR_UNKNOWN, false, PS_ERRORTEXT_psStats_STATS_FAILED);
             psFree(stats);
+            psFree(inF32);
+            psFree(errorsF32);
             psTrace(__func__, 3,"---- %s(NULL) end  ----\n", __func__);
             return(NULL);
@@ -2560 +2095 @@
     // ************************************************************************
     if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
-        psS32 rc = p_psVectorClippedStats(inF32, errorsF32, mask, maskVal, stats);
-        if (rc < 0) {
+        if (!vectorClippedStats(inF32, errorsF32, maskU8, maskVal, stats)) {
             psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics for input psVector.\n");
             stats->clippedMean = NAN;
@@ -2569 +2103 @@
 
     // ************************************************************************
-    if (stats->options & PS_STAT_MAX) {
-        if (0 != p_psVectorMax(inF32, mask, maskVal, stats)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector maximum");
-            stats->max = NAN;
-        }
-    }
-
-    // ************************************************************************
-    if (stats->options & PS_STAT_MIN) {
-        if (0 != p_psVectorMin(inF32, mask, maskVal, stats)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum");
-            stats->min = NAN;
-        }
-    }
-
-    if (mustFreeVectorIn == 1) {
-        psFree(inF32);
-    }
-    if (mustFreeVectorErrors == 1) {
-        psFree(errorsF32);
-    }
+    if (stats->options & (PS_STAT_MAX | PS_STAT_MIN)) {
+        if (vectorMinMax(inF32, maskU8, maskVal, stats) == 0) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum and maximum");
+        }
+    }
+
+    psFree(inF32);
+    psFree(errorsF32);
+    psFree(maskU8);
     psTrace(__func__, 3,"---- %s() end  ----\n", __func__);
     return (stats);