Changeset 27838 for branches/tap_branches/psModules/src/astrom/pmAstrometryWCS.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)
• psModules (modified) (1 prop)
• psModules/src/astrom/pmAstrometryWCS.c (modified) (6 diffs)

branches/tap_branches/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;
@@ -866 +865 @@
     }
 
+#define noCOMPARE_TRANS
 #ifdef COMPARE_TRANS
     // compare the computed coordintes from this transform with the original
     psPlane *new = psPlaneAlloc();
+    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);
@@ -875 +876 @@
         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.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" : "");
     }
     psFree(new);
@@ -887 +892 @@
 }
 
-bool pmAstromLinearizeTransforms(pmFPA *fpa, pmChip *chip)
-{
-    psRegion    *chipBounds = pmChipPixels(chip);
-
-    psPlaneTransform *newToFPA = linearFitToTransform(chip->toFPA, chipBounds);
-    if (!newToFPA) {
-        psFree(chipBounds);
-        psError(PS_ERR_UNKNOWN, false, "linear fit for toFPA failed");
-        return false;
-    }
-
-    psFree(chip->toFPA);
-    chip->toFPA = newToFPA;
-
-    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;
-    }
-
-    fpa->toSky->type = PS_PROJ_TAN;
-
-    psFree(chipBounds);
-    psFree(fpaBounds);
+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) {
+        psError(PS_ERR_UNKNOWN, false, "failed to create chipToTPA");
+        return false;
+    }
+
+    // Next do the linear fit within the output boundary pixels
+    psPlaneTransform *chipToFPA = linearFitToTransform(chipToTPA, outputBounds);
+    psFree(chipToTPA);
+    if (!chipToFPA) {
+        psError(PS_ERR_UNKNOWN, false, "linear fit of chip to TPA transform failed");
+        return false;
+    }
+
+    // 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 output chip toFPA failed");
+        return false;
+    }
+
+    // Success. Now set the fpa's toTPA and fromTPA to identity and replace the chip's transforms.
+    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;