Changeset 3598 for trunk/psLib/src/astronomy/psCoord.c

Timestamp:

Mar 31, 2005, 1:01:46 PM (21 years ago)

Author:

desonia

Message:

cosmetic tweaks.

File:

• trunk/psLib/src/astronomy/psCoord.c (modified) (5 diffs)

trunk/psLib/src/astronomy/psCoord.c

@@ -10 +10 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-03-29 19:41:56 $
+*  @version $Revision: 1.61 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-03-31 23:01:46 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -26 +26 @@
 #include "psLogMsg.h"
 #include "psAstronomyErrors.h"
+#include "psAstrometry.h"
+#include "psMatrix.h"
 #include <math.h>
 #include <float.h>
@@ -59 +61 @@
 /*  FUNCTION IMPLEMENTATION - LOCAL                                          */
 /*****************************************************************************/
+/*****************************************************************************
+multiplyDPoly2D(trans1, trans2): Takes two 2-D polynomials as input and
+multiplies them.  Basically, for each non-zero coeff in the trans1 coeff[][]
+array, you must multiply by all non-zero coeffs in trans2.
+ 
+XXX: Inefficient in that the out polynomial is allocated every time.
+ *****************************************************************************/
+
+psDPolynomial2D *multiplyDPoly2D(psDPolynomial2D *trans1,
+                                 psDPolynomial2D *trans2)
+{
+    psS32 orderX = (trans1->nX + trans2->nX) - 1;
+    psS32 orderY = (trans1->nX + trans2->nX) - 1;
+
+    psDPolynomial2D *out = psDPolynomial2DAlloc(orderX, orderY, PS_POLYNOMIAL_ORD);
+    for (psS32 i = 0 ; i < out->nX; i++) {
+        for (psS32 j = 0 ; j < out->nY; j++) {
+            out->coeff[i][j] = 0.0;
+            out->mask[i][j] = 0;
+        }
+    }
+
+    for (psS32 t1x = 0 ; t1x < trans1->nX ; t1x++) {
+        for (psS32 t1y = 0 ; t1y < trans1->nY ; t1y++) {
+            if (0.0 != trans1->coeff[t1x][t1y]) {
+                for (psS32 t2x = 0 ; t2x < trans2->nX ; t2x++) {
+                    for (psS32 t2y = 0 ; t2y < trans2->nY ; t2y++) {
+                        out->coeff[t1x+t2x][t1y+t2y]+= (trans1->coeff[t1x][t1y] * trans2->coeff[t2x][t2y]);
+                    }
+                }
+            }
+        }
+    }
+    return(out);
+}
 
 /*****************************************************************************/
