Changeset 6186

Timestamp:

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

Author:

gusciora

Message:

Misc code cleaning

Location:

trunk/psLib/src

Files:

• astro/psCoord.c (modified) (2 diffs)
• imageops/psImageStats.c (modified) (3 diffs)
• math/psConstants.h (modified) (4 diffs)
• math/psMinimizePolyFit.c (modified) (96 diffs)
• math/psMinimizePowell.c (modified) (4 diffs)
• math/psPolynomial.c (modified) (17 diffs)
• math/psPolynomial.h (modified) (6 diffs)
• math/psSpline.c (modified) (4 diffs)
• math/psStats.c (modified) (5 diffs)

trunk/psLib/src/astro/psCoord.c

@@ -10 +10 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.99 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-12-24 00:33:13 $
+*  @version $Revision: 1.100 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-23 22:25:31 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -925 +925 @@
         return(p_psPlaneTransformLinearInvert((psPlaneTransform *) in));
     }
-    PS_INT_COMPARE(1, nSamples, NULL);
+    PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(nSamples, 1, NULL);
 
     // Ensure that the input transformation is symmetrical.

trunk/psLib/src/imageops/psImageStats.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.86 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-24 00:33:17 $
+ *  @version $Revision: 1.87 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-23 22:25:31 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -212 +212 @@
 
     return (scalingFactors);
-}
-
-// XXX: Use a static array of Chebyshev polynomials.
-psPolynomial1D **p_psCreateChebyshevPolys(psS32 maxChebyPoly)
-{
-    PS_ASSERT_INT_POSITIVE(maxChebyPoly, NULL);
-    psPolynomial1D **chebPolys = NULL;
-    psS32 i = 0;
-    psS32 j = 0;
-
-    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
-    for (i = 0; i < maxChebyPoly; i++) {
-        // XXX: verify this, poly nOrder/nTerms change.
-        chebPolys[i] = psPolynomial1DAlloc(i, PS_POLYNOMIAL_ORD);
-    }
-
-    // Create the Chebyshev polynomials.
-    // Polynomial i has i-th order.
-    chebPolys[0]->coeff[0] = 1;
-    chebPolys[1]->coeff[1] = 1;
-    for (i = 2; i < maxChebyPoly; i++) {
-        for (j = 0; j < (1 + chebPolys[i - 1]->nX); j++) {
-            chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j];
-        }
-        for (j = 0; j < (1 + chebPolys[i - 2]->nX); j++) {
-            chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
-        }
-    }
-
-    return (chebPolys);
 }
 
@@ -302 +272 @@
     // Determine how many Chebyshev polynomials
     // are needed, then create them.
-    // XXX: recorde or verify the poly order/nterm change
+    // XXX: record or verify the poly order/nterm change
     maxChebyPoly = coeffs->nX;
     if (coeffs->nY > coeffs->nX) {

trunk/psLib/src/math/psConstants.h

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.81 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-11-23 23:54:43 $
+ *  @version $Revision: 1.82 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-23 22:25:31 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -14 +14 @@
  *       called with complex expressions.
  *
- *  XXX: All functions which use the PS_CHECK macros must be scrutinized so
+ *  XXX: All functions which use the PS_ASSERT macros must be scrutinized so
  *  that we ensure that an argument which is expected to be output is
  *  psFree'ed before reurning NULL.
@@ -21 +21 @@
  *  throw a psError if the CONDITION is true.  However, some throw the error
  *  if the CONDITION is false.  This should be consistant.
- *
- *  XXX: rename all these with form PS_ASSERT_CONDITION().
  *
  */
@@ -156 +154 @@
 }
 
