Changeset 10848 for trunk/psLib/src/math/psMinimizePolyFit.c

Timestamp:

Dec 28, 2006, 6:38:42 PM (19 years ago)

Author:

magnier

Message:

changed return value for psVectorFit functions to bool; now using stats->option to set FitClip

File:

• trunk/psLib/src/math/psMinimizePolyFit.c (modified) (72 diffs)

trunk/psLib/src/math/psMinimizePolyFit.c

@@ -10 +10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-12-17 09:43:48 $
+ *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-12-29 04:38:42 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -55 +55 @@
     }\
 }\
+
+// free a local temporary F64 vector (TEMP) which is a copy of a non-F64 vector (ORIG)
+# define PS_FREE_TEMP_F64_VECTOR(ORIG, TEMP) \
+if ((ORIG != NULL) && (ORIG->type.type != PS_TYPE_F64)) { psFree(TEMP); }
 
 /*****************************************************************************/
@@ -281 +285 @@
  *****************************************************************************/
 
-
 /******************************************************************************
  ******************************************************************************
@@ -288 +291 @@
  *****************************************************************************/
 
-static psPolynomial1D* VectorFitPolynomial1DOrd(
-    psPolynomial1D* myPoly,
-    const psVector *mask,
-    psMaskType maskValue,
-    const psVector *f,
-    const psVector *fErr,
-    const psVector *x);
-
 /******************************************************************************
 vectorFitPolynomial1DCheb():  This routine will fit a Chebyshev
@@ -301 +296 @@
 coefficients of that polynomial.
 *****************************************************************************/
-static psPolynomial1D *vectorFitPolynomial1DCheb(
+static bool vectorFitPolynomial1DCheb(
     psPolynomial1D* myPoly,
     const psVector *mask,
@@ -404 +399 @@
         if (!psMatrixGJSolve(A, B)) {
             psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
-            psFree(myPoly);
-            myPoly = NULL;
+            goto escape_GJ;
         } else {
             // the first nTerm entries in B correspond directly to the desired
@@ -427 +421 @@
         if (ALUD == NULL) {
             psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
-            psFree(myPoly);
-            myPoly = NULL;
+            goto escape_LUD;
         } else {
             coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
             if (coeffs == NULL) {
                 psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
-                psFree(myPoly);
-                myPoly = NULL;
+                goto escape_LUD;
             } else {
                 for (psS32 k = 0; k < numTerms; k++) {
@@ -441 +433 @@
             }
         }
-
-
         psFree(ALUD);
         psFree(coeffs);
         psFree(outPerm);
     }
-
     psFree(A);
     psFree(B);
-
-    return(myPoly);
+    return true;
+
+escape_LUD:
+    // XXX drop the LUD version!
+    psFree(A);
+    psFree(B);
+    return false;
+
+escape_GJ:
+    psFree(A);
+    psFree(B);
+    return false;
 }
 
@@ -459 +458 @@
 x and fErr vectors may be NULL.  All non-NULL vectors must be of type
 PS_TYPE_F64.
- *****************************************************************************/
-static psPolynomial1D* VectorFitPolynomial1DOrd(
+ 
+XXX EAM : since this is a private function, can we drop the ASSERTS?
+XXX EAM : can we drop the LUD version? it does not calculate coeffErr values!!
+*****************************************************************************/
+static bool VectorFitPolynomial1DOrd(
     psPolynomial1D* myPoly,
     const psVector *mask,
@@ -469 +471 @@
 {
     psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL);
+    PS_ASSERT_POLY_NON_NULL(myPoly, false);
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false);
     if (mask) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
     }
     if (x) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-        PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false);
+        PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false);
     }
     if (fErr) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, false);
     }
 
@@ -510 +512 @@
     if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) {
         psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B.  Returning NULL.\n");
-        psFree(myPoly);
         psFree(A);
         psFree(B);
         psTrace("psLib.math", 4, "---- %s() End ----\n", __func__);
-        return(NULL);
+        return false;
     }
 
@@ -596 +597 @@
         if (!psMatrixGJSolve(A, B)) {
             psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
-            psFree(myPoly);
-            myPoly = NULL;
+            goto escape_GJ;
         } else {
             // the first nTerm entries in B correspond directly to the desired
@@ -617 +617 @@
         if (ALUD == NULL) {
             psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
-            psFree(myPoly);
-            myPoly = NULL;
+            goto escape_LUD;
         } else {
             coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
             if (coeffs == NULL) {
                 psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
-                psFree(myPoly);
-                myPoly = NULL;
+                goto escape_LUD;
             } else {
                 for (psS32 k = 0; k < nTerm; k++) {
@@ -632 +630 @@
             }
         }
-
         psFree(ALUD);
         psFree(coeffs);
         psFree(outPerm);
     }
-
-
     psFree(A);
     psFree(B);
 
     psTrace("psLib.math", 4, "---- %s() End ----\n", __func__);
-    return (myPoly);
+    return true;
+
+escape_LUD:
+    psFree(A);
+    psFree(B);
+    return false;
+
+escape_GJ:
+    psFree(A);
+    psFree(B);
+    return false;
 }
 
@@ -653 +658 @@
 conversion only.
  *****************************************************************************/
-psPolynomial1D *psVectorFitPolynomial1D(
+bool psVectorFitPolynomial1D(
     psPolynomial1D *poly,
     const psVector *mask,
@@ -661 +666 @@
     const psVector *x)
 {
-    // Internal pointers for possibly NULL or mis-typed vectors.
+    PS_ASSERT_POLY_NON_NULL(poly, false);
+    PS_ASSERT_INT_NONNEGATIVE(poly->nX, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_NON_EMPTY(f, false);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
+    }
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false);
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, false);
+    }
+    if (x != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false);
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, false);
+    }
+
+    // Convert input vectors to F64 if necessary.
+    psVector *f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy (NULL, f, PS_TYPE_F64);
     psVector *x64 = NULL;
-    psVector *f64 = NULL;
+    if (x != NULL) {
+        x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy (NULL, x, PS_TYPE_F64);
+    }
     psVector *fErr64 = NULL;
