Changeset 750

Timestamp:

May 20, 2004, 3:04:45 PM (22 years ago)

Author:

gusciora

Message:

...

Location:

trunk/psLib/src

Files:

• dataManip/psFunctions.c (modified) (12 diffs)
• dataManip/psFunctions.h (modified) (1 diff)
• dataManip/psMinimize.c (modified) (2 diffs)
• math/psMinimize.c (modified) (2 diffs)
• math/psPolynomial.c (modified) (12 diffs)
• math/psPolynomial.h (modified) (1 diff)
• math/psSpline.c (modified) (12 diffs)
• math/psSpline.h (modified) (1 diff)

trunk/psLib/src/dataManip/psFunctions.c

@@ -24 +24 @@
 psGaussian(float x,
            float mean,
-           float stddev)
-{
-    return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+           float stddev,
+           int normal)
+{
+    float tmp = 0.0;
+
+    if (nomal == true) {
+        tmp = 1.0 / sqrtf(2.0 * PI * (stddev * stddev));
+        return(tmp * exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+    } else {
+        return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+    }
+}
+
+psVector *psGaussianDev(float mean,
+                        float sigma,
+                        int Npts)
+{
+    psVector *gauss = NULL;
+
+    gauss = psVectorAlloc(Npts, PS_TYPE_FLOAT);
 }
 
@@ -204 +221 @@
 
 /*****************************************************************************
-    Polynomial coefficients will be accessed in (x,y,z) fashion.
+    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
  *****************************************************************************/
 float
@@ -212 +229 @@
     int loop_x = 0;
     float polySum = 0.0;
+    float xSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->n;loop_x++) {
-        polySum+= x + myPoly->coeff[loop_x];
+        polySum+= xSum * myPoly->coeff[loop_x];
+        xSum*=x;
     }
 
@@ -229 +248 @@
     int loop_y = 0;
     float polySum = 0.0;
+    float xSum = 1.0;
+    float ySum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        xSum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
-            polySum+= x + myPoly->coeff[loop_x][loop_y];
-        }
+            polySum+= ySum * myPoly->coeff[loop_x][loop_y];
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -249 +273 @@
     int loop_z = 0;
     float polySum = 0.0;
+    float xSum = 1.0;
+    float ySum = 1.0;
+    float zSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        ySum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+            zSum = ySum;
             for (loop_z=0;loop_z<myPoly->nZ;loop_z++) {
-                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
-            }
-        }
+                polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z];
+                zSum*=z;
+            }
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -273 +305 @@
     int loop_z = 0;
     float polySum = 0.0;
+    float wSum = 1.0;
+    float xSum = 1.0;
+    float ySum = 1.0;
+    float zSum = 1.0;
 
     for (loop_w=0;loop_w<myPoly->nW;loop_w++) {
+        xSum = wSum;
         for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+            ySum = xSum;
             for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+                zSum = ySum;
                 for (loop_z=0;loop_z<myPoly->nZ;loop_z++) {
-                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    polySum+= zSum *
+                              myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    zSum*=z;
                 }
-            }
-        }
+                ySum*=y;
+            }
+            xSum*=x;
+        }
+        wSum*=w;
     }
 
@@ -461 +505 @@
 }
 
+/*****************************************************************************
+    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
+ *****************************************************************************/
 double
-psevalDPolynomial1D(double x,
-                    const psPolynomial1D *myPoly)
+psDEvalPolynomial1D(double x,
+                    const psDPolynomial1D *myPoly)
 {
     int loop_x = 0;
     double polySum = 0.0;
+    double xSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->n;loop_x++) {
-        polySum+= x + myPoly->coeff[loop_x];
+        polySum+= xSum * myPoly->coeff[loop_x];
+        xSum*=x;
     }
 
@@ -477 +526 @@
 
 double
-psevalDPolynomial2D(double x,
+psDEvalPolynomial2D(double x,
                     double y,
-                    const psPolynomial2D *myPoly)
+                    const psDPolynomial2D *myPoly)
 {
     int loop_x = 0;
     int loop_y = 0;
     double polySum = 0.0;
+    double xSum = 1.0;
+    double ySum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        xSum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
-            polySum+= x + myPoly->coeff[loop_x][loop_y];
-        }
+            polySum+= ySum * myPoly->coeff[loop_x][loop_y];
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -495 +549 @@
 
 double
-psevalDPolynomial3D(double x,
+psDEvalPolynomial3D(double x,
                     double y,
                     double z,
-                    const psPolynomial3D *myPoly)
+                    const psDPolynomial3D *myPoly)
 {
     int loop_x = 0;
@@ -504 +558 @@
     int loop_z = 0;
     double polySum = 0.0;
+    double xSum = 1.0;
+    double ySum = 1.0;
+    double zSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        ySum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+            zSum = ySum;
             for (loop_z=0;loop_z<myPoly->nZ;loop_z++) {
-                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
-            }
-        }
+                polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z];
+                zSum*=z;
+            }
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -517 +579 @@
 
 double