-// Produce an error if (NAME1 > NAME2)
-// XXX: Get rid of this, use above macros.
-#define PS_INT_COMPARE(NAME1, NAME2, RVAL) \
-if ((NAME1) > (NAME2)) { \
-    psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
-            "Error: (%s > %s) (%d %d).", \
-            #NAME1, #NAME2, NAME1, NAME2); \
-    return(RVAL); \
-}
-
-// Produce an error if ((NAME1 > NAME2)
-// XXX: Get rid of this, use above macros.
-#define PS_FLOAT_COMPARE(NAME1, NAME2, RVAL) \
-if ((NAME1) > (NAME2)) { \
-    psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
-            "Error: (%s > %s) (%f %f)", \
-            #NAME1, #NAME2, NAME1, NAME2); \
-    return(RVAL); \
-}
-
 #define PS_ASSERT_FLOAT_NON_EQUAL(NAME1, NAME2, RVAL) \
 if (fabs((NAME2) - (NAME1)) < FLT_EPSILON) { \

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;

trunk/psLib/src/math/psMinimizePowell.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-21 02:43:31 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-23 22:25:31 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -32 +32 @@
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* TYPE DEFINITIONS                                                          */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* GLOBAL VARIABLES                                                          */
-/* XXX: Do these conform to code standard?         */
-/*****************************************************************************/
-static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL;
-static psVector *myValue;
-static psVector *myError;
-/*****************************************************************************/
-/* FILE STATIC VARIABLES                                                     */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* FUNCTION IMPLEMENTATION - LOCAL                                           */
 /*****************************************************************************/
 
@@ -82 +63 @@
 } \
 
+/*****************************************************************************/
+/* TYPE DEFINITIONS                                                          */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* GLOBAL VARIABLES                                                          */
+/* XXX: Do these conform to code standard?         */
+/*****************************************************************************/
+static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL;
+static psVector *myValue;
+static psVector *myError;
+/*****************************************************************************/
+/* FILE STATIC VARIABLES                                                     */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* FUNCTION IMPLEMENTATION - LOCAL                                           */
+/*****************************************************************************/
 
 /******************************************************************************
@@ -402 +401 @@
     }
 
-    // XXX:Use psTraceGetLevel()
-    for (i=0;i<params->n;i++) {
-        psTrace(__func__, 6,
-                "(params, paramMask, line)[%d] is (%f %d %f)\n", i,
-                params->data.F32[i],
-                paramMask->data.U8[i],
-                line->data.F32[i]);
+    if (6 <= psTraceGetLevel(__func__)) {
+        for (i=0;i<params->n;i++) {
+            psTrace(__func__, 6, "(params, paramMask, line)[%d] is (%f %d %f)\n", i,
+                    params->data.F32[i], paramMask->data.U8[i], line->data.F32[i]);
+        }
     }
 

trunk/psLib/src/math/psPolynomial.c

@@ -7 +7 @@
 *  polynomials.  It also contains a Gaussian functions.
 *
-*  @version $Revision: 1.135 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-01-21 02:43:31 $
+*  @version $Revision: 1.136 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-23 22:25:31 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -14 +14 @@
 *  XXX: Should the "coeffErr[]" be used as well?  Bug ???.  Ignore coeffErr
 *
-*  XXX: In the various polyAlloc(n) functions, n is really the order of the
-*  polynomial plus 1.  To create a 2nd-order polynomial, n == 3.
 */
 /*****************************************************************************/
@@ -154 +152 @@
 
 /*****************************************************************************
-createChebyshevPolys(n): this routine takes as input the required order n,
+p_psCreateChebyshevPolys(n): this routine takes as input the required order n,
 and returns as output as a pointer to an array of n psPolynomial1D
 structures, corresponding to the first n Chebyshev polynomials.
@@ -162 +160 @@
 outer coefficients of the Chebyshev polynomials.
  
-XXX: From the poly nOrder/nTerm changem the maxChebyPoly here is still nTerms,
-not nOrder.
  *****************************************************************************/
-psPolynomial1D **createChebyshevPolys(psS32 numPolys)
+psPolynomial1D **p_psCreateChebyshevPolys(psS32 numPolys)
 {
     PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(numPolys, 1, NULL);
@@ -202 +198 @@
 }
 