-
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    //PS_ASSERT_INT_NONNEGATIVE(poly->nX, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_NON_EMPTY(f, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
-    if (mask != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    if (x != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
-    }
     if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
-    }
-
-    f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy(NULL, f, PS_TYPE_F64);
-
-    if (x != NULL) {
-        x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy(NULL, x, PS_TYPE_F64);
-    }
-
-    if (fErr != NULL) {
-        fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy(NULL, fErr, PS_TYPE_F64);
-    }
+        fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy (NULL, fErr, PS_TYPE_F64);
+    }
+
+    bool result = true;
 
     switch (poly->type) {
     case PS_POLYNOMIAL_ORD:
-        poly = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64);
-        if (poly == NULL) {
+        result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64);
+        if (!result) {
             psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning NULL.\n");
         }
@@ -703 +707 @@
     case PS_POLYNOMIAL_CHEB:
         if (mask != NULL) {
-            //            psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n");
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n");
         }
         if (fErr != NULL) {
@@ -713 +717 @@
         }
 
-        poly = vectorFitPolynomial1DCheb(poly, mask, maskValue, f64, fErr64, x64);
+        result = vectorFitPolynomial1DCheb(poly, mask, maskValue, f64, fErr64, x64);
+        if (!result) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning NULL.\n");
+        }
 
         if (x == NULL) {
@@ -721 +728 @@
     default:
         psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type (%d).  Returning NULL.\n", poly->type);
-        poly = NULL;
+        result = false;
         break;
     }
 
     // Free psVectors that were created for NULL arguments.
-    if (f->type.type != PS_TYPE_F64) {
-        psFree(f64);
-    }
-
-    if ((x != NULL) && (x->type.type != PS_TYPE_F64)) {
-        psFree(x64);
-    }
-
-    if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
-        psFree(fErr64);
-    }
-
-    return(poly);
+    PS_FREE_TEMP_F64_VECTOR (f, f64);
+    PS_FREE_TEMP_F64_VECTOR (x, x64);
+    PS_FREE_TEMP_F64_VECTOR (fErr, fErr64);
+
+    return result;
 }
 
 // This function accepts F32 and F64 input vectors.
-psPolynomial1D *psVectorClipFitPolynomial1D(
+bool psVectorClipFitPolynomial1D(
     psPolynomial1D *poly,
     psStats *stats,
@@ -752 +751 @@
 {
     psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_PTR_NON_NULL(stats, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(mask, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, NULL);
-    PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    PS_ASSERT_POLY_NON_NULL(poly, false);
+    PS_ASSERT_PTR_NON_NULL(stats, false);
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
+    PS_ASSERT_VECTOR_NON_NULL(mask, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, false);
+    PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
+
     if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, false);
+        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false);
     }
 
@@ -767 +767 @@
     psVector *x = NULL;
     if (xIn != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, NULL);
-        PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, false);
+        PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, false);
         x = (psVector *) xIn;
     } else {
@@ -781 +781 @@
         } else {
             psError(PS_ERR_UNKNOWN, true, "Error, bad poly type.\n");
-            return(NULL);
-        }
+            return false;
+        }
+    }
+
+    // the user supplies one of various stats option pairs,
+    // determine the desired mean and stdev STATS options:
+    // XXX enforce consistency?
+    // XXX psStatsGetValue() probably has inverted precedence
+    psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2);
+    psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2);
+    if (!meanOption) {
+        psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");
+        return false;
+    }
+    if (!stdevOption) {
+        psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected");
+        return false;
     }
 
@@ -798 +813 @@
         minClipSigma = fabs(stats->clipSigma);
     }
-
-    psVector *fit   = NULL;
     psVector *resid = psVectorAlloc(f->n, PS_TYPE_F64);
 
-    // eventual expansion: user supplies one of various stats option pairs,
-    // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to
-    // evaluate the clipping sigma
-    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
-    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-    stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter);
     psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
@@ -822 +829 @@
             }
         }
-        poly = psVectorFitPolynomial1D(poly, mask, maskValue, f, fErr, x);
-        if (poly == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning NULL.\n");
+        if (!psVectorFitPolynomial1D(poly, mask, maskValue, f, fErr, x)) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
             if (xIn == NULL) {
                 psFree(x);
             }
-            return(NULL);
-        }
-
-        fit = psPolynomial1DEvalVector(poly, x);
+            return false;
+        }
+
+        psVector *fit = psPolynomial1DEvalVector(poly, x);
         if (fit == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector().  Returning NULL.\n");
+            psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector().  Returning false.\n");
             psFree(resid);
-            return(NULL);
+            return false;
         }
         for (psS32 i = 0 ; i < f->n ; i++) {
@@ -856 +862 @@
 
         if (!psVectorStats(stats, resid, NULL, mask, maskValue)) {
-            psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector.  Returning NULL.\n");
-            psFree(resid)
-            psFree(fit)
-            return(NULL);
-        }
-        # if (USE_ROBUST_STATS_FOR_CLIPPING)
-            psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian);
-        psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev);
-        psF32 minClipValue = -minClipSigma*stats->robustStdev;
-        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
-        psF32 clipMedian = stats->robustMedian;
-        # else
-
-            psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian);
-        psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev);
-        psF32 minClipValue = -minClipSigma*stats->robustStdev;
-        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
-        psF32 clipMedian = stats->sampleMedian;
-        # endif
+            psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector.  Returning false.\n");
+            psFree(resid);
+            psFree(fit);
+            return false;
+        }
+
+        double meanValue = psStatsGetValue (stats, meanOption);
+        double stdevValue = psStatsGetValue (stats, stdevOption);
+
+        psTrace("psLib.math", 5, "Mean is %f\n", meanValue);
+        psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue);
+        psF32 minClipValue = -minClipSigma*stdevValue;
+        psF32 maxClipValue = +maxClipSigma*stdevValue;
 
         // set mask if pts are not valid
@@ -884 +884 @@
             }
 
