Changeset 10513

Timestamp:

Dec 6, 2006, 2:49:49 PM (20 years ago)

Author:

eugene

Message:

working on non-linear wcs terms

Location:

trunk/psastro/src

Files:

• psastro.h (modified) (1 diff)
• psastroChipAstrom.c (modified) (1 diff)
• psastroWCS.c (modified) (5 diffs)

trunk/psastro/src/psastro.h

@@ -43 +43 @@
 
 bool              pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic);
-bool              pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header, double plateScale);
+bool              pmAstromWriteWCS (psPlaneTransform *toFPA, psPlaneDistort *toTPA, psProjection *toSky, psMetadata *header, double plateScale);
 psPlaneDistort   *psPlaneDistortIdentity ();
+bool psPlaneDistortIsIdentity (psPlaneDistort *distort);
 
 psArray          *psastroLoadRefstars (pmConfig *config);

trunk/psastro/src/psastroChipAstrom.c

@@ -55 +55 @@
                     return false;
                 }
-                pmAstromWriteWCS (chip->toFPA, fpa->toSky, updates, plateScale);
+                pmAstromWriteWCS (chip->toFPA, fpa->toTangentPlane, fpa->toSky, updates, plateScale);
                 psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_DATA_METADATA, "psastro header stats", updates);
                 psFree (updates);

trunk/psastro/src/psastroWCS.c

@@ -208 +208 @@
 // convert toFPA / toSky components to traditional WCS
 // plateScale is nominal physical scale on tangent plane (microns / arcsecond)
-bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header, double plateScale) { 
+// this requires toTP to be the identity transformation
+bool pmAstromWriteWCS (psPlaneTransform *toFPA, psPlaneDistort *toTPA, psProjection *toSky, psMetadata *header, double plateScale) { 
+
+    // techinically, we can have a plate scale here (toTPA:dx,dy != 1)
+    if (!psPlaneDistortIsIdentity (toTPA)) psAbort ("psastro", "invalid TPA transformation");
 
     switch (toSky->type) {
@@ -232 +236 @@
     }
 
-    // XXX not really right: needs to deal with non-identity toTP coeffs
+# if (0)
+    // XXX not really right: needs to deal with non-identity to coeffs
     // XXX actually, totally wrong.  fix the conversions
     // XXX need to handle the plateScale
+    
+    // solve for CDELT1,2:
+    cdelt1 = toSky->Xs*DEG_RAD*toTPA->x->coeff[1][0][0][0];
+    cdelt2 = toSky->Ys*DEG_RAD*toTPA->y->coeff[0][1][0][0];
+
+    // L,M = toFPA(X,Y)
+    // solve for CRPIX1,2 (Xo,Yo) : L,M(Xo,Yo) = 0,0
+    // linear solution for Xo,Yo:
+    R  = (xsum[1][0]*ysum[0][1] - xsum[0][1]*ysum[1][0]);
+    Xo = det*(ysum[0][0]*xsum[0][1] - xsum[0][0]*ysum[0][1]);
+    Yo = det*(xsum[0][0]*ysum[1][0] - ysum[0][0]*xsum[1][0]);
+
+    if (
+# endif
+
     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",   PS_META_REPLACE, "", toSky->R*DEG_RAD);
     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",   PS_META_REPLACE, "", toSky->D*DEG_RAD);
@@ -362 +382 @@
 }
 
-// XXX for the moment, we mimic the limited Elixir mosastro orders
+// construct a psPlaneDistort which is the identify transformation
 psPlaneDistort *psPlaneDistortIdentity (int order) {
 
@@ -376 +396 @@
         for (int j = 0; j <= order; j++) {
             if (i + j > order) {
-                distort->x->mask [1][1][0][0] = 1;
-                distort->y->mask [1][1][0][0] = 1;
+                distort->x->mask [i][j][0][0] = 1;
+                distort->y->mask [i][j][0][0] = 1;
             } 
         }
@@ -386 +406 @@
     return distort;
 }
+
+// check that the given psPlaneDistort is the identity
+bool psPlaneDistortIsIdentity (psPlaneDistort *distort) {
+
+    int order;
+    bool status;
+
+    // we currently only support up to 3rd order polynomials
+    if (distort->x->nX < 1) return false;
+    if (distort->x->nY < 1) return false;
+    if (distort->y->nX < 1) return false;
+    if (distort->y->nY < 1) return false;
+
+    if (distort->x->nX > 3) return false;
+    if (distort->x->nY > 3) return false;
+    if (distort->y->nX > 3) return false;
+    if (distort->y->nY > 3) return false;
+    
+    if (distort->x->nZ > 0) return false;
+    if (distort->x->nT > 0) return false;
+    if (distort->y->nZ > 0) return false;
+    if (distort->y->nT > 0) return false;
+    
+    if (distort->x->nX != distort->x->nY) return false;
+    if (distort->y->nX != distort->y->nY) return false;
+
+    status = true;
+    order = distort->x->nX;
+    for (int i = 0; i <= order; i++) {
+        for (int j = 0; j <= order; j++) {
+            if (i + j > order) {
+                // high-order cross terms must be masked (eg, x^3 y^2)
+                status &= !distort->x->mask[i][j][0][0];
+            } else {
+                if (i + j != 1) {
+                    // non-linear and off-diagonal terms must be 0 (eg, x^2, x y)
+                    status &= (distort->x->coeff[i][j][0][0] == 0.0);
+                } else {
+                    // linear, diagonal terms must be 1.0
+                    status &= (distort->x->coeff[i][j][0][0] == 1.0);
+                }
+            }
+        }
+    }
+
+    order = distort->y->nX;
+    for (int i = 0; i <= order; i++) {
+        for (int j = 0; j <= order; j++) {
+            if (i + j > order) {
+                // high-order cross terms must be masked (eg, x^3 y^2)
+                status &= !distort->y->mask[i][j][0][0];
+            } else {
+                if (i + j != 1) {
+                    // non-linear and off-diagonal terms must be 0 (eg, x^2, x y)
+                    status &= (distort->y->coeff[i][j][0][0] == 0.0);
+                } else {
+                    // linear, diagonal terms must be 1.0
+                    status &= (distort->y->coeff[i][j][0][0] == 1.0);
+                }
+            }
+        }
+    }
+    return status;
+}