Changeset 12483 for trunk/psModules/src/astrom/pmAstrometryWCS.c

Timestamp:

Mar 18, 2007, 12:05:43 PM (19 years ago)

Author:

magnier

Message:

fixed adjustment to coordinate center when applying a common f.p. scale; get chip and fpa dimensions from concepts

File:

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

trunk/psModules/src/astrom/pmAstrometryWCS.c

@@ -7 +7 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-02-07 23:58:17 $
+ *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-03-18 22:05:43 $
  *
  *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -20 +20 @@
 #include <pslib.h>
 #include "pmFPA.h"
+#include "pmFPAExtent.h"
 #include "pmAstrometryWCS.h"
 #include "pmAstrometryUtils.h"
@@ -82 +83 @@
 
 // convert toFPA / toSky components to traditional WCS
-bool pmAstromReadBilevelMosaic (pmFPA *fpa, const psMetadata *header)
+bool pmAstromReadBilevelMosaic (pmFPA *fpa, const psMetadata *header, double pixelScale)
 {
     pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
@@ -89 +90 @@
     }
 
-    bool status = pmAstromWCSBileveltoFPA (fpa, wcs);
+    bool status = pmAstromWCSBileveltoFPA (fpa, wcs, pixelScale);
 
     psFree (wcs);
@@ -404 +405 @@
         double rX = toSky->Xs / fpa->toSky->Xs;
         double rY = toSky->Ys / fpa->toSky->Ys;
+
         for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
             for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
@@ -411 +413 @@
         }
 
-        // adjust for common toSky, toTPA for mosaic:
+        // adjust reference pixel for new toSky reference coordinate
         // find the FPA coordinate of 0,0 for this chip.
-        psPlane *chip = psPlaneAlloc();
-        psPlane *fp = psPlaneAlloc();
+        psPlane *fpOld = psPlaneAlloc();
+        psPlane *fpNew = psPlaneAlloc();
         psPlane *tp = psPlaneAlloc();
         psSphere *sky = psSphereAlloc();
-        chip->x = chip->y = 0;
-
-        psPlaneTransformApply (fp, toFPA, chip); // find the focal-plane coordinate of this chip's 0,0 coordinate
-        psPlaneTransformApply (tp, fpa->toTPA, fp);
-        psDeproject (sky, tp, toSky); // find the RA,DEC coord of the focal-plane coordinate
+
+        sky->r = toSky->R;
+        sky->d = toSky->D;
         psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
-        psPlaneTransformApply (fp, fpa->fromTPA, tp);
-
-        toFPA->x->coeff[0][0] = fp->x;
-        toFPA->y->coeff[0][0] = fp->y;
+        psPlaneTransformApply (fpOld, fpa->fromTPA, tp);
+
+        sky->r = fpa->toSky->R;
+        sky->d = fpa->toSky->D;
+        psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
+        psPlaneTransformApply (fpNew, fpa->fromTPA, tp);
+
+        toFPA->x->coeff[0][0] -= fpNew->x - fpOld->x;
+        toFPA->y->coeff[0][0] -= fpNew->y - fpOld->y;
 
         psFree (sky);
         psFree (tp);
-        psFree (fp);
-        psFree (chip);
+        psFree (fpOld);
+        psFree (fpNew);
+
         psFree (toSky);
     }
@@ -437 +443 @@
     // free an existing toFPA structure
     chip->toFPA = toFPA;
-    // XXX this needs to perform the full (non-linear) inversion
-    // XXX we need to pull the region from the chip metadata
-    // is chip trimsec defined?  do we need to ensure the
-    // XXX should this function be defined for the CELL, not the CHIP?
-    // psRegion region = psMetadataLookupXXX (chip->concepts, "CHIP.TRIMSEC");
-    psRegion region = psRegionSet (0, 4000, 0, 4000);
-    chip->fromFPA = psPlaneTransformInvert(NULL, toFPA, region, 50);
+    // NOTE that region specifies the FP region in PIXELS covered by the chip
+    psRegion *region = pmChipExtent (chip);
+    region->x1 -= region->x0;
+    region->y1 -= region->y0;
+    region->x0 = 0;
+    region->y0 = 0;
+    psFree (chip->fromFPA);
+    chip->fromFPA = psPlaneTransformInvert(NULL, toFPA, *region, 50);
+    psFree (region);
 
     // this can take a very long time...
@@ -475 +483 @@
     // create transformation with 0,0 reference pixel and units of microns/pixel
     chip->toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
-
-    // XXX this needs to perform the full (non-linear) inversion
-    // XXX we need to pull the region from the chip metadata
-    // is chip trimsec defined?  do we need to ensure the
-    // XXX should this function be defined for the CELL, not the CHIP?
-    // psRegion region = psMetadataLookupXXX (chip->concepts, "CHIP.TRIMSEC");
-    psRegion region = psRegionSet (0, 4000, 0, 4000);
-    chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, region, 50);
+    psRegion *region = pmChipExtent (chip);
+    region->x1 -= region->x0;
+    region->y1 -= region->y0;
+    region->x0 = 0;
+    region->y0 = 0;
+    psFree (chip->fromFPA);
+    chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50);
+    psFree (region);
 
     return true;
@@ -488 +496 @@
 
 // convert a pmAstromWCS structure representing a bilevel mosaic into corresponding fpa elements
-bool pmAstromWCSBileveltoFPA (pmFPA *fpa, const pmAstromWCS *wcs)
+bool pmAstromWCSBileveltoFPA (pmFPA *fpa, const pmAstromWCS *wcs, double pixelScale)
 {
     // projection from TPA (microns) to SKY (radians)
@@ -505 +513 @@
     }
 
-    // XXX this needs to perform the full (non-linear) inversion
-    // XXX we need to pull the region from the chip metadata
-    // is chip trimsec defined?  do we need to ensure the
-    // XXX should this function be defined for the CELL, not the CHIP?
-    // psRegion region = psMetadataLookupXXX (chip->concepts, "CHIP.TRIMSEC");
-    psRegion region = psRegionSet (0, 4000, 0, 4000);
-    fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, region, 50);
+    // XXX this is still wrong: the pmFPAExtent function returns the FP dimensions in pixel units
+    // relative to the 0,0 corner of chip 0,0.  for now, just recenter...
+    // XXX I need to have a function which loops over all cells, determines the pixel dimensions 
+    // for each cell, converts them to chip pixels, then uses the fpa astrometry structures
+    // to get the total dimensions of the fpa in fpa units.  equivalent to pmFPAExtent
+    psRegion *region = pmFPAExtent (fpa);
+    region->x0 *= pixelScale;
+    region->x1 *= pixelScale;
+    region->y0 *= pixelScale;
+    region->y1 *= pixelScale;
+    int xCenter = 0.5*(region->x0 + region->x1);
+    int yCenter = 0.5*(region->y0 + region->y1);
+    region->x0 -= xCenter;
+    region->x1 -= xCenter;
+    region->y0 -= yCenter;
+    region->y1 -= yCenter;
+    psFree (fpa->fromTPA);
+    fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, *region, 50);
+    psFree (region);
 
     return true;