Index: unk/archive/scripts/src/phase2/pmChipMosaic.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmChipMosaic.c (revision 5798) +++ (revision ) @@ -1,175 +1,0 @@ -#include -#include - -#include "pslib.h" -#include "pmFPA.h" -#include "pmChipMosaic.h" - -// Compare a value with a maximum and minimum -#define COMPARE(value,min,max) \ - if ((value) < (min)) { \ - (min) = (value); \ - } \ - if ((value) > (max)) { \ - (max) = (value); \ - } - -// Mosaic multiple images, with flips, binning and offsets -psImage *p_pmImageMosaic(const psArray *source, // Images to splice in - const psVector *xFlip, const psVector *yFlip, // Need to flip x and y? - const psVector *xBinSource, const psVector *yBinSource, // Binning in x and y of - // source images - int xBinTarget, int yBinTarget, // Binning in x and y of target images - const psVector *x0, const psVector *y0 // Offsets for source images on target - ) -{ - assert(source); - assert(xFlip && xFlip->type.type == PS_TYPE_U8); - assert(yFlip && yFlip->type.type == PS_TYPE_U8); - assert(xBinSource && xBinSource->type.type == PS_TYPE_S32); - assert(yBinSource && yBinSource->type.type == PS_TYPE_S32); - assert(x0 && x0->type.type == PS_TYPE_S32); - assert(y0 && y0->type.type == PS_TYPE_S32); - - // Get the maximum extent of the mosaic image - int xMin = INT_MAX; - int xMax = 0; - int yMin = INT_MAX; - int yMax = 0; - for (int i = 0; i < source->n; i++) { - psImage *image = source->data[i]; // The image of interest - - assert(image->type.type == PS_TYPE_F32); // Only implemented for F32 images so far. - - // Size of cell in x and y - int xParity = xFlip->data.U8[i] ? -1 : 1; - int yParity = yFlip->data.U8[i] ? -1 : 1; - psTrace(__func__, 5, "Extent of cell %d: %d -> %d , %d -> %d\n", i, x0->data.S32[i], - x0->data.S32[i] + xParity * xBinSource->data.S32[i] * image->numCols, y0->data.S32[i], - y0->data.S32[i] + yParity * yBinSource->data.S32[i] * image->numRows); - - COMPARE(x0->data.S32[i], xMin, xMax); - COMPARE(y0->data.S32[i], yMin, yMax); - // Subtract the parity to get the inclusive limit (not exclusive) - COMPARE(x0->data.S32[i] + xParity * xBinSource->data.S32[i] * image->numCols - xParity, xMin, xMax); - COMPARE(y0->data.S32[i] + yParity * yBinSource->data.S32[i] * image->numRows - yParity, yMin, yMax); - } - - // Set up the image - // Since both upper and lower values are inclusive, we need to add one to the size - float xSize = (float)(xMax - xMin + 1) / (float)xBinTarget; - if (xSize - (int)xSize > 0) { - xSize += 1; - } - float ySize = (float)(yMax - yMin + 1) / (float)yBinTarget; - if (ySize - (int)ySize > 0) { - ySize += 1; - } - - psTrace(__func__, 3, "Spliced image will be %dx%d\n", (int)xSize, (int)ySize); - psImage *mosaic = psImageAlloc((int)xSize, (int)ySize, PS_TYPE_F32); // The mosaic image - psImageInit(mosaic, 0.0); - - // Next pass through the images to do the mosaicking - for (int i = 0; i < source->n; i++) { - psImage *image = source->data[i]; // The image of interest - if (xBinSource->data.S32[i] == xBinTarget && yBinSource->data.S32[i] == yBinTarget && - xFlip->data.U8[i] == 0 && yFlip->data.U8[i] == 0) { - // Let someone else do the hard work; useful to test psImageOverlaySection if no other reason - psImageOverlaySection(mosaic, image, x0->data.S32[i], y0->data.S32[i], "+"); - } else { - // We have to do the hard work ourself - for (int y = 0; y < image->numRows; y++) { - int yParity = yFlip->data.U8[i] ? -1 : 1; - float yTargetBase = (y0->data.S32[i] + yParity * yBinSource->data.S32[i] * y) / yBinTarget; - for (int x = 0; x < image->numCols; x++) { - int xParity = xFlip->data.U8[i] ? -1 : 1; - float xTargetBase = (x0->data.S32[i] + xParity * xBinSource->data.S32[i] * x) / - xBinTarget; - - // In case the original image is binned but the mosaic is not, we need to fill in the - // values in the mosaic. - for (int j = 0; j < yBinSource->data.S32[i]; j++) { - int yTarget = (int)(yTargetBase + yParity * (float)j / (float)yBinTarget); - for (int i = 0; i < xBinSource->data.S32[i]; i++) { - int xTarget = (int)(xTargetBase + xParity * (float)i / (float)xBinTarget); - - mosaic->data.F32[yTarget][xTarget] += image->data.F32[y][x]; - } - } // Iterating over mosaic image for binned input image - } - } // Iterating over input image - } - } - - return mosaic; -} - -psImage *pmChipMosaic(pmChip *chip, // Chip to mosaic - int xBinChip, int yBinChip // Binning of mosaic image in x and y - ) -{ - psArray *cells = chip->cells; // The array of cells - psArray *images = psArrayAlloc(cells->n); // Array of images that will be mosaicked - psVector *x0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin x coordinates - psVector *y0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin y coordinates - psVector *xBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in x - psVector *yBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in y - psVector *xFlip = psVectorAlloc(cells->n, PS_TYPE_U8); // Flip in x? - psVector *yFlip = psVectorAlloc(cells->n, PS_TYPE_U8); // Flip in y? - - // Set up the required inputs - for (int i = 0; i < cells->n; i++) { - pmCell *cell = cells->data[i]; // The cell of interest - x0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); - y0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); - xBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); - yBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); - int xParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); - int yParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); - if (xParity == 1) { - xFlip->data.U8[i] = 0; - } else if (xParity == -1) { - xFlip->data.U8[i] = 1; - } else { - psLogMsg(__func__, PS_LOG_WARN, "The x parity of cell %d is not +/- 1 (it's %d) --- " - "assuming +1.\n", i, xParity); - xFlip->data.U8[i] = 0; - } - if (yParity == 1) { - yFlip->data.U8[i] = 0; - } else if (yParity == -1) { - yFlip->data.U8[i] = 1; - } else { - psLogMsg(__func__, PS_LOG_WARN, "The y parity of cell %d is not +/- 1 (it's %d) --- " - "assuming +1.\n", i, yParity); - yFlip->data.U8[i] = 0; - } - - // Trim the image to get rid of the overscan - psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); - psTrace(__func__, 7, "Cell %d trimsec: [%.0f:%.0f,%.0f:%.0f]\n", i, trimsec->x0, trimsec->x1, - trimsec->y0, trimsec->y1); - psArray *readouts = cell->readouts; // The array of readouts - if (readouts->n > 1) { - psLogMsg(__func__, PS_LOG_WARN, "Cell %d contains more than one readout --- only the first will " - "be mosaicked.\n", i); - } - psImage *image = ((pmReadout*)readouts->data[0])->image; // The image to put into the mosaic - images->data[i] = psImageSubset(image, *trimsec); // Trimmed image - } - - // Mosaic the images together and we're done - psImage *mosaic = p_pmImageMosaic(images, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0); - - // Clean up - psFree(x0); - psFree(y0); - psFree(xBin); - psFree(yBin); - psFree(xFlip); - psFree(yFlip); - psFree(images); - - return mosaic; -} Index: unk/archive/scripts/src/phase2/pmChipMosaic.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmChipMosaic.h (revision 5798) +++ (revision ) @@ -1,22 +1,0 @@ -#ifndef PM_CHIP_MOSAIC_H -#define PM_CHIP_MOSAIC_H - -#include "pslib.h" -#include "pmFPA.h" - -// Mosaic multiple images, with flips, binning and offsets -psImage *p_pmImageMosaic(const psArray *source, // Images to splice in - const psVector *xFlip, const psVector *yFlip, // Need to flip x and y? - const psVector *xBinSource, const psVector *yBinSource, // Binning in x and y of - // source images - int xBinTarget, int yBinTarget, // Binning in x and y of target images - const psVector *x0, const psVector *y0 // Offsets for source images on target - ); - -// Mosaic all the cells in a chip together (neglecting the overscans) -psImage *pmChipMosaic(pmChip *chip, // Chip to mosaic - int xBinChip, int yBinChip // Binning of mosaic image in x and y - ); - - -#endif Index: unk/archive/scripts/src/phase2/pmConfig.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmConfig.c (revision 5798) +++ (revision ) @@ -1,308 +1,0 @@ -#include -#include -#include -#include "pslib.h" -#include "psAdditionals.h" -#include "pmConfig.h" - -#define PS_SITE "PS_SITE" // Name of the environment variable containing the site config file -#define DEFAULT_SITE "ipprc.config" // Default site config file - -static bool readConfig(psMetadata **config, // Config to output - const char *name,// Name of file - const char *description // Description of file - ) -{ - int numBadLines = 0; // Number of bad lines in config file - psLogMsg(__func__, PS_LOG_INFO, "Loading %s configuration from file %s\n", description, name); - *config = psMetadataConfigParse(NULL, &numBadLines, name, true); - if (numBadLines > 0) { - psLogMsg(__func__, PS_LOG_WARN, "%d bad lines in %s configuration file (%s)\n", description, - name); - } - if (! *config) { - psError(PS_ERR_IO, false, "Unable to read %s configuration from %s\n", description, name); - return false; - } - - return true; -} - - -bool pmConfigRead(psMetadata **site, psMetadata **camera, psMetadata **recipe, - int *argc, char **argv, const char *recipeName) -{ - // Make sure we've been given the correct inputs, and that we won't leak memory - assert(site && *site == NULL); - assert(camera && *camera == NULL); - assert(recipe && *recipe == NULL); - assert(*argc > 0); - assert(argv); - - char *siteName = NULL; // Name of the site configuration file - - // First, try command line - int argNum = 0; // Number of the site argument - if (argNum = psArgumentGet(*argc, argv, "-site")) { - (void)psArgumentRemove(argNum, argc, argv); - if (argNum >= *argc) { - psLogMsg(__func__, PS_LOG_WARN, - "-site command-line switch provided without the required filename --- ignored.\n"); - } else { - siteName = argv[argNum]; - (void)psArgumentRemove(argNum, argc, argv); - } - } - // Next, try environment variable - if (! siteName) { - siteName = getenv(PS_SITE); - } - // Last chance is ~/.ipprc - bool cleanupSiteName = false; // Do I have to psFree siteName? - if (! siteName) { - siteName = psStringCopy(DEFAULT_SITE); - cleanupSiteName = true; - } - - if (! readConfig(site, siteName, "site")) { - if (cleanupSiteName) { - psFree(siteName); - } - return false; - } - - // Next is the camera configuration - if (argNum = psArgumentGet(*argc, argv, "-camera")) { - (void)psArgumentRemove(argNum, argc, argv); - if (argNum >= *argc) { - psLogMsg(__func__, PS_LOG_WARN, - "-camera command-line switch provided without the required filename --- ignored.\n"); - } else { - (void)psArgumentRemove(argNum, argc, argv); - (void)readConfig(camera, argv[argNum], "camera"); - } - } - - // And then the recipe configuration - if (argNum = psArgumentGet(*argc, argv, "-recipe")) { - (void)psArgumentRemove(argNum, argc, argv); - if (argNum >= *argc) { - psLogMsg(__func__, PS_LOG_WARN, - "-recipe command-line switch provided without the required filename --- ignored.\n"); - } else { - (void)psArgumentRemove(argNum, argc, argv); - (void)readConfig(recipe, argv[argNum], "recipe"); - } - } - // Or, load the recipe from the camera file, if appropriate - if (! *recipe && *camera && recipeName) { - *recipe = pmConfigRecipeFromCamera(*camera, recipeName); - } - - - // Now we can look into the site configuration and do the required stuff - bool mdok = true; // Status of MD lookup result - psString timeName = psMetadataLookupString(&mdok, *site, "TIME"); // Name of time file - if (mdok && timeName) { - psTrace(__func__, 7, "Initialising psTime with file %s\n", timeName); -#ifdef PRODUCTION - psTimeInitialize(timeName); -#else - psLibInit(timeName); -#endif - } - - int logLevel = psMetadataLookupS32(&mdok, *site, "LOGLEVEL"); // Logging level - if (mdok && logLevel >= 0) { - psTrace(__func__, 7, "Setting log level to %d\n", logLevel); - psLogSetLevel(logLevel); - } - - psString logFormat = psMetadataLookupString(&mdok, *site, "LOGLEVEL"); // Log format - if (mdok && logFormat) { - psTrace(__func__, 7, "Setting log format to %s\n", logFormat); - psLogSetFormat(logFormat); - } - - psString logDest = psMetadataLookupString(&mdok, *site, "LOGDEST"); // Log destination - if (mdok && logDest) { - // XXX: Only stdout is provided for now; this section should be expanded in the future to do files, - // and perhaps even sockets. - if (strcasecmp(logDest, "STDOUT") != 0) { - psLogMsg(__func__, PS_LOG_WARN, "Only STDOUT is currently supported as a log destination.\n"); - } - psTrace(__func__, 7, "Setting log destination to STDOUT.\n"); - psLogSetDestination(PS_LOG_TO_STDOUT); - } - - psMetadata *trace = psMetadataLookupMD(&mdok, *site, "TRACE"); // Trace levels - if (mdok && trace) { - psMetadataIterator *traceIter = psMetadataIteratorAlloc(trace, PS_LIST_HEAD, NULL); // Iterator - psMetadataItem *traceItem = NULL; // Item from MD iteration - while (traceItem = psMetadataGetAndIncrement(traceIter)) { - if (traceItem->type != PS_DATA_S32) { - psLogMsg(__func__, PS_LOG_WARN, "The level for trace component %s is not of type S32 (%x)\n", - traceItem->name, traceItem->type); - continue; - } - psTrace(__func__, 7, "Setting trace level for %s to %d\n", traceItem->name, traceItem->data.S32); - psTraceSetLevel(traceItem->name, traceItem->data.S32); - } - psFree(traceIter); - } - - if (cleanupSiteName) { - psFree(siteName); - } - return true; -} - -bool pmConfigValidateCamera(const psMetadata *camera, const psMetadata *header) -{ - // Read the rule for that camera - bool mdStatus = true; // Status of MD lookup - psMetadata *rule = psMetadataLookupMD(&mdStatus, camera, "RULE"); - if (! mdStatus || ! rule) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to read rule for camera.\n"); - return false; - } - - // Apply the rules - psMetadataIterator *ruleIter = psMetadataIteratorAlloc(rule, PS_LIST_HEAD, NULL); // Rule iterator - psMetadataItem *ruleItem = NULL; // Item from the metadata - bool match = true; // Does it match? - while ((ruleItem = psMetadataGetAndIncrement(ruleIter)) && match) { - // Check for the existence of the rule - psMetadataItem *headerItem = psMetadataLookup((psMetadata*)header, ruleItem->name); - if (! headerItem || headerItem->type != ruleItem->type) { - match = false; - break; - } - - // Check to see if the rule works - switch (ruleItem->type) { - case PS_DATA_STRING: - psTrace(__func__, 8, "Matching %s: '%s' vs '%s'\n", ruleItem->name, - ruleItem->data.V, headerItem->data.V); - if (strncmp(ruleItem->data.V, headerItem->data.V, - strlen(ruleItem->data.V)) != 0) { - match = false; - } - break; - case PS_TYPE_S32: - case PS_TYPE_BOOL: - psTrace(__func__, 8, "Matching %s: %d vs %d\n", ruleItem->name, - ruleItem->data.S32, headerItem->data.S32); - if (ruleItem->data.S32 != headerItem->data.S32) { - match = false; - } - break; - case PS_TYPE_F32: - psTrace(__func__, 8, "Matching %s: %f vs %f\n", ruleItem->name, - ruleItem->data.F32, headerItem->data.F32); - if (ruleItem->data.F32 != headerItem->data.F32) { - match = false; - } - break; - case PS_TYPE_F64: - psTrace(__func__, 8, "Matching %s: %g vs %g\n", ruleItem->name, - ruleItem->data.F64, headerItem->data.F64); - if (ruleItem->data.F64 != headerItem->data.F64) { - match = false; - } - break; - default: - psLogMsg(__func__, PS_LOG_WARN, "Ignoring invalid type in metadata: %x\n", - ruleItem->type); - } - } // Iterating through the RULEs - - psFree(ruleIter); - - return match; -} - - - -// Work out what camera we have, based on the FITS header and a set of rules specified in the IPP -// configuration; return the camera configuration -psMetadata *pmConfigCameraFromHeader(const psMetadata *ipprc, // The IPP configuration - const psMetadata *header // The FITS header - ) -{ - bool mdStatus = false; // Metadata lookup status - psMetadata *cameras = psMetadataLookupMD(&mdStatus, ipprc, "CAMERAS"); - if (! mdStatus) { - psError(PS_ERR_IO, false, "Unable to find CAMERAS in the configuration.\n"); - return NULL; - } - - psMetadata *winner = NULL; // The camera configuration whose rule first matches the supplied header - - // Iterate over the cameras - psMetadataIterator *iterator = psMetadataIteratorAlloc(cameras, PS_LIST_HEAD, NULL); // MD Iterator - psMetadataItem *cameraItem = NULL; // Item from the metadata - while (cameraItem = psMetadataGetAndIncrement(iterator)) { - // Open the camera information - psTrace(__func__, 3, "Inspecting camera %s (%s)\n", cameraItem->name, - cameraItem->comment); - psMetadata *camera = NULL; // The camera metadata - if (cameraItem->type == PS_DATA_METADATA) { - camera = psMemIncrRefCounter(cameraItem->data.md); - } else if (cameraItem->type == PS_DATA_STRING) { - psTrace(__func__, 5, "Reading camera configuration for %s...\n", cameraItem->name); - int badLines = 0; // Number of bad lines in reading camera configuration - camera = psMetadataConfigParse(NULL, &badLines, cameraItem->data.V, true); - if (badLines > 0) { - psLogMsg(__func__, PS_LOG_WARN, "%d bad lines encountered while reading camera" - "configuration %s\n", badLines, cameraItem->name); - } - } - - if (! camera) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to interpret camera configuration for %s (%s)\n", - cameraItem->name, cameraItem->comment); - continue; - } - - if (pmConfigValidateCamera(camera, header)) { - if (! winner) { - // This is the first match - winner = psMemIncrRefCounter(camera); - psLogMsg(__func__, PS_LOG_INFO, "FITS header matches camera %s\n", - cameraItem->name); - } else { - // We have a duplicate match - psLogMsg(__func__, PS_LOG_WARN, "Additional camera found that matches the rules: %s\n", - cameraItem->name); - } - } // Done inspecting the camera - - psFree(camera); - - } // Done looking at all cameras - if (! winner) { - psError(PS_ERR_IO, true, "Unable to find a camera that matches input FITS header!\n"); - } - - psFree(iterator); - return winner; -} - -psMetadata *pmConfigRecipeFromCamera(const psMetadata *camera, const char *recipeName) -{ - assert(camera); - assert(recipeName); - - psMetadata *recipe = NULL; // Recipe to read - bool mdok = true; // Status of MD lookup - psMetadata *recipes = psMetadataLookupMD(&mdok, camera, "RECIPES"); // The list of recipes - if (! mdok || ! recipes) { - psLogMsg(__func__, PS_LOG_WARN, "RECIPES in the camera configuration file is not of type METADATA\n"); - } else { - psString recipeFileName = psMetadataLookupString(&mdok, recipes, recipeName); - (void)readConfig(&recipe, recipeFileName, "recipe"); - } - - return recipe; -} Index: unk/archive/scripts/src/phase2/pmConfig.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmConfig.h (revision 5798) +++ (revision ) @@ -1,12 +1,0 @@ -#ifndef PM_CONFIG_H -#define PM_CONFIG_H - -#include "pslib.h" - -bool pmConfigRead(psMetadata **site, psMetadata **camera, psMetadata **recipe, - int *argc, char **argv, const char *recipeName); -bool pmConfigValidateCamera(const psMetadata *camera, const psMetadata *header); -psMetadata *pmConfigCameraFromHeader(const psMetadata *site, const psMetadata *header); -psMetadata *pmConfigRecipeFromCamera(const psMetadata *camera, const char *recipeName); - -#endif Index: unk/archive/scripts/src/phase2/pmFPA.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPA.c (revision 5798) +++ (revision ) @@ -1,293 +1,0 @@ -#include -#include -#include -#include "pslib.h" -#include "psAdditionals.h" - -#include "pmFPA.h" - -// Get the bias images for a readout, using the CELL.BIASSEC -psList *pmReadoutGetBias(pmReadout *readout // Readout for which to get the bias sections - ) -{ - pmCell *cell = readout->parent; // The parent cell - psList *sections = (psList*)psMetadataLookupPtr(NULL, cell->concepts, "CELL.BIASSEC"); // CELL.BIASSEC - - psImage *image = readout->image; // The image that contains both the biassec and the trimsec - - psList *images = psListAlloc(NULL); // List of images from bias sections - psListIterator *sectionsIter = psListIteratorAlloc(sections, PS_LIST_HEAD, true); // Iterator - psRegion *region = NULL; // Bias region from list - while (region = psListGetAndIncrement(sectionsIter)) { -#if 0 - // Convert from FITS standard to PS standard - // We don't touch the x1 and y1 values because they aren't included by psImageSubset. - region->x0 -= 1; - region->y0 -= 1; -#endif - - psImage *bias = psImageSubset(image, *region); // Image from bias section - psListAdd(images, PS_LIST_TAIL, bias); - psFree(bias); - } - psFree(sectionsIter); - - return images; -} - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Allocators -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -p_pmHDU *p_pmHDUAlloc(const char *extname) -{ - p_pmHDU *pd = psAlloc(sizeof(p_pmHDU)); - psMemSetDeallocator(pd, (psFreeFunc)p_pmHDUFree); - - pd->extname = extname; - pd->header = NULL; - pd->images = NULL; - pd->masks = NULL; - pd->weights = NULL; - - return pd; -} - -void p_pmHDUFree(p_pmHDU *hdu) -{ - psFree(hdu->header); - psFree(hdu->images); - psFree(hdu->masks); - psFree(hdu->weights); -} - -pmFPA *pmFPAAlloc(const psMetadata *camera // Camera configuration - ) -{ - pmFPA *fpa = psAlloc(sizeof(pmFPA));// The FPA - psMemSetDeallocator(fpa, (psFreeFunc)p_pmFPAFree); - - // Fill in the components - fpa->fromTangentPlane = NULL; - fpa->toTangentPlane = NULL; - fpa->projection = NULL; - - fpa->concepts = psMetadataAlloc(); - fpa->camera = psMemIncrRefCounter((psPtr)camera); - fpa->chips = psArrayAlloc(0); - - fpa->phu = NULL; - fpa->hdu = NULL; - - return fpa; -} - -void p_pmFPAFree(pmFPA *fpa) -{ - psFree(fpa->fromTangentPlane); - psFree(fpa->toTangentPlane); - psFree(fpa->projection); - - psFree(fpa->concepts); - psFree(fpa->camera); - psFree(fpa->chips); - - psFree(fpa->phu); - psFree(fpa->hdu); -} - -pmChip *pmChipAlloc(pmFPA *fpa, // FPA to which the chip belongs - psString name // Chip name - ) -{ - pmChip *chip = psAlloc(sizeof(pmChip)); // The chip - psMemSetDeallocator(chip, (psFreeFunc)p_pmChipFree); - - // Push onto the array of chips - fpa->chips = psArrayAdd(fpa->chips, 0, chip); - - // Fill in the components - *(int*)&chip->col0 = 0; // Good enough for now - *(int*)&chip->row0 = 0; // Good enough for now - - chip->toFPA = NULL; - chip->fromFPA = NULL; - - chip->concepts = psMetadataAlloc(); - chip->cells = psArrayAlloc(0); - chip->parent = fpa; // We don't increment the reference counter on this --- it's a - // "hidden" link. If we increment the reference counter, we get stuck - // in a circle. - chip->valid = true; - - chip->hdu = NULL; - - psMetadataAddStr(chip->concepts, PS_LIST_HEAD, "CHIP.NAME", 0, "Chip name added at pmChipAlloc", name); - - return chip; -} - -void p_pmChipFree(pmChip *chip) -{ - psFree(chip->toFPA); - psFree(chip->fromFPA); - - psFree(chip->concepts); - psFree(chip->cells); - - psFree(chip->hdu); - - // We don't free the parent member, since that would generate a circular call. We don't increment the - // reference counter when we add it, anyway, so that's OK. -} - -pmCell *pmCellAlloc(pmChip *chip, // Chip to which the cell belongs - psMetadata *cameraData, // Camera data - psString name // Name of cell - ) -{ - pmCell *cell = psAlloc(sizeof(pmCell)); // The cell - psMemSetDeallocator(cell, (psFreeFunc)p_pmCellFree); - - // Push onto the array of chips - chip->cells = psArrayAdd(chip->cells, 0, cell); - - // Fill in components - *(int*)&cell->col0 = 0; // Good enough for now - *(int*)&cell->row0 = 0; // Good enough for now - - cell->toChip = NULL; - cell->toFPA = NULL; - cell->toSky = NULL; - - cell->concepts = psMetadataAlloc(); - psMetadataAddStr(cell->concepts, PS_LIST_HEAD, "CELL.NAME", 0, "Cell name added at pmCellAlloc", name); - cell->camera = psMemIncrRefCounter(cameraData); - - cell->readouts = psArrayAlloc(0); - cell->parent = chip; // We don't increment the reference counter on this --- it's a - // "hidden" link. If we increment the reference counter, we get stuck - // in a circle. - cell->valid = true; - - cell->hdu = NULL; - - return cell; -} - -void p_pmCellFree(pmCell *cell) -{ - psFree(cell->toChip); - psFree(cell->toFPA); - psFree(cell->toSky); - - psFree(cell->concepts); - psFree(cell->camera); - psFree(cell->readouts); - - psFree(cell->hdu); - - // We don't free the parent member, since that would generate a circular call. We don't increment the - // reference counter when we add it, anyway, so that's OK. -} - -pmReadout *pmReadoutAlloc(pmCell *cell, // Cell to which the readout belongs - psImage *image, // The pixels - psImage *mask,// The mask pixels - int col0, int row0, int colBin, int rowBin // Data - ) -{ - pmReadout *readout = psAlloc(sizeof(pmReadout)); - psMemSetDeallocator(readout, (psFreeFunc)p_pmReadoutFree); - if (cell) { - cell->readouts = psArrayAdd(cell->readouts, 0, readout); - } - - // Set the components - readout->image = psMemIncrRefCounter(image); - readout->mask = psMemIncrRefCounter(mask); - readout->weight = NULL; - - //readout->concepts = psMetadataAlloc(); - - *(int*)&readout->col0 = col0; - *(int*)&readout->row0 = row0; - *(int*)&readout->colBins = colBin; - *(int*)&readout->rowBins = rowBin; - - readout->parent = cell; - - return readout; -} - -void p_pmReadoutFree(pmReadout *readout) -{ - psFree(readout->image); - psFree(readout->mask); - psFree(readout->weight); - //psFree(readout->concepts); -} - - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Select and Exclude chips -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -bool pmFPASelectChip(pmFPA *fpa, int chipNum) -{ - assert(fpa); - - if (chipNum < 0 || chipNum > fpa->chips->n) { - return false; - } - psArray *chips = fpa->chips; // Component chips - for (int i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // The chip of interest - if (i == chipNum) { - psTrace(__func__, 5, "Marking chip %d valid.\n", i); - chip->valid = true; - psArray *cells = chip->cells; // Component cells - for (int j = 0; j < cells->n; j++) { - pmCell *cell = cells->data[j]; // Cell of interest - cell->valid = true; - } - } else { - psTrace(__func__, 5, "Marking chip %d invalid.\n", i); - chip->valid = false; - psArray *cells = chip->cells; // Component cells - for (int j = 0; j < cells->n; j++) { - pmCell *cell = cells->data[j]; // Cell of interest - cell->valid = false; - } - } - } - return true; -} - -int pmFPAExcludeChip(pmFPA *fpa, int chipNum) -{ - assert(fpa); - - if (chipNum < 0 || chipNum > fpa->chips->n) { - psLogMsg(__func__, PS_LOG_WARN, "Invalid chip number: %d\n", chipNum); - } - int numValid = 0; // Number of valid chips - psArray *chips = fpa->chips; // Component chips - for (int i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // The chip of interest - if (i == chipNum) { - psTrace(__func__, 5, "Marking chip %d invalid.\n", i); - chip->valid = false; - psArray *cells = chip->cells; // Component cells - for (int j = 0; j < cells->n; j++) { - pmCell *cell = cells->data[j]; // Cell of interest - cell->valid = false; - } - } else if (chip->valid) { - numValid++; - } - } - - psTrace(__func__, 5, "%d valid chips in fpa.\n", numValid); - return numValid; -} Index: unk/archive/scripts/src/phase2/pmFPA.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPA.h (revision 5798) +++ (revision ) @@ -1,135 +1,0 @@ -#ifndef PM_FPA_H -#define PM_FPA_H - -#include "pslib.h" -#include "psAdditionals.h" - -// Temporary metadata types -//#define PS_META_CHIP PS_DATA_UNKNOWN -//#define PS_META_CELL PS_DATA_UNKNOWN - -typedef struct { - const char *extname; // Extension name, if it corresponds to this level - psMetadata *header; // The FITS header, if it corresponds to this level - psArray *images; // The pixel data, if it corresponds to this level - psArray *masks; // The mask data, if it corresponds to this level - psArray *weights; // The weight data, if it corresponds to this level -} p_pmHDU; - -typedef struct { - // Astrometric transformations - psPlaneDistort* fromTangentPlane; // Transformation from tangent plane to focal plane - psPlaneDistort* toTangentPlane; // Transformation from focal plane to tangent plane - psProjection *projection; // Projection from tangent plane to sky - // Information - psMetadata *concepts; // Cache for PS concepts - psMetadata *analysis; // FPA-level analysis metadata - const psMetadata *camera; // Camera configuration - psArray *chips; // The chips - p_pmHDU *hdu; // FITS data - psMetadata *phu; // Primary Header -} pmFPA; - -typedef struct { - // Offset specifying position on focal plane - int col0; // Offset from the left of FPA. - int row0; // Offset from the bottom of FPA. - // Astrometric transformations - psPlaneTransform* toFPA; // Transformation from chip to FPA coordinates - psPlaneTransform* fromFPA; // Transformation from FPA to chip coordinates - // Information - psMetadata *concepts; // Cache for PS concepts - psMetadata *analysis; // Chip-level analysis metadata - psArray *cells; // The cells (referred to by name) - pmFPA *parent; // Parent FPA - bool valid; // Do we bother about reading and working with this chip? - p_pmHDU *hdu; // FITS data -} pmChip; - -typedef struct { - // Offset specifying position on chip - int col0; // Offset from the left of chip. - int row0; // Offset from the bottom of chip. - // Astrometric transformations - psPlaneTransform* toChip; // Transformations from cell to chip coordinates - psPlaneTransform* toFPA; // Transformations from cell to FPA coordinates - psPlaneTransform* toSky; // Transformations from cell to sky coordinates - // Information - psMetadata *concepts; // Cache for PS concepts - psMetadata *camera; // Camera information - psMetadata *analysis; // Cell-level analysis metadata - psArray *readouts; // The readouts (referred to by number) - pmChip *parent; // Parent chip - bool valid; // Do we bother about reading and working with this cell? - p_pmHDU *hdu; // FITS data -} pmCell; - -typedef struct { - // Position on the cell - int col0; // Offset from the left of cell. - int row0; // Offset from the bottom of cell. - int colBins; // Amount of binning in x-dimension and parity (from sign) - int rowBins; // Amount of binning in y-dimension and parity (from sign) - // Information - psImage *image; // Imaging area of readout - psImage *mask; // Mask for image - psImage *weight; // Weight for image -#if 0 - psList *bias; // List of bias section (sub-)images -#endif - psMetadata *analysis; // Readout-level analysis metadata - pmCell *parent; // Parent cell -} pmReadout; - -#if 0 -typedef struct { - // Details for position on the cell - const int col0; // Offset from the left of cell. - const int row0; // Offset from the bottom of cell. - const int colParity; // Readout Direction X - const int rowParity; // Readout Direction Y - const unsigned int colBins; // Amount of binning in x-dimension - const unsigned int rowBins; // Amount of binning in y-dimension - // Information - psMetadata *concepts; // Concepts for readouts - pmCell *parent; // Parent cell - psImage *image; // The pixels - psImage *mask; // Mask image - psImage *weight; // Weight image -} pmReadout; -#endif - -psList *pmReadoutGetBias(pmReadout *readout // Readout for which to get the bias sections - ); - - -// Allocators and deallocators -p_pmHDU *p_pmHDUAlloc(const char *extname); -void p_pmHDUFree(p_pmHDU *hdu); -pmFPA *pmFPAAlloc(const psMetadata *camera // Camera configuration - ); -void p_pmFPAFree(pmFPA *fpa); - -pmChip *pmChipAlloc(pmFPA *fpa, // FPA to which the chip belongs - psString name // Name of chip - ); -void p_pmChipFree(pmChip *chip); - -pmCell *pmCellAlloc(pmChip *chip, // Chip to which the cell belongs - psMetadata *cameraData, // Camera data - psString name // Name of cell - ); -void p_pmCellFree(pmCell *cell); - -pmReadout *pmReadoutAlloc(pmCell *cell, // Cell to which the readout belongs - psImage *image, // The pixels - psImage *mask,// The mask pixels - int col0, int row0, int colBin, int rowBin // Data - ); -void p_pmReadoutFree(pmReadout *readout); - -// Select and exclude chips -bool pmFPASelectChip(pmFPA *fpa, int chipNum); -int pmFPAExcludeChip(pmFPA *fpa, int chipNum); - -#endif Index: unk/archive/scripts/src/phase2/pmFPAConceptsGet.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAConceptsGet.c (revision 5798) +++ (revision ) @@ -1,1054 +1,0 @@ -#include -#include -#include "pslib.h" - -#include "pmFPA.h" - -#include "papStuff.h" - -#include "pmFPAConceptsGet.h" - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Private functions -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -psMetadataItem *p_pmFPAConceptGetFromCamera(pmCell *cell, // The cell - const char *concept // Name of concept - ) -{ - if (cell) { - psMetadata *camera = cell->camera; // Camera data - // Need the CELL.NAME first, which should be in the cell->concepts since creation of the cell - bool mdStatus = true; // Status of MD lookup - // Finally, the info that we want - psMetadataItem *item = psMetadataLookup(camera, concept); - return item; - } - return NULL; -} - - -psMetadataItem *p_pmFPAConceptGetFromHeader(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - const char *concept // Name of concept - ) -{ - bool mdStatus = true; // Status of MD lookup - psMetadata *translation = psMetadataLookupMD(&mdStatus, fpa->camera, "TRANSLATION"); // FITS translation - if (! mdStatus) { - psError(PS_ERR_IO, false, "Unable to find TRANSLATION in camera configuration.\n"); - return NULL; - } - - // Look for how to translate the concept into a FITS header name - const char *keyword = psMetadataLookupString(&mdStatus, translation, concept); - if (mdStatus && strlen(keyword) > 0) { - // We have a FITS header to look up --- search each level - if (cell && cell->hdu) { - psMetadataItem *cellItem = psMetadataLookup(cell->hdu->header, keyword); - if (cellItem) { - // XXX: Need to clean up before returning - return cellItem; - } - } - - if (chip && chip->hdu) { - psMetadataItem *chipItem = psMetadataLookup(chip->hdu->header, keyword); - if (chipItem) { - // XXX: Need to clean up before returning - return chipItem; - } - } - - if (fpa->hdu) { - psMetadataItem *fpaItem = psMetadataLookup(fpa->hdu->header, keyword); - if (fpaItem) { - // XXX: Need to clean up before returning - return fpaItem; - } - } - - if (fpa->phu) { - psMetadataItem *fpaItem = psMetadataLookup(fpa->phu, keyword); - if (fpaItem) { - // XXX: Need to clean up before returning - return fpaItem; - } - } - } - - // No header value - return NULL; -} - - -// Look for a default -psMetadataItem *p_pmFPAConceptGetFromDefault(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - const char *concept // Name of concept - ) -{ - bool mdOK = true; // Status of MD lookup - psMetadata *defaults = psMetadataLookupMD(&mdOK, fpa->camera, "DEFAULTS"); - if (! mdOK) { - psError(PS_ERR_IO, false, "Unable to find DEFAULTS in camera configuration.\n"); - return NULL; - } - - psMetadataItem *defItem = psMetadataLookup(defaults, concept); - if (defItem) { - if (defItem->type == PS_DATA_METADATA) { - // A dependent default - psTrace(__func__, 7, "Evaluating dependent default....\n"); - psMetadata *dependents = defItem->data.V; // The list of dependents - // Find out what it depends on - psString dependName = psStringCopy(concept); - psStringAppend(&dependName, ".DEPEND"); - psString dependsOn = psMetadataLookupString(&mdOK, defaults, dependName); - if (! mdOK) { - psError(PS_ERR_IO, false, "Unable to find %s in camera configuration for dependent default" - " --- ignored\n", dependName); - // XXX: Need to clean up before returning - return NULL; - } - psFree(dependName); - // Find the value of the dependent concept - psMetadataItem *depItem = p_pmFPAConceptGetFromHeader(fpa, chip, cell, dependsOn); - if (! depItem) { - psError(PS_ERR_IO, true, "Unable to find value for %s (required for %s)\n", dependsOn, - concept); - return NULL; - } - if (depItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "Value of %s is not of type string, as required for dependency" - " --- ignored.\n", dependsOn); - } - - defItem = psMetadataLookup(dependents, depItem->data.V); // This is now what we were after - } - } - - // XXX: Need to clean up before returning - return defItem; // defItem is either NULL or points to what was desired -} - - -// Look for a database lookup -// XXX: Not tested -psMetadataItem *p_pmFPAConceptGetFromDB(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - psDB *db, // DB handle - const char *concept // Name of concept - ) -{ - if (! db) { - // No database initialised - return NULL; - } - - bool mdStatus = true; // Status of MD lookup - psMetadata *database = psMetadataLookupMD(&mdStatus, fpa->camera, "DATABASE"); - if (! mdStatus) { - // No error, because not everyone needs to use the DB - return NULL; - } - - psMetadata *dbLookup = psMetadataLookupMD(&mdStatus, database, concept); - if (dbLookup) { - const char *tableName = psMetadataLookupString(&mdStatus, dbLookup, "TABLE"); // Name of the table - const char *colName = psMetadataLookupString(&mdStatus, dbLookup, "COLUMN"); // Name of the column - const char *givenCols = psMetadataLookupString(&mdStatus, dbLookup, "GIVENDBCOL"); // Name of "where" - // columns - const char *givenPS = psMetadataLookupString(&mdStatus, dbLookup, "GIVENPS"); // Values for "where" - // columns - - // Now, need to get the "given"s - if (strlen(givenCols) || strlen(givenPS)) { - psList *cols = papSplit(givenCols, ",;"); // List of column names - psList *values = papSplit(givenPS, ",;"); // List of value names for the columns - psMetadata *selection = psMetadataAlloc(); // The stuff to select in the DB - if (cols->n != values->n) { - psLogMsg(__func__, PS_LOG_WARN, "The GIVENDBCOL and GIVENPS entries for %s do not have " - "the same number of entries --- ignored.\n", concept); - } else { - // Iterators for the lists - psListIterator *colsIter = psListIteratorAlloc(cols, PS_LIST_HEAD, false); - psListIterator *valuesIter = psListIteratorAlloc(values, PS_LIST_HEAD, false); - char *column = NULL; // Name of the column - while (column = psListGetAndIncrement(colsIter)) { - char *name = psListGetAndIncrement(valuesIter); // Name for the value - if (!strlen(column) || !strlen(name)) { - psLogMsg(__func__, PS_LOG_WARN, "One of the columns or value names for %s is " - " empty --- ignored.\n", concept); - } else { - // Search for the value name - psMetadataItem *item = p_pmFPAConceptGetFromHeader(fpa, chip, cell, name); - if (! item) { - item = p_pmFPAConceptGetFromDefault(fpa, chip, cell, name); - } - if (! item) { - psLogMsg(__func__, PS_LOG_ERROR, "Unable to find the value name %s for DB " - " lookup on %s --- ignored.\n", name, concept); - } else { - // We need to create a new psMetadataItem. I don't think we can't simply hack - // the existing one, since that could conceivably cause memory leaks - psMetadataItem *newItem = psMetadataItemAlloc(concept, item->type, - item->comment, item->data.V); - psMetadataAddItem(selection, newItem, PS_LIST_TAIL, PS_META_REPLACE); - psFree(newItem); - } - } - psFree(name); - psFree(column); - } // Iterating through the columns - psFree(colsIter); - psFree(valuesIter); - - psArray *dbResult = psDBSelectRows(db, tableName, selection, 2); // Lookup result - // Note that we use limit=2 in order to test if there are multiple rows returned - - psMetadataItem *result = NULL; // The final result of the DB lookup - if (dbResult->n == 0) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to find any rows in DB for %s --- ignored\n", - concept); - } else { - if (dbResult-> n > 1) { - psLogMsg(__func__, PS_LOG_WARN, "Multiple rows returned in DB lookup for %s --- " - " using the first one only.\n", concept); - } - result = (psMetadataItem*)dbResult->data[0]; - } - // XXX: Need to clean up before returning - return result; - } - psFree(cols); - psFree(values); - } - } // Doing the "given"s. - - psAbort(__func__, "Shouldn't ever get here.\n"); -} - - -// Concept lookup -psMetadataItem *p_pmFPAConceptGet(pmFPA *fpa, // The FPA - pmChip *chip,// The chip - pmCell *cell, // The cell - psDB *db, // DB handle - const char *concept // Concept name - ) -{ - // Try headers, database, defaults in order - psMetadataItem *item = p_pmFPAConceptGetFromCamera(cell, concept); - if (! item) { - item = p_pmFPAConceptGetFromHeader(fpa, chip, cell, concept); - } - if (! item) { - item = p_pmFPAConceptGetFromDB(fpa, chip, cell, db, concept); - } - if (! item) { - item = p_pmFPAConceptGetFromDefault(fpa, chip, cell, concept); - } - return item; // item is either NULL, or points to what was desired -} - - -void p_pmFPAConceptGetF32(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept - ) -{ - psMetadataItem *item = p_pmFPAConceptGet(fpa, chip, cell, db, name); - float value = NAN; - if (item) { - switch (item->type) { - case PS_DATA_F32: - value = item->data.F32; - break; - case PS_DATA_F64: - // Assume it's OK to truncate to floating point from double - value = (float)item->data.F64; - break; - case PS_DATA_S32: - // Promote to float - value = (float)item->data.S32; - break; - default: - psError(PS_ERR_IO, true, "Concept %s (%s) is not of floating point type (%x) --- treating as " - "undefined.\n", name, comment, item->type); - } - } else { - psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); - } - psTrace(__func__, 7, "Adding %s (%s): %f\n", name, comment, value); - - psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_F32 | PS_META_REPLACE, comment, value); -} - -void p_pmFPAConceptGetF64(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept - ) -{ - psMetadataItem *item = p_pmFPAConceptGet(fpa, chip, cell, db, name); - double value = NAN; - if (item) { - switch (item->type) { - case PS_TYPE_F64: - value = item->data.F64; - break; - case PS_TYPE_F32: - // Promote to double - value = (double)item->data.F32; - break; - case PS_TYPE_S32: - // Promote to double - value = (double)item->data.S32; - break; - default: - psError(PS_ERR_IO, true, "Concept %s (%s) is not of double-precision floating point type (%x) " - "--- treating as undefined.\n", name, comment, item->type); - } - } else { - psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); - } - psTrace(__func__, 7, "Adding %s (%s): %f\n", name, comment, value); - - psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_F64 | PS_META_REPLACE, comment, value); -} - -void p_pmFPAConceptGetS32(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept - ) -{ - psMetadataItem *item = p_pmFPAConceptGet(fpa, chip, cell, db, name); - int value = 0; - if (item) { - switch (item->type) { - case PS_TYPE_S32: - value = item->data.S32; - break; - case PS_TYPE_F32: - psLogMsg(__func__, PS_LOG_WARN, "Concept %s (%s) should be S32, but is F32 --- converting.\n", - name, comment); - value = (int)item->data.F32; - break; - case PS_TYPE_F64: - psLogMsg(__func__, PS_LOG_WARN, "Concept %s (%s) should be S32, but is F64 --- converting.\n", - name, comment); - value = (int)item->data.F64; - break; - default: - psError(PS_ERR_IO, true, "Concept %s (%s) is not of integer type (%x) --- treating as " - "undefined.\n", name, comment, item->type); - } - } else { - psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); - } - psTrace(__func__, 7, "Adding %s (%s): %d\n", name, comment, value); - - psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_S32 | PS_META_REPLACE, comment, value); -} - -void p_pmFPAConceptGetString(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept - ) -{ - psMetadataItem *item = p_pmFPAConceptGet(fpa, chip, cell, db, name); - psString value = NULL; - if (item) { - switch (item->type) { - case PS_DATA_STRING: - value = psMemIncrRefCounter(item->data.V); - break; - case PS_DATA_F32: - psStringAppend(&value, "%f", item->data.F32); - break; - case PS_DATA_S32: - psStringAppend(&value, "%d", item->data.S32); - break; - default: - psError(PS_ERR_IO, true, "Concept %s (%s) is not of string type (%x) --- treating as " - "undefined.\n", name, comment, item->type); - } - } else { - psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); - value = psStringCopy(""); - } - psTrace(__func__, 7, "Adding %s (%s): %s\n", name, comment, value); - - psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_DATA_STRING | PS_META_REPLACE, comment, value); - psFree(value); -} - - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Public functions -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -// Ingest concepts for the FPA -void pmFPAIngestConcepts(pmFPA *fpa, // The FPA - psDB *db // DB handle - ) -{ - if (! fpa->concepts) { - fpa->concepts = psMetadataAlloc(); - } - - // FPA.NAME - p_pmFPAConceptGetString(fpa, NULL, NULL, db, fpa->concepts, "FPA.NAME", "Name of FPA"); - - // FPA.AIRMASS - p_pmFPAConceptGetF32(fpa, NULL, NULL, db, fpa->concepts, "FPA.AIRMASS", "Airmass at boresight"); - - // FPA.FILTER - p_pmFPAConceptGetString(fpa, NULL, NULL, db, fpa->concepts, "FPA.FILTER", "Filter used"); - - // FPA.POSANGLE - p_pmFPAConceptGetF32(fpa, NULL, NULL, db, fpa->concepts, "FPA.POSANGLE", "Position angle for instrument"); - - // FPA.RADECSYS - p_pmFPAConceptGetString(fpa, NULL, NULL, db, fpa->concepts, "FPA.RADECSYS", "Celestial coordinate system"); - - // These take some extra work - - // FPA.RA - { - double ra = NAN; // The RA - psMetadataItem *raItem = p_pmFPAConceptGet(fpa, NULL, NULL, db, "FPA.RA"); // The FPA.RA item - if (raItem) { - switch (raItem->type) { - case PS_TYPE_F32: - ra = raItem->data.F32; - break; - case PS_TYPE_F64: - ra = raItem->data.F64; - break; - case PS_DATA_STRING: - // Sexagesimal format - { - int big, medium; - float small; - // XXX: Upgrade path is to allow dd:mm.mmm - if (sscanf(raItem->data.V, "%d:%d:%f", &big, &medium, &small) != 3 && - sscanf(raItem->data.V, "%d %d %f", &big, &medium, &small) != 3) { - psError(PS_ERR_IO, true, "Cannot interpret FPA.RA: %s\n", raItem->data.V); - break; - } - ra = abs(big) + (float)medium/60.0 + small/3600.0; - if (big < 0) { - ra *= -1.0; - } - } - break; - default: - psError(PS_ERR_IO, true, "FPA.RA is of an unexpected type: %x\n", raItem->type); - } - - // How to interpret the RA - const psMetadata *camera = fpa->camera; // Camera configuration data - bool mdok = true; // Status of MD lookup - psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); - if (mdok && formats) { - psString raFormat = psMetadataLookupString(&mdok, formats, "FPA.RA"); - if (mdok && strlen(raFormat) > 0) { - if (strcasecmp(raFormat, "HOURS") == 0) { - ra *= M_PI / 12.0; - } else if (strcasecmp(raFormat, "DEGREES") == 0) { - ra *= M_PI / 180.0; - } else if (strcasecmp(raFormat, "RADIANS") == 0) { - // No action required - } else { - psLogMsg(__func__, PS_LOG_WARN, "Don't understand FPA.RA in FORMATS (%s) --- assuming" - " HOURS.\n"); - ra *= M_PI / 12.0; - } - } else { - psError(PS_ERR_IO, false, "Unable to find FPA.RA in FORMATS --- assuming HOURS.\n"); - ra *= M_PI / 12.0; - } - } else { - psError(PS_ERR_IO, false, "Unable to find FORMAT metadata in camera configuration --- " - "assuming format for FPA.RA is HOURS.\n"); - ra *= M_PI / 12.0; - } - } else { - psError(PS_ERR_IO, false, "Couldn't find FPA.RA.\n"); - } - - psMetadataAdd(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_DATA_F64 | PS_META_REPLACE, - "Right Ascension of the boresight (radians)", ra); - - } - - // FPA.DEC - { - double dec = NAN; // The DEC - psMetadataItem *decItem = p_pmFPAConceptGet(fpa, NULL, NULL, db, "FPA.DEC"); // The FPA.DEC item - if (decItem) { - switch (decItem->type) { - case PS_TYPE_F32: - dec = decItem->data.F32; - break; - case PS_TYPE_F64: - dec = decItem->data.F64; - break; - case PS_DATA_STRING: - // Sexagesimal format - { - int big, medium; - float small; - // XXX: Upgrade path is to allow dd:mm.mmm - if (sscanf(decItem->data.V, "%d:%d:%f", &big, &medium, &small) != 3 && - sscanf(decItem->data.V, "%d %d %f", &big, &medium, &small) != 3) { - psError(PS_ERR_IO, true, "Cannot interpret FPA.DEC: %s\n", decItem->data.V); - break; - } - dec = abs(big) + (float)medium/60.0 + small/3600.0; - if (big < 0) { - dec *= -1.0; - } - } - break; - default: - psError(PS_ERR_IO, true, "FPA.DEC is of an unexpected type: %x\n", decItem->type); - } - - // How to interpret the DEC - const psMetadata *camera = fpa->camera; // Camera configuration data - bool mdok = true; // Status of MD lookup - psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); - if (mdok && formats) { - psString decFormat = psMetadataLookupString(&mdok, formats, "FPA.DEC"); - if (mdok && strlen(decFormat) > 0) { - if (strcasecmp(decFormat, "HOURS") == 0) { - dec *= M_PI / 12.0; - } else if (strcasecmp(decFormat, "DEGREES") == 0) { - dec *= M_PI / 180.0; - } else if (strcasecmp(decFormat, "RADIANS") == 0) { - // No action required - } else { - psLogMsg(__func__, PS_LOG_WARN, "Don't understand FPA.DEC in FORMATS (%s) --- " - "assuming DEGREES.\n"); - dec *= M_PI / 180.0; - } - } else { - psError(PS_ERR_IO, false, "Unable to find FPA.DEC in FORMATS --- assuming DEGREES.\n"); - dec *= M_PI / 180.0; - } - } else { - psError(PS_ERR_IO, false, "Unable to find FORMATS metadata in camera configuration --- " - "assuming format for FPA.DEC is DEGREES.\n"); - dec *= M_PI / 180.0; - } - } else { - psError(PS_ERR_IO, false, "Couldn't find FPA.DEC.\n"); - } - - psMetadataAdd(fpa->concepts, PS_LIST_TAIL, "FPA.DEC", PS_DATA_F64 | PS_META_REPLACE, - "Declination of the boresight (radians)", dec); - - } - - // Pau. -} - - -// Ingest concepts for the chip -bool pmChipIngestConcepts(pmChip *chip, // The chip - psDB *db // DB handle - ) -{ - pmFPA *fpa = chip->parent; // The parent FPA - - if (! chip->concepts) { - chip->concepts = psMetadataAlloc(); - } - - // CHIP.NAME --- added by pmFPAConstruct - - // Pau. -} - - -// Add corrective to a position --- in case the user wants FORTRAN indexing (like the FITS standard...) -static bool correctPosition(pmFPA *fpa, // FPA, contains the camera configuration - pmCell *cell, // Cell containing the concept to correct - const char *conceptName // Name of concept to correct - ) -{ - bool mdok = false; // Result of MD lookup - psMetadata *formats = psMetadataLookupMD(&mdok, fpa->camera, "FORMATS"); - if (mdok && formats) { - psString format = psMetadataLookupString(&mdok, formats, conceptName); - if (mdok && strlen(format) > 0 && strcasecmp(format, "FORTRAN") == 0) { - psMetadataItem *valueItem = psMetadataLookup(cell->concepts, conceptName); - valueItem->data.S32 -= 1; - } - } - return true; -} - -bool p_pmCellIngestConcept(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - const char *concept // Name of the concept - ) -{ - if (strcmp(concept, "CELL.GAIN") == 0) { - p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.GAIN", "CCD gain (e/count)"); - } else if (strcmp(concept, "CELL.READNOISE") == 0) { - p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.READNOISE", "CCD read noise (e)"); - } else if (strcmp(concept, "CELL.SATURATION") == 0) { - p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.SATURATION", - "Saturation level (ADU)"); - } else if (strcmp(concept, "CELL.BAD") == 0) { - p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.BAD", "Bad level (ADU)"); - } else if (strcmp(concept, "CELL.XPARITY") == 0) { - p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.XPARITY", - "Orientation in x compared to the rest of the FPA"); - } else if (strcmp(concept, "CELL.YPARITY") == 0) { - p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.YPARITY", - "Orientation in y compared to the rest of the FPA"); - } else if (strcmp(concept, "CELL.READDIR") == 0) { - p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.READDIR", - "Read direction: 1=row, 2=col"); - } else if (strcmp(concept, "CELL.EXPOSURE") == 0) { // used to be READOUT.EXPOSURE - p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.EXPOSURE", "Exposure time (sec)"); - } else if (strcmp(concept, "CELL.DARKTIME") == 0) { // used to be READOUT.DARKTIME - p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.DARKTIME", - "Time since CCD flush (sec)"); - // These take some extra work - } else if (strcmp(concept, "CELL.TRIMSEC") == 0) { - psRegion *trimsec = psAlloc(sizeof(psRegion)); // Make space for a psRegion (usually passed by value) - - psMetadataItem *secItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TRIMSEC"); - if (! secItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.TRIMSEC.\n"); - *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); - } else if (secItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "CELL.TRIMSEC is not of type STR (%x)\n", secItem->type); - *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); - } else { - psString section = secItem->data.V; // The section string - - psMetadataItem *sourceItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TRIMSEC.SOURCE"); - if (! sourceItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.TRIMSEC.SOURCE.\n"); - *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); - } else if (sourceItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "CELL.TRIMSEC.SOURCE is not of type STR (%x)\n", sourceItem->type); - *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); - } else { - psString source = sourceItem->data.V; // The source string - - if (strcasecmp(source, "VALUE") == 0) { - *trimsec = psRegionFromString(section); - } else if (strcasecmp(source, "HEADER") == 0) { - psMetadata *header = NULL; // The FITS header - if (cell->hdu) { - header = cell->hdu->header; - } else if (chip->hdu) { - header = chip->hdu->header; - } else if (fpa->hdu) { - header = fpa->hdu->header; - } - if (! header) { - psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); - *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); - } else { - bool mdok = true; // Status of MD lookup - psString secValue = psMetadataLookupString(&mdok, header, section); - if (! mdok || ! secValue) { - psError(PS_ERR_IO, false, "Unable to locate header %s\n", section); - *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); - } else { - *trimsec = psRegionFromString(secValue); - } - } - } else { - psError(PS_ERR_IO, true, "CELL.TRIMSEC.SOURCE (%s) is not HEADER or VALUE --- trying " - "VALUE.\n", source); - *trimsec = psRegionFromString(section); - } // Value of CELL.TRIMSEC.SOURCE - } // Looking up CELL.TRIMSEC.SOURCE - } // Looking up CELL.TRIMSEC - - psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TRIMSEC", PS_DATA_UNKNOWN | PS_META_REPLACE, - "Trim section", trimsec); - psFree(trimsec); - - } else if (strcmp(concept, "CELL.BIASSEC") == 0) { - psList *biassecs = psListAlloc(NULL); // List of bias sections - - psMetadataItem *secItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.BIASSEC"); - if (! secItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.BIASSEC.\n"); - } else if (secItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "CELL.BIASSEC is not of type STR (%x)\n", secItem->type); - } else { - psString sections = secItem->data.V; // The section string - - psMetadataItem *sourceItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.BIASSEC.SOURCE"); - if (! sourceItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.BIASSEC.SOURCE.\n"); - } else if (sourceItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "CELL.BIASSEC.SOURCE is not of type STR (%x)\n", sourceItem->type); - } else { - psString source = sourceItem->data.V; // The source string - - psList *secList = papSplit(sections, " ;"); // List of sections - psListIterator *secIter = psListIteratorAlloc(secList, PS_LIST_HEAD, false); // Iterator over - // sections - psString aSection = NULL; // A section from the list - while (aSection = psListGetAndIncrement(secIter)) { - psRegion *region = psAlloc(sizeof(psRegion)); // Make space for a psRegion (usually passed - // by value) - - if (strcasecmp(source, "VALUE") == 0) { - *region = psRegionFromString(aSection); - } else if (strcasecmp(source, "HEADER") == 0) { - psMetadata *header = NULL; // The FITS header - if (cell->hdu) { - header = cell->hdu->header; - } else if (chip->hdu) { - header = chip->hdu->header; - } else if (fpa->hdu) { - header = fpa->hdu->header; - } - if (! header) { - psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); - *region = psRegionSet(0.0,0.0,0.0,0.0); - } else { - bool mdok = true; // Status of MD lookup - psString secValue = psMetadataLookupString(&mdok, header, aSection); - if (! mdok || ! secValue) { - psError(PS_ERR_IO, false, "Unable to locate header %s\n", aSection); - *region = psRegionSet(0.0,0.0,0.0,0.0); - } else { - *region = psRegionFromString(secValue); - } - } - } else { - psError(PS_ERR_IO, true, "CELL.BIASSEC.SOURCE (%s) is not HEADER or VALUE --- trying " - "VALUE.\n", source); - *region = psRegionFromString(aSection); - } // Value of CELL.BIASSEC.SOURCE - - psListAdd(biassecs, PS_LIST_TAIL, region); - psFree(region); - } // Iterating over multiple sections - psFree(secIter); - psFree(secList); - } // Looking up CELL.BIASSEC.SOURCE - } // Looking up CELL.BIASSEC - - psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.BIASSEC", PS_DATA_LIST | PS_META_REPLACE, - "Bias sections", biassecs); - psFree(biassecs); - - } else if (strcmp(concept, "CELL.XBIN") == 0) { - int xBin = 1; // Binning factor in x - psMetadataItem *binItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.XBIN"); - if (! binItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.XBIN.\n"); - } else if (binItem->type == PS_DATA_STRING) { - psString binString = binItem->data.V; // The string containing the binning - if (sscanf(binString, "%d %*d", &xBin) != 1 && - sscanf(binString, "%d,%*d", &xBin) != 1) { - psError(PS_ERR_IO, true, "Unable to read string to get x binning: %s\n", binString); - } - } else if (binItem->type == PS_TYPE_S32) { - xBin = binItem->data.S32; - } else { - psError(PS_ERR_IO, true, "Note sure how to interpret CELL.XBIN of type %x --- assuming 1.\n", - binItem->type); - } - - psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_TYPE_S32 | PS_META_REPLACE, - "Binning in x", xBin); - - } else if (strcmp(concept, "CELL.YBIN") == 0) { - int yBin = 1; // Binning factor in y - psMetadataItem *binItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.YBIN"); - if (! binItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.YBIN.\n"); - } else if (binItem->type == PS_DATA_STRING) { - psString binString = binItem->data.V; // The string containing the binning - if (sscanf(binString, "%*d %d", &yBin) != 1 && - sscanf(binString, "%*d,%d", &yBin) != 1) { - psError(PS_ERR_IO, true, "Unable to read string to get y binning: %s\n", binString); - } - } else if (binItem->type == PS_TYPE_S32) { - yBin = binItem->data.S32; - } else { - psError(PS_ERR_IO, true, "Note sure how to interpret CELL.YBIN of type %x --- assuming 1.\n", - binItem->type); - } - - psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_TYPE_S32 | PS_META_REPLACE, - "Binning in y", yBin); - - } else if (strcmp(concept, "CELL.TIME") == 0 || strcmp(concept, "CELL.TIMESYS") == 0) { - // Do CELL.TIME and CELL.TIMESYS together - psTime *time = NULL; // The time - psTimeType timeSys = PS_TIME_UTC; // The time system - - // CELL.TIMESYS - psMetadataItem *sysItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TIMESYS"); - if (! sysItem) { - psError(PS_ERR_IO, true, "Couldn't find CELL.TIMESYS --- assuming UTC.\n"); - } else if (sysItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "CELL.TIMESYS isn't of type STRING --- assuming UTC.\n"); - } else { - psString sys = sysItem->data.V; // The time system string - if (strcasecmp(sys, "TAI") == 0) { - timeSys = PS_TIME_TAI; - } else if (strcasecmp(sys, "UTC") == 0) { - timeSys = PS_TIME_UTC; - } else if (strcasecmp(sys, "UT1") == 0) { - timeSys = PS_TIME_UT1; - } else if (strcasecmp(sys, "TT") == 0) { - timeSys = PS_TIME_TT; - } else { - psError(PS_ERR_IO, true, "Can't interpret CELL.TIMESYS --- assuming UTC.\n"); - } - } - psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TIMESYS", PS_TYPE_S32 | PS_META_REPLACE, - "Time system", timeSys); - - psMetadataItem *timeItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TIME"); - if (! timeItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.TIME.\n"); - } else { - // Get format - const psMetadata *camera = fpa->camera; // The camera configuration data - bool mdok = true; // Status of MD lookup - psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); - if (mdok && formats) { - psString timeFormat = psMetadataLookupString(&mdok, formats, "CELL.TIME"); - if (mdok && strlen(timeFormat) > 0) { - switch (timeItem->type) { - case PS_DATA_STRING: - { - psString timeString = timeItem->data.V; // String with the time - if (strcasecmp(timeFormat, "ISO") == 0) { - // timeString contains an ISO time - time = psTimeFromISO(timeString, timeSys); - } else if (strstr(timeFormat, "SEPARATE")) { - // timeString contains headers for the date and time - psMetadata *header = NULL; // The FITS header - if (cell->hdu) { - header = cell->hdu->header; - } else if (chip->hdu) { - header = chip->hdu->header; - } else if (fpa->hdu) { - header = fpa->hdu->header; - } - if (! header) { - psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); - } else { - // Get the headers - char *stuff1 = strpbrk(timeString, " ,;"); - psString dateName = psStringNCopy(timeString, - strlen(timeString) - strlen(stuff1)); - char *stuff2 = strpbrk(stuff1, "abcdefghijklmnopqrstuvwxyz" - "ABCDEFGHIJKLMNOPQRSTUVWXYZ"); - psString timeName = psStringCopy(stuff2); - - bool mdok = true; // Status of MD lookup - psString dateString = psMetadataLookupString(&mdok, header, dateName); - psFree(dateName); - int day = 0, month = 0, year = 0; - if (sscanf(dateString, "%d-%d-%d", &day, &month, &year) != 3 && - sscanf(dateString, "%d/%d/%d", &day, &month, &year) != 3) { - psError(PS_ERR_IO, true, "Unable to read date: %s\n", dateString); - } else { - if (strstr(timeFormat, "BACKWARDS")) { - int temp = day; - day = year; - year = temp; - } - if (strstr(timeFormat, "PRE2000") || year < 2000) { - year += 2000; - } - - psMetadataItem *timeItem = psMetadataLookup(header, timeName); - if (! timeItem) { - psError(PS_ERR_IO, false, "Unable to find time header: %s\n", - timeName); - } else if (timeItem->type == PS_DATA_STRING) { - // Time is a string, in the usual way: - psStringAppend(&dateString, "T%s", timeItem->data.V); - } else { - // Assume that time is specified in Second of Day - double seconds = NAN; - switch (timeItem->type) { - case PS_TYPE_S32: - seconds = timeItem->data.S32; - break; - case PS_TYPE_F32: - seconds = timeItem->data.F32; - break; - case PS_TYPE_F64: - seconds = timeItem->data.F64; - break; - default: - psError(PS_ERR_IO, true, "Time header (%s) is not of an " - "expected type: %x\n", timeName, timeItem->type); - } - // Now print to timeString as "hh:mm:ss.ss" - int hours = seconds / 3600; - seconds -= (double)hours * 3600.0; - int minutes = seconds / 60; - seconds -= (double)minutes * 60.0; - psStringAppend(&dateString, "T%02d:%02d:%02f", hours, minutes, - seconds); - } - time = psTimeFromISO(dateString, timeSys); - } // Reading date and time - psFree(timeName); - } // Reading headers - } else { - psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%s) --- trying " - "ISO\n", timeString); - time = psTimeFromISO(timeString, timeSys); - } // Interpreting the time string - } - break; - case PS_TYPE_F32: - { - double timeValue = (double)timeItem->data.F32; - if (strcasecmp(timeFormat, "JD") == 0) { - time = psTimeFromJD(timeValue); - } else if (strcasecmp(timeFormat, "MJD") == 0) { - time = psTimeFromMJD(timeValue); - } else { - psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%f) --- trying " - "JD\n", timeValue); - time = psTimeFromJD(timeValue); - } - } - break; - case PS_TYPE_F64: - { - double timeValue = (double)timeItem->data.F64; - if (strcasecmp(timeFormat, "JD") == 0) { - time = psTimeFromJD(timeValue); - } else if (strcasecmp(timeFormat, "MJD") == 0) { - time = psTimeFromMJD(timeValue); - } else { - psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%f) --- trying " - "JD\n", timeValue); - time = psTimeFromJD(timeValue); - } - } - break; - default: - psError(PS_ERR_IO, true, "Unable to parse CELL.TIME.\n"); - } - } else { - psError(PS_ERR_IO, false, "Unable to find CELL.TIME in FORMATS.\n"); - } // Getting the format - } else { - psError(PS_ERR_IO, false, "Unable to find FORMATS in camera configuration.\n"); - } // Getting the formats - } // Getting CELL.TIME - - psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TIME", PS_DATA_TIME | PS_META_REPLACE, - "Time of exposure", time); - psFree(time); - - // These are new and experimental concepts: CELL.X0 and CELL.Y0 - } else if (strcmp(concept, "CELL.X0") == 0) { - p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.X0", "Position of (0,0) on the chip"); - correctPosition(fpa, cell, "CELL.X0"); - } else if (strcmp(concept, "CELL.Y0") == 0) { - p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.Y0", "Position of (0,0) on the chip"); - correctPosition(fpa, cell, "CELL.Y0"); - } - - return true; -} - - - -// Ingest concepts for the cell -bool pmCellIngestConcepts(pmCell *cell, // The cell - psDB *db // DB handle - ) -{ - pmChip *chip = cell->parent; // The parent chip - pmFPA *fpa = chip->parent; // The parent FPA - - if (! cell->concepts) { - cell->concepts = psMetadataAlloc(); - } - - // CELL.NAME --- added by pmFPAConstruct - - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.GAIN"); - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.READNOISE"); - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.SATURATION"); - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.BAD"); - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.XPARITY"); - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.YPARITY"); - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.READDIR"); - - // These used to be pmReadoutGetExposure and pmReadoutGetDarkTime, but that doesn't really make sense at - // the moment. Maybe we need to add a "parent" link to the readouts. But then how are the exposure times - // REALLY derived? They're not in the FITS headers, because a readout is a plane in a 3D image. We'll - // have to dream up some additional suffix to specify these, but for now.... - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.EXPOSURE"); - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.DARKTIME"); - - - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.TRIMSEC"); - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.BIASSEC"); - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.TIME"); // This also does CELL.TIMESYS - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.XBIN"); - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.YBIN"); - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.X0"); - p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.Y0"); - -} - - -// Retrieve a concept -psMetadataItem *pmCellGetConcept(pmCell *cell, // The cell - const char *concept // The concept - ) -{ - psMetadataItem *item = psMetadataLookup(cell->concepts, concept); - if (! item) { - pmChip *chip = cell->parent; - item = psMetadataLookup(chip->concepts, concept); - if (! item) { - pmFPA *fpa = chip->parent; - item = psMetadataLookup(fpa->concepts, concept); - } - } - return item; // item is either NULL or is what we want. -} Index: unk/archive/scripts/src/phase2/pmFPAConceptsGet.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAConceptsGet.h (revision 5798) +++ (revision ) @@ -1,94 +1,0 @@ -#ifndef PM_FPA_CONCEPTS_GET_H -#define PM_FPA_CONCEPTS_GET_H - -#include "pmFPA.h" -#include "pslib.h" -#include "psAdditionals.h" - -psMetadataItem *p_pmFPAConceptGetFromCamera(pmCell *cell, // The cell - const char *concept // Name of concept - ); -psMetadataItem *p_pmFPAConceptGetFromHeader(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - const char *concept // Name of concept - ); -psMetadataItem *p_pmFPAConceptGetFromDefault(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - const char *concept // Name of concept - ); -psMetadataItem *p_pmFPAConceptGetFromDB(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - psDB *db, // DB handle - const char *concept // Name of concept - ); -psMetadataItem *p_pmFPAConceptGet(pmFPA *fpa, // The FPA - pmChip *chip,// The chip - pmCell *cell, // The cell - psDB *db, // DB handle - const char *concept // Concept name - ); -void p_pmFPAConceptGetF32(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept - ); -void p_pmFPAConceptGetF64(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept - ); -void p_pmFPAConceptGetS32(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept - ); -void p_pmFPAConceptGetString(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept - ); - - -// Ingest a single cell concept; this is the workhorse for ingesting cell concepts -bool p_pmCellIngestConcept(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - const char *concept // Name of the concept - ); - -// Ingest concepts for the FPA -void pmFPAIngestConcepts(pmFPA *fpa, // The FPA - psDB *db // DB handle - ); -// Ingest concepts for the chip -bool pmChipIngestConcepts(pmChip *chip, // The chip - psDB *db // DB handle - ); -// Ingest concepts for the cell -bool pmCellIngestConcepts(pmCell *cell, // The cell - psDB *db // DB handle - ); - -// Retrieve a concept -psMetadataItem *pmCellGetConcept(pmCell *cell, // The cell - const char *concept // The concept - ); - - -#endif Index: unk/archive/scripts/src/phase2/pmFPAConceptsSet.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAConceptsSet.c (revision 5798) +++ (revision ) @@ -1,887 +1,0 @@ -#include -#include -#include "pslib.h" - -#include "pmFPA.h" - -#include "papStuff.h" - -#include "pmFPAConceptsGet.h" -#include "pmFPAConceptsSet.h" - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// File-static functions -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -#define COMPARE_REGIONS(a,b) (((a)->x0 == (b)->x0 && \ - (a)->x1 == (b)->x1 && \ - (a)->y0 == (b)->y0 && \ - (a)->y1 == (b)->y1) ? true : false) - - -static bool compareConcepts(psMetadataItem *item1, // First item to compare - psMetadataItem *item2 // Second item to compare - ) -{ - // First order checks - if (! item1 || ! item2) { - return false; - } - if (strcasecmp(item1->name, item2->name) != 0) { - return false; - } - - // Check the more boring types - switch (item1->type) { - case PS_TYPE_S32: - switch (item2->type) { - case PS_TYPE_S32: - return (item1->data.S32 == item2->data.S32) ? true : false; - case PS_TYPE_F32: - return (item1->data.S32 == (int)item2->data.F32) ? true : false; - case PS_TYPE_F64: - return (item1->data.S32 == (int)item2->data.F64) ? true : false; - default: - return false; - } - case PS_TYPE_F32: - switch (item2->type) { - case PS_TYPE_S32: - return (item1->data.F32 == (float)item2->data.S32) ? true : false; - case PS_TYPE_F32: - return (item1->data.F32 == item2->data.F32) ? true : false; - case PS_TYPE_F64: - return (item1->data.F32 == (float)item2->data.F64) ? true : false; - default: - return false; - } - case PS_TYPE_F64: - switch (item2->type) { - case PS_TYPE_S32: - return (item1->data.F64 == (double)item2->data.S32) ? true : false; - case PS_TYPE_F32: - return (item1->data.F64 == (double)item2->data.F32) ? true : false; - case PS_TYPE_F64: - return (item1->data.F64 == item2->data.F64) ? true : false; - default: - return false; - } - return (item1->data.F64 == item2->data.F64) ? true : false; - case PS_DATA_STRING: - if (item2->type != PS_DATA_STRING) { - return false; - } - return (strcasecmp(item1->data.V, item2->data.V) == 0) ? true : false; - default: - return false; - } - psAbort(__func__, "Should never get here.\n"); -} - - -// Well, not really "set", but check to make sure it's there and matches -static bool setConceptInCamera(pmCell *cell, // The cell - psMetadataItem *concept // Concept - ) -{ - if (cell) { - psMetadataItem *item = psMetadataLookup(cell->camera, concept->name); // Info we want - return compareConcepts(item, concept); - } -} - - -static bool setConceptInHeader(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - psMetadataItem *concept // Concept - ) -{ - bool mdok = true; // Status of MD lookup - bool status = false; // Status of setting header - psMetadata *translation = psMetadataLookupMD(&mdok, fpa->camera, "TRANSLATION"); // FITS translation - if (! mdok) { - psError(PS_ERR_IO, false, "Unable to find TRANSLATION in camera configuration.\n"); - return false; - } - - // Look for how to translate the concept into a FITS header name - const char *keyword = psMetadataLookupString(&mdok, translation, concept->name); - if (mdok && strlen(keyword) > 0) { - psTrace(__func__, 6, "It's in keyword %s\n", keyword); - psMetadataItem *headerItem = NULL; // Item to add to header - // XXX: Need to expand range of types - switch (concept->type) { - case PS_DATA_STRING: - headerItem = psMetadataItemAllocStr(keyword, concept->comment, concept->data.V); - break; - case PS_DATA_S32: - headerItem = psMetadataItemAllocS32(keyword, concept->comment, concept->data.S32); - break; - case PS_DATA_F32: - headerItem = psMetadataItemAllocF32(keyword, concept->comment, concept->data.F32); - break; - case PS_DATA_F64: - headerItem = psMetadataItemAllocF64(keyword, concept->comment, concept->data.F64); - break; - default: - headerItem = psMetadataItemAlloc(keyword, concept->type, concept->comment, - concept->data.V); // Item for the header - } - - // We have a FITS header to look up --- search each level - if (cell && cell->hdu) { - psTrace(__func__, 7, "Adding to the cell level header...\n"); - psMetadataAddItem(cell->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); - status = true; - } else if (chip && chip->hdu) { - psTrace(__func__, 7, "Adding to the chip level header...\n"); - psMetadataAddItem(chip->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); - status = true; - } else if (fpa->hdu) { - psTrace(__func__, 7, "Adding to the FPA level header...\n"); - psMetadataAddItem(fpa->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); - status = true; - } else { - // In desperation, add to the PHU --- it HAS to be in the header somewhere! - if (! fpa->phu) { - fpa->phu = psMetadataAlloc(); - } - psTrace(__func__, 7, "Adding to the PHU...\n"); - psMetadataAddItem(fpa->phu, headerItem, PS_LIST_TAIL, PS_META_REPLACE); - status = true; - } - psFree(headerItem); - } - - // No header value - return status; -} - - -// Well, not really "set", but check to see if it's there, and matches -static bool setConceptInDefault(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - psMetadataItem *concept // Concept - ) -{ - bool mdOK = true; // Status of MD lookup - psMetadata *defaults = psMetadataLookupMD(&mdOK, fpa->camera, "DEFAULTS"); - if (! mdOK || ! defaults) { - psError(PS_ERR_IO, false, "Unable to find DEFAULTS in camera configuration.\n"); - return false; - } - - psMetadataItem *defItem = psMetadataLookup(defaults, concept->name); - bool status = false; // Result of checking the database - if (defItem) { - if (defItem->type == PS_DATA_METADATA) { - // A dependent default - psTrace(__func__, 7, "Evaluating dependent default....\n"); - psMetadata *dependents = defItem->data.V; // The list of dependents - // Find out what it depends on - psString dependName = psStringCopy(concept->name); - psStringAppend(&dependName, ".DEPEND"); - psString dependsOn = psMetadataLookupString(&mdOK, defaults, dependName); - if (! mdOK) { - psError(PS_ERR_IO, false, "Unable to find %s in camera configuration for dependent default" - " --- ignored\n", dependName); - // XXX: Need to clean up before returning - return false; - } - psFree(dependName); - // Find the value of the dependent concept - psMetadataItem *depItem = p_pmFPAConceptGetFromHeader(fpa, chip, cell, dependsOn); - if (! depItem) { - psError(PS_ERR_IO, true, "Unable to find value for %s (required for %s)\n", dependsOn, - concept->name); - return false; - } - if (depItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "Value of %s is not of type string, as required for dependency" - " --- ignored.\n", dependsOn); - } - - defItem = psMetadataLookup(dependents, depItem->data.V); // This is now what we were after - } - - status = compareConcepts(defItem, concept); - if (! status) { - psError(PS_ERR_IO, true, "Concept %s is specified by default in the camera configuration, " - "but doesn't match the actual value.\n", concept->name); - } - } - - // XXX: Need to clean up before returning - return status; -} - - -// XXX: Not tested at all -static bool setConceptInDB(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadataItem *concept // Concept - ) -{ - if (! db) { - // No database initialised - return false; - } - - bool mdStatus = true; // Status of MD lookup - psMetadata *database = psMetadataLookupMD(&mdStatus, fpa->camera, "DATABASE"); - if (! mdStatus) { - // No error, because not everyone needs to use the DB - return NULL; - } - - psMetadata *dbLookup = psMetadataLookupMD(&mdStatus, database, concept->name); - if (dbLookup) { - const char *tableName = psMetadataLookupString(&mdStatus, dbLookup, "TABLE"); // Name of the table - const char *colName = psMetadataLookupString(&mdStatus, dbLookup, "COLUMN"); // Name of the column - const char *givenCols = psMetadataLookupString(&mdStatus, dbLookup, "GIVENDBCOL"); // Name of "where" - // columns - const char *givenPS = psMetadataLookupString(&mdStatus, dbLookup, "GIVENPS"); // Values for "where" - // columns - - // Now, need to get the "given"s - if (strlen(givenCols) || strlen(givenPS)) { - psList *cols = papSplit(givenCols, ",;"); // List of column names - psList *values = papSplit(givenPS, ",;"); // List of value names for the columns - psMetadata *selection = psMetadataAlloc(); // The stuff to select in the DB - if (cols->n != values->n) { - psLogMsg(__func__, PS_LOG_WARN, "The GIVENDBCOL and GIVENPS entries for %s do not have " - "the same number of entries --- ignored.\n", concept); - } else { - // Iterators for the lists - psListIterator *colsIter = psListIteratorAlloc(cols, PS_LIST_HEAD, false); - psListIterator *valuesIter = psListIteratorAlloc(values, PS_LIST_HEAD, false); - char *column = NULL; // Name of the column - while (column = psListGetAndIncrement(colsIter)) { - char *name = psListGetAndIncrement(valuesIter); // Name for the value - if (!strlen(column) || !strlen(name)) { - psLogMsg(__func__, PS_LOG_WARN, "One of the columns or value names for %s is " - " empty --- ignored.\n", concept); - } else { - // Search for the value name - psMetadataItem *item = p_pmFPAConceptGetFromHeader(fpa, chip, cell, name); - if (! item) { - item = p_pmFPAConceptGetFromDefault(fpa, chip, cell, name); - } - if (! item) { - psLogMsg(__func__, PS_LOG_ERROR, "Unable to find the value name %s for DB " - " lookup on %s --- ignored.\n", name, concept); - } else { - // We need to create a new psMetadataItem. I don't think we can't simply hack - // the existing one, since that could conceivably cause memory leaks - psMetadataItem *newItem = psMetadataItemAlloc(concept->name, item->type, - item->comment, item->data.V); - psMetadataAddItem(selection, newItem, PS_LIST_TAIL, PS_META_REPLACE); - psFree(newItem); - } - } - psFree(name); - psFree(column); - } // Iterating through the columns - psFree(colsIter); - psFree(valuesIter); - - // Check first to make sure we're only going to touch one row - psArray *dbResult = psDBSelectRows(db, tableName, selection, 2); // Lookup result - // Note that we use limit=2 in order to test if there are multiple rows returned - - psMetadataItem *result = NULL; // The final result of the DB lookup - if (! dbResult || dbResult->n == 0) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to find any rows in DB for %s --- ignored\n", - concept->name); - return false; - } else { - if (dbResult->n > 1) { - psLogMsg(__func__, PS_LOG_WARN, "Multiple rows returned in DB lookup for %s --- " - " ignored.\n", concept->name); - } - // Update the DB - psMetadata *update = psMetadataAlloc(); - psMetadataAddItem(update, concept, PS_LIST_HEAD, 0); - psDBUpdateRows(db, tableName, selection, update); - psFree(update); - return true; - } - } - psFree(cols); - psFree(values); - } - } // Doing the "given"s. - - psAbort(__func__, "Shouldn't ever get here?\n"); -} - - -// Concept set from item -static bool setConceptItem(pmFPA *fpa, // The FPA - pmChip *chip,// The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadataItem *concept // Concept item - ) -{ - // Try headers, database, defaults in order - psTrace(__func__, 3, "Trying to set concept %s...\n", concept->name); - bool status = setConceptInCamera(cell, concept); // Status for return - if (! status) { - psTrace(__func__, 5, "Trying header....\n"); - status = setConceptInHeader(fpa, chip, cell, concept); - } - if (! status) { - psTrace(__func__, 5, "Trying database....\n"); - status = setConceptInDB(fpa, chip, cell, db, concept); - } - if (! status) { - psTrace(__func__, 5, "Checking defaults....\n"); - status = setConceptInDefault(fpa, chip, cell, concept); - } - - if (! status) { - psError(PS_ERR_IO, true, "Unable to set %s (%s).\n", concept->name, concept->comment); - } - - return status; -} - - -// Concept set -static bool setConcept(pmFPA *fpa, // The FPA - pmChip *chip,// The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // Concepts MD from which to set - const char *name // Name of the concept - ) -{ - psMetadataItem *concept = psMetadataLookup(concepts, name); - if (! concept) { - psError(PS_ERR_IO, true, "No such concept as %s\n", name); - return false; - } - return setConceptItem(fpa, chip, cell, db, concept); -} - - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Public functions -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -// Outget concepts from the FPA -void pmFPAOutgestConcepts(pmFPA *fpa, // The FPA - psDB *db // DB handle - ) -{ - // FPA.NAME - setConcept(fpa, NULL, NULL, db, fpa->concepts, "FPA.NAME"); - - // FPA.AIRMASS - setConcept(fpa, NULL, NULL, db, fpa->concepts, "FPA.AIRMASS"); - - // FPA.FILTER - setConcept(fpa, NULL, NULL, db, fpa->concepts, "FPA.FILTER"); - - // FPA.POSANGLE - setConcept(fpa, NULL, NULL, db, fpa->concepts, "FPA.POSANGLE"); - - // FPA.RADECSYS - setConcept(fpa, NULL, NULL, db, fpa->concepts, "FPA.RADECSYS"); - - // These take some extra work - - // FPA.RA - { - double ra = psMetadataLookupF64(NULL, fpa->concepts, "FPA.RA"); // The RA - - // How to interpret the RA - const psMetadata *camera = fpa->camera; // Camera configuration data - bool mdok = true; // Status of MD lookup - psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); - if (mdok && formats) { - psString raFormat = psMetadataLookupString(&mdok, formats, "FPA.RA"); - if (mdok && strlen(raFormat) > 0) { - if (strcasecmp(raFormat, "HOURS") == 0) { - ra /= M_PI / 12.0; - } else if (strcasecmp(raFormat, "DEGREES") == 0) { - ra /= M_PI / 180.0; - } else if (strcasecmp(raFormat, "RADIANS") == 0) { - // No action required - } else { - psLogMsg(__func__, PS_LOG_WARN, "Don't understand FPA.RA in FORMATS (%s) --- assuming" - " HOURS.\n"); - ra /= M_PI / 12.0; - } - } else { - psError(PS_ERR_IO, false, "Unable to find FPA.RA in FORMATS --- assuming HOURS.\n"); - ra /= M_PI / 12.0; - } - } else { - psError(PS_ERR_IO, false, "Unable to find FORMAT metadata in camera configuration --- " - "assuming format for FPA.RA is HOURS.\n"); - ra /= M_PI / 12.0; - } - - // We choose to write sexagesimal format - int big, medium; - float small; - big = (int)ra; - medium = (int)(60.0*(ra - (double)big)); - small = 3600.0*(ra - (double)big - 60.0 * (double)medium); - psString raString = psStringCopy(""); - psStringAppend(&raString, "%d:%d:%.2f", big, medium, small); - psMetadataItem *raItem = psMetadataItemAllocStr("FPA.RA", "Right Ascension of the boresight", - raString); - setConceptItem(fpa, NULL, NULL, db, raItem); - psFree(raItem); - psFree(raString); - } - - // FPA.DEC - { - double dec = psMetadataLookupF64(NULL, fpa->concepts, "FPA.DEC"); // The DEC - - // How to interpret the DEC - const psMetadata *camera = fpa->camera; // Camera configuration data - bool mdok = true; // Status of MD lookup - psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); - if (mdok && formats) { - psString decFormat = psMetadataLookupString(&mdok, formats, "FPA.DEC"); - if (mdok && strlen(decFormat) > 0) { - if (strcasecmp(decFormat, "HOURS") == 0) { - dec /= M_PI / 12.0; - } else if (strcasecmp(decFormat, "DEGREES") == 0) { - dec /= M_PI / 180.0; - } else if (strcasecmp(decFormat, "RADIANS") == 0) { - // No action required - } else { - psLogMsg(__func__, PS_LOG_WARN, "Don't understand FPA.DEC in FORMATS (%s) --- assuming" - " DEGREES.\n"); - dec /= M_PI / 180.0; - } - } else { - psError(PS_ERR_IO, false, "Unable to find FPA.DEC in FORMATS --- assuming DEGREES.\n"); - dec /= M_PI / 180.0; - } - } else { - psError(PS_ERR_IO, false, "Unable to find FORMAT metadata in camera configuration --- " - "assuming format for FPA.DEC is HOURS.\n"); - dec /= M_PI / 12.0; - } - - // We choose to write sexagesimal format - int big, medium; - float small; - big = (int)dec; - medium = (int)(60.0*(dec - (double)big)); - small = 3600.0*(dec - (double)big - 60.0 * (double)medium); - psString decString = psStringCopy(""); - psStringAppend(&decString, "%d:%d:%.2f", big, medium, small); - psMetadataItem *decItem = psMetadataItemAllocStr("FPA.DEC", "Right Ascension of the boresight", - decString); - setConceptItem(fpa, NULL, NULL, db, decItem); - psFree(decItem); - psFree(decString); - - } - - // Pau. -} - - -// Outgest concepts for the chip -void pmChipOutgestConcepts(pmChip *chip, // The chip - psDB *db // DB handle - ) -{ - pmFPA *fpa = chip->parent; // The parent FPA - - // CHIP.NAME --- no need to do anything - - // Pau. -} - - -// Outgest concepts for the cell -void pmCellOutgestConcepts(pmCell *cell, // The cell - psDB *db // DB handle - ) -{ - pmChip *chip = cell->parent; // The parent chip - pmFPA *fpa = chip->parent; // The parent FPA - - // CELL.NAME --- no need to do anything - - // CELL.GAIN - setConcept(fpa, chip, cell, db, cell->concepts, "CELL.GAIN"); - - // CELL.READNOISE - setConcept(fpa, chip, cell, db, cell->concepts, "CELL.READNOISE"); - - // CELL.SATURATION - setConcept(fpa, chip, cell, db, cell->concepts, "CELL.SATURATION"); - - // CELL.BAD - setConcept(fpa, chip, cell, db, cell->concepts, "CELL.BAD"); - - // CELL.XPARITY - setConcept(fpa, chip, cell, db, cell->concepts, "CELL.XPARITY"); - - // CELL.YPARITY - setConcept(fpa, chip, cell, db, cell->concepts, "CELL.YPARITY"); - - // CELL.READDIR - setConcept(fpa, chip, cell, db, cell->concepts, "CELL.READDIR"); - - // These used to be pmReadoutGetExposure and pmReadoutGetDarkTime, but that doesn't really make sense at - // the moment. Maybe we need to add a "parent" link to the readouts. But then how are the exposure times - // REALLY derived? They're not in the FITS headers, because a readout is a plane in a 3D image. We'll - // have to dream up some additional suffix to specify these, but for now.... - - // CELL.EXPOSURE (used to be READOUT.EXPOSURE) - setConcept(fpa, chip, cell, db, cell->concepts, "CELL.EXPOSURE"); - - // CELL.DARKTIME (used to be READOUT.DARKTIME) - setConcept(fpa, chip, cell, db, cell->concepts, "CELL.DARKTIME"); - - // These take some extra work - - // CELL.TRIMSEC - { - psMetadataItem *trimsecItem = psMetadataLookup(cell->concepts, "CELL.TRIMSEC"); - psRegion *trimsec = trimsecItem->data.V; // The trimsec region - psMetadataItem *sourceItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TRIMSEC.SOURCE"); - if (! sourceItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.TRIMSEC.SOURCE.\n"); - } else if (sourceItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "CELL.TRIMSEC.SOURCE is not of type STR (%x)\n", sourceItem->type); - } else { - psString source = sourceItem->data.V; // The source string - - if (strcasecmp(source, "VALUE") == 0) { - // Check that it's the same value as stored in the camera - psString checkString = psMetadataLookupString(NULL, cell->camera, "CELL.TRIMSEC"); - psRegion checkRegion = psRegionFromString(checkString); - if (! COMPARE_REGIONS(&checkRegion, trimsec)) { - psError(PS_ERR_IO, true, "Target CELL.TRIMSEC is specified by value, and values don't " - "match.\n"); - } - } else if (strcasecmp(source, "HEADER") == 0) { - psString keyword = psMetadataLookupString(NULL, cell->camera, "CELL.TRIMSEC"); - psMetadata *header = NULL; // The FITS header - if (cell->hdu) { - header = cell->hdu->header; - } else if (chip->hdu) { - header = chip->hdu->header; - } else if (fpa->hdu) { - header = fpa->hdu->header; - } - if (! header) { - psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); - } else { - psMetadataAddItem(header, trimsecItem, PS_LIST_TAIL, PS_META_REPLACE); - } - } else { - psError(PS_ERR_IO, true, "CELL.TRIMSEC.SOURCE (%s) is not HEADER or VALUE.\n", source); - } - } - } - - // CELL.BIASSEC - { - psMetadataItem *biassecItem = psMetadataLookup(cell->concepts, "CELL.BIASSEC"); - psList *biassecs = biassecItem->data.V; // The biassecs region list - psMetadataItem *sourceItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.BIASSEC.SOURCE"); - if (! sourceItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.BIASSEC.SOURCE.\n"); - } else if (sourceItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "CELL.BIASSEC.SOURCE is not of type STR (%x)\n", sourceItem->type); - } else { - psString source = sourceItem->data.V; // The source string - - if (strcasecmp(source, "VALUE") == 0) { - // Check that it's the same value as stored in the camera - psString checkString = psMetadataLookupString(NULL, cell->camera, "CELL.BIASSEC"); - psList *checkList = papSplit(checkString, " ;"); - if (biassecs->n != checkList->n) { - psError(PS_ERR_IO, true, "Target CELL.BIASSEC is specified by value, but number of " - "entries doesn't match.\n"); - } else { - // We don't care if the order matches or not - psVector *check = psVectorAlloc(biassecs->n, PS_TYPE_U8); // Vector to mark regions off - for (int i = 0; i < check->n; i++) { - check->data.U8[i] = 0; - } - psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, - false); // Iterator - psListIterator *checkListIter = psListIteratorAlloc(checkList, PS_LIST_HEAD, false); - psRegion *biassec = NULL; // Region from iteration - while (biassec = psListGetAndIncrement(biassecsIter)) { - psListIteratorSet(checkListIter, PS_LIST_HEAD); - psString checkRegionString = NULL; // Region string from iteration - int i = 0; // Counter - while (checkRegionString = psListGetAndIncrement(checkListIter)) { - psRegion checkRegion = psRegionFromString(checkRegionString); - psTrace(__func__, 7, "Checking [%.0f:%.0f,%.0f:%.0f] against " - "[%.0f:%.0f,%.0f:%.0f]\n", biassec->x0, biassec->x1, biassec->y0, - biassec->y1, checkRegion.x0, checkRegion.x1, checkRegion.y0, - checkRegion.y1); - if (COMPARE_REGIONS(biassec, &checkRegion)) { - check->data.U8[i] = 1; - i++; - break; - } - i++; - } - } - for (int i = 0; i < check->n; i++) { - if (check->data.U8[i] == 0) { - psError(PS_ERR_IO, true, "Target CELL.BIASSEC is specified by value, but values " - "don't match.\n"); - break; - } - } - psFree(checkListIter); - psFree(biassecsIter); - psFree(check); - } - psFree(checkList); - } else if (strcasecmp(source, "HEADER") == 0) { - psString keywordsString = psMetadataLookupString(NULL, cell->camera, "CELL.BIASSEC"); - psList *keywords = papSplit(keywordsString, " ,;"); - if (biassecs->n != keywords->n) { - psError(PS_ERR_IO, true, "Target CELL.BIASSEC is sepcified by headers, but the number " - "of headers doesn't match.\n"); - } else { - psMetadata *header = NULL; // The FITS header - if (cell->hdu) { - header = cell->hdu->header; - } else if (chip->hdu) { - header = chip->hdu->header; - } else if (fpa->hdu) { - header = fpa->hdu->header; - } - if (! header) { - psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); - } else { - psListIterator *keywordsIter = psListIteratorAlloc(keywords, PS_LIST_HEAD, false); - psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); - psString keyword = NULL; // Header keyword from list - while (keyword = psListGetAndIncrement(keywordsIter)) { - // Update the header - psRegion *biassec = psListGetAndIncrement(biassecsIter); - psString biassecString = psRegionToString(*biassec); - psMetadataAdd(header, PS_LIST_TAIL, keyword, PS_DATA_STRING | PS_META_REPLACE, - "Bias section", biassecString); - psFree(biassecString); - } - psFree(keywordsIter); - psFree(biassecsIter); - } - } - psFree(keywords); - } else { - psError(PS_ERR_IO, true, "CELL.BIASSEC.SOURCE (%s) is not HEADER or VALUE.\n", source); - } - } - } - - // CELL.XBIN and CELL.YBIN - { - psMetadataItem *xBinItem = psMetadataLookup(cell->concepts, "CELL.XBIN"); // Binning factor in x - psMetadataItem *yBinItem = psMetadataLookup(cell->concepts, "CELL.YBIN"); // Binning factor in y - - // If these are specified by the header, we need to check to see if it's the same header keyword - const psMetadata *camera = fpa->camera; - psMetadata *translation = psMetadataLookupMD(NULL, camera, "TRANSLATION"); - psString xKeyword = psMetadataLookupString(NULL, translation, "CELL.XBIN"); - psString yKeyword = psMetadataLookupString(NULL, translation, "CELL.YBIN"); - if (strlen(xKeyword) > 0 && strcasecmp(xKeyword, yKeyword) != 0) { - setConceptInHeader(fpa, chip, cell, xBinItem); - xBinItem = NULL; - } - if (strlen(yKeyword) > 0 && strcasecmp(xKeyword, yKeyword) != 0) { - setConceptInHeader(fpa, chip, cell, yBinItem); - yBinItem = NULL; - } - if (strlen(xKeyword) > 0 && strlen(yKeyword) > 0 && strcasecmp(xKeyword, yKeyword) == 0) { - psString binString = psStringCopy(""); - psStringAppend(&binString, "%d,%d", xBinItem->data.S32, yBinItem->data.S32); - psMetadataItem *binItem = psMetadataItemAllocStr(xKeyword, "Binning factor in x and y", - binString); - psFree(binString); - psMetadata *header = NULL; // The FITS header - if (cell->hdu) { - header = cell->hdu->header; - } else if (chip->hdu) { - header = chip->hdu->header; - } else if (fpa->hdu) { - header = fpa->hdu->header; - } - if (! header) { - psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); - } else { - psMetadataAddItem(header, binItem, PS_LIST_TAIL, PS_META_REPLACE); - xBinItem = NULL; - yBinItem = NULL; - } - psFree(binItem); - } - - // Do it in the usual way if we have to - if (xBinItem) { - setConceptItem(fpa, chip, cell, db, xBinItem); - } - if (yBinItem) { - setConceptItem(fpa, chip, cell, db, yBinItem); - } - } - - // CELL.TIMESYS - { - psMetadataItem *sysItem = psMetadataLookup(cell->concepts, "CELL.TIMESYS"); - psString sys = NULL; // String to store - switch (sysItem->data.S32) { - case PS_TIME_TAI: - sys = psStringCopy("TAI"); - break; - case PS_TIME_UTC: - sys = psStringCopy("UTC"); - break; - case PS_TIME_UT1: - sys = psStringCopy("UT1"); - break; - case PS_TIME_TT: - sys = psStringCopy("TT"); - break; - default: - sys = psStringCopy("Unknown"); - } - psMetadataItem *newItem = psMetadataItemAllocStr("CELL.TIMESYS", "Time system", sys); - setConceptItem(fpa, chip, cell, db, newItem); - psFree(newItem); - psFree(sys); - } - - // CELL.TIME - { - psMetadataItem *timeItem = psMetadataLookup(cell->concepts, "CELL.TIME"); - psTime *time = timeItem->data.V; // The time - psString dateTimeString = psTimeToISO(time); // String representation - - psMetadata *formats = psMetadataLookupMD(NULL, fpa->camera, "FORMATS"); - psString format = psMetadataLookupString(NULL, formats, "CELL.TIME"); - - if (strlen(format) == 0) { - // No format specified --> do it in the usual way (maybe it's a DB lookup) - psMetadataItem *newTimeItem = psMetadataItemAllocStr("CELL.TIME", "Time of observation", - dateTimeString); - setConceptItem(fpa, chip, cell, db, newTimeItem); - psFree(newTimeItem); - } else { - if (strcasecmp(format, "ISO") == 0) { - // dateTimeString contains an ISO time - psMetadataItem *newTimeItem = psMetadataItemAllocStr("CELL.TIME", "Time of observation", - dateTimeString); - setConceptItem(fpa, chip, cell, db, newTimeItem); - psFree(newTimeItem); - } else if (strstr(format, "SEPARATE")) { - // We're working with two separate headers - psMetadata *header = NULL; // The FITS header - if (cell->hdu) { - header = cell->hdu->header; - } else if (chip->hdu) { - header = chip->hdu->header; - } else if (fpa->hdu) { - header = fpa->hdu->header; - } - if (! header) { - psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); - } else { - // Get the headers - const psMetadata *camera = fpa->camera; - psMetadata *translation = psMetadataLookupMD(NULL, camera, "TRANSLATION"); - psString keywords = psMetadataLookupString(NULL, translation, "CELL.TIME"); - psList *dateTimeKeywords = papSplit(keywords, " ,;"); - psString dateKeyword = psListGet(dateTimeKeywords, PS_LIST_HEAD); - psString timeKeyword = psListGet(dateTimeKeywords, PS_LIST_TAIL); - - psList *dateTime = papSplit(dateTimeString, " T"); // Find the middle T - psString dateString = psListGet(dateTime, PS_LIST_HEAD); - psString timeString = psListGet(dateTime, PS_LIST_TAIL); - - // XXX: Couldn't be bothered doing these right now - if (strstr(format, "PRE2000")) { - psError(PS_ERR_IO, true, "Don't you realise it's the twenty-first century?\n"); - } - if (strstr(format, "BACKWARDS")) { - psError(PS_ERR_IO, true, "You want it BACKWARDS? Not right now, thanks.\n"); - } - - psMetadataItem *dateItem = psMetadataItemAllocStr(dateKeyword, "Date of observation", - dateString); - psMetadataItem *timeItem = psMetadataItemAllocStr(timeKeyword, "Time of observation", - timeString); - psMetadataAddItem(header, dateItem, PS_LIST_TAIL, PS_META_REPLACE); - psMetadataAddItem(header, timeItem, PS_LIST_TAIL, PS_META_REPLACE); - - psFree(dateTimeKeywords); - psFree(dateTime); - psFree(dateItem); - psFree(timeItem); - } - } else if (strcasecmp(format, "MJD") == 0) { - double mjd = psTimeToMJD(time); - psMetadataItem *newTimeItem = psMetadataItemAllocF64("CELL.TIME", "MJD of observation", mjd); - setConceptItem(fpa, chip, cell, db, newTimeItem); - psFree(newTimeItem); - } else if (strcasecmp(format, "JD") == 0) { - double jd = psTimeToMJD(time); - psMetadataItem *newTimeItem = psMetadataItemAllocF64("CELL.TIME", "JD of observation", jd); - setConceptItem(fpa, chip, cell, db, newTimeItem); - psFree(newTimeItem); - } else { - psError(PS_ERR_IO, true, "Not sure how to outgest CELL.TIME (%s) --- trying " - "ISO\n", dateTimeString); - psMetadataItem *newTimeItem = psMetadataItemAllocStr("CELL.TIME", "Time of observation", - dateTimeString); - setConceptItem(fpa, chip, cell, db, newTimeItem); - psFree(newTimeItem); - } - } - - psFree(dateTimeString); - } - - // Remove corrective - { - const psMetadata *camera = fpa->camera; - bool mdok = false; // Result of MD lookup - psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); - if (mdok && formats) { - psString format = psMetadataLookupString(&mdok, formats, "CELL.X0"); - if (mdok && strlen(format) > 0 && strcasecmp(format, "FORTRAN") == 0) { - psMetadataItem *cellx0 = psMetadataLookup(cell->concepts, "CELL.X0"); - cellx0->data.S32 += 1; - } - format = psMetadataLookupString(&mdok, formats, "CELL.Y0"); - if (mdok && strlen(format) > 0 && strcasecmp(format, "FORTRAN") == 0) { - psMetadataItem *celly0 = psMetadataLookup(cell->concepts, "CELL.Y0"); - celly0->data.S32 += 1; - } - } - } - // CELL.X0 - setConcept(fpa, chip, cell, db, cell->concepts, "CELL.X0"); - // CELL.Y0 - setConcept(fpa, chip, cell, db, cell->concepts, "CELL.Y0"); - - // Pau. -} - Index: unk/archive/scripts/src/phase2/pmFPAConceptsSet.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAConceptsSet.h (revision 5798) +++ (revision ) @@ -1,17 +1,0 @@ -#ifndef PM_FPA_CONCEPTS_SET_H -#define PM_FPA_CONCEPTS_SET_H - -void pmFPAOutgestConcepts(pmFPA *fpa, // The FPA - psDB *db // DB handle - ); -void pmChipOutgestConcepts(pmChip *chip, // The chip - psDB *db // DB handle - ); -void pmCellOutgestConcepts(pmCell *cell, // The cell - psDB *db // DB handle - ); - - - - -#endif Index: unk/archive/scripts/src/phase2/pmFPAConstruct.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAConstruct.c (revision 5798) +++ (revision ) @@ -1,320 +1,0 @@ -#include -#include -#include "pslib.h" - -#include "pmFPA.h" -#include "pmFPAConstruct.h" - -#include "papStuff.h" - -// NOTE: Need to deal with header inheritance - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// File-static functions -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -// Read data for a particular cell from the camera configuration -static psMetadata *getCellData(const psMetadata *camera, // The camera configuration - const char *cellName // The name of the cell - ) -{ - bool status = true; // Result of MD lookup - psMetadata *cells = psMetadataLookupMD(&status, camera, "CELLS"); // The CELLS - if (! status) { - psError(PS_ERR_IO, false, "Unable to determine CELLS of camera.\n"); - return NULL; - } - - psMetadata *cellData = psMetadataLookupMD(&status, cells, cellName); // The data for the particular cell - if (! status) { - psError(PS_ERR_IO, false, "Unable to find specs for cell %s: ignored\n", cellName); - } - -#if 0 - // Need to create a new instance, so that each cell can work with its own - psMetadata *copy = psMetadataAlloc(); - psMetadataIterator *iter = psMetadataIteratorAlloc(cellData, PS_LIST_HEAD, NULL); // Iterator - psMetadataItem *item = NULL; // Item from iteration - while (item = psMetadataGetAndIncrement(iter)) { - if (item->type == PS_DATA_METADATA_MULTI || item->type == PS_DATA_METADATA) { - psLogMsg(__func__, PS_LOG_WARN, "PS_DATA_METADATA_MULTI and PS_DATA_METADATA are not supported " - "in a cell definition --- %s ignored.\n", item->name); - continue; - } - if (! psMetadataAdd(copy, PS_LIST_TAIL, item->name, item->type, item->comment, item->data.V)) { - psAbort(__func__, "Should never reach here!\n"); - } - } - psFree(iter); - - return copy; -#else - return cellData; -#endif - -} - - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Public functions -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -pmFPA *pmFPAConstruct(const psMetadata *camera // The camera configuration - ) -{ - pmFPA *fpa = pmFPAAlloc(camera); // The FPA to fill out - bool mdStatus = true; // Status from metadata lookups - const char *phuType = psMetadataLookupString(&mdStatus, camera, "PHU"); // What is the PHU? - const char *extType = psMetadataLookupString(&mdStatus, camera, "EXTENSIONS"); // What's in the extns? - - if (strcasecmp(phuType, "FPA") == 0) { - // The FITS file contains an entire FPA - - psMetadata *contents = psMetadataLookupMD(&mdStatus, camera, "CONTENTS"); // The CONTENTS - if (! mdStatus) { - psError(PS_ERR_IO, false, "Unable to determine CONTENTS of camera.\n"); - psFree(fpa); - return NULL; - } - - // Set up iteration over the contents - psMetadataIterator *contentsIter = psMetadataIteratorAlloc(contents, PS_LIST_HEAD, NULL); - psMetadataItem *contentItem = NULL; // Item from the metadata - - if (strcasecmp(extType, "CHIP") == 0) { - // Extensions are chips; Content contains a list of cells - while (contentItem = psMetadataGetAndIncrement(contentsIter)) { - psString extName = contentItem->name; // The name of the extension - pmChip *chip = pmChipAlloc(fpa, extName); // The chip - chip->hdu = p_pmHDUAlloc(extName); // Prepare chip to receive FITS data - if (contentItem->type != PS_DATA_STRING) { - psLogMsg(__func__, PS_LOG_WARN, "Type of content item (%x) is not string: ignored\n", - contentItem->type); - continue; - } - - const char *content = contentItem->data.V; // The content of the extension - psTrace(__func__, 7, "Content of %s is: %s\n", extName, content); - psList *cellNames = papSplit(content, " ,"); // A list of the component cells - psListIterator *cellNamesIter = psListIteratorAlloc(cellNames, PS_LIST_HEAD, NULL); - char *cellName = NULL; // The name of a cell - while (cellName = psListGetAndIncrement(cellNamesIter)) { - // Get the cell data - psMetadata *cellData = getCellData(camera, cellName); - pmCell *cell = pmCellAlloc(chip, cellData, cellName); // The cell -// psFree(cellData); - } - psFree(cellNamesIter); - psFree(cellNames); - } - - } else if (strcasecmp(extType, "CELL") == 0) { - // Extensions are cells; Content contains a chip name and cell type - psMetadata *chips = psMetadataAlloc(); // Given a chip name, holds the chip number - while (contentItem = psMetadataGetAndIncrement(contentsIter)) { - psString extName = contentItem->name; // The name of the extension - psTrace(__func__, 1, "Getting %s....\n", extName); - - if (contentItem->type != PS_DATA_STRING) { - psLogMsg(__func__, PS_LOG_WARN, "Type of content item (%x) is not string: ignored\n", - contentItem->type); - continue; - } - const char *content = contentItem->data.V; // The content of the extension - psList *contents = papSplit(content, ": "); // Split the name from the type - if (contents->n != 2) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to read contents of %s: ignored.\n", extName); - } else { - psString chipName = psListGet(contents, 0); // The name of the chip - psString cellType = psListGet(contents, 1); // The type of cell - psTrace(__func__, 7, "Extension is cell of type %s, from chip %s\n", cellType, - chipName); - - psMetadataItem *chipItem = psMetadataLookup(chips, chipName); // Item containing the chip - pmChip *chip = NULL; // The chip - if (! chipItem) { - chip = pmChipAlloc(fpa, chipName); - psMetadataAdd(chips, PS_LIST_TAIL, chipName, PS_DATA_UNKNOWN, "", chip); - } else { - chip = psMemIncrRefCounter(chipItem->data.V); - } - // The cell - psMetadata *cellData = getCellData(camera, cellType); - pmCell *cell = pmCellAlloc(chip, cellData, extName); // The cell -// psFree(cellData); - cell->hdu = p_pmHDUAlloc(extName); // Prepare cell to receive FITS data - - psFree(chip); - psFree(cell); - } - psFree(contents); - } - psFree(chips); - - } else if (strcasecmp(extType, "NONE") == 0) { - // No extensions; Content contains metadata, each entry is a chip with its component cells - fpa->hdu = p_pmHDUAlloc("PHU"); // Prepare FPA to receive FITS data - while (contentItem = psMetadataGetAndIncrement(contentsIter)) { - psString chipName = contentItem->name; // The name of the chip - - if (contentItem->type != PS_DATA_STRING) { - psLogMsg(__func__, PS_LOG_WARN, "Type of content item (%x) is not string: ignored\n", - contentItem->type); - continue; - } - psString content = contentItem->data.V; // The content of the extension - psTrace(__func__, 5, "Component cells are: %s\n", content); - pmChip *chip = pmChipAlloc(fpa, chipName); // The chip - psList *cellNames = papSplit(content, ", "); // Split the list of cells - psListIterator *cellNamesIter = psListIteratorAlloc(cellNames, PS_LIST_HEAD, false); - char *cellName = NULL; // Name of the cell - while (cellName = psListGetAndIncrement(cellNamesIter)) { - psTrace(__func__, 7, "Processing cell %s....\n", cellName); - psMetadata *cellData = getCellData(camera, cellName); - pmCell *cell = pmCellAlloc(chip, cellData, cellName); // The cell -// psFree(cellData); - } - psFree(cellNamesIter); - } - - } else { - psError(PS_ERR_IO, false, "EXTENSIONS in camera definition is not CHIP, CELL or NONE.\n"); - psFree(fpa); - return NULL; - } // Type of extension - - psFree(contentsIter); - - } else if (strcasecmp(phuType, "CHIP") == 0) { - // The FITS file contains a single chip only - psString chipName = psStringCopy("CHIP"); // Name for chip - pmChip *chip = pmChipAlloc(fpa, chipName); // The chip - - if (strcasecmp(extType, "NONE") == 0) { - // There are no extensions --- only the PHU - chip->hdu = p_pmHDUAlloc("PHU"); - - const char *contents = psMetadataLookupString(&mdStatus, camera, "CONTENTS"); - if (! mdStatus) { - psError(PS_ERR_IO, false, "Unable to determine CONTENTS of camera.\n"); - psFree(fpa); - return NULL; - } - psList *cellNames = papSplit(contents, " ,"); // Names of cells - psListIterator *cellIter = psListIteratorAlloc(cellNames, PS_LIST_HEAD, false); // Iterator - psString cellName = NULL; - while (cellName = psListGetAndIncrement(cellIter)) { - psMetadata *cellData = getCellData(camera, cellName); - pmCell *cell = pmCellAlloc(chip, cellData, cellName); // The cell -// psFree(cellData); - } - psFree(cellIter); - - } else if (strcasecmp(extType, "CELL") == 0) { - // Extensions are cells - psMetadata *contents = psMetadataLookupMD(&mdStatus, camera, "CONTENTS"); // The CONTENTS - if (! mdStatus) { - psError(PS_ERR_IO, false, "Unable to determine CONTENTS of camera.\n"); - psFree(fpa); - return NULL; - } - - if (strcasecmp(extType, "CELL") != 0) { - psLogMsg(__func__, PS_LOG_WARN, "EXTENSIONS in camera definition is %s, but PHU is CHIP.\n" - "EXTENSIONS assumed to be CELL.\n", extType); - } - - // Iterate through the contents - psMetadataIterator *contentsIter = psMetadataIteratorAlloc(contents, PS_LIST_HEAD, NULL); - psMetadataItem *contentItem = NULL; // Item from metadata - while (contentItem = psMetadataGetAndIncrement(contentsIter)) { - psString extName = contentItem->name; // The name of the extension - // Content is a cell type - if (contentItem->type != PS_DATA_STRING) { - psLogMsg(__func__, PS_LOG_WARN, - "CONTENT metadata for extension %s is not of type string, but %x --- ignored\n", - extName, contentItem->type); - continue; - } - const char *cellType = contentItem->data.V; // The type of cell - psTrace(__func__, 5, "Cell type is %s\n", cellType); - psMetadata *cellData = getCellData(camera, cellType); - pmCell *cell = pmCellAlloc(chip, cellData, extName); // The cell -// psFree(cellData); - cell->hdu = p_pmHDUAlloc(extName); // Prepare cell to receive FITS data - - psFree(cell); - } // Iterating through contents - psFree(contentsIter); - - } else { - psError(PS_ERR_IO, false, "EXTENSIONS in camera definition is neither CELL or NONE.\n"); - psFree(fpa); - return NULL; - } - psFree(chipName); - psFree(chip); - - } else { - psError(PS_ERR_IO, true, - "The PHU type specified in the camera configuration (%s) is not FPA or CHIP.\n", - phuType); - psFree(fpa); - return NULL; - } - - return fpa; -} - -// Print out the focal plane structure -void pmFPAPrint(pmFPA *fpa // FPA to print - ) -{ - psTrace(__func__, 1, "FPA:\n"); - if (fpa->hdu) { - psTrace(__func__, 2, "---> FPA is extension %s.\n", fpa->hdu->extname); - if (! fpa->hdu->images) { - psTrace(__func__, 2, "---> NO PIXELS for extension %s\n", fpa->hdu->extname); - } - } - psMetadataPrint(fpa->concepts, 2); - - psArray *chips = fpa->chips; // Array of chips - // Iterate over the FPA - for (int i = 0; i < chips->n; i++) { - psTrace(__func__, 3, "Chip: %d\n", i); - pmChip *chip = chips->data[i]; // The chip - if (chip->hdu) { - psTrace(__func__, 4, "---> Chip is extension %s.\n", chip->hdu->extname); - if (! chip->hdu->images) { - psTrace(__func__, 4, "---> NO PIXELS for extension %s\n", chip->hdu->extname); - } - } - psMetadataPrint(chip->concepts, 4); - - // Iterate over the chip - psArray *cells = chip->cells; // Array of cells - for (int j = 0; j < cells->n; j++) { - psTrace(__func__, 5, "Cell: %d\n", j); - pmCell *cell = cells->data[j]; // The cell - if (cell->hdu) { - psTrace(__func__, 6, "---> Cell is extension %s.\n", cell->hdu->extname); - if (! cell->hdu->images) { - psTrace(__func__, 6, "---> NO PIXELS for extension %s\n", cell->hdu->extname); - } - } - psMetadataPrint(cell->concepts, 6); - - psTrace(__func__, 7, "Readouts:\n"); - psArray *readouts = cell->readouts; // Array of readouts - for (int k = 0; k < readouts->n; k++) { - pmReadout *readout = readouts->data[k]; // The readout - psImage *image = readout->image; // The image - psTrace(__func__, 8, "Image: [%d:%d,%d:%d] (%dx%d)\n", image->col0, image->col0 + - image->numCols, image->row0, image->row0 + image->numRows, image->numCols, - image->numRows); - } // Iterating over cell - } // Iterating over chip - } // Iterating over FPA - -} Index: unk/archive/scripts/src/phase2/pmFPAConstruct.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAConstruct.h (revision 5798) +++ (revision ) @@ -1,16 +1,0 @@ -#ifndef PM_FPA_CONSTRUCT_H -#define PM_FPA_CONSTRUCT_H - -#include "pslib.h" -#include "pmFPA.h" - -// Read the contents of a FITS file (format specified by the camera configuration) into memory -pmFPA *pmFPAConstruct(const psMetadata *camera // The camera configuration - ); - -// Print out the FPA -void pmFPAPrint(pmFPA *fpa // FPA to print - ); - - -#endif Index: unk/archive/scripts/src/phase2/pmFPARead.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPARead.c (revision 5798) +++ (revision ) @@ -1,644 +1,0 @@ -#include -#include -#include "pslib.h" - -#include "pmFPA.h" -#include "pmFPAConceptsGet.h" -#include "pmFPARead.h" - -#include "papStuff.h" // For "split" - -// NOTE: Need to deal with header inheritance - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// File-static functions -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - - -// Read a FITS extension into a chip -static bool readExtension(p_pmHDU *hdu, // Pixel data into which to read - psFits *fits // The FITS file from which to read - ) -{ - const char *extName = hdu->extname; // Extension name - - psTrace(__func__, 7, "Moving to extension %s...\n", extName); - if (strncmp(extName, "PHU", 3) == 0) { - if (!psFitsMoveExtNum(fits, 0, false)) { - psError(PS_ERR_IO, false, "Unable to find PHU in FITS file!\n"); - return false; - } - } else if (! psFitsMoveExtName(fits, extName)) { - psError(PS_ERR_IO, false, "Unable to find extension %s in FITS file!\n", extName); - return false; - } - psTrace(__func__, 7, "Reading header....\n"); - psMetadata *header = psFitsReadHeader(NULL, fits); // Header - if (! header) { - psError(PS_ERR_IO, false, "Unable to read FITS header!\n"); - return false; - } - - psTrace(__func__, 7, "Checking NAXIS....\n"); - bool mdStatus = false; - int nAxis = psMetadataLookupS32(&mdStatus, header, "NAXIS"); - if (!mdStatus) { - psLogMsg(__func__, PS_LOG_WARN, "There is no NAXIS keyword in the FITS header of extension %s!\n", - extName); - } - if (nAxis != 2 && nAxis != 3) { - psLogMsg(__func__, PS_LOG_WARN, "Image is not 2- or 3-dimensional --- reading into a single image " - "anyway.\n"); - } - psTrace(__func__, 9, "NAXIS = %d\n", nAxis); - - int numPlanes = 1; // Number of planes - if (nAxis == 3) { - numPlanes = psMetadataLookupS32(&mdStatus, header, "NAXIS3"); - if (!mdStatus) { - psError(PS_ERR_IO, false, "Unable to read NAXIS3 for 3-dimensional image!\n"); - return false; - } - } - psRegion region = {0, 0, 0, 0}; // Region to read is everything - - // Read each plane into the array - psTrace(__func__, 7, "Reading %d planes into array....\n", numPlanes); - psArray *pixels = psArrayAlloc(numPlanes); // Array of images - for (int i = 0; i < numPlanes; i++) { - psTrace(__func__, 9, "Reading plane %d\n", i); - psMemCheckCorruption(true); - psImage *image = psFitsReadImage(NULL, fits, region, i); - - // XXX: Type conversion here to support the modules, which don't have multiple type support yet - if (image->type.type != PS_TYPE_F32) { - pixels->data[i] = psImageCopy(NULL, image, PS_TYPE_F32); - -#ifndef PRODUCTION - // XXX: Temporary fix for writing images, until psFits gets cleaned up - psMetadataItem *bitpixItem = psMetadataLookup(header, "BITPIX"); - bitpixItem->data.S32 = -32; - psMetadataRemove(header, 0, "BZERO"); - psMetadataRemove(header, 0, "BSCALE"); -#endif - - psFree(image); - } else { - pixels->data[i] = image; - } - psTrace(__func__, 10, "Done\n"); - if (! pixels->data[i]) { - psError(PS_ERR_IO, false, "Unable to read image for extension %s, plane %d\n", extName, i); - return false; - } - } - - // Update the HDU with the new data - hdu->images = pixels; - hdu->header = header; - - // XXX: Insert mask and weight reading here - - return true; -} - -// Portion out an image into the cell -static bool generateReadouts(pmCell *cell, // The cell that gets its bits - p_pmHDU *hdu // Pixel data, containing image, mask, weights - ) -{ - psArray *images = hdu->images; // Array of images (each of which is a readout) - psArray *masks = hdu->masks; // Array of masks (one for each readout) - // Iterate over each of the image planes - for (int i = 0; i < images->n; i++) { - psImage *image = images->data[i]; // The i-th plane - psImage *mask = NULL; // The mask - if (masks) { - mask = masks->data[i]; - } - - pmReadout *readout = pmReadoutAlloc(cell, image, mask, 0,0,1,1); - psFree(readout); // This seems silly, but the alloc registers it with the cell, so we - // don't have to do anything with it. Perhaps we should change - // pmReadoutAlloc to pmReadoutAdd? - } - - return true; -} - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Public functions -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -bool pmFPARead(pmFPA *fpa, // FPA to read into - psFits *fits, // FITS file from which to read - psMetadata *phu, // Primary header - psDB *db // Database handle, for concept ingest - ) -{ - p_pmHDU *hdu = NULL; // Pixel data from FITS file - - // Read the PHU, if required - if (! phu) { - if (! psFitsMoveExtNum(fits, 0, false)) { - psError(PS_ERR_IO, false, "Unable to find PHU in FITS file!\n"); - return false; - } - psTrace(__func__, 7, "Reading PHU....\n"); - fpa->phu = psFitsReadHeader(NULL, fits); // Primary header - } else { - fpa->phu = psMemIncrRefCounter(phu); - } - - // Read in.... - psTrace(__func__, 1, "Working on FPA...\n"); - if (fpa->hdu) { - hdu = fpa->hdu; - psTrace(__func__, 3, "Reading pixels from extension %s into FPA.\n", hdu->extname); - if (! readExtension(hdu, fits)) { - psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", hdu->extname); - return false; - } - } - pmFPAIngestConcepts(fpa, db); - - psArray *chips = fpa->chips; // Array of chips - // Iterate over the FPA - for (int i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // The chip - - // Only read chips marked "valid" - if (! chip->valid) { - psTrace(__func__, 2, "Ignoring chip %d...\n", i); - continue; - } - psTrace(__func__, 2, "Reading in chip %d...\n", i); - - if (chip->hdu) { - hdu = chip->hdu; - psTrace(__func__, 3, "Reading pixels from extension %s into chip %d.\n", hdu->extname, i); - if (! readExtension(hdu, fits)) { - psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", hdu->extname); - return false; - } - } - pmChipIngestConcepts(chip, db); - - // Iterate over the chip - psArray *cells = chip->cells; // Array of cells - for (int j = 0; j < cells->n; j++) { - pmCell *cell = cells->data[j]; // The cell - - // Only read cells marked "valid" - if (! cell->valid) { - psTrace(__func__, 3, "Ignoring chip %d...\n", i); - continue; - } - psTrace(__func__, 3, "Reading in cell %d...\n", j); - - if (cell->hdu) { - hdu = cell->hdu; - psTrace(__func__, 5, "Reading pixels from extension %s into cell %d.\n", hdu->extname, - j); - if (! readExtension(hdu, fits)) { - psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", - hdu->extname); - return false; - } - } - pmCellIngestConcepts(cell, db); - - psTrace(__func__, 5, "Allocating readouts for chip %d cell %d...\n", - i, j); - generateReadouts(cell, hdu); - } - } - - return true; -} - -// Translate a name from the configuration file, containing something like "%a_%d" to a real value -psString p_pmFPATranslateName(const psString name, // The name to translate - const pmCell *cell // The cell for which to translate - ) -{ - // %a is the FPA.NAME - // %b is the CHIP.NAME - // %c is the CELL.NAME - // %d is the chip number - // %e if the cell number - // %f is the extension name - - psString translation = psMemIncrRefCounter(name); // The translated string - char *temp = NULL; // Temporary string - - pmChip *chip = cell->parent; // Chip of interest - pmFPA *fpa = chip->parent; // FPA of interest - - // FPA.NAME - if (temp = strstr(translation, "%a")) { - // This is not particularly friendly to the memory allocation, but the cache should make it OK, - // and there's not much of it anyway... - psFree(translation); - translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string - psString fpaName = psMetadataLookupString(NULL, fpa->concepts, "FPA.NAME"); // The value of FPA.NAME - psStringAppend(&translation, "%s%s", fpaName, temp + 2); - // So "translation" now contains the first part, the replaced string, and the last part. - } - - // CHIP.NAME - if (temp = strstr(translation, "%b")) { - // This is not particularly friendly to the memory allocation, but the cache should make it OK, - // and there's not much of it anyway... - psFree(translation); - translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string - psString chipName = psMetadataLookupString(NULL, chip->concepts, "CHIP.NAME"); // CHIP.NAME's value - psStringAppend(&translation, "%s%s", chipName, temp + 2); - // So "translation" now contains the first part, the replaced string, and the last part. - } - - // CELL.NAME - if (temp = strstr(translation, "%c")) { - // This is not particularly friendly to the memory allocation, but the cache should make it OK, - // and there's not much of it anyway... - psFree(translation); - translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string - psString cellName = psMetadataLookupString(NULL, cell->concepts, "CELL.NAME"); // CELL.NAME's value - psStringAppend(&translation, "%s%s", cellName, temp + 2); - // So "translation" now contains the first part, the replaced string, and the last part. - } - - // Chip number - if (temp = strstr(translation, "%d")) { - // This is not particularly friendly to the memory allocation, but the cache should make it OK, - // and there's not much of it anyway... - psFree(translation); - translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string - // Search for the pointer to get the chip number - int chipNum = -1; - psArray *chips = fpa->chips; // The array of chips - for (int i = 0; i < chips->n && chipNum < 0; i++) { - if (chips->data[i] == chip) { - chipNum = i; - } - } - if (chipNum < 0) { - psError(PS_ERR_IO, true, "Unable to find chip to get name: %s\n", name); - // Try to muddle on by leaving the number out - psStringAppend(&translation, "%s", temp + 2); - } else { - psStringAppend(&translation, "%d%s", chipNum, temp + 2); - } - // So "translation" now contains the first part, the replaced string, and the last part. - } - - // Cell number - if (temp = strstr(translation, "%e")) { - // This is not particularly friendly to the memory allocation, but the cache should make it OK, - // and there's not much of it anyway... - psFree(translation); - translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string - // Search for the pointer to get the cell number - int cellNum = -1; - psArray *cells = chip->cells; // The array of cells - for (int i = 0; i < cells->n && cellNum < 0; i++) { - if (cells->data[i] == cell) { - cellNum = i; - } - } - if (cellNum < 0) { - psError(PS_ERR_IO, true, "Unable to find cell to get name: %s\n", name); - // Try to muddle on by leaving the number out - psStringAppend(&translation, "%s", temp + 2); - } else { - psStringAppend(&translation, "%d%s", cellNum, temp + 2); - } - // So "translation" now contains the first part, the replaced string, and the last part. - } - - // Extension name - if (temp = strstr(translation, "%f")) { - // This is not particularly friendly to the memory allocation, but the cache should make it OK, - // and there's not much of it anyway... - psFree(translation); - translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string - const char *extname = NULL; - if (cell->hdu) { - extname = cell->hdu->extname; - } else if (chip->hdu) { - extname = chip->hdu->extname; - } else if (fpa->hdu) { - extname = fpa->hdu->extname; - } - psStringAppend(&translation, "%s%s", extname, temp + 2); - // So "translation" now contains the first part, the replaced string, and the last part. - } - - return translation; -} - -// Get filename and extension out of "somefile.fits:ext" -psString p_pmFPATranslateFileExt(psString *extName, // Extension name, to be returned - const psString name, // The string to be translated and then parsed - const pmCell *cell // The cell - ) -{ - psString filenameExt = p_pmFPATranslateName(name, cell); - const char *colon = strchr(filenameExt, ':'); // Pointer to a colon in the filename-extn - psString filename = NULL; // The filename - if (extName && *extName) { - psFree(*extName); - } - if (extName) { - *extName = NULL; - } - - if (colon) { - filename = psStringNCopy(filenameExt, strlen(filenameExt) - strlen(colon)); - if (strlen(colon) > 1) { - if (extName) { - *extName = psStringCopy(colon + 1); - } - } - } else { - filename = psMemIncrRefCounter(filenameExt); - } - - return filename; -} - -// Read a mask into the FPA -bool pmFPAReadMask(pmFPA *fpa, // FPA to read into - psFits *source // Source FITS file (for the original data) - ) -{ - const psMetadata *camera = fpa->camera; // Camera configuration for FPA - bool mdok = false; // Status of MD lookup - - // Get the required information from the camera configuration - psMetadata *supps = psMetadataLookupMD(&mdok, camera, "SUPPLEMENTARY"); // Rules for supplementary data - if (! mdok || ! supps) { - psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n"); - return false; - } - psString sourceType = psMetadataLookupString(&mdok, supps, "MASK.SOURCE"); // Type of source: EXT | FILE - if (! mdok || strlen(sourceType) <= 0) { - psError(PS_ERR_IO, false, "Unable to find MASK.SOURCE in SUPPLEMENTARY section of camera " - "configuration!\n"); - return false; - } - psString name = psMetadataLookupString(&mdok, supps, "MASK.NAME"); // Name of mask - if (! mdok || strlen(sourceType) <= 0) { - psError(PS_ERR_IO, false, "Unable to find MASK.NAME in SUPPLEMENTARY section of camera " - "configuration!\n"); - return false; - } - - // Go through the FPA to each cell/readout to get the mask - p_pmHDU *hdu = fpa->hdu; // The HDU into which we will read the mask - psArray *chips = fpa->chips; // Array of chips - for (int chipNum = 0; chipNum < chips->n; chipNum++) { - pmChip *chip = chips->data[chipNum]; // The current chip of interest - if (chip->valid) { - if (chip->hdu) { - hdu = chip->hdu; - } - psArray *cells = chip->cells; // Array of cells - for (int cellNum = 0; cellNum < cells->n; cellNum++) { - pmCell *cell = cells->data[cellNum]; // The current cell of interest - if (cell->valid) { - if (cell->hdu) { - hdu = cell->hdu; - } - - // Now, need to find out where to get the pixels - psFits *maskSource = psMemIncrRefCounter(source); // Source of mask image - if (strcasecmp(sourceType, "FILE") == 0) { - // Source is a file (with optional extension, e.g., "myMaskFile.fits:thisExt" - psString extname = NULL; // Extension name - psString filename = p_pmFPATranslateFileExt(&extname, name, cell); // The filename - - psFree(maskSource); - maskSource = psFitsOpen(filename, "r"); - if (extname) { - if (! psFitsMoveExtName(maskSource, extname)) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to find extension %s to read mask.\n", - extname); - return false; - } - } - psFree(filename); - psFree(extname); - } else if (strncasecmp(sourceType, "EXT", 3) == 0) { - // Source is an extension in the original file - psString extname = p_pmFPATranslateName(name, cell); - psFitsMoveExtName(maskSource, extname); - } - - // We've arrived where the pixels are. Now we need to read them in. - psMetadata *header = psFitsReadHeader(NULL, maskSource); // The header - bool mdStatus = false; - int nAxis = psMetadataLookupS32(&mdStatus, header, "NAXIS"); - if (!mdStatus) { - psLogMsg(__func__, PS_LOG_WARN, "There is no NAXIS keyword in the FITS header for " - "mask (%s)!\n", name); - } - if (nAxis != 2 && nAxis != 3) { - psLogMsg(__func__, PS_LOG_WARN, "Image is not 2- or 3-dimensional --- reading into " - "a single image anyway.\n"); - } - - int numPlanes = 1; // Number of planes - if (nAxis == 3) { - numPlanes = psMetadataLookupS32(&mdStatus, header, "NAXIS3"); - if (!mdStatus) { - psError(PS_ERR_IO, false, "Unable to read NAXIS3 for 3-dimensional image!\n"); - // Try to proceed by taking only the first plane - numPlanes = 1; - } - if (numPlanes != 1 && numPlanes != cell->readouts->n) { - psError(PS_ERR_IO, false, "Number of masks (%d) does not match number of " - "readouts (%d)\n", numPlanes, cell->readouts->n); - // Try to proceed by taking only the first plane - numPlanes = 1; - } - } - - hdu->masks = psArrayAlloc(hdu->images->n); - - // Read each plane into the array - psArray *readouts = cell->readouts; // The array of readouts - for (int i = 0; i < hdu->masks->n; i++) { - psImage *mask = NULL; // The mask to be added - if (i < numPlanes) { - // Read the mask from the file - psTrace(__func__, 9, "Reading plane %d\n", i); - psRegion region = {0, 0, 0, 0}; - mask = psFitsReadImage(NULL, maskSource, region, i); - } else { - // One mask in the file is provided for all planes in the original image - psTrace(__func__, 9, "Copying plane 0 into plane %d\n", i); - psImage *original = hdu->masks->data[0]; - mask = psImageCopy(NULL, original, original->type.type); - } - hdu->masks->data[0] = mask; - pmReadout *readout = readouts->data[i]; - readout->mask = mask; - // Check the dimensions - // XXX: Perhaps this check should be against the CELL.TRIMSEC, not the image size? - if (mask->numCols < readout->image->numCols || - mask->numRows < readout->image->numRows) - { - psError(PS_ERR_IO, false, "Mask size (%dx%d) not compatible with image (%dx%d)\n", - mask->numCols, mask->numRows, readout->image->numCols, - readout->image->numRows); - return false; - } - } // Iterating over readouts - psFree(maskSource); - } // Valid cells - } // Iterating over cells - } // Valid chips - } // Iterating over chips - - return true; -} - - -// Read a mask into the FPA -// This is just a copy of the above pmFPAReadMask, replacing "mask" with "weight" throughout. -bool pmFPAReadWeight(pmFPA *fpa, // FPA to read into - psFits *source // Source FITS file (for the original data) - ) -{ - const psMetadata *camera = fpa->camera; // Camera configuration for FPA - bool mdok = false; // Status of MD lookup - - // Get the required information from the camera configuration - psMetadata *supps = psMetadataLookupMD(&mdok, camera, "SUPPLEMENTARY"); // Rules for supplementary data - if (! mdok || ! supps) { - psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n"); - return false; - } - psString sourceType = psMetadataLookupString(&mdok, supps, "WEIGHT.SOURCE"); // Type of source: EXT | FILE - if (! mdok || strlen(sourceType) <= 0) { - psError(PS_ERR_IO, false, "Unable to find WEIGHT.SOURCE in SUPPLEMENTARY section of camera " - "configuration!\n"); - return false; - } - psString name = psMetadataLookupString(&mdok, supps, "WEIGHT.NAME"); // Name of weight - if (! mdok || strlen(sourceType) <= 0) { - psError(PS_ERR_IO, false, "Unable to find WEIGHT.NAME in SUPPLEMENTARY section of camera " - "configuration!\n"); - return false; - } - - // Go through the FPA to each cell/readout to get the weight - p_pmHDU *hdu = fpa->hdu; // The HDU into which we will read the weight - psArray *chips = fpa->chips; // Array of chips - for (int chipNum = 0; chipNum < chips->n; chipNum++) { - pmChip *chip = chips->data[chipNum]; // The current chip of interest - if (chip->valid) { - if (chip->hdu) { - hdu = chip->hdu; - } - psArray *cells = chip->cells; // Array of cells - for (int cellNum = 0; cellNum < cells->n; cellNum++) { - pmCell *cell = cells->data[cellNum]; // The current cell of interest - if (cell->valid) { - if (cell->hdu) { - hdu = cell->hdu; - } - - // Now, need to find out where to get the pixels - psFits *weightSource = psMemIncrRefCounter(source); // Source of weight image - if (strcasecmp(sourceType, "FILE") == 0) { - // Source is a file (with optional extension, e.g., "myWeightFile.fits:thisExt" - psString extname = NULL; // Extension name - psString filename = p_pmFPATranslateFileExt(&extname, name, cell); // The filename - - psFree(weightSource); - weightSource = psFitsOpen(filename, "r"); - if (extname) { - if (! psFitsMoveExtName(weightSource, extname)) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to find extension %s to read " - "weight.\n", extname); - return false; - } - } - psFree(filename); - psFree(extname); - } else if (strncasecmp(sourceType, "EXT", 3) == 0) { - // Source is an extension in the original file - psString extname = p_pmFPATranslateName(name, cell); - psFitsMoveExtName(weightSource, extname); - } - - // We've arrived where the pixels are. Now we need to read them in. - psMetadata *header = psFitsReadHeader(NULL, weightSource); // The header - bool mdStatus = false; - int nAxis = psMetadataLookupS32(&mdStatus, header, "NAXIS"); - if (!mdStatus) { - psLogMsg(__func__, PS_LOG_WARN, "There is no NAXIS keyword in the FITS header for " - "weight (%s)!\n", name); - } - if (nAxis != 2 && nAxis != 3) { - psLogMsg(__func__, PS_LOG_WARN, "Image is not 2- or 3-dimensional --- reading into " - "a single image anyway.\n"); - } - - int numPlanes = 1; // Number of planes - if (nAxis == 3) { - numPlanes = psMetadataLookupS32(&mdStatus, header, "NAXIS3"); - if (!mdStatus) { - psError(PS_ERR_IO, false, "Unable to read NAXIS3 for 3-dimensional image!\n"); - // Try to proceed by taking only the first plane - numPlanes = 1; - } - if (numPlanes != 1 && numPlanes != cell->readouts->n) { - psError(PS_ERR_IO, false, "Number of weights (%d) does not match number of " - "readouts (%d)\n", numPlanes, cell->readouts->n); - // Try to proceed by taking only the first plane - numPlanes = 1; - } - } - - hdu->weights = psArrayAlloc(hdu->images->n); - - // Read each plane into the array - psArray *readouts = cell->readouts; // The array of readouts - for (int i = 0; i < hdu->weights->n; i++) { - psImage *weight = NULL; // The weight to be added - if (i < numPlanes) { - // Read the weight from the file - psTrace(__func__, 9, "Reading plane %d\n", i); - psRegion region = {0, 0, 0, 0}; - weight = psFitsReadImage(NULL, weightSource, region, i); - } else { - // One weight in the file is provided for all planes in the original image - psTrace(__func__, 9, "Copying plane 0 into plane %d\n", i); - psImage *original = hdu->weights->data[0]; - weight = psImageCopy(NULL, original, original->type.type); - } - hdu->weights->data[0] = weight; - pmReadout *readout = readouts->data[i]; - readout->weight = weight; - // Check the dimensions - // XXX: Perhaps this check should be against the CELL.TRIMSEC, not the image size? - if (weight->numCols < readout->image->numCols || - weight->numRows < readout->image->numRows) - { - psError(PS_ERR_IO, false, "Weight size (%dx%d) not compatible with image " - "(%dx%d)\n", weight->numCols, weight->numRows, readout->image->numCols, - readout->image->numRows); - return false; - } - } // Iterating over readouts - psFree(weightSource); - } // Valid cells - } // Iterating over cells - } // Valid chips - } // Iterating over chips - - return true; -} Index: unk/archive/scripts/src/phase2/pmFPARead.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPARead.h (revision 5798) +++ (revision ) @@ -1,30 +1,0 @@ -#ifndef PM_FPA_READ_H -#define PM_FPA_READ_H - -#include "pmFPA.h" - -bool pmFPARead(pmFPA *fpa, // FPA to read into - psFits *fits, // FITS file from which to read - psMetadata *phu, // Primary header - psDB *db // Database handle, for concept ingest - ); - -psString p_pmFPATranslateName(const psString name, // The name to translate - const pmCell *cell // The cell for which to translate - ); - -psString p_pmFPATranslateFileExt(psString *extName, // Extension name, to be returned - const psString filenameExt, // The string to parse into filename and ext - const pmCell *cell // The cell - ); - -bool pmFPAReadMask(pmFPA *fpa, // FPA to read into - psFits *source // Source FITS file (for the original data) - ); - -bool pmFPAReadWeight(pmFPA *fpa, // FPA to read into - psFits *source // Source FITS file (for the original data) - ); - - -#endif Index: unk/archive/scripts/src/phase2/pmFPAWrite.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAWrite.c (revision 5798) +++ (revision ) @@ -1,257 +1,0 @@ -#include -#include -#include "pslib.h" - -#include "pmFPA.h" -#include "pmFPAConceptsSet.h" -#include "pmFPARead.h" - -static bool writeHDU(psFits *fits, // FITS file to which to write - p_pmHDU *hdu // Pixel data to write - ) -{ - bool status = true; // Status of write, to return - for (int i = 0; i < hdu->images->n; i++) { - status &= psFitsWriteImage(fits, hdu->header, hdu->images->data[i], i); - // XXX: Insert here the writing on mask and weight images - } - - return status; -} - - -bool pmFPAWrite(psFits *fits, // FITS file to which to write - pmFPA *fpa, // FPA to write - psDB *db // Database to update - ) -{ - bool status = true; // Status of writing, to return - - psTrace(__func__, 7, "Outgesting FPA concepts...\n"); - pmFPAOutgestConcepts(fpa, db); - - // Write the primary header - if (fpa->phu) { - status &= psFitsMoveExtNum(fits, 0, false); - status &= psFitsWriteHeader(fpa->phu, fits); - } - - psArray *chips = fpa->chips; // Array of component chips - for (int i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // The component chip - if (chip->valid) { - psTrace(__func__, 1, "Writing out chip %d...\n", i); - - pmChipOutgestConcepts(chip, db); - - psArray *cells = chip->cells; // Array of component cells - for (int j = 0; j < cells->n; j++) { - pmCell *cell = cells->data[j]; // The component cell - if (cell->valid) { - psTrace(__func__, 2, "Writing out cell, %d...\n", j); - - pmCellOutgestConcepts(cell, db); - - if (cell->hdu && strlen(cell->hdu->extname) > 0) { - status &= writeHDU(fits, cell->hdu); - } - } - } - - if (chip->hdu && strlen(chip->hdu->extname) > 0) { - status &= writeHDU(fits, chip->hdu); - } - } - - } - - if (fpa->hdu && strlen(fpa->hdu->extname) > 0) { - status &= writeHDU(fits, fpa->hdu); - } - - - - return status; -} - - -bool pmFPAWriteMask(pmFPA *fpa, // FPA containing mask to write - psFits *fits // FITS file for image - ) -{ - const psMetadata *camera = fpa->camera; // Camera configuration for FPA - bool mdok = false; // Status of MD lookup - - // Get the required information from the camera configuration - psMetadata *supps = psMetadataLookupMD(&mdok, camera, "SUPPLEMENTARY"); // Rules for supplementary data - if (! mdok || ! supps) { - psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n"); - return false; - } - psString sourceType = psMetadataLookupString(&mdok, supps, "MASK.SOURCE"); // Type of source: EXT | FILE - if (! mdok || strlen(sourceType) <= 0) { - psError(PS_ERR_IO, false, "Unable to find MASK.SOURCE in SUPPLEMENTARY section of camera " - "configuration!\n"); - return false; - } - psString name = psMetadataLookupString(&mdok, supps, "MASK.NAME"); // Name of mask - if (! mdok || strlen(sourceType) <= 0) { - psError(PS_ERR_IO, false, "Unable to find MASK.NAME in SUPPLEMENTARY section of camera " - "configuration!\n"); - return false; - } - - // Go through the FPA to each cell/readout to get the mask - p_pmHDU *hdu = fpa->hdu; // The HDU into which we will read the mask - psArray *chips = fpa->chips; // Array of chips - for (int chipNum = 0; chipNum < chips->n; chipNum++) { - pmChip *chip = chips->data[chipNum]; // The current chip of interest - if (chip->valid) { - if (chip->hdu) { - hdu = chip->hdu; - } - psArray *cells = chip->cells; // Array of cells - for (int cellNum = 0; cellNum < cells->n; cellNum++) { - pmCell *cell = cells->data[cellNum]; // The current cell of interest - if (cell->valid) { - if (cell->hdu) { - hdu = cell->hdu; - } - - // Now, need to find out where to write the pixels - psFits *maskDest = psMemIncrRefCounter(fits); // Destination of mask image - psMetadata *header = psMetadataAlloc(); // A dummy header, containing the extension name - if (strcasecmp(sourceType, "FILE") == 0) { - // Source is a file (with optional extension, e.g., "myMaskFile.fits:thisExt" - psString extname = NULL; // Extension name - psString filename = p_pmFPATranslateFileExt(&extname, name, cell); // The filename - psFree(maskDest); - maskDest = psFitsOpen(filename, "rw"); - if (extname) { - psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); - } - psFree(filename); - psFree(extname); - } else if (strncasecmp(sourceType, "EXT", 3) == 0) { - // Source is an extension in the original file - psString extname = p_pmFPATranslateName(name, cell); - psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); - psFree(extname); - } - - // We've arrived where the pixels are. Now we need to write them out. - psArray *readouts = cell->readouts; // The array of readouts - for (int readNum = 0; readNum < readouts->n; readNum++) { - pmReadout *readout = readouts->data[readNum]; // The readout of interest - if (! readout->mask) { - psLogMsg(__func__, PS_LOG_WARN, "No mask to write out in %d,%d,%d\n", - chipNum, cellNum, readNum); - } else { - // XXX: Need to add the extname to the existing header - if (! psFitsWriteImage(maskDest, header, readout->mask, readNum)) { - psError(PS_ERR_IO, false, "Unable to write mask plane %d in extension %s\n", - readNum, hdu->extname); - return false; - } - } - } // Iterating over readouts - psFree(header); - psFree(maskDest); - } // Valid cells - } // Iterating over cells - } // Valid chips - } // Iterating over chips - - return true; -} - - -bool pmFPAWriteWeight(pmFPA *fpa, // FPA containing mask to write - psFits *fits // FITS file for image - ) -{ - const psMetadata *camera = fpa->camera; // Camera configuration for FPA - bool mdok = false; // Status of MD lookup - - // Get the required information from the camera configuration - psMetadata *supps = psMetadataLookupMD(&mdok, camera, "SUPPLEMENTARY"); // Rules for supplementary data - if (! mdok || ! supps) { - psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n"); - return false; - } - psString sourceType = psMetadataLookupString(&mdok, supps, "WEIGHT.SOURCE"); // Type of source: EXT | FILE - if (! mdok || strlen(sourceType) <= 0) { - psError(PS_ERR_IO, false, "Unable to find WEIGHT.SOURCE in SUPPLEMENTARY section of camera " - "configuration!\n"); - return false; - } - psString name = psMetadataLookupString(&mdok, supps, "WEIGHT.NAME"); // Name of weight - if (! mdok || strlen(sourceType) <= 0) { - psError(PS_ERR_IO, false, "Unable to find WEIGHT.NAME in SUPPLEMENTARY section of camera " - "configuration!\n"); - return false; - } - - // Go through the FPA to each cell/readout to get the weight - p_pmHDU *hdu = fpa->hdu; // The HDU into which we will read the weight - psArray *chips = fpa->chips; // Array of chips - for (int chipNum = 0; chipNum < chips->n; chipNum++) { - pmChip *chip = chips->data[chipNum]; // The current chip of interest - if (chip->valid) { - if (chip->hdu) { - hdu = chip->hdu; - } - psArray *cells = chip->cells; // Array of cells - for (int cellNum = 0; cellNum < cells->n; cellNum++) { - pmCell *cell = cells->data[cellNum]; // The current cell of interest - if (cell->valid) { - if (cell->hdu) { - hdu = cell->hdu; - } - - // Now, need to find out where to write the pixels - psFits *weightDest = psMemIncrRefCounter(fits); // Destination of weight image - psMetadata *header = psMetadataAlloc(); // A dummy header, containing the extension name - if (strcasecmp(sourceType, "FILE") == 0) { - // Source is a file (with optional extension, e.g., "myWeightFile.fits:thisExt" - psString extname = NULL; // Extension name - psString filename = p_pmFPATranslateFileExt(&extname, name, cell); // The filename - - psFree(weightDest); - weightDest = psFitsOpen(filename, "rw"); - if (extname) { - psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); - } - psFree(filename); - psFree(extname); - } else if (strncasecmp(sourceType, "EXT", 3) == 0) { - // Source is an extension in the original file - psString extname = p_pmFPATranslateName(name, cell); - psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); - psFree(extname); - } - - // We've arrived where the pixels are. Now we need to write them out. - psArray *readouts = cell->readouts; // The array of readouts - for (int readNum = 0; readNum < readouts->n; readNum++) { - pmReadout *readout = readouts->data[readNum]; // The readout of interest - if (! readout->weight) { - psLogMsg(__func__, PS_LOG_WARN, "No weight image to write out in %d,%d,%d\n", - chipNum, cellNum, readNum); - } else { - if (! psFitsWriteImage(weightDest, header, readout->weight, readNum)) { - psError(PS_ERR_IO, false, "Unable to write weight plane %d in extension %s\n", - readNum, hdu->extname); - return false; - } - } - } // Iterating over readouts - psFree(header); - psFree(weightDest); - } // Valid cells - } // Iterating over cells - } // Valid chips - } // Iterating over chips - - return true; -} Index: unk/archive/scripts/src/phase2/pmFPAWrite.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAWrite.h (revision 5798) +++ (revision ) @@ -1,20 +1,0 @@ -#ifndef PM_FPA_WRITE_H -#define PM_FPA_WRITE_H - -#include "pmFPA.h" - -bool pmFPAWrite(psFits *fits, // FITS file to which to write - pmFPA *fpa, // FPA to write - psDB *db // Database to update - ); - -bool pmFPAWriteMask(pmFPA *fpa, // FPA containing mask to write - psFits *fits // FITS file for image - ); - -bool pmFPAWriteWeight(pmFPA *fpa, // FPA containing mask to write - psFits *fits // FITS file for image - ); - - -#endif Index: unk/archive/scripts/src/phase2/pmFlatField.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFlatField.c (revision 5798) +++ (revision ) @@ -1,182 +1,0 @@ -/** @file pmFlatField.c - * - * @brief Given an input image and a flat field image, pmFlatField shall divide the input image by the flat - * field image. - * - * The input image, in, and the flat field image, flat, need not be the same size, since the input image may - * already have been trimmed (following overscan subtraction), but the function shall use the offsets in the - * image (in->x0 and in->y0) to determine the appropriate offsets to obtain the correct pixel on the flat - * field. In the event that the flat image is too small (i.e., pixels on the input image refer to pixels - * outside the range of the flat image), the function shall generate an error. Pixels which are negative or - * zero in the flat shall be masked in the input image with the value PM_MASK_FLAT. Negative pixels in the - * flat may be set to zero so that they are treated identically to zeroes. Any pixels masked in the flat - * shall be masked with corresponding values in the output. The function shall not normalize the flat; this - * responsibility is left to the caller. This function is basically equivalent to a divide (with psImageOp), - * but with care for the region that is divided, checking for negative pixels, and copying of the mask from - * the flat to the output. - * - * @author Ross Harman, MHPCC - * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-30 21:46:15 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - */ - -#if HAVE_CONFIG_H -#include -#endif - -#include -#include -#include - -#include "pslib.h" -#include "pmFlatField.h" -#include "pmMaskBadPixels.h" -#include "pmFlatFieldErrors.h" - - -bool pmFlatField(pmReadout *in, pmReadout *mask, const pmReadout *flat) -{ - // XXX: Not sure if this is correct. Must consult with IfA. - PS_ASSERT_PTR_NON_NULL(mask, false); - int i = 0; - int j = 0; - int totOffCol = 0; - int totOffRow = 0; - psElemType inType; - psElemType flatType; - psElemType maskType; - psImage *inImage = NULL; - psImage *inMask = NULL; - psImage *flatImage = NULL; - - - // Check for nulls - if (in == NULL) { - return true; // Readout may not have data in it - } else if(flat==NULL) { - psError( PS_ERR_BAD_PARAMETER_NULL, true, - PS_ERRORTEXT_pmFlatField_NULL_FLAT_READOUT); - return false; - } - - inImage = in->image; - flatImage = flat->image; - if (inImage == NULL) { - psError( PS_ERR_BAD_PARAMETER_NULL, true, - PS_ERRORTEXT_pmFlatField_NULL_INPUT_IMAGE); - return false; - } else if(flatImage == NULL) { - psError( PS_ERR_BAD_PARAMETER_NULL, true, - PS_ERRORTEXT_pmFlatField_NULL_FLAT_IMAGE); - return false; - } - inMask = mask->image; - - // Check input image and its mask are not larger than flat image - - if (inImage->numRows>flatImage->numRows || inImage->numCols>flatImage->numCols) { - psError( PS_ERR_BAD_PARAMETER_SIZE, true, - PS_ERRORTEXT_pmFlatField_SIZE_INPUT_IMAGE, - inImage->numRows, inImage->numCols, flatImage->numRows, flatImage->numCols); - return false; - } - if (inMask->numRows > flatImage->numRows || inMask->numCols > flatImage->numCols) { - psError( PS_ERR_BAD_PARAMETER_SIZE, true, - PS_ERRORTEXT_pmFlatField_SIZE_MASK_IMAGE, - inMask->numRows, inMask->numCols, flatImage->numRows, flatImage->numCols); - return false; - } - - // Determine total offset based on image offset with chip offset - totOffCol = inImage->col0 + in->col0; - totOffRow = inImage->row0 + in->row0; - - // Check that offsets are within image limits - if(totOffRow>=flatImage->numRows || totOffCol>=flatImage->numCols) { - psError( PS_ERR_BAD_PARAMETER_SIZE, true, - PS_ERRORTEXT_pmFlatField_OFFSET_FLAT_IMAGE, - totOffRow, totOffCol, flatImage->numRows, flatImage->numCols); - return false; - } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) { - psError( PS_ERR_BAD_PARAMETER_SIZE, true, - PS_ERRORTEXT_pmFlatField_OFFSET_INPUT_IMAGE, - totOffRow, totOffCol, inImage->numRows, inImage->numCols); - return false; - } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) { - psError( PS_ERR_BAD_PARAMETER_SIZE, true, - PS_ERRORTEXT_pmFlatField_OFFSET_MASK_IMAGE, - totOffRow, totOffCol, inMask->numRows, inMask->numCols); - return false; - } - - // Check for incorrect types - inType = inImage->type.type; - flatType = flatImage->type.type; - maskType = inMask->type.type; - if(PS_IS_PSELEMTYPE_COMPLEX(inType)) { - psError( PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_pmFlatField_TYPE_INPUT_IMAGE, - inType); - return false; - } else if(PS_IS_PSELEMTYPE_COMPLEX(flatType)) { - psError( PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_pmFlatField_TYPE_FLAT_IMAGE, - flatType); - return false; - } else if(maskType != PS_TYPE_MASK) { - psError( PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_pmFlatField_TYPE_MASK_IMAGE, - maskType); - return false; - } else if(inType != flatType) { - psError( PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_pmFlatField_TYPE_MISMATCH, - inType, flatType); - return false; - } - - // Macro for all PS types - #define PM_FLAT_DIVISION(TYPE) \ -case PS_TYPE_##TYPE: \ - /* Per Eugene's request, use two sets of loops: first to fill mask, second to avoid div with bad pix */ \ - for(j = totOffRow; j < inImage->numRows; j++) { \ - for(i = totOffCol; i < inImage->numCols; i++) { \ - if(flatImage->data.TYPE[j][i] <= 0.0) { \ - /* Negative or zero flat pixels shall be masked in input image as PM_MASK_FLAT */ \ - inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_FLAT; \ - flatImage->data.TYPE[j][i] = 0.0; \ - } \ - } \ - } \ - for(j = totOffRow; j < inImage->numRows; j++) { \ - for(i = totOffCol; i < inImage->numCols; i++) { \ - if(!inMask->data.PS_TYPE_MASK_DATA[j][i]) { \ - /* Module shall divide the input image by the flat-fielded image */ \ - inImage->data.TYPE[j][i] /= flatImage->data.TYPE[j][i]; \ - } \ - } \ - } \ - break; - - switch(inType) { - PM_FLAT_DIVISION(U8); - PM_FLAT_DIVISION(U16); - PM_FLAT_DIVISION(U32); - PM_FLAT_DIVISION(U64); - PM_FLAT_DIVISION(S8); - PM_FLAT_DIVISION(S16); - PM_FLAT_DIVISION(S32); - PM_FLAT_DIVISION(S64); - PM_FLAT_DIVISION(F32); - PM_FLAT_DIVISION(F64); - default: - psError( PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_pmFlatField_TYPE_UNSUPPORTED, - inType); - } - - return true; -} Index: unk/archive/scripts/src/phase2/pmFlatField.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmFlatField.h (revision 5798) +++ (revision ) @@ -1,42 +1,0 @@ -/** @file pmFlatField.h - * - * @brief Given an input image and a flat field image, pmFlatField shall divide the input image by the flat - * field image. - * - * The input image, in, and the flat field image, flat, need not be the same size, since the input image may - * already have been trimmed (following overscan subtraction), but the function shall use the offsets in the - * image (in->x0 and in->y0) to determine the appropriate offsets to obtain the correct pixel on the flat - * field. In the event that the flat image is too small (i.e., pixels on the input image refer to pixels - * outside the range of the flat image), the function shall generate an error. Pixels which are negative or - * zero in the flat shall be masked in the input image with the value PM_MASK_FLAT. Negative pixels in the - * flat may be set to zero so that they are treated identically to zeroes. Any pixels masked in the flat - * shall be masked with corresponding values in the output. The function shall not normalize the flat; this - * responsibility is left to the caller. This function is basically equivalent to a divide (with psImageOp), - * but with care for the region that is divided, checking for negative pixels, and copying of the mask from - * the flat to the output. - * - * @author Ross Harman, MHPCC - * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-03 01:30:32 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - */ - -#include "pslib.h" -#include "pmFPA.h" // #include "pmAstrometry.h" - - -/** Execute flat field module. - * - * Given an input image and a flat-field image, pmFlatField shall divide the input image by the flat field - * image. - * - * @return bool: True or false for success or failure - */ -bool pmFlatField( - pmReadout *in, ///< Readout with input image - pmReadout *mask, ///< Input image mask - const pmReadout *flat ///< Readout with flat image -); - Index: unk/archive/scripts/src/phase2/pmFlatFieldErrors.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmFlatFieldErrors.h (revision 5798) +++ (revision ) @@ -1,47 +1,0 @@ -/** @file pmFlatFieldErrors.h - * - * @brief Contains the error text for the flat field module - * - * @ingroup ErrorHandling - * - * @author Ross Harman, MHPCC - * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-03 01:30:32 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - */ - -#ifndef PM_FLATFIELD_ERRORS_H -#define PM_FLATFIELD_ERRORS_H - -/* N.B., lines between '//~Start' and '//~End' are automatic generated from - * the template following the '//~Start'. The template is used to generate - * the other lines by, for each error text in psDataManipErrors.dat, the following - * substitutions are made: - * $1 The error text macro name (first word in the psFlatFieldErrors.h lines) - * $2 The error text (rest of the line in psFlatFieldErrors.h) - * $n The order of the source line in psFlatFieldErrors.h (comments excluded) - * - * DO NOT EDIT THE LINES BETWEEN //~Start and //~End! ANY CHANGES WILL BE OVERWRITTEN. - */ - -#define PS_ERRORNAME_DOMAIN "psModule.src." - -//~Start #define PS_ERRORTEXT_$1 "$2" -#define PS_ERRORTEXT_pmFlatField_NULL_FLAT_READOUT "Null not allowed for flat readout." -#define PS_ERRORTEXT_pmFlatField_NULL_INPUT_IMAGE "Null not allowed for input image." -#define PS_ERRORTEXT_pmFlatField_NULL_FLAT_IMAGE "Null not allowed for flat image." -#define PS_ERRORTEXT_pmFlatField_SIZE_INPUT_IMAGE "Input image size exceeds that of flat image: (%d, %d) vs (%d, %d)" -#define PS_ERRORTEXT_pmFlatField_SIZE_MASK_IMAGE "Input image mask size exceeds that of flat image: (%d, %d) vs (%d, %d)" -#define PS_ERRORTEXT_pmFlatField_OFFSET_FLAT_IMAGE "Total offset >= flat image size: (%d, %d) vs (%d, %d)" -#define PS_ERRORTEXT_pmFlatField_OFFSET_INPUT_IMAGE "Total offset >= input image: (%d, %d) vs (%d, %d)" -#define PS_ERRORTEXT_pmFlatField_OFFSET_MASK_IMAGE "Total offset >= input image mask: (%d, %d) vs (%d, %d)" -#define PS_ERRORTEXT_pmFlatField_TYPE_INPUT_IMAGE "Complex types not allowed for input image. Type: %d" -#define PS_ERRORTEXT_pmFlatField_TYPE_FLAT_IMAGE "Complex types not allowed for flat image. Type: %d" -#define PS_ERRORTEXT_pmFlatField_TYPE_MASK_IMAGE "Mask must be PS_TYPE_MASK type. Type: %d" -#define PS_ERRORTEXT_pmFlatField_TYPE_MISMATCH "Input and flat image types differ: (%d vs %d)" -#define PS_ERRORTEXT_pmFlatField_TYPE_UNSUPPORTED "Unsupported image datatype. Type: %d" -//~End - -#endif Index: unk/archive/scripts/src/phase2/pmMaskBadPixels.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmMaskBadPixels.c (revision 5798) +++ (revision ) @@ -1,183 +1,0 @@ -/** @file pmMaskBadPixels.c - * - * @brief Given an input image, a bad pixel mask, a corresponding value in the bad pixel mask to mask, a - * saturation level, and a growing radius, mask in the input image those pixels in the bad pixel mask that - * match the value to mask. - * - * Given an input image, in, a bad pixel mask, a corresponding value in the bad pixel mask to mask in the - * input image, maskVal, a saturation level, and a growing radius, pmMaskBadPixels shall mask in the input - * image those pixels in the bad pixel mask that match the value to mask. Note that the input image, in, is - * modified in-place. All pixels in the mask which satisfy the maskVal shall have their corresponding pixels - * masked in the input image, in. All pixels which satisfy the growVal shall have their corresponding - * pixels, along with all pixels within the grow radius masked. Pixels which have flux greater than sat shall - * also be masked, but not grown. Note that the input image, in, and the mask need not be the same size, since - * the input image may already have been trimmed (following overscan subtraction), but the function shall use - * the offsets in the image (in->x0 and in->y0) to determine the appropriate offsets to obtain the correct - * pixel on the mask. In the event that the mask image is too small (i.e., pixels on the input image refer to - * pixels outside the range of the mask image), the function shall generate an error. - - * @author Ross Harman, MHPCC - * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-30 21:46:15 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - */ - -#if HAVE_CONFIG_H -#include -#endif - -#include -#include -#include - -#include "pmMaskBadPixels.h" -#include "pmMaskBadPixelsErrors.h" - -bool pmMaskBadPixels(pmReadout *in, const psImage *mask, unsigned int maskVal, float sat, - unsigned int growVal, int grow) -{ - int i = 0; - int j = 0; - int jj = 0; - int ii = 0; - int rRound = 0; - int rowMin = 0; - int rowMax = 0; - int colMin = 0; - int colMax = 0; - int totOffCol = 0; - int totOffRow = 0; - float r = 0.0f; - psElemType inType; - psElemType maskType; - psImage *inImage = NULL; - psImage *inMask = NULL; - - - // Check for nulls - if (in == NULL) { - return true; // Readout may not have data in it - } else if(mask==NULL) { - psError( PS_ERR_BAD_PARAMETER_NULL, true, - PS_ERRORTEXT_pmMaskBadPixels_NULL_MASK_IMAGE); - return false; - } - - inImage = in->image; - if (inImage == NULL) { - psError( PS_ERR_BAD_PARAMETER_NULL, true, - PS_ERRORTEXT_pmMaskBadPixels_NULL_INPUT_IMAGE); - return false; - } else if(in->mask == NULL) { - in->mask = psImageAlloc(inImage->numCols, inImage->numRows, PS_TYPE_MASK); - memset(in->mask->data.V[0], 0, inImage->numCols*inImage->numRows*PSELEMTYPE_SIZEOF(PS_TYPE_MASK)); - } - inMask = in->mask; - - // Check input image and its mask are not larger than mask - if(inImage->numRows > mask->numRows || inImage->numCols > mask->numCols) { - psError( PS_ERR_BAD_PARAMETER_SIZE, true, - PS_ERRORTEXT_pmMaskBadPixels_SIZE_INPUT_IMAGE, - inImage->numRows, inImage->numCols, mask->numRows, mask->numCols); - return false; - } else if(inMask->numRows>mask->numRows || inMask->numCols>mask->numCols) { - psError( PS_ERR_BAD_PARAMETER_SIZE, true, - PS_ERRORTEXT_pmMaskBadPixels_SIZE_MASK_IMAGE, - inMask->numRows, inMask->numCols, mask->numRows, mask->numCols); - return false; - } - - // Determine total offset based on image offset with chip offset - totOffCol = inImage->col0 + in->col0; - totOffRow = inImage->row0 + in->row0; - - // Check that offsets are within image limits - if(totOffRow>=mask->numRows || totOffCol>=mask->numCols) { - psError( PS_ERR_BAD_PARAMETER_SIZE, true, - PS_ERRORTEXT_pmMaskBadPixels_OFFSET_MASK_IMAGE, - totOffRow, totOffCol, mask->numRows, mask->numCols); - return false; - } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) { - psError( PS_ERR_BAD_PARAMETER_SIZE, true, - PS_ERRORTEXT_pmMaskBadPixels_OFFSET_INPUT_IMAGE, - totOffRow, totOffCol, inImage->numRows, inImage->numCols); - return false; - } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) { - psError( PS_ERR_BAD_PARAMETER_SIZE, true, - PS_ERRORTEXT_pmMaskBadPixels_OFFSET_INPUT_IMAGE_MASK, - totOffRow, totOffCol, inMask->numRows, inMask->numCols); - return false; - } - - // Check for incorrect types - inType = inImage->type.type; - maskType = mask->type.type; - if(PS_IS_PSELEMTYPE_COMPLEX(inType)) { - psError( PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_pmMaskBadPixels_TYPE_INPUT_IMAGE, - inType); - return false; - } else if(maskType!=PS_TYPE_MASK) { - psError( PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_pmMaskBadPixels_TYPE_MASK_IMAGE, - maskType); - return false; - } - - // Macro for all PS types - #define PM_BAD_PIXELS(TYPE) \ -case PS_TYPE_##TYPE: \ - for(j=totOffRow; jnumRows; j++) { \ - for(i=totOffCol; inumCols; i++) { \ - \ - /* Pixels with flux greater than sat shall be masked */ \ - if(inImage->data.TYPE[j][i] > sat) { \ - inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_SAT; \ - } \ - \ - /* Pixels which satisfy maskVal shall be masked */ \ - inMask->data.PS_TYPE_MASK_DATA[j][i] |= (mask->data.PS_TYPE_MASK_DATA[j][i]&maskVal); \ - \ - /* Pixels which satisfy growVal and within the grow radius shall be masked */ \ - if(mask->data.PS_TYPE_MASK_DATA[j][i] & growVal) { \ - rowMin = MAX(j-grow, 0); \ - rowMax = MIN(j+grow+1, inImage->numRows); \ - colMin = MAX(i-grow, 0); \ - colMax = MIN(i+grow+1, inImage->numCols); \ - for(jj=rowMin; jjdata.PS_TYPE_MASK_DATA[jj][ii] |= \ - (mask->data.PS_TYPE_MASK_DATA[j][i]&growVal); \ - } \ - } \ - } \ - } \ - } \ - } \ - break; - - // Switch to call bad pixel masking macro defined above - switch(inType) { - PM_BAD_PIXELS(U8); - PM_BAD_PIXELS(U16); - PM_BAD_PIXELS(U32); - PM_BAD_PIXELS(U64); - PM_BAD_PIXELS(S8); - PM_BAD_PIXELS(S16); - PM_BAD_PIXELS(S32); - PM_BAD_PIXELS(S64); - PM_BAD_PIXELS(F32); - PM_BAD_PIXELS(F64); - default: - psError( PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_pmMaskBadPixels_TYPE_UNSUPPORTED, - inType); - } - - return false; -} Index: unk/archive/scripts/src/phase2/pmMaskBadPixels.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmMaskBadPixels.h (revision 5798) +++ (revision ) @@ -1,61 +1,0 @@ -/** @file pmMaskBadPixels.h - * - * @brief Given an input image, a bad pixel mask, a corresponding value in the bad pixel mask to mask, a - * saturation level, and a growing radius, mask in the input image those pixels in the bad pixel mask that - * match the value to mask. - * - * Given an input image, in, a bad pixel mask, a corresponding value in the bad pixel mask to mask in the - * input image, maskVal, a saturation level, and a growing radius, pmMaskBadPixels shall mask in the input - * image those pixels in the bad pixel mask that match the value to mask. Note that the input image, in, is - * modified in-place. All pixels in the mask which satisfy the maskVal shall have their corresponding pixels - * masked in the input image, in. All pixels which satisfy the growVal shall have their corresponding - * pixels, along with all pixels within the grow radius masked. Pixels which have flux greater than sat shall - * also be masked, but not grown. Note that the input image, in, and the mask need not be the same size, - * since the input image may already have been trimmed (following overscan subtraction), but the function - * shall use the offsets in the image (in->x0 and in->y0) to determine the appropriate offsets to obtain the - * correct pixel on the mask. In the event that the mask image is too small (i.e., pixels on the input image - * refer to pixels outside the range of the mask image), the function shall generate an error. - * - * @author Ross Harman, MHPCC - * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-23 02:58:30 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - */ - -#include "pslib.h" -#include "pmFPA.h" - -/** Mask values */ -typedef enum { - PM_MASK_TRAP = 0x0001, ///< The pixel is a charge trap. - PM_MASK_BADCOL = 0x0002, ///< The pixel is a bad column. - PM_MASK_SAT = 0x0004, ///< The pixel is saturated. - PM_MASK_FLAT = 0x0008 ///< The pixel is non-positive in the flat-field. -} pmMaskValue; - -/** Macro to find maximum of two numbers */ -#define MAX(A,B)((A)>=(B)?(A):(B)) - -/** Macro to find minimum of two numbers */ -#define MIN(A,B)((A)<=(B)?(A):(B)) - - -/** Execute bad pixels module. - * - * Given an input image, a bad pixel mask, a corresponding value in the bad pixel mask to mask, a - * saturation level, and a growing radius, mask in the input image those pixels in the bad pixel mask that - * match the value to mask. - * - * @return bool: True or false for success or failure - */ -bool pmMaskBadPixels( - pmReadout *in, ///< Readout containing input image data. - const psImage *mask, ///< Mask data to be added to readout mask data. - unsigned int maskVal, ///< Mask value to determine what to add to input mask. - float sat, ///< Saturation limit to mask bad pixels. - unsigned int growVal, ///< Mask data to determine if a circurlar area should be masked. - int grow ///< Radius of mask to apply around pixel. -); - Index: unk/archive/scripts/src/phase2/pmMaskBadPixelsErrors.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmMaskBadPixelsErrors.h (revision 5798) +++ (revision ) @@ -1,45 +1,0 @@ -/** @file pmMaskBadPixelsErrors.h - * - * @brief Contains the error text for the mask bad pixels module - * - * @ingroup ErrorHandling - * - * @author Ross Harman, MHPCC - * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-23 02:58:30 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - */ - -#ifndef PM_FLATFIELD_ERRORS_H -#define PM_FLATFIELD_ERRORS_H - -/* N.B., lines between '//~Start' and '//~End' are automatic generated from - * the template following the '//~Start'. The template is used to generate - * the other lines by, for each error text in pmMaskBadPixelsErrors.dat, the following - * substitutions are made: - * $1 The error text macro name (first word in the pmMaskBadPixelsErrors.h lines) - * $2 The error text (rest of the line in pmMaskBadPixelsErrors.h) - * $n The order of the source line in pmMaskBadPixelsErrors.h (comments excluded) - * - * DO NOT EDIT THE LINES BETWEEN //~Start and //~End! ANY CHANGES WILL BE OVERWRITTEN. - */ - -#define PS_ERRORNAME_DOMAIN "psModule.src." - -//~Start #define PS_ERRORTEXT_$1 "$2" -#define PS_ERRORTEXT_pmMaskBadPixels_NULL_MASK_IMAGE "Null not allowed for mask image." -#define PS_ERRORTEXT_pmMaskBadPixels_NULL_INPUT_IMAGE "Null not allowed for input image." -#define PS_ERRORTEXT_pmMaskBadPixels_SIZE_INPUT_IMAGE "Input image size exceeds that of mask image: (%d, %d) vs (%d, %d)" -#define PS_ERRORTEXT_pmMaskBadPixels_SIZE_MASK_IMAGE "Input image mask size exceeds that of mask image: (%d, %d) vs (%d, %d)" -#define PS_ERRORTEXT_pmMaskBadPixels_OFFSET_MASK_IMAGE "Total offset >= mask image: (%d, %d) vs (%d, %d)" -#define PS_ERRORTEXT_pmMaskBadPixels_OFFSET_INPUT_IMAGE "Total offset >= input image: (%d, %d) vs (%d, %d)" -#define PS_ERRORTEXT_pmMaskBadPixels_OFFSET_INPUT_IMAGE_MASK "Total offset >= input image mask: (%d, %d) vs (%d, %d)" -#define PS_ERRORTEXT_pmMaskBadPixels_TYPE_INPUT_IMAGE "Complex types not allowed for input image. Type: %d" -#define PS_ERRORTEXT_pmMaskBadPixels_TYPE_MASK_IMAGE "Mask must be PS_TYPE_MASK type. Type: %d" -#define PS_ERRORTEXT_pmMaskBadPixels_TYPE_MISMATCH "Input and flat image types differ: (%d vs %d)" -#define PS_ERRORTEXT_pmMaskBadPixels_TYPE_UNSUPPORTED "Unsupported image datatype. Type: %d" -//~End - -#endif Index: unk/archive/scripts/src/phase2/pmNonLinear.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmNonLinear.c (revision 5798) +++ (revision ) @@ -1,123 +1,0 @@ -/** @file pmNonLinear.c - * - * Provides polynomial or table lookup non-linearity corrections to readouts. - * - * @author GLG, MHPCC - * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-03 03:21:36 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - * XXX: The SDR is silent about image types. Only F32 was implemented. - * - */ - -#if HAVE_CONFIG_H -#include -#endif - -#include -#include - -#include "pmNonLinear.h" - -/****************************************************************************** -pmNonLinearityLookup(): This routine will take an pmReadout image as input -and a 1-D polynomial. For each pixel in the input image, the polynomial will -be evaluated at that pixels value, and the image pixel will then be set to -that value. - *****************************************************************************/ - -pmReadout *pmNonLinearityPolynomial(pmReadout *inputReadout, - const psPolynomial1D *input1DPoly) -{ - PS_ASSERT_PTR_NON_NULL(inputReadout, NULL); - PS_ASSERT_PTR_NON_NULL(inputReadout->image, NULL); - PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL); - PS_ASSERT_PTR_NON_NULL(input1DPoly, NULL); - - psS32 i; - psS32 j; - - for (i=0;iimage->numRows;i++) { - for (j=0;jimage->numCols;j++) { - inputReadout->image->data.F32[i][j] = psPolynomial1DEval(input1DPoly, inputReadout->image->data.F32[i][j]); - } - } - return(inputReadout); -} - - -/****************************************************************************** -pmNonLinearityLookup(): This routine will take an pmReadout image as input -and two input vectors, which constitute a lookup table. For each pixel in -the input image, that pixels value will be determined in the input vector -inFluxe, and the corresponding value in outFlux. The image pixel will then -be set to the value from outFlux. - *****************************************************************************/ -pmReadout *pmNonLinearityLookup(pmReadout *inputReadout, - const psVector *inFlux, - const psVector *outFlux) -{ - PS_ASSERT_PTR_NON_NULL(inputReadout,NULL); - PS_ASSERT_PTR_NON_NULL(inputReadout->image,NULL); - PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL); - PS_ASSERT_PTR_NON_NULL(inFlux,NULL); - if (inFlux->n < 2) { - psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): input vector less than 2 elements. Returning inputReadout image."); - return(inputReadout); - } - PS_ASSERT_PTR_NON_NULL(outFlux,NULL); - psS32 tableSize = inFlux->n; - if (inFlux->n != outFlux->n) { - tableSize = PS_MIN(inFlux->n, outFlux->n); - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmNonLinear.c: pmNonLinearityLookup(): input vectors have different sizes (%d, %d)\n", inFlux->n, outFlux->n); - } - PS_ASSERT_VECTOR_TYPE(inFlux, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_TYPE(outFlux, PS_TYPE_F32, NULL); - - psS32 i; - psS32 j; - psS32 binNum; - psScalar x; - psS32 numPixels = 0; - psF32 slope; - - x.type.type = PS_TYPE_F32; - for (i=0;iimage->numRows;i++) { - for (j=0;jimage->numCols;j++) { - x.data.F32 = inputReadout->image->data.F32[i][j]; - binNum = p_psVectorBinDisect((psVector *)inFlux, &x); - - if (binNum == -2) { - // We get here if x is below the table lookup range. - inputReadout->image->data.F32[i][j] = outFlux->data.F32[0]; - numPixels++; - - } else if (binNum == -1) { - // We get here if x is above the table lookup range. - inputReadout->image->data.F32[i][j] = outFlux->data.F32[tableSize-1]; - numPixels++; - - } else if (binNum < -2) { - // We get here if there was some other problem. - psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): Could not perform p_psVectorBinDisect(). Returning inputReadout image."); - return(inputReadout); - numPixels++; - } else { - // Perform linear interpolation. - slope = (outFlux->data.F32[binNum+1] - outFlux->data.F32[binNum]) / - (inFlux->data.F32[binNum+1] - inFlux->data.F32[binNum]); - inputReadout->image->data.F32[i][j] = outFlux->data.F32[binNum] + - ((x.data.F32 - inFlux->data.F32[binNum]) * slope); - } - } - } - if (numPixels > 0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmNonLinear.c: pmNonLinearityLookup(): %d pixels outside table.", numPixels); - } - return(inputReadout); -} Index: unk/archive/scripts/src/phase2/pmNonLinear.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmNonLinear.h (revision 5798) +++ (revision ) @@ -1,27 +1,0 @@ -/** @file pmNonLinear.h - * - * Provides polynomial or table lookup non-linearity corrections to readouts. - * - * @author GLG, MHPCC - * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-03 03:21:36 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#if !defined(PM_NON_LINEAR_H) -#define PM_NON_LINEAR_H - -#include "pslib.h" -#include "pmFPA.h" // #include "pmAstrometry.h" - -pmReadout *pmNonLinearityPolynomial(pmReadout *in, - const psPolynomial1D *coeff); - -pmReadout *pmNonLinearityLookup(pmReadout *in, - const psVector *inFlux, - const psVector *outFlux); - -#endif Index: unk/archive/scripts/src/phase2/pmSubtractBias.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmSubtractBias.c (revision 5798) +++ (revision ) @@ -1,685 +1,0 @@ -/** @file pmSubtractBias.c - * - * This file will contain a module which will subtract the detector bias - * in place from an input image. - * - * @author GLG, MHPCC - * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-12-14 03:47:35 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#if HAVE_CONFIG_H -#include -#endif - -#include "pmSubtractBias.h" - -#define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2 -#define PM_SUBTRACT_BIAS_SPLINE_ORDER 3 - -// XXX: put these in psConstants.h -void PS_POLY1D_PRINT(psPolynomial1D *poly) -{ - printf("-------------- PS_POLY1D_PRINT() --------------\n"); - printf("poly->nX is %d\n", poly->nX); - for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) { - printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]); - } -} - -void PS_PRINT_SPLINE(psSpline1D *mySpline) -{ - printf("-------------- PS_PRINT_SPLINE() --------------\n"); - printf("mySpline->n is %d\n", mySpline->n); - for (psS32 i = 0 ; i < mySpline->n ; i++) { - PS_POLY1D_PRINT(mySpline->spline[i]); - } - PS_VECTOR_PRINT_F32(mySpline->knots); -} - -#define PS_IMAGE_PRINT_F32_HIDEF(NAME) \ -printf("======== printing %s ========\n", #NAME); \ -for (int i = 0 ; i < (NAME)->numRows ; i++) { \ - for (int j = 0 ; j < (NAME)->numCols ; j++) { \ - printf("%.5f ", (NAME)->data.F32[i][j]); \ - } \ - printf("\n"); \ -}\ - -/****************************************************************************** -psSubtractFrame(): this routine will take as input a readout for the input -image and a readout for the bias image. The bias image is subtracted in -place from the input image. -*****************************************************************************/ -static pmReadout *SubtractFrame(pmReadout *in, - const pmReadout *bias) -{ - psS32 i; - psS32 j; - - if (bias == NULL) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.c: SubtractFrame(): bias frame is NULL. Returning original image.\n"); - return(in); - } - - - if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) { - psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows. Returning in image\n"); - return(in); - } - if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) { - psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns. Returning in image\n"); - return(in); - } - - for (i=0;iimage->numRows;i++) { - for (j=0;jimage->numCols;j++) { - in->image->data.F32[i][j]-= - bias->image->data.F32[i+in->row0-bias->row0][j+in->col0-bias->col0]; - if ((in->mask != NULL) && (bias->mask != NULL)) { - (in->mask->data.U8[i][j])|= - bias->mask->data.U8[i+in->row0-bias->row0][j+in->col0-bias->col0]; - } - } - } - - return(in); -} - -/****************************************************************************** -ImageSubtractScalar(): subtract a scalar from the input image. - -XXX: Use a psLib function for this. - -XXX: This should - *****************************************************************************/ -static psImage *ImageSubtractScalar(psImage *image, - psF32 scalar) -{ - for (psS32 i=0;inumRows;i++) { - for (psS32 j=0;jnumCols;j++) { - image->data.F32[i][j]-= scalar; - } - } - return(image); -} - -/****************************************************************************** -GenNewStatOptions(): this routine will take as input the options member of the -stat data structure, determine if multiple options have been specified, issue -a warning message if so, and return the highest priority option (according to -the order of the if-statements in this code). The higher priority options are -listed lower in the code. - *****************************************************************************/ -static psStatsOptions GenNewStatOptions(const psStats *stat) -{ - psS32 numOptions = 0; - psStatsOptions opt = 0; - - if (stat->options & PS_STAT_ROBUST_MODE) { - if (numOptions == 0) { - opt = PS_STAT_ROBUST_MODE; - } - numOptions++; - } - - if (stat->options & PS_STAT_ROBUST_MEDIAN) { - if (numOptions == 0) { - opt = PS_STAT_ROBUST_MEDIAN; - } - numOptions++; - } - - if (stat->options & PS_STAT_ROBUST_MEAN) { - if (numOptions == 0) { - opt = PS_STAT_ROBUST_MEAN; - } - numOptions++; - } - - if (stat->options & PS_STAT_CLIPPED_MEAN) { - if (numOptions == 0) { - opt = PS_STAT_CLIPPED_MEAN; - } - numOptions++; - } - - if (stat->options & PS_STAT_SAMPLE_MEDIAN) { - if (numOptions == 0) { - opt = PS_STAT_SAMPLE_MEDIAN; - } - numOptions++; - } - - if (stat->options & PS_STAT_SAMPLE_MEAN) { - numOptions++; - opt = PS_STAT_SAMPLE_MEAN; - } - - - if (numOptions == 0) { - psError(PS_ERR_UNKNOWN,true, "No statistics options have been specified.\n"); - } - if (numOptions != 1) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.c: GenNewStatOptions(): Too many statistics options have been specified\n"); - } - return(opt); -} - - - -/****************************************************************************** -ScaleOverscanVector(): this routine takes as input an arbitrary vector, -creates a new vector of length n, and fills the new vector with the -interpolated values of the old vector. The type of interpolation is: - PM_FIT_POLYNOMIAL: fit a polynomial to the entire input vector data. - PM_FIT_SPLINE: fit splines to the input vector data. -XXX: Doesn't it make more sense to do polynomial interpolation on a few -elements of the input vector, rather than fit a polynomial to the entire -vector? - *****************************************************************************/ -static psVector *ScaleOverscanVector(psVector *overscanVector, - psS32 n, - void *fitSpec, - pmFit fit) -{ - psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4, - "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n); - // PS_VECTOR_PRINT_F32(overscanVector); - - if (NULL == overscanVector) { - return(overscanVector); - } - - // Allocate the new vector. - psVector *newVec = psVectorAlloc(n, PS_TYPE_F32); - - // - // If the new vector is the same size as the old, simply copy the data. - // - if (n == overscanVector->n) { - for (psS32 i = 0 ; i < n ; i++) { - newVec->data.F32[i] = overscanVector->data.F32[i]; - } - return(newVec); - } - psPolynomial1D *myPoly; - psSpline1D *mySpline; - psF32 x; - psS32 i; - if (fit == PM_FIT_POLYNOMIAL) { - // Fit a polynomial to the old overscan vector. - myPoly = (psPolynomial1D *) fitSpec; - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL); - if (myPoly == NULL) { - psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(1): Could not fit a polynomial to the psVector.\n"); - return(NULL); - } - - // For each element of the new vector, convert the x-ordinate to that - // of the old vector, use the fitted polynomial to determine the - // interpolated value at that point, and set the new vector. - for (i=0;in) / ((psF32) n); - newVec->data.F32[i] = psPolynomial1DEval(myPoly, x); - } - } else if (fit == PM_FIT_SPLINE) { - psS32 mustFreeSpline = 0; - // Fit a spline to the old overscan vector. - mySpline = (psSpline1D *) fitSpec; - // XXX: Does it make any sense to have a psSpline argument? - if (mySpline == NULL) { - mustFreeSpline = 1; - } - - // - // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL, - // splines endpoints will be from 0.0 to overscanVector->n-1. Must scale - // properly when doing the spline eval. - // - // mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL); - mySpline = psVectorFitSpline1D(NULL, overscanVector); - if (mySpline == NULL) { - psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(2): Could not fit a spline to the psVector.\n"); - return(NULL); - } - // PS_PRINT_SPLINE(mySpline); - - // For each element of the new vector, convert the x-ordinate to that - // of the old vector, use the fitted polynomial to determine the - // interpolated value at that point, and set the new vector. - for (i=0;in - 1] - x = ((psF32) i) * ((psF32) (overscanVector->n-1)) / ((psF32) n); - newVec->data.F32[i] = psSpline1DEval(mySpline, x); - } - if (mustFreeSpline ==1) { - psFree(mySpline); - } - // PS_VECTOR_PRINT_F32(newVec); - - - } else { - psError(PS_ERR_UNKNOWN, true, "unknown fit type. Returning NULL.\n"); - psFree(newVec); - return(NULL); - } - - psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4, - "---- ScaleOverscanVector() exit ----\n"); - return(newVec); -} - -/****************************************************************************** -XXX: The SDRS does not specify type support. F32 is implemented here. - *****************************************************************************/ -pmReadout *pmSubtractBias(pmReadout *in, - void *fitSpec, - const psList *overscans, - pmOverscanAxis overScanAxis, - psStats *stat, - psS32 nBinOrig, - pmFit fit, - const pmReadout *bias) -{ - psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4, - "---- pmSubtractBias() begin ----\n"); - PS_ASSERT_READOUT_NON_NULL(in, NULL); - PS_ASSERT_READOUT_NON_EMPTY(in, NULL); - PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL); - - // - // If the overscans != NULL, then check the type of each image. - // - if (overscans != NULL) { - psListElem *tmpOverscan = (psListElem *) overscans->head; - while (NULL != tmpOverscan) { - psImage *myOverscanImage = (psImage *) tmpOverscan->data; - PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL); - tmpOverscan = tmpOverscan->next; - } - } - - if ((overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE)) { - psError(PS_ERR_UNKNOWN,true, "(overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE). Returning in image\n"); - return(in); - } - - // Check for an unallowable pmFit. - if ((fit != PM_OVERSCAN_NONE) && - (fit != PM_OVERSCAN_ROWS) && - (fit != PM_OVERSCAN_COLUMNS) && - (fit != PM_OVERSCAN_ALL)) { - psError(PS_ERR_UNKNOWN, true, "fit is unallowable (%d). Returning in image.\n", fit); - return(in); - } - // Check for an unallowable pmOverscanAxis. - if ((overScanAxis != PM_OVERSCAN_NONE) && - (overScanAxis != PM_OVERSCAN_ROWS) && - (overScanAxis != PM_OVERSCAN_COLUMNS) && - (overScanAxis != PM_OVERSCAN_ALL)) { - psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d). Returning in image.\n", overScanAxis); - return(in); - } - psS32 i; - psS32 j; - psS32 numBins = 0; - static psVector *overscanVector = NULL; - psVector *tmpRow = NULL; - psVector *tmpCol = NULL; - psVector *myBin = NULL; - psVector *binVec = NULL; - psListElem *tmpOverscan = NULL; - double statValue; - psImage *myOverscanImage = NULL; - psPolynomial1D *myPoly = NULL; - psSpline1D *mySpline = NULL; - psS32 nBin; - - // - // Create a static stats data structure and determine the highest - // priority stats option. - // - static psStats *myStats = NULL; - if (myStats == NULL) { - myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); - p_psMemSetPersistent(myStats, true); - } - if (stat != NULL) { - myStats->options = GenNewStatOptions(stat); - } - - - if (overScanAxis == PM_OVERSCAN_NONE) { - if (fit != PM_FIT_NONE) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): overScanAxis equals NONE, and fit does not equal NONE. Proceeding to full fram subtraction.\n"); - } - - if (overscans != NULL) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): overScanAxis equals NONE and overscans does not equal NULL. Proceeding to full fram subtraction.\n"); - } - return(SubtractFrame(in, bias)); - } - - if ((overScanAxis == PM_OVERSCAN_ALL) && (fit != PM_FIT_NONE)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): overScanAxis equals ALL, and fit does not equal NONE. Proceeding with the rest of the module.\n"); - } - - - // - // We subtract each overscan region from the image data. - // If we get here we know that overscans != NULL. - // - - if (overScanAxis == PM_OVERSCAN_ALL) { - tmpOverscan = (psListElem *) overscans->head; - while (NULL != tmpOverscan) { - myOverscanImage = (psImage *) tmpOverscan->data; - - PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL); - psStats *rc = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff); - if (rc == NULL) { - psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation. Returning in image.\n"); - return(in); - } - if (false == p_psGetStatValue(myStats, &statValue)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); - return(in); - } - ImageSubtractScalar(in->image, statValue); - - tmpOverscan = tmpOverscan->next; - } - return(in); - } - - // This check is redundant with above code. - if (!((overScanAxis == PM_OVERSCAN_ROWS) || (overScanAxis == PM_OVERSCAN_COLUMNS))) { - psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d).\nReturning in image.\n", overScanAxis); - return(in); - } - - tmpOverscan = (psListElem *) overscans->head; - while (NULL != tmpOverscan) { - // PS_IMAGE_PRINT_F32_HIDEF(in->image); - myOverscanImage = (psImage *) tmpOverscan->data; - - if (overScanAxis == PM_OVERSCAN_ROWS) { - if (myOverscanImage->numCols != (in->image)->numCols) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): overscan image has %d columns, input image has %d columns\n", - myOverscanImage->numCols, in->image->numCols); - } - - // We create a row vector and subtract this vector from image. - // XXX: Is there a better way to extract a psVector from a psImage without - // having to copy every element in that vector? - overscanVector = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32); - for (i=0;in;i++) { - overscanVector->data.F32[i] = 0.0; - } - tmpRow = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32); - - // For each column of the input image, loop through every row, - // collect the pixel in that row, then performed the specified - // statistical op on those pixels. Store this in overscanVector. - for (i=0;inumCols;i++) { - for (j=0;jnumRows;j++) { - tmpRow->data.F32[j] = myOverscanImage->data.F32[j][i]; - } - psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0); - if (rc == NULL) { - psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); - return(in); - } - if (false == p_psGetStatValue(rc, &statValue)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); - return(in); - } - overscanVector->data.F32[i] = statValue; - } - psFree(tmpRow); - - // Scale the overscan vector to the size of the input image. - if (overscanVector->n != in->image->numCols) { - if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) { - psVector *newVec = ScaleOverscanVector(overscanVector, - in->image->numCols, - fitSpec, fit); - if (newVec == NULL) { - psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector. Returning in image.\n"); - return(in); - } - psFree(overscanVector); - overscanVector = newVec; - } else { - psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector. Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL. Returning in image.\n"); - psFree(overscanVector); - return(in); - } - } - } - - if (overScanAxis == PM_OVERSCAN_COLUMNS) { - if (myOverscanImage->numRows != (in->image)->numRows) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): overscan image has %d rows, input image has %d rows\n", - myOverscanImage->numRows, in->image->numRows); - } - - // We create a column vector and subtract this vector from image. - overscanVector = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32); - for (i=0;in;i++) { - overscanVector->data.F32[i] = 0.0; - } - tmpCol = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32); - - // For each row of the input image, loop through every column, - // collect the pixel in that row, then performed the specified - // statistical op on those pixels. Store this in overscanVector. - for (i=0;inumRows;i++) { - for (j=0;jnumCols;j++) { - tmpCol->data.F32[j] = myOverscanImage->data.F32[i][j]; - } - psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0); - if (rc == NULL) { - psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); - return(in); - } - if (false == p_psGetStatValue(rc, &statValue)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); - return(in); - } - overscanVector->data.F32[i] = statValue; - } - psFree(tmpCol); - - // Scale the overscan vector to the size of the input image. - if (overscanVector->n != in->image->numRows) { - if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) { - psVector *newVec = ScaleOverscanVector(overscanVector, - in->image->numRows, - fitSpec, fit); - if (newVec == NULL) { - psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector. Returning in image.\n"); - return(in); - } - psFree(overscanVector); - overscanVector = newVec; - } else { - psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector. Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL. Returning in image.\n"); - psFree(overscanVector); - return(in); - } - } - } - - // - // Re-bin the overscan vector (change its length). - // - // Only if nBinOrig > 1. - if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) { - numBins = 1+((overscanVector->n)/nBinOrig); - myBin = psVectorAlloc(numBins, PS_TYPE_F32); - binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32); - - for (i=0;in > ((i*nBinOrig)+j)) { - binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j]; - } else { - // XXX: we get here if nBinOrig does not evenly divide - // the overscanVector vector. This is the last bin. Should - // we change the binVec->n to acknowledge that? - binVec->n = j; - } - } - psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0); - if (rc == NULL) { - psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); - return(in); - } - if (false == p_psGetStatValue(rc, &statValue)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); - return(in); - } - myBin->data.F32[i] = statValue; - } - - // Change the effective size of overscanVector. - overscanVector->n = numBins; - for (i=0;idata.F32[i] = myBin->data.F32[i]; - } - psFree(binVec); - psFree(myBin); - nBin = nBinOrig; - } else { - nBin = 1; - } - - // At this point the number of data points in overscanVector should be - // equal to the number of rows/columns (whatever is appropriate) in the - // image divided by numBins. - // - - - // - // This doesn't seem right. The only way to do a spline fit is if, - // by SDRS requirements, fitSpec is not-NULL> But in order for it - // to be non-NULL, someone must have called psSpline1DAlloc() with - // the min, max, and number of splines. - // - if (!((fitSpec == NULL) || (fit == PM_FIT_NONE))) { - // - // Fit a polynomial or spline to the overscan vector. - // - if (fit == PM_FIT_POLYNOMIAL) { - myPoly = (psPolynomial1D *) fitSpec; - myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL); - if (myPoly == NULL) { - psError(PS_ERR_UNKNOWN, false, "(3) Could not fit a polynomial to overscan vector. Returning in image.\n"); - psFree(overscanVector); - return(in); - } - } else if (fit == PM_FIT_SPLINE) { - // XXX: This makes no sense - // XXX: must free mySpline? - mySpline = (psSpline1D *) fitSpec; - mySpline = psVectorFitSpline1D(NULL, overscanVector); - if (mySpline == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector. Returning in image.\n"); - psFree(overscanVector); - return(in); - } - } - - // - // Subtract fitted overscan vector row-wise from the image. - // - if (overScanAxis == PM_OVERSCAN_ROWS) { - for (i=0;i<(in->image)->numCols;i++) { - psF32 tmpF32 = 0.0; - if (fit == PM_FIT_POLYNOMIAL) { - tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin)); - } else if (fit == PM_FIT_SPLINE) { - tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin)); - } - for (j=0;j<(in->image)->numRows;j++) { - (in->image)->data.F32[j][i]-= tmpF32; - } - } - } - - // - // Subtract fitted overscan vector column-wise from the image. - // - if (overScanAxis == PM_OVERSCAN_COLUMNS) { - for (i=0;i<(in->image)->numRows;i++) { - psF32 tmpF32 = 0.0; - if (fit == PM_FIT_POLYNOMIAL) { - tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin)); - } else if (fit == PM_FIT_SPLINE) { - tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin)); - } - - for (j=0;j<(in->image)->numCols;j++) { - (in->image)->data.F32[i][j]-= tmpF32; - } - } - } - } else { - // - // If we get here, then no polynomials were fit to the overscan - // vector. We simply subtract it, taking into account binning, - // from the image. - // - - // - // Subtract overscan vector row-wise from the image. - // - if (overScanAxis == PM_OVERSCAN_ROWS) { - for (i=0;i<(in->image)->numCols;i++) { - for (j=0;j<(in->image)->numRows;j++) { - (in->image)->data.F32[j][i]-= overscanVector->data.F32[i/nBin]; - } - } - } - - // - // Subtract overscan vector column-wise from the image. - // - if (overScanAxis == PM_OVERSCAN_COLUMNS) { - for (i=0;i<(in->image)->numRows;i++) { - for (j=0;j<(in->image)->numCols;j++) { - (in->image)->data.F32[i][j]-= overscanVector->data.F32[i/nBin]; - } - } - } - } - - psFree(overscanVector); - - tmpOverscan = tmpOverscan->next; - } - - psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4, - "---- pmSubtractBias() exit ----\n"); - - if (bias != NULL) { - return(SubtractFrame(in, bias)); - } - return(in); -} - - Index: unk/archive/scripts/src/phase2/pmSubtractBias.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmSubtractBias.h (revision 5798) +++ (revision ) @@ -1,50 +1,0 @@ -/** @file pmSubtractBias.h - * - * This file will contain a module which will subtract the detector bias - * in place from an input image. - * - * @author GLG, MHPCC - * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-03 01:30:32 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#if !defined(PM_SUBTRACT_BIAS_H) -#define PM_SUBTRACT_BIAS_H - -#if HAVE_CONFIG_H -#include -#endif - -#include -#include -#include "pslib.h" - -#include "pmFPA.h" //#include "pmAstrometry.h" - -typedef enum { - PM_OVERSCAN_NONE, ///< No overscan subtraction - PM_OVERSCAN_ROWS, ///< Subtract rows - PM_OVERSCAN_COLUMNS, ///< Subtract columns - PM_OVERSCAN_ALL ///< Subtract the statistic of all pixels in overscan region -} pmOverscanAxis; - -typedef enum { - PM_FIT_NONE, ///< No fit - PM_FIT_POLYNOMIAL, ///< Fit polynomial - PM_FIT_SPLINE ///< Fit cubic splines -} pmFit; - -pmReadout *pmSubtractBias(pmReadout *in, ///< The input pmReadout image - void *fitSpec, ///< A polynomial or spline, defining the fit type. - const psList *overscans, ///< A psList of overscan images - pmOverscanAxis overScanAxis, ///< Defines how overscans are applied - psStats *stat, ///< The statistic to be used in combining overscan data - int nBin, ///< The amount of binning to be done image pixels. - pmFit fit, ///< PM_FIT_SPLINE, PM_FIT_POLYNOMIAL, or PM_FIT_NONE - const pmReadout *bias); ///< A possibly NULL bias pmReadout which is to be subtracted - -#endif