-            if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) {
+            if ((resid->data.F64[i] - meanValue > maxClipValue) || (resid->data.F64[i] - meanValue < minClipValue)) {
                 if (f->type.type == PS_TYPE_F64) {
                     psTrace("psLib.math", 6, "Masking element %d (%f).  resid->data.F64[%d] is %f\n",
@@ -906 +906 @@
         //
         psTrace("psLib.math", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n);
+        stats->clippedNvalues = Nkeep;
         psFree(fit);
     }
@@ -917 +918 @@
 
     psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
-    return (poly);
+    return true;
 }
 
@@ -933 +934 @@
  
  *****************************************************************************/
-static psPolynomial2D* VectorFitPolynomial2DOrd(
+static bool VectorFitPolynomial2DOrd(
     psPolynomial2D* myPoly,
     const psVector* mask,
@@ -943 +944 @@
 {
     psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL);
+    PS_ASSERT_POLY_NON_NULL(myPoly, false);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, false);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, false);
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false);
     if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL);
-    }
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, false);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, false);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(x, false);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false);
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false);
     if (mask != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, false);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
     }
 
@@ -974 +975 @@
     if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) {
         psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B.  Returning NULL.\n");
-        psFree(myPoly);
         psFree(A);
         psFree(B);
         psTrace("psLib.math", 4, "---- %s() End ----\n", __func__);
-        return(NULL);
+        return false;
     }
 
@@ -1040 +1040 @@
     if (!psMatrixGJSolve(A, B)) {
         psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
-        psFree(myPoly);
-        myPoly = NULL;
-    } else {
-        // select the appropriate solution entries
-        for (int i = 0; i < nTerm; i++) {
-            int l = i / nYterm;         // x index
-            int m = i % nYterm;         // y index
-            myPoly->coeff[l][m] = B->data.F64[i];
-            myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]);
-        }
-    }
-
+        psFree(A);
+        psFree(B);
+        return false;
+    }
+
+    // select the appropriate solution entries
+    for (int i = 0; i < nTerm; i++) {
+        int l = i / nYterm;         // x index
+        int m = i % nYterm;         // y index
+        myPoly->coeff[l][m] = B->data.F64[i];
+        myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]);
+    }
     psFree(A);
     psFree(B);
 
     psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
-    return (myPoly);
+    return true;
 }
 
@@ -1065 +1065 @@
 vector conversion only.
  *****************************************************************************/
-psPolynomial2D *psVectorFitPolynomial2D(
+bool psVectorFitPolynomial2D(
     psPolynomial2D *poly,
     const psVector *mask,
@@ -1074 +1074 @@
     const psVector *y)
 {
-    // Internal pointers for possibly NULL or mis-typed vectors.
-    psVector *x64 = NULL;
-    psVector *y64 = NULL;
-    psVector *f64 = NULL;
+    PS_ASSERT_POLY_NON_NULL(poly, false);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
+    PS_ASSERT_VECTOR_NON_NULL(x, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false);
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
+    }
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false);
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, false);
+    }
+
+    // Convert input vectors to F64 if necessary.
+    psVector *f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy(NULL, f, PS_TYPE_F64);
+    psVector *x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy(NULL, x, PS_TYPE_F64);
+    psVector *y64 = (y->type.type == PS_TYPE_F64) ? (psVector *) y : psVectorCopy(NULL, y, PS_TYPE_F64);
+
     psVector *fErr64 = NULL;
-
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
-
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
-
-    if (mask != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
     if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
-    }
-
-    //
-    // Convert input vectors to F64 if necessary.
-    //
-    if (f->type.type != PS_TYPE_F64) {
-        f64 = psVectorCopy(NULL, f, PS_TYPE_F64);
-    } else {
-        f64 = (psVector *) f;
-    }
-
-    if (x->type.type != PS_TYPE_F64) {
-        x64 = psVectorCopy(NULL, x, PS_TYPE_F64);
-    } else {
-        x64 = (psVector *) x;
-    }
-
-    if (y->type.type != PS_TYPE_F64) {
-        y64 = psVectorCopy(NULL, y, PS_TYPE_F64);
-    } else {
-        y64 = (psVector *) y;
-    }
-
-    //
-    // fErr
-    //
-    if (fErr != NULL) {
-        if (fErr->type.type != PS_TYPE_F64) {
-            fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64);
-        } else {
-            fErr64 = (psVector *) fErr;
-        }
-    }
-
-    if (poly->type == PS_POLYNOMIAL_ORD) {
-        poly = VectorFitPolynomial2DOrd(poly, mask, maskValue, f64, fErr64, x64, y64);
-        if (poly == NULL) {
+        fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy(NULL, fErr, PS_TYPE_F64);
+    }
+
+    bool result = true;
+
+    switch (poly->type) {
+    case PS_POLYNOMIAL_ORD:
+        result = VectorFitPolynomial2DOrd(poly, mask, maskValue, f64, fErr64, x64, y64);
+        if (!result) {
             psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial.  Returning NULL.\n");
-            // Free psVectors that were created for NULL arguments.
-            if (f->type.type != PS_TYPE_F64) {
-                psFree(f64);
-            }
-
-            if (x->type.type != PS_TYPE_F64) {
-                psFree(x64);
-            }
-
-            if (y->type.type != PS_TYPE_F64) {
-                psFree(y64);
-            }
-
-            if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
-                psFree(fErr64);
-            }
-            return(NULL);
-        }
-    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        }
+        break;
+    case PS_POLYNOMIAL_CHEB:
         if (mask != NULL) {
             psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n");
         }
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
-        psFree(poly);
-        poly = NULL;
-    } else {
-        // Free psVectors that were created for NULL arguments.
-        if (f->type.type != PS_TYPE_F64) {
-            psFree(f64);
-        }
-
-        if (x->type.type != PS_TYPE_F64) {
-            psFree(x64);
-        }
-
-        if (y->type.type != PS_TYPE_F64) {
-            psFree(y64);
-        }
-
-        if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
-            psFree(fErr64);
-        }
+        psError(PS_ERR_UNKNOWN, true, "2-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
+        result = false;
+        break;
+    default:
         psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type.  Returning NULL.\n");
-    }
-
+        result = false;
+        break;
+    }
 
     // Free psVectors that were created for NULL arguments.