-psevalDPolynomial4D(double w,
+psDEvalPolynomial4D(double w,
                     double x,
                     double y,
                     double z,
-                    const psPolynomial4D *myPoly)
+                    const psDPolynomial4D *myPoly)
 {
     int loop_w = 0;
@@ -528 +590 @@
     int loop_z = 0;
     double polySum = 0.0;
+    double wSum = 1.0;
+    double xSum = 1.0;
+    double ySum = 1.0;
+    double zSum = 1.0;
 
     for (loop_w=0;loop_w<myPoly->nW;loop_w++) {
+        xSum = wSum;
         for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+            ySum = xSum;
             for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+                zSum = ySum;
                 for (loop_z=0;loop_z<myPoly->nZ;loop_z++) {
-                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    polySum+= zSum *
+                              myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    zSum*=z;
                 }
-            }
-        }
-    }
-
-    return(polySum);
-}
+                ySum*=y;
+            }
+            xSum*=x;
+        }
+        wSum*=w;
+    }
+
+    return(polySum);
+}

trunk/psLib/src/dataManip/psFunctions.h

@@ -91 +91 @@
 void psPolynomial4DFree(psPolynomial4D *myPoly ///< Polynomial to destroy
                        );
+
 
 /** Evaluate 1D polynomial */

trunk/psLib/src/dataManip/psMinimize.c

@@ -13 +13 @@
 #include "psFunctions.h"
 #include "psSort.h"
-
 #include "float.h"
 #include <math.h>
