Changeset 750 for trunk/psLib/src/math/psMinimize.c

Timestamp:

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

Author:

gusciora

Message:

...

File:

• trunk/psLib/src/math/psMinimize.c (modified) (2 diffs)

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 *