-    if (f->type.type != PS_TYPE_F64) {
-        psFree(f64);
-    }
-
-    if (x->type.type != PS_TYPE_F64) {
-        psFree(x64);
-    }
-
-    if (y->type.type != PS_TYPE_F64) {
-        psFree(y64);
-    }
-
-    if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
-        psFree(fErr64);
-    }
-
-    return(poly);
+    PS_FREE_TEMP_F64_VECTOR (f, f64);
+    PS_FREE_TEMP_F64_VECTOR (x, x64);
+    PS_FREE_TEMP_F64_VECTOR (y, y64);
+    PS_FREE_TEMP_F64_VECTOR (fErr, fErr64);
+
+    return result;
 }
 
-psPolynomial2D *psVectorClipFitPolynomial2D(
+bool psVectorClipFitPolynomial2D(
     psPolynomial2D *poly,
     psStats *stats,
@@ -1212 +1144 @@
 {
     psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
-    PS_ASSERT_PTR_NON_NULL(stats, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(mask, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
-    if (mask != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-    PS_ASSERT_VECTOR_TYPE(x, f->type.type, NULL);
-
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_TYPE(y, f->type.type, NULL);
+    PS_ASSERT_POLY_NON_NULL(poly, false);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false);
+    PS_ASSERT_PTR_NON_NULL(stats, false);
+    PS_ASSERT_VECTOR_NON_NULL(mask, false);
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(x, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false);
+    PS_ASSERT_VECTOR_TYPE(x, f->type.type, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false);
+    PS_ASSERT_VECTOR_TYPE(y, f->type.type, false);
+
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false);
+    PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
 
     if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false);
+        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false);
+    }
+
+    // the user supplies one of various stats option pairs,
+    // determine the desired mean and stdev STATS options:
+    // XXX enforce consistency?
+    // XXX psStatsGetValue() probably has inverted precedence
+    psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2);
+    psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2);
+    if (!meanOption) {
+        psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");
+        return false;
+    }
+    if (!stdevOption) {
+        psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected");
+        return false;
     }
 
@@ -1250 +1197 @@
     psVector *resid = psVectorAlloc(f->n, PS_TYPE_F64);
 
-    // eventual expansion: user supplies one of various stats option pairs,
-    // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to
-    // evaluate the clipping sigma
-    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
-    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-    stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter);
     psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
@@ -1270 +1211 @@
         }
 
-        poly = psVectorFitPolynomial2D(poly, mask, maskValue, f, fErr, x, y);
-        if (poly == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data.  Returning NULL.\n");
+        if (!psVectorFitPolynomial2D(poly, mask, maskValue, f, fErr, x, y)) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data.  Returning false.\n");
             psFree(resid)
-            return(NULL);
+            return false;
         }
 
@@ -1281 +1221 @@
             psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector().  Returning NULL.\n");
             psFree(resid)
-            return(NULL);
+            return false;
         }
 
@@ -1307 +1247 @@
             psFree(resid)
             psFree(fit)
-            return(NULL);
-        }
-        # if (USE_ROBUST_STATS_FOR_CLIPPING)
-            psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian);
-        psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev);
-        psTrace("psLib.math", 5, "Sample Median is %f\n", stats->sampleMedian);
-        psTrace("psLib.math", 5, "Sample Stdev is %f\n", stats->sampleStdev);
-        psF32 minClipValue = -minClipSigma*stats->robustStdev;
-        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
-        psF32 clipMedian = stats->robustMedian;
-        # else
-
-            psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian);
-        psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev);
-        psTrace("psLib.math", 5, "Robust Median is %f\n", stats->robustMedian);
-        psTrace("psLib.math", 5, "Robust Stdev is %f\n", stats->robustStdev);
-        psF32 minClipValue = -minClipSigma*stats->robustStdev;
-        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
-        psF32 clipMedian = stats->sampleMedian;
-        # endif
+            return false;
+        }
+
+        double meanValue = psStatsGetValue (stats, meanOption);
+        double stdevValue = psStatsGetValue (stats, stdevOption);
+
+        psTrace("psLib.math", 5, "Mean is %f\n", meanValue);
+        psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue);
+        psF32 minClipValue = -minClipSigma*stdevValue;
+        psF32 maxClipValue = +maxClipSigma*stdevValue;
 
         // set mask if pts are not valid
@@ -1336 +1266 @@
             }
 
-            if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) {
+            if ((resid->data.F64[i] - meanValue > maxClipValue) || (resid->data.F64[i] - meanValue < minClipValue)) {
                 if (fit->type.type == PS_TYPE_F64) {
                     psTrace("psLib.math", 6, "Masking element %d (%f).  resid->data.F64[%d] is %f\n",
@@ -1352 +1282 @@
             Nkeep++;
         }
-
         psTrace("psLib.math", 4, "keeping %d of %ld pts for fit\n", Nkeep, x->n);
+        stats->clippedNvalues = Nkeep;
         psFree(fit);
     }
@@ -1359 +1289 @@
     psFree(resid);
 
-    if (poly == NULL) {
-        psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data.  Returning NULL.\n");
-        return(NULL);
-    }
-
     psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
-    return(poly);
+    return true;
 }
 
@@ -1381 +1306 @@
  
  *****************************************************************************/
-static psPolynomial3D* VectorFitPolynomial3DOrd(
+static bool VectorFitPolynomial3DOrd(
     psPolynomial3D* myPoly,
     const psVector* mask,
@@ -1392 +1317 @@
 {
     psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL);
-
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL);
+    PS_ASSERT_POLY_NON_NULL(myPoly, false);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, false);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, false);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false);
     if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL);
-    }
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
-    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, false);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, false);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(x, false);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false);
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false);
+    PS_ASSERT_VECTOR_NON_NULL(z, false);
+    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false);
     if (mask != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
+    }
 
     int nXterm = 1 + myPoly->nX;        // Number of x terms