-
+// GUS: rewrite so there is no maximum order for the polynomials.
+#define MAX_POLY_ORDER 10
+#define MAX_POLYNOMIAL_TERMS  (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2)
+int MyInfoLevel = 0;
 /******************************************************************************
     This routine must minimize an arbitrary function.
@@ -43 +45 @@
 }
 
+/** @brief This procedure calculates various combinations of powers of x and y
+ *   and stores them in the data structure sums[][].  After it completes:
+ *          sums[i][j] == x^i * y^j
+ */
+void buildSums(double x,
+               double y,
+               /*@out@*/double sums[MAX_POLY_ORDER+1][MAX_POLY_ORDER+1],
+               int polyOrder)
+{
+    int         i = 0;          // loop index variable
+    int         j = 0;          // loop index variable
+    double       xSum = 0.0;    // The running sum of X terms
+    double       ySum = 0.0;    // The running sum of Y terms
+
+    xSum = 1.0;
+    ySum = 1.0;
+    for(i=0;i<=polyOrder;i++) {
+        ySum = xSum;
+        for(j=0;j<=polyOrder;j++) {
+            sums[i][j] = ySum;
+            ySum*= y;
+        }
+        xSum*= x;
+    }
+}
+
+/** @brief The coefficients of the matrix in equation (7) from the ADD will
+ * be very large if the x and y values are in the 0-511 range (ie: the sum y^7
+ * for all 0<y<512).  In order to avoid potential numerical instability, we
+ * added ability to scale those x,y values arbitrarily.  The following code
+ * creates a 1-D matrix imageScalingFactors[] which holds the scaled down
+ * values of x,y: the i-th element of imageScalingFactors[] contains the scaled
+ * down value for x=i, or y=i.
+ *
+ *     Input:
+ *     <ul>
+ *         <li>height
+ *         <li>width
+ *     </ul>
+ *
+ *     Output:
+ *     <ul>
+ *         <li>imageScalingFactors
+ *     </ul>
+ *
+ * @return error status (PsError) indicating error information, or NULL on
+ * success.
+ */
+void buildImageScalingFactors(int height,
+                              int width,
+                              float **imageScalingFactors)
+{
+    int maxDim = 0;             // The largest dimension of the image.
+    int i = 0;                  // loop index variable.
+
+    // Calculate the maximum dimensional extent of the image.
+    if (height > width) {
+        maxDim = height;
+    } else {
+        maxDim = width;
+    }
+
+
+    // Allocate memory for the output array.
+    *imageScalingFactors = (float *) psAlloc((maxDim+10) * sizeof(float));
+
+    // This code is somewhat arbitrary.  For an image with a height/width
+    // of 512x512, the scaling factors will be between 0.0-1.0.
+    for (i=0;i<maxDim;i++) {
+        (*imageScalingFactors)[i] = (((float) i) / ((float) maxDim)) - 0.5;
+        //        (*imageScalingFactors)[i] = ((float) i);
+    }
+}
+
+
+/** @brief buildPolyTerms(): this routine computes a 3-D array polyTerms[] that
+ *         holds terms for the polynomial that is used to model the sky
+ *         background.  We use this array primarily for convenience in many
+ *         computations involving that sky model polynomials. It is defined as:
+ * 
+ *             polyTerms[poly][i][0] = the power to which X is raised in the
+ *         i-th term of in an poly-order sky
+ *         background polynomial</P>.
+ *             polyTerms[poly][i][1] = the power to which Y is raised in the
+ *         i-th term of in an poly-order sky
+ *         background polynomial</P>.
+ * 
+ *    NOTE: the C_0 term defined in the ADD begins at i=2 in our data
+ *        structures (ie. the x/y powers of the i-th term in the sky model
+ *        polynomial are actually stored at polyTerms[][i+2][].  There are two
+ *        reasons for this.  First, there is a term prior to C_0 in equation
+ *        (7) of the ADD.  Second, our linear algebra codes assume data is
+ *        stored offset from index 1.
+ * 
+ *     Input:
+ *     <ul>
+ *         <li>polyTerms[][][]
+ *     </ul>
+ *
+ *     Output:
+ *     <ul>
+ *         <li>polyTerms[][][]
+ *     </ul>
+ *
+ * @return error status (PsError) indicating error information, or NULL on
+ * success.
+ */
+void buildPolyTerms(/*@out@*/ int polyTerms[MAX_POLY_ORDER+1][(MAX_POLYNOMIAL_TERMS+2)][2])
+{
+    int polyOrder=0;                    // loop index variable.
+    int i=0;                            // loop index variable.
+    int term = 0;                       // loop index variable.
+    int num=0;                          // loop index variable.
+
+    for(polyOrder=0;polyOrder<=MAX_POLY_ORDER;polyOrder++) {
+        // The following 4 terms should not be used in any of the subsequent
+        // computation.  We initialize them to zero in order to produce stable
+        // results for debugging purposes should they mistakenly be used.
+        polyTerms[polyOrder][0][0] = 0;
+        polyTerms[polyOrder][0][1] = 0;
+        polyTerms[polyOrder][1][0] = 0;
+        polyTerms[polyOrder][1][1] = 0;
+
+        // This code segment loops through each term i in the polynomial and
+        // calculates the power to which x/y are raised in that i-th term.
+        i=2;
+        for (term=0;term<=polyOrder;term++) {
+            for (num=0;num<=term;num++) {
+                polyTerms[polyOrder][i][0] = term-num;
+                polyTerms[polyOrder][i][1] = num;
+                if (MyInfoLevel > 2) {
+                    printf("%d-th order Sky polynomial term %d is x^%d y^%d\n",
+                           polyOrder, i,
+                           polyTerms[polyOrder][i][0], polyTerms[polyOrder][i][1]);
+                }
+                i++;
+            }
+        }
+    }
+}
+
+
+/** @brief This routine checks if all polyOrder-th terms in the polyOrder-th
+ * order sky background polynomial defined by the coefficients in the array B[]
+ * are consistent with zero.  If true, then *flag is set to 1.  Otherwise,
+ * *flag is set to 0.  The matrix inversion code in the middle of this
+ * procedure draws from Numerical Recipes in C page 48.
+ * 
+ *     Input:
+ *     <ul>
+ *         <li> A       This is the LUD decomposition of the original matrix A.
+ *         <li> N       The size of the matrix (plus 1, actually, since offset 1).
+ *         <li> indx    misc Numerical Recipes data structure.
+ *         <li> B       The coefficients of the sky polynomial.
+ *         <li> polyOrder The degree of the sky polynomial.
+ *     </ul>
+ *     Output:
+ *     <ul>
+ *         <li> *flag   Set this to 1 if we must recalculate the coefficients.
+ *     </ul>
+ *
+ * @return error status (PsError) indicating error information, or NULL on
+ * success.
+ */
+void polyOrderCheck(float **A,
+                    int N,
+                    int *indx,
+                    float *B,
+                    int polyOrder,
+                    int *flag)
+{
+    float     **y = NULL;  // This 2-D matrix will hold A^-1
+    float      *col = NULL;             // misc NumerRecipes data structure
+    float      *error=NULL;             // will hold the sqrt() of the
+    // diagonal of y[][].
+    int         i=0;                    // loop-index variable
+    int         j=0;                    // loop-index variable
+    int         numPolyTerms = 0;       // The number of terms in the
+    // polynomial.
+    int         lastTerm = 0;           // The index location of the first
+    // n-th order term in array B[].
+    int         firstTerm = 0;          // Index location of last such term.
+
+    // Allocate the necessary data structures for this procedure...
+    error = (float *) psAlloc((N + 1) * sizeof(float));
+    col = (float *) psAlloc((N + 1) * sizeof(float));
+    y = (float **) psAlloc((N + 1) * sizeof(float *));
+    for(i=1;i<=N;i++) {
+        y[i] = (float *) psAlloc((N + 1) * sizeof(float));
+    }
+
+    // Invert the matrix A and put the result in y[][].  This code is taken
+    // from Numerical Recipes in C page 48.
+    for(j=1;j<=N;j++) {
+        for(i=1;i<=N;i++) {
+            col[i] = 0.0;
+        }
+        col[j] = 1.0;
+        // GUS: substitue the LUD rotine
+        //        lubksb(A, N, indx, col);
+        for(i=1;i<=N;i++) {
+            y[i][j] = col[i];
+        }
+    }
+
+    // Determine where the first n-th order (in this comment, n equals
+    // polyOrder) polynomial term is stored in the matrix B[], and also were
+    // the last n-order term is stored.  Then we loop over all the n-order
+    // terms and check if they are consistent with zero.
+
+    numPolyTerms = (((polyOrder+1) * (polyOrder + 2)) / 2);
+    lastTerm = numPolyTerms + 1;
+    firstTerm = lastTerm - polyOrder;
+    *flag = 1;
+    for (i=firstTerm; i<=lastTerm; i++) {
+        error[i] = sqrtf(y[i][i]);
+        if (!((B[i]  <= (2.0f * error[i])) &&
+                ((-2.0f * error[i]) <= B[i]))) {
+            *flag = 0;
+        }
+    }
+
+    // Free all memory allocated in this routine.
+    psFree(error);
+    psFree(col);
+    for(j=1;j<=N;j++) {
+        psFree(y[j]);
+    }
+    psFree(y);
+}
+
+
+
+
+
+
+
+
+
 /******************************************************************************
-    This routine must minimize an arbotrary function.
+    This routine must fit a polynomial of degree myPoly to the data points
+    (x, y) and return the coefficients of that polynomial, as well as the
+    error for each data poiny (yErr).
  *****************************************************************************/
 psPolynomial1D *

