Index: /trunk/psModules/src/astrom/pmAstrometryTable.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryTable.c (revision 15357) +++ /trunk/psModules/src/astrom/pmAstrometryTable.c (revision 15357) @@ -0,0 +1,272 @@ +/** @file pmAstrometryObjects.c +* +* @brief This file defines the basic types for matching objects +* based on their astrometry. +* +* @ingroup AstroImage +* +* @author EAM, IfA +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2007-10-23 20:52:42 $ +* +* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +*/ + +#ifdef HAVE_CONFIG_H +#include +#endif + +/******************************************************************************/ +/* INCLUDE FILES */ +/******************************************************************************/ +#include +#include +#include +#include +#include +#include // for unlink +#include + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmAstrometryObjects.h" + +/********************* CheckDataStatus functions *****************************/ + +bool pmAstromCheckDataStatusForView (const pmFPAview *view, pmFPAfile *file) { + + pmFPA *fpa = file->fpa; + + if (view->chip == -1) { + bool exists = pmAstromCheckDataStatusForFPA (fpa); + return exists; + } + if (view->chip >= fpa->chips->n) { + psError(PS_ERR_IO, true, "Requested chip == %d >= fpa->chips->n == %ld", view->chip, fpa->chips->n); + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + bool exists = pmAstromCheckDataStatusForChip (chip); + return exists; + } + if (view->cell >= chip->cells->n) { + psError(PS_ERR_IO, true, "Requested cell == %d >= chip->cells->n == %ld", view->cell, chip->cells->n); + return false; + } + psError(PS_ERR_IO, false, "Astrometry only valid at the chip level"); + return false; +} + +bool pmAstromCheckDataStatusForFPA (const pmFPA *fpa) { + + if (!fpa->toTPA) return false; + if (!fpa->fromTPA) return false; + if (!fpa->toSky) return false; + + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + if (!chip) continue; + if (pmAstromCheckDataStatusForChip (chip)) return true; + } + return false; +} + +bool pmAstromCheckDataStatusForChip (const pmChip *chip) { + + if (!chip->toFPA) return false; + if (!chip->fromFPA) return false; // XXX not strictly needed? + return true; +} + +/********************* WritePHU function *****************************/ + +bool pmAstromWritePHU (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { + + // not needed if already written + if (file->wrote_phu) return true; + + // not needed if not FPA + // XXX this prevents us from defining a SPLIT/MEF CMF file... + if (file->fileLevel != PM_FPA_LEVEL_FPA) return true; + + // find the FPA phu + // XXX do we need to supply a header from the image? + pmHDU *phu = pmFPAviewThisPHU (view, file->fpa); + + // if there is no PHU, this is a single header+image (extension-less) file. This could be + // the case for an input SPLIT set of files being written out as a MEF. if there is a PHU, + // write it out as a 'blank' + psMetadata *outhead = psMetadataAlloc(); + if (phu) { + psMetadataCopy (outhead, phu->header); + } else { + pmConfigConformHeader (outhead, file->format); + + // XXX probably need other info to conform the header + psMetadata *fileData = psMetadataLookupMetadata(NULL, file->format, "FILE"); // File information + const char *fpaNameHdr = psMetadataLookupStr(NULL, fileData, "FPA.NAME"); + if (fpaNameHdr && strlen(fpaNameHdr) > 0) { + const char *fpaName = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.NAME"); + psMetadataAddStr(outhead, PS_LIST_TAIL, fpaNameHdr, PS_META_REPLACE, "FPA name", fpaName); + } + } + + psMetadataAddBool (outhead, PS_LIST_TAIL, "EXTEND", PS_META_REPLACE, "this file has extensions", true); + psFitsWriteBlank (file->fits, outhead, ""); + file->wrote_phu = true; + + psTrace ("pmFPAfile", 5, "wrote phu %s (type: %d)\n", file->filename, file->type); + psFree (outhead); + + return true; +} + +/********************* Write Data functions *****************************/ + +bool pmAstromWriteForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { + + pmFPA *fpa = file->fpa; + + if (view->chip == -1) { + if (!pmAstromWriteFPA (fpa, view, file, config)) { + psError(PS_ERR_IO, false, "Failed to write Astrometry for fpa"); + return false; + } + return true; + } + if (view->chip >= fpa->chips->n) { + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + if (!pmAstromWriteChip (chip, view, file, config)) { + psError(PS_ERR_IO, false, "Failed to write Astrometry for chip"); + return false; + } + return true; + } + + psError(PS_ERR_IO, false, "Astrometry must be written at the chip level"); + return false; +} + +// write out all chip-level Astrometry data for this FPA +bool pmAstromWriteFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +{ + pmFPAview *thisView = pmFPAviewAlloc (view->nRows); + *thisView = *view; + + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + thisView->chip = i; + if (!pmAstromWriteChip (chip, thisView, file, config)) { + psError(PS_ERR_IO, false, "Failed to write Astrometry for %dth chip", i); + psFree(thisView); + return false; + } + } + psFree(thisView); + return true; +} + +// write out chip-level Astrometry data for this chip +// 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?) +bool pmAstromWriteChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +{ + // 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; + } + extname = pmFPAfileNameFromRule (rule, file, view); + } + + // write the chip elements in the following form: + // chipID, ...? I forget the form: look at ICD + + return true; +} + +/********************* Read Data functions *****************************/ + +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? + 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) +{ + // 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; + } + extname = pmFPAfileNameFromRule (rule, file, view); + } + + // read the chip elements in the following form: + // chipID, ...? I forget the form: look at ICD + + return true; +} + Index: /trunk/psastro/src/psastroEnforceChip.c =================================================================== --- /trunk/psastro/src/psastroEnforceChip.c (revision 15357) +++ /trunk/psastro/src/psastroEnforceChip.c (revision 15357) @@ -0,0 +1,75 @@ +# include "psastroInternal.h" + +bool psastroEnforceChips (pmFPA *fpa, psMetadata *recipe) { + + pmChip *chip = NULL; + // identify reference astrometry table, stored as an fpa? better format? + + pmFPAfile *refAstrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.REF.ASTROM"); + + pmFPAview *view = pmFPAviewAlloc (0); + + // physical pixel scale in microns per pixel + double pixelTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.TOLERANCE"); + if (!status) { + psError(PS_ERR_IO, true, "Failed to lookup pixel tolerance"); + return false; + } + double angleTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.ANGLE.TOLERANCE"); + if (!status) { + psError(PS_ERR_IO, true, "Failed to lookup angle tolerance"); + return false; + } + + // loop over all chips + while ((obsChip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process); + if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } + + // set the chip astrometry using the refAstrom file + pmChip *refChip = pmFPAviewThisChip (view, refAstrom->fpa); + + // bad Astrometry test: ref pixel or angle outside nominal + + psPlane refPixel = {0.0, 0.0, 0.0, 0.0}; + psPlane obsCoord, refCoord; + + // find location of 0,0 pixel in focal plane coords for this chip + psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel); + + // find location of 0,0 pixel in focal plane coords for ref chip + psPlaneTransformApply (&refCoord, refChip->toFPA, &refPixel); + + psPlane offPixel = {0.0, 0.0, 0.0, 0.0}; + psPlane obsOffPt, refOffPt; + + // find location of 0,0 pixel in focal plane coords for this chip + psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel); + + // find location of 0,0 pixel in focal plane coords for ref chip + psPlaneTransformApply (&refOffPt, refChip->toFPA, &offPixel); + + double obsAngle = atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x); + + badAstrom = false; + badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol; + badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol; + badAstrom |= fabs(obsAngle - refAngle) > angleTol; + + if (!badAstrom) continue; + + psFree (obsChip->toFPA); + psFree (obsChip->fromFPA); + + // XXX can we do a straight copy? I think so, but we may need to apply a scaling factor based + // on the toSky or toTPA values + obsChip->toFPA = psMemIncrRefCounter (refChip->toFPA); + obsChip->fromFPA = psMemIncrRefCounter (refChip->fromFPA); + } + return true; +} + +// XXX for this function to work, I need to: +// 1 *) define badAstrom (ref pixel too far from nominal?) (ref angle too far off?) +// 2) load the refAstrom fpa in the pmFPAfileIO stage +// 3) allow for this operation to be optional Index: /trunk/psphot/src/psphotRadialProfile.c =================================================================== --- /trunk/psphot/src/psphotRadialProfile.c (revision 15357) +++ /trunk/psphot/src/psphotRadialProfile.c (revision 15357) @@ -0,0 +1,45 @@ +# include "psphotInternal.h" + +bool psphotRadialProfile (pmSource *source, psMetadata *recipe, psMaskType maskVal) { + + // allocate pmSourceExtendedParameters, if not already defined + if (!source->extpars) { + source->extpars = pmSourceExtendedParametersAlloc (); + } + + if (!source->extpars->profile) { + source->extpars->profile = pmSourceRadialProfile (); + } + + int nPts = source->pixels->numRows * source->pixels->numCols; + source->extpars->profile->radius = psVectorAllocEmpty (nPts, PS_TYPE_F32); + source->extpars->profile->flux = psVectorAllocEmpty (nPts, PS_TYPE_F32); + source->extpars->profile->weight = psVectorAllocEmpty (nPts, PS_TYPE_F32); + + psVector *radius = source->extpars->profile->radius; + psVector *flux = source->extpars->profile->flux; + psVector *weight = source->extpars->profile->radius; + + // XXX use the extended source model here for Xo, Yo? + // XXX define a radius scaled to the elliptical contour? + + int n = 0; + float Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0; + float Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + if (source->maskObj->data.U8[iy][ix]) continue; + radius->data.F32[n] = hypot (ix - Xo, iy - Yo) ; + flux->data.F32[n] = source->pixels->data.F32[iy][ix]; + weight->data.F32[n] = source->weight->data.F32[iy][ix]; + n++; + } + } + radius->n = n; + weight->n = n; + flux->n = n; + + SortVectorsByRadius (radius, flux, weight); + + return true; +}