@@ -1429 +1353 @@
     if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) {
         psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B.  Returning NULL.\n");
-        psFree(myPoly);
         psFree(A);
         psFree(B);
         psTrace("psLib.math", 4, "---- %s() End ----\n", __func__);
-        return(NULL);
+        return false;
     }
 
@@ -1511 +1434 @@
         // The matrices were overflowing, so I switched to LUD.
         if (!psMatrixGJSolve(A, B)) {
-            psFree(A);
-            psFree(B);
             psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n");
-            return(NULL);
+            goto escape_GJ;
         }
         // select the appropriate solution entries
@@ -1534 +1455 @@
         if (ALUD == NULL) {
             psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
-            psFree(myPoly);
-            myPoly = NULL;
+            goto escape_LUD;
         } else {
             coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
             if (coeffs == NULL) {
                 psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
-                psFree(myPoly);
-                myPoly = NULL;
+                goto escape_LUD;
             } else {
                 // select the appropriate solution entries
@@ -1565 +1484 @@
 
     psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
-    return (myPoly);
+    return true;
+
+escape_LUD:
+    psFree(A);
+    psFree(B);
+    return false;
+
+escape_GJ:
+
+    psFree(A);
+    psFree(B);
+    return false;
 }
 
@@ -1574 +1504 @@
 vector conversion only.
  *****************************************************************************/
-psPolynomial3D *psVectorFitPolynomial3D(
+bool psVectorFitPolynomial3D(
     psPolynomial3D *poly,
     const psVector *mask,
@@ -1584 +1514 @@
     const psVector *z)
 {
-    // Internal pointers for possibly NULL or mis-typed vectors.
-    psVector *x64 = NULL;
-    psVector *y64 = NULL;
-    psVector *z64 = NULL;
-    psVector *f64 = NULL;
+    PS_ASSERT_POLY_NON_NULL(poly, false);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
+    PS_ASSERT_VECTOR_NON_NULL(x, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false);
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false);
+    PS_ASSERT_VECTOR_NON_NULL(z, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
+    }
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false);
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, false);
+    }
+
+    // Convert input vectors to F64 if necessary.
+    psVector *f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy(NULL, f, PS_TYPE_F64);
+    psVector *x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy(NULL, x, PS_TYPE_F64);
+    psVector *y64 = (y->type.type == PS_TYPE_F64) ? (psVector *) y : psVectorCopy(NULL, y, PS_TYPE_F64);
+    psVector *z64 = (z->type.type == PS_TYPE_F64) ? (psVector *) z : psVectorCopy(NULL, z, PS_TYPE_F64);
+
     psVector *fErr64 = NULL;
-
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
-
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
-    if (mask != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
     if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
-        //        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL);
-    }
-
-    //
-    // Convert input vectors to F64 if necessary.
-    //
-    if (f->type.type != PS_TYPE_F64) {
-        f64 = psVectorCopy(NULL, f, PS_TYPE_F64);
-    } else {
-        f64 = (psVector *) f;
-    }
-    if (x->type.type != PS_TYPE_F64) {
-        x64 = psVectorCopy(NULL, x, PS_TYPE_F64);
-    } else {
-        x64 = (psVector *) x;
-    }
-    if (y->type.type != PS_TYPE_F64) {
-        y64 = psVectorCopy(NULL, y, PS_TYPE_F64);
-    } else {
-        y64 = (psVector *) y;
-    }
-
-    if (z->type.type != PS_TYPE_F64) {
-        z64 = psVectorCopy(NULL, z, PS_TYPE_F64);
-    } else {
-        z64 = (psVector *) z;
-    }
-
-    if (fErr != NULL) {
-        if (fErr->type.type != PS_TYPE_F64) {
-            fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64);
-        } else {
-            fErr64 = (psVector *) fErr;
-        }
-    }
-
-    if (poly->type == PS_POLYNOMIAL_ORD) {
-        poly = VectorFitPolynomial3DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64);
-        if (poly == NULL) {
+        fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy(NULL, fErr, PS_TYPE_F64);
+    }
+
+    bool result = true;
+
+    switch (poly->type) {
+    case PS_POLYNOMIAL_ORD:
+        result = VectorFitPolynomial3DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64);
+        if (!result) {
             psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial.  Returning NULL.\n");
-            // Free psVectors that were created for NULL arguments.
-            if (f->type.type != PS_TYPE_F64) {
-                psFree(f64);
-            }
-
-            if (x->type.type != PS_TYPE_F64) {
-                psFree(x64);
-            }
-
-            if (y->type.type != PS_TYPE_F64) {
-                psFree(y64);
-            }
-
-            if (z->type.type != PS_TYPE_F64) {
-                psFree(z64);
-            }
-
-            if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
-                psFree(fErr64);
-            }
-            return(NULL);
-        }
-    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        }
+        break;
+    case PS_POLYNOMIAL_CHEB:
         if (mask != NULL) {
             psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n");
         }
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
-        psFree(poly);
-        poly = NULL;
-    } else {
-        // Free psVectors that were created for NULL arguments.
-        if (f->type.type != PS_TYPE_F64) {
-            psFree(f64);
-        }
-
-        if (x->type.type != PS_TYPE_F64) {
-            psFree(x64);
-        }
-
-        if (y->type.type != PS_TYPE_F64) {
-            psFree(y64);
-        }
-
-        if (z->type.type != PS_TYPE_F64) {
-            psFree(z64);
-        }
-
-        if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
-            psFree(fErr64);
-        }
+        psError(PS_ERR_UNKNOWN, true, "3-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
+        result = false;
+        break;
+    default:
         psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type.  Returning NULL.\n");
-    }
-
+        result = false;
+        break;
+    }
 
     // Free psVectors that were created for NULL arguments.
-    if (f->type.type != PS_TYPE_F64) {
-        psFree(f64);
-    }
-
-    if (x->type.type != PS_TYPE_F64) {
-        psFree(x64);
-    }
-
-    if (y->type.type != PS_TYPE_F64) {
-        psFree(y64);
-    }
-
-    if (z->type.type != PS_TYPE_F64) {
-        psFree(z64);
-    }
-
-    if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
-        psFree(fErr64);
-    }
-
-    return(poly);
+    PS_FREE_TEMP_F64_VECTOR (f, f64);
+    PS_FREE_TEMP_F64_VECTOR (x, x64);
+    PS_FREE_TEMP_F64_VECTOR (y, y64);
+    PS_FREE_TEMP_F64_VECTOR (z, z64);
+    PS_FREE_TEMP_F64_VECTOR (fErr, fErr64);
+
+    return result;
 }
 