-/*
-static psPolynomial1D **createChebyshevPolysOld(psS32 maxChebyPoly)
-{
-    PS_ASSERT_INT_NONNEGATIVE(maxChebyPoly, NULL);
- 
-    psPolynomial1D **chebPolys = NULL;
- 
-    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
-    for (psS32 i = 0; i < maxChebyPoly; i++) {
-        chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
-    }
- 
-    // Create the Chebyshev polynomials.
-    // Polynomial i has i-th order.
-    chebPolys[0]->coeff[0] = 1;
- 
-    // XXX: Bug 296
-    if (maxChebyPoly > 1) {
-        chebPolys[1]->coeff[1] = 1;
- 
-        for (psS32 i = 2; i < maxChebyPoly; i++) {
-            for (psS32 j = 0; j < chebPolys[i - 1]->nX; j++) {
-                chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j];
-            }
-            for (psS32 j = 0; j < chebPolys[i - 2]->nX; j++) {
-                chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
-            }
-        }
-    } else {
-        // XXX: Code this.
-        printf("WARNING: %u-order chebyshev polynomials not correctly implemented.\n", maxChebyPoly);
-    }
- 
-    return (chebPolys);
-}
-*/
 
 /*****************************************************************************
     Polynomial coefficients will be accessed in [w][x][y][z] fashion.
  *****************************************************************************/
-static psF64 ordPolynomial1DEval(psF64 x,
-                                 const psPolynomial1D* poly)
+static psF64 ordPolynomial1DEval(
+    psF64 x,
+    const psPolynomial1D* poly)
 {
     unsigned int loop_x = 0;
@@ -249 +210 @@
     psF64 xSum = 1.0;
 
-    psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4,
-            "---- Calling ordPolynomial1DEval(%lf)\n", x);
-    psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4,
-            "Polynomial order is %u\n", poly->nX);
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
+    psTrace(__func__, 4, "Polynomial order is %u\n", poly->nX);
     for (loop_x = 0; loop_x < poly->nX+1; loop_x++) {
-        psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4,
-                "Polynomial coeff[%u] is %lf\n", loop_x, poly->coeff[loop_x]);
+        psTrace(__func__, 4, "Polynomial coeff[%u] is %lf\n", loop_x, poly->coeff[loop_x]);
     }
 
@@ -262 +220 @@
             // XXX: If you set the tracelevel to 10 here, and later set the tracelevel to
             // 2 or higher in the test code, you get seg faults.
-            psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 8,
+            psTrace(__func__, 8,
                     "polysum+= sum*coeff [%lf+= (%lf * %lf)\n", polySum, xSum, poly->coeff[loop_x]);
             polySum += xSum * poly->coeff[loop_x];
@@ -269 +227 @@
     }
 
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
     return(polySum);
 }
@@ -284 +243 @@
     psVector *d;
 
-    // XXX: n should be nTerms here (for clarity).  Or get rid of the variable.
     unsigned int nTerms = 1 + poly->nX;
     unsigned int i;
@@ -338 +296 @@
     } else {
         // This is old code that does not use Clenshaw's formula.  Get rid of it.
-        psPolynomial1D **chebPolys = createChebyshevPolys(1 + poly->nX);
+        psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(1 + poly->nX);
 
         tmp = 0.0;
@@ -402 +360 @@
         maxChebyPoly = poly->nY;
     }
-    chebPolys = createChebyshevPolys(maxChebyPoly + 1);
+    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly + 1);
 
     for (loop_x = 0; loop_x < (1 + poly->nX); loop_x++) {
@@ -476 +434 @@
         maxChebyPoly = poly->nZ;
     }
-    chebPolys = createChebyshevPolys(maxChebyPoly + 1);
+    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly + 1);
 
     for (loop_x = 0; loop_x < (1 + poly->nX); loop_x++) {
@@ -567 +525 @@
         maxChebyPoly = poly->nT;
     }
