Index: /trunk/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 15664) +++ /trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 15665) @@ -8,6 +8,6 @@ * @author EAM, IfA * -* @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ -* @date $Date: 2007-09-29 02:28:04 $ +* @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ +* @date $Date: 2007-11-21 07:02:55 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -269,5 +269,5 @@ /****************************************************************************** -pmAstromRotateObj(old, center, angle): rotate the focal-plane coordinates +pmAstromRotateObj(old, center, angle, angle): rotate & scale the focal-plane coordinates about the center coordinate angle specified in radians ******************************************************************************/ @@ -275,5 +275,6 @@ const psArray *old, psPlane center, - double angle) + double angle, + double scale) { PS_ASSERT_PTR_NON_NULL(old, NULL); @@ -284,6 +285,6 @@ psArray *new = psArrayAlloc (old->n); - double cs = cos(angle); - double sn = sin(angle); + double cs = scale*cos(angle); + double sn = scale*sin(angle); double xCenter = center.x; double yCenter = center.y; @@ -326,4 +327,5 @@ // stats->offset = {0, 0, 0, 0}; stats->angle = 0.0; + stats->scale = 1.0; stats->minMetric = 0.0; stats->minVar = 0.0; @@ -703,4 +705,8 @@ center.y = 0.5*(yMin + yMax); + double minScale = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.SCALE"); + double maxScale = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MAX.SCALE"); + double delScale = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.DEL.SCALE"); + double minAngle = RAD_DEG*psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.ANGLE"); double maxAngle = RAD_DEG*psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MAX.ANGLE"); @@ -709,24 +715,27 @@ minStat->minMetric = 1e10; - for (double angle = minAngle; angle <= maxAngle; angle += delAngle) { - rot = pmAstromRotateObj (raw, center, angle); - newStat = pmAstromGridAngle (rot, ref, config); - newStat->angle = angle; - newStat->center = center; - - if (newStat->minMetric < minStat->minMetric) { - *minStat = *newStat; - psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric) *", - minStat->offset.x, minStat->offset.y, DEG_RAD*minStat->angle, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric)); - } else { - psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)", - newStat->offset.x, newStat->offset.y, DEG_RAD*newStat->angle, newStat->nMatch, newStat->nSigma, newStat->minVar, log10(newStat->minMetric)); - - } - psFree (newStat); - psFree (rot); - } - psLogMsg ("psModule.astrom.grid.match", 4, "grid best - offset: %7.2f,%7.2f @ %6.1f deg (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)", - minStat->offset.x, minStat->offset.y, DEG_RAD*minStat->angle, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric)); + for (double scale = minScale; scale <= maxScale; scale += delScale) { + for (double angle = minAngle; angle <= maxAngle; angle += delAngle) { + rot = pmAstromRotateObj (raw, center, angle, scale); + newStat = pmAstromGridAngle (rot, ref, config); + newStat->angle = angle; + newStat->scale = scale; + newStat->center = center; + + if (newStat->minMetric < minStat->minMetric) { + *minStat = *newStat; + psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric) *", + minStat->offset.x, minStat->offset.y, DEG_RAD*minStat->angle, minStat->scale, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric)); + } else { + psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)", + newStat->offset.x, newStat->offset.y, DEG_RAD*newStat->angle, newStat->scale, newStat->nMatch, newStat->nSigma, newStat->minVar, log10(newStat->minMetric)); + + } + psFree (newStat); + psFree (rot); + } + } + psLogMsg ("psModule.astrom.grid.match", 4, "grid best - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)", + minStat->offset.x, minStat->offset.y, DEG_RAD*minStat->angle, minStat->scale, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric)); // I need to decide if a solution is likely to be a good solution or just a mis-match @@ -756,5 +765,5 @@ int nBin, xBin, yBin; - rot = pmAstromRotateObj (raw, stats->center, stats->angle); + rot = pmAstromRotateObj (raw, stats->center, stats->angle, stats->scale); // sampling scale of the grid @@ -846,6 +855,6 @@ PS_ASSERT_POLY_NON_NULL(map->y, NULL); - double cs = cos (stat->angle); - double sn = sin (stat->angle); + double cs = stat->scale * cos (stat->angle); + double sn = stat->scale * sin (stat->angle); double dx = (map->x->coeff[0][0] - stat->center.x); Index: /trunk/psModules/src/astrom/pmAstrometryObjects.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 15664) +++ /trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 15665) @@ -4,6 +4,6 @@ * @author EAM, IfA * - * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-23 23:44:05 $ + * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-21 07:02:55 $ * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii */ @@ -68,4 +68,5 @@ psPlane center; ///< psPlane offset; ///< + double scale; ///< double angle; ///< double minMetric; ///< @@ -142,5 +143,6 @@ const psArray *old, psPlane center, - double angle + double angle, + double scale ); Index: /trunk/psModules/src/astrom/pmAstrometryTable.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryTable.c (revision 15664) +++ /trunk/psModules/src/astrom/pmAstrometryTable.c (revision 15665) @@ -8,6 +8,6 @@ * @author EAM, IfA * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-11-10 01:09:20 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-21 07:02:55 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -36,6 +36,8 @@ #include "pmFPAview.h" #include "pmFPAfile.h" +#include "pmFPAExtent.h" +#include "pmAstrometryWCS.h" +#include "pmAstrometryRegions.h" #include "pmAstrometryTable.h" - /********************* CheckDataStatus functions *****************************/ @@ -181,36 +183,172 @@ // write out chip-level Astrometry data for this chip +// XXX this function writes out astrometry elements for a chip->fpa->tp->sky layer set // XXX we need to ensure the chip and fpa level data are written // XXX can we write the rows one at a time (or groups at a time?) +// XXX allow this to be chip level, or just force to be fpa? is it valid at the chip level?? bool pmAstromWriteChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { - bool status; - // *** define the EXTNAME values used for the table data, and residual image segments *** - - // 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, "ASTROMETRY"); - if (!rule) { - psError(PS_ERR_UNKNOWN, false, "missing entry for ASTROMETRY in EXTNAME.RULES in camera.config"); - return false; - } - - // XXX need to finish this: uncomment when used - # if (0) - char *extname = pmFPAfileNameFromRule (rule, file, view); - # endif - - // write the chip elements in the following form: - // chipID, ...? I forget the form: look at ICD - + psArray *table; + psMetadata *row; + psMetadata *header; + + // write out the blank PHU to start the file + psFitsWriteBlank (file->fits, NULL, NULL); + + // first layer is the sky + header = psMetadataAlloc(); + psMetadataAddStr(header, PS_LIST_TAIL, "COORD", PS_META_REPLACE, "name of this layer", "SKY"); + psMetadataAddStr(header, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "NONE"); + psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region", "NONE"); + psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent", "NONE"); + + table = psArrayAllocEmpty (1); + row = psMetadataAlloc (); + psMetadataAddStr(row, PS_LIST_TAIL, "SEGMENT", PS_META_REPLACE, "name of this segment", "SKY"); + psMetadataAddStr(row, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "NONE"); + + psArrayAdd (table, 100, row); + psFree (row); + + if (!psFitsWriteTable (file->fits, header, table, "SKY")) { + psError(PS_ERR_IO, false, "writing sky data\n"); + psFree(table); + return false; + } + + psFree (table); + psFree (header); + + // second layer is the tangent plane + header = psMetadataAlloc(); + psMetadataAddStr(header, PS_LIST_TAIL, "COORD", PS_META_REPLACE, "name of this layer", "TANGENT_PLANE"); + psMetadataAddStr(header, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "SKY"); + psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region", "RECTANGLE"); + psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent", "PROJECTION"); + + table = psArrayAllocEmpty (1); + row = psMetadataAlloc (); + psMetadataAddStr(row, PS_LIST_TAIL, "SEGMENT", PS_META_REPLACE, "name of this segment", "TANGENT_PLANE"); + psMetadataAddStr(row, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "SKY"); + + psRegion *region = pmAstromFPAExtent (file->fpa); + psMetadataAddF32(row, PS_LIST_TAIL, "MINX", PS_META_REPLACE, "range", region->x0); + psMetadataAddF32(row, PS_LIST_TAIL, "MAXX", PS_META_REPLACE, "range", region->x1); + psMetadataAddF32(row, PS_LIST_TAIL, "MINY", PS_META_REPLACE, "range", region->y0); + psMetadataAddF32(row, PS_LIST_TAIL, "MAXY", PS_META_REPLACE, "range", region->y1); + + 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); + + psArrayAdd (table, 100, row); + psFree (row); + + if (!psFitsWriteTable (file->fits, header, table, "TP")) { + psError(PS_ERR_IO, false, "writing sky data\n"); + psFree(table); + return false; + } + + psFree (table); + psFree (header); + + // third layer is the focal plane + header = psMetadataAlloc(); + psMetadataAddStr(header, PS_LIST_TAIL, "COORD", PS_META_REPLACE, "name of this layer", "FOCAL_PLANE"); + psMetadataAddStr(header, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "TANGENT_PLANE"); + psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region", "RECTANGLE"); + psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent", "POLYNOMIAL"); + + table = psArrayAllocEmpty (1); + + // XXX is this or the tpa region correct? + region = pmAstromFPAExtent (file->fpa); + + for (int i = 0; i <= file->fpa->toTPA->x->nX; i++) { + for (int j = 0; j <= file->fpa->toTPA->x->nY; j++) { + row = psMetadataAlloc (); + psMetadataAddStr(row, PS_LIST_TAIL, "SEGMENT", PS_META_REPLACE, "name of this segment", "FOCAL_PLANE"); + psMetadataAddStr(row, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "TANGENT_PLANE"); + psMetadataAddF32(row, PS_LIST_TAIL, "MINX", PS_META_REPLACE, "range", region->x0); + psMetadataAddF32(row, PS_LIST_TAIL, "MAXX", PS_META_REPLACE, "range", region->x1); + psMetadataAddF32(row, PS_LIST_TAIL, "MINY", PS_META_REPLACE, "range", region->y0); + psMetadataAddF32(row, PS_LIST_TAIL, "MAXY", PS_META_REPLACE, "range", region->y1); + + psMetadataAddS32(row, PS_LIST_TAIL, "NX", PS_META_REPLACE, "", file->fpa->toTPA->x->nX); + psMetadataAddS32(row, PS_LIST_TAIL, "NY", PS_META_REPLACE, "", file->fpa->toTPA->x->nY); + psMetadataAddF32(row, PS_LIST_TAIL, "POLY_X", PS_META_REPLACE, "", file->fpa->toTPA->x->coeff[i][j]); + psMetadataAddF32(row, PS_LIST_TAIL, "POLY_Y", PS_META_REPLACE, "", file->fpa->toTPA->y->coeff[i][j]); + psMetadataAddF32(row, PS_LIST_TAIL, "ERROR_X", PS_META_REPLACE, "", file->fpa->toTPA->x->coeffErr[i][j]); + psMetadataAddF32(row, PS_LIST_TAIL, "ERROR_Y", PS_META_REPLACE, "", file->fpa->toTPA->y->coeffErr[i][j]); + + psArrayAdd (table, 100, row); + psFree (row); + } + } + + if (!psFitsWriteTable (file->fits, header, table, "FP")) { + psError(PS_ERR_IO, false, "writing sky data\n"); + psFree(table); + return false; + } + + psFree (table); + psFree (header); + + // fourth layer holds the chips + header = psMetadataAlloc(); + psMetadataAddStr(header, PS_LIST_TAIL, "COORD", PS_META_REPLACE, "name of this layer", "FOCAL_PLANE"); + psMetadataAddStr(header, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "TANGENT_PLANE"); + psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region", "RECTANGLE"); + psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent", "POLYNOMIAL"); + + table = psArrayAllocEmpty (1); + + pmFPAview *view = pmFPAviewAlloc (0); + + pmChip *chip = NULL; + while ((chip = pmFPAviewNextChip (view, file->fpa, 1)) != NULL) { + + psRegion *region = pmChipPixels (chip); + + char *chiprule = psStringCopy ("CHIP.{CHIP.NAME}"); + char *chipname = pmFPAfileNameFromRule (chiprule, file, view); + + for (int i = 0; i <= chip->toFPA->x->nX; i++) { + for (int j = 0; j <= chip->toFPA->x->nY; j++) { + row = psMetadataAlloc (); + // xXXX set chip name based on metadata + psMetadataAddStr(row, PS_LIST_TAIL, "SEGMENT", PS_META_REPLACE, "name of this segment", chipname); + psMetadataAddStr(row, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "TANGENT_PLANE"); + psMetadataAddF32(row, PS_LIST_TAIL, "MINX", PS_META_REPLACE, "range", region->x0); + psMetadataAddF32(row, PS_LIST_TAIL, "MAXX", PS_META_REPLACE, "range", region->x1); + psMetadataAddF32(row, PS_LIST_TAIL, "MINY", PS_META_REPLACE, "range", region->y0); + 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); + 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]); + psMetadataAddF32(row, PS_LIST_TAIL, "ERROR_X", PS_META_REPLACE, "", chip->toFPA->x->coeffErr[i][j]); + psMetadataAddF32(row, PS_LIST_TAIL, "ERROR_Y", PS_META_REPLACE, "", chip->toFPA->y->coeffErr[i][j]); + psArrayAdd (table, 100, row); + psFree (row); + } + } + psFree (chiprule); + psFree (chipname); + } + psFree (view); + + if (!psFitsWriteTable (file->fits, header, table, "CHIPS")) { + psError(PS_ERR_IO, false, "writing sky data\n"); + psFree(table); + return false; + } + + psFree (table); + psFree (header); return true; } @@ -267,7 +405,7 @@ // EXTNAME for image header - rule = psMetadataLookupStr(&status, menu, "ASTROMETRY"); + rule = psMetadataLookupStr(&status, menu, "REF.ASTROM"); if (!rule) { - psError(PS_ERR_UNKNOWN, false, "missing entry for ASTROMETRY in EXTNAME.RULES in camera.config"); + psError(PS_ERR_UNKNOWN, false, "missing entry for REF.ASTROM in EXTNAME.RULES in camera.config"); return false; } Index: /trunk/psModules/src/camera/pmFPAfile.c =================================================================== --- /trunk/psModules/src/camera/pmFPAfile.c (revision 15664) +++ /trunk/psModules/src/camera/pmFPAfile.c (revision 15665) @@ -439,4 +439,7 @@ if (!strcasecmp (type, "HEADER")) { return PM_FPA_FILE_HEADER; + } + if (!strcasecmp (type, "ASTROM")) { + return PM_FPA_FILE_ASTROM; } @@ -473,4 +476,6 @@ case PM_FPA_FILE_HEADER: return ("HEADER"); + case PM_FPA_FILE_ASTROM: + return ("ASTROM"); default: return ("NONE");