Index: /trunk/psLib/src/astro/psCoord.c
===================================================================
--- /trunk/psLib/src/astro/psCoord.c	(revision 5541)
+++ /trunk/psLib/src/astro/psCoord.c	(revision 5542)
@@ -10,6 +10,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.88 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-10-12 21:02:20 $
+*  @version $Revision: 1.89 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-11-18 19:39:29 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -57,56 +57,33 @@
 "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 there is now a general psPlaneTransformInvert() function, we
 should rename this.
  
-XXX: Use the ADD version which is based on determinants.
+To derive this transformation, start with the simple 2-by-2 matrix inversion
+based on discriminants.  This will invert the
+    (x_2, y_2) = matrix(a b c d) * vector (x, y)
+where there are no constant terms.  Then you substitute 
+    x_2 = x_1 - e
+    y_2 = y_1 - f
+for (x_2, y_2) to get the desired inverse.
  *****************************************************************************/
 psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform)
 {
-    PS_ASSERT_PTR_NON_NULL(transform, 0);
-    PS_ASSERT_PTR_NON_NULL(transform->x, 0);
-    PS_ASSERT_PTR_NON_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;
-
-    A = transform->x->coeff[0][0];
-    if (transform->x->nX >= 1) {
-        B = transform->x->coeff[1][0];
-    }
-    if (transform->x->nY >= 1) {
-        C = transform->x->coeff[0][1];
-    }
-    D = transform->y->coeff[0][0];
-    if (transform->y->nX >= 1) {
-        E = transform->y->coeff[1][0];
-    }
-    if (transform->y->nY >= 1) {
-        F = transform->y->coeff[0][1];
-    }
-
-    psPlaneTransform *out = psPlaneTransformAlloc(2, 2);
-
-    /* This is sample code from IfA.  It didn't work initially, and I did not
-       spend any time debugging it.
-
+    PS_ASSERT_PTR_NON_NULL(transform, NULL);
+    PS_ASSERT_PTR_NON_NULL(transform->x, NULL);
+    PS_ASSERT_PTR_NON_NULL(transform->y, NULL);
+    if ((transform->x->nX < 1) ||
+            (transform->x->nY < 1) ||
+            (transform->y->nX < 1) ||
+            (transform->y->nY < 1)) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: input transform is not invertible.");
+        return(NULL);
+    }
+
+    psPlaneTransform *out = psPlaneTransformAlloc(1, 1);
+
+    if (0) {
+        //This is sample code from IfA.  It didn't work initially, and I did not
+        //spend any time debugging it.
         psF64 a = transform->x->coeff[1][0];
         psF64 b = transform->x->coeff[0][1];
@@ -118,5 +95,5 @@
         psF64 invDet = 1.0 / (a * d - b * c); // Inverse of the determinant
 
-        // Not entirely sure why this works, but it appears to do so....................................!
+        // Not entirely sure why this works, but it appears to do so
         out->x->coeff[1][0] = invDet * a;
         out->x->coeff[0][1] = - invDet * b;
@@ -126,11 +103,29 @@
         out->x->coeff[0][0] = - invDet * (d * e + c * f);
         out->y->coeff[0][0] = - invDet * (b * e + a * f);
-    */
-    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));
+    }
+
+    if (1) {
+        psF64 A = transform->x->coeff[1][0];
+        psF64 B = transform->x->coeff[0][1];
+        psF64 C = transform->y->coeff[1][0];
+        psF64 D = transform->y->coeff[0][1];
+        psF64 E = transform->x->coeff[0][0];
+        psF64 F = transform->y->coeff[0][0];
+        psF64 invDet = 1.0 / (A * D - B * C);
+
+        out->x->coeff[1][0] = invDet *  D;
+        out->x->coeff[0][1] = invDet * -B;
+        out->y->coeff[1][0] = invDet * -C;
+        out->y->coeff[0][1] = invDet *  A;
+
+        //
+        // This comes from substituting the (x_2 = x_1 - e) and (y_2 = y_1 - f) equations:
+        //
+        out->x->coeff[0][0] = invDet * (B * F - D * E);
+        out->y->coeff[0][0] = invDet * (E * C - A * F);
+
+    }
+    //    printf("HMMM: out->x: (%f %f %f)\n", out->x->coeff[0][0], out->x->coeff[1][0], out->x->coeff[0][1]);
+    //    printf("HMMM: out->y: (%f %f %f)\n", out->y->coeff[0][0], out->y->coeff[1][0], out->y->coeff[0][1]);
 
     return(out);
