Changeset 26492

Timestamp:

Dec 29, 2009, 11:06:35 AM (16 years ago)

Author:

bills

Message:

combine toFPA and toTPA transforms and do a single linear fit

File:

• trunk/psModules/src/astrom/pmAstrometryWCS.c (modified) (7 diffs)

trunk/psModules/src/astrom/pmAstrometryWCS.c

@@ -378 +378 @@
     // XXX make it optional to write out CDi_j terms, or other versions
     // apply CDELT1,2 (degrees / pixel) to yield PCi,j terms of order unity
-    if (!(wcs->wcsCDkeys)) {
+    if (!wcs->wcsCDkeys) { 
 
       double cdelt1 = wcs->cdelt1;
@@ -419 +419 @@
         psMetadataRemoveKey(header, "CD2_2");
       }
-    }
-    if (wcs->wcsCDkeys) {
+    } else {
 
       psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_1", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0]);
@@ -845 +844 @@
     int k=0;
     for (int j=0; j<nSamples; j++) {
-        double y = j * deltaY / nSamples;
+        double y = bounds->y0 + (j * deltaY / nSamples);
         for (int i=0; i<nSamples; i++) {
             psPlane *s = psPlaneAlloc();
-            s->x = i * deltaX / nSamples;
+            s->x = bounds->x0 + (i * deltaX / nSamples);
             s->y = y;
             psArraySet(src, k, s);
             psPlane *d = psPlaneTransformApply(NULL, trans, s);
             psArraySet(dst, k, d);
-            psFree(s);
+            psFree(s);  // drop our refs to s and d
             psFree(d);
             ++k;
@@ -869 +868 @@
     // compare the computed coordintes from this transform with the original
     psPlane *new = psPlaneAlloc();
+    printf("   i     chip_x  fpa_x     fpa_x_fit     dx         chip_y    fpa_y     fpa_y_fit     dy     dx > 0.5 || dy > 0.5\n");
     for (int i=0; i<psArrayLength(dst); i++) {
         psPlane *d = (psPlane *) psArrayGet(dst, i);
@@ -875 +875 @@
         new = psPlaneTransformApply(new, newTrans, s);
 
-        printf("%4d %f %f\n", i, 100.*(new->x - d->x)/d->x, 100.*(new->y - d->y)/d->y);
+        double xerr = new->x - d->x;
+        double yerr = new->y - d->y;
+        bool bigerr = (fabs(xerr) > .5) || (fabs(yerr) > .5);
+        printf("%4d %9.2f %9.2f %9.2f %9.2f     %9.2f %9.2f %9.2f %9.2f   %s\n"
+        , i, s->x, new->x, d->x, xerr, s->y, new->y, d->y, yerr, bigerr ? "BIGERR" : "");
     }
     psFree(new);
@@ -891 +895 @@
     psRegion    *chipBounds = pmChipPixels(chip);
 
-    psPlaneTransform *newToFPA = linearFitToTransform(chip->toFPA, chipBounds);
+#ifdef TESTING_CAMRUN_14728
+    chipBounds->y0 = 1874.;
+    chipBounds->x1 = 2387.;
+#endif
+
+    // First combine the chip to FPA and FPA to TPA into a single transformation
+    psPlaneTransform *chipToTPA = psPlaneTransformCombine(NULL, chip->toFPA, fpa->toTPA, *chipBounds, 50);
+    if (!chipToTPA) {
+        psFree(chipBounds);
+        psError(PS_ERR_UNKNOWN, false, "failed to create chipToTPA");
+        return false;
+    }
+
+    // Next do a linear fit
+    psPlaneTransform *newToFPA = linearFitToTransform(chipToTPA, chipBounds);
+    psFree(chipToTPA);
     if (!newToFPA) {
         psFree(chipBounds);
-        psError(PS_ERR_UNKNOWN, false, "linear fit for toFPA failed");
-        return false;
-    }
+        psError(PS_ERR_UNKNOWN, false, "linear fit of chip to TPA transform failed");
+        return false;
+    }
+
+    psPlaneTransform *newFromFPA = psPlaneTransformInvert(NULL, newToFPA, *chipBounds, 50);
+    psFree(chipBounds);
+    if (!newFromFPA) {
+        psFree(newToFPA);
+        psError(PS_ERR_UNKNOWN, false, "inversion of fit of chip to TPA transform failed");
+        return false;
+    }
+
+    // Success. Now set the fpa's toTPA and fromTPA to identity and replace the chip's transforms.
+
+    psPlaneTransform *newToTPA   = psPlaneTransformIdentity(1);
+    psFree(fpa->toTPA);
+    fpa->toTPA = newToTPA;
+
+    psPlaneTransform *newFromTPA = psPlaneTransformIdentity(1);
+    psFree(fpa->fromTPA);
+    fpa->fromTPA = newFromTPA;
 
     psFree(chip->toFPA);
@@ -902 +939 @@
 
     psFree(chip->fromFPA);
-    chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *chipBounds, 50);
-    if (!chip->fromFPA) {
-        psError(PS_ERR_UNKNOWN, false, "failed to invert linear fit for toFPA");
-        return false;
-    }
-
-    psPlane *chip0 = psPlaneAlloc();
-    chip0->x = 0;
-    chip0->y = 0;
-    psPlane *chip1 = psPlaneAlloc();
-    chip1->x = chipBounds->x1;
-    chip1->y = chipBounds->y1;
-
-    // compute bounding region for fpa
-    psPlane *fpa0 = psPlaneTransformApply(NULL, newToFPA, chip0);
-    psPlane *fpa1 = psPlaneTransformApply(NULL, newToFPA, chip1);
-
-    psRegion *fpaBounds = psRegionAlloc(fpa0->x, fpa1->x, fpa0->y, fpa1->y);
-    psFree(chip0);
-    psFree(chip1);
-    psFree(fpa0);
-    psFree(fpa1);
-
-    psPlaneTransform *newToTPA = linearFitToTransform(fpa->toTPA, fpaBounds);
-    if (!newToTPA) {
-        psError(PS_ERR_UNKNOWN, false, "failed to perform linear fit to toTPA");
-        psFree(fpaBounds);
-        return false;
-    }
-    psFree(fpa->toTPA);
-    fpa->toTPA = newToTPA;
-
-    // XXX: is this region ok?
-    psFree(fpa->fromTPA);
-    fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, *fpaBounds, 50);
-    if (!fpa->fromTPA) {
-        psError(PS_ERR_UNKNOWN, false, "failed to invert linear fit to toTPA");
-        return false;
-    }
-
+    chip->fromFPA = newFromFPA;
+
+    // Finally change the type for the projection.
     fpa->toSky->type = PS_PROJ_TAN;
-
-    psFree(chipBounds);
-    psFree(fpaBounds);
 
     return true;