-psPolynomial3D *psVectorClipFitPolynomial3D(
+bool psVectorClipFitPolynomial3D(
     psPolynomial3D *poly,
     psStats *stats,
@@ -1738 +1589 @@
 {
     psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
-    PS_ASSERT_PTR_NON_NULL(stats, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(mask, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
-    if (mask != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-    PS_ASSERT_VECTOR_TYPE(x, f->type.type, NULL);
-
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_TYPE(y, f->type.type, NULL);
-
-    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
-    PS_ASSERT_VECTOR_TYPE(z, f->type.type, NULL);
+    PS_ASSERT_POLY_NON_NULL(poly, false);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false);
+    PS_ASSERT_PTR_NON_NULL(stats, false);
+    PS_ASSERT_VECTOR_NON_NULL(mask, false);
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(x, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false);
+    PS_ASSERT_VECTOR_TYPE(x, f->type.type, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false);
+    PS_ASSERT_VECTOR_TYPE(y, f->type.type, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(z, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false);
+    PS_ASSERT_VECTOR_TYPE(z, f->type.type, false);
+
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false);
+    PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
 
     if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false);
+        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false);
+    }
+
+    // the user supplies one of various stats option pairs,
+    // determine the desired mean and stdev STATS options:
+    // XXX enforce consistency?
+    // XXX psStatsGetValue() probably has inverted precedence
+    psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2);
+    psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2);
+    if (!meanOption) {
+        psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");
+        return false;
+    }
+    if (!stdevOption) {
+        psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected");
+        return false;
     }
 
@@ -1778 +1644 @@
         minClipSigma = fabs(stats->clipSigma);
     }
-    psVector *fit   = NULL;
     psVector *resid = psVectorAlloc(f->n, PS_TYPE_F64);
 
-    // eventual expansion: user supplies one of various stats option pairs,
-    // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to
-    // evaluate the clipping sigma
-    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
-    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-    stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter);
     psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
@@ -1801 +1660 @@
         }
 
-        poly = psVectorFitPolynomial3D(poly, mask, maskValue, f, fErr, x, y, z);
-        if (poly == NULL) {
+        if (!psVectorFitPolynomial3D(poly, mask, maskValue, f, fErr, x, y, z)) {
             psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data.  Returning NULL.\n");
             psFree(resid)
-            psFree(fit)
-            return(NULL);
-        }
-        fit = psPolynomial3DEvalVector(poly, x, y, z);
+            return false;
+        }
+        psVector *fit = psPolynomial3DEvalVector(poly, x, y, z);
         if (fit == NULL) {
             psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector().  Returning NULL.\n");
             psFree(resid)
-            return(NULL);
+            return false;
         }
         for (psS32 i = 0 ; i < f->n ; i++) {
@@ -1837 +1694 @@
             psFree(resid)
             psFree(fit)
-            return(NULL);
-        }
-
-        # if (USE_ROBUST_STATS_FOR_CLIPPING)
-            psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian);
-        psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev);
-        psF32 minClipValue = -minClipSigma*stats->robustStdev;
-        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
-        psF32 clipMedian = stats->robustMedian;
-        # else
-
-            psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian);
-        psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev);
-        psF32 minClipValue = -minClipSigma*stats->robustStdev;
-        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
-        psF32 clipMedian = stats->sampleMedian;
-        # endif
+            return false;
+        }
+
+        double meanValue = psStatsGetValue (stats, meanOption);
+        double stdevValue = psStatsGetValue (stats, stdevOption);
+
+        psTrace("psLib.math", 5, "Mean is %f\n", meanValue);
+        psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue);
+        psF32 minClipValue = -minClipSigma*stdevValue;
+        psF32 maxClipValue = +maxClipSigma*stdevValue;
 
         // set mask if pts are not valid
@@ -1863 +1713 @@
             }
 
-            if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue))  {
+            if ((resid->data.F64[i] - meanValue > maxClipValue) || (resid->data.F64[i] - meanValue < minClipValue))  {
                 if (f->type.type == PS_TYPE_F64) {
                     psTrace("psLib.math", 6, "Masking element %d (%f).  resid->data.F64[%d] is %f\n",
@@ -1879 +1729 @@
             Nkeep++;
         }
-
         psTrace("psLib.math", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n);
+        stats->clippedNvalues = Nkeep;
         psFree(fit);
     }
@@ -1886 +1736 @@
     psFree(resid);
 
-    if (poly == NULL) {
-        psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data.  Returning NULL.\n");
-        return(NULL);
-    }
-
     psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
-    return(poly);
+    return true;
 }
 
@@ -1906 +1751 @@
  
  *****************************************************************************/
-static psPolynomial4D* VectorFitPolynomial4DOrd(
+static bool VectorFitPolynomial4DOrd(
     psPolynomial4D* myPoly,
     const psVector* mask,
@@ -1918 +1763 @@
 {
     psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(myPoly->nT, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL);
+    PS_ASSERT_POLY_NON_NULL(myPoly, false);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, false);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, false);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, false);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nT, false);
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false);
     if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL);
-    }
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
-    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(t, NULL);
-    PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F64, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, false);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, false);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(x, false);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false);
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false);
+    PS_ASSERT_VECTOR_NON_NULL(z, false);
+    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false);
+    PS_ASSERT_VECTOR_NON_NULL(t, false);
+    PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F64, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, false);
     if (mask) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, false);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
     }
 
@@ -1959 +1804 @@
     if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) {
         psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B.  Returning NULL.\n");
-        psFree(myPoly);
         psFree(A);
         psFree(B);
         psTrace("psLib.math", 4, "---- %s() End ----\n", __func__);
-        return(NULL);
+        return false;
     }
 