trunk/psLib/src/math/psMinimize.c

@@ -13 +13 @@
 #include "psFunctions.h"
 #include "psSort.h"
-
 #include "float.h"
 #include <math.h>
-
+// GUS: rewrite so there is no maximum order for the polynomials.
+#define MAX_POLY_ORDER 10
+#define MAX_POLYNOMIAL_TERMS  (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2)
+int MyInfoLevel = 0;
 /******************************************************************************
     This routine must minimize an arbitrary function.
@@ -43 +45 @@
 }
 
+/** @brief This procedure calculates various combinations of powers of x and y
+ *   and stores them in the data structure sums[][].  After it completes:
+ *          sums[i][j] == x^i * y^j
+ */
+void buildSums(double x,
+               double y,
+               /*@out@*/double sums[MAX_POLY_ORDER+1][MAX_POLY_ORDER+1],
+               int polyOrder)
+{
+    int         i = 0;          // loop index variable
+    int         j = 0;          // loop index variable
+    double       xSum = 0.0;    // The running sum of X terms
+    double       ySum = 0.0;    // The running sum of Y terms
+
+    xSum = 1.0;
+    ySum = 1.0;
+    for(i=0;i<=polyOrder;i++) {
+        ySum = xSum;
+        for(j=0;j<=polyOrder;j++) {
+            sums[i][j] = ySum;
+            ySum*= y;
+        }
+        xSum*= x;
+    }
+}
+
+/** @brief The coefficients of the matrix in equation (7) from the ADD will
+ * be very large if the x and y values are in the 0-511 range (ie: the sum y^7
+ * for all 0<y<512).  In order to avoid potential numerical instability, we
+ * added ability to scale those x,y values arbitrarily.  The following code
+ * creates a 1-D matrix imageScalingFactors[] which holds the scaled down
+ * values of x,y: the i-th element of imageScalingFactors[] contains the scaled
+ * down value for x=i, or y=i.
+ *
+ *     Input:
+ *     <ul>
+ *         <li>height
+ *         <li>width
+ *     </ul>
+ *
+ *     Output:
+ *     <ul>
+ *         <li>imageScalingFactors
+ *     </ul>
+ *
+ * @return error status (PsError) indicating error information, or NULL on
+ * success.
+ */
+void buildImageScalingFactors(int height,
+                              int width,
+                              float **imageScalingFactors)
+{
+    int maxDim = 0;             // The largest dimension of the image.
+    int i = 0;                  // loop index variable.
+
+    // Calculate the maximum dimensional extent of the image.
+    if (height > width) {
+        maxDim = height;
+    } else {
+        maxDim = width;
+    }
+
+
+    // Allocate memory for the output array.
+    *imageScalingFactors = (float *) psAlloc((maxDim+10) * sizeof(float));
+
+    // This code is somewhat arbitrary.  For an image with a height/width
+    // of 512x512, the scaling factors will be between 0.0-1.0.
+    for (i=0;i<maxDim;i++) {
+        (*imageScalingFactors)[i] = (((float) i) / ((float) maxDim)) - 0.5;
+        //        (*imageScalingFactors)[i] = ((float) i);
+    }
+}
+
+
+/** @brief buildPolyTerms(): this routine computes a 3-D array polyTerms[] that
+ *         holds terms for the polynomial that is used to model the sky
+ *         background.  We use this array primarily for convenience in many
+ *         computations involving that sky model polynomials. It is defined as:
+ * 
+ *             polyTerms[poly][i][0] = the power to which X is raised in the
+ *         i-th term of in an poly-order sky
+ *         background polynomial</P>.
+ *             polyTerms[poly][i][1] = the power to which Y is raised in the
+ *         i-th term of in an poly-order sky
+ *         background polynomial</P>.
+ * 
+ *    NOTE: the C_0 term defined in the ADD begins at i=2 in our data
+ *        structures (ie. the x/y powers of the i-th term in the sky model
+ *        polynomial are actually stored at polyTerms[][i+2][].  There are two
+ *        reasons for this.  First, there is a term prior to C_0 in equation
+ *        (7) of the ADD.  Second, our linear algebra codes assume data is
+ *        stored offset from index 1.
+ * 
+ *     Input:
+ *     <ul>
+ *         <li>polyTerms[][][]
+ *     </ul>
+ *
+ *     Output:
+ *     <ul>
+ *         <li>polyTerms[][][]
+ *     </ul>
+ *
+ * @return error status (PsError) indicating error information, or NULL on
+ * success.
+ */
+void buildPolyTerms(/*@out@*/ int polyTerms[MAX_POLY_ORDER+1][(MAX_POLYNOMIAL_TERMS+2)][2])
+{
+    int polyOrder=0;                    // loop index variable.
+    int i=0;                            // loop index variable.
+    int term = 0;                       // loop index variable.
+    int num=0;                          // loop index variable.
+
+    for(polyOrder=0;polyOrder<=MAX_POLY_ORDER;polyOrder++) {
+        // The following 4 terms should not be used in any of the subsequent
+        // computation.  We initialize them to zero in order to produce stable
+        // results for debugging purposes should they mistakenly be used.
+        polyTerms[polyOrder][0][0] = 0;
+        polyTerms[polyOrder][0][1] = 0;
+        polyTerms[polyOrder][1][0] = 0;
+        polyTerms[polyOrder][1][1] = 0;
+
+        // This code segment loops through each term i in the polynomial and
+        // calculates the power to which x/y are raised in that i-th term.
+        i=2;
+        for (term=0;term<=polyOrder;term++) {
+            for (num=0;num<=term;num++) {
+                polyTerms[polyOrder][i][0] = term-num;
+                polyTerms[polyOrder][i][1] = num;
+                if (MyInfoLevel > 2) {
+                    printf("%d-th order Sky polynomial term %d is x^%d y^%d\n",
+                           polyOrder, i,
+                           polyTerms[polyOrder][i][0], polyTerms[polyOrder][i][1]);
+                }
+                i++;
+            }
+        }
+    }
+}
+
+
+/** @brief This routine checks if all polyOrder-th terms in the polyOrder-th
+ * order sky background polynomial defined by the coefficients in the array B[]
+ * are consistent with zero.  If true, then *flag is set to 1.  Otherwise,
+ * *flag is set to 0.  The matrix inversion code in the middle of this
+ * procedure draws from Numerical Recipes in C page 48.
+ * 
+ *     Input:
+ *     <ul>
+ *         <li> A       This is the LUD decomposition of the original matrix A.
+ *         <li> N       The size of the matrix (plus 1, actually, since offset 1).
+ *         <li> indx    misc Numerical Recipes data structure.
+ *         <li> B       The coefficients of the sky polynomial.
+ *         <li> polyOrder The degree of the sky polynomial.
+ *     </ul>
+ *     Output:
+ *     <ul>
+ *         <li> *flag   Set this to 1 if we must recalculate the coefficients.
+ *     </ul>
+ *
+ * @return error status (PsError) indicating error information, or NULL on
+ * success.
+ */
+void polyOrderCheck(float **A,
+                    int N,
+                    int *indx,
+                    float *B,
+                    int polyOrder,
+                    int *flag)
+{
+    float     **y = NULL;  // This 2-D matrix will hold A^-1
+    float      *col = NULL;             // misc NumerRecipes data structure
+    float      *error=NULL;             // will hold the sqrt() of the
+    // diagonal of y[][].
+    int         i=0;                    // loop-index variable
+    int         j=0;                    // loop-index variable
+    int         numPolyTerms = 0;       // The number of terms in the
+    // polynomial.
+    int         lastTerm = 0;           // The index location of the first
+    // n-th order term in array B[].
+    int         firstTerm = 0;          // Index location of last such term.
+
+    // Allocate the necessary data structures for this procedure...
+    error = (float *) psAlloc((N + 1) * sizeof(float));
+    col = (float *) psAlloc((N + 1) * sizeof(float));
+    y = (float **) psAlloc((N + 1) * sizeof(float *));
+    for(i=1;i<=N;i++) {
+        y[i] = (float *) psAlloc((N + 1) * sizeof(float));
+    }
+
+    // Invert the matrix A and put the result in y[][].  This code is taken
+    // from Numerical Recipes in C page 48.
+    for(j=1;j<=N;j++) {
+        for(i=1;i<=N;i++) {
+            col[i] = 0.0;
+        }
+        col[j] = 1.0;
+        // GUS: substitue the LUD rotine
+        //        lubksb(A, N, indx, col);
+        for(i=1;i<=N;i++) {
+            y[i][j] = col[i];
+        }
+    }
+
+    // Determine where the first n-th order (in this comment, n equals
+    // polyOrder) polynomial term is stored in the matrix B[], and also were
+    // the last n-order term is stored.  Then we loop over all the n-order
+    // terms and check if they are consistent with zero.
+
+    numPolyTerms = (((polyOrder+1) * (polyOrder + 2)) / 2);
+    lastTerm = numPolyTerms + 1;
+    firstTerm = lastTerm - polyOrder;
+    *flag = 1;
+    for (i=firstTerm; i<=lastTerm; i++) {
+        error[i] = sqrtf(y[i][i]);
+        if (!((B[i]  <= (2.0f * error[i])) &&
+                ((-2.0f * error[i]) <= B[i]))) {
+            *flag = 0;
+        }
+    }
+
+    // Free all memory allocated in this routine.
+    psFree(error);
+    psFree(col);
+    for(j=1;j<=N;j++) {
+        psFree(y[j]);
+    }
+    psFree(y);
+}
+
+
+
+
+
+
+
+
+
 /******************************************************************************
-    This routine must minimize an arbotrary function.
+    This routine must fit a polynomial of degree myPoly to the data points
+    (x, y) and return the coefficients of that polynomial, as well as the
+    error for each data poiny (yErr).
  *****************************************************************************/
 psPolynomial1D *

