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

Timestamp:

Jul 6, 2004, 8:07:53 PM (22 years ago)

Author:

gusciora

Message:

/

File:

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

trunk/psLib/src/math/psMinimize.c

@@ -600 +600 @@
 
 
-// NOTE: 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;
-/** @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
-
+/******************************************************************************
+p_psBuildSums1D(x, sums, polyOrder): this routine calculates the powers of
+input parameter "x" between 0 and input parameter polyOrder.  The result is
+returned as a psVector.
+ *****************************************************************************/
+psVector *p_psBuildSums1D(double x,
+                          int polyOrder)
+{
+    int       i = 0;
+    double    xSum = 0.0;
+    psVector *sums = NULL;
+
+    sums = psVectorAlloc(polyOrder+1, PS_TYPE_F32);
     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;
-        }
+        sums->data.F32[i] = xSum;
         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.
+    return(sums);
+}
+
+
+/******************************************************************************
+p_psBuildSums1D(x): this routine returns a psVector with "x" elements.  The
+values of the vector will be scaled uniformly between -1.0 and 1.0.
+ *****************************************************************************/
+psVector *psBuildImageScalingFactors(int x)
+
+{
     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++;
-            }
-        }
-    }
-}
-
+    psVector *imageScalingFactors = NULL;
+
+
+    imageScalingFactors = psVectorAlloc(x, PS_TYPE_F32);
+
+    for (i=0;i<x;i++) {
+        imageScalingFactors->data.F32[i] = (((float) 2*i) / ((float) x)) - 1.0;
+    }
+    return(imageScalingFactors);
+}
 
 /** @brief This routine checks if all polyOrder-th terms in the polyOrder-th
@@ -863 +758 @@
                      const psVector *restrict yErr)
 {
+    /*
+        int polyOrder = myPoly->n;
+        float **A;
+     
+     
+        // Numerical Recipes routines are all index offset 1.
+        B = (float *) psAlloc((polyOrder+2) * sizeof(float));
+        ludIndex = (int *) psAlloc((polyOrder+2) * sizeof(int));
+     
+        A = (float **) psAlloc((polyOrder+2) * sizeof(float *));
+        for(i=0;i<(polyOrder+2);i++) {
+            A[i] = (float *) psAlloc((polyOrder+2) * sizeof(float));
+        }
+     
+        // Initialize data structures.
+        for(i=0;i<(polyOrder+2);i++) {
+            B[i] = 0.0;
+            ludIndex[i] = 0;
+            for(j=0;j<(polyOrder+2);j++) {
+                A[i][j] = 0.0;
+            }
+        }
+     
+        for(k=1;k<(polyOrder+2);k++) {
+            for (i=0;i<x->n;i++) {
+                B[k]+= y->data.F32[i] * (pow(x->data.F32[i], k));
+            }
+        }
+        for(k=1;k<(polyOrder+2);k++) {
+            for(i=1;i<(polyOrder+2);i++) {
+                for (i=0;i<x->n;i++) {
+                    A[k][j]+= (pow(y->data.F32[i], k) * pow(x->data.F32[i], j));
+                }
+            }
+        }
+     
+     
+    */
     return(NULL);
 }
-
-
-
-
-
-