@@ -143,14 +138,12 @@
  
 Returns:
-    1: if linear
-    0: otherwise
- 
-XXX: This should be a psBool
+    true: if linear
+    false: otherwise
  *****************************************************************************/
-psS32 p_psIsProjectionLinear(psPlaneTransform *transform)
-{
-    PS_ASSERT_PTR_NON_NULL(transform, 0);
-    PS_ASSERT_PTR_NON_NULL(transform->x, 0);
-    PS_ASSERT_PTR_NON_NULL(transform->y, 0);
+psBool p_psIsProjectionLinear(psPlaneTransform *transform)
+{
+    PS_ASSERT_PTR_NON_NULL(transform, false);
+    PS_ASSERT_PTR_NON_NULL(transform->x, false);
+    PS_ASSERT_PTR_NON_NULL(transform->y, false);
 
     for (psS32 i=0;i<(1 + transform->x->nX);i++) {
@@ -160,5 +153,5 @@
                         ((i == 0) && (j == 1)) ||
                         ((i == 1) && (j == 0)))) {
-                    return(0);
+                    return(false);
                 }
             }
@@ -172,5 +165,5 @@
                         ((i == 0) && (j == 1)) ||
                         ((i == 1) && (j == 0)))) {
-                    return(0);
+                    return(false);
                 }
             }
@@ -178,5 +171,5 @@
     }
 
-    return(1);
+    return(true);
 }
 
@@ -215,6 +208,6 @@
 
     psPlaneTransform *pt = psAlloc(sizeof(psPlaneTransform));
-    pt->x = psPolynomial2DAlloc(n1-1, n2-1, PS_POLYNOMIAL_ORD);
-    pt->y = psPolynomial2DAlloc(n1-1, n2-1, PS_POLYNOMIAL_ORD);
+    pt->x = psPolynomial2DAlloc(n1, n2, PS_POLYNOMIAL_ORD);
+    pt->y = psPolynomial2DAlloc(n1, n2, PS_POLYNOMIAL_ORD);
 
     psMemSetDeallocator(pt, (psFreeFunc) planeTransformFree);
@@ -243,7 +236,8 @@
 }
 
-psPlane* psPlaneTransformApply(psPlane* out,
-                               const psPlaneTransform* transform,
-                               const psPlane* coords)
+psPlane* psPlaneTransformApply(
+    psPlane* out,
+    const psPlaneTransform* transform,
+    const psPlane* coords)
 {
     PS_ASSERT_PTR_NON_NULL(transform, NULL);
@@ -256,15 +250,6 @@
     }
 
-    out->x = psPolynomial2DEval(
-                 transform->x,
-                 coords->x,
-                 coords->y
-             );
-
-    out->y = psPolynomial2DEval(
-                 transform->y,
-                 coords->x,
-                 coords->y
-             );
+    out->x = psPolynomial2DEval(transform->x, coords->x, coords->y);
+    out->y = psPolynomial2DEval(transform->y, coords->x, coords->y);
 
     return (out);
@@ -383,6 +368,8 @@
 
 
-psPlane* psProject(const psSphere* coord,
-                   const psProjection* projection)
+psPlane* p_psProject(
+    psPlane *outPlane,
+    const psSphere* coord,
+    const psProjection* projection)
 {
     PS_ASSERT_PTR_NON_NULL(coord, NULL);
@@ -393,5 +380,10 @@
 
     // Allocate return value
-    psPlane* out = psPlaneAlloc();
+    psPlane* out = NULL;
+    if (outPlane == NULL) {
+        out = psPlaneAlloc();
+    } else {
+        out = outPlane;
+    }
 
     // Convert to projection spherical coordinate system
@@ -435,6 +427,14 @@
 }
 
