Changeset 27838 for branches/tap_branches/pstamp/src/ppstampMakeStamp.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• pstamp/src/ppstampMakeStamp.c (modified) (4 diffs)

branches/tap_branches/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;
         }
@@ -393 +422 @@
     double dx = 0.5 * options->roip.dRA  / fpa->toSky->Xs;
     double dy = 0.5 * options->roip.dDEC / fpa->toSky->Ys;
+
+    if (dx > 5000) {
+        fprintf(stderr, "requested width %f too large reducing to 5000\n", dx);
+        dx = 5000;
+    }
+    if (dy > 5000) {
+        fprintf(stderr, "requested height %f too large reducing to 5000\n", dy);
+        dy = 5000;
+    }
 
     pt->TP->xErr = 0;
@@ -497 +535 @@
             int width  = options->roip.dX;
             int height = options->roip.dY;
+            if (width > 5000) {
+                fprintf(stderr, "requested width %d too large reducing to 5000\n", width);
+                width = 5000;
+            }
+            if (height > 5000) {
+                fprintf(stderr, "requested height %d too large reducing to 5000\n", height);
+                height = 5000;
+            }
 
             // calculate the ROI in chip coordinates