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

Timestamp:

Jan 23, 2006, 12:25:31 PM (20 years ago)

Author:

gusciora

Message:

Misc code cleaning

File:

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

trunk/psLib/src/math/psMinimizePolyFit.c

@@ -10 +10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-23 20:44:29 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-23 22:25:31 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
  *
- *  XXX: must follow coding name standards on local functions.
- *  XXX: put local functions in front.
  *  XXX: For clip-fit functions, what should we do if the mask is NULL?
- *  XXX: Generate macros for
- * PS_ASSERT_VECTOR_TYPES_MATCH()
- * PS_ASSERT_VECTOR_SIZES_MATCH()
- *  They should have better error messages.
  *
  */
@@ -42 +36 @@
 /*****************************************************************************/
 
+#define PS_VECTOR_GEN_CHEBY_INDEX(VEC, SIZE, TYPE) \
+VEC = psVectorAlloc(SIZE, TYPE); \
+if (TYPE == PS_TYPE_F64) { \
+    for (psS32 i = 0 ; i < SIZE ; i++) { \
+        VEC->data.F64[i] = ((2.0 / ((psF64) (SIZE - 1))) * ((psF64) i)) - 1.0; \
+    }\
+} else if (TYPE == PS_TYPE_F32){ \
+    for (psS32 i = 0 ; i < SIZE ; i++) { \
+        VEC->data.F32[i] = ((2.0 / ((psF32) (SIZE - 1))) * ((psF32) i)) - 1.0; \
+    }\
+}\
+
 /*****************************************************************************/
 /* TYPE DEFINITIONS                                                          */
@@ -58 +64 @@
 /*****************************************************************************/
 
