Index: trunk/psModules/src/astrom/pmAstrometryTable.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryTable.c (revision 15865) +++ trunk/psModules/src/astrom/pmAstrometryTable.c (revision 15884) @@ -1,15 +1,17 @@ -/** @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,6 +40,9 @@ #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,5 +101,5 @@ // 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,6 +227,6 @@ 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,5 +318,5 @@ // set the chip name - char *chiprule = psStringCopy ("CHIP.{CHIP.NAME}"); + char *chiprule = psStringCopy ("{CHIP.NAME}"); char *chipname = pmFPAfileNameFromRule (chiprule, file, view); @@ -327,6 +332,6 @@ 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,67 +362,292 @@ 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); +} + Index: trunk/psModules/src/astrom/pmAstrometryTable.h =================================================================== --- trunk/psModules/src/astrom/pmAstrometryTable.h (revision 15865) +++ trunk/psModules/src/astrom/pmAstrometryTable.h (revision 15884) @@ -4,6 +4,6 @@ * @author EAM, IfA * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-11-21 09:59:40 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-12-19 18:57:05 $ * Copyright 2006 Institute for Astronomy, University of Hawaii */ @@ -27,7 +27,14 @@ bool pmAstromWriteChips (pmFPAfile *file); +int pmConceptsChipNumberFromName (pmFPA *fpa, char *name); +pmChip *pmConceptsChipFromName (pmFPA *fpa, char *name); + bool pmAstromReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config); -bool pmAstromReadFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config); -bool pmAstromReadChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config); +bool pmAstromReadFPA (pmFPAfile *file); +bool pmAstromReadPHU (pmFPAfile *file); +bool pmAstromReadChips (pmFPAfile *file); +bool pmAstromReadFP (pmFPAfile *file); +bool pmAstromReadTP (pmFPAfile *file); +bool pmAstromReadSky (pmFPAfile *file); /// @} Index: trunk/psModules/src/astrom/pmAstrometryUtils.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryUtils.c (revision 15865) +++ trunk/psModules/src/astrom/pmAstrometryUtils.c (revision 15884) @@ -7,6 +7,6 @@ * @author EAM, IfA * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-09 19:27:04 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-12-19 18:57:05 $ * * Copyright 2006 Institute for Astronomy, University of Hawaii @@ -129,4 +129,41 @@ } +// rotate a transformation L(x,y) by theta +psPlaneTransform *psPlaneTransformRotate (psPlaneTransform *output, psPlaneTransform *input, double theta) +{ + // validate fit order + if (output == NULL) { + output = psPlaneTransformAlloc(input->x->nX, input->x->nY); + } + + /* given the polynomial transformations: + * L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j and + * M(x,y) = \sum_i \sum_j B_{i,j} x^i y^j + * we can rotate L,M to L',M' by applying the rotation matrix (c,s),(-s,c). + * the resulting terms of L and M are: + * A'_{i,j} = c A_{i,j} + s B_{i,j} + * B'_{i,j} = c B_{i,j} - s A_{i,j} + */ + + float cs = cos(theta); + float sn = sin(theta); + + for (int i = 0; i <= input->x->nX; i++) { + for (int j = 0; j <= input->x->nY; j++) { + // XXX what about inconsistent x and y masking? + output->x->coeffMask[i][j] = input->x->coeffMask[i][j]; + output->y->coeffMask[i][j] = input->y->coeffMask[i][j]; + if (output->x->coeffMask[i][j]) { + output->x->coeff[i][j] = 0.0; + output->y->coeff[i][j] = 0.0; + } else { + output->x->coeff[i][j] = cs*input->x->coeff[i][j] + sn*input->y->coeff[i][j]; + output->y->coeff[i][j] = cs*input->y->coeff[i][j] - sn*input->x->coeff[i][j]; + } + } + } + return output; +} + // construct a psPlaneTransform which is the identify transformation psPlaneTransform *psPlaneTransformIdentity (int order) Index: trunk/psModules/src/astrom/pmAstrometryUtils.h =================================================================== --- trunk/psModules/src/astrom/pmAstrometryUtils.h (revision 15865) +++ trunk/psModules/src/astrom/pmAstrometryUtils.h (revision 15884) @@ -4,6 +4,6 @@ * @author EAM, IfA * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-01-24 02:54:14 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-12-19 18:57:05 $ * Copyright 2006 Institute for Astronomy, University of Hawaii */ @@ -17,4 +17,5 @@ psPlane *psPlaneTransformGetCenter (psPlaneTransform *trans, double tol); psPlaneTransform *psPlaneTransformSetCenter (psPlaneTransform *output, psPlaneTransform *input, double Xo, double Yo); +psPlaneTransform *psPlaneTransformRotate (psPlaneTransform *output, psPlaneTransform *input, double theta); psPlaneTransform *psPlaneTransformIdentity (int order); Index: trunk/psModules/src/psmodules.h =================================================================== --- trunk/psModules/src/psmodules.h (revision 15865) +++ trunk/psModules/src/psmodules.h (revision 15884) @@ -72,4 +72,5 @@ #include #include +#include #include