Changeset 3540 for trunk/psLib/src/astronomy/psAstrometry.c

Timestamp:

Mar 29, 2005, 9:41:56 AM (21 years ago)

Author:

gusciora

Message:

Plane transform invert, combine, and fit.

File:

• trunk/psLib/src/astronomy/psAstrometry.c (modified) (6 diffs)

trunk/psLib/src/astronomy/psAstrometry.c

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-03-24 22:36:16 $
+ *  @version $Revision: 1.61 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-03-29 19:41:56 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -913 +913 @@
 
 
-
-
-
-/*****************************************************************************
-psPlaneTransformInvertHEY(out, in, region, nSamples): this is an earlier
-version of this function which doesnot separate the code for building and
-solving the matrixequations.
- 
-XXX: Delete this codeonce you know the other version works.
- *****************************************************************************/
-psPlaneTransform *psPlaneTransformInvertHEY(psPlaneTransform *out,
-        const psPlaneTransform *in,
-        psRegion *region,
-        int nSamples)
-{
-    PS_PTR_CHECK_NULL(in, NULL);
-    if (isProjectionLinear((psPlaneTransform *) in)) {
-        return(invertPlaneTransform((psPlaneTransform *) in));
-    }
-    psPlaneTransform *myPT = NULL;
-    psPlane *inCoord = psPlaneAlloc();
-    psPlane *outCoord = psPlaneAlloc();
-
-    PS_PTR_CHECK_NULL(region, NULL);
-    PS_INT_COMPARE(0, nSamples, NULL);
-
-    // XXX: Is this correct?
-    psS32 order = PS_MAX(in->x->nX, in->x->nY);
-
-    //
-    // 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))) {
-            psFree(out);
-            myPT = psPlaneTransformAlloc(order, order);
-        } else {
-            myPT = out;
-        }
-    }
-    // XXX: Initialize myPT?
-
-    //
-    // 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; // Set all coeffecients to 1
-            fakePoly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate
-        }
-    }
-
-    //
-    // Create a grid of xin,yin --> xout,yout
-    //
-    psVector *xIn = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
-    psVector *yIn = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
-    psVector *xOut = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
-    psVector *yOut = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
-
-    //
-    // Initialize the grid of points
-    //
-    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);
-            xOut->data.F32[yint*nSamples + xint] = inCoord->x;
-            yOut->data.F32[yint*nSamples + xint] = inCoord->y;
-            xIn->data.F32[yint*nSamples + xint] = outCoord->x;
-            yIn->data.F32[yint*nSamples + xint] = outCoord->y;
-        }
-    }
-
-    //
-    // Initialise 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 < nSamples*nSamples; 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 ijPoly = psDPolynomial2DEval(fakePoly, (psF64) xIn->data.F32[g], (psF64) yIn->data.F32[g]);
-                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, (psF64) xIn->data.F32[g], (psF64) yIn->data.F32[g]);
-                        fakePoly->mask[m][n] = 1;
-
-                        matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
-                    }
-                }
-
-                xVector->data.F64[ijIndex] += ijPoly * (psF64)xOut->data.F32[g];
-                yVector->data.F64[ijIndex] += ijPoly * (psF64)yOut->data.F32[g];
-            }
-        }
-    }
-
-    //
-    // 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
-
-    //
-    // Stuff coefficients into transformation
-    //
-    for (psS32 i = 0, ijIndex = 0; i < order; i++) {
-        for (psS32 j = 0; j < order - i; j++, ijIndex++) {
-            myPT->x->coeff[i][j] = xSolution->data.F64[ijIndex];
-            myPT->y->coeff[i][j] = ySolution->data.F64[ijIndex];
-        }
-    }
-
-    return(myPT);
-}
-
-
 /*****************************************************************************
 multiplyCoeffs(trans1, trans2): Takes two 2-D polynomials as input and
@@ -1100 +952 @@
 
 /*****************************************************************************
-psPlaneTransformCombineTmp(out, trans1, trans2)
+psPlaneTransformCombine(out, trans1, trans2)
  
 XXX: Much room for optimization.  Currently, we call the polyMultiply
 routine far too many times.
  *****************************************************************************/
-psPlaneTransform *psPlaneTransformCombineTmp(psPlaneTransform *out,
+psPlaneTransform *psPlaneTransformCombine(psPlaneTransform *out,
         const psPlaneTransform *trans1,
         const psPlaneTransform *trans2)
@@ -1240 +1092 @@
 
 /*****************************************************************************
-psPlaneTranformFitTmp(trans, source, dest, nRejIter, sigmaClip)
+psPlaneTranformFit(trans, source, dest, nRejIter, sigmaClip)
  
 XXX: What about nRejIter?  Iterations?
 XXX: Use static vectors for internal data.
  *****************************************************************************/
-bool psPlaneTranformFitTmp(psPlaneTransform *trans,
-                           const psArray *source,
-                           const psArray *dest,
-                           int nRejIter,
-                           float sigmaClip)
+bool psPlaneTranformFit(psPlaneTransform *trans,
+                        const psArray *source,
+                        const psArray *dest,
+                        int nRejIter,
+                        float sigmaClip)
 {
     PS_PTR_CHECK_NULL(trans, NULL);
@@ -1360 +1212 @@
 
 /*****************************************************************************
-psPlaneTransformInvertTmp(out, in, region, nSamples)
+psPlaneTransformInvert(out, in, region, nSamples)
  
 // XXX: Use static data structures.
  *****************************************************************************/
-psPlaneTransform *psPlaneTransformInvertTmp(psPlaneTransform *out,
+psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out,
         const psPlaneTransform *in,
         psRegion *region,
@@ -1448 +1300 @@
         }
     }
-    bool rc = psPlaneTranformFitTmp(myPT, inData, outData, 10, 100.0);
+    bool rc = psPlaneTranformFit(myPT, inData, outData, 10, 100.0);
 
     psFree(inCoord);