-    // XXX: Add 1 since createChebyshevPolys() takes nTerms, not nOrder.
-    chebPolys = createChebyshevPolys(maxChebyPoly + 1);
+    // Add 1 since p_psCreateChebyshevPolys() takes nTerms, not nOrder.
+    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly + 1);
 
     for (loop_x = 0; loop_x < (1 + poly->nX); loop_x++) {
@@ -648 +606 @@
 /*****************************************************************************
     This routine must allocate memory for the polynomial structures.
-    XXX: replaces nterms variables with nOrder.  Lots of potential for bugs
-         here.  Probably should create separately named private variables
-         in these functions.
  *****************************************************************************/
-psPolynomial1D* psPolynomial1DAlloc(unsigned int n,
-                                    psPolynomialType type)
+psPolynomial1D* psPolynomial1DAlloc(
+    unsigned int n,
+    psPolynomialType type)
 {
     PS_ASSERT_INT_NONNEGATIVE(n, NULL);
-
-    unsigned int i = 0;
-    psPolynomial1D* newPoly = NULL;
-
-    newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D));
+    psU32 nOrder = n;
+    psPolynomial1D *newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D));
     psMemSetDeallocator(newPoly, (psFreeFunc) polynomial1DFree);
 
     newPoly->type = type;
-    newPoly->nX = n;
+    newPoly->nX = nOrder;
     newPoly->p_min = NAN;
     newPoly->p_max = NAN;