trunk/psLib/src/math/psPolynomial.c

@@ -24 +24 @@
 psGaussian(float x,
            float mean,
-           float stddev)
-{
-    return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+           float stddev,
+           int normal)
+{
+    float tmp = 0.0;
+
+    if (nomal == true) {
+        tmp = 1.0 / sqrtf(2.0 * PI * (stddev * stddev));
+        return(tmp * exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+    } else {
+        return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+    }
+}
+
+psVector *psGaussianDev(float mean,
+                        float sigma,
+                        int Npts)
+{
+    psVector *gauss = NULL;
+
+    gauss = psVectorAlloc(Npts, PS_TYPE_FLOAT);
 }
 
@@ -204 +221 @@
 
 /*****************************************************************************
-    Polynomial coefficients will be accessed in (x,y,z) fashion.
+    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
  *****************************************************************************/
 float
@@ -212 +229 @@
     int loop_x = 0;
     float polySum = 0.0;
+    float xSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->n;loop_x++) {
-        polySum+= x + myPoly->coeff[loop_x];
+        polySum+= xSum * myPoly->coeff[loop_x];
+        xSum*=x;
     }
 
@@ -229 +248 @@
     int loop_y = 0;
     float polySum = 0.0;
+    float xSum = 1.0;
+    float ySum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        xSum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
-            polySum+= x + myPoly->coeff[loop_x][loop_y];
-        }
+            polySum+= ySum * myPoly->coeff[loop_x][loop_y];
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -249 +273 @@
     int loop_z = 0;
     float polySum = 0.0;
