Changeset 4375 for trunk/psphot/src/psPolynomials.c

Timestamp:

Jun 24, 2005, 5:38:37 AM (21 years ago)

Author:

eugene

Message:

adding model tests

File:

• trunk/psphot/src/psPolynomials.c (modified) (5 diffs)

trunk/psphot/src/psPolynomials.c

@@ -1 +1 @@
 # include "psphot.h"
 
-// XXX EAM : this version uses myPoly->nX as Norder, not Nterms
-psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly,
-                                 const psVector *x,
-                                 const psVector *y)
-
-{
-    PS_POLY_CHECK_NULL(myPoly, NULL);
-    PS_VECTOR_CHECK_NULL(x, NULL);
-    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
-    PS_VECTOR_CHECK_NULL(y, NULL);
-    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
-
-    psVector *tmp;
-    psS32 vecLen=x->n;
-
-    // Determine the length of the output vector to by the minimum of the x,y vectors
-    if (y->n < vecLen) {
-        vecLen = y->n;
-    }
-
-    // Create output vector to return
-    tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
-
-    // Evaluate the polynomial at the specified points
-    for (psS32 i=0; i<vecLen; i++) {
-        tmp->data.F32[i] = Polynomial2DEval(myPoly,x->data.F32[i],y->data.F32[i]);
-    }
-
-    // Return output vector
-    return(tmp);
-}
-
-// XXX EAM : this version uses the F64 vectors
-psVector *Polynomial2DEvalVectorD(const psPolynomial2D *myPoly,
-                                  const psVector *x,
-                                  const psVector *y)
-
+// write out the terms of the given 1D polynomial
+void psPolynomial1DDump (psPolynomial1D *poly) {
+
+    for (int i = 0; i < poly->n + 1; i++) {
+        fprintf (stderr, "x^%d : %g +/- %g\n", i, poly->coeff[i], poly->coeffErr[i]);
+    }
+}    
+
+psF32 Polynomial1DEval_EAM(psF32 x, const psPolynomial1D* myPoly)
+{
+    psS32 loop_x = 0;
+    psF32 polySum = 0.0;
+    psF32 xSum = 1.0;
+
+    for (loop_x = 0; loop_x < myPoly->n + 1; loop_x++) {
+        if (myPoly->mask[loop_x] == 0) {
+            polySum += xSum * myPoly->coeff[loop_x];
+        }
+        xSum *= x;
+    }
+
+    return(polySum);
+}
+
+psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly,
+                                       const psVector *x)
 {
     PS_POLY_CHECK_NULL(myPoly, NULL);
     PS_VECTOR_CHECK_NULL(x, NULL);
     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
-    PS_VECTOR_CHECK_NULL(y, NULL);
-    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
 
     psVector *tmp;
-    psS32 vecLen=x->n;
-
-    // Determine the length of the output vector to by the minimum of the x,y vectors
-    if (y->n < vecLen) {
-        vecLen = y->n;
-    }
-
-    // Create output vector to return
-    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
-
-    // Evaluate the polynomial at the specified points
-    for (psS32 i=0; i<vecLen; i++) {
-        tmp->data.F64[i] = Polynomial2DEval(myPoly,x->data.F64[i],y->data.F64[i]);
-    }
-
-    // Return output vector
+
+    tmp = psVectorAlloc(x->n, PS_TYPE_F64);
+    for (psS32 i=0;i<x->n;i++) {
+        tmp->data.F64[i] = Polynomial1DEval_EAM(x->data.F64[i], myPoly);
+    }
+
     return(tmp);
 }
@@ -93 +70 @@
 }
 
-// XXX EAM : use Nterm = Norder + 1 definition  
-// the user requests a polynomial of order Norder
-psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder,
-                                    psPolynomialType type)
-{
-    PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL);
-    PS_INT_CHECK_NON_NEGATIVE(nYorder, NULL);
-
-    psS32 x = 0;
-    psS32 y = 0;
-    psS32 nXterm = nXorder + 1;
-    psS32 nYterm = nYorder + 1;
-    psPolynomial2D* newPoly = NULL;
-
-    newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D));
-    // p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree);
-    // XXX EAM : me, being lazy
-
-    newPoly->type = type;
-    newPoly->nX = nXorder;
-    newPoly->nY = nYorder;
-
-    newPoly->coeff = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
-    newPoly->coeffErr = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
-    newPoly->mask = (psU8 **)psAlloc(nXterm * sizeof(psU8 *));
-    for (x = 0; x < nXterm; x++) {
-        newPoly->coeff[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
-        newPoly->coeffErr[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
-        newPoly->mask[x] = (psU8 *)psAlloc(nYterm * sizeof(psU8));
-    }
-    for (x = 0; x < nXterm; x++) {
-        for (y = 0; y < nYterm; y++) {
-            newPoly->coeff[x][y] = 0.0;
-            newPoly->coeffErr[x][y] = 0.0;
-            newPoly->mask[x][y] = 0;
-        }
-    }
-    return(newPoly);
-}
-
 // XXX EAM : my alternate BuildSums1D
 static psVector *BuildSums1D(psVector* sums, 
@@ -153 +90 @@
         sums->data.F64[i] = xSum;
         xSum *= x;
-    }
-    return (sums);
-}
-
-// XXX EAM : BuildSums2D in analogy with BuildSums1D
-static psImage *BuildSums2D(psImage* sums,
-                            psF64 x,      psF64 y,
-                            psS32 nXterm, psS32 nYterm)
-{
-    psS32 nXsum = 0;
-    psS32 nYsum = 0;
-    psF64 xSum = 1.0;
-    psF64 ySum = 1.0;
-
-    nXsum = 2*nXterm;
-    nYsum = 2*nYterm;
-    if (sums == NULL) {
-        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
-    }
-    if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
-        psFree (sums);
-        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
-    }
-
-    ySum = 1.0;
-    for (int j = 0; j < nYsum; j++) {
-        xSum = ySum;
-        for (int i = 0; i < nXsum; i++) {
-            sums->data.F64[i][j] = xSum;
-            xSum *= x;
-        }
-        ySum *= y;
     }
     return (sums);
