Index: trunk/psLib/src/astronomy/psAstrometry.c
===================================================================
--- trunk/psLib/src/astronomy/psAstrometry.c	(revision 3559)
+++ trunk/psLib/src/astronomy/psAstrometry.c	(revision 3598)
@@ -8,6 +8,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.62 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-03-30 02:21:14 $
+ *  @version $Revision: 1.63 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-03-31 23:01:46 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -51,106 +51,4 @@
 }
 
-/*****************************************************************************
-isProjectionLinear(): 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.
- *****************************************************************************/
-static psS32 isProjectionLinear(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);
-}
-
-/*****************************************************************************
-invertPlaneTransform(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.
- 
- *****************************************************************************/
-static psPlaneTransform *invertPlaneTransform(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);
-}
 
 static void FPAFree(psFPA* fpa)
@@ -914,5 +812,5 @@
 
     // generate an error if cell->toTP is not linear.
-    if (0 == isProjectionLinear(cell->toTP)) {
+    if (0 == p_psIsProjectionLinear(cell->toTP)) {
         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                 PS_ERRORTEXT_psAstrometry_NONLINEAR_TRANSFORM,
@@ -920,5 +818,5 @@
     }
 
-    TPtoCell = invertPlaneTransform(cell->toTP);
+    TPtoCell = p_psPlaneTransformLinearInvert(cell->toTP);
     cellCoord = psPlaneTransformApply(cellCoord, TPtoCell, tpCoord);
 
@@ -930,397 +828,3 @@
 
 
-/*****************************************************************************
-multiplyCoeffs(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);
-}
-
-
-
-
-/*****************************************************************************
-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 (isProjectionLinear((psPlaneTransform *) in)) {
-        printf("COOL: is linear\n");
-        return(invertPlaneTransform((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);
-}
+