Changeset 15884 for trunk/psModules/src/astrom/pmAstrometryTable.c

Timestamp:

Dec 19, 2007, 8:57:05 AM (18 years ago)

Author:

eugene

Message:

substantial work on astrometry model tables

File:

• trunk/psModules/src/astrom/pmAstrometryTable.c (modified) (7 diffs)

trunk/psModules/src/astrom/pmAstrometryTable.c

@@ -1 @@
-/** @file  pmAstrometryObjects.c
+/** @file  pmAstrometryTable.c
  *
- *  @brief This file defines the basic types for matching objects
- *  based on their astrometry.
+ *  @brief Functions to read and write astrometric information using the generic astrometry
+ *  model.
+ * 
+ *  The generic model does not specify the location of the boresite on the sky, and it includes
+ *  a model for the rotator and motion of the boresite.
  *
  *  @ingroup AstroImage
  *
  *  @author EAM, IfA
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-12-19 18:57:05 $
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-11-21 12:01:10 $
- *
- *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
+ *  Copyright 2007 Institute for Astronomy, University of Hawaii
  */
 
@@ -38 +40 @@
 #include "pmFPAExtent.h"
 #include "pmAstrometryWCS.h"
+#include "pmAstrometryUtils.h"
 #include "pmAstrometryRegions.h"
 #include "pmAstrometryTable.h"
+
+# define REQUIRE(TEST,MESSAGE){ if (!(TEST)) { psAbort (MESSAGE); }}
 
 /********************* CheckDataStatus functions *****************************/
@@ -96 +101 @@
     // write the full table in one pass: require the level to be FPA
     if (view->chip != -1) {
-        psError(PS_ERR_IO, false, "Astrometry must be written at the chip level");
+        psError(PS_ERR_IO, false, "Astrometry must be written at the FPA level");
         return false;
     }
@@ -222 +227 @@
     psMetadataAddF32(row,    PS_LIST_TAIL, "XSCALE",   PS_META_REPLACE, "", file->fpa->toSky->Xs * PM_DEG_RAD);
     psMetadataAddF32(row,    PS_LIST_TAIL, "YSCALE",   PS_META_REPLACE, "", file->fpa->toSky->Ys * PM_DEG_RAD);
-    psMetadataAddF32(row,    PS_LIST_TAIL, "XREF",     PS_META_REPLACE, "", file->fpa->toSky->R * PM_DEG_RAD);
-    psMetadataAddF32(row,    PS_LIST_TAIL, "YREF",     PS_META_REPLACE, "", file->fpa->toSky->D * PM_DEG_RAD);
+    psMetadataAddF32(row,    PS_LIST_TAIL, "XREF",     PS_META_REPLACE, "", file->fpa->toSky->R  * PM_DEG_RAD);
+    psMetadataAddF32(row,    PS_LIST_TAIL, "YREF",     PS_META_REPLACE, "", file->fpa->toSky->D  * PM_DEG_RAD);
 
     psArrayAdd (table, 100, row);
@@ -313 +318 @@
 
         // set the chip name
-        char *chiprule = psStringCopy ("CHIP.{CHIP.NAME}"); 
+        char *chiprule = psStringCopy ("{CHIP.NAME}"); 
         char *chipname = pmFPAfileNameFromRule (chiprule, file, view);
 
@@ -327 +332 @@
                 psMetadataAddF32(row,    PS_LIST_TAIL, "MAXY",     PS_META_REPLACE, "range", region->y1);
     
-                psMetadataAddS32(row,    PS_LIST_TAIL, "NX",       PS_META_REPLACE, "", chip->toFPA->x->nX);
-                psMetadataAddS32(row,    PS_LIST_TAIL, "NY",       PS_META_REPLACE, "", chip->toFPA->x->nY);
+                psMetadataAddS32(row,    PS_LIST_TAIL, "NX",       PS_META_REPLACE, "", i);
+                psMetadataAddS32(row,    PS_LIST_TAIL, "NY",       PS_META_REPLACE, "", j);
                 psMetadataAddF32(row,    PS_LIST_TAIL, "POLY_X",   PS_META_REPLACE, "", chip->toFPA->x->coeff[i][j]);
                 psMetadataAddF32(row,    PS_LIST_TAIL, "POLY_Y",   PS_META_REPLACE, "", chip->toFPA->y->coeff[i][j]);