+    float xSum = 1.0;
+    float ySum = 1.0;
+    float zSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        ySum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+            zSum = ySum;
             for (loop_z=0;loop_z<myPoly->nZ;loop_z++) {
-                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
-            }
-        }
+                polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z];
+                zSum*=z;
+            }
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -273 +305 @@
     int loop_z = 0;
     float polySum = 0.0;
+    float wSum = 1.0;
+    float xSum = 1.0;
+    float ySum = 1.0;
+    float zSum = 1.0;
 
     for (loop_w=0;loop_w<myPoly->nW;loop_w++) {
+        xSum = wSum;
         for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+            ySum = xSum;
             for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+                zSum = ySum;
                 for (loop_z=0;loop_z<myPoly->nZ;loop_z++) {
-                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    polySum+= zSum *
+                              myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    zSum*=z;
                 }
-            }
-        }
+                ySum*=y;
+            }
+            xSum*=x;
+        }
+        wSum*=w;
     }
 
@@ -461 +505 @@
 }
 
+/*****************************************************************************
+    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
+ *****************************************************************************/
 double
-psevalDPolynomial1D(double x,
-                    const psPolynomial1D *myPoly)
+psDEvalPolynomial1D(double x,
+                    const psDPolynomial1D *myPoly)
 {
     int loop_x = 0;
     double polySum = 0.0;
+    double xSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->n;loop_x++) {
-        polySum+= x + myPoly->coeff[loop_x];
+        polySum+= xSum * myPoly->coeff[loop_x];
+        xSum*=x;
     }
 