-/******************************************************************************
- ******************************************************************************
- Analytical 1-D fitting routines.
- ******************************************************************************
- *****************************************************************************/
-// XXX: Make this a general type conversion macro, or function
-#define PS_VECTOR_GEN_F64_FROM_F32(VECF64, VECF32) \
-VECF64 = psVectorAlloc(VECF32->n, PS_TYPE_F64); \
-for (psS32 i = 0 ; i < VECF32->n ; i++) { \
-    VECF64->data.F64[i] = (psF64) VECF32->data.F32[i]; \
-} \
-
-// XXX: Consolidate these
-#define PS_VECTOR_GEN_CHEBY_INDEX(VECF64, SIZE) \
-VECF64 = psVectorAlloc(SIZE, PS_TYPE_F64); \
-for (psS32 i = 0 ; i < SIZE ; i++) { \
-    VECF64->data.F64[i] = ((2.0 / ((psF64) (SIZE - 1))) * ((psF64) i)) - 1.0; \
-}\
-
-#define PS_VECTOR_GEN_CHEBY_INDEX_F64(VECF64, SIZE) \
-VECF64 = psVectorAlloc(SIZE, PS_TYPE_F64); \
-for (psS32 i = 0 ; i < SIZE ; i++) { \
-    VECF64->data.F64[i] = ((2.0 / ((psF64) (SIZE - 1))) * ((psF64) i)) - 1.0; \
-}\
-
-#define PS_VECTOR_GEN_CHEBY_INDEX_F32(VECF32, SIZE) \
-VECF32 = psVectorAlloc(SIZE, PS_TYPE_F32); \
-for (psS32 i = 0 ; i < SIZE ; i++) { \
-    VECF32->data.F32[i] = ((2.0 / ((psF32) (SIZE - 1))) * ((psF32) i)) - 1.0; \
-}\
 /******************************************************************************
 BuildSums1D(sums, x, polyOrder, sums): this routine calculates the powers of
@@ -113 +89 @@
 
     xSum = 1.0;
-    for (int i = 0; i < nSum; i++) {
+    for (psS32 i = 0; i < nSum; i++) {
         sums->data.F64[i] = xSum;
         xSum *= x;
@@ -148 +124 @@
 
     xSum = 1.0;
-    for (int i = 0; i < nXsum; i++) {
+    for (psS32 i = 0; i < nXsum; i++) {
         ySum = xSum;
-        for (int j = 0; j < nYsum; j++) {
+        for (psS32 j = 0; j < nYsum; j++) {
             sums->data.F64[i][j] = ySum;
             ySum *= y;
         }
         xSum *= x;
-    }
-
-    if (0) {
-        printf("--------------------- BuildSums2D(%.2f %.2f) ---------------------\n", x, y);
-        for (int i = 0; i < nXsum; i++) {
-            for (int j = 0; j < nYsum; j++) {
-                printf("(%.2f) ", sums->data.F64[i][j]);
-            }
-            printf("\n");
-        }
     }
 
@@ -196 +162 @@
     if (sums == NULL) {
         sums = (psF64 ***) psAlloc (nXsum*sizeof(psF64));
-        for (int i = 0; i < nXsum; i++) {
+        for (psS32 i = 0; i < nXsum; i++) {
             sums[i] = (psF64 **) psAlloc (nYsum*sizeof(psF64));
-            for (int j = 0; j < nYsum; j++) {
+            for (psS32 j = 0; j < nYsum; j++) {
                 sums[i][j] = (psF64 *) psAlloc (nZsum*sizeof(psF64));
             }
@@ -207 +173 @@
     if (1) {
         zSum = 1.0;
-        for (int k = 0; k < nZsum; k++) {
+        for (psS32 k = 0; k < nZsum; k++) {
             ySum = zSum;
-            for (int j = 0; j < nYsum; j++) {
+            for (psS32 j = 0; j < nYsum; j++) {
                 xSum = ySum;
-                for (int i = 0; i < nXsum; i++) {
+                for (psS32 i = 0; i < nXsum; i++) {
                     sums[i][j][k] = xSum;
                     xSum *= x;
@@ -221 +187 @@
     } else {
         xSum = 1.0;
-        for (int i = 0; i < nXsum; i++) {
+        for (psS32 i = 0; i < nXsum; i++) {
             ySum = xSum;
-            for (int j = 0; j < nYsum; j++) {
+            for (psS32 j = 0; j < nYsum; j++) {
                 zSum = ySum;
-                for (int k = 0; k < nZsum; k++) {
+                for (psS32 k = 0; k < nZsum; k++) {
                     sums[i][j][k] = zSum;
                     zSum *= z;
@@ -240 +206 @@
 /******************************************************************************
     BuildSums4D(sums, x, y, z, t, nXterm, nYterm, nZterm, nTterm). equiv to
-    BuildSums2D(). The result is returned as a double ****
+    BuildSums2D(). The result is returned as a psF64 ****
 *****************************************************************************/
 static psF64 ****BuildSums4D(
@@ -268 +234 @@
     if (sums == NULL) {
         sums = (psF64 ****) psAlloc (nXsum*sizeof(psF64));
-        for (int i = 0; i < nXsum; i++) {
+        for (psS32 i = 0; i < nXsum; i++) {
             sums[i] = (psF64 ***) psAlloc (nYsum*sizeof(psF64));
-            for (int j = 0; j < nYsum; j++) {
+            for (psS32 j = 0; j < nYsum; j++) {
                 sums[i][j] = (psF64 **) psAlloc (nZsum*sizeof(psF64));
-                for (int k = 0; k < nZsum; k++) {
+                for (psS32 k = 0; k < nZsum; k++) {
                     sums[i][j][k] = (psF64 *) psAlloc (nTsum*sizeof(psF64));
                 }
@@ -281 +247 @@
 
     tSum = 1.0;
-    for (int m = 0; m < nTsum; m++) {
+    for (psS32 m = 0; m < nTsum; m++) {
         zSum = tSum;
-        for (int k = 0; k < nZsum; k++) {
+        for (psS32 k = 0; k < nZsum; k++) {
             ySum = zSum;
-            for (int j = 0; j < nYsum; j++) {
+            for (psS32 j = 0; j < nYsum; j++) {
                 xSum = ySum;
-                for (int i = 0; i < nXsum; i++) {
+                for (psS32 i = 0; i < nXsum; i++) {
                     sums[i][j][k][m] = xSum;
                     xSum *= x;
@@ -299 +265 @@
     return (sums);
 }
+
+/******************************************************************************
+ ******************************************************************************
+ Analytical 1-D fitting routines.
+ ******************************************************************************
+ *****************************************************************************/
+
 
 /******************************************************************************
@@ -356 +329 @@
         NUM_DATA = x->n;
     }
-    // XX: psTrace    printf("vectorFitPolynomial1DChebySlow(): NUM_DATA is %d\n", NUM_DATA);
-
-    psPolynomial1D **chebPolys = createChebyshevPolys(NUM_POLY);
+    // psTrace    printf("vectorFitPolynomial1DChebySlow(): NUM_DATA is %d\n", NUM_DATA);
+
+    psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(NUM_POLY);
     if (0) {
         for (psS32 j = 0; j < NUM_POLY; j++) {
@@ -416 +389 @@
     }
 
-    // GaussJordan version
     if (0) {
-        // does the solution in place
-        // XXX: Check error codes!
-        psGaussJordan (A, B);
-
-        // the first nTerm entries in B correspond directly to the desired
-        // polynomial coefficients.  this is only true for the 1D case
-        for (psS32 k = 0; k < NUM_POLY; k++) {
-            myPoly->coeff[k] = B->data.F64[k];
+        // GaussJordan version
+        if (false == psGaussJordan(A, B)) {
+            psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
+            psFree(myPoly);
+            myPoly = NULL;
+        } else {
+            // the first nTerm entries in B correspond directly to the desired
+            // polynomial coefficients.  this is only true for the 1D case
+            for (psS32 k = 0; k < NUM_POLY; k++) {
+                myPoly->coeff[k] = B->data.F64[k];
+            }
         }
     } else {
@@ -433 +408 @@
         psVector* coeffs = NULL;
 
+        // XXX: Check return codes.
         ALUD = psImageAlloc(NUM_POLY, NUM_POLY, PS_TYPE_F64);
         ALUD = psMatrixLUD(ALUD, &outPerm, A);
@@ -468 +444 @@
     the [-1:1] interval.
  
-XXX: mask, maskValue, yErr are currently ignored.
+    NOTE: mask, maskValue, yErr are ignored with this function.
 *****************************************************************************/
 static psPolynomial1D *vectorFitPolynomial1DChebyFast(
@@ -501 +477 @@
     tmpScalar.type.type = PS_TYPE_F64;
 
-    // XXX: These assignments appear too simple to warrant code and
+    // These assignments appear too simple to warrant code and
     // variable declarations.  I retain them here to maintain coherence
     // with the NR code.
@@ -538 +514 @@
             fScalar = p_psVectorInterpolate((psVector *) x, (psVector *) y,
                                             3, &tmpScalar);
+            if (fScalar == NULL) {
+                psError(PS_ERR_UNKNOWN, false, "Could not perform vector interpolation.  Returning NULL.\n");
+                psFree(myPoly)
+                return(NULL);
+            }
+
             f->data.F64[i] = fScalar->data.F64;
             psFree(fScalar);
         }
 
-        psTrace(".psLib.dataManip.vectorFitPolynomial1DCheby", 6,
-                "(x, X, y, f(X)) is (%f, %f, %f, %f)\n",
+        psTrace(__func__, 6, "(x, X, y, f(X)) is (%f, %f, %f, %f)\n",
                 x->data.F64[i], X, y->data.F64[i], f->data.F64[i]);
     }
@@ -551 +532 @@
 
     fac = 2.0/((psF32) n);
-    // XXX: is this loop bound correct?
     for (j=0;j<myPoly->nX+1;j++) {
         sum = 0.0;
@@ -573 +553 @@
 PS_TYPE_F64.
  *****************************************************************************/
-psPolynomial1D* VectorFitPolynomial1DOrd(
+static psPolynomial1D* VectorFitPolynomial1DOrd(
     psPolynomial1D* myPoly,
     const psVector *mask,
@@ -581 +561 @@
     const psVector *x)
 {
-    psTrace(__func__, 4, "---- VectorFitPolynomial1DOrd() begin ----\n");
-    // XXX: these ASSERTS are redundant.
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
     PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
@@ -600 +579 @@
     }
 
-    psImage*     A = NULL;
-    psVector*    B = NULL;
+    psImage* A = NULL;
+    psVector* B = NULL;
     psVector* xSums = NULL;
     psS32 nTerm;
@@ -607 +586 @@
     psF64 wt;
 
-    if (psTraceGetLevel (__func__) >= 6) {
+    if (psTraceGetLevel(__func__) >= 6) {
         psTrace(__func__, 6, "VectorFitPolynomial1D()\n");
         for (psS32 i = 0; i < f->n; i++) {
@@ -627 +606 @@
     nOrder = nTerm - 1;
 
-    A     = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
-    B     = psVectorAlloc(nTerm, PS_TYPE_F64);
-
-    //
+    A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
+    B = psVectorAlloc(nTerm, PS_TYPE_F64);
     // Initialize data structures.
-    // XXX: Use psLib function.
-    //
-    PS_VECTOR_SET_F64(B, 0.0);
-    PS_IMAGE_SET_F64(A, 0.0);
+    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(__func__, 4, "---- %s() End ----\n", __func__);
+        return(NULL);
+    }
 
     // xSums look like: 1, x, x^2, ... x^(2n+1)
@@ -641 +622 @@
     // XXX EAM : use temp pointers eg vB = B->data.F64 to save redirects
     // XXX EAM : this function is only valid for data vectors of F64
-    for (int k = 0; k < f->n; k++) {
+    for (psS32 k = 0; k < f->n; k++) {
         if ((mask != NULL) && (mask->data.U8[k] && maskValue)) {
             continue;
@@ -657 +638 @@
             wt = (fErr->data.F64[k] == 0) ? 0.0 : 1.0 / PS_SQR(fErr->data.F64[k]);
         }
-        for (int i = 0; i < nTerm; i++) {
+        for (psS32 i = 0; i < nTerm; i++) {
             B->data.F64[i] += f->data.F64[k] * xSums->data.F64[i] * wt;
         }
@@ -663 +644 @@
         // we could skip half of the array and assign at the end
         // we must handle masked orders
-        for (int i = 0; i < nTerm; i++) {
-            for (int j = 0; j < nTerm; j++) {
+        for (psS32 i = 0; i < nTerm; i++) {
+            for (psS32 j = 0; j < nTerm; j++) {
                 A->data.F64[i][j] += xSums->data.F64[i + j] * wt;
             }
@@ -670 +651 @@
     }
 
-    // GaussJordan version
     if (0) {
-        // does the solution in place
-        // XXX: Check error codes!
-        psGaussJordan (A, B);
-
-        // the first nTerm entries in B correspond directly to the desired
-        // polynomial coefficients.  this is only true for the 1D case
-        for (int k = 0; k < nTerm; k++) {
-            myPoly->coeff[k] = B->data.F64[k];
-        }
+        // GaussJordan version
+        if (false == psGaussJordan(A, B)) {
+            psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
+            psFree(myPoly);
+            myPoly = NULL;
+        } else {
+            // the first nTerm entries in B correspond directly to the desired
+            // polynomial coefficients.  this is only true for the 1D case
+            for (psS32 k = 0; k < nTerm; k++) {
+                myPoly->coeff[k] = B->data.F64[k];
+            }
+        }
+
     } else {
         // LUD version of the fit
@@ -687 +671 @@
         psVector* coeffs = NULL;
 
+        // XXX: Check return codes.
         ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
         ALUD = psMatrixLUD(ALUD, &outPerm, A);
         coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-        for (int k = 0; k < nTerm; k++) {
+        for (psS32 k = 0; k < nTerm; k++) {
             myPoly->coeff[k] = coeffs->data.F64[k];
         }
@@ -703 +688 @@
     psFree(xSums);
 
-    psTrace(__func__, 4, "---- VectorFitPolynomial1DOrd() End ----\n");
+    psTrace(__func__, 4, "---- %s() End ----\n", __func__);
     return (myPoly);
 }
@@ -713 +698 @@
 Types F32 and F64 are supported, however, type F32 is done via vector
 conversion only.
- 
-XXX: Put the asserts at the beginning.  Why is this not failing tests?
  *****************************************************************************/
 psPolynomial1D *psVectorFitPolynomial1D(
@@ -748 +731 @@
 
     if (f->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(f64, f);
+        f64 = psVectorCopy(NULL, f, PS_TYPE_F64);
     } else {
         f64 = (psVector *) f;
@@ -755 +738 @@
     if (x != NULL) {
         if (x->type.type != PS_TYPE_F64) {
-            PS_VECTOR_GEN_F64_FROM_F32(x64, x);
+            x64 = psVectorCopy(NULL, x, PS_TYPE_F64);
         } else {
             x64 = (psVector *) x;
@@ -763 +746 @@
     if (fErr != NULL) {
         if (fErr->type.type != PS_TYPE_F64) {
-            PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr);
+            fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64);
         } else {
             fErr64 = (psVector *) fErr;
@@ -784 +767 @@
         if (x == NULL) {
             // If x==NULL, create an x64 vector with x values set to (-1:1).
-            PS_VECTOR_GEN_CHEBY_INDEX(x64, f64->n);
+            PS_VECTOR_GEN_CHEBY_INDEX(x64, f64->n, PS_TYPE_F64);
         }
 
@@ -850 +833 @@
         } else if (poly->type == PS_POLYNOMIAL_CHEB) {
             if (f->type.type == PS_TYPE_F32) {
-                PS_VECTOR_GEN_CHEBY_INDEX_F32(x, f->n);
+                PS_VECTOR_GEN_CHEBY_INDEX(x, f->n, PS_TYPE_F32);
             } else if (f->type.type == PS_TYPE_F64) {
-                PS_VECTOR_GEN_CHEBY_INDEX_F64(x, f->n);
+                PS_VECTOR_GEN_CHEBY_INDEX(x, f->n, PS_TYPE_F64);
             }
         } else {
-            printf("XXX: error, bad poly type.\n");
+            psError(PS_ERR_UNKNOWN, false, "Error, bad poly type.\n");
+            return(NULL);
         }
     }
@@ -885 +869 @@
 
     //
-    for (int N = 0; N < stats->clipIter; N++) {
+    for (psS32 N = 0; N < stats->clipIter; N++) {
         psTrace(__func__, 6, "Loop iteration %d.  Calling psVectorFitPolynomial1D()\n");
-        int Nkeep = 0;
+        psS32 Nkeep = 0;
         if (psTraceGetLevel(__func__) >= 6) {
             if (mask != NULL) {
@@ -897 +881 @@
         poly = psVectorFitPolynomial1D(poly, mask, maskValue, f, fErr, x);
         if (poly == NULL) {
-            printf("XXX: psVectorFitPolynomial1D() returned NULL: Generate error, free data.\n");
+            psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial.  Returning NULL.\n");
+            if (xIn == NULL) {
+                psFree(x);
+            }
             return(NULL);
         }
 
-        // XXX: Check error codes
         fit = psPolynomial1DEvalVector(poly, x);
+        if (fit == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector().  Returning NULL.\n");
+            psFree(resid)
+            return(NULL);
+        }
         for (psS32 i = 0 ; i < f->n ; i++) {
             if (f->type.type == PS_TYPE_F64) {
@@ -930 +921 @@
         // we are masking out any point which is out of range
         // recovery is not allowed with this scheme
-        for (int i = 0; i < resid->n; i++) {
+        for (psS32 i = 0; i < resid->n; i++) {
             if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
                 continue;
@@ -954 +945 @@
 
         //
-        // XXX: We should probably exit this loop if no new elements were masked
+        // We should probably exit this loop if no new elements were masked
         // since the polynomial fit won't change.
         //
@@ -984 +975 @@
 pairs.  All non-NULL vectors must be of type PS_TYPE_F64.
  
-// XXX: What about the maskValue?
  *****************************************************************************/
-psPolynomial2D* VectorFitPolynomial2DOrd(
+static psPolynomial2D* VectorFitPolynomial2DOrd(
     psPolynomial2D* myPoly,
     const psVector* mask,
@@ -995 +985 @@
     const psVector *y)
 {
-    psTrace(__func__, 4, "---- VectorFitPolynomial2DOrd() begin ----\n");
-    // These ASSERTS are redundant.
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
     PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
@@ -1025 +1014 @@
     psS32 nTerm;
 
-    // XXX:Watch for changes to the psPolys: nTerm != nOrder.
     psS32 nXterm = 1 + myPoly->nX;
     psS32 nYterm = 1 + myPoly->nY;
@@ -1032 +1020 @@
     A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
     B = psVectorAlloc(nTerm, PS_TYPE_F64);
-
-    //
     // Initialize data structures.
-    // XXX: Use psLib function.
-    //
-    PS_VECTOR_SET_F64(B, 0.0);
-    PS_IMAGE_SET_F64(A, 0.0);
+    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(__func__, 4, "---- %s() End ----\n", __func__);
+        return(NULL);
+    }
 
     // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)
 
     // Build the B and A data structs.
-    for (int k  = 0; k < x->n; k++) {
+    for (psS32 k  = 0; k < x->n; k++) {
         if ((mask != NULL) && (mask->data.U8[k] & maskValue)) {
             continue;
@@ -1059 +1049 @@
         // we could skip half of the array and assign at the end
         // we must handle masked orders
-        for (int n = 0; n < nXterm; n++) {
-            for (int m = 0; m < nYterm; m++) {
+        for (psS32 n = 0; n < nXterm; n++) {
+            for (psS32 m = 0; m < nYterm; m++) {
                 B->data.F64[n+m*nXterm] += f->data.F64[k] * Sums->data.F64[n][m] * wt;
             }
         }
 
-        for (int i = 0; i < nXterm; i++) {
-            for (int j = 0; j < nYterm; j++) {
-                for (int n = 0; n < nXterm; n++) {
-                    for (int m = 0; m < nYterm; m++) {
+        for (psS32 i = 0; i < nXterm; i++) {
+            for (psS32 j = 0; j < nYterm; j++) {
+                for (psS32 n = 0; n < nXterm; n++) {
+                    for (psS32 m = 0; m < nYterm; m++) {
                         A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt;
                     }
@@ -1076 +1066 @@
     }
 
-    // does the solution in place
-    // XXX: Check error codes!
-    psGaussJordan (A, B);
-
-    // select the appropriate solution entries
-    for (int n = 0; n < nXterm; n++) {
-        for (int m = 0; m < nYterm; m++) {
-            myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
+    if (false == psGaussJordan(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 (psS32 n = 0; n < nXterm; n++) {
+            for (psS32 m = 0; m < nYterm; m++) {
+                myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
+            }
         }
     }
@@ -1091 +1083 @@
     psFree(Sums);
 
-    psTrace(__func__, 4, "---- VectorFitPolynomial2DOrd() end ----\n");
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
     return (myPoly);
 }
@@ -1139 +1131 @@
     //
     if (f->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(f64, f);
+        f64 = psVectorCopy(NULL, f, PS_TYPE_F64);
     } else {
         f64 = (psVector *) f;
@@ -1145 +1137 @@
 
     if (x->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(x64, x);
+        x64 = psVectorCopy(NULL, x, PS_TYPE_F64);
     } else {
         x64 = (psVector *) x;
@@ -1151 +1143 @@
 
     if (y->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(y64, y);
+        y64 = psVectorCopy(NULL, y, PS_TYPE_F64);
     } else {
         y64 = (psVector *) y;
@@ -1161 +1153 @@
     if (fErr != NULL) {
         if (fErr->type.type != PS_TYPE_F64) {
-            PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr);
+            fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64);
         } else {
             fErr64 = (psVector *) fErr;
@@ -1272 +1264 @@
 
     // clipping range defined by min and max and/or clipSigma
-    float minClipSigma;
-    float maxClipSigma;
+    psF32 minClipSigma;
+    psF32 maxClipSigma;
     if (isfinite(stats->max)) {
         maxClipSigma = fabs(stats->max);
@@ -1294 +1286 @@
     psTrace(__func__, 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
 
-    for (int N = 0; N < stats->clipIter; N++) {
+    for (psS32 N = 0; N < stats->clipIter; N++) {
         psTrace(__func__, 6, "Loop iteration %d.  Calling psVectorFitPolynomial1D()\n");
-        int Nkeep = 0;
+        psS32 Nkeep = 0;
         if (psTraceGetLevel(__func__) >= 6) {
             if (mask != NULL) {
@@ -1306 +1298 @@
 
         poly = psVectorFitPolynomial2D(poly, mask, maskValue, f, fErr, x, y);
-        // XXX: Check error codes
+        if (poly == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data.  Returning NULL.\n");
+            psFree(resid)
+            return(NULL);
+        }
+
         psVector *fit = psPolynomial2DEvalVector(poly, x, y);
+        if (fit == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector().  Returning NULL.\n");
+            psFree(resid)
+            return(NULL);
+        }
 
         for (psS32 i = 0 ; i < f->n ; i++) {
@@ -1328 +1330 @@
         }
 
-        // XXX: Check error codes
         stats  = psVectorStats(stats, resid, NULL, mask, maskValue);
+        if (stats == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector.  Returning NULL.\n");
+            psFree(resid)
+            psFree(fit)
+            return(NULL);
+        }
         psTrace(__func__, 6, "Median is %f\n", stats->sampleMedian);
         psTrace(__func__, 6, "Stdev is %f\n", stats->sampleStdev);
-        float minClipValue = -minClipSigma*stats->sampleStdev;
-        float maxClipValue = +maxClipSigma*stats->sampleStdev;
+        psF32 minClipValue = -minClipSigma*stats->sampleStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
 
         // set mask if pts are not valid
         // we are masking out any point which is out of range
         // recovery is not allowed with this scheme
-        for (int i = 0; i < resid->n; i++) {
+        for (psS32 i = 0; i < resid->n; i++) {
             if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
                 continue;
@@ -1388 +1395 @@
 y, z)-(f) pairs.  All non-NULL vectors must be of type PS_TYPE_F64.
  
-XXX: This routine needs to be written.  Currently, this is simply a shell.  We
-can assume that all vectors have been converted to F64, that (f, x, y, z) are
-non-null and F64.  fErr may be NULL, but will be F64 is not.
  *****************************************************************************/
-psPolynomial3D* VectorFitPolynomial3DOrd(
+static psPolynomial3D* VectorFitPolynomial3DOrd(
     psPolynomial3D* myPoly,
     const psVector* mask,
@@ -1402 +1406 @@
     const psVector *z)
 {
-    psTrace(__func__, 4, "---- VectorFitPolynomial3DOrd() begin ----\n");
-    // These ASSERTS are redundant.
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
     PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
@@ -1435 +1438 @@
     psS32 nTerm;
 
-    // XXX:Watch for changes to the psPolys: nTerm != nOrder.
     psS32 nXterm = 1 + myPoly->nX;
     psS32 nYterm = 1 + myPoly->nY;
@@ -1443 +1445 @@
     A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
     B = psVectorAlloc(nTerm, PS_TYPE_F64);
-
     // Initialize data structures.
-    psVectorInit (B, 0.0);
-    psImageInit (A, 0.0);
+    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(__func__, 4, "---- %s() End ----\n", __func__);
+        return(NULL);
+    }
 
     // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)*y, ...
 
     // Build the B and A data structs.
-    for (int k  = 0; k < x->n; k++) {
+    for (psS32 k  = 0; k < x->n; k++) {
         if ((mask != NULL) && (mask->data.U8[k] & maskValue)) {
             continue;
@@ -1467 +1474 @@
         // we could skip half of the array and assign at the end
         // we must handle masked orders
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
                     if (myPoly->mask[ix][iy][iz]) {
                         continue;
                     } else {
-                        int nx = ix + iy*nXterm + iz*nXterm*nYterm;
+                        psS32 nx = ix + iy*nXterm + iz*nXterm*nYterm;
                         B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz] * wt;
                     }
@@ -1480 +1487 @@
         }
 
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
                     if (myPoly->mask[ix][iy][iz])
                         continue;
-                    int nx = ix+iy*nXterm+iz*nXterm*nYterm;
-                    for (int jx = 0; jx < nXterm; jx++) {
-                        for (int jy = 0; jy < nYterm; jy++) {
-                            for (int jz = 0; jz < nZterm; jz++) {
+                    psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
+                    for (psS32 jx = 0; jx < nXterm; jx++) {
+                        for (psS32 jy = 0; jy < nYterm; jy++) {
+                            for (psS32 jz = 0; jz < nZterm; jz++) {
                                 if (myPoly->mask[jx][jy][jz])
                                     continue;
-                                int ny = jx+jy*nXterm+jz*nXterm*nYterm;
+                                psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm;
                                 A->data.F64[nx][ny] += Sums[ix+jx][iy+jy][iz+jz] * wt;
                             }
@@ -1501 +1508 @@
     }
 
-    for (int ix = 0; ix < nXterm; ix++) {
-        for (int iy = 0; iy < nYterm; iy++) {
-            for (int iz = 0; iz < nZterm; iz++) {
+    for (psS32 ix = 0; ix < nXterm; ix++) {
+        for (psS32 iy = 0; iy < nYterm; iy++) {
+            for (psS32 iz = 0; iz < nZterm; iz++) {
                 if (!myPoly->mask[ix][iy][iz])
                     continue;
-                int nx = ix+iy*nXterm+iz*nXterm*nYterm;
+                psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
                 B->data.F64[nx] = 0;
-                for (int jx = 0; jx < nXterm; jx++) {
-                    for (int jy = 0; jy < nYterm; jy++) {
-                        for (int jz = 0; jz < nZterm; jz++) {
-                            int ny = jx+jy*nXterm+jz*nXterm*nYterm;
+                for (psS32 jx = 0; jx < nXterm; jx++) {
+                    for (psS32 jy = 0; jy < nYterm; jy++) {
+                        for (psS32 jz = 0; jz < nZterm; jz++) {
+                            psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm;
                             A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
                         }
@@ -1523 +1530 @@
         // does the solution in place
         // The matrices were overflowing, so I switched to LUD.
-        if (false == psGaussJordan (A, B)) {
+        if (false == psGaussJordan(A, B)) {
             psFree(A);
             psFree(B);
 
-            for (int ix = 0; ix < 2*nXterm; ix++) {
-                for (int iy = 0; iy < 2*nYterm; iy++) {
+            for (psS32 ix = 0; ix < 2*nXterm; ix++) {
+                for (psS32 iy = 0; iy < 2*nYterm; iy++) {
                     psFree(Sums[ix][iy]);
                 }
@@ -1538 +1545 @@
         }
         // select the appropriate solution entries
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
-                    int nx = ix+iy*nXterm+iz*nXterm*nYterm;
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
+                    psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
                     myPoly->coeff[ix][iy][iz] = B->data.F64[nx];
                     myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
@@ -1553 +1560 @@
         psVector* coeffs = NULL;
 
+        // XXX: Check return codes.
         ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
         ALUD = psMatrixLUD(ALUD, &outPerm, A);
@@ -1558 +1566 @@
 
         // select the appropriate solution entries
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
-                    int nx = ix+iy*nXterm+iz*nXterm*nYterm;
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
+                    psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
                     myPoly->coeff[ix][iy][iz] = coeffs->data.F64[nx];
                     myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
@@ -1575 +1583 @@
     psFree(B);
 
-    for (int ix = 0; ix < 2*nXterm; ix++) {
-        for (int iy = 0; iy < 2*nYterm; iy++) {
+    for (psS32 ix = 0; ix < 2*nXterm; ix++) {
+        for (psS32 iy = 0; iy < 2*nYterm; iy++) {
             psFree(Sums[ix][iy]);
         }
@@ -1583 +1591 @@
     psFree(Sums);
 
-    psTrace(__func__, 4, "---- VectorFitPolynomial3DOrd() end ----\n");
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
     return (myPoly);
 }
@@ -1621 +1629 @@
     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);
     if (fErr != NULL) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
@@ -1637 +1642 @@
     //
     if (f->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(f64, f);
+        f64 = psVectorCopy(NULL, f, PS_TYPE_F64);
     } else {
         f64 = (psVector *) f;
     }
     if (x->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(x64, x);
+        x64 = psVectorCopy(NULL, x, PS_TYPE_F64);
     } else {
         x64 = (psVector *) x;
     }
     if (y->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(y64, y);
+        y64 = psVectorCopy(NULL, y, PS_TYPE_F64);
     } else {
         y64 = (psVector *) y;
@@ -1653 +1658 @@
 
     if (z->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(z64, z);
+        z64 = psVectorCopy(NULL, z, PS_TYPE_F64);
     } else {
         z64 = (psVector *) z;
@@ -1660 +1665 @@
     if (fErr != NULL) {
         if (fErr->type.type != PS_TYPE_F64) {
-            PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr);
+            fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64);
         } else {
             fErr64 = (psVector *) fErr;
@@ -1788 +1793 @@
 
     // clipping range defined by min and max and/or clipSigma
-    float minClipSigma;
-    float maxClipSigma;
+    psF32 minClipSigma;
+    psF32 maxClipSigma;
     if (isfinite(stats->max)) {
         maxClipSigma = fabs(stats->max);
@@ -1811 +1816 @@
     psTrace(__func__, 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
 
-    for (int N = 0; N < stats->clipIter; N++) {
+    for (psS32 N = 0; N < stats->clipIter; N++) {
         psTrace(__func__, 6, "Loop iteration %d.  Calling psVectorFitPolynomial1D()\n");
-        int Nkeep = 0;
+        psS32 Nkeep = 0;
         if (psTraceGetLevel(__func__) >= 6) {
             if (mask != NULL) {
@@ -1823 +1828 @@
 
         poly = psVectorFitPolynomial3D(poly, mask, maskValue, f, fErr, x, y, z);
-        // XXX: Check error codes
+        if (poly == NULL) {
+            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);
+        if (fit == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector().  Returning NULL.\n");
+            psFree(resid)
+            return(NULL);
+        }
         for (psS32 i = 0 ; i < f->n ; i++) {
             if (f->type.type == PS_TYPE_F64) {
@@ -1844 +1859 @@
         }
 
-        // XXX: Check error codes
         stats  = psVectorStats(stats, resid, NULL, mask, maskValue);
+        if (stats == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector.  Returning NULL.\n");
+            psFree(resid)
+            psFree(fit)
+            return(NULL);
+        }
+
         psTrace(__func__, 6, "Median is %f\n", stats->sampleMedian);
         psTrace(__func__, 6, "Stdev is %f\n", stats->sampleStdev);
-        float minClipValue = -minClipSigma*stats->sampleStdev;
-        float maxClipValue = +maxClipSigma*stats->sampleStdev;
+        psF32 minClipValue = -minClipSigma*stats->sampleStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
 
         // set mask if pts are not valid
         // we are masking out any point which is out of range
         // recovery is not allowed with this scheme
-        for (int i = 0; i < resid->n; i++) {
+        for (psS32 i = 0; i < resid->n; i++) {
             if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
                 continue;
@@ -1902 +1923 @@
 y, z, t)-(f) pairs.  All non-NULL vectors must be of type PS_TYPE_F64.
  
-XXX: This routine needs to be written.  Currently, this is simply a shell.  We
-can assume that all vectors have been converted to F64, that (f, x, y, z) are
-non-null and F64.  fErr may be NULL, but will be F64 is not.
  *****************************************************************************/
-psPolynomial4D* VectorFitPolynomial4DOrd(
+static psPolynomial4D* VectorFitPolynomial4DOrd(
     psPolynomial4D* myPoly,
     const psVector* mask,
@@ -1918 +1936 @@
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    // These ASSERTS are redundant.
     PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
     PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
@@ -1954 +1971 @@
     psS32 nTerm;
 
-    // XXX:Watch for changes to the psPolys: nTerm != nOrder.
     psS32 nXterm = 1 + myPoly->nX;
     psS32 nYterm = 1 + myPoly->nY;
@@ -1963 +1979 @@
     A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
     B = psVectorAlloc(nTerm, PS_TYPE_F64);
-
     // Initialize data structures.
-    psVectorInit (B, 0.0);
-    psImageInit (A, 0.0);
+    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(__func__, 4, "---- %s() End ----\n", __func__);
+        return(NULL);
+    }
 
     // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)*y, ...
 
     // Build the B and A data structs.
-    for (int k  = 0; k < x->n; k++) {
+    for (psS32 k  = 0; k < x->n; k++) {
         if ((mask != NULL) && (mask->data.U8[k] & maskValue)) {
             continue;
@@ -1987 +2008 @@
         // we could skip half of the array and assign at the end
         // we must handle masked orders
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
-                    for (int it = 0; it < nTterm; it++) {
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
+                    for (psS32 it = 0; it < nTterm; it++) {
                         if (myPoly->mask[ix][iy][iz][it])
                             continue;
-                        int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+                        psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
                         B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz][it] * wt;
                     }
@@ -2000 +2021 @@
         }
 
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
-                    for (int it = 0; it < nTterm; it++) {
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
+                    for (psS32 it = 0; it < nTterm; it++) {
                         if (myPoly->mask[ix][iy][iz][it])
                             continue;
-                        int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
-                        for (int jx = 0; jx < nXterm; jx++) {
-                            for (int jy = 0; jy < nYterm; jy++) {
-                                for (int jz = 0; jz < nZterm; jz++) {
-                                    for (int jt = 0; jt < nTterm; jt++) {
+                        psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+                        for (psS32 jx = 0; jx < nXterm; jx++) {
+                            for (psS32 jy = 0; jy < nYterm; jy++) {
+                                for (psS32 jz = 0; jz < nZterm; jz++) {
+                                    for (psS32 jt = 0; jt < nTterm; jt++) {
                                         if (myPoly->mask[jx][jy][jz][jt])
                                             continue;
-                                        int ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
+                                        psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
                                         A->data.F64[nx][ny]+= Sums[ix+jx][iy+jy][iz+jz][it+jt] * wt;
                                     }
@@ -2025 +2046 @@
     }
 
-    for (int ix = 0; ix < nXterm; ix++) {
-        for (int iy = 0; iy < nYterm; iy++) {
-            for (int iz = 0; iz < nZterm; iz++) {
-                for (int it = 0; it < nTterm; it++) {
+    for (psS32 ix = 0; ix < nXterm; ix++) {
+        for (psS32 iy = 0; iy < nYterm; iy++) {
+            for (psS32 iz = 0; iz < nZterm; iz++) {
+                for (psS32 it = 0; it < nTterm; it++) {
                     if (!myPoly->mask[ix][iy][iz][it])
                         continue;
-                    int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+                    psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
                     B->data.F64[nx] = 0;
-                    for (int jx = 0; jx < nXterm; jx++) {
-                        for (int jy = 0; jy < nYterm; jy++) {
-                            for (int jz = 0; jz < nZterm; jz++) {
-                                for (int jt = 0; jt < nTterm; jt++) {
-                                    int ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
+                    for (psS32 jx = 0; jx < nXterm; jx++) {
+                        for (psS32 jy = 0; jy < nYterm; jy++) {
+                            for (psS32 jz = 0; jz < nZterm; jz++) {
+                                for (psS32 jt = 0; jt < nTterm; jt++) {
+                                    psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
                                     A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
                                 }
@@ -2051 +2072 @@
     if (0) {
         // does the solution in place
-        // XXX: The GaussJordan version was overflowing, so I'm using LUD.
+        // The GaussJordan version was overflowing, so I'm using LUD.
         if (false == psGaussJordan(A, B)) {
             psFree(A);
             psFree(B);
-            for (int ix = 0; ix < 2*nXterm; ix++) {
-                for (int iy = 0; iy < 2*nYterm; iy++) {
-                    for (int iz = 0; iz < 2*nZterm; iz++) {
+            for (psS32 ix = 0; ix < 2*nXterm; ix++) {
+                for (psS32 iy = 0; iy < 2*nYterm; iy++) {
+                    for (psS32 iz = 0; iz < 2*nZterm; iz++) {
                         psFree(Sums[ix][iy][iz]);
                     }
@@ -2070 +2091 @@
 
         // select the appropriate solution entries
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
-                    for (int it = 0; it < nTterm; it++) {
-                        int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
+                    for (psS32 it = 0; it < nTterm; it++) {
+                        psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
                         myPoly->coeff[ix][iy][iz][it] = B->data.F64[nx];
                         myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
@@ -2087 +2108 @@
         psVector* coeffs = NULL;
 
+        // XXX: Check return codes.
         ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
         ALUD = psMatrixLUD(ALUD, &outPerm, A);
@@ -2092 +2114 @@
 
         // select the appropriate solution entries
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
-                    for (int it = 0; it < nTterm; it++) {
-                        int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
+                    for (psS32 it = 0; it < nTterm; it++) {
+                        psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
                         myPoly->coeff[ix][iy][iz][it] = coeffs->data.F64[nx];
                         myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
@@ -2113 +2135 @@
     psFree(B);
 
-    for (int ix = 0; ix < 2*nXterm; ix++) {
-        for (int iy = 0; iy < 2*nYterm; iy++) {
-            for (int iz = 0; iz < 2*nZterm; iz++) {
+    for (psS32 ix = 0; ix < 2*nXterm; ix++) {
+        for (psS32 iy = 0; iy < 2*nYterm; iy++) {
+            for (psS32 iz = 0; iz < 2*nZterm; iz++) {
                 psFree(Sums[ix][iy][iz]);
             }
@@ -2180 +2202 @@
     //
     if (f->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(f64, f);
+        f64 = psVectorCopy(NULL, f, PS_TYPE_F64);
     } else {
         f64 = (psVector *) f;
     }
     if (x->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(x64, x);
+        x64 = psVectorCopy(NULL, x, PS_TYPE_F64);
     } else {
         x64 = (psVector *) x;
     }
     if (y->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(y64, y);
+        y64 = psVectorCopy(NULL, y, PS_TYPE_F64);
     } else {
         y64 = (psVector *) y;
     }
     if (z->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(z64, z);
+        z64 = psVectorCopy(NULL, z, PS_TYPE_F64);
     } else {
         z64 = (psVector *) z;
     }
     if (t->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(t64, t);
+        t64 = psVectorCopy(NULL, t, PS_TYPE_F64);
     } else {
         t64 = (psVector *) t;
@@ -2209 +2231 @@
     if (fErr != NULL) {
         if (fErr->type.type != PS_TYPE_F64) {
-            PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr);
+            fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64);
         } else {
             fErr64 = (psVector *) fErr;
@@ -2310 +2332 @@
 
 
-// XXX: Add F64 support.
 psPolynomial4D *psVectorClipFitPolynomial4D(
     psPolynomial4D *poly,
@@ -2356 +2377 @@
 
     // clipping range defined by min and max and/or clipSigma
-    float minClipSigma;
-    float maxClipSigma;
+    psF32 minClipSigma;
+    psF32 maxClipSigma;
     if (isfinite(stats->max)) {
         maxClipSigma = fabs(stats->max);
@@ -2379 +2400 @@
     psTrace(__func__, 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
 
-    for (int N = 0; N < stats->clipIter; N++) {
+    for (psS32 N = 0; N < stats->clipIter; N++) {
         psTrace(__func__, 6, "Loop iteration %d.  Calling psVectorFitPolynomial1D()\n");
-        int Nkeep = 0;
+        psS32 Nkeep = 0;
         if (psTraceGetLevel(__func__) >= 6) {
             if (mask != NULL) {
@@ -2391 +2412 @@
 
         poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t);
-        // XXX: Check error codes
+        if (poly == NULL) {
+            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);
-        // XXX: This is coded differently for 1D, 2D.  Consolidate.
+        if (fit == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector().  Returning NULL.\n");
+            psFree(resid)
+            return(NULL);
+        }
         for (psS32 i = 0 ; i < f->n ; i++) {
             if (f->type.type == PS_TYPE_F64) {
@@ -2401 +2432 @@
             }
         }
-        // XXX: Check error codes
 
         if (psTraceGetLevel(__func__) >= 6) {
@@ -2414 +2444 @@
         }
 
-        // XXX: Check error codes
         stats  = psVectorStats (stats, resid, NULL, mask, maskValue);
+        if (stats == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector.  Returning NULL.\n");
+            psFree(resid)
+            psFree(fit)
+            return(NULL);
+        }
         psTrace(__func__, 6, "Median is %f\n", stats->sampleMedian);
         psTrace(__func__, 6, "Stdev is %f\n", stats->sampleStdev);
-        float minClipValue = -minClipSigma*stats->sampleStdev;
-        float maxClipValue = +maxClipSigma*stats->sampleStdev;
+        psF32 minClipValue = -minClipSigma*stats->sampleStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
 
         // set mask if pts are not valid
         // we are masking out any point which is out of range
         // recovery is not allowed with this scheme
-        for (int i = 0; i < resid->n; i++) {
+        for (psS32 i = 0; i < resid->n; i++) {
             if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
                 continue;