@@ -2051 +1895 @@
         // The GaussJordan version was overflowing, so I'm using LUD.
         if (!psMatrixGJSolve(A, B)) {
-            psFree(A);
-            psFree(B);
             psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n");
-            return(NULL);
+            goto escape_GJ;
         }
 
@@ -2076 +1918 @@
         if (ALUD == NULL) {
             psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
-            psFree(myPoly);
-            myPoly = NULL;
+            goto escape_LUD;
         } else {
             coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
             if (coeffs == NULL) {
                 psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
-                psFree(myPoly);
-                myPoly = NULL;
+                goto escape_LUD;
             } else {
                 // select the appropriate solution entries
@@ -2100 +1940 @@
             }
         }
-
         psFree(ALUD);
         psFree(coeffs);
         psFree(outPerm);
-
-    }
-
+    }
     psFree(A);
     psFree(B);
 
     psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
-    return (myPoly);
+    return true;
+
+escape_LUD:
+    psFree(A);
+    psFree(B);
+    return false;
+
+escape_GJ:
+    psFree(A);
+    psFree(B);
+    return false;
 }
 
@@ -2120 +1967 @@
 via vector conversion only.
  *****************************************************************************/
-psPolynomial4D *psVectorFitPolynomial4D(
+bool psVectorFitPolynomial4D(
     psPolynomial4D *poly,
     const psVector *mask,
@@ -2131 +1978 @@
     const psVector *t)
 {
-    // Internal pointers for possibly NULL or mis-typed vectors.
-    psVector *x64 = NULL;
-    psVector *y64 = NULL;
-    psVector *z64 = NULL;
-    psVector *t64 = NULL;
-    psVector *f64 = NULL;
+    PS_ASSERT_POLY_NON_NULL(poly, false);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
+    PS_ASSERT_VECTOR_NON_NULL(x, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false);
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false);
+    PS_ASSERT_VECTOR_NON_NULL(z, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false);
+    PS_ASSERT_VECTOR_NON_NULL(t, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, false);
+    if (mask) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
+    }
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false);
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, false);
+    }
+
+    // Convert input vectors to F64 if necessary.
+    psVector *f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy(NULL, f, PS_TYPE_F64);
+    psVector *x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy(NULL, x, PS_TYPE_F64);
+    psVector *y64 = (y->type.type == PS_TYPE_F64) ? (psVector *) y : psVectorCopy(NULL, y, PS_TYPE_F64);
+    psVector *z64 = (z->type.type == PS_TYPE_F64) ? (psVector *) z : psVectorCopy(NULL, z, PS_TYPE_F64);
+    psVector *t64 = (t->type.type == PS_TYPE_F64) ? (psVector *) t : psVectorCopy(NULL, t, PS_TYPE_F64);
+
     psVector *fErr64 = NULL;
-
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
-
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
-    if (mask) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(t, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
     if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
-    }
-
-
-    //
-    // Convert input vector to F64 if necessary.
-    //
-    if (f->type.type != PS_TYPE_F64) {
-        f64 = psVectorCopy(NULL, f, PS_TYPE_F64);
-    } else {
-        f64 = (psVector *) f;
-    }
-    if (x->type.type != PS_TYPE_F64) {
-        x64 = psVectorCopy(NULL, x, PS_TYPE_F64);
-    } else {
-        x64 = (psVector *) x;
-    }
-    if (y->type.type != PS_TYPE_F64) {
-        y64 = psVectorCopy(NULL, y, PS_TYPE_F64);
-    } else {
-        y64 = (psVector *) y;
-    }
-    if (z->type.type != PS_TYPE_F64) {
-        z64 = psVectorCopy(NULL, z, PS_TYPE_F64);
-    } else {
-        z64 = (psVector *) z;
-    }
-    if (t->type.type != PS_TYPE_F64) {
-        t64 = psVectorCopy(NULL, t, PS_TYPE_F64);
-    } else {
-        t64 = (psVector *) t;
-    }
-    //
-    // fErr
-    //
-    if (fErr != NULL) {
-        if (fErr->type.type != PS_TYPE_F64) {
-            fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64);
-        } else {
-            fErr64 = (psVector *) fErr;
-        }
-    }
-
-    if (poly->type == PS_POLYNOMIAL_ORD) {
-        poly = VectorFitPolynomial4DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64, t64);
-        if (poly == NULL) {
+        fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy(NULL, fErr, PS_TYPE_F64);
+    }
+
+    bool result = true;
+
+    switch (poly->type) {
+    case PS_POLYNOMIAL_ORD:
+        result = VectorFitPolynomial4DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64, t64);
+        if (!result) {
             psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial.  Returning NULL.\n");
-            // Free psVectors that were created for NULL arguments.
-            if (f->type.type != PS_TYPE_F64) {
-                psFree(f64);
-            }
-
-            if (x->type.type != PS_TYPE_F64) {
-                psFree(x64);
-            }
-
-            if (y->type.type != PS_TYPE_F64) {
-                psFree(y64);
-            }
-
-            if (z->type.type != PS_TYPE_F64) {
-                psFree(z64);
-            }
-
-            if (t->type.type != PS_TYPE_F64) {
-                psFree(t64);
-            }
-
-            if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
-                psFree(fErr64);
-            }
-            return(NULL);
-        }
-    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        }
+        break;
+    case PS_POLYNOMIAL_CHEB:
         if (mask != NULL) {
             psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n");
         }
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
-        psFree(poly);
-        poly = NULL;
-    } else {
-        // Free psVectors that were created for NULL arguments.
-        if (f->type.type != PS_TYPE_F64) {
-            psFree(f64);
-        }
-
-        if (x->type.type != PS_TYPE_F64) {
-            psFree(x64);
-        }
-
-        if (y->type.type != PS_TYPE_F64) {
-            psFree(y64);
-        }
-
-        if (z->type.type != PS_TYPE_F64) {
-            psFree(z64);
-        }
-
-        if (t->type.type != PS_TYPE_F64) {
-            psFree(t64);
-        }
-
-        if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
-            psFree(fErr64);
-        }
+        psError(PS_ERR_UNKNOWN, true, "4-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
+        result = false;
+        break;
+    default:
         psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type.  Returning NULL.\n");
-    }
-
+        result = false;
+        break;
+    }
 
     // Free psVectors that were created for NULL arguments.
