Changeset 27840 for branches/simtest_nebulous_branches/psModules/src/astrom/pmAstrometryWCS.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/astrom/pmAstrometryWCS.c (modified) (13 diffs)

branches/simtest_nebulous_branches/psModules/src/astrom/pmAstrometryWCS.c

@@ -289 +289 @@
     // test the CDELTi varient
     if (pcKeys) {
+        wcs->wcsCDkeys = 0;
         wcs->cdelt1 = psMetadataLookupF64 (&status, header, "CDELT1");
         wcs->cdelt2 = psMetadataLookupF64 (&status, header, "CDELT2");
@@ -334 +335 @@
     // test the CDi_j varient
     if (cdKeys) {
+        wcs->wcsCDkeys = 1;
+
         wcs->trans->x->coeff[1][0] = psMetadataLookupF64 (&status, header, "CD1_1"); // == PC1_1
         wcs->trans->x->coeff[0][1] = psMetadataLookupF64 (&status, header, "CD1_2"); // == PC1_2
@@ -375 +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
-    double cdelt1 = wcs->cdelt1;
-    double cdelt2 = wcs->cdelt2;
-    psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1);
-    psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2);
-
-    // test the PC00i00j varient:
-    psMetadataAddF64 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1
-    psMetadataAddF64 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1] / cdelt2); // == PC1_2
-    psMetadataAddF64 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1
-    psMetadataAddF64 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2
-
-    // Elixir-style polynomial terms
-    // XXX currently, Elixir/DVO cannot accept mixed orders
-    // XXX need to respect the masks
-    // XXX is wcs->cdelt1,2 always consistent?
-    int fitOrder = wcs->trans->x->nX;
-    if (fitOrder > 1) {
+    if (!wcs->wcsCDkeys) { 
+
+      double cdelt1 = wcs->cdelt1;
+      double cdelt2 = wcs->cdelt2;
+      psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1);
+      psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2);
+      
+      // test the PC00i00j varient:
+      psMetadataAddF64 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1
+      psMetadataAddF64 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1] / cdelt2); // == PC1_2
+      psMetadataAddF64 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1
+      psMetadataAddF64 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2
+      
+      // Elixir-style polynomial terms
+      // XXX currently, Elixir/DVO cannot accept mixed orders
+      // XXX need to respect the masks
+      // XXX is wcs->cdelt1,2 always consistent?
+      int fitOrder = wcs->trans->x->nX;
+      if (fitOrder > 1) {
         for (int i = 0; i <= fitOrder; i++) {
-            for (int j = 0; j <= fitOrder; j++) {
-                if (i + j < 2)
-                    continue;
-                if (i + j > fitOrder)
-                    continue;
-                sprintf (name, "PCA1X%1dY%1d", i, j);
-                psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
-                sprintf (name, "PCA2X%1dY%1d", i, j);
-                psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
-            }
+          for (int j = 0; j <= fitOrder; j++) {
+            if (i + j < 2)
+              continue;
+            if (i + j > fitOrder)
+              continue;
+            sprintf (name, "PCA1X%1dY%1d", i, j);
+            psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
+            sprintf (name, "PCA2X%1dY%1d", i, j);
+            psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
+          }
         }
         psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", fitOrder);
-    }
-
-    // remove any existing 'CDi_j style' wcs keywords
-    if (psMetadataLookup(header, "CD1_1")) {
+      }
+      
+      // remove any existing 'CDi_j style' wcs keywords
+      if (psMetadataLookup(header, "CD1_1")) {
         psMetadataRemoveKey(header, "CD1_1");
         psMetadataRemoveKey(header, "CD1_2");
         psMetadataRemoveKey(header, "CD2_1");
         psMetadataRemoveKey(header, "CD2_2");
+      }
+    } else {
+
+      psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_1", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0]);
+      psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_2", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1]);
+      psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_1", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0]);
+      psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_2", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1]);
+
+      if (psMetadataLookup(header, "PC001001")) {
+        psMetadataRemoveKey(header, "PC001001");
+        psMetadataRemoveKey(header, "PC001002");
+        psMetadataRemoveKey(header, "PC002001");
+        psMetadataRemoveKey(header, "PC002002");
+      }
     }
 
@@ -535 +554 @@
         fpa->toSky->R -= 2.0*M_PI;
 
+    fpa->wcsCDkeys = wcs->wcsCDkeys;
+
     psTrace ("psastro", 5, "toFPA: %f %f  (%f,%f),(%f,%f)\n",
              chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0],
@@ -648 +669 @@
     wcs->crval2 = fpa->toSky->D*PS_DEG_RAD;
 
+    // generate a transform that has 0.0 rotation:
+    // get the current posangle of the ref chip
+    // XXX average angles for x and y...
+    float angle = atan2 (toTPA->y->coeff[1][0], toTPA->x->coeff[1][0]);
+    // fprintf (stderr, "angle: %f\n", angle*PS_DEG_RAD);
+    psPlaneTransform *tpa1 = psPlaneTransformRotate (NULL, toTPA, angle);
+
     // given transformation, solve for coordinates which yields output coordinates of 0,0
-    psPlane *center = psPlaneTransformGetCenter (toTPA, tol);
+    psPlane *center = psPlaneTransformGetCenter (tpa1, tol);
     if (!center) {
         psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
@@ -657 +685 @@
     }
 
+    // generate transform with the original orientation (does this rotate about 'center'?) 
+    psPlaneTransform *tpa2 = psPlaneTransformRotate (NULL, tpa1, -1.0*angle);
+
+    // prove that the center coordinates give 0,0:
+    // float Xo = psPolynomial2DEval (tpa1->x, center->x, center->y);
+    // float Yo = psPolynomial2DEval (tpa1->x, center->x, center->y);
+    // fprintf (stderr, "tpa1: Xo, Yo: %f, %f\n", Xo, Yo);
+
+    // prove that the center coordinates give 0,0:
+    // Xo = psPolynomial2DEval (tpa2->x, center->x, center->y);
+    // Yo = psPolynomial2DEval (tpa2->x, center->x, center->y);
+    // fprintf (stderr, "tpa2: Xo, Yo: %f, %f\n", Xo, Yo);
+
     // create wcs transform from toFPA, resulting transformation has units of microns/pixel
     // adjust wcs transform to use center as reference coordinate
-    psPlaneTransformSetCenter (wcs->trans, toTPA, center->x, center->y);
+    psPlaneTransformSetCenter (wcs->trans, tpa2, center->x, center->y);
 
     // calculated center is crpix1,2
@@ -665 +706 @@
     wcs->crpix2 = center->y;
     psFree (center);
+    psFree (tpa1);
+    psFree (tpa2);
 
     // pdelt1,2 has units of degrees/micron
@@ -682 +725 @@
     wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
 
+    wcs->wcsCDkeys = fpa->wcsCDkeys;
     psFree (toTPA);
 
@@ -800 +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;
@@ -821 +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);
@@ -830 +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);
@@ -842 +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;
@@ -922 +1084 @@
     wcs->trans = psPlaneTransformAlloc (nXorder, nYorder);
     wcs->toSky = NULL;
+    wcs->wcsCDkeys = 0;
 
     memset (wcs->ctype1, 0, PM_ASTROM_WCS_TYPE_SIZE);