Changeset 27461

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

Location:

trunk

Files:

• psModules/src/astrom/pmAstrometryWCS.c (modified) (3 diffs)
• psModules/src/astrom/pmAstrometryWCS.h (modified) (1 diff)
• pstamp/src/ppstampMakeStamp.c (modified) (2 diffs)

trunk/psModules/src/astrom/pmAstrometryWCS.c

@@ -865 +865 @@
     }
 
+#define noCOMPARE_TRANS
 #ifdef COMPARE_TRANS
     // 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");
+    printf("   i     chip_x  tpa_x     tpa_x_fit     dx         chip_y    tpa_y     tpa_y_fit     dy     dx > 0.5 || dy > 0.5\n");
     for (int i=0; i<psArrayLength(dst); i++) {
         psPlane *d = (psPlane *) psArrayGet(dst, i);
@@ -878 +879 @@
         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"
+        printf("%4d %9.2f %9.2f %9.2f %9.4f     %9.2f %9.2f %9.2f %9.4f   %s\n"
         , i, s->x, new->x, d->x, xerr, s->y, new->y, d->y, yerr, bigerr ? "BIGERR" : "");
     }
@@ -891 +892 @@
 }
 
-bool pmAstromLinearizeTransforms(pmFPA *fpa, pmChip *chip)
-{
-    psRegion    *chipBounds = pmChipPixels(chip);
-
-#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);
+bool pmAstromLinearizeTransforms(pmFPA *inFPA, pmChip *inChip, pmFPA *outFPA, pmChip *outChip, psRegion *outputBounds, double offset_x, double offset_y)
+{
+    PS_ASSERT_PTR_NON_NULL(inFPA, NULL);
+    PS_ASSERT_PTR_NON_NULL(inChip, NULL);
+
+    if (outFPA == NULL) {
+        outFPA = inFPA;
+    }
+    if (outChip == NULL) {
+        outChip = inChip;
+    }
+    if (outputBounds == NULL) {
+        outputBounds = pmChipPixels(outChip);
+    }
+
+    // First combine the "chip to FPA" and "FPA to TPA" into a single transformation
+    psPlaneTransform *chipToTPA = psPlaneTransformCombine(NULL, inChip->toFPA, inFPA->toTPA, *outputBounds, 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);
+    // Next do the linear fit within the output boundary pixels
+    psPlaneTransform *chipToFPA = linearFitToTransform(chipToTPA, outputBounds);
     psFree(chipToTPA);
-    if (!newToFPA) {
-        psFree(chipBounds);
+    if (!chipToFPA) {
         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 requested,  change the center
+    psPlaneTransform *outToFPA;
+    if (offset_x != 0. && offset_y != 0.) {
+        outToFPA = psPlaneTransformSetCenter(NULL, chipToFPA, offset_x, offset_y);
+        psFree(chipToFPA);
+    } else {
+        outToFPA = chipToFPA;
+    }
+
+    psPlaneTransform *outFromFPA = psPlaneTransformInvert(NULL, outToFPA, *outputBounds, 50);
+    if (!outFromFPA) {
+        psFree(outToFPA);
+        psError(PS_ERR_UNKNOWN, false, "inversion of fit of output chip toFPA failed");
+        return false;
+    }
+
+    // Success. Now set the fpa's toTPA and fromTPA to identity and replace the chip's transforms.
+
+    psFree(outFPA->toTPA);
+    outFPA->toTPA =  psPlaneTransformIdentity(1);
+
+    psFree(outFPA->fromTPA);
+    outFPA->fromTPA = psPlaneTransformIdentity(1);
+
+    psFree(outChip->toFPA);
+    outChip->toFPA = outToFPA;
+
+    psFree(outChip->fromFPA);
+    outChip->fromFPA = outFromFPA;
+
+    // Finally, change the type for the projection.
+    outFPA->toSky->type = PS_PROJ_TAN;
+
+    return true;
+}
+
+bool pmAstromLinearizeToSky(pmFPA *inFPA, pmChip *inChip, pmFPA *outFPA, pmChip *outChip, psRegion *bounds)
+{
+    PS_ASSERT_PTR_NON_NULL(inFPA, NULL);
+    PS_ASSERT_PTR_NON_NULL(inChip, NULL);
+    PS_ASSERT_PTR_NON_NULL(outFPA, NULL);
+    PS_ASSERT_PTR_NON_NULL(outChip, NULL);
+    PS_ASSERT_PTR_NON_NULL(bounds, NULL);
+
+    // outFPA projection must be defined as the goal
+    
+    // the output transformations are:
+    // chip -> FPA : standard linear trans with needed rotation, etc
+    // FPA  -> TPA : identidy
+
+    int nSamples = 10;  // 10 samples in each dimension
+
+    double deltaX = (bounds->x1 - bounds->x0);
+    double deltaY = (bounds->y1 - bounds->y0);
+
+    psArray *src = psArrayAllocEmpty(nSamples * nSamples);
+    psArray *dst = psArrayAllocEmpty(nSamples * nSamples);
+
+    psPlane srcFP, srcTP;
+
+    for (int j = 0; j < nSamples; j++) {
+        double y = bounds->y0 + (j * deltaY / nSamples);
+        for (int i =  0; i < nSamples; i++) {
+
+            psSphere srcSky;
+            psPlane *srcChip = psPlaneAlloc();
+            psPlane *dstTP = psPlaneAlloc();
+
+            srcChip->x = bounds->x0 + (i * deltaX / nSamples);
+            srcChip->y = y;
+
+            psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip);
+            psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP);
+            psDeproject (&srcSky, &srcTP, inFPA->toSky);
+            
+            // fprintf (stderr, "%f %f | %f %f | %f %f | %f %f\n", srcChip->x, srcChip->y, srcFP.x, srcFP.y, srcTP.x, srcTP.y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD);
+
+            psProject (dstTP, &srcSky, outFPA->toSky);
+
+            srcChip->x -= bounds->x0;
+            srcChip->y -= bounds->y0;
+            psArrayAdd (src, 100, srcChip);
+            psArrayAdd (dst, 100, dstTP);
+
+            psFree(srcChip);  // drop our refs to s and d
+            psFree(dstTP);
+        }
+    }
+
+    psPlaneTransform *newToFPA = psPlaneTransformAlloc(1, 1);
+    newToFPA->x->coeffMask[1][1] = 1;
+    newToFPA->y->coeffMask[1][1] = 1;
+
+    if (!psPlaneTransformFit(newToFPA, src, dst, 0, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed");
+        psFree(src);
+        psFree(dst);
+        return NULL;
+    }
+    
+# if (0)
+    for (int i = 0; i < src->n; i++) {
+        
+        psSphere srcSky, dstSky;
+        psPlane *srcChip = src->data[i];
+        psPlane *dstTP   = dst->data[i];
+
+        psPlaneTransformApply (&srcFP, newToFPA, srcChip);
+        psDeproject (&srcSky, &srcFP, outFPA->toSky);
+        psDeproject (&dstSky, dstTP, outFPA->toSky);
+
+        double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0;
+        double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0;
+        fprintf (stderr, "%f %f | %f %f | %f %f | %f %f | %f %f | %f %f\n", dX, dY, srcChip->x, srcChip->y, srcFP.x, srcFP.y, dstTP->x, dstTP->y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD, dstSky.r*PS_DEG_RAD, dstSky.d*PS_DEG_RAD);
+
+    }
+# endif
+
+    psFree(src);
+    psFree(dst);
+
+    // this is a linear transformation
+    psPlaneTransform *newFromFPA = psPlaneTransformInvert(NULL, newToFPA, *bounds, 1);
     if (!newFromFPA) {
         psFree(newToFPA);
-        psError(PS_ERR_UNKNOWN, false, "inversion of fit of chip to TPA transform failed");
+        psError(PS_ERR_UNKNOWN, false, "inversion of fit of output chip toFPA 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);
-    chip->toFPA = newToFPA;
-
-    psFree(chip->fromFPA);
-    chip->fromFPA = newFromFPA;
-
-    // Finally change the type for the projection.
-    fpa->toSky->type = PS_PROJ_TAN;
+    psFree(outChip->toFPA);
+    outChip->toFPA = newToFPA;
+
+    psFree(outChip->fromFPA);
+    outChip->fromFPA = newFromFPA;
+
+    psFree(outFPA->toTPA);
+    outFPA->toTPA =  psPlaneTransformIdentity(1);
+
+    psFree(outFPA->fromTPA);
+    outFPA->fromTPA = psPlaneTransformIdentity(1);
 
     return true;

trunk/psModules/src/astrom/pmAstrometryWCS.h

@@ -63 +63 @@
 bool pmAstromWriteBilevelMosaic (psMetadata *header, const pmFPA *fpa, double tol);
 
-bool pmAstromLinearizeTransforms(pmFPA *fpa, pmChip *);
+bool pmAstromLinearizeTransforms(pmFPA *inFPA, pmChip *inChip, pmFPA *outFPA, pmChip *outChip, psRegion *outputBounds, double offset_x, double offset_y);
+bool pmAstromLinearizeToSky(pmFPA *inFPA, pmChip *inChip, pmFPA *outFPA, pmChip *outChip, psRegion *bounds);
 
 // move to pslib

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;
         }