Changeset 3375 for trunk/stac/src/stacInvertMaps.c

Timestamp:

Mar 4, 2005, 11:37:01 AM (21 years ago)

Author:

Paul Price

Message:

Updated to psLib rel5alpha1

File:

• trunk/stac/src/stacInvertMaps.c (modified) (6 diffs)

trunk/stac/src/stacInvertMaps.c

@@ -1 +1 @@
 #include <stdio.h>
+#include <assert.h>
 #include "pslib.h"
 #include "stac.h"
@@ -24 +25 @@
         psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[mapNum]; // Uninverted map
         // Check input
-        if (oldMap->x->nX != oldMap->x->nY) {
-            psError("stac.invertMaps", "Polynomial order in x and y don't match for map %d\n", mapNum);
-            return NULL;
-        }
-        if (oldMap->y->nX != oldMap->y->nY) {
-            psError("stac.invertMaps", "Polynomial order in x and y don't match for map %d\n", mapNum);
-            return NULL;
-        }
-        if (oldMap->x->nX != oldMap->y->nX) {
-            psError("stac.invertMaps", "Polynomial order in x and y don't match for map %d\n", mapNum);
-            return NULL;
-        }
+        assert(oldMap->x->nX == oldMap->x->nY && oldMap->y->nX == oldMap->y->nY &&
+               oldMap->x->nX == oldMap->y->nX);
         int order = oldMap->x->nX;      // Polynomial order
         psTrace("stac.invertMaps", 4, "Generating order %d polynomial inverse transformation.\n", order);
@@ -91 +82 @@
 
                     fakePoly->mask[i][j] = 0;
-                    double ijPoly = psDPolynomial2DEval(xIn->data.F32[g], yIn->data.F32[g], fakePoly);
+                    double ijPoly = psDPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
                     fakePoly->mask[i][j] = 1;
 
@@ -98 +89 @@
                             
                             fakePoly->mask[m][n] = 0;
-                            double mnPoly = psDPolynomial2DEval(xIn->data.F32[g], yIn->data.F32[g], fakePoly);
+                            double mnPoly = psDPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
                             fakePoly->mask[m][n] = 1;
 
@@ -113 +104 @@
         // 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
+        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
@@ -178 +169 @@
 #if 0
         // Can't handle higher order than linear yet
-        if (((psPlaneTransform*)maps->data[i])->x->nX != 2 ||
-            ((psPlaneTransform*)maps->data[i])->x->nY != 2 ||
-            ((psPlaneTransform*)maps->data[i])->y->nX != 2 ||
-            ((psPlaneTransform*)maps->data[i])->y->nY != 2) {
-            psError("stac.invertMaps",
-                    "STAC cannot currently support orders other than linear.\n");
-            psFree(inverted);
-            return NULL;
-        }
-
+        assert(((psPlaneTransform*)maps->data[i])->x->nX == 2 &&
+               ((psPlaneTransform*)maps->data[i])->x->nY == 2 &&
+               ((psPlaneTransform*)maps->data[i])->y->nX == 2 &&
+               ((psPlaneTransform*)maps->data[i])->y->nY == 2);
         psPlaneTransform *newMap = psPlaneTransformAlloc(2, 2); // Inverted map
         psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[i]; // Uninverted map