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

Timestamp:

Mar 21, 2007, 12:00:49 PM (19 years ago)

Author:

magnier

Message:

fixed the source of the focal-plane region; added astrometric extent functions

File:

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

trunk/psModules/src/astrom/pmAstrometryWCS.c

@@ -7 +7 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-03-18 22:05:43 $
+ *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-03-21 22:00:49 $
  *
  *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -23 +23 @@
 #include "pmAstrometryWCS.h"
 #include "pmAstrometryUtils.h"
+#include "pmAstrometryRegions.h"
 
 // the following functions support coordinate transformations direcly related to the FITS WCS
@@ -83 +84 @@
 
 // convert toFPA / toSky components to traditional WCS
-bool pmAstromReadBilevelMosaic (pmFPA *fpa, const psMetadata *header, double pixelScale)
+// we require the header to have NAXIS1,NAXIS2, the field of the FPA 
+// the center of the TPA/Sky projection is 0.5*(NAXIS1,NAXIS2)
+bool pmAstromReadBilevelMosaic (pmFPA *fpa, const psMetadata *header)
 {
     pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
     if (!wcs) {
-        return false;
-    }
-
-    bool status = pmAstromWCSBileveltoFPA (fpa, wcs, pixelScale);
+        psError(PS_ERR_UNKNOWN, false, "failure to determine WCS terms from header");
+        return false;
+    }
+
+    bool status1 = false;
+    bool status2 = false;
+    int Nx = psMetadataLookupS32 (&status1, header, "NAXIS1");
+    int Ny = psMetadataLookupS32 (&status2, header, "NAXIS2");
+    if (!status1 || !status2) {
+        psFree (wcs);
+        psError(PS_ERR_UNKNOWN, false, "missing required FPA size in header");
+        return false;
+    }
+
+    psRegion region = psRegionSet (-0.5*Nx, +0.5*Nx, -0.5*Ny, +0.5*Ny);
+    bool status = pmAstromWCSBileveltoFPA (fpa, wcs, region);
 
     psFree (wcs);
@@ -113 +128 @@
     pmAstromWCS *wcs = pmAstromWCSBilevelMosaicFromFPA (fpa, tol);
 
+    // we need to specify the dimensions of the FPA
+    // if we have chips defined, we can do
+    psRegion *region = pmAstromFPAExtent (fpa);
+    int Nx = region->x1 - region->x0;
+    int Ny = region->y1 - region->y0;
+    psMetadataAddS32 (header, PS_LIST_TAIL, "NAXIS1", PS_META_REPLACE, "Mosaic Dimensions", Nx);
+    psMetadataAddS32 (header, PS_LIST_TAIL, "NAXIS2", PS_META_REPLACE, "Mosaic Dimensions", Ny);
+
     pmAstromWCStoHeader (header, wcs);
 
+    psFree (region);
     psFree (wcs);
     return true;
@@ -442 +466 @@
 
     // free an existing toFPA structure
+    psFree (chip->toFPA);
     chip->toFPA = toFPA;
-    // 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;
+
+    // determine the inverse transformation: we need the chip pixels covered by this transform
+    psRegion *region = pmChipPixels (chip);
+
     psFree (chip->fromFPA);
-    chip->fromFPA = psPlaneTransformInvert(NULL, toFPA, *region, 50);
+    chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50);
     psFree (region);
 
@@ -482 +505 @@
 
     // create transformation with 0,0 reference pixel and units of microns/pixel
+    psFree (chip->toFPA);
     chip->toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
-    psRegion *region = pmChipExtent (chip);
-    region->x1 -= region->x0;
-    region->y1 -= region->y0;
-    region->x0 = 0;
-    region->y0 = 0;
+
+    // determine the inverse transformation: we need the chip pixels covered by this transform
+    psRegion *region = pmChipPixels (chip);
+
     psFree (chip->fromFPA);
     chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50);
@@ -496 +519 @@
 
 // convert a pmAstromWCS structure representing a bilevel mosaic into corresponding fpa elements
-bool pmAstromWCSBileveltoFPA (pmFPA *fpa, const pmAstromWCS *wcs, double pixelScale)
+bool pmAstromWCSBileveltoFPA (pmFPA *fpa, const pmAstromWCS *wcs, psRegion region)
 {
     // projection from TPA (microns) to SKY (radians)
@@ -513 +536 @@
     }
 
-    // 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;
+    // the transformation used the region to define the inversion grid
+    // the region defines the FPA pixels covered by the tranformation
     psFree (fpa->fromTPA);
-    fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, *region, 50);
-    psFree (region);
-
+    fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, region, 50);
     return true;
 }