@@ -477 +526 @@
 
 double
-psevalDPolynomial2D(double x,
+psDEvalPolynomial2D(double x,
                     double y,
-                    const psPolynomial2D *myPoly)
+                    const psDPolynomial2D *myPoly)
 {
     int loop_x = 0;
     int loop_y = 0;
     double polySum = 0.0;
+    double xSum = 1.0;
+    double ySum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        xSum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
-            polySum+= x + myPoly->coeff[loop_x][loop_y];
-        }
+            polySum+= ySum * myPoly->coeff[loop_x][loop_y];
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -495 +549 @@
 
 double
-psevalDPolynomial3D(double x,
+psDEvalPolynomial3D(double x,
                     double y,
                     double z,
-                    const psPolynomial3D *myPoly)
+                    const psDPolynomial3D *myPoly)
 {
     int loop_x = 0;
@@ -504 +558 @@
     int loop_z = 0;
     double polySum = 0.0;
+    double xSum = 1.0;
+    double ySum = 1.0;
+    double zSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        ySum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+            zSum = ySum;
             for (loop_z=0;loop_z<myPoly->nZ;loop_z++) {
-                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
-            }
-        }
+                polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z];
+                zSum*=z;
+            }
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -517 +579 @@
 
 double
-psevalDPolynomial4D(double w,
+psDEvalPolynomial4D(double w,
                     double x,
                     double y,
                     double z,
-                    const psPolynomial4D *myPoly)
+                    const psDPolynomial4D *myPoly)
 {
     int loop_w = 0;
@@ -528 +590 @@
     int loop_z = 0;
     double polySum = 0.0;
+    double wSum = 1.0;
+    double xSum = 1.0;
+    double ySum = 1.0;
+    double zSum = 1.0;
 
     for (loop_w=0;loop_w<myPoly->nW;loop_w++) {
+        xSum = wSum;
         for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+            ySum = xSum;
             for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+                zSum = ySum;
                 for (loop_z=0;loop_z<myPoly->nZ;loop_z++) {
-                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    polySum+= zSum *
+                              myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    zSum*=z;
                 }
-            }
-        }
-    }
-
-    return(polySum);
-}
+                ySum*=y;
+            }
+            xSum*=x;
+        }
+        wSum*=w;
+    }
+
+    return(polySum);
+}

trunk/psLib/src/math/psPolynomial.h

@@ -91 +91 @@
 void psPolynomial4DFree(psPolynomial4D *myPoly ///< Polynomial to destroy
                        );
+
 
 /** Evaluate 1D polynomial */

trunk/psLib/src/math/psSpline.c

@@ -24 +24 @@
 psGaussian(float x,
            float mean,
-           float stddev)
-{
-    return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+           float stddev,
+           int normal)
+{
+    float tmp = 0.0;
+
+    if (nomal == true) {
+        tmp = 1.0 / sqrtf(2.0 * PI * (stddev * stddev));
+        return(tmp * exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+    } else {
+        return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+    }
+}
+
+psVector *psGaussianDev(float mean,
+                        float sigma,
+                        int Npts)
+{
+    psVector *gauss = NULL;
+
+    gauss = psVectorAlloc(Npts, PS_TYPE_FLOAT);
 }
 
@@ -204 +221 @@
 
 /*****************************************************************************
-    Polynomial coefficients will be accessed in (x,y,z) fashion.
+    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
  *****************************************************************************/
 float
@@ -212 +229 @@
     int loop_x = 0;
     float polySum = 0.0;
+    float xSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->n;loop_x++) {
-        polySum+= x + myPoly->coeff[loop_x];
+        polySum+= xSum * myPoly->coeff[loop_x];
+        xSum*=x;
     }
 
