Changeset 2897 for trunk/stac/src/stacRead.c

Timestamp:

Jan 4, 2005, 3:55:40 PM (22 years ago)

Author:

Paul Price

Message:

Updated to handle higher polynomial orders than linear, and to calculate the output size itself

File:

• trunk/stac/src/stacRead.c (modified) (2 diffs)

trunk/stac/src/stacRead.c

@@ -3 +3 @@
 #include "pslib.h"
 #include "stac.h"
-
-#define ORDER 2
 
 psArray *stacReadImages(stacConfig *config)
@@ -54 +52 @@
         // Read the file
         psTrace("stac.read.maps",5,"Reading map file %s....\n",mapfile);
-        psPlaneTransform *map = psPlaneTransformAlloc(ORDER,ORDER);
 
-        // Format is as following:
-        // A0  A1
-        // B00 B01
-        // B10 B11
-        //
-        // where:
-        // x_sky = A0 + B00*x + B10*y
-        // y_sky = A1 + B01*x + B11*y
+        // Format is now:
+        // order
+        // x coefficients
+        // y coefficients
+        // 
+        // where the order is 1 for linear, 2 for quadratic, 3 for cubic.
+        // and the coefficients are read by the following pseudo-code:
+        // 
+        //      for (int k = 0; k < order + 1; k++)
+        //          for (int j = 0; j < k; j++)
+        //              int i = k - j;
+        //              read coefficient of x^i y^j
+        // 
+        // This produces the ordering:
+        // x^0 y^0, x^1 y^0, x^0 y^1, x^2 y^0, x^1 y^1, x^0 y^2, x^3 y^0, x^2 y^1, x^1 y^2, x^0 y^3
+        // 
+        // This is, of course, for third order polynomials.
+        // For lower orders, the list is truncated at the appropriate level.
 
-        double a0, a1, b00, b10, b01, b11;
-
-        if (fscanf(mapfp,"%lf %lf %lf %lf %lf %lf", &a0, &a1, &b00, &b01, &b10, &b11) != 6) {
-            psLogMsg("stac.read.maps", PS_LOG_ERROR, "Couldn't read all 6 transformation values from %s\n",
-                     mapfile);
+        int order = 0;                  // Polynomial order
+        if (fscanf(mapfp, "%d", &order) != 1) {
+            psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
             exit(EXIT_FAILURE);
         }
-        fclose(mapfp);
-        map->x->coeff[0][0] = a0;
-        map->x->coeff[1][0] = b00;
-        map->x->coeff[0][1] = b10;
-        map->y->coeff[0][0] = a1;
-        map->y->coeff[1][0] = b01;
-        map->y->coeff[0][1] = b11;
 
-        // Mask out the cross-terms
-        map->x->mask[1][1] = 1;
-        map->y->mask[1][1] = 1;
-        // Set the cross-terms to zero, just in case...
-        map->x->coeff[1][1] = 0.0;
-        map->y->coeff[1][1] = 0.0;
+        psTrace("stac.read.maps", 5, "Polynomial order: %d\n", order);
+
+        psPlaneTransform *map = psPlaneTransformAlloc(order + 1, order + 1); // The transformation
+        // Set coefficient masks
+        for (int i = 0; i < order + 1; i++) {
+            for (int j = 0; j < order + 1; j++) {
+                if (i + j > order + 1) {
+                    map->x->mask[i][j] = 1;
+                    map->y->mask[i][j] = 1;
+                } else {
+                    map->x->mask[i][j] = 0;
+                    map->y->mask[i][j] = 0;
+                }
+            }
+        }
+
+        // Read x coefficients
+        psTrace("stac.read.maps", 7, "x' = \n");
+        for (int k = 0; k < order + 1; k++) {
+            for (int j = 0; j < k + 1; j++) {
+                int i = k - j;
+                if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
+                    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
+                    exit(EXIT_FAILURE);
+                }
+                psTrace("stac.read.maps", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
+            }
+        }
+        // Read y coefficients
+        psTrace("stac.read.maps", 7, "y' = \n");
+        for (int k = 0; k < order + 1; k++) {
+            for (int j = 0; j < k + 1; j++) {
+                int i = k - j;
+                if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
+                    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
+                    exit(EXIT_FAILURE);
+                }
+                psTrace("stac.read.maps", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
+            }
+        }
+
+        fclose(mapfp);
 
         // Plug it into the array
         maps->data[i] = map;
-
-#ifdef TESTING
-        // Check psLib's answer is the same as I would get with old version
-        double detx = 256.0, dety = 256.0;
-        // My code
-        double myx = a0 + b00*detx + b10*dety;
-        double myy = a1 + b01*detx + b11*dety;
-        // psLib
-        psPlane *det = psAlloc(sizeof(psPlane));
-        det->x = 256.0;
-        det->y = 256.0;
-        psPlane *pslib = psPlaneTransformApply(NULL, map, det);
-
-        psTrace("stac.read.maps.test",0,"me: %lf,%lf\n",myx,myy);
-        psTrace("stac.read.maps.test",0,"ps: %lf,%lf\n",pslib->x,pslib->y);
-        psFree(det);
-        psFree(pslib);
-#endif
 
     }