@@ -291 +196 @@
 }
 
+// ********************** 2D polynomial functions ******************
+
+// XXX EAM : this version uses myPoly->nX as Norder, not Nterms
+psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly,
+                                 const psVector *x,
+                                 const psVector *y)
+
+{
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+
+    psVector *tmp;
+    psS32 vecLen=x->n;
+
+    // Determine the length of the output vector to by the minimum of the x,y vectors
+    if (y->n < vecLen) {
+        vecLen = y->n;
+    }
+
+    // Create output vector to return
+    tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
+
+    // Evaluate the polynomial at the specified points
+    for (psS32 i=0; i<vecLen; i++) {
+        tmp->data.F32[i] = Polynomial2DEval(myPoly,x->data.F32[i],y->data.F32[i]);
+    }
+
+    // Return output vector
+    return(tmp);
+}
+
+// XXX EAM : this version uses the F64 vectors
+psVector *Polynomial2DEvalVectorD(const psPolynomial2D *myPoly,
+                                  const psVector *x,
+                                  const psVector *y)
+
+{
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+
+    psVector *tmp;
+    psS32 vecLen=x->n;
+
+    // Determine the length of the output vector to by the minimum of the x,y vectors
+    if (y->n < vecLen) {
+        vecLen = y->n;
+    }
+
+    // Create output vector to return
+    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
+
+    // Evaluate the polynomial at the specified points
+    for (psS32 i=0; i<vecLen; i++) {
+        tmp->data.F64[i] = Polynomial2DEval(myPoly,x->data.F64[i],y->data.F64[i]);
+    }
+
+    // Return output vector
+    return(tmp);
+}
+
+// XXX EAM : use Nterm = Norder + 1 definition  
+// the user requests a polynomial of order Norder
+psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder,
+                                    psPolynomialType type)
+{
+    PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nYorder, NULL);
+
+    psS32 x = 0;
+    psS32 y = 0;
+    psS32 nXterm = nXorder + 1;
+    psS32 nYterm = nYorder + 1;
+    psPolynomial2D* newPoly = NULL;
+
+    newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D));
+    // p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree);
+    // XXX EAM : me, being lazy
+
+    newPoly->type = type;
+    newPoly->nX = nXorder;
+    newPoly->nY = nYorder;
+
+    newPoly->coeff = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
+    newPoly->coeffErr = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
+    newPoly->mask = (psU8 **)psAlloc(nXterm * sizeof(psU8 *));
+    for (x = 0; x < nXterm; x++) {
+        newPoly->coeff[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
+        newPoly->coeffErr[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
+        newPoly->mask[x] = (psU8 *)psAlloc(nYterm * sizeof(psU8));
+    }
+    for (x = 0; x < nXterm; x++) {
+        for (y = 0; y < nYterm; y++) {
+            newPoly->coeff[x][y] = 0.0;
+            newPoly->coeffErr[x][y] = 0.0;
+            newPoly->mask[x][y] = 0;
+        }
+    }
+    return(newPoly);
+}
+
+// XXX EAM : BuildSums2D in analogy with BuildSums1D
+static psImage *BuildSums2D(psImage* sums,
+                            psF64 x,      psF64 y,
+                            psS32 nXterm, psS32 nYterm)
+{
+    psS32 nXsum = 0;
+    psS32 nYsum = 0;
+    psF64 xSum = 1.0;
+    psF64 ySum = 1.0;
+
+    nXsum = 2*nXterm;
+    nYsum = 2*nYterm;
+    if (sums == NULL) {
+        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
+    }
+    if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
+        psFree (sums);
+        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
+    }
+
+    ySum = 1.0;
+    for (int j = 0; j < nYsum; j++) {
+        xSum = ySum;
+        for (int i = 0; i < nXsum; i++) {
+            sums->data.F64[i][j] = xSum;
+            xSum *= x;
+        }
+        ySum *= y;
+    }
+    return (sums);
+}
+
 // XXX EAM : test version of 2d fitting
 psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly,
@@ -387 +430 @@
 }    
 
-// write out the terms of the given 1D polynomial
-void psPolynomial1DDump (psPolynomial1D *poly) {
-
-    for (int i = 0; i < poly->n + 1; i++) {
-        fprintf (stderr, "x^%d : %g +/- %g\n", i, poly->coeff[i], poly->coeffErr[i]);
-    }
-}    
-
 psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly,
                             const psVector* x,