@@ -357 +362 @@
 
 bool pmAstromReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    pmFPA *fpa = file->fpa;
-
-    if (view->chip == -1) {
-        return pmAstromReadFPA(fpa, view, file, config);
-    }
-
-    if (view->chip >= fpa->chips->n) {
-        psAbort("Programming error: view does not apply to FPA.");
-    }
-    pmChip *chip = fpa->chips->data[view->chip];
-
-    if (view->cell == -1) {
-        return pmAstromReadChip(chip, view, file, config);
-    }
-
-    psError(PS_ERR_IO, false, "Astrometry must be read at the chip level");
-    return false;
-}
-
-// read in fpa-level Astrometry files for this FPA, and all chips
-bool pmAstromReadFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-    bool success = true;                // Was everything successful?
+    {
+
+        // write the full table in one pass: require the level to be FPA
+        if (view->chip != -1) {
+            psError(PS_ERR_IO, false, "Astrometry must be read at the FPA level");
+            return false;
+        }
+
+        if (!pmAstromReadFPA (file)) {
+            psError(PS_ERR_IO, false, "Failed to read Astrometry for fpa");
+            return false;
+        }
+        return true;
+    }
+
+// read out all chip-level Astrometry data for this FPA
+bool pmAstromReadFPA (pmFPAfile *file) {
+
+    if (!pmAstromReadPHU (file)) {
+        psError(PS_ERR_IO, false, "Failed to read PHU for Astrometry table");
+        return false;
+    }
+
+    if (!pmAstromReadChips (file)) {
+        psError(PS_ERR_IO, false, "Failed to read Astrometry for chips");
+        return false;
+    }
+
+    if (!pmAstromReadFP (file)) {
+        psError(PS_ERR_IO, false, "Failed to read Sky for Astrometry table");
+        return false;
+    }
+
+    // NOTE : TP must come after FP as it applies the POS, ROT boresite corrections to the
+    // transformation determined in FP
+    if (!pmAstromReadTP (file)) {
+        psError(PS_ERR_IO, false, "Failed to read Sky for Astrometry table");
+        return false;
+    }
+
+    if (!pmAstromReadSky (file)) {
+        psError(PS_ERR_IO, false, "Failed to read Sky for Astrometry table");
+        return false;
+    }
+
+    return true;
+}
+
+bool pmAstromReadPHU (pmFPAfile *file) { 
+
+    fprintf (stderr, "not sure we need to read PHU\n");
+    return true;
+}
+
+int pmConceptsChipNumberFromName (pmFPA *fpa, char *name) {
+
     for (int i = 0; i < fpa->chips->n; i++) {
-        pmChip *chip = fpa->chips->data[i];
-        success &= pmAstromReadChip(chip, view, file, config);
-    }
-    return success;
-}
-
-// read in chip-level Astrometry data for this chip
-bool pmAstromReadChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
+        pmChip *chip = fpa->chips->data[i];
+        if (!chip) continue;
+        char *thisone = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME");
+        if (!thisone) continue;
+        if (!strcmp (name, thisone)) return (i);
+    }
+    return -1;
+}
+
+pmChip *pmConceptsChipFromName (pmFPA *fpa, char *name) {
+
+    for (int i = 0; i < fpa->chips->n; i++) {
+        pmChip *chip = fpa->chips->data[i];
+        if (!chip) continue;
+        char *thisone = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME");
+        if (!thisone) continue;
+        if (!strcmp (name, thisone)) return (chip);
+    }
+    return NULL;
+}
+
+// first layer converts Chip to Focal Plane
+bool pmAstromReadChips (pmFPAfile *file) {
+
     bool status;
 
-    // lookup the EXTNAME values used for table data and image header segments
-    char *rule = NULL;
-
-    // Menu of EXTNAME rules
-    psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
-    if (!menu) {
-        psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
-        return false;
-    }
-
-    // EXTNAME for image header
-    rule = psMetadataLookupStr(&status, menu, "REF.ASTROM");
-    if (!rule) {
-        psError(PS_ERR_UNKNOWN, false, "missing entry for REF.ASTROM in EXTNAME.RULES in camera.config");
-        return false;
-    }
-
-    // XXX finish this: uncomment when used
-# if (0)
-    char *extname = pmFPAfileNameFromRule (rule, file, view);
-# endif
-
-    // read the chip elements in the following form:
-    // chipID, ...?  I forget the form: look at ICD
-
-    return true;
-}
-
+    // set FITS cursor
+    if (!psFitsMoveExtName (file->fits, "CHIPS")) {
+        psError(PS_ERR_IO, false, "missing CHIPS extension in astrometry table\n");
+        return false;
+    }
+
+    // free exising tranformations in prep for new alloc below
+    for (int i = 0; i < file->fpa->chips->n; i++) {
+        pmChip *chip = file->fpa->chips->data[i];
+        psFree (chip->toFPA);
+        chip->toFPA = NULL;
+    }
+
+    // XXX do I need anything from the header?
+    psMetadata *header = psFitsReadHeader(NULL, file->fits); // The FITS header
+    if (!header) psAbort("cannot read table header");
+
+    // load the full table in one shot
+    psArray *table = psFitsReadTable (file->fits);
+    if (!table) psAbort("cannot read table");
+    fprintf (stderr, "read %ld rows from FP\n", table->n);
+
+    // parse the table entries
+    for (int i = 0; i < table->n; i++) {
+        psMetadata *row = table->data[i];
+
+        // name of the chip for this row. 
+        char *chipname = psMetadataLookupStr (&status, row, "SEGMENT");
+
+        // get chip from name
+        pmChip *chip = pmConceptsChipFromName (file->fpa, chipname);
+        REQUIRE (chip, "invalid chip name");
+
+        // define the toFPA transform if not already defined
+        int nX = psMetadataLookupS32(&status, row, "NXORDER"); REQUIRE (status, "missing NXORDER");
+        int nY = psMetadataLookupS32(&status, row, "NYORDER"); REQUIRE (status, "missing NYORDER");
+        if (chip->toFPA == NULL) {
+            chip->toFPA = psPlaneTransformAlloc(nX, nY);
+        } else {
+            REQUIRE (chip->toFPA->x->nX == nX, "mismatch in chip order");
+            REQUIRE (chip->toFPA->x->nY == nY, "mismatch in chip order");
+            REQUIRE (chip->toFPA->y->nX == nX, "mismatch in chip order");
+            REQUIRE (chip->toFPA->y->nY == nY, "mismatch in chip order");
+        }
+
+        int ix = psMetadataLookupS32(&status, row, "XORDER");  REQUIRE (status, "missing XORDER");
+        int iy = psMetadataLookupS32(&status, row, "YORDER");  REQUIRE (status, "missing YORDER");
+
+        // XXX do I need to use this or can i rely on the camera concepts?
+        // minX = psMetadataLookupF32(&status, row, "MINX");
+        // maxX = psMetadataLookupF32(&status, row, "MAXX");
+        // minY = psMetadataLookupF32(&status, row, "MINY");
+        // maxY = psMetadataLookupF32(&status, row, "MAXY");
+
+        // XXX need to include mask values
+        chip->toFPA->x->coeff[ix][iy]    = psMetadataLookupF32(&status, row, "POLY_X");
+        chip->toFPA->y->coeff[ix][iy]    = psMetadataLookupF32(&status, row, "POLY_Y");
+        chip->toFPA->x->coeffErr[ix][iy] = psMetadataLookupF32(&status, row, "ERROR_X");
+        chip->toFPA->y->coeffErr[ix][iy] = psMetadataLookupF32(&status, row, "ERROR_Y");
+    }
+
+    // convert the toFPA transfomations to fromFPA transformations
+    for (int i = 0; i < file->fpa->chips->n; i++) {
+        pmChip *chip = file->fpa->chips->data[i];
+        psRegion *region = pmChipPixels (chip);
+        psFree (chip->fromFPA);
+        chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50);
+        psFree (region);
+    }
+
+    psFree (table);
+    psFree (header);
+    return true;
+}
+
+// second layer converts Focal Plane to Tangent Plane (unrotated)
+bool pmAstromReadFP (pmFPAfile *file) {
+
+    bool status;
+
+    if (!psFitsMoveExtName (file->fits, "FP")) {
+        psError(PS_ERR_IO, false, "missing FP extension in astrometry table\n");
+        return false;
+    }
+
+    psMetadata *header = psFitsReadHeader(NULL, file->fits); // The FITS header
+    if (!header) psAbort("cannot read table header");
+
+    // there is only one transformation in this table; the order is defined in the header
+    int nX = psMetadataLookupS32(&status, header, "NXORDER"); REQUIRE (status, "missing NXORDER");
+    int nY = psMetadataLookupS32(&status, header, "NYORDER"); REQUIRE (status, "missing NYORDER");
+
+    // allocate the new transformation
+    psFree (file->fpa->toTPA);
+    file->fpa->toTPA = psPlaneTransformAlloc(nX, nY);
+
+    // free & DON'T allocate the new transformation fromTPA; this is created after the boresite
+    // is applied in pmAstromReadTP
+    psFree (file->fpa->fromTPA);
+    file->fpa->fromTPA = NULL;
+
+    // read the complete table data at one shot
+    psArray *table = psFitsReadTable (file->fits);
+    fprintf (stderr, "read %ld rows from FP\n", table->n);
+
+    // parse the table
+    for (int i = 0; i < table->n; i++) {
+        psMetadata *row = table->data[i];
+        int ix = psMetadataLookupS32(&status, row, "XORDER"); REQUIRE (status, "missing XORDER");
+        int iy = psMetadataLookupS32(&status, row, "YORDER"); REQUIRE (status, "missing YORDER");
+        file->fpa->toTPA->x->coeff[ix][iy]    = psMetadataLookupF32(&status, row, "POLY_X"); REQUIRE (status, "missing POLY_X");
+        file->fpa->toTPA->y->coeff[ix][iy]    = psMetadataLookupF32(&status, row, "POLY_Y"); REQUIRE (status, "missing POLY_Y");
+        file->fpa->toTPA->x->coeffErr[ix][iy] = psMetadataLookupF32(&status, row, "ERROR_X"); REQUIRE (status, "missing ERROR_X");
+        file->fpa->toTPA->y->coeffErr[ix][iy] = psMetadataLookupF32(&status, row, "ERROR_Y"); REQUIRE (status, "missing ERROR_Y");
+    }
+
+    psFree (table);
+    psFree (header);
+    return true;
+}
+
+// third layer applies boresite corrections and converts tangent plane to sky
+bool pmAstromReadTP (pmFPAfile *file) {
+    
+    bool status;
+
+    // these externally supplied values are used to set the final transformation terms
+    double RA  = psMetadataLookupF64 (&status, file->fpa->concepts, "FPA.RA"); REQUIRE (status, "missing FPA.RA");
+    double DEC = psMetadataLookupF64 (&status, file->fpa->concepts, "FPA.DEC"); REQUIRE (status, "missing FPA.DEC");
+    double POS = psMetadataLookupF64 (&status, file->fpa->concepts, "FPA.POSANGLE"); REQUIRE (status, "missing FPA.POSANGLE");
+    double ROT = psMetadataLookupF64 (&status, file->fpa->concepts, "FPA.ROTANGLE"); REQUIRE (status, "missing FPA.ROTANGLE");
+
+    if (!psFitsMoveExtName (file->fits, "TP")) {
+        psError(PS_ERR_IO, false, "missing TP extension in astrometry table\n");
+        return false;
+    }
+
+    psMetadata *header = psFitsReadHeader(NULL, file->fits); // The FITS header
+    if (!header) psAbort("cannot read table header");
+
+    psArray *table = psFitsReadTable (file->fits);
+    fprintf (stderr, "read %ld rows from TP\n", table->n);
+    if (table->n != 1) psAbort("invalid number of rows in TP table (%ld)", table->n);
+
+    psMetadata *row = table->data[0];
+
+    // get projection scale; center is supplied
+    float Xs = psMetadataLookupF32(&status, row, "XSCALE") * PM_RAD_DEG; REQUIRE (status, "missing XSCALE");
+    float Ys = psMetadataLookupF32(&status, row, "YSCALE") * PM_RAD_DEG; REQUIRE (status, "missing YSCALE");
+
+    // allocate a new toSky projection
+    psFree (file->fpa->toSky);
+    file->fpa->toSky = psProjectionAlloc (RA, DEC, Xs, Ys, PS_PROJ_WRP);
+
+    // get boresite correction terms.  L,M,T define an ellipse along which the boresite travels
+    double R_L = psMetadataLookupF32(&status, row, "ROTATOR_L"); REQUIRE (status, "missing ROTATOR_L");
+    double R_M = psMetadataLookupF32(&status, row, "ROTATOR_M"); REQUIRE (status, "missing ROTATOR_M");
+    double R_T = psMetadataLookupF32(&status, row, "ROTATOR_T"); REQUIRE (status, "missing ROTATOR_T");
+
+    // the position of the boresite on the ellipse is the ROTANGLE - ROT_ZERO
+    double R_0 = psMetadataLookupF32(&status, row, "ROT_ZERO"); REQUIRE (status, "missing ROT_ZERO");
+
+    // find the current position of the boresite for this image
+    double Lo = R_L*cos(ROT - R_0)*cos(R_T) + R_M*sin(ROT - R_0)*sin(R_T);
+    double Mo = R_M*sin(ROT - R_0)*cos(R_T) - R_L*cos(ROT - R_0)*sin(R_T);
+
+    // apply new reference pixel to transformation
+    psPlaneTransform *tmpToTPA = psPlaneTransformSetCenter (NULL, file->fpa->toTPA, Lo, Mo);
+    psFree (file->fpa->toTPA);
+    
+    // apply POSANGLE 
+    double Po = psMetadataLookupF32(&status, row, "POS_ZERO"); REQUIRE (status, "missing POS_ZERO");
+    file->fpa->toTPA = psPlaneTransformRotate (NULL, tmpToTPA, POS - Po);
+    psFree (tmpToTPA);
+
+    psRegion *region = pmAstromFPAExtent (file->fpa);
+    psFree (file->fpa->fromTPA);
+    file->fpa->fromTPA = psPlaneTransformInvert(NULL, file->fpa->toTPA, *region, 50);
+
+    psFree (region);
+    psFree (table);
+    psFree (header);
+    return (true);
+}
+
+// first layer is the sky
+bool pmAstromReadSky (pmFPAfile *file) {
+
+    if (!psFitsMoveExtName (file->fits, "SKY")) {
+        psError(PS_ERR_IO, false, "missing SKY extension in astrometry table\n");
+        return false;
+    }
+
+    psMetadata *header = psFitsReadHeader(NULL, file->fits); // The FITS header
+    if (!header) psAbort("cannot read table header");
+
+    psArray *table = psFitsReadTable (file->fits);
+    fprintf (stderr, "read %ld rows from SKY\n", table->n);
+    if (table->n != 1) psAbort("invalid number of rows in SKY table (%ld)", table->n);
+
+    // XXX not much information of interest in this table...
+
+    psFree (table);
+    psFree (header);
+    return (true);
+}
+