@@ -229 +248 @@
     int loop_y = 0;
     float polySum = 0.0;
+    float xSum = 1.0;
+    float ySum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        xSum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
-            polySum+= x + myPoly->coeff[loop_x][loop_y];
-        }
+            polySum+= ySum * myPoly->coeff[loop_x][loop_y];
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -249 +273 @@
     int loop_z = 0;
     float polySum = 0.0;
+    float xSum = 1.0;
+    float ySum = 1.0;
+    float zSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        ySum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+            zSum = ySum;
             for (loop_z=0;loop_z<myPoly->nZ;loop_z++) {
-                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
-            }
-        }
+                polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z];
+                zSum*=z;
+            }
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -273 +305 @@
     int loop_z = 0;
     float polySum = 0.0;
+    float wSum = 1.0;
+    float xSum = 1.0;
+    float ySum = 1.0;
+    float zSum = 1.0;
 
     for (loop_w=0;loop_w<myPoly->nW;loop_w++) {
+        xSum = wSum;
         for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+            ySum = xSum;
             for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+                zSum = ySum;
                 for (loop_z=0;loop_z<myPoly->nZ;loop_z++) {
-                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    polySum+= zSum *
+                              myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    zSum*=z;
                 }
-            }
-        }
+                ySum*=y;
+            }
+            xSum*=x;
+        }
+        wSum*=w;
     }
 
@@ -461 +505 @@
 }
 
+/*****************************************************************************
+    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
+ *****************************************************************************/
 double
-psevalDPolynomial1D(double x,
-                    const psPolynomial1D *myPoly)
+psDEvalPolynomial1D(double x,
+                    const psDPolynomial1D *myPoly)
 {
     int loop_x = 0;
     double polySum = 0.0;
+    double xSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->n;loop_x++) {
-        polySum+= x + myPoly->coeff[loop_x];
+        polySum+= xSum * myPoly->coeff[loop_x];
+        xSum*=x;
     }
 
@@ -477 +526 @@
 
 double
-psevalDPolynomial2D(double x,
+psDEvalPolynomial2D(double x,
                     double y,
-                    const psPolynomial2D *myPoly)
+                    const psDPolynomial2D *myPoly)
 {
     int loop_x = 0;
     int loop_y = 0;
     double polySum = 0.0;
+    double xSum = 1.0;
+    double ySum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        xSum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
-            polySum+= x + myPoly->coeff[loop_x][loop_y];
-        }
+            polySum+= ySum * myPoly->coeff[loop_x][loop_y];
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -495 +549 @@
 
 double
-psevalDPolynomial3D(double x,
+psDEvalPolynomial3D(double x,
                     double y,
                     double z,
-                    const psPolynomial3D *myPoly)
+                    const psDPolynomial3D *myPoly)
 {
     int loop_x = 0;
@@ -504 +558 @@
     int loop_z = 0;
     double polySum = 0.0;
+    double xSum = 1.0;
+    double ySum = 1.0;
+    double zSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        ySum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+            zSum = ySum;
             for (loop_z=0;loop_z<myPoly->nZ;loop_z++) {
-                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
-            }
-        }
+                polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z];
+                zSum*=z;
+            }
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -517 +579 @@
 
 double
-psevalDPolynomial4D(double w,
+psDEvalPolynomial4D(double w,
                     double x,
                     double y,
                     double z,
-                    const psPolynomial4D *myPoly)
+                    const psDPolynomial4D *myPoly)
 {
     int loop_w = 0;
@@ -528 +590 @@
     int loop_z = 0;
     double polySum = 0.0;
+    double wSum = 1.0;
+    double xSum = 1.0;
+    double ySum = 1.0;
+    double zSum = 1.0;
 
     for (loop_w=0;loop_w<myPoly->nW;loop_w++) {
+        xSum = wSum;
         for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+            ySum = xSum;
             for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+                zSum = ySum;
                 for (loop_z=0;loop_z<myPoly->nZ;loop_z++) {
-                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    polySum+= zSum *
+                              myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    zSum*=z;
                 }
-            }
-        }
-    }
-
-    return(polySum);
-}
+                ySum*=y;
+            }
+            xSum*=x;
+        }
+        wSum*=w;
+    }
+
+    return(polySum);
+}

trunk/psLib/src/math/psSpline.h

@@ -91 +91 @@
 void psPolynomial4DFree(psPolynomial4D *myPoly ///< Polynomial to destroy
                        );
+
 
 /** Evaluate 1D polynomial */