Index: /branches/eam_branch_20070830/psModules/src/config/pmConfigCamera.c =================================================================== --- /branches/eam_branch_20070830/psModules/src/config/pmConfigCamera.c (revision 14729) +++ /branches/eam_branch_20070830/psModules/src/config/pmConfigCamera.c (revision 14729) @@ -0,0 +1,849 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include /* for strn?casecmp */ +#include + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmVersion.h" +#include "pmConcepts.h" +#include "pmConfigCamera.h" + +#define TABLE_OF_CONTENTS "CONTENTS" // Name for camera format metadata containing the contents +#define CHIP_TYPES "CHIPS" // Name for camera format metadata containing the chip types +#define CELL_TYPES "CELLS" // Name for camera format metadata containing the cell types + +// local helper functions defined below +static void removeCellConceptsSources(psMetadata *source); +static void removeChipConceptsSources(psMetadata *source); + +// Generate the skycell version of a named camera configuration +bool pmConfigCameraSkycellVersion(psMetadata *site, // The site configuration + const char *name // Name of the un-mosaicked camera + ) +{ + PS_ASSERT_METADATA_NON_NULL(site, false); + PS_ASSERT_STRING_NON_EMPTY(name, false); + + bool mdok; // Status of MD lookup + psMetadata *cameras = psMetadataLookupMetadata(&mdok, site, "CAMERAS"); // List of cameras + if (!mdok || !cameras) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find CAMERAS in the site configuration.\n"); + return false; + } + if (!pmConfigGenerateSkycellVersion(cameras, cameras, name, site)) { + psError(PS_ERR_UNKNOWN, true, "Failed to build skycell camera description for %s\n", name); + return false; + } + return true; +} + + +bool pmConfigCameraSkycellVersionsAll(psMetadata *site) +{ + PS_ASSERT_METADATA_NON_NULL(site, false); + + bool mdok; // Status of MD lookup + psMetadata *cameras = psMetadataLookupMetadata(&mdok, site, "CAMERAS"); // List of cameras + if (!mdok || !cameras) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find CAMERAS in the site configuration.\n"); + return false; + } + + psMetadataIterator *camerasIter = psMetadataIteratorAlloc(cameras, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *camerasItem = NULL; // Item from iteration + psMetadata *new = psMetadataAlloc();// New cameras to add + while ((camerasItem = psMetadataGetAndIncrement(camerasIter))) { + assert(camerasItem->type == PS_DATA_METADATA); // Only metadata are allowed here! + if (!pmConfigGenerateSkycellVersion(cameras, new, camerasItem->name, site)) { + psError(PS_ERR_UNKNOWN, true, "Failed to build skycell camera description for %s\n", + camerasItem->name); + return false; + } + } + psFree(camerasIter); + + // Now put the new cameras at the top of the list of cameras, so they get recognised first + // Note: going from the top, and putting everything to the top as we get there, so that the last one on + // goes to the top. This preserves the original order of the cameras, putting the skycell versions + // before the originals. + camerasIter = psMetadataIteratorAlloc(new, PS_LIST_HEAD, NULL); // Iterator + while ((camerasItem = psMetadataGetAndIncrement(camerasIter))) { + psMetadataAddItem(cameras, camerasItem, PS_LIST_HEAD, 0); + } + psFree(camerasIter); + psFree(new); + + return true; +} + +// Don't update these skycell concepts; last one MUST be 0 (i.e., NULL). +const static char *skycellConceptsCell[] = { "CELL.BIASSEC", "CELL.TRIMSEC", "CELL.READDIR", "CELL.XPARITY", + "CELL.YPARITY", "CELL.XWINDOW", "CELL.YWINDOW", "CELL.X0", + "CELL.Y0", "CELL.XSIZE", "CELL.YSIZE", 0 }; +const static char *skycellConceptsChip[] = { "CHIP.XPARITY", "CHIP.YPARITY", "CHIP.XSIZE", "CHIP.YSIZE", 0 }; +const static char *skycellConceptsFPA[] = { 0 }; + +// What do we call the skycell concept in the FITS header? +static const char *skycellConceptName(const char *name, // Name of concept + const char **concepts, // List of concepts NOT to update + const psMetadata *site // Site configuration + ) +{ + for (int i = 0; concepts[i]; i++) { + if (strcmp(name, concepts[i]) == 0) { + return NULL; + } + } + + if (!site) { + return name; + } + bool mdok; // Status of MD lookup + psMetadata *skycells = psMetadataLookupMetadata(&mdok, site, "SKYCELLS"); // Skycell concept headers + if (!skycells) { + return name; + } + const char *keyword = psMetadataLookupStr(&mdok, skycells, name); // Keyword to use for this concept + if (!mdok || !keyword || strlen(keyword) == 0) { + return name; + } + return keyword; +} + + +// Generate a skycell version of a camera configuration +bool pmConfigGenerateSkycellVersion(psMetadata *oldCameras, // Old list of camera configurations + psMetadata *newCameras, // New list of camera configurations + const char *name, // Name of original camera configuration + const psMetadata *site // Site configuration + ) +{ + assert(oldCameras); + assert(newCameras); + assert(name); + + // See if the old one is there + psMetadata *camera = psMetadataLookupMetadata(NULL, oldCameras, name); // The camera configuration + if (!camera) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find camera to be skycelled in camera list."); + return false; + } + + // See if the new one is already there + psString newName = NULL; // Name of skycelled camera + psStringAppend(&newName, "_%s-SKYCELL", name); + if (psMetadataLookup(oldCameras, newName)) { + return true; + } + + psMetadata *new = psMetadataCopy(NULL, camera); // Copy of the camera description + + // Fix the FPA description to contain a single chip with single cell + { + psMetadata *fpa = psMetadataAlloc();// The FPA description + psMetadataAddStr(fpa, PS_LIST_HEAD, "SkyChip", 0, "Single chip with single cell", "SkyCell"); + psMetadataAddMetadata(new, PS_LIST_TAIL, "FPA", PS_META_REPLACE, "Description of FPA hierarchy", fpa); + psFree(fpa); + } + + // Clear out the formats, replace them with the One True Format + psMetadata *formats = psMetadataLookupMetadata(NULL, new, "FORMATS"); // The list of formats + if (!formats) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find FORMATS within camera configuration."); + psFree(new); + return false; + } + while (psListLength(formats->list) > 0) { + psMetadataRemoveIndex(formats, PS_LIST_HEAD); + } + psMetadata *format = psMetadataAlloc(); // The One True Format + + { + psMetadata *rule = psMetadataAlloc(); // The RULE --- how to recognise the camera + psMetadataAddStr(rule, PS_LIST_TAIL, "PSCAMERA", 0, "Camera name", name); + psMetadataAddStr(rule, PS_LIST_TAIL, "PSFORMAT", 0, "Camera format", "SKYCELL"); + psString version = psModulesVersion(); + psMetadataAddStr(rule, PS_LIST_TAIL, "ORIGIN", 0, "File origin", version); + psFree(version); + psMetadataAddMetadata(format, PS_LIST_TAIL, "RULE", 0, "How to recognise this type of file", rule); + psFree(rule); + } + + { + psMetadata *file = psMetadataAlloc(); // The FILE --- how to read the data + psMetadataAddStr(file, PS_LIST_TAIL, "PHU", 0, "What level the FITS file represents", "FPA"); + psMetadataAddStr(file, PS_LIST_TAIL, "EXTENSIONS", 0, "What level the extensions represent", "NONE"); + psMetadataAddStr(file, PS_LIST_TAIL, "FPA.NAME", 0, "PHU keyword for unique identifier", "FPA.NAME"); + psMetadataAddMetadata(format, PS_LIST_TAIL, "FILE", 0, "How to read this type of file", file); + psFree(file); + } + + psMetadataAddStr(format, PS_LIST_TAIL, "CONTENTS", 0, "What's in this type of file", + "SkyChip:SkyCell:_skycell"); + + { + psMetadata *cells = psMetadataAlloc(); // The CELLS --- how to read the cells + psMetadata *skycell = psMetadataAlloc(); // How to read the skycell + psMetadataAddStr(skycell, PS_LIST_TAIL, "CELL.TRIMSEC", 0, "Trim section", "CELL.TRIMSEC"); + psMetadataAddStr(skycell, PS_LIST_TAIL, "CELL.BIASSEC", 0, "Bias section", "CELL.BIASSEC"); + psMetadataAddStr(skycell, PS_LIST_TAIL, "CELL.TRIMSEC.SOURCE", 0, "Source for trim section", + "HEADER"); + psMetadataAddStr(skycell, PS_LIST_TAIL, "CELL.BIASSEC.SOURCE", 0, "Source for bias section", + "HEADER"); + psMetadataAddMetadata(cells, PS_LIST_TAIL, "_skycell", 0, "Skycell specification", skycell); + psFree(skycell); + psMetadataAddMetadata(format, PS_LIST_TAIL, "CELLS", 0, "How to read the cells", cells); + psFree(cells); + } + + // Stuffing all concepts into the header, by their PS concept name (e.g., "FPA.AIRMASS"). + // (HIERARCH will take care of the long names, implemented in psLib.) + // Some people may not like this, but it's quick and easy and will do for now. + // An alternative may be provided later. + { + psMetadata *translation = psMetadataAlloc(); // The TRANSLATION --- how to read the FITS headers + + psList *concepts; // List of concepts for each level + psListIterator *iter; // Iterator for concepts + psString name; // Concept name, from iteration + + concepts = pmConceptsList(PM_FPA_LEVEL_FPA); // FPA-level concepts + iter = psListIteratorAlloc(concepts, PS_LIST_HEAD, false); + while ((name = psListGetAndIncrement(iter))) { + const char *new = skycellConceptName(name, skycellConceptsFPA, site); // Name for skycell + if (new) { + psMetadataAddStr(translation, PS_LIST_TAIL, name, 0, NULL, new); + } + } + psFree(iter); + psFree(concepts); + + concepts = pmConceptsList(PM_FPA_LEVEL_CHIP); + iter = psListIteratorAlloc(concepts, PS_LIST_HEAD, false); + while ((name = psListGetAndIncrement(iter))) { + const char *new = skycellConceptName(name, skycellConceptsChip, site); // Name for skycell + if (new) { + psMetadataAddStr(translation, PS_LIST_TAIL, name, 0, NULL, new); + } + } + psFree(iter); + psFree(concepts); + + concepts = pmConceptsList(PM_FPA_LEVEL_CELL); + iter = psListIteratorAlloc(concepts, PS_LIST_HEAD, false); + while ((name = psListGetAndIncrement(iter))) { + const char *new = skycellConceptName(name, skycellConceptsCell, site); // Name for skycell + if (new) { + psMetadataAddStr(translation, PS_LIST_TAIL, name, 0, NULL, new); + } + } + psFree(iter); + psFree(concepts); + + psMetadataAddStr(translation, PS_LIST_TAIL, "CELL.XSIZE", 0, NULL, "NAXIS1"); + psMetadataAddStr(translation, PS_LIST_TAIL, "CELL.YSIZE", 0, NULL, "NAXIS2"); + psMetadataAddStr(translation, PS_LIST_TAIL, "CELL.XWINDOW", 0, NULL, "NAXIS1"); + psMetadataAddStr(translation, PS_LIST_TAIL, "CELL.YWINDOW", 0, NULL, "NAXIS2"); + psMetadataAddStr(translation, PS_LIST_TAIL, "CHIP.XSIZE", 0, NULL, "NAXIS1"); + psMetadataAddStr(translation, PS_LIST_TAIL, "CHIP.YSIZE", 0, NULL, "NAXIS2"); + + psMetadataAddMetadata(format, PS_LIST_TAIL, "TRANSLATION", 0, "How to translate the FITS headers", + translation); + psFree(translation); + } + + { + psMetadata *defaults = psMetadataAlloc(); // Default values for concepts + + psMetadataAddS32(defaults, PS_LIST_TAIL, "CELL.XPARITY", 0, NULL, 1); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CELL.YPARITY", 0, NULL, 1); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CHIP.XPARITY", 0, NULL, 1); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CHIP.YPARITY", 0, NULL, 1); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CELL.Y0", 0, NULL, 0); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CELL.X0", 0, NULL, 0); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CHIP.Y0", 0, NULL, 0); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CHIP.X0", 0, NULL, 0); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CELL.READDIR", 0, "Read direction (rows)", 1); + + psMetadataAddMetadata(format, PS_LIST_TAIL, "DEFAULTS", 0, "Default values for concepts", defaults); + psFree(defaults); + + } + + { + psMetadata *database = psMetadataAlloc(); // Database values for concepts + psMetadataAddMetadata(format, PS_LIST_TAIL, "DATABASE", 0, "Database values for concepts", database); + psFree(database); + } + + { + psMetadata *conceptFormats = psMetadataAlloc(); // Format peculiarities for various concepts + // These are the only essential formats + psMetadataAddStr(conceptFormats, PS_LIST_TAIL, "FPA.RA", 0, "Units for RA", "HOURS"); + psMetadataAddStr(conceptFormats, PS_LIST_TAIL, "FPA.DEC", 0, "Units for RA", "DEGREES"); + psMetadataAddStr(conceptFormats, PS_LIST_TAIL, "FPA.TIME", 0, "Format for time", "MJD"); + psMetadataAddStr(conceptFormats, PS_LIST_TAIL, "CELL.TIME", 0, "Format for time", "MJD"); + + psMetadataAddMetadata(format, PS_LIST_TAIL, "FORMATS", 0, "Formats for various concepts", + conceptFormats); + psFree(conceptFormats); + } + + psMetadataAddMetadata(formats, PS_LIST_TAIL, "SKYCELL", 0, "The One True Format for skycells", format); + psFree(format); + + // New camera MUST go to the head of the metadata, so that it will be recognised first + // The old camera doesn't contain the PSCAMERA and PSFORMAT headers, so it will match anything! + psMetadataAddMetadata(newCameras, PS_LIST_HEAD, newName, PS_META_REPLACE, + "Automatically generated", new); + psFree(newName); + psFree(new); + + return true; +} + +// Generate the Chip and FPA mosaicked version of a named camera configuration +bool pmConfigCameraMosaickedVersions(psMetadata *site, // The site configuration + const char *name // Name of the un-mosaicked camera + ) +{ + PS_ASSERT_METADATA_NON_NULL(site, false); + PS_ASSERT_STRING_NON_EMPTY(name, false); + + bool mdok; // Status of MD lookup + psMetadata *cameras = psMetadataLookupMetadata(&mdok, site, "CAMERAS"); // List of cameras + if (!mdok || !cameras) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find CAMERAS in the site configuration.\n"); + return false; + } + if (!pmConfigGenerateMosaickedVersion(cameras, cameras, name, PM_FPA_LEVEL_CHIP)) { + psError(PS_ERR_UNKNOWN, true, "Failed to build Chip mosaic camera description for %s\n", name); + return false; + } + if (!pmConfigGenerateMosaickedVersion(cameras, cameras, name, PM_FPA_LEVEL_FPA)) { + psError(PS_ERR_UNKNOWN, true, "Failed to build FPA mosaic camera description for %s\n", name); + return false; + } + return true; +} + +// the operation putting the new entries first is now implemented in pmConfigGenerateMosaickedVersion +// Generate the Chip and FPA mosaicked version of a named camera configuration +bool pmConfigCameraMosaickedVersionsAll(psMetadata *site) +{ + PS_ASSERT_METADATA_NON_NULL(site, false); + + bool mdok; // Status of MD lookup + psMetadata *cameras = psMetadataLookupMetadata(&mdok, site, "CAMERAS"); // List of cameras + if (!mdok || !cameras) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find CAMERAS in the site configuration.\n"); + return false; + } + + psMetadataIterator *camerasIter = psMetadataIteratorAlloc(cameras, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *camerasItem = NULL; // Item from iteration + psMetadata *new = psMetadataAlloc();// New cameras to add + while ((camerasItem = psMetadataGetAndIncrement(camerasIter))) { + assert(camerasItem->type == PS_DATA_METADATA); // Only metadata are allowed here! + if (!pmConfigGenerateMosaickedVersion(cameras, new, camerasItem->name, PM_FPA_LEVEL_CHIP)) { + psError(PS_ERR_UNKNOWN, true, "Failed to build Chip mosaic camera description for %s\n", camerasItem->name); + return false; + } + if (!pmConfigGenerateMosaickedVersion(cameras, new, camerasItem->name, PM_FPA_LEVEL_FPA)) { + psError(PS_ERR_UNKNOWN, true, "Failed to build FPA mosaic camera description for %s\n", camerasItem->name); + return false; + } + } + psFree(camerasIter); + + // Now put the new cameras at the top of the list of cameras, so they get recognised first + // Note: going from the top, and putting everything to the top as we get there, so that the last one on + // goes to the top. This preserves the original order of the cameras, putting the mosaicked versions + // before the originals. + camerasIter = psMetadataIteratorAlloc(new, PS_LIST_HEAD, NULL); // Iterator + while ((camerasItem = psMetadataGetAndIncrement(camerasIter))) { + psMetadataAddItem(cameras, camerasItem, PS_LIST_HEAD, 0); + } + psFree(camerasIter); + psFree(new); + + return true; +} + +// Generate a mosaicked version of a camera configuration +bool pmConfigGenerateMosaickedVersion(psMetadata *oldCameras, // Old list of camera configurations + psMetadata *newCameras, // New list of camera configurations + const char *name, // Name of original camera configuration + pmFPALevel mosaicLevel // Level to which we are mosaicking + ) +{ + assert(oldCameras); + assert(newCameras); + assert(name); + assert(mosaicLevel == PM_FPA_LEVEL_CHIP || mosaicLevel == PM_FPA_LEVEL_FPA); + + // Don't generate a mosaicked version of a skycell + if (name[0] == '_' && strlen(name) > 9 && strcmp(name + strlen(name) - 8, "-SKYCELL") == 0) { + return true; + } + + // See if the old one is there + psMetadata *camera = psMetadataLookupMetadata(NULL, oldCameras, name); // The camera configuration + if (!camera) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find camera to be mosaicked in camera list."); + return false; + } + + // See if the new one is already there + psString newName = NULL; // Name of mosaicked camera + psStringAppend(&newName, "_%s-%s", name, mosaicLevel == PM_FPA_LEVEL_CHIP ? "CHIP" : "FPA"); + if (psMetadataLookup(oldCameras, newName)) { + return true; + } + + psMetadata *new = psMetadataCopy(NULL, camera); // Copy of the camera description + bool mdok; // Status of MD lookups + + // ** Fix up the contents of the FPA description to match the mosaicked camera ** + // select the FPA description + psMetadata *fpa = psMetadataLookupMetadata(NULL, new, "FPA"); // FPA in the camera configuration + if (!fpa) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find FPA within camera configuration."); + psFree(new); + return false; + } + switch (mosaicLevel) { + // For CHIP mosaic, replace the contents of each chip with a single cell + case PM_FPA_LEVEL_CHIP: { + psMetadataIterator *fpaIter = psMetadataIteratorAlloc(fpa, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *fpaItem = NULL; // Item from iteration + while ((fpaItem = psMetadataGetAndIncrement(fpaIter))) { + if (fpaItem->type != PS_DATA_STRING) { + psError(PS_ERR_UNKNOWN, true, "Element %s within FPA in camera configuration is not of type STR.", + fpaItem->name); + psFree(new); + return false; + } + + psFree(fpaItem->data.str); + fpaItem->data.str = psStringCopy("MosaickedCell"); + psFree(fpaItem->comment); + fpaItem->comment = psStringCopy("Mosaicked cell; automatically generated"); + } + psFree(fpaIter); + break; + } + // For FPA mosaic, replace the contents of the FPA with a single chip containing a single cell + case PM_FPA_LEVEL_FPA: { + while (psListLength(fpa->list) > 0) { + psMetadataRemoveIndex(fpa, PS_LIST_TAIL); + } + + psMetadataAddStr(fpa, PS_LIST_HEAD, "MosaickedChip", 0, + "Mosaicked chip with mosaicked cell; automatically generated", + "MosaickedCell"); + break; + } + default: + psAbort("Should never get here.\n"); + } + + // ** Update the camera formats : add a new (mosaicked) format for each existing camera format ** + // select the list of all camera formats + psMetadata *formats = psMetadataLookupMetadata(NULL, new, "FORMATS"); // FORMATS in the configuration + assert(formats); // It had better be there --- we've already read them in + // loop over each of the formats + psMetadataIterator *formatsIter = psMetadataIteratorAlloc(formats, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *formatsItem = NULL; // Item from iteration + while ((formatsItem = psMetadataGetAndIncrement(formatsIter))) { + assert(formatsItem->type == PS_DATA_METADATA); // We should have read it by now! + psMetadata *format = formatsItem->data.V; // The camera format + + // Add a new RULE which uniquely describes the mosaicked format. this is needed so + // that when a mosaic is written to a FITS file, it can be recognised again when read. + psMetadata *rule = psMetadataLookupMetadata(NULL, format, "RULE"); // Way to identify format from PHU + if (!rule) { + // a camera format without a rule is not allowed. + psError(PS_ERR_UNKNOWN, false, "Camera format %s has no RULE", formatsItem->name); + return false; + } + + // the new rule is supplemented by the mosaicLevel + switch (mosaicLevel) { + case PM_FPA_LEVEL_CHIP: + psMetadataAddStr(rule, PS_LIST_TAIL, "PSMOSAIC", 0, "Mosaicked level", "CHIP"); + break; + case PM_FPA_LEVEL_FPA: + psMetadataAddStr(rule, PS_LIST_TAIL, "PSMOSAIC", 0, "Mosaicked level", "FPA"); + break; + default: + psAbort("Should never get here.\n"); + } + + // Fix the FILE information: need to fix the levels for the PHU and EXTENSIONS. + // both of these elements are required in the format; we raise an error if they are not found + // If EXTENSIONS is NONE, then we need to change the CONTENT specifier to point to the chip name. + psMetadata *file = psMetadataLookupMetadata(NULL, format, "FILE"); // File information + if (!file) { + psError(PS_ERR_UNKNOWN, false, "Camera format %s has no FILE", formatsItem->name); + return false; + } + psMetadataItem *phuItem = psMetadataLookup(file, "PHU"); // PHU level + if (!phuItem || phuItem->type != PS_DATA_STRING) { + psError(PS_ERR_UNKNOWN, false, "Camera format %s is missing PHU in the FILE information", formatsItem->name); + return false; + } + psMetadataItem *extensionsItem = psMetadataLookup(file, "EXTENSIONS"); // Extensions level + if (!extensionsItem || extensionsItem->type != PS_DATA_STRING) { + psError(PS_ERR_UNKNOWN, false, "Camera format %s is missing EXTENSIONS in the FILE information", formatsItem->name); + return false; + } + + // mosaicLevel == CHIP: + // Case PHU EXTENSIONS Modifications + // ==== === ========== =========== + // 1. FPA CHIP NONE + // 2. FPA CELL EXT->CHIP + // 3. FPA NONE NONE + // 4. CHIP CELL EXT->NONE + // 5. CHIP NONE NONE + // 6. CELL NONE PHU->CHIP + // possible outcomes: + // FPA CHIP + // FPA NONE + // CHIP NONE + + // mosaicLevel == FPA: + // Case PHU EXTENSIONS Modifications + // ==== === ========== =========== + // 1. FPA CHIP EXT->NONE + // 2. FPA CELL EXT->NONE + // 3. FPA NONE NONE + // 4. CHIP CELL PHU->FPA, EXT->NONE + // 5. CHIP NONE PHU->FPA + // 6. CELL NONE PHU->FPA + // possible outcomes: + // FPA NONE + + // modify the values of phuItem and extensionsItem + switch (mosaicLevel) { + case PM_FPA_LEVEL_CHIP: + if (!strcasecmp(phuItem->data.str, "FPA") && !strcasecmp(extensionsItem->data.str, "CHIP")) { + break; + } + if (!strcasecmp(phuItem->data.str, "FPA") && !strcasecmp(extensionsItem->data.str, "CELL")) { + psFree(extensionsItem->data.str); + extensionsItem->data.str = psStringCopy("CHIP"); + break; + } + if (!strcasecmp(phuItem->data.str, "FPA") && !strcasecmp(extensionsItem->data.str, "NONE")) { + break; + } + if (!strcasecmp(phuItem->data.str, "CHIP") && !strcasecmp(extensionsItem->data.str, "CELL")) { + psFree(extensionsItem->data.str); + extensionsItem->data.str = psStringCopy("NONE"); + break; + } + if (!strcasecmp(phuItem->data.str, "CHIP") && !strcasecmp(extensionsItem->data.str, "NONE")) { + break; + } + if (!strcasecmp(phuItem->data.str, "CELL") && !strcasecmp(extensionsItem->data.str, "NONE")) { + psFree(phuItem->data.str); + phuItem->data.str = psStringCopy("CHIP"); + break; + } + psAbort ("should not reach here"); + + case PM_FPA_LEVEL_FPA: + if (!strcasecmp(phuItem->data.str, "FPA") && !strcasecmp(extensionsItem->data.str, "CHIP")) { + psFree(extensionsItem->data.str); + extensionsItem->data.str = psStringCopy("NONE"); + break; + } + if (!strcasecmp(phuItem->data.str, "FPA") && !strcasecmp(extensionsItem->data.str, "CELL")) { + psFree(extensionsItem->data.str); + extensionsItem->data.str = psStringCopy("NONE"); + break; + } + if (!strcasecmp(phuItem->data.str, "FPA") && !strcasecmp(extensionsItem->data.str, "NONE")) { + break; + } + if (!strcasecmp(phuItem->data.str, "CHIP") && !strcasecmp(extensionsItem->data.str, "CELL")) { + psFree(phuItem->data.str); + phuItem->data.str = psStringCopy("FPA"); + psFree(extensionsItem->data.str); + extensionsItem->data.str = psStringCopy("NONE"); + break; + } + if (!strcasecmp(phuItem->data.str, "CHIP") && !strcasecmp(extensionsItem->data.str, "NONE")) { + psFree(phuItem->data.str); + phuItem->data.str = psStringCopy("FPA"); + break; + } + if (!strcasecmp(phuItem->data.str, "CELL") && !strcasecmp(extensionsItem->data.str, "NONE")) { + psFree(phuItem->data.str); + phuItem->data.str = psStringCopy("FPA"); + break; + } + psAbort ("should not reach here"); + + default: + psAbort("Should never get here.\n"); + } + +#if 0 + // XXXX when do I need to adjust the value of CONTENT, CELL.NAME, CHIP.NAME? + + // Don't need a "CONTENT" to identify the content! + if (psMetadataLookup(file, "CONTENT")) { + psMetadataRemoveKey(file, "CONTENT"); + } + psMetadataAddStr(file, PS_LIST_TAIL, "CONTENT", PS_META_REPLACE, "Key to CONTENTS menu", + "{CHIP.NAME}"); + // Don't need CELL.NAME for chip-mosaicked camera + if (psMetadataLookup(file, "CELL.NAME")) { + psMetadataRemoveKey(file, "CELL.NAME"); + } + // Don't need CHIP.NAME or CELL.NAME for fpa-mosaicked camera + if (psMetadataLookup(file, "CHIP.NAME")) { + psMetadataRemoveKey(file, "CHIP.NAME"); + } + if (psMetadataLookup(file, "CELL.NAME")) { + psMetadataRemoveKey(file, "CELL.NAME"); + } +#endif + + // Fix up the CONTENTS to contain only the mosaicked cell for each chip + switch (mosaicLevel) { + case PM_FPA_LEVEL_FPA: + psMetadataAddStr(format, PS_LIST_TAIL, TABLE_OF_CONTENTS, PS_META_REPLACE, NULL, + "MosaickedChip:MosaickedCell:_mosaic"); + break; + case PM_FPA_LEVEL_CHIP: + if (!strcasecmp(phuItem->data.str, "FPA") && !strcasecmp(extensionsItem->data.str, "CHIP")) { + // ensure the value of CONTENT in the FILE section has the right value + psMetadataAddStr(file, PS_LIST_TAIL, "CONTENT", PS_META_REPLACE, "Key to CONTENTS menu", "{CHIP.NAME}"); + + // List the chipName:chipType for each chip. + psMetadata *contents = psMetadataAlloc(); // List of contents, with chipName:chipType + + psMetadataIterator *fpaIter = psMetadataIteratorAlloc(fpa, PS_LIST_HEAD, NULL); // Iteratr + psMetadataItem *fpaItem; // Item from iteration + while ((fpaItem = psMetadataGetAndIncrement(fpaIter))) { + assert (fpaItem->type == PS_DATA_STRING); + psString content = NULL; // Content to add + psStringAppend(&content, "%s:_mosaicChip ", fpaItem->name); + psMetadataAddStr(contents, PS_LIST_TAIL, fpaItem->name, 0, NULL, content); + psFree(content); + } + psFree(fpaIter); + psMetadataAddMetadata(format, PS_LIST_TAIL, TABLE_OF_CONTENTS, PS_META_REPLACE, + "List of contents", contents); + psFree(contents); + + psMetadata *chips = psMetadataAlloc(); // List of chip types, with cellName:cellType + psMetadataAddStr(chips, PS_LIST_TAIL, "_mosaicChip", 0, NULL, + "MosaickedCell:_mosaic"); + psMetadataAddMetadata(format, PS_LIST_TAIL, CHIP_TYPES, PS_META_REPLACE, + "List of chip types", chips); + psFree(chips); + break; + } + if (!strcasecmp(phuItem->data.str, "FPA") && !strcasecmp(extensionsItem->data.str, "NONE")) { + // List the contents on a single line + psString contentsLine = NULL; // Contents of the PHU + psMetadataIterator *fpaIter = psMetadataIteratorAlloc(fpa, PS_LIST_HEAD, NULL); // Iteratr + psMetadataItem *fpaItem; // Item from iteration + while ((fpaItem = psMetadataGetAndIncrement(fpaIter))) { + assert (fpaItem->type == PS_DATA_STRING); + psStringAppend(&contentsLine, "%s:MosaickedCell:_mosaic ", fpaItem->name); + } + psFree(fpaIter); + psMetadataAddStr(format, PS_LIST_TAIL, TABLE_OF_CONTENTS, PS_META_REPLACE, + NULL, contentsLine); + psFree(contentsLine); + break; + } + if (!strcasecmp(phuItem->data.str, "CHIP") && !strcasecmp(extensionsItem->data.str, "NONE")) { + // XXX recode this to match the structure above (FPA/CHIP)? + // select and remove the old contents + psMetadata *contents = psMetadataLookupMetadata(NULL, format, TABLE_OF_CONTENTS); // File contents + if (!contents) { + psError(PS_ERR_UNKNOWN, false, "Couldn't find %s in the camera format %s.\n", + TABLE_OF_CONTENTS, formatsItem->name); + return false; + } + while (psListLength(contents->list) > 0) { + psMetadataRemoveIndex(contents, PS_LIST_TAIL); + } + + // update with the new contents + psMetadataIterator *fpaIter = psMetadataIteratorAlloc(fpa, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *fpaItem; // Item from iteration + while ((fpaItem = psMetadataGetAndIncrement(fpaIter))) { + assert (fpaItem->type == PS_DATA_STRING); + psMetadataAddStr(contents, PS_LIST_TAIL, fpaItem->name, 0, NULL, "_mosaicChip"); + } + psFree(fpaIter); + + psMetadata *chips = psMetadataAlloc(); // List of chip types, with cellName:cellType + psMetadataAddStr(chips, PS_LIST_TAIL, "_mosaicChip", 0, NULL, + "MosaickedCell:_mosaic"); + psMetadataAddMetadata(format, PS_LIST_TAIL, CHIP_TYPES, PS_META_REPLACE, + "List of chip types", chips); + psFree(chips); + break; + } + default: + psAbort("Should never get here.\n"); + } + + // Fix the cell type + psMetadata *cells = psMetadataLookupMetadata(NULL, format, CELL_TYPES); // CELLS information + if (!cells) { + psError(PS_ERR_UNKNOWN, false, "Couldn't find CELLS of type METADATA in the camera format %s.\n", formatsItem->name); + return false; + } + psMetadata *cell = psMetadataAlloc(); // Cell information + psMetadataAddStr(cell, PS_LIST_TAIL, "CELL.TRIMSEC", 0, "Trim section", "TRIMSEC"); + psMetadataAddStr(cell, PS_LIST_TAIL, "CELL.BIASSEC", 0, "Bias section", "BIASSEC"); + psMetadataAddStr(cell, PS_LIST_TAIL, "CELL.TRIMSEC.SOURCE", 0, "Trim section source", "HEADER"); + psMetadataAddStr(cell, PS_LIST_TAIL, "CELL.BIASSEC.SOURCE", 0, "Bias section source", "HEADER"); + psMetadataAddMetadata(cells, PS_LIST_HEAD, "_mosaic", PS_META_REPLACE, "Mosaic cell information", cell); + psFree(cell); // Drop reference + + // Update the concepts, so that they are all stored in the FITS headers, under headers of the same + // name as the concept + psMetadata *database = psMetadataLookupMetadata(&mdok, format, "DATABASE"); // DATABASE concepts + psMetadata *defaults = psMetadataLookupMetadata(&mdok, format, "DEFAULTS"); // DEFAULTS concepts + if (!mdok || !defaults) { + psWarning("Couldn't find DEFAULTS of type METADATA in the camera format %s.\n", + formatsItem->name); + continue; + } + psMetadata *translation = psMetadataLookupMetadata(&mdok, format, "TRANSLATION"); // TRANSLATION info + if (!mdok || !translation) { + psWarning("Couldn't find TRANSLATION of type METADATA in the camera format %s.\n", + formatsItem->name); + continue; + } + + removeCellConceptsSources(translation); + removeCellConceptsSources(database); + removeCellConceptsSources(defaults); + + // Add in the positioning concepts + psMetadataAddS32(defaults, PS_LIST_TAIL, "CELL.XPARITY", 0, NULL, 1); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CELL.YPARITY", 0, NULL, 1); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CELL.X0", 0, NULL, 0); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CELL.Y0", 0, NULL, 0); + if (formats) { + if (psMetadataLookup(formats, "CELL.X0")) { + psMetadataRemoveKey(formats, "CELL.X0"); + } + if (psMetadataLookup(formats, "CELL.Y0")) { + psMetadataRemoveKey(formats, "CELL.Y0"); + } + } + + if (mosaicLevel == PM_FPA_LEVEL_FPA) { + removeChipConceptsSources(translation); + removeChipConceptsSources(database); + removeChipConceptsSources(defaults); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CHIP.XPARITY", 0, NULL, 1); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CHIP.YPARITY", 0, NULL, 1); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CHIP.X0", 0, NULL, 0); + psMetadataAddS32(defaults, PS_LIST_TAIL, "CHIP.Y0", 0, NULL, 0); + if (formats) { + if (psMetadataLookup(formats, "CHIP.X0")) { + psMetadataRemoveKey(formats, "CHIP.X0"); + } + if (psMetadataLookup(formats, "CHIP.Y0")) { + psMetadataRemoveKey(formats, "CHIP.Y0"); + } + } + } + + } + psFree(formatsIter); + + // New camera MUST go to the head of the metadata, so that it will be recognised first + // The old camera doesn't contain the PSMOSAIC header, so it will match anything! + psMetadataAddMetadata(newCameras, PS_LIST_HEAD, newName, PS_META_REPLACE, + "Automatically generated", new); + psFree(newName); + psFree(new); + + return true; +} + +/*** Helper Functions ***/ + +// Remove a concept from the list of sources. Need to check to see if it exists first, to avoid a warning. +static void removeConcept(psMetadata *source, // Source from which to remove concept + const char *concept // Concept name to remove + ) +{ + assert(source); + assert(concept && strlen(concept) > 0); + + if (psMetadataLookup(source, concept)) { + psMetadataRemoveKey(source, concept); + } + + return; +} + +// Remove certain concepts from the list of sources. These concepts are important in the mosaicking process, +// and are added explicitly to the defaults (elsewhere) so that the user can't get them wrong. +static void removeCellConceptsSources(psMetadata *source // Source for concepts + ) +{ + if (!source) { + return; + } + + removeConcept(source, "CELL.BIASSEC"); + removeConcept(source, "CELL.TRIMSEC"); + removeConcept(source, "CELL.XPARITY"); + removeConcept(source, "CELL.YPARITY"); + removeConcept(source, "CELL.X0"); + removeConcept(source, "CELL.Y0"); + + // For the sake of the defaults, include the .DEPEND + removeConcept(source, "CELL.XPARITY.DEPEND"); + removeConcept(source, "CELL.YPARITY.DEPEND"); + removeConcept(source, "CELL.X0.DEPEND"); + removeConcept(source, "CELL.Y0.DEPEND"); + + return; +} + +// Remove certain concepts from the list of sources. These concepts are important in the mosaicking process, +// and are added explicitly to the defaults (elsewhere) so that the user can't get them wrong. +static void removeChipConceptsSources(psMetadata *source // Source for concepts + ) +{ + if (!source) { + return; + } + + removeConcept(source, "CHIP.XPARITY"); + removeConcept(source, "CHIP.YPARITY"); + removeConcept(source, "CHIP.X0"); + removeConcept(source, "CHIP.Y0"); + + // For the sake of the defaults, include the .DEPEND + removeConcept(source, "CHIP.XPARITY.DEPEND"); + removeConcept(source, "CHIP.YPARITY.DEPEND"); + removeConcept(source, "CHIP.X0.DEPEND"); + removeConcept(source, "CHIP.Y0.DEPEND"); + + return; +} Index: /branches/eam_branch_20070830/psphot/doc/2Dmodels.txt =================================================================== --- /branches/eam_branch_20070830/psphot/doc/2Dmodels.txt (revision 14729) +++ /branches/eam_branch_20070830/psphot/doc/2Dmodels.txt (revision 14729) @@ -0,0 +1,29 @@ + +I am testing a new concept for the 2D representation of various values +measured by psphot. I have been using 2D polynomials to fit the +variations of, eg, the psf parameters as a function of position in the +images. this has been somewhat successful, but it has also been +fraught with a number of problems. the worst issue has been the fact +that the polynomial form is not physically motivated. in some cases, +eg the variation of the psf shape parameters, the true variations for +real distributions are not well represented by a polynomial. I have +done some ad hoc tricks to get the functional form to look roughly +like a polynomial, but this is really quite a hack. the other problem +is that it is difficult to decide if the polynomial does an acceptable +job representing the variations. + +My new concept for modeling 2D variations in some parameter is to +define a low-resolution 2D image to represent the variation. The +advantage is that the image may be define with whatever resolution may +be sampled by the input data, and the functional form need not be +constrained by the polynomial basis function. + +** is this a poor substitute for using a chebychev polynomial as a + better basis function? + +psphot parameters which vary by position: + +pmPSF.params +pmPSF.ApTrend +pmPSF.FluxScale (flux given normalization of 1.0) + Index: /branches/eam_branch_20070830/psphot/src/psphotImageMedian.c =================================================================== --- /branches/eam_branch_20070830/psphot/src/psphotImageMedian.c (revision 14729) +++ /branches/eam_branch_20070830/psphot/src/psphotImageMedian.c (revision 14729) @@ -0,0 +1,387 @@ +# include "psphotInternal.h" +static int npass = 0; + +// we have 4 possibilities: (INTERNAL or I/O file) and (exists or not) +// select model pixels (from output background model file, or create internal file) +static pmReadout *get_model_readout(const char *name, // name of internal/external file + const pmConfig *config, // configuration information + const pmFPAview *view, + pmFPA *inFPA, + const psImageBinning *binning) { + pmReadout *model = NULL; + + bool status = true; + pmFPAfile *file = psMetadataLookupPtr(&status, config->files, name); + if (file == NULL) { + // we are not using PSPHOT.BACKMDL as an I/O file: define an internal version + model = pmFPAfileDefineInternal (config->files, name, binning->nXruff, binning->nYruff, PS_TYPE_F32); + } else { + if (file->mode == PM_FPA_MODE_INTERNAL) { + // we are not using PSPHOT.BACKMDL as an I/O file: already defined above + model = file->readout; + } else { + // we are using PSPHOT.BACKMDL as an I/O file: select readout or create + model = pmFPAviewThisReadout (view, file->fpa); + if (model == NULL) { + // readout does not yet exist: create from input + // XXX we have an inconsistency in this calculation here and in pmFPACopy + // XXX use the psImageBinning functions to set the output image size + pmFPAfileCopyStructureView (file->fpa, inFPA, binning->nXbin, binning->nYbin, view); + model = pmFPAviewThisReadout (view, file->fpa); + PS_ASSERT (binning->nXruff == model->image->numCols, false); + PS_ASSERT (binning->nYruff == model->image->numRows, false); + } + } + } + + return model; +} + +// generate the median in NxN boxes, clipping heavily +// linear interpolation to generate full-scale model +bool psphotImageMedian (pmConfig *config, const pmFPAview *view, psMaskType maskVal) +{ + bool status = true; + static char *defaultStatsName = "FITTED_MEAN"; + + pmReadout *readout = NULL; + pmReadout *background = NULL; + pmReadout *backSub = NULL; + + psTimerStart ("psphot"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + + // user supplied seed, if available + unsigned long seed = psMetadataLookupS32 (&status, recipe, "IMSTATS_SEED"); + if (!status) { + seed = 0; + } + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, seed); + + // subtract this amount extra from the sky + float SKY_BIAS = psMetadataLookupF32 (&status, recipe, "SKY_BIAS"); + if (!status) { + SKY_BIAS = 0; + } + + // supply the sky background statistics options + char *statsName = psMetadataLookupStr (&status, recipe, "SKY_STAT"); + if (statsName == NULL) { + statsName = defaultStatsName; + } + psStatsOptions statsOptionLocation = psStatsOptionFromString(statsName); + if (!(statsOptionLocation & (PS_STAT_SAMPLE_MEAN | + PS_STAT_SAMPLE_MEDIAN | + PS_STAT_ROBUST_MEDIAN | + PS_STAT_ROBUST_QUARTILE | + PS_STAT_CLIPPED_MEAN | + PS_STAT_FITTED_MEAN | + PS_STAT_FITTED_MEAN_V2 | + PS_STAT_FITTED_MEAN_V3))) { + statsOptionLocation = PS_STAT_FITTED_MEAN; + } + + psStatsOptions statsOptionWidth = PS_STAT_NONE; + if (statsOptionLocation & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN)) { + statsOptionWidth = PS_STAT_SAMPLE_STDEV; + } else if (statsOptionLocation & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_QUARTILE)) { +#if 1 + statsOptionWidth = PS_STAT_ROBUST_STDEV; // not set; => NaN +#else + statsOptionWidth = PS_STAT_FITTED_STDEV; +#endif + } else if (statsOptionLocation & PS_STAT_FITTED_MEAN) { + statsOptionWidth = PS_STAT_FITTED_STDEV; + } else if (statsOptionLocation & PS_STAT_CLIPPED_MEAN) { + statsOptionWidth = PS_STAT_CLIPPED_STDEV; + } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V2) { + statsOptionWidth = PS_STAT_FITTED_STDEV_V2; + } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V3) { + statsOptionWidth = PS_STAT_FITTED_STDEV_V3; + } else { + psAbort("Unable to estimate variance of selected statsOptionLocations 0x%x", statsOptionLocation); + } + + const psStatsOptions statsOption = statsOptionLocation | statsOptionWidth; + psStats *stats = psStatsAlloc (statsOption); + + // set range for old-version of sky statistic + if (statsOptionLocation & PS_STAT_ROBUST_QUARTILE) { + stats->min = 0.25; + stats->max = 0.75; + } + + // set user-option for number of pixels per region + stats->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX"); + if (!status) { + stats->nSubsample = 1000; + } + + // optionally set the binsize + stats->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS"); + if (status) { + stats->options |= PS_STAT_USE_BINSIZE; + } + + // optionally set the sigma clipping + stats->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA"); + if (!status) { + if ((stats->options & PS_STAT_FITTED_MEAN) || (stats->options & PS_STAT_FITTED_MEAN_V2)) { + stats->clipSigma = 1.0; + } else { + stats->clipSigma = 3.0; + } + } + + // stats is not initialized by psStats??? use this to save the input options + // XXX re-work this to use the new psStatsInit function + psStats *statsDefaults = psStatsAlloc (statsOption); + *statsDefaults = *stats; + + // find the currently selected readout + pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + pmFPA *inFPA = file->fpa; + readout = pmFPAviewThisReadout (view, inFPA); + + psImage *image = readout->image; + psImage *mask = readout->mask; + + // I have the fine image size, I know the binning factor, determine the ruff image size + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = image->numCols; + binning->nYfine = image->numRows; + binning->nXbin = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); + binning->nYbin = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); + + psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER); + psImageBinningSetSkip(binning, image); + status = psMetadataAddPtr(recipe, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning); + PS_ASSERT (status, false); + + // we save the binning structure for use in psphotMagnitudes + pmReadout *model = get_model_readout("PSPHOT.BACKMDL", config, view, inFPA, binning); + pmReadout *modelStdev = get_model_readout("PSPHOT.BACKMDL.STDEV", config, view, inFPA, binning); + + psF32 **modelData = model->image->data.F32; + psF32 **modelStdevData = modelStdev->image->data.F32; + + // measure clipped median for subimages + psRegion ruffRegion = {0,0,0,0}; + psRegion fineRegion = {0,0,0,0}; + for (int iy = 0; iy < model->image->numRows; iy++) { + for (int ix = 0; ix < model->image->numCols; ix++) { + + // convert the ruff grid cell to the equivalent fine grid cell + // XXX we need to watch out for row0,col0 + ruffRegion = psRegionSet (ix, ix + 2, iy, iy + 2); + fineRegion = psImageBinningSetFineRegion (binning, ruffRegion); + fineRegion = psRegionForImage (image, fineRegion); + + psImage *subset = psImageSubset (image, fineRegion); + if (!subset->numCols || !subset->numRows) { + psFree (subset); + continue; + } + psImage *submask = psImageSubset (mask, fineRegion); + + // reset the default values + *stats = *statsDefaults; + + // Use the selected background statistic for the first pass + // If it fails, fall back on the "ROBUST_MEDIAN" version + // If both fail, set the pixel to NAN and (below) interpolate + // XXX psImageBackground will probably be renamed psImageStats + // XXX don't bother trying if there are no valid pixels... + + if (psImageBackground(stats, subset, submask, maskVal, rng)) { + if (stats->options & PS_STAT_ROBUST_QUARTILE) { + modelData[iy][ix] = stats->robustMedian; + } else { + modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation); + } + modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth); + } else { + psStatsOptions currentOptions = stats->options; + stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV; + if (!psImageBackground(stats, subset, submask, maskVal, rng)) { + psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for " + "(%dx%d, (row0,col0) = (%d,%d)", + subset->numRows, subset->numCols, subset->row0, subset->col0); + modelData[iy][ix] = modelStdevData[iy][ix] = NAN; + } else { + modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN); + modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV); + } + // drop errors caused by psImageBackground failures + // XXX we probably should trap and exit on serious failures + psErrorClear(); + stats->options = currentOptions; + } + modelData[iy][ix] += SKY_BIAS; + psFree (subset); + psFree (submask); + } + } + + // patch over bad regions (use average of 8 possible neighbor pixels) + // XXX consider testing all pixels against the 8 neighbors and replacing outliers... + double Count = 0; // number of good pixels + double Value = 0; // sum of good pixel's value + double ValueStdev = 0; // sum of good pixel's standard deviations + for (int iy = 0; iy < model->image->numRows; iy++) { + for (int ix = 0; ix < model->image->numCols; ix++) { + if (!isnan(modelData[iy][ix])) { + Value += modelData[iy][ix]; + ValueStdev += modelStdevData[iy][ix]; + Count++; + continue; + } + + double value = 0; + double count = 0; + for (int jy = iy - 1; jy <= iy + 1; jy++) { + if (jy < 0) continue; + if (jy >= model->image->numRows) continue; + for (int jx = ix - 1; jx <= ix + 1; jx++) { + if (!jx && !jy) continue; + if (jx < 0) continue; + if (jx >= model->image->numCols) continue; + value += modelData[jy][jx]; + count += 1.0; + } + } + if (count > 0) modelData[iy][ix] = value / count; + } + } + if (Count == 0) { + psError (PSPHOT_ERR_DATA, true, "failed to build background image"); + psFree(stats); + psFree(statsDefaults); + psFree(binning); + psFree(rng); + return false; + } + + Value /= Count; + ValueStdev /= Count; + + // patch over remaining bad regions (use global average) + for (int iy = 0; iy < model->image->numRows; iy++) { + for (int ix = 0; ix < model->image->numCols; ix++) { + if (!isnan(modelData[iy][ix])) continue; + modelData[iy][ix] = Value; + modelStdevData[iy][ix] = ValueStdev; + } + } + + psLogMsg ("psphot", PS_LOG_MINUTIA, "build median image: %f sec\n", psTimerMark ("psphot")); + + psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value); + psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev); + psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f stdev %f", Value, ValueStdev); + + // measure image and background stats and save for later output + psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN | + PS_STAT_SAMPLE_STDEV | + PS_STAT_MIN | + PS_STAT_MAX); + psImageStats (statsBck, model->image, NULL, 0); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MN", + PS_META_REPLACE, "sky model mean", statsBck->sampleMean); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_SIG", + PS_META_REPLACE, "sky model stdev", statsBck->sampleStdev); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MAX", + PS_META_REPLACE, "sky model maximum value", statsBck->max); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MIN", + PS_META_REPLACE, "sky model minimum value", statsBck->min); + psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NX", + PS_META_REPLACE, "sky model size (x)", model->image->numCols); + psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NY", + PS_META_REPLACE, "sky model size (y)", model->image->numRows); + psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f", + statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev); + psFree (statsBck); + + // select background pixels, from output background file, or create + file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND"); + if (file) { + // we are using PSPHOT.BACKGND as an I/O file: select readout or create + if (file->mode == PM_FPA_MODE_INTERNAL) { + background = file->readout; + } else { + background = pmFPAviewThisReadout (view, file->fpa); + } + if (background == NULL) { + // readout does not yet exist: create from input + pmFPAfileCopyStructureView (file->fpa, inFPA, 1, 1, view); + background = pmFPAviewThisReadout (view, file->fpa); + if ((image->numCols != background->image->numCols) || (image->numRows != background->image->numRows)) { + psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for background dimensions"); + return false; + } + } + } else { + background = pmFPAfileDefineInternal (config->files, "PSPHOT.BACKGND", image->numCols, image->numRows, PS_TYPE_F32); + } + psF32 **backData = background->image->data.F32; + + // linear interpolation to full-scale + if (!psImageUnbin (background->image, model->image, binning)) { + psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning"); + return false; + } + + psLogMsg ("psphot", PS_LOG_MINUTIA, "build resampled image: %f sec\n", psTimerMark ("psphot")); + + // back-sub image pixels, from output background file (don't create if not requested) + file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKSUB"); + if (file) { + // we are using PSPHOT.BACKSUB as an I/O file: select readout or create + backSub = pmFPAviewThisReadout (view, file->fpa); + if (backSub == NULL) { + // readout does not yet exist: create from input + pmFPAfileCopyStructureView (file->fpa, inFPA, 1, 1, view); + backSub = pmFPAviewThisReadout (view, file->fpa); + } + } + + if (psTraceGetLevel("psphot") > 5) { + char name[256]; + sprintf (name, "image.%02d.fits", npass); + psphotSaveImage (NULL, image, name); + sprintf (name, "back.%02d.fits", npass); + psphotSaveImage (NULL, background->image, name); + sprintf (name, "mask.%02d.fits", npass); + psphotSaveImage (NULL, mask, name); + sprintf (name, "backmdl.%02d.fits", npass); + psphotSaveImage (NULL, model->image, name); + } + + // subtract the background model (save in backSub, if requested) + for (int j = 0; j < image->numRows; j++) { + for (int i = 0; i < image->numCols; i++) { + image->data.F32[j][i] -= backData[j][i]; + if (backSub) { + backSub->image->data.F32[j][i] = image->data.F32[j][i]; + } + } + } + + if (psTraceGetLevel("psphot") > 5) { + char name[256]; + sprintf (name, "backsub.%02d.fits", npass); + psphotSaveImage (NULL, image, name); + } + npass ++; + + psLogMsg ("psphot", PS_LOG_INFO, "subtracted background model: %f sec\n", psTimerMark ("psphot")); + psFree(stats); + psFree(statsDefaults); + psFree(binning); + psFree(rng); + + // the pmReadout selected in this function are all view on entries in config->files + return true; +}