Changeset 6389

Timestamp:

Feb 8, 2006, 3:18:53 PM (20 years ago)

Author:

gusciora

Message:

Complete rewrite of psPlaneTransformInvert(). I now call the vectorFit
poly routines.

File:

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

trunk/psLib/src/astro/psCoord.c

@@ -10 +10 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.110 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-02-09 00:52:18 $
+*  @version $Revision: 1.111 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-02-09 01:18:53 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -766 +766 @@
 psPlaneTransformFit(trans, source, dest, nRejIter, sigmaClip)
  
-XXX: This code ignores nRejIter and sigmaClip.
+XXX: This code ignores nRejIter and sigmaClip.  We must call the ClipFit
+routines instead.
  *****************************************************************************/
 bool psPlaneTransformFit(
@@ -824 +825 @@
     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)) {
+    PS_ASSERT_INT_LARGER_THAN(nSamples, 0, NULL);
+
+    // Reject a trivially non-invertible case.
+    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);
     }
+
+    // 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.");
+    }
+
     //
     // If the transform is linear, then invert it exactly and return.
@@ -839 +846 @@
     PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(nSamples, 1, 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.");
-    }
-    // XXX: recode or verify after poly changes
-    psS32 order = 1 + in->x->nX;
-
+    //
+    // Allocate a new psPlaneTransform if "out" is NULL, or has the wrong size.
+    //
+    psS32 order = in->x->nX;
     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) {
-        // XXX: Verify this (poly order/nterm change)
         myPT = psPlaneTransformAlloc(order, order);
     } else {
-        // XXX: recode or verify after poly changes
-        if (((1 + out->x->nX) == order) && ((1 + out->x->nY) == order) &&
-                ((1 + out->y->nX) == order) && ((1 + out->y->nX) == order)) {
+        if ((out->x->nX == order) && (out->x->nY == order) &&
+                (out->y->nX == order) && (out->y->nX == order)) {
             myPT = out;
         } else {
@@ -870 +864 @@
 
     //
-    // Copy the input transform to myPT.
-    //
-    for (psS32 i = 0 ; i < (1 + in->x->nX) ; i++) {
-        for (psS32 j = 0 ; j < (1 + in->x->nY) ; j++) {
-            myPT->x->coeff[i][j] = in->x->coeff[i][j];
-        }
-    }
-    for (psS32 i = 0 ; i < (1 + in->y->nX) ; i++) {
-        for (psS32 j = 0 ; j < (1 + 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.  Since we want the inverse of the transformation, the
     // inCoords are written to the outData vector, and the outCoords are written
     // to the inData vector.
     //
+    psVector *xIn = psVectorAlloc(nSamples*nSamples, PS_TYPE_F64);
+    psVector *yIn = psVectorAlloc(nSamples*nSamples, PS_TYPE_F64);
+    psVector *xOut = psVectorAlloc(nSamples*nSamples, PS_TYPE_F64);
+    psVector *yOut = psVectorAlloc(nSamples*nSamples, PS_TYPE_F64);
+    psPlane *inCoord = psPlaneAlloc();
+    psPlane *outCoord = psPlaneAlloc();
     psS32 cnt = 0;
     for (int yint = 0; yint < nSamples; yint++) {
@@ -904 +880 @@
             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;
-
+            xOut->data.F64[cnt] = inCoord->x;
+            yOut->data.F64[cnt] = inCoord->y;
+            xIn->data.F64[cnt] = outCoord->x;
+            yIn->data.F64[cnt] = outCoord->y;
             cnt++;
         }
     }
-    // XXX: what values should be used here?
-    bool rc = psPlaneTransformFit(myPT, inData, outData, 10, 100.0);
+
+    myPT->x = psVectorFitPolynomial2D(myPT->x, NULL, 0, xOut, NULL, xIn, yIn);
+    myPT->y = psVectorFitPolynomial2D(myPT->y, NULL, 0, yOut, NULL, xIn, yIn);
 
     psFree(inCoord);
     psFree(outCoord);
-    psFree(inData);
-    psFree(outData);
-
-    if (rc == true) {
-        return(myPT);
-    } else {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not fit a transform to coordinates.");
+    psFree(xIn);
+    psFree(yIn);
+    psFree(xOut);
+    psFree(yOut);
+
+    if ((myPT->x == NULL) || (myPT->y == NULL)) {
+        psError( PS_ERR_UNKNOWN, true, "psVectorFitPolynomial2D() returned NULL: could not fit a 2-D polynomial to the data.\n");
+        psFree(out);
         return(NULL);
     }
-}
+
+    return(myPT);
+}
+
 
 psPlane *psPlaneTransformDeriv(