-    newPoly->coeff = psAlloc((n + 1) * sizeof(psF64));
-    newPoly->coeffErr = psAlloc((n + 1) * sizeof(psF64));
-    newPoly->mask = (psMaskType *)psAlloc((n + 1) * sizeof(psMaskType));
-    for (i = 0; i < (n + 1); i++) {
+    newPoly->coeff = psAlloc((nOrder + 1) * sizeof(psF64));
+    newPoly->coeffErr = psAlloc((nOrder + 1) * sizeof(psF64));
+    newPoly->mask = (psMaskType *)psAlloc((nOrder + 1) * sizeof(psMaskType));
+    for (psU32 i = 0; i < (nOrder + 1); i++) {
         newPoly->coeff[i] = 0.0;
         newPoly->coeffErr[i] = 0.0;
@@ -967 +920 @@
 // XXX: The output of this routine is always psF64 while 1D and 2D are
 // dependent on the input vectors.
-psVector *psPolynomial3DEvalVector(const psPolynomial3D *poly,
-                                   const psVector *x,
-                                   const psVector *y,
-                                   const psVector *z)
-
+psVector *psPolynomial3DEvalVector(
+    const psPolynomial3D *poly,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z)
 {
     PS_ASSERT_POLY_NON_NULL(poly, NULL);
@@ -1036 +989 @@
 }
 
-psF64 psPolynomial4DEval(const psPolynomial4D* poly,
-                         psF64 x,
-                         psF64 y,
-                         psF64 z,
-                         psF64 t)
+psF64 psPolynomial4DEval(
+    const psPolynomial4D* poly,
+    psF64 x,
+    psF64 y,
+    psF64 z,
+    psF64 t)
 {
     PS_ASSERT_POLY_NON_NULL(poly, NAN);
@@ -1056 +1010 @@
 }
 
-psVector *psPolynomial4DEvalVector(const psPolynomial4D *poly,
-                                   const psVector *x,
-                                   const psVector *y,
-                                   const psVector *z,
-                                   const psVector *t)
+psVector *psPolynomial4DEvalVector(
+    const psPolynomial4D *poly,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z,
+    const psVector *t)
 {
     PS_ASSERT_POLY_NON_NULL(poly, NULL);

trunk/psLib/src/math/psPolynomial.h

@@ -11 +11 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-19 23:58:47 $
+ *  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-23 22:25:31 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -69 +69 @@
 psPolynomialType;
 
-// XXX: These are incorrect names for the order of the polynomial.  We
-// keep them here temporarily so we can later sed replace them with the
-// correct names.
 /** One-dimensional polynomial */
 typedef struct
@@ -85 +82 @@
 psPolynomial1D;
 
-// XXX: These are incorrect names for the order of the polynomial.  We
-// keep them here temporarily so we can later sed replace them with the
-// correct names.
 /** Two-dimensional polynomial */
 typedef struct
@@ -100 +94 @@
 psPolynomial2D;
 
-// XXX: These are incorrect names for the order of the polynomial.  We
-// keep them here temporarily so we can later sed replace them with the
-// correct names.
 /** Three-dimensional polynomial */
 typedef struct
@@ -116 +107 @@
 psPolynomial3D;
 
-// XXX: These are incorrect names for the order of the polynomial.  We
-// keep them here temporarily so we can later sed replace them with the
-// correct names.
 /** Four-dimensional polynomial */
 typedef struct
@@ -301 +289 @@
 
 
-// XXX: Coding Standard
-psPolynomial1D **createChebyshevPolys(psS32 numPolys);
+
+
+/** Creates the specified number of chebyshev polys.
+ *
+ *  @return psPolynomial1D** The chebyshev polys.
+ *
+ */
+psPolynomial1D **p_psCreateChebyshevPolys(
+    psS32 numPolys
+);
 
 typedef struct

trunk/psLib/src/math/psSpline.c

@@ -7 +7 @@
 *  splines.
 *
-*  @version $Revision: 1.132 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-11-15 20:10:32 $
+*  @version $Revision: 1.133 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-23 22:25:31 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -90 +90 @@
 (F32).
  
-XXX: This algorithm is derived from the Numerical Recipes.
 XXX: use recycled vectors for internal data.
 XXX: do an F64 version?
@@ -111 +110 @@
     psF32 *Y = (psF32 *) & (y->data.F32[0]);
     //
-    // XXX: The second derivatives at the endpoints, undefined in the SDR,
+    // The second derivatives at the endpoints, undefined in the SDR,
     // are set in psConstants.h: PS_LEFT_SPLINE_DERIV, PS_RIGHT_SPLINE_DERIV.
     //
@@ -335 +334 @@
 LaGrange interpolation.
  *****************************************************************************/
-static psF32 interpolate1DF32(psF32 *domain,
-                              psF32 *range,
-                              unsigned int n,
-                              unsigned int order,
-                              psF32 x)
+static psF32 interpolate1DF32(
+    psF32 *domain,
+    psF32 *range,
+    unsigned int n,
+    unsigned int order,
+    psF32 x)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);

trunk/psLib/src/math/psStats.c

@@ -17 +17 @@
  *
  *
- *  @version $Revision: 1.159 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-23 20:44:29 $
+ *  @version $Revision: 1.160 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-23 22:25:31 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -2317 +2317 @@
     psTrace(__func__, 5, "(lower, upper, n) is (%f, %f, %d)\n", lower, upper, n);
     PS_ASSERT_INT_POSITIVE(n, NULL);
-    PS_FLOAT_COMPARE(lower, upper, NULL);
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(upper, lower, NULL);
 
     psS32 i = 0;                  // Loop index variable
@@ -2369 +2369 @@
     PS_ASSERT_VECTOR_NON_NULL(bounds, NULL);
     PS_ASSERT_VECTOR_TYPE(bounds, PS_TYPE_F32, NULL);
-    PS_INT_COMPARE(2, bounds->n, NULL);
+    PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(bounds->n, 2, NULL);
 
     psHistogram* newHist = NULL;        // The new histogram structure
@@ -2437 +2437 @@
     PS_ASSERT_PTR_NON_NULL(out->nums, -1);
     PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, ((out->nums->n)-1), -2);
-    PS_FLOAT_COMPARE(0.0, error, -3);
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(error, 0.0, -3);
+
     PS_ASSERT_FLOAT_WITHIN_RANGE(data, out->bounds->data.F32[0], out->bounds->data.F32[(out->bounds->n)-1], -4);
 
@@ -2756 +2757 @@
 
     if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) {
-        PS_FLOAT_COMPARE(stats->min, stats->max, stats);
+        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->max, stats->min, stats);
     }