Changeset 6388

Timestamp:

Feb 8, 2006, 2:52:18 PM (20 years ago)

Author:

gusciora

Message:

The psPlaneTransformFit() was completely rewritten and simplified. We now
simply call the psVectorPolynomial2DFit() routines.

File:

• trunk/psLib/src/astro/psCoord.c (modified) (8 diffs)

trunk/psLib/src/astro/psCoord.c

@@ -10 +10 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.109 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-02-08 21:52:05 $
+*  @version $Revision: 1.110 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-02-09 00:52:18 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -28 +28 @@
 #include "psErrorText.h"
 #include "psMatrix.h"
+#include "psMinimizePolyFit.h"
 #include <math.h>
 #include <float.h>
@@ -461 +462 @@
         out = psSphereAlloc();
     } else {
+        if ( !psMemCheckType(PS_DATA_SPHERE, outSphere) ) {
+            psError(PS_ERR_UNKNOWN, true, "outSphere must be of type *psSphere.\n");
+            return(NULL);
+        }
         out = outSphere;
-        // XXX: Do a memory checkpointer.
     }
 
@@ -565 +569 @@
                 for (psS32 t2x = 0 ; t2x < (1 + trans2->nX) ; t2x++) {
                     for (psS32 t2y = 0 ; t2y < (1 + trans2->nY) ; t2y++) {
-                        /* Possible debug-only macro which checks these coords?
-                        if ((t1x+t2x) >= orderX)
-                            printf("BAD 1\n");
-                        if ((t1y+t2y) >= orderY)
-                            printf("BAD 2\n");
-                        */
                         out->coeff[t1x+t2x][t1y+t2y]+= (trans1->coeff[t1x][t1y] * trans2->coeff[t2x][t2y]);
                     }
@@ -625 +623 @@
     // Allocate the new psPlaneTransform, if necessary.
     //
-    // XXX: rename, or verify, or recode, after the poly norder/nterm change.
 
     psPlaneTransform *myPT = NULL;
@@ -639 +636 @@
             // Initialize the new psPlaneTransform, if necessary.
             //
-            for (psS32 i = 0 ; i < orderX ; i++) {
-                for (psS32 j = 0 ; j < orderY ; j++) {
+            for (psS32 i = 0 ; i < orderX+1 ; i++) {
+                for (psS32 j = 0 ; j < orderY+1 ; j++) {
                     myPT->x->coeff[i][j] = 0.0;
                     myPT->x->mask[i][j] = 0;
@@ -769 +766 @@
 psPlaneTransformFit(trans, source, dest, nRejIter, sigmaClip)
  
-XXX: What about nRejIter?  Iterations?
-XXX: This code has problems with data that corresponds to a non-linear fit.
+XXX: This code ignores nRejIter and sigmaClip.
  *****************************************************************************/
 bool psPlaneTransformFit(
@@ -783 +779 @@
     PS_ASSERT_PTR_NON_NULL(dest, NULL);
 
-    psBool rc = true;
+    //
+    // Create the x and y vectors for the psVectorFitPolynomial2D() function.
+    //
     psS32 numCoords = PS_MIN(source->n, dest->n);
-    psS32 order = PS_MAX(trans->x->nX, trans->x->nY);
-    order = PS_MAX(order, trans->y->nX);
-    order = PS_MAX(order, trans->y->nY);
-    // XXX: Verify, or recode, after the poly norder/nterm change.
-    order++;
-
-    //
-    // Create fake polynomial to use in evaluation
-    //
-    // XXX: Verify this (poly order/nterm change)
-    psPolynomial2D *fakePoly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, order-1, order-1);
-    for (int i = 0; i < order; i++) {
-        for (int j = 0; j < order; j++) {
-            fakePoly->coeff[i][j] = 1.0;
-            fakePoly->mask[i][j] = 1;       // Mask all coefficients; unmask to evaluate
-        }
-    }
-
-    //
-    // Initialize the matrix and vectors
-    //
-    psS32 nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
-    psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
-    psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
-    psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
-    for (psS32 i = 0; i < nCoeff; i++) {
-        for (psS32 j = 0; j < nCoeff; j++) {
-            matrix->data.F64[i][j] = 0.0;
-        }
-        xVector->data.F64[i] = 0.0;
-        yVector->data.F64[i] = 0.0;
-    }
-
-    //
-    // Iterate over the grid points
-    //
+    psVector *xIn = psVectorAlloc(numCoords, PS_TYPE_F64);
+    psVector *yIn = psVectorAlloc(numCoords, PS_TYPE_F64);
+    psVector *xOut = psVectorAlloc(numCoords, PS_TYPE_F64);
+    psVector *yOut = psVectorAlloc(numCoords, PS_TYPE_F64);
     for (psS32 g = 0; g < numCoords; g++) {
-        // Iterate over the polynomial coefficients, accumulating the matrix and vectors
-
-        for (psS32 i = 0, ijIndex = 0; i < order; i++) {
-            for (psS32 j = 0; j < order - i; j++, ijIndex++) {
-                fakePoly->mask[i][j] = 0;
-                psF64 xIn = ((psPlane *) source->data[g])->x;
-                psF64 yIn = ((psPlane *) source->data[g])->y;
-                psF64 xOut = ((psPlane *) dest->data[g])->x;
-                psF64 yOut = ((psPlane *) dest->data[g])->y;
-                psF64 ijPoly = psPolynomial2DEval(fakePoly, xIn, yIn);
-                fakePoly->mask[i][j] = 1;
-
-                for (psS32 m = 0, mnIndex = 0; m < order; m++) {
-                    for (psS32 n = 0; n < order - m; n++, mnIndex++) {
-                        fakePoly->mask[m][n] = 0;
-                        psF64 mnPoly = psPolynomial2DEval(fakePoly, xIn, yIn);
-                        fakePoly->mask[m][n] = 1;
-
-                        matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
-                    }
-                }
-
-                xVector->data.F64[ijIndex] += ijPoly * xOut;
-                yVector->data.F64[ijIndex] += ijPoly * yOut;
-            }
-        }
-    }
-
-    //
-    // Solution via LU Decomposition
-    //
-    psVector *xSolution;
-    psVector *ySolution;
-    psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
-    psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
-    if (luMatrix == NULL) {
-        rc = false;
-        psError(PS_ERR_UNKNOWN, true, "psMatrixLUD() returned NULL.  Returning FALSE.\n");
-    } else {
-        xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
-        ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
-
-        if ((xSolution == NULL) || (ySolution == NULL)) {
-            rc = false;
-            psError(PS_ERR_UNKNOWN, true, "psMatrixLUSolve() returned NULL.  Returning FALSE.\n");
-        } else {
-            //
-            // Stuff coefficients into transformation
-            //
-            for (psS32 i = 0, ijIndex = 0; i < order; i++) {
-                for (psS32 j = 0; j < order - i; j++, ijIndex++) {
-                    trans->x->coeff[i][j] = xSolution->data.F64[ijIndex];
-                    trans->y->coeff[i][j] = ySolution->data.F64[ijIndex];
-                }
-            }
-        }
-    }
-
-    psFree(fakePoly);
-    psFree(permutation);
-    psFree(luMatrix);
-    psFree(xSolution);
-    psFree(ySolution);
-    psFree(matrix);
-    psFree(xVector);
-    psFree(yVector);
-
-    return(rc);
-}
+        xIn->data.F64[g] = ((psPlane *) source->data[g])->x;
+        yIn->data.F64[g] = ((psPlane *) source->data[g])->y;
+        xOut->data.F64[g] = ((psPlane *) dest->data[g])->x;
+        yOut->data.F64[g] = ((psPlane *) dest->data[g])->y;
+    }
+
+    trans->x = psVectorFitPolynomial2D(trans->x, NULL, 0, xOut, NULL, xIn, yIn);
+    trans->y = psVectorFitPolynomial2D(trans->y, NULL, 0, yOut, NULL, xIn, yIn);
+    psFree(xIn);
+    psFree(yIn);
+    psFree(xOut);
+    psFree(yOut);
+
+    if ((trans->x == NULL) || (trans->y == NULL)) {
+        psError( PS_ERR_UNKNOWN, true, "psVectorFitPolynomial2D() returned NULL: could not fit a 2-D polynomial to the data.\n");
+        return(false);
+    }
+
+    return(true);
+}
+