Changeset 15896

Timestamp:

Dec 22, 2007, 7:52:57 AM (18 years ago)

Author:

eugene

Message:

apply toTPA to toFPA when writing out single chip version

File:

• trunk/psModules/src/astrom/pmAstrometryWCS.c (modified) (9 diffs)

trunk/psModules/src/astrom/pmAstrometryWCS.c

@@ -7 +7 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-10-09 19:27:04 $
+ *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-12-22 17:52:57 $
  *
  *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -56 +56 @@
 bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip, double tol)
 {
-    // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY
-    // XXX require chip->toFPA->y->nX == chip->toFPA->y->nY
-    // XXX require chip->toFPA->x->nX == chip->toFPA->y->nX
-    // XXX require chip->toFPA->nX == 1,2,3
-
-    // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
-    if (!psPlaneTransformIsDiagonal (fpa->toTPA))
-        psAbort("invalid TPA transformation");
-
     pmAstromWCS *wcs = pmAstromWCSfromFPA(fpa, chip, tol);
+    if (!wcs) return false;
+
     pmAstromWCStoHeader (header, wcs);
 
@@ -118 +111 @@
 {
     pmAstromWCS *wcs = pmAstromWCSBilevelChipFromFPA (chip, tol);
-
     pmAstromWCStoHeader (header, wcs);
 
@@ -337 +329 @@
     char *type;
 
+    if (!wcs) return false;
+
     type = psProjectTypeToString (wcs->toSky->type, "RA--");
     psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", type);
@@ -385 +379 @@
     }
 
-    return (true);
+    return true;
 }
 
@@ -550 +544 @@
 pmAstromWCS *pmAstromWCSfromFPA (const pmFPA *fpa, const pmChip *chip, double tol)
 {
-
-    // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
-    // XXX not really: toTPA needs to have unity scale for distortion fitting function
-    if (!psPlaneTransformIsDiagonal (fpa->toTPA))
-        psAbort("invalid TPA transformation");
-
     // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY
     // XXX require chip->toFPA->y->nX == chip->toFPA->y->nY
@@ -561 +549 @@
     // XXX require chip->toFPA->nX == 1,2,3
 
-    pmAstromWCS *wcs = pmAstromWCSAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY);
+    // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
+    // XXX not really: toTPA needs to have unity scale for distortion fitting function
+    // if (!psPlaneTransformIsDiagonal (fpa->toTPA))
+    // psAbort("invalid TPA transformation");
+
+    // create a temporary transform which combines toTPA and toFPA.  allow toTPA to have 0th
+    // and 1st order terms
+    
+    // XXX require fpa->toTPA->x->nX == fpa->toTPA->x->nY
+    // XXX require fpa->toTPA->y->nX == fpa->toTPA->y->nY
+    // XXX require fpa->toTPA->x->nX == fpa->toTPA->y->nX
+    // XXX require fpa->toTPA->x->coeffMask[1][1]
+    // XXX require fpa->toTPA->y->coeffMask[1][1]
+    // XXX require fpa->toTPA->nX == 1
+
+    psPlaneTransform *toTPA = psPlaneTransformAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY);
+
+    for (int i = 0; i <= toTPA->x->nX; i++) {
+        for (int j = 0; j <= toTPA->x->nY; j++) {
+            double f1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[1][0]*chip->toFPA->x->coeff[i][j];
+            double f2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[0][1]*chip->toFPA->y->coeff[i][j];
+            toTPA->x->coeff[i][j] = f1 + f2;
+
+            double g1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[1][0]*chip->toFPA->x->coeff[i][j];
+            double g2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[0][1]*chip->toFPA->y->coeff[i][j];
+            toTPA->y->coeff[i][j] = g1 + g2;
+        }
+    }
+    toTPA->x->coeff[0][0] += fpa->toTPA->x->coeff[0][0];
+    toTPA->y->coeff[0][0] += fpa->toTPA->y->coeff[0][0];
+
+    pmAstromWCS *wcs = pmAstromWCSAlloc(toTPA->x->nX, toTPA->x->nY);
 
     // convert projection from FPA to SKY into wcs projection (degrees to radians)
@@ -569 +588 @@
 
     // given transformation, solve for coordinates which yields output coordinates of 0,0
-    psPlane *center = psPlaneTransformGetCenter (chip->toFPA, tol);
+    psPlane *center = psPlaneTransformGetCenter (toTPA, tol);
 
     // create wcs transform from toFPA, resulting transformation has units of microns/pixel
     // adjust wcs transform to use center as reference coordinate
-    psPlaneTransformSetCenter (wcs->trans, chip->toFPA, center->x, center->y);
+    psPlaneTransformSetCenter (wcs->trans, toTPA, center->x, center->y);
 
     // calculated center is crpix1,2
@@ -595 +614 @@
     wcs->cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
     wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
+
+    psFree (toTPA);
 
     return wcs;