Changeset 27461 for trunk/pstamp/src

Timestamp:

Mar 26, 2010, 11:21:55 AM (16 years ago)

Author:

bills

Message:

Improve fit to bi-level astrometry for postage stamps by making astrometry terms that are centered on the output chip center

File:

• trunk/pstamp/src/ppstampMakeStamp.c (modified) (2 diffs)

trunk/pstamp/src/ppstampMakeStamp.c

@@ -17 +17 @@
 
 static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip);
+static void chipToSky(pmAstromObj *pt, pmFPA *fpa, pmChip *chip);
 static bool setMaskedToNAN(pmConfig *config, psImage *image, psImage *mask, psImage *variance);
 
 // convert the input chip's transforms to the output
-static bool convertWCS(pmHDU *outHDU, pmFPA *outFPA, pmChip *outChip, pmChip *inChip, psRegion *roi)
-{
-    pmHDU *hdu = pmHDUFromChip(inChip);
-    PS_ASSERT_PTR_NON_NULL(hdu, 1)
-    PS_ASSERT_PTR_NON_NULL(hdu->header, 1)
-
-    // Change the reference pixel to account for the change in origin between the stamp and
-    // the original image
-    // XXX: shouldn't we be using the center of the stamp instead of the corner?
-    outChip->toFPA   = psPlaneTransformSetCenter(NULL, inChip->toFPA, (int) roi->x0, (int) roi->y0);
-
-    outChip->fromFPA = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50);
+static bool convertWCS(pmHDU *outHDU, pmFPA *inFPA, pmChip *inChip, pmFPA *outFPA, pmChip *outChip, psRegion *roi)
+{
+    PS_ASSERT_PTR_NON_NULL(outHDU, 1);
+    PS_ASSERT_PTR_NON_NULL(outHDU->header, 1);
+
+    // Center of stamp in chip coordinates
+    double center_x = (roi->x0 + roi->x1) / 2;
+    double center_y = (roi->y0 + roi->y1) / 2;
+    pmAstromObj *check = pmAstromObjAlloc();
+    check->chip->x = center_x;
+    check->chip->y = center_y;
+    chipToSky(check, inFPA, inChip);
+    double r_before = check->sky->r;
+    double d_before = check->sky->d;
+    psLogMsg("ppstampMakeStamp", 2, "Before fit  Center Pixel RA: %f   DEC: %f (degrees)\n",
+             r_before * PS_DEG_RAD, d_before * PS_DEG_RAD);
 
     if (outFPA->toSky->type != PS_PROJ_TAN) {
-        if (!pmAstromLinearizeTransforms(outFPA, outChip)) {
+        // we need to make astrometry terms for the output which are centered on the output chip center,
+        // but we keep the original plate scale
+        psFree(outFPA->toSky);
+        outFPA->toSky = psProjectionAlloc (r_before, d_before, inFPA->toSky->Xs, inFPA->toSky->Ys, PS_PROJ_TAN);
+
+        if (!pmAstromLinearizeToSky(inFPA, inChip, outFPA, outChip, roi)) {
+            psFree(check);
             psError(PS_ERR_UNKNOWN, false, "Failed to linearize astrometry\n");
             return false;
         }
-    }
+        check->chip->x = (roi->x1 - roi->x0) / 2;
+        check->chip->y = (roi->y1 - roi->y0) / 2;
+        chipToSky(check, outFPA, outChip);
+        double r_after = check->sky->r;
+        double d_after = check->sky->d;
+
+        psLogMsg("ppstampMakeStamp", 2, "After fit:  Center Pixel RA: %f   DEC: %f (degrees)\n",
+                 r_after * PS_DEG_RAD, d_after * PS_DEG_RAD);
+        psLogMsg("ppstampMakeStamp", 2, "Error in fit to astrometry %.2f    %.2f (arcseconds)\n", 
+                 (r_after - r_before) *  PS_DEG_RAD * 3600, (d_after - d_before) * PS_DEG_RAD * 3600);
+
+        // XXX: should we fail if the fit is bad ??
+
+    } else {
+        outChip->toFPA     = psPlaneTransformSetCenter(NULL, inChip->toFPA, (int) roi->x0, (int) roi->y0);
+        outChip->fromFPA   = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50);
+        if (!outChip->fromFPA) {
+            psError(PS_ERR_UNKNOWN, false, "inversion of toFPA transform failed");
+            return false;
+        }
+    }
+    psFree(check);
 
     if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_NONLIN_TOL)) {
@@ -92 +124 @@
     // If input had WCS convert it for stamp
     if (input->fpa->toSky) {
-        // steal the input fpa's transforms
-        output->fpa->toTPA   = input->fpa->toTPA;
-        output->fpa->fromTPA = input->fpa->fromTPA;
-        output->fpa->toSky   = input->fpa->toSky;
-
-        // drop the references
-        input->fpa->toTPA   = 0;
-        input->fpa->fromTPA = 0;
-        input->fpa->toSky   = 0;
-
-        if (!convertWCS(outHDU, output->fpa, outChip, inChip, &options->roi)) {
+        // copy the input fpa's transforms
+        // XXX: if we change to making multiple stamps per process we'll need to copy
+        // the transformations
+        output->fpa->toTPA   = psMemIncrRefCounter(input->fpa->toTPA);
+        output->fpa->fromTPA = psMemIncrRefCounter(input->fpa->fromTPA);
+        output->fpa->toSky   = psMemIncrRefCounter(input->fpa->toSky);
+
+        if (!convertWCS(outHDU, input->fpa, inChip, output->fpa, outChip, &options->roi)) {
             return false;
         }