-    if (f->type.type != PS_TYPE_F64) {
-        psFree(f64);
-    }
-
-    if (x->type.type != PS_TYPE_F64) {
-        psFree(x64);
-    }
-
-    if (y->type.type != PS_TYPE_F64) {
-        psFree(y64);
-    }
-
-    if (z->type.type != PS_TYPE_F64) {
-        psFree(z64);
-    }
-
-    if (t->type.type != PS_TYPE_F64) {
-        psFree(t64);
-    }
-
-    if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
-        psFree(fErr64);
-    }
-
-    return(poly);
+    PS_FREE_TEMP_F64_VECTOR (f, f64);
+    PS_FREE_TEMP_F64_VECTOR (x, x64);
+    PS_FREE_TEMP_F64_VECTOR (y, y64);
+    PS_FREE_TEMP_F64_VECTOR (z, z64);
+    PS_FREE_TEMP_F64_VECTOR (t, t64);
+    PS_FREE_TEMP_F64_VECTOR (fErr, fErr64);
+
+    return result;
 }
 
 
-psPolynomial4D *psVectorClipFitPolynomial4D(
+bool psVectorClipFitPolynomial4D(
     psPolynomial4D *poly,
     psStats *stats,
@@ -2309 +2059 @@
 {
     psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
-    PS_ASSERT_PTR_NON_NULL(stats, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(mask, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
-    if (mask) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-    PS_ASSERT_VECTOR_TYPE(x, f->type.type, NULL);
-
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_TYPE(y, f->type.type, NULL);
-
-    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
-    PS_ASSERT_VECTOR_TYPE(z, f->type.type, NULL);
-
-    PS_ASSERT_VECTOR_NON_NULL(t, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
-    PS_ASSERT_VECTOR_TYPE(t, f->type.type, NULL);
+    PS_ASSERT_POLY_NON_NULL(poly, false);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false);
+    PS_ASSERT_PTR_NON_NULL(stats, false);
+    PS_ASSERT_VECTOR_NON_NULL(mask, false);
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(x, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false);
+    PS_ASSERT_VECTOR_TYPE(x, f->type.type, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false);
+    PS_ASSERT_VECTOR_TYPE(y, f->type.type, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(z, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false);
+    PS_ASSERT_VECTOR_TYPE(z, f->type.type, false);
+
+    PS_ASSERT_VECTOR_NON_NULL(t, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, false);
+    PS_ASSERT_VECTOR_TYPE(t, f->type.type, false);
+
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false);
+    PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
 
     if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false);
+        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false);
+    }
+
+    // the user supplies one of various stats option pairs,
+    // determine the desired mean and stdev STATS options:
+    // XXX enforce consistency?
+    // XXX psStatsGetValue() probably has inverted precedence
+    psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2);
+    psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2);
+    if (!meanOption) {
+        psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");
+        return false;
+    }
+    if (!stdevOption) {
+        psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected");
+        return false;
     }
 
@@ -2353 +2118 @@
         minClipSigma = fabs(stats->clipSigma);
     }
-    psVector *fit   = NULL;
     psVector *resid = psVectorAlloc(f->n, PS_TYPE_F64);
 
-    // eventual expansion: user supplies one of various stats option pairs,
-    // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to
-    // evaluate the clipping sigma
-    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
-    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-    stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter);
     psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
@@ -2376 +2134 @@
         }
 
-        poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t);
-        if (poly == NULL) {
+        if (!psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t)) {
             psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data.  Returning NULL.\n");
             psFree(resid)
-            psFree(fit)
-            return(NULL);
-        }
-
-        fit = psPolynomial4DEvalVector (poly, x, y, z, t);
+            return false;
+        }
+
+        psVector *fit = psPolynomial4DEvalVector (poly, x, y, z, t);
         if (fit == NULL) {
             psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial4DEvalVector().  Returning NULL.\n");
             psFree(resid)
-            return(NULL);
+            return false;
         }
         for (psS32 i = 0 ; i < f->n ; i++) {
@@ -2413 +2169 @@
             psFree(resid)
             psFree(fit)
-            return(NULL);
-        }
-        # if (USE_ROBUST_STATS_FOR_CLIPPING)
-            psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian);
-        psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev);
-        psF32 minClipValue = -minClipSigma*stats->robustStdev;
-        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
-        psF32 clipMedian = stats->robustMedian;
-        # else
-
-            psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian);
-        psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev);
-        psF32 minClipValue = -minClipSigma*stats->robustStdev;
-        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
-        psF32 clipMedian = stats->sampleMedian;
-        # endif
+            return false;
+        }
+
+        double meanValue = psStatsGetValue (stats, meanOption);
+        double stdevValue = psStatsGetValue (stats, stdevOption);
+
+        psTrace("psLib.math", 5, "Mean is %f\n", meanValue);
+        psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue);
+        psF32 minClipValue = -minClipSigma*stdevValue;
+        psF32 maxClipValue = +maxClipSigma*stdevValue;
 
         // set mask if pts are not valid
@@ -2438 +2188 @@
             }
 
-            if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) {
+            if ((resid->data.F64[i] - meanValue > maxClipValue) || (resid->data.F64[i] - meanValue < minClipValue)) {
                 if (f->type.type == PS_TYPE_F64) {
                     psTrace("psLib.math", 6, "Masking element %d (%f).  resid->data.F64[%d] is %f\n",
@@ -2452 +2202 @@
                 continue;
             }
-
             Nkeep++;
         }
-
         psTrace("psLib.math", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n);
+        stats->clippedNvalues = Nkeep;
         psFree (fit);
     }
@@ -2462 +2211 @@
     psFree (resid);
 
-    if (poly == NULL) {
-        psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data.  Returning NULL.\n");
-        return(NULL);
-    }
-
     psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
-    return(poly);
+    return true;
 }