@@ -66 +103 @@
 {
     // There are non dynamic allocated items
+}
+
+/*****************************************************************************
+p_psPlaneTransformLinearInvert(transform): : this is a private function which
+simply inverts the supplied psPlaneTransform transform.  It assumes that
+"transform" is linear.
+ 
+This program assumes that the inverse of the following linear equations:
+        X2 = A + (B * X1) + (C * Y1);
+        Y2 = D + (E * X1) + (F * Y1);
+is
+        Y1 = (Y2 - ((E/B) * X2) - D + ((E*A)/B)) / (F - ((C*E)/B));
+        X1 = (Y2 - ((F/C) * X2) - D + ((F*A)/C)) / (E - ((F*B)/C));
+or
+ X1 = (-D + ((F*A)/C)) / (E - ((F*B)/C)) +
+      (X2 * -((F/C) / (E - ((F*B)/C)))) +
+             (Y2 * (1.0 / (E - ((F*B)/C))));
+        Y1 = (-D + ((E*A)/B))/(F - ((C*E)/B)) +
+             (X2 * -((E/B) / (F - ((C*E)/B)))) +
+             (Y2 * (1.0 / (F - ((C*E)/B))));
+ 
+XXX: Since thre is now a general psPlaneTransformInvert() function, we
+should rename this.
+ 
+ *****************************************************************************/
+psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform)
+{
+    PS_PTR_CHECK_NULL(transform, 0);
+    PS_PTR_CHECK_NULL(transform->x, 0);
+    PS_PTR_CHECK_NULL(transform->y, 0);
+
+    psF64 A = 0.0;
+    psF64 B = 0.0;
+    psF64 C = 0.0;
+    psF64 D = 0.0;
+    psF64 E = 0.0;
+    psF64 F = 0.0;
+
+    // XXX: Test this for correctness.
+    A = transform->x->coeff[0][0];
+    if (transform->x->nX >= 2) {
+        B = transform->x->coeff[1][0];
+    }
+    if (transform->x->nY >= 2) {
+        C = transform->x->coeff[0][1];
+    }
+    D = transform->y->coeff[0][0];
+    if (transform->y->nX >= 2) {
+        E = transform->y->coeff[1][0];
+    }
+    if (transform->y->nY >= 2) {
+        F = transform->y->coeff[0][1];
+    }
+
+    // XXX: Use the constructor here.
+    psPlaneTransform *out = psPlaneTransformAlloc(2, 2);
+
+    out->x->coeff[0][0] = -D + ((F*A)/C) / (E - ((F*B)/C));
+    out->x->coeff[1][0] = -(F/C) / (E - ((F*B)/C));
+    out->x->coeff[0][1] =  1.0 / (E - ((F*B)/C));
+    out->y->coeff[0][0] = -D + ((E*A)/B) / (F - ((C*E)/B));
+    out->y->coeff[1][0] = -(E/B) / (F - ((C*E)/B));
+    out->y->coeff[0][1] =  1.0 / (F - ((C*E)/B));
+
+    return(out);
+}
+
+/*****************************************************************************
+p_psIsProjectionLinear(): this is a private function which simply determines
+if the supplied psPlaneTransform transform is linear: if any of the
+cooefficients of order 2 are higher are non-zero, then it is not linear.
+ *****************************************************************************/
+psS32 p_psIsProjectionLinear(psPlaneTransform *transform)
+{
+    PS_PTR_CHECK_NULL(transform, 0);
+    PS_PTR_CHECK_NULL(transform->x, 0);
+    PS_PTR_CHECK_NULL(transform->y, 0);
+
+    for (psS32 i=0;i<(transform->x->nX);i++) {
+        for (psS32 j=0;j<(transform->x->nY);j++) {
+            if (transform->x->coeff[i][j] != 0.0) {
+                if (!(((i == 0) && (j == 0)) ||
+                        ((i == 0) && (j == 1)) ||
+                        ((i == 1) && (j == 0)))) {
+                    return(0);
+                }
+            }
+        }
+    }
+
+    for (psS32 i=0;i<(transform->y->nX);i++) {
+        for (psS32 j=0;j<(transform->y->nY);j++) {
+            if (transform->y->coeff[i][j] != 0.0) {
+                if (!(((i == 0) && (j == 0)) ||
+                        ((i == 0) && (j == 1)) ||
+                        ((i == 1) && (j == 0)))) {
+                    return(0);
+                }
+            }
+        }
+    }
+
+    return(1);
 }
 
@@ -721 +861 @@
 }
 