-psSphere* psDeproject(const psPlane* coord,
-                      const psProjection* projection)
+psPlane* psProject(const psSphere* coord,
+                   const psProjection* projection)
+{
+    return(p_psProject(NULL, coord, projection));
+}
+
+psSphere* p_psDeproject(
+    psSphere *outSphere,
+    const psPlane* coord,
+    const psProjection* projection)
 {
     PS_ASSERT_PTR_NON_NULL(coord, NULL);
@@ -445,5 +445,11 @@
 
     // Allocate return sphere structure
-    psSphere* out = psSphereAlloc();
+    psSphere *out = NULL;
+    if (outSphere == NULL) {
+        out = psSphereAlloc();
+    } else {
+        out = outSphere;
+        // XXX: Do a memory checkpointer.
+    }
 
     // Remove plate scales
@@ -490,4 +496,11 @@
 }
 
+psSphere* psDeproject(const psPlane* coord,
+                      const psProjection* projection)
+{
+    return(p_psDeproject(NULL, coord, projection));
+}
+
+
 /*****************************************************************************
 multiplyDPoly2D(trans1, trans2): Takes two 2-D polynomials as input and
@@ -813,4 +826,13 @@
 {
     PS_ASSERT_PTR_NON_NULL(in, NULL);
+    PS_ASSERT_PTR_NON_NULL(in->x, NULL);
+    PS_ASSERT_PTR_NON_NULL(in->y, NULL);
+    if ((in->x->nX < 1) ||
+            (in->x->nY < 1) ||
+            (in->y->nX < 1) ||
+            (in->y->nY < 1)) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: input transform is not invertible.");
+        return(NULL);
+    }
     //
     // If the transform is linear, then invert it exactly and return.
Index: /trunk/psLib/src/astro/psCoord.h
===================================================================
--- /trunk/psLib/src/astro/psCoord.h	(revision 5541)
+++ /trunk/psLib/src/astro/psCoord.h	(revision 5542)
@@ -11,6 +11,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-09-07 23:51:20 $
+*  @version $Revision: 1.47 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-11-18 19:39:29 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -290,4 +290,10 @@
 );
 
+psPlane* p_psProject(
+    psPlane *out,
+    const psSphere* coord,             ///< coordinate to project
+    const psProjection* projection     ///< parameters of the projection
+);
+
 /** Reverse projection of a linear coordinate to a spherical coordinate system
  *
@@ -295,4 +301,14 @@
  */
 psSphere* psDeproject(
+    const psPlane* coord,              ///< coordinate to project
+    const psProjection* projection     ///< parameters of the projection
+);
+
+/** Private reverse projection of a linear coordinate to a spherical coordinate system
+ *
+ *  @return psPlane*    projected coordinate
+ */
+psSphere* p_psDeproject(
+    psSphere *outSphere,
     const psPlane* coord,              ///< coordinate to project
     const psProjection* projection     ///< parameters of the projection
@@ -314,7 +330,8 @@
  *  the order of the projection
 */
-psS32 p_psIsProjectionLinear(
+psBool p_psIsProjectionLinear(
     psPlaneTransform *transform        ///< transform to test for linearity
 );
+
 
 /** inverts a given transformation.
Index: /trunk/psLib/src/math/psPolynomial.c
===================================================================
--- /trunk/psLib/src/math/psPolynomial.c	(revision 5541)
+++ /trunk/psLib/src/math/psPolynomial.c	(revision 5542)
@@ -7,6 +7,6 @@
 *  polynomials.  It also contains a Gaussian functions.
 *
-*  @version $Revision: 1.130 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-10-12 21:02:20 $
+*  @version $Revision: 1.131 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-11-18 19:39:29 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -644,6 +644,6 @@
                                      psPolynomialType type)
 {
-    PS_ASSERT_INT_POSITIVE(nX, NULL);
-    PS_ASSERT_INT_POSITIVE(nY, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
 
     unsigned int x = 0;
@@ -682,7 +682,7 @@
                                      psPolynomialType type)
 {
-    PS_ASSERT_INT_POSITIVE(nX, NULL);
-    PS_ASSERT_INT_POSITIVE(nY, NULL);
-    PS_ASSERT_INT_POSITIVE(nZ, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nZ, NULL);
 
     unsigned int x = 0;
@@ -731,8 +731,8 @@
                                      psPolynomialType type)
 {
-    PS_ASSERT_INT_POSITIVE(nX, NULL);
-    PS_ASSERT_INT_POSITIVE(nY, NULL);
-    PS_ASSERT_INT_POSITIVE(nZ, NULL);
-    PS_ASSERT_INT_POSITIVE(nT, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nZ, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nT, NULL);
 
     unsigned int x = 0;