+
+
+/*****************************************************************************
+psPlaneTransformCombine(out, trans1, trans2)
+ 
+XXX: Much room for optimization.  Currently, we call the polyMultiply
+routine far too many times.
+ *****************************************************************************/
+psPlaneTransform *psPlaneTransformCombine(psPlaneTransform *out,
+        const psPlaneTransform *trans1,
+        const psPlaneTransform *trans2)
+{
+    PS_PTR_CHECK_NULL(trans1, NULL);
+    PS_PTR_CHECK_NULL(trans2, NULL);
+
+    //
+    // Determine the size of the new psPlaneTransform.
+    //
+    // PS_MAX(  Number of x terms in T2->x * number of x terms in T1->x,
+    //          Number of y terms in T2->x * number of x terms in T1->y,
+    psS32 orderXnX = PS_MAX((trans2->x->nX * trans1->x->nX),
+                            (trans2->x->nY * trans1->y->nX));
+    psS32 orderXnY = PS_MAX((trans2->x->nX * trans1->x->nY),
+                            (trans2->x->nY * trans1->y->nY));
+
+    psS32 orderYnX = PS_MAX((trans2->y->nX * trans1->x->nX),
+                            (trans2->y->nY * trans1->y->nX));
+    psS32 orderYnY = PS_MAX((trans2->y->nX * trans1->x->nY),
+                            (trans2->y->nY * trans1->y->nY));
+    psS32 orderX = PS_MAX(orderXnX, orderYnX);
+    psS32 orderY = PS_MAX(orderXnY, orderYnY);
+
+    //
+    // Allocate the new psPlaneTransform, if necessary.
+    //
+    psPlaneTransform *myPT = NULL;
+    if (out == NULL) {
+        myPT = psPlaneTransformAlloc(orderX, orderY);
+    } else {
+        if ((out->x->nX == orderX) && (out->x->nY == orderY) &&
+                (out->y->nX == orderX) && (out->y->nY == orderY)) {
+            myPT = out;
+        } else {
+            psFree(out);
+            myPT = psPlaneTransformAlloc(orderX, orderY);
+        }
+    }
+
+    //
+    // Initialize the new psPlaneTransform, if necessary.
+    //
+    for (psS32 i = 0 ; i < orderX ; i++) {
+        for (psS32 j = 0 ; j < orderY ; j++) {
+            myPT->x->coeff[i][j] = 0.0;
+            myPT->x->mask[i][j] = 0;
+            myPT->y->coeff[i][j] = 0.0;
+            myPT->y->mask[i][j] = 0;
+        }
+    }
+
+    //
+    // For each term (a * x^i * y^j) in trans2, we substitute the appropriate
+    // equation from trans1, and raise it to the appropriate power.  This is
+    // done via the multiplyDPoly2D().  The result is a polynomial (currPoly)
+    // and its coefficients are added into the myPT coeff matrix.
+    //
+    // XXX: This is horribly inefficient in that the trans1 polys are repeatedly
+    // multiplied against themselves.  This can easily be improved.
+    //
+    for (psS32 t2x = 0 ; t2x < trans2->x->nX ; t2x++) {
+        for (psS32 t2y = 0 ; t2y < trans2->x->nY ; t2y++) {
+            psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
+            currPoly->coeff[0][0] = 1.0;
+            currPoly->mask[0][0] = 0;
+            psDPolynomial2D *newPoly = NULL;
+
+            if (trans2->x->mask[t2x][t2y] == 0) {
+
+                // Must raise trans1->y to the t2y-power.
+                for (psS32 c = 0 ; c < t2y; c++) {
+                    newPoly = multiplyDPoly2D(currPoly, trans1->y);
+                    psFree(currPoly);
+                    currPoly = newPoly;
+                }
+
+                // Must raise trans1->x to the t2x-power.
+                for (psS32 c = 0 ; c < t2x; c++) {
+                    newPoly = multiplyDPoly2D(currPoly, trans1->x);
+                    psFree(currPoly);
+                    currPoly = newPoly;
+                }
+
+                // Set the appropriate coeffs in myPT->x
+                for (psS32 i = 0 ; i < currPoly->nX ; i++) {
+                    for (psS32 j = 0 ; j < currPoly->nY ; j++) {
+                        myPT->x->coeff[i][j]+= currPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y];
+                    }
+                }
+            }
+            psFree(currPoly);
+        }
+    }
+
+
+
+    for (psS32 t2x = 0 ; t2x < trans2->y->nX ; t2x++) {
+        for (psS32 t2y = 0 ; t2y < trans2->y->nY ; t2y++) {
+            psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
+            currPoly->coeff[0][0] = 1.0;
+            currPoly->mask[0][0] = 0;
+            psDPolynomial2D *newPoly = NULL;
+
+            if (trans2->y->mask[t2x][t2y] == 0) {
+
+                // Must raise trans1->y to the t2y-power.
+                for (psS32 c = 0 ; c < t2y; c++) {
+                    newPoly = multiplyDPoly2D(currPoly, trans1->y);
+                    psFree(currPoly);
+                    currPoly = newPoly;
+                }
+
+                // Must raise trans1->x to the t2x-power.
+                for (psS32 c = 0 ; c < t2x; c++) {
+                    newPoly = multiplyDPoly2D(currPoly, trans1->x);
+                    psFree(currPoly);
+                    currPoly = newPoly;
+                }
+
+                // Set the appropriate coeffs in myPT->x
+                for (psS32 i = 0 ; i < currPoly->nX ; i++) {
+                    for (psS32 j = 0 ; j < currPoly->nY ; j++) {
+                        myPT->y->coeff[i][j]+= currPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y];
+                    }
+                }
+            }
+            psFree(currPoly);
+        }
+    }
+
+    return(myPT);
+}
+
+/*****************************************************************************
+psPlaneTransformFit(trans, source, dest, nRejIter, sigmaClip)
+ 
+XXX: What about nRejIter?  Iterations?
+XXX: Use static vectors for internal data.
+ *****************************************************************************/
+bool psPlaneTransformFit(psPlaneTransform *trans,
+                         const psArray *source,
+                         const psArray *dest,
+                         int nRejIter,
+                         float sigmaClip)
+{
+    PS_PTR_CHECK_NULL(trans, NULL);
+    PS_PTR_CHECK_NULL(source, NULL);
+    PS_PTR_CHECK_NULL(dest, NULL);
+
+    psS32 numCoords = PS_MIN(source->n, dest->n);
+    // This is not really necessary because of above conditionals.
+    psS32 order = PS_MAX(trans->x->nX, trans->x->nY);
+
+    //
+    // Create fake polynomial to use in evaluation
+    //
+    psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
+    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
+    //
+    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 = psDPolynomial2DEval(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 = psDPolynomial2DEval(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 *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
+    psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
+    psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
+    psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
+
+    //
+    // XXX: Should check the output of the matrix routines and return false if bad.
+    //
+
+    //
+    // 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(true);
+}
+
+
+/*****************************************************************************
+psPlaneTransformInvert(out, in, region, nSamples)
+ 
+// XXX: Use static data structures.
+ *****************************************************************************/
+psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out,
+        const psPlaneTransform *in,
+        psRegion *region,
+        int nSamples)
+{
+    PS_PTR_CHECK_NULL(in, NULL);
+    //
+    // If the transform is linear, then invert it exactly and return.
+    //
+    if (p_psIsProjectionLinear((psPlaneTransform *) in)) {
+        printf("COOL: is linear\n");
+        return(p_psPlaneTransformLinearInvert((psPlaneTransform *) in));
+    }
+    PS_PTR_CHECK_NULL(region, NULL);
+    PS_INT_COMPARE(1, nSamples, NULL);
+
+    // Ensure that the input transformation is symmetrical.
+    if ((in->x->nX != in->x->nY) ||
+            (in->y->nX != in->y->nY) ||
+            (in->x->nX != in->y->nX)) {
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Input transformation must have same nX==nY.");
+    }
+    psS32 order = PS_MAX(in->x->nX, in->x->nY);
+
+    psPlaneTransform *myPT = NULL;
+    psPlane *inCoord = psPlaneAlloc();
+    psPlane *outCoord = psPlaneAlloc();
+
+    //
+    // Allocate a new psPlaneTransform if "out" is NULL, or has the wrong size.
+    //
+    if (out == NULL) {
+        myPT = psPlaneTransformAlloc(order, order);
+    } else {
+        if ((out->x->nX == order) && (out->x->nY == order) &&
+                (out->y->nX == order) && (out->y->nY == order)) {
+            myPT = out;
+        } else {
+            psFree(out);
+            myPT = psPlaneTransformAlloc(order, order);
+        }
+    }
+
+    //
+    // Copy the input transform to myPT.
+    //
+    for (psS32 i = 0 ; i < in->x->nX ; i++) {
+        for (psS32 j = 0 ; j < in->x->nY ; j++) {
+            myPT->x->coeff[i][j] = in->x->coeff[i][j];
+        }
+    }
+    for (psS32 i = 0 ; i < in->y->nX ; i++) {
+        for (psS32 j = 0 ; j < in->y->nY ; j++) {
+            myPT->y->coeff[i][j] = in->y->coeff[i][j];
+        }
+    }
+
+    //
+    // Create a grid of xin,yin --> xout,yout
+    //
+    psArray *inData = psArrayAlloc(nSamples * nSamples);
+    psArray *outData = psArrayAlloc(nSamples * nSamples);
+    for (psS32 i = 0 ; i < inData->n; i++) {
+        inData->data[i] = (psPtr *) psPlaneAlloc();
+        outData->data[i] = (psPtr *) psPlaneAlloc();
+    }
+
+    //
+    // Initialize the grid.
+    //
+    psS32 cnt = 0;
+    for (int yint = 0; yint < nSamples; yint++) {
+        inCoord->y = region->y0 + ((psF32) yint) * ((region->y1 - region->y0) / ((psF32) nSamples));
+        for (int xint = 0; xint < nSamples; xint++) {
+            inCoord->x = region->x0 + ((psF32) xint) * ((region->x1 - region->x0) / ((psF32) nSamples));
+            (void)psPlaneTransformApply(outCoord, in, inCoord);
+
+            ((psPlane *) outData->data[cnt])->x = inCoord->x;
+            ((psPlane *) outData->data[cnt])->y = inCoord->y;
+            ((psPlane *) inData->data[cnt])->x = outCoord->x;
+            ((psPlane *) inData->data[cnt])->y = outCoord->y;
+
+            cnt++;
+        }
+    }
+    bool rc = psPlaneTransformFit(myPT, inData, outData, 10, 100.0);
+
+    psFree(inCoord);
+    psFree(outCoord);
+    psFree(inData);
+    psFree(outData);
+
+    if (rc == true) {
+        return(myPT);
+    }
+
+    // XXX: Generate an error message, or warning message.
+    return(NULL);
+}