Index: trunk/psModules/src/camera/Makefile.am =================================================================== --- trunk/psModules/src/camera/Makefile.am (revision 6872) +++ trunk/psModules/src/camera/Makefile.am (revision 7017) @@ -4,19 +4,39 @@ libpsmodulecamera_la_LDFLAGS = -release $(PACKAGE_VERSION) libpsmodulecamera_la_SOURCES = \ + pmFPA.c \ + pmFPAConstruct.c \ + pmFPACopy.c \ + pmFPAHeader.c \ + pmFPARead.c \ + pmFPAUtils.c \ + pmFPAWrite.c \ + pmHDU.c \ + pmHDUUtils.c \ + pmReadout.c \ pmChipMosaic.c \ - pmFPAConceptsGet.c \ - pmFPAConceptsSet.c \ - pmFPAConstruct.c \ - pmFPAMorph.c \ - pmFPARead.c \ - pmFPAWrite.c + pmFPA_JPEG.c \ + pmFPAview.c \ + pmFPAfile.c + +# pmFPAMaskWeight.c \ +# pmChipMosaic.c psmoduleincludedir = $(includedir) psmoduleinclude_HEADERS = \ + pmFPA.h \ + pmFPAConstruct.h \ + pmFPACopy.h \ + pmFPAHeader.h \ + pmFPARead.h \ + pmFPAUtils.h \ + pmFPAWrite.h \ + pmHDU.h \ + pmHDUUtils.h \ + pmReadout.h \ pmChipMosaic.h \ - pmFPAConceptsGet.h \ - pmFPAConceptsSet.h \ - pmFPAConstruct.h \ - pmFPAMorph.h \ - pmFPARead.h \ - pmFPAWrite.h + pmFPA_JPEG.h \ + pmFPAview.h \ + pmFPAfile.h + +# pmFPAMaskWeight.h \ +# pmChipMosaic.h Index: trunk/psModules/src/camera/pmChipMosaic.c =================================================================== --- trunk/psModules/src/camera/pmChipMosaic.c (revision 7017) +++ trunk/psModules/src/camera/pmChipMosaic.c (revision 7017) @@ -0,0 +1,634 @@ +#include +#include +#include "pslib.h" +#include "pmFPA.h" +#include "pmHDU.h" +#include "psRegionIsBad.h" + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// File-static (private) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Do two regions overlap? +#define REGIONS_OVERLAP(region1, region2) \ +((region1->x0 > region2->x0 && region1->x0 < region2->x1) || \ + (region1->x1 > region2->x0 && region1->x1 < region2->x1) || \ + (region1->y0 > region2->y0 && region1->y0 < region2->y1) || \ + (region1->y1 > region2->y0 && region1->y1 < region2->y1)) + +// Compare a value with a maximum and minimum +#define COMPARE(value,min,max) \ +if ((value) < (min)) { \ + (min) = (value); \ +} \ +if ((value) > (max)) { \ + (max) = (value); \ +} + +// Update a metadata entry directly +#define MD_UPDATE(MD, NAME, TYPE, VALUE) \ +{ \ + psMetadataItem *item = psMetadataLookup(MD, NAME); \ + item->data.TYPE = VALUE; \ +} + +// Are the pixels for the chip contiguous on the HDU? +// Work this out by examining all the CELL.TRIMSEC and CELL.BIASSEC regions for the component cells +static bool chipContiguous(psRegion *bounds, // The bounds of the image, altered if primary==true + pmChip *chip, // The chip to examine for contiguity + bool primary // Is this the primary chip of interest? + ) +{ + if (primary) { + *bounds = psRegionSet(INFINITY, 0, INFINITY, 0); + } + psArray *cells = chip->cells; // The array of cells + bool mdok = true; // Status of MD lookup + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // Cell of interest + psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section + if (!mdok || !trimsec || psRegionIsBad(*trimsec)) { + psError(PS_ERR_IO, true, "CELL.TRIMSEC hasn't been set for cell %d.\n", i); + return false; + } + + if (primary) { + if (trimsec->x0 < bounds->x0) { + bounds->x0 = trimsec->x0; + } + if (trimsec->x1 > bounds->x1) { + bounds->x1 = trimsec->x1; + } + if (trimsec->y0 < bounds->y0) { + bounds->y0 = trimsec->y0; + } + if (trimsec->y1 > bounds->y1) { + bounds->y1 = trimsec->y1; + } + } else if (REGIONS_OVERLAP(trimsec, bounds)) { + return false; + } + + psList *biassecs = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.BIASSEC"); // Bias sections + if (!mdok || !biassecs) { + psError(PS_ERR_IO, true, "CELL.BIASSEC hasn't been set for cell %d.\n", i); + return false; + } + if (biassecs->n == 0) { + // No point allocating an iterator if there's nothing there to iterate on + continue; + } + psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); // Iterator + psRegion *biassec = NULL; // Bias section from iteration + while ((biassec = psListGetAndIncrement(biassecsIter))) { + if (psRegionIsBad(*biassec)) { + continue; + } + if (REGIONS_OVERLAP(biassec, bounds)) { + psFree(biassecsIter); + return false; + } + } + psFree(biassecsIter); + } + + // If we've gotten this far, everything is fine. + return true; +} + + +// Is the chip "nice"? If so, return the region containing the chip pixels +static psRegion *niceChip(int *xBinChip, int *yBinChip, // Binning for chip, to be returned + pmChip *chip // Chip to examine for "niceness". + ) +{ + // Check that we've got the HDU in the chip or the FPA + if ((!chip->hdu || !chip->hdu->images) && (!chip->parent->hdu || !chip->parent->hdu->images)) { + return NULL; + } + + // Check parity and binning for component cells + bool mdok = true; // Status of MD lookup + *xBinChip = 0; + *yBinChip = 0; + for (int i = 0; i < chip->cells->n; i++) { + pmCell *cell = chip->cells->data[i]; // The cell of interest + + // A "nice" chip must have only a single readout + if (cell->readouts->n != 1) { + return NULL; + } + + // A "nice" chip must have parity == 1 + int xParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XPARITY"); // Parity in x + if (!mdok || xParity == 0) { + psError(PS_ERR_IO, true, "CELL.XPARITY hasn't been set for cell %d.\n", i); + return NULL; + } + if (xParity != 1) { + return NULL; + } + int yParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YPARITY"); // Parity in y + if (!mdok || yParity == 0) { + psError(PS_ERR_IO, true, "CELL.YPARITY hasn't been set for cell %d.\n", i); + return NULL; + } + if (yParity != 1) { + return NULL; + } + + // A "nice" chip must have consistent binning + int xBin = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XBIN"); // Binning in x + if (!mdok || xBin <= 0) { + psError(PS_ERR_IO, true, "CELL.XPARITY hasn't been set for cell %d.\n", i); + return NULL; + } + int yBin = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YBIN"); // Binning in y + if (!mdok || yBin <= 0) { + psError(PS_ERR_IO, true, "CELL.YPARITY hasn't been set for cell %d.\n", i); + return NULL; + } + if (*xBinChip == 0 || *yBinChip == 0) { + *xBinChip = xBin; + *yBinChip = yBin; + } else if (xBin != *xBinChip || yBin != *yBinChip) { + return NULL; + } + } + + // Now check that the pixels are all contiguous + psRegion *imageBounds = psRegionAlloc(0, 0, 0, 0); // Bound of image on HDU + if (!chipContiguous(imageBounds, chip, true)) { + psTrace(__func__, 5, "Image isn't contiguous.\n"); + psFree(imageBounds); + return NULL; + } + + psString region = psRegionToString(*imageBounds); + psTrace(__func__, 7, "Image bounds: %s\n", region); + psFree(region); + + for (int i = 0; i < chip->cells->n; i++) { + pmCell *cell = chip->cells->data[i]; // The cell of interest + + // A "nice" chip must have the (0,0) pixel at CELL.X0,CELL.Y0 + int x0 = psMetadataLookupS32(&mdok, cell->concepts, "CELL.X0"); // Position of (0,0) on chip + if (!mdok) { + psError(PS_ERR_IO, true, "CELL.X0 hasn't been set for cell %d.\n", i); + return NULL; + } + int y0 = psMetadataLookupS32(&mdok, cell->concepts, "CELL.Y0"); // Position of (0,0) on chip + if (!mdok) { + psError(PS_ERR_IO, true, "CELL.Y0 hasn't been set for cell %d.\n", i); + return NULL; + } + pmReadout *readout = cell->readouts->data[0]; // A representative readout + if (!readout) { + return NULL; // Nothing here + } + if (x0 != readout->col0 + readout->image->col0 - (int)imageBounds->x0 || + y0 != readout->row0 + readout->image->row0 - (int)imageBounds->y0) { + psTrace(__func__, 5, "CELL.X0,Y0 don't match: %d,%d vs %d,%d\n", x0, y0, + readout->col0 + readout->image->col0 - (int)imageBounds->x0, + readout->row0 + readout->image->row0 - (int)imageBounds->y0); + psFree(imageBounds); + return NULL; + } + } + + // Need to check all the other chips if the HDU is in the FPA + pmFPA *fpa = chip->parent; // The parent FPA + if (fpa->hdu && fpa->hdu->images) { + psArray *chips = fpa->chips; // Array of chips + for (int i = 0; i < chips->n; i++) { + pmChip *testChip = chips->data[i]; // The chip of interest + if (testChip == chip) { + // Already done this one + continue; + } + if (!chipContiguous(imageBounds, testChip, false)) { + psTrace(__func__, 5, "Image isn't contiguous.\n"); + psFree(imageBounds); + return NULL; + } + } + } + + return imageBounds; +} + +// Mosaic multiple images, with flips, binning and offsets +static psImage *imageMosaic(const psArray *source, // Images to splice in + const psVector *xFlip, const psVector *yFlip, // Need to flip x and y? + const psVector *xBinSource, // Binning in x of source images + const psVector *yBinSource, // Binning in 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); + assert(xFlip->n == source->n); + assert(yFlip->n == source->n); + assert(xBinSource->n == source->n); + assert(yBinSource->n == source->n); + assert(x0->n == source->n); + assert(y0->n == source->n); + + if (source->n == 0) { + return NULL; + } + + // Get the maximum extent of the mosaic image + int xMin = INT_MAX; + int xMax = - INT_MAX; + int yMin = INT_MAX; + int yMax = - INT_MAX; + psElemType type = 0; + int numImages = 0; // Number of images + for (int i = 0; i < source->n; i++) { + psImage *image = source->data[i]; // The image of interest + if (!image) { + continue; + } + numImages++; + + // Only implemented for F32 and U8 images so far. + assert(image->type.type == PS_TYPE_F32 || image->type.type == PS_TYPE_U8); + // All input types must be the same + if (type == 0) { + type = image->type.type; + } + assert(type == image->type.type); + + // 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); + } + if (numImages == 0) { + return NULL; + } + + // 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, type); // 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 (!image) { + continue; + } + 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. + #define FILL_IN(TYPE) \ + 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.TYPE[yTarget][xTarget] += image->data.TYPE[y][x]; \ + } \ + } + + switch (type) { + case PS_TYPE_F32: + FILL_IN(F32); + break; + case PS_TYPE_U8: + FILL_IN(U8); + break; + default: + psAbort(__func__, "Should never get here.\n"); + } + + } + } // Iterating over input image + } + } + + return mosaic; +} + + +// Set the concepts in the new cell, based on the values in the old one +static bool cellConcepts(pmCell *target,// Target cell + psArray *sources, // Source cells + int xBin, int yBin, // Binning + psRegion *trimsec // The trim section + ) +{ + bool success = true; // Result of setting everything + float gain = 0.0; // Gain + float readnoise = 0.0; // Read noise + float saturation = INFINITY; // Saturation level + float bad = -INFINITY; // Bad level + float exposure = 0.0; // Exposure time + float darktime = 0.0; // Dark time + double time = 0.0; // Time of observation + psTimeType timeSys = 0; // Time system + + int nCells = 0; // Number of cells; + for (int i = 0; i < sources->n; i++) { + pmCell *cell = sources->data[i];// The cell of interest + if (!cell) { + continue; + } + + nCells++; + gain += psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); + readnoise += psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE"); + exposure += psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); + darktime += psMetadataLookupF32(NULL, cell->concepts, "CELL.DARKTIME"); + psTime *cellTime = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TIME"); + time += psTimeToMJD(cellTime); + if (i == 0) { + timeSys = psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS"); + } else if (timeSys != psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS")) { + psLogMsg(__func__, PS_LOG_ERROR, "Differing time systems in use: %d vs %d\n", timeSys, + psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS")); + success = false; + } + + float cellSaturation = psMetadataLookupF32(NULL, cell->concepts, "CELL.SATURATION"); + if (cellSaturation < saturation) { + saturation = cellSaturation; + } + float cellBad = psMetadataLookupF32(NULL, cell->concepts, "CELL.BAD"); + if (cellBad < bad) { + bad = cellBad; + } + } + gain /= (float)nCells; + readnoise /= (float)nCells; + exposure /= (float)nCells; + darktime /= (float)nCells; + time /= (double)nCells; + + MD_UPDATE(target->concepts, "CELL.GAIN", F32, gain); + MD_UPDATE(target->concepts, "CELL.READNOISE", F32, readnoise); + MD_UPDATE(target->concepts, "CELL.SATURATION", F32, saturation); + MD_UPDATE(target->concepts, "CELL.BAD", F32, bad); + MD_UPDATE(target->concepts, "CELL.EXPOSURE", F32, exposure); + MD_UPDATE(target->concepts, "CELL.DARKTIME", F32, darktime); + MD_UPDATE(target->concepts, "CELL.TIMESYS", S32, timeSys); + MD_UPDATE(target->concepts, "CELL.X0", S32, 0); + MD_UPDATE(target->concepts, "CELL.Y0", S32, 0); + MD_UPDATE(target->concepts, "CELL.XPARITY", S32, 1); + MD_UPDATE(target->concepts, "CELL.YPARITY", S32, 1); + MD_UPDATE(target->concepts, "CELL.XBIN", S32, xBin); + MD_UPDATE(target->concepts, "CELL.YBIN", S32, yBin); + + // CELL.TIME needs special care + { + psMetadataItem *timeItem = psMetadataLookup(target->concepts, "CELL.TIME"); + psFree(timeItem->data.V); + timeItem->data.V = psTimeFromMJD(time); + } + + // CELL.TRIMSEC needs special care + { + psMetadataItem *trimsecItem = psMetadataLookup(target->concepts, "CELL.TRIMSEC"); + psFree(trimsecItem->data.V); + trimsecItem->data.V = psMemIncrRefCounter(trimsec); + } + + + return success; +} + + +// Mosaic together the cells in a chip +bool chipMosaic(psImage **mosaicImage, // The mosaic image, to be returned + psImage **mosaicMask, // The mosaic mask, to be returned + psImage **mosaicWeights, // The mosaic weights, to be returned + int *xBinChip, int *yBinChip, // The binning in x and y, to be returned + pmChip *chip // Chip to mosaic + ) +{ + psArray *cells = chip->cells; // The array of cells + int numCells = cells->n; // Number of cells + psArray *images = psArrayAlloc(numCells); // Array of images that will be mosaicked + psArray *weights = psArrayAlloc(numCells); // Array of weight images to be mosaicked + psArray *masks = psArrayAlloc(numCells); // Array of mask images to be mosaicked + psVector *x0 = psVectorAlloc(numCells, PS_TYPE_S32); // Origin x coordinates + psVector *y0 = psVectorAlloc(numCells, PS_TYPE_S32); // Origin y coordinates + psVector *xBin = psVectorAlloc(numCells, PS_TYPE_S32); // Binning in x + psVector *yBin = psVectorAlloc(numCells, PS_TYPE_S32); // Binning in y + psVector *xFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Flip in x? + psVector *yFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Flip in y? + images->n = weights->n = masks->n = numCells; + x0->n = y0->n = xBin->n = yBin->n = xFlip->n = yFlip->n = numCells; + + // Binning for the mosaicked chip is the minimum binning allowed by the cells + *xBinChip = INT_MAX; + *yBinChip = INT_MAX; + + // Set up the required inputs + psTrace(__func__, 1, "Mosaicking %d cells...\n", numCells); + bool good = true; // Is everything good, well-behaved? + for (int i = 0; i < numCells && good; i++) { + pmCell *cell = cells->data[i]; // The cell of interest + if (!cell) { + continue; + } + bool mdok = true; // Status of MD lookup + x0->data.S32[i] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.X0"); + if (!mdok) { + psError(PS_ERR_IO, true, "CELL.X0 for cell %d is not set.\n", i); + good = false; + } + y0->data.S32[i] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.Y0"); + if (!mdok) { + psError(PS_ERR_IO, true, "CELL.Y0 for cell %d is not set.\n", i); + good = false; + } + psTrace(__func__, 5, "Cell %d: x0=%d y0=%d\n", i, x0->data.S32[i], y0->data.S32[i]); + xBin->data.S32[i] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XBIN"); + if (!mdok || xBin->data.S32[i] == 0) { + psError(PS_ERR_IO, true, "CELL.XBIN for cell %d is not set.\n", i); + good = false; + } + if (xBin->data.S32[i] < *xBinChip) { + *xBinChip = xBin->data.S32[i]; + } + yBin->data.S32[i] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YBIN"); + if (!mdok || yBin->data.S32[i] == 0) { + psError(PS_ERR_IO, true, "CELL.YBIN for cell %d is not set.\n", i); + good = false; + } + if (yBin->data.S32[i] < *yBinChip) { + *yBinChip = yBin->data.S32[i]; + } + int xParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XPARITY"); + if (!mdok || (xParity != 1 && xParity != -1)) { + psError(PS_ERR_IO, true, "CELL.XPARITY for cell %d is not set.\n", i); + good = false; + } + int yParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YPARITY"); + if (!mdok || (yParity != 1 && yParity != -1)) { + psError(PS_ERR_IO, true, "CELL.YPARITY for cell %d is not set.\n", i); + good = false; + } + if (xParity == -1) { + xFlip->data.U8[i] = 1; + } else { + xFlip->data.U8[i] = 0; + } + if (yParity == -1) { + yFlip->data.U8[i] = 1; + } else { + yFlip->data.U8[i] = 0; + } + + 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); + } + pmReadout *readout = readouts->data[0]; // The only readout we'll bother with + + // The images to put into the mosaic + images->data[i] = psMemIncrRefCounter(readout->image); + weights->data[i] = psMemIncrRefCounter(readout->weight); + masks->data[i] = psMemIncrRefCounter(readout->mask); + } + + // Mosaic the images together and we're done + if (good) { + *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0); + *mosaicWeights = imageMosaic(weights, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0); + *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0); + } + + // Clean up + psFree(images); + psFree(weights); + psFree(masks); + psFree(xFlip); + psFree(yFlip); + psFree(xBin); + psFree(yBin); + psFree(x0); + psFree(y0); + + return good; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Mosaic all the cells in a chip together. +// +// It is desirable to do this without using psImageOverlay (or similar) if it can be at all avoided (because +// it's really really slow in that case). There are therefore two cases: +// +// 1. The HDU is at the Chip or FPA level. This is the fast case, and only works if the HDU is "nice", by +// which I mean: +// +// - the CELL.TRIMSECs are contiguous on the HDU image +// - the CELL.PARITYs are identically +1 +// - the CELL.XBIN and CELL.YBIN are all identical +// +// Then we can just use psImageSubset to get the "mosaicked" chip. +// +// +// 2. The HDU is at the cell level, or the above requirements are not met, in which case we mosaic the cells. +// This is the slow case. We need to: +// +// - Throw away the bias regions +// - Convert all cells to common parity +// - Mosaic the cells into an HDU image using CELL.X0 and CELL.Y0 +// - Update CELL.TRIMSECs +// +// Once the demands of case 1 have been met, or case 2 has been performed, then we can create a cell to hold +// the mosaic image. +bool pmChipMosaic(pmChip *chip // Chip whose cells will be mosaicked + ) +{ + psImage *mosaicImage = NULL; // The mosaic image + psImage *mosaicMask = NULL; // The mosaic mask + psImage *mosaicWeights = NULL; // The mosaic weights + + // Find the HDU + psRegion *chipRegion = NULL; // Region on the HDU that corresponds to the chip + int xBin = 0, yBin = 0; // Binning for the chip mosaic + if ((chipRegion = niceChip(&xBin, &yBin, chip))) { + // Case 1 --- we need only cut out the region + psTrace(__func__, 1, "Case 1 mosaicking: simple cut-out.\n"); + pmHDU *hdu = chip->hdu; // The HDU that has the pixels + if (!hdu || !hdu->images) { + hdu = chip->parent->hdu; + } + mosaicImage = psMemIncrRefCounter(psImageSubset(hdu->images->data[0], *chipRegion)); + if (hdu->masks) { + mosaicMask = psMemIncrRefCounter(psImageSubset(hdu->masks->data[0], *chipRegion)); + } + if (hdu->weights) { + mosaicWeights = psMemIncrRefCounter(psImageSubset(hdu->weights->data[0], *chipRegion)); + } + } else { + // Case 2 --- we need to mosaic by cut and paste + psTrace(__func__, 1, "Case 2 mosaicking: cut and paste.\n"); + if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicWeights, &xBin, &yBin, chip)) { + psError(PS_ERR_IO, false, "Unable to mosaic cells.\n"); + return false; + } + chipRegion = psRegionAlloc(NAN, NAN, NAN, NAN); // We've cut and paste, so there's no valid trimsec + } + + // Construct a new cell, set the concepts, and add the mosaic in + pmCell *newCell = pmCellAlloc(NULL, "MOSAIC"); // New cell + cellConcepts(newCell, chip->cells, xBin, yBin, chipRegion); + psFree(chipRegion); + chip->mosaic = (struct pmCell*)newCell; + + // Now make a new readout to go in the new cell + pmReadout *newReadout = pmReadoutAlloc(newCell); // New readout + newReadout->image = mosaicImage; + newReadout->mask = mosaicMask; + newReadout->weight = mosaicWeights; + psFree(newReadout); // Drop reference + + return true; +} Index: trunk/psModules/src/camera/pmChipMosaic.h =================================================================== --- trunk/psModules/src/camera/pmChipMosaic.h (revision 7017) +++ trunk/psModules/src/camera/pmChipMosaic.h (revision 7017) @@ -0,0 +1,11 @@ +#ifndef PM_CHIP_MOSAIC_H +#define PM_CHIP_MOSAIC_H + +#include "pslib.h" +#include "pmFPA.h" + +// Mosaic all cells within a chip +bool pmChipMosaic(pmChip *chip // Chip whose cells will be mosaicked + ); + +#endif Index: trunk/psModules/src/camera/pmFPA.c =================================================================== --- trunk/psModules/src/camera/pmFPA.c (revision 7017) +++ trunk/psModules/src/camera/pmFPA.c (revision 7017) @@ -0,0 +1,575 @@ +/** @file pmFPA.c +* +* @brief This file defines the basic types for the FPA hierarchy +* +* @ingroup AstroImage +* +* @author GLG, MHPCC +* +* XXX: We should review the extent of the warning messages on these functions +* when the transformations are not successful. +* +* XXX: Should we implement non-linear cell->chip transforms? +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-05-01 01:55:43 $ +* +* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +*/ + +/******************************************************************************/ +/* INCLUDE FILES */ +/******************************************************************************/ +#include +#include +#include +#include "pslib.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmConcepts.h" +#include "pmMaskBadPixels.h" + +/****************************************************************************** + *****************************************************************************/ + +static void readoutFree(pmReadout *readout) +{ + // if this readout has a parent, drop that instance + if (readout->parent) { + psTrace(__func__, 9, "Removing readout %lx from cell %lx...\n", (size_t)readout, (size_t)readout->parent); + psArray *readouts = readout->parent->readouts; + for (int i = 0; i < readouts->n; i++) { + if (readouts->data[i] == readout) { + readouts->data[i] = NULL; + } + } + } + psTrace(__func__, 9, "Freeing readout %lx\n", (size_t) readout); + psFree(readout->image); + psFree(readout->mask); + psFree(readout->weight); + psFree(readout->analysis); + psFree(readout->bias); +} + +static void cellFree(pmCell *cell) +{ + + // if this cell has a parent, drop that instance + if (cell->parent) { + psTrace(__func__, 9, "Removing cell %lx from chip %lx...\n", (size_t)cell, (size_t)cell->parent); + psArray *cells = cell->parent->cells; + for (int i = 0; i < cells->n; i++) { + if (cells->data[i] == cell) { + cells->data[i] = NULL; + } + } + } + psTrace(__func__, 9, "Freeing cell %lx\n", (size_t)cell); + pmCellFreeReadouts(cell); + psFree(cell->readouts); + + psFree(cell->concepts); + psFree(cell->analysis); + psFree(cell->config); + psFree(cell->hdu); +} + +static void chipFree(pmChip* chip) +{ + // if this chip has a parent, drop that instance + if (chip->parent) { + psTrace(__func__, 9, "Removing chip %lx from fpa %lx...\n", (size_t)chip, (size_t)chip->parent); + psArray *chips = chip->parent->chips; + for (int i = 0; i < chips->n; i++) { + if (chips->data[i] == chip) { + chips->data[i] = NULL; + } + } + } + + psTrace(__func__, 9, "Freeing chip %lx\n", (size_t)chip); + pmChipFreeCells(chip); + psFree(chip->cells); + + psFree(chip->concepts); + psFree(chip->analysis); + psFree(chip->hdu); + psFree(chip->mosaic); + + # if FPA_ASTROM + + psFree(chip->toFPA); + psFree(chip->fromFPA); + # endif +} + + +static void FPAFree(pmFPA *fpa) +{ + psTrace(__func__, 9, "Freeing fpa %lx\n", (size_t)fpa); + + // NULL the parent pointers + psArray *chips = fpa->chips; + for (int i = 0 ; i < chips->n ; i++) { + pmChip *tmpChip = chips->data[i]; + if (! tmpChip) { + continue; + } + tmpChip->parent = NULL; + } + psFree(fpa->chips); + psFree(fpa->concepts); + psFree(fpa->analysis); + psFree(fpa->camera); + psFree(fpa->hdu); + + # if FPA_ASTROM + + psFree(fpa->fromTangentPlane); + psFree(fpa->toTangentPlane); + psFree(fpa->projection); + # endif +} + +void pmCellFreeReadouts(pmCell *cell) +{ + // + // Set the parent to NULL in all cell->readouts before psFree(cell->readouts) + // in order to avoid memory reference counter problems. + // + psArray *readouts = cell->readouts; + for (psS32 i = 0 ; i < readouts->n ; i++) { + pmReadout *tmpReadout = readouts->data[i]; + if (! tmpReadout) { + continue; + } + tmpReadout->parent = NULL; + psTrace(__func__, 9, "Will now free readout %lx...\n", (size_t)tmpReadout); + } + cell->readouts = psArrayRealloc(cell->readouts, 0); + cell->readouts->n = 0; +} + + +void pmChipFreeCells(pmChip *chip) +{ + // + // Set the parent to NULL in all chip->cells before psFree(chip->cells) + // in order to avoid memory reference counter problems. + // + psArray *cells = chip->cells; + for (int i = 0 ; i < cells->n ; i++) { + pmCell *tmpCell = cells->data[i]; + if (! tmpCell) { + continue; + } + tmpCell->parent = NULL; + pmCellFreeReadouts(tmpCell);// Drop all readouts the cell holds + } + chip->cells = psArrayRealloc(chip->cells, 0); + chip->cells->n = 0; +} + + +pmReadout *pmReadoutAlloc(pmCell *cell) +{ + pmReadout *tmpReadout = (pmReadout *) psAlloc(sizeof(pmReadout)); + psMemSetDeallocator(tmpReadout, (psFreeFunc) readoutFree); + + tmpReadout->image = NULL; + tmpReadout->mask = NULL; + tmpReadout->weight = NULL; + tmpReadout->bias = psListAlloc(NULL); + tmpReadout->analysis = psMetadataAlloc(); + tmpReadout->parent = cell; + if (cell != NULL) { + cell->readouts = psArrayAdd(cell->readouts, 1, (psPtr) tmpReadout); + } + + tmpReadout->process = true; // All cells are processed by default + tmpReadout->file_exists = false; // file not yet identified + tmpReadout->data_exists = false; // data yet read in + + tmpReadout->row0 = 0; + tmpReadout->col0 = 0; + + return(tmpReadout); +} + +pmCell *pmCellAlloc( + pmChip *chip, + const char *name +) +{ + pmCell *tmpCell = (pmCell *) psAlloc(sizeof(pmCell)); + psMemSetDeallocator(tmpCell, (psFreeFunc) cellFree); + + tmpCell->config = NULL; + tmpCell->analysis = psMetadataAlloc(); + tmpCell->readouts = psArrayAlloc(0); + tmpCell->parent = chip; + if (chip != NULL) { + chip->cells = psArrayAdd(chip->cells, 1, (psPtr) tmpCell); + } + tmpCell->hdu = NULL; + tmpCell->process = true; // All cells are processed by default + tmpCell->file_exists = false; // Not yet identified + tmpCell->data_exists = false; // Not yet read in + + tmpCell->concepts = psMetadataAlloc(); + tmpCell->conceptsRead = PM_CONCEPT_SOURCE_NONE; + if (!psMetadataAddStr(tmpCell->concepts, PS_LIST_HEAD, "CELL.NAME", 0, NULL, name)) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not add CELL.NAME to metadata.\n"); + } + pmConceptsBlankCell(tmpCell); + + return(tmpCell); +} + +pmChip *pmChipAlloc( + pmFPA *fpa, + const char *name) +{ + pmChip *tmpChip = (pmChip *) psAlloc(sizeof(pmChip)); + psMemSetDeallocator(tmpChip, (psFreeFunc) chipFree); + + # if FPA_ASTROM + + tmpChip->col0 = 0; +tmpChip->row0 = 0; +tmpChip->toFPA = NULL; +tmpChip->fromFPA = NULL; +# endif + +tmpChip->analysis = psMetadataAlloc(); + tmpChip->cells = psArrayAlloc(0); + tmpChip->parent = fpa; + if (fpa != NULL) { + psArrayAdd(fpa->chips, 1, (psPtr) tmpChip); + } + tmpChip->hdu = NULL; + tmpChip->mosaic = NULL; + tmpChip->process = true; // Work on all chips, by default + tmpChip->file_exists = false; // Not yet identified + tmpChip->data_exists = false; // Not yet read in + + tmpChip->concepts = psMetadataAlloc(); + tmpChip->conceptsRead = PM_CONCEPT_SOURCE_NONE; + if (!psMetadataAddStr(tmpChip->concepts, PS_LIST_HEAD, "CHIP.NAME", 0, NULL, name)) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not add CHIP.NAME %s to concepts.\n", name); + } + pmConceptsBlankChip(tmpChip); + return(tmpChip); +} + +pmFPA *pmFPAAlloc(const psMetadata *camera) +{ + pmFPA *tmpFPA = (pmFPA *) psAlloc(sizeof(pmFPA)); + psMemSetDeallocator(tmpFPA, (psFreeFunc) FPAFree); + + # if FPA_ASTROM + + tmpFPA->fromTangentPlane = NULL; +tmpFPA->toTangentPlane = NULL; +tmpFPA->projection = NULL; +# endif + +tmpFPA->analysis = NULL; +tmpFPA->camera = psMemIncrRefCounter((psPtr)camera); + tmpFPA->chips = psArrayAlloc(0); + tmpFPA->hdu = NULL; + + tmpFPA->concepts = psMetadataAlloc(); + tmpFPA->conceptsRead = PM_CONCEPT_SOURCE_NONE; + pmConceptsBlankFPA(tmpFPA); + + return(tmpFPA); +} + +static psBool cellCheckParents(pmCell *cell) +{ + if (!cell) { + return(true); + } + psBool flag = true; + + for (psS32 i = 0 ; i < cell->readouts->n ; i++) { + pmReadout *tmpReadout = (pmReadout *) cell->readouts->data[i]; + PS_ASSERT_PTR_NON_NULL(tmpReadout, false); + if (tmpReadout->parent != cell) { + tmpReadout->parent = cell; + flag = false; + } + } + return(flag); +} + +static psBool chipCheckParents(pmChip *chip) +{ + if (!chip) { + return(true); + } + psBool flag = true; + + for (psS32 i = 0 ; i < chip->cells->n ; i++) { + pmCell *tmpCell = (pmCell *) chip->cells->data[i]; + PS_ASSERT_PTR_NON_NULL(tmpCell, false); + if (tmpCell->parent != chip) { + tmpCell->parent = chip; + flag = false; + } + + flag &= cellCheckParents(tmpCell); + } + return(flag); +} + +psBool pmFPACheckParents(pmFPA *fpa) +{ + if (!fpa) { + return(true); + } + psBool flag = true; + + for (psS32 i = 0 ; i < fpa->chips->n ; i++) { + pmChip *tmpChip = (pmChip *) fpa->chips->data[i]; + PS_ASSERT_PTR_NON_NULL(tmpChip, false); + if (tmpChip->parent != fpa) { + tmpChip->parent = fpa; + flag = false; + } + + flag &= chipCheckParents(tmpChip); + } + return(flag); +} + +/** functions to turn on/off the file_exists flag **/ +bool pmFPASetFileStatus (pmFPA *fpa, bool status) +{ + + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + pmChipSetFileStatus (chip, status); + } + return true; +} + +bool pmChipSetFileStatus (pmChip *chip, bool status) +{ + + chip->file_exists = status; + for (int i = 0; i < chip->cells->n; i++) { + pmCell *cell = chip->cells->data[i]; + pmCellSetFileStatus (cell, status); + } + return true; +} + +bool pmCellSetFileStatus (pmCell *cell, bool status) +{ + + cell->file_exists = status; + for (int i = 0; i < cell->readouts->n; i++) { + pmReadout *readout = cell->readouts->data[i]; + readout->file_exists = status; + } + return true; +} + +/** functions to turn on/off the data_exists flag **/ +bool pmFPASetDataStatus (pmFPA *fpa, bool status) +{ + + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + pmChipSetDataStatus (chip, status); + } + return true; +} + +bool pmChipSetDataStatus (pmChip *chip, bool status) +{ + + chip->data_exists = status; + for (int i = 0; i < chip->cells->n; i++) { + pmCell *cell = chip->cells->data[i]; + pmCellSetDataStatus (cell, status); + } + return true; +} + +bool pmCellSetDataStatus (pmCell *cell, bool status) +{ + + cell->data_exists = status; + for (int i = 0; i < cell->readouts->n; i++) { + pmReadout *readout = cell->readouts->data[i]; + readout->data_exists = status; + } + return true; +} + +/***************************************************************************** + *****************************************************************************/ + +// Set cells within a chip to be processed or not +static bool setCellsProcess(const pmChip *chip, // Chip of interest + bool process // Process this chip? + ) +{ + PS_ASSERT_PTR_NON_NULL(chip, false); + + psArray *cells = chip->cells; // Component cells + if (! cells) { + return false; + } + for (int i = 0; i < cells->n; i++) { + pmCell *tmpCell = cells->data[i]; // Cell of interest + if (tmpCell) { + tmpCell->process = process; + } + } + + return true; +} + +/***************************************************************************** +XXX EAM : I've added the 'exclusive' option. if true all other chips are de-selected +XXX EAM : a negative value is valid and, in combinations with exclusive, de-selects all chips + *****************************************************************************/ +bool pmFPASelectChip(pmFPA *fpa, int chipNum, bool exclusive) +{ + PS_ASSERT_PTR_NON_NULL(fpa, false); + + psArray *chips = fpa->chips; // Component chips + if ((chips == NULL) || (chipNum >= chips->n)) { + return(false); + } + + for (int i = 0 ; i < chips->n ; i++) { + pmChip *tmpChip = (pmChip *) chips->data[i]; + if (tmpChip == NULL) { + continue; + } + if (i == chipNum) { + tmpChip->process = true; + setCellsProcess(tmpChip, true); + } else { + if (exclusive) { + tmpChip->process = false; + setCellsProcess(tmpChip, false); + } + } + + } + + return true; +} + +/***************************************************************************** +XXX EAM : I've added the 'exclusive' option. if true all other chips are de-selected +XXX EAM : a negative value is valid and, in combinations with exclusive, de-selects all cells +XXX this function should probably be re-defined to merge with 'setCellsProcess' + *****************************************************************************/ +bool pmChipSelectCell(pmChip *chip, int cellNum, bool exclusive) +{ + assert(chip); + + psArray *cells = chip->cells; // Component cells + if (!cells || cellNum > cells->n) { + return false; + } + + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; + if (!cell) { + continue; + } + if (i == cellNum) { + cell->process = true; + } else { + if (exclusive) { + cell->process = false; + } + } + } + return true; +} + +/***************************************************************************** +XXX: The SDRS is ambiguous on a few things: + Whether or not the other chips should be set process=true. [PAP: No] + Should we return the number of chip process=true before or after they're set, [PAP: After] + *****************************************************************************/ +/** + * + * pmFPAExcludeChip shall set process to false only for the specified chip + * number (chipNum). In the event that the specified chip number does not exist + * within the fpa, the function shall generate a warning, and perform no action. + * The function shall return the number of chips within the fpa that have process + * set to true. + * + */ +int pmFPAExcludeChip( + pmFPA *fpa, + int chipNum) +{ + PS_ASSERT_PTR_NON_NULL(fpa, false); + + psArray *chips = fpa->chips; // Component chips + if (chips == NULL) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: fpa->chips == NULL\n"); + return(0); + } + if ((chipNum >= chips->n) || (NULL == (pmChip *) chips->data[chipNum])) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: the specified chip (%d) does not exist.\n", chipNum); + return(0); + } + + int numChips = 0; // Number of chips to be processed + for (int i = 0 ; i < chips->n ; i++) { + pmChip *tmpChip = (pmChip *) chips->data[i]; // Chip of interest + if (tmpChip != NULL) { + if (i == chipNum) { + tmpChip->process = false; + setCellsProcess(tmpChip, false); // Wipe out the cell as well + } else if (tmpChip->process) { + numChips++; + } + } + } + + return(numChips); +} + +int pmChipExcludeCell(pmChip *chip, + int cellNum + ) +{ + assert(chip); + + psArray *cells = chip->cells; // The component cells + if (!cells || cellNum > cells->n) { + return 0; + } + + int numCells = 0; // Number of cells to be processed + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; + if (!cell) { + continue; + } + if (i == cellNum) { + cell->process = false; + } else { + numCells++; + } + } + + return numCells; +} + + Index: trunk/psModules/src/camera/pmFPA.h =================================================================== --- trunk/psModules/src/camera/pmFPA.h (revision 7017) +++ trunk/psModules/src/camera/pmFPA.h (revision 7017) @@ -0,0 +1,286 @@ +/** @file pmFPA.h +* +* @brief This file defines the basic types the focal plane hierarchy. +* +* @ingroup AstroImage +* +* @author GLG, MHPCC +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-05-01 01:55:43 $ +* +* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +*/ + +#ifndef PM_FPA_H +#define PM_FPA_H +#if HAVE_CONFIG_H +#include +#endif + +#include "pslib.h" +#include "pmHDU.h" + +# define FPA_ASTROM 1 + +/// @addtogroup AstroImage +/// @{ + +/** Focal plane data structure + * + * A focal plane consists of one or more chips (according to the number of + * pieces of contiguous silicon). It contains metadata containers for the + * concepts and analysis, a link to the parent, and pointers to the FITS header, + * if that corresponds to this level (the FPA may be the PHU, but will not ever + * contain pixels). For astrometry, it contains a transformation from the focal + * plane to the tangent plane and the fixed pattern residuals. It is expected + * that the transformation will consist of two 4D polynomials (i.e. a function + * of two coordinates in position, the magnitude of the object, and the color of + * the object) in order to correct for optical distortions and the effects of + * the atmosphere; hence we think that it is prudent to include a reverse + * transformation which will be derived from numerically inverting the forward + * transformation. + * + */ +typedef struct +{ + # if FPA_ASTROM + // 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 +# endif +// Information +psMetadata *concepts; ///< Cache for PS concepts +unsigned int conceptsRead; ///< Which concepts have been read +psMetadata *analysis; ///< FPA-level analysis metadata +const psMetadata *camera; ///< Camera configuration +psArray *chips; ///< The chips +pmHDU *hdu; ///< FITS data +} +pmFPA; + +/** Chip data structure + * + * A chip consists of one or more cells (according to the number of amplifiers + * on the device). The chip contains metadata containers for the concepts and + * analysis, a link to the parent, and pointers to the pointers to the various + * FITS data, if that corresponds to this level. For astrometry, in addition to + * the rough positioning information, it contains a coordinate transform from + * the chip to the focal plane. It is expected that this transform will consist + * of two second-order 2D polynomials; hence we think that it is prudent to + * include a reverse transformation which will be derived from numerically + * inverting the forward transformation. A boolean indicates whether the chip is + * of interest, allowing it to be excluded from analysis. + * + */ +typedef struct +{ +# if FPA_ASTROM + // 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 +# endif +// Information +psMetadata *concepts; ///< Cache for PS concepts +unsigned int conceptsRead; ///< Which concepts have been read +psMetadata *analysis; ///< Chip-level analysis metadata +psArray *cells; ///< The cells (referred to by name) +pmFPA *parent; ///< Parent FPA +bool process; ///< Do we bother about reading and working with this chip? +bool file_exists; ///< Does the file for this chip exist (read case only)? +bool data_exists; ///< Does the data for this chip exist (read case only)? +pmHDU *hdu; ///< FITS data +struct pmCell *mosaic; ///< A mosaic cell +} +pmChip; + +/** Cell data structure + * + * A cell consists of one or more readouts. It also contains a pointer to the + * cell's metadata, and its parent chip. On the astrometry side, it also + * contains coordinate transforms from the cell to chip, from the cell to + * focal-plane, as well as a "quick and dirty" tranform from the cell to + * sky coordinates. + * + */ +typedef struct +{ +psMetadata *concepts; ///< Cache for PS concepts +unsigned int conceptsRead; ///< Which concepts have been read +psMetadata *config; ///< Cell configuration info +psMetadata *analysis; ///< Cell-level analysis metadata +psArray *readouts; ///< The readouts (referred to by number) +pmChip *parent; ///< Parent chip +bool process; ///< Do we bother about reading and working with this cell? +bool file_exists; ///< Does the file for this cell exist (read case only)? +bool data_exists; ///< Does the data for this cell exist (read case only)? +pmHDU *hdu; ///< FITS data +} +pmCell; + +/** Readout data structure. + * + * A readout is the result of a single read of a cell (or a portion thereof). + * It contains the offset from the lower-left corner of the chip, in the case + * that the CCD was windowed, as well as the binning factors and parity (if the + * binning value is negative, then the parity is reversed). It also contains the + * pixel data, metadata containers for the concepts and analysis, and a link to + * the parent. + * + */ +typedef struct +{ +int col0; ///< Column offset; non-zero if reading in columns bit by bit +int row0; ///< Row offset; non-zero if reading in rows bit by bit +psImage *image; ///< Imaging area of readout +psImage *mask; ///< Mask of input image +psImage *weight; ///< Weight of input image +psList *bias; ///< Overscan images +psMetadata *analysis; ///< Readout-level analysis metadata +pmCell *parent; ///< Parent cell +bool process; ///< Do we bother about reading and working with this readout? +bool file_exists; ///< Does the file for this readout exist (read case only)? +bool data_exists; ///< Does the data for this readout exist (read case only)? +} +pmReadout; + +void pmCellFreeReadouts(pmCell *cell); +void pmChipFreeCells(pmChip *chip); + +/** Allocates a pmReadout + * + * The constructor shall make an empty pmReadout. If the parent cell is not + * NULL, the parent link is made and the readout shall be placed in the + * parents array of readouts. The metadata containers shall be allocated. All + * other pointers in the structure shall be initialized to NULL. + * + * @return pmReadout* newly allocated pmReadout with all internal pointers set to NULL + */ +pmReadout *pmReadoutAlloc( + pmCell *cell ///< Parent cell +); + +/** Allocates a pmCell + * + * The constructor shall make an empty pmCell. If the parent chip is not NULL, + * the parent link is made and the cell shall be placed in the parents array of + * cells. The readouts array shall be allocated with a zero size, and the + * metadata containers constructed. All other pointers in the structure shall be + * initialized to NULL. + * + * @return pmCell* newly allocated pmCell + */ +pmCell *pmCellAlloc( + pmChip *chip, ///< Parent chip + const char *name ///< Name of cell +); + +/** Allocates a pmChip + * + * The constructor shall make an empty pmChip. If the parent fpa is not NULL, + * the parent link is made and the chip shall be placed in the parent's array + * of chips. The cells array shall be allocated with a zero size, and the + * metadata containers constructed. All other pointers in the structure shall be + * initialized to NULL. + * + * @return pmChip* newly allocated pmChip + */ +pmChip *pmChipAlloc( + pmFPA *fpa, ///< FPA to which the chip belongs + const char *name ///< Name of chip +); + +/** Allocates a pmFPA + * + * The constructor shall make an empty pmFPA. The chips array shall be + * allocated with a zero size, the camera and db pointers set to the values + * provided, and the concepts metadata constructed. All other pointers in the + * structure shall be initialized to NULL. + * + */ +pmFPA *pmFPAAlloc( + const psMetadata *camera ///< Camera configuration +); + + +/** Verify parent links. + * + * This function checks the validity of the parent links in the FPA hierarchy. + * If a parent link is not set (or not set correctly), it is corrected, and the + * function shall return false. If all the parent pointers were correct, the + * function shall return true. + * + */ +bool pmFPACheckParents( + pmFPA *fpa +); + +/* functions to turn on/off the file_exists flags */ +bool pmFPASetFileStatus (pmFPA *fpa, bool status); +bool pmChipSetFileStatus (pmChip *chip, bool status); +bool pmCellSetFileStatus (pmCell *cell, bool status); + +/* functions to turn on/off the data_exists flags */ +bool pmFPASetDataStatus (pmFPA *fpa, bool status); +bool pmChipSetDataStatus (pmChip *chip, bool status); +bool pmCellSetDataStatus (pmCell *cell, bool status); + +/** Specify the level for an operation. + */ +typedef enum { + PM_FPA_LEVEL_NONE, ///< No particular level specified + PM_FPA_LEVEL_FPA, ///< Level corresponds to an FPA + PM_FPA_LEVEL_CHIP, ///< Level corresponds to a Chip + PM_FPA_LEVEL_CELL, ///< Level corresponds to a Cell + PM_FPA_LEVEL_READOUT ///< Level corresponds to a Readout +} pmFPALevel; + + +/** + * + * pmFPASelectChip shall set valid to true for the specified chip number + * (chipNum), and all other chips shall have valid set to false. In the event + * that the specified chip number does not exist within the fpa, the function + * shall return false. + * + */ +bool pmFPASelectChip( + pmFPA *fpa, + int chipNum, + bool exclusive +); + +bool pmChipSelectCell(pmChip *chip, + int cellNum, + bool exclusive + ); + +/** + * + * pmFPAExcludeChip shall set valid to false only for the specified chip + * number (chipNum). In the event that the specified chip number does not exist + * within the fpa, the function shall generate a warning, and perform no action. + * The function shall return the number of chips within the fpa that have valid + * set to true. + * + */ +int pmFPAExcludeChip( + pmFPA *fpa, + int chipNum +); + +int pmChipExcludeCell(pmChip *chip, + int cellNum + ); + +// Set the weights and masks within a cell, based on the gain and RN +bool pmCellSetWeights(pmCell *cell // Cell for which to set weights + ); + + + +#endif // #ifndef PM_FPA_H Index: trunk/psModules/src/camera/pmFPAConstruct.c =================================================================== --- trunk/psModules/src/camera/pmFPAConstruct.c (revision 7017) +++ trunk/psModules/src/camera/pmFPAConstruct.c (revision 7017) @@ -0,0 +1,781 @@ +#include +#include +#include +#include "pslib.h" +#include "psMetadataItemParse.h" + +#include "pmFPA.h" +#include "pmConcepts.h" +#include "pmFPAConstruct.h" +#include "pmFPAview.h" +#include "pmFPAUtils.h" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// File-static functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Read data for a particular cell from the camera format description +static psMetadata *getCellData(const psMetadata *format, // The camera format description + const char *cellName // The name of the cell + ) +{ + bool status = true; // Result of MD lookup + psMetadata *cells = psMetadataLookupMD(&status, format, "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 || !cellData) { + psLogMsg(__func__, PS_LOG_WARN, "Unable to find specs for cell %s: ignored\n", cellName); + } + + return cellData; +} + +// Parse a list of first:second:third pairs in a string +static bool parseContent(psArray **first, // Array of the first values + psArray **second, // Array of the second values + psArray **third, // Array of the third values + const char *string // The string to parse + ) +{ + bool allOK = true; // Everything was OK? + psList *values = psStringSplit(string, " ,;", true); // List of the parts + *first = psArrayAlloc(values->n); + *second = psArrayAlloc(values->n); + *third = psArrayAlloc(values->n); + (*first)->n = (*second)->n = (*third)->n = values->n; + int num = 0; + psListIterator *valuesIter = psListIteratorAlloc(values, PS_LIST_HEAD, false); // Iterator for values + psString value = NULL; // "first:second:third" string + while ((value = psListGetAndIncrement(valuesIter))) { + psList *firstSecondThird = psStringSplit(value, ":", true); // List containing the first, second, third + psArray *fst = psListToArray(firstSecondThird); // An array representation + psFree(firstSecondThird); + psString firstPart = NULL; // The first part + psString secondPart = NULL; // The second part + psString thirdPart = NULL; // The third part + switch (firstSecondThird->n) { + case 1: + thirdPart = fst->data[0]; + break; + case 2: + secondPart = fst->data[0]; + thirdPart = fst->data[1]; + break; + case 3: + firstPart = fst->data[0]; + secondPart = fst->data[1]; + thirdPart = fst->data[2]; + break; + default: + psLogMsg(__func__, PS_LOG_WARN, "Badly formated specifier: %s --- ignored.\n", value); + allOK = false; + psFree(fst); + continue; + } + psArraySet(*first, num, firstPart); + psArraySet(*second, num, secondPart); + psArraySet(*third, num, thirdPart); + num++; + psFree(fst); + } + psFree(valuesIter); + psFree(values); + + return allOK; +} + +// Get the name of a PHU chip or cell from the header +static psString phuNameFromHeader(const char *name, // The name to lookup: "CELL.NAME" or "CHIP.NAME" + const psMetadata *fileInfo, // FILE within the camera format description + const psMetadata *header // Primary header + ) +{ + bool mdok = true; // Result of MD lookup + psString keyword = psMetadataLookupStr(&mdok, fileInfo, name); + if (!mdok || strlen(keyword) == 0) { + return false; + } + psMetadataItem *resultItem = psMetadataLookup(header, keyword); + if (! resultItem) { + psError(PS_ERR_IO, false, "Unable to find %s in primary header to identify %s.\n", keyword, name); + return NULL; + } + return psMetadataItemParseString(resultItem); +} + +// Add an HDU to the FPA +static bool addHDUtoFPA(pmFPA *fpa, // FPA to which to add + pmHDU *hdu // HDU to be added + ) +{ + assert(fpa); + assert(hdu); + + if (fpa->hdu) { + // Something's already here + if (fpa->hdu != hdu) { + psError(PS_ERR_IO, true, "Unable to add HDU since FPA already has one.\n"); + } + return false; + } + fpa->hdu = psMemIncrRefCounter(hdu); + if (hdu->header) { + pmConceptsReadFPA(fpa, PM_CONCEPT_SOURCE_HEADER, NULL); + } + pmFPASetFileStatus(fpa, true); + + return true; +} + +// Add an HDU to the chip +static bool addHDUtoChip(pmChip *chip, // Chip to which to add + pmHDU *hdu // HDU to be added + ) +{ + assert(chip); + assert(hdu); + + if (chip->hdu) { + // Something's already here + if (chip->hdu != hdu) { + psError(PS_ERR_IO, true, "Unable to add HDU since chip already has one.\n"); + } + return false; + } + chip->hdu = psMemIncrRefCounter(hdu); + if (hdu->header) { + pmConceptsReadChip(chip, PM_CONCEPT_SOURCE_HEADER, true, NULL); + } + pmChipSetFileStatus(chip, true); + + return true; +} + +// Add an HDU to the cell +static bool addHDUtoCell(pmCell *cell, // Cell to which to add + pmHDU *hdu // HDU to be added + ) +{ + assert(cell); + assert(hdu); + + if (cell->hdu) { + // Something's already here + if (cell->hdu != hdu) { + psError(PS_ERR_IO, true, "Unable to add HDU since cell already has one.\n"); + } + return false; + } + cell->hdu = psMemIncrRefCounter(hdu); + if (hdu->header) { + pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_HEADER, true, NULL); + } + pmCellSetFileStatus(cell, true); + + return true; +} + + +// Looks up the particular content based on the header +static const char *getContent(psMetadata *fileInfo, // The FILE from the camera format configuration + psMetadata *contents, // The CONTENTS from the camera format configuration + psMetadata *header // The (primary) header + ) +{ + bool mdok = true; // Status of MD lookup + const char *contentHeaders = psMetadataLookupStr(&mdok, fileInfo, "CONTENT"); // Headers for content + if (!mdok || !contentHeaders || strlen(contentHeaders) == 0) { + psError(PS_ERR_IO, true, "Unable to find CONTENT in FILE within camera format configuration.\n"); + return NULL; + } + + psList *keywords = psStringSplit(contentHeaders, " ,;", true); // List of keywords + psListIterator *keywordsIter = psListIteratorAlloc(keywords, PS_LIST_HEAD, false); // Iterator + psString keyword = NULL; // Keyword, from iteration + psString contentsKey = NULL; // Key to the CONTENTS menu + bool first = true; // Is it the first value? + while ((keyword = psListGetAndIncrement(keywordsIter))) { + const char *value = psMetadataLookupStr(&mdok, header, keyword); + if (first) { + psStringAppend(&contentsKey, "%s", value); + first = false; + } else { + psStringAppend(&contentsKey, "_%s", value); + } + } + psFree(keywordsIter); + psFree(keywords); + psTrace(__func__, 5, "Looking up %s in the CONTENTS.\n", contentsKey); + const char *content = psMetadataLookupStr(&mdok, contents, contentsKey); + if (!mdok || !content || strlen(content) == 0) { + psError(PS_ERR_IO, true, "Unable to find %s in the CONTENTS.\n", contentsKey); + return NULL; + } + psFree(contentsKey); + + return content; +} + + +// Given a (string) list of contents "chip:cell:type chip:cell:type", put the HDU in the correct place and +// plug in the cell configuration information +static int processContents(pmFPA *fpa, // The FPA + pmChip *chip, // The chip, or NULL + pmCell *cell, // The cell, or NULL + pmHDU *hdu, // The HDU to be added + pmFPALevel level, // The level at which to add the HDU + const char *contents, // The contents line, consisting of a list of chip:cell + psMetadata *format // Camera format configuration + ) +{ + assert(fpa); + assert(contents && strlen(contents) > 0); + assert(!cell || (cell && chip)); // Need both chip and cell if given a cell + + if (hdu && level == PM_FPA_LEVEL_FPA) { + addHDUtoFPA(fpa, hdu); + } + + // Parse the list of chip:cell:type + psArray *chips = NULL; // The chips (first bits) + psArray *cells = NULL; // The cells (second bits) + psArray *types = NULL; // The cell types (third bits) + int numCells = 0; // Number of cells processed + parseContent(&chips, &cells, &types, contents); + for (int i = 0; i < types->n; i++) { + psString chipName = chips->data[i]; // The name of the chip + psString cellName = cells->data[i]; // The name of the cell + psString cellType = types->data[i]; // The type of the cell + + // Get the chip + pmChip *newChip = NULL; // The chip specified + if (chip) { + newChip = chip; + } else if (chipName) { + // Find the chip + int chipNum = pmFPAFindChip(fpa, chipName); // The chip we're looking for + if (chipNum == -1) { + psLogMsg(__func__, PS_LOG_WARN, "Unable to find chip %s in fpa --- ignored.\n", chipName); + continue; + } + newChip = fpa->chips->data[chipNum]; + } + + if (!newChip) { + psLogMsg(__func__, PS_LOG_WARN, "Unable to determine chip for entry %d of content (follows) --- " + "ignored.\n\t%s", i, contents); + continue; + } + + // Get the cell + pmCell *newCell = NULL; // The cell specified + if (cell) { + newCell = cell; + } else if (cellName) { + // Find the cell + int cellNum = pmChipFindCell(newChip, cellName); // The cell we're looking for + if (cellNum == -1) { + psLogMsg(__func__, PS_LOG_WARN, "Unable to find cell %s in chip %s --- ignored.\n", + cellName, chipName); + continue; + } + newCell = newChip->cells->data[cellNum]; + } + if (!newCell) { + psLogMsg(__func__, PS_LOG_WARN, "Unable to determine cell for entry %d of content (follows) --- " + "ignored.\n\t%s", i, contents); + continue; + } + + // Get the type + if (!cellType) { + psLogMsg(__func__, PS_LOG_WARN, "Unable to determine cell type for entry %d of content (follows) " + "--- ignored.\n\t%s", i, contents); + continue; + } + psMetadata *cellData = getCellData(format, cellType); // Data for this cell + + // Put in the HDU + if (hdu && level == PM_FPA_LEVEL_CHIP) { + addHDUtoChip(newChip, hdu); + } + if (hdu && level == PM_FPA_LEVEL_CELL) { + addHDUtoCell(newCell, hdu); + } + + // Put in the cell data + if (newCell->config) { + psLogMsg(__func__, PS_LOG_WARN, "Overwriting cell data in chip %s, cell %s\n", chipName, + cellName); + psFree(newCell->config); // Make way! + } + newCell->config = psMemIncrRefCounter(cellData); + pmConceptsReadCell(newCell, PM_CONCEPT_SOURCE_CAMERA | PM_CONCEPT_SOURCE_DEFAULTS, false, NULL); + numCells++; + } + psFree(chips); + psFree(cells); + psFree(types); + + return numCells; +} + +static pmFPALevel hduLevel(psMetadata *format // The camera format configuration + ) +{ + bool mdok = true; // Status of MD lookup + psMetadata *file = psMetadataLookupMD(&mdok, format, "FILE"); // File information + if (!mdok || !file) { + psError(PS_ERR_IO, true, "Unable to find FILE information in camera format configuration.\n"); + return PM_FPA_LEVEL_NONE; + } + const char *extType = psMetadataLookupStr(&mdok, file, "EXTENSIONS"); + if (!mdok || !extType || strlen(extType) == 0) { + psError(PS_ERR_IO, true, "Unable to find EXTENSIONS in the FILE information in the camera format" + " configuration.\n"); + return PM_FPA_LEVEL_NONE; + } + + // Where do we stick in the HDUs? + pmFPALevel level = PM_FPA_LEVEL_NONE; // Level for HDU insertion + if (strcasecmp(extType, "CHIP") == 0) { + level = PM_FPA_LEVEL_CHIP; + } else if (strcasecmp(extType, "CELL") == 0) { + level = PM_FPA_LEVEL_CELL; + } else if (strcasecmp(extType, "NONE") != 0) { + psError(PS_ERR_IO, true, "EXTENSIONS is not CHIP or CELL or NONE.\n"); + } + + return level; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Construct an FPA instance on the basis of a camera configuration +pmFPA *pmFPAConstruct(const psMetadata *camera // The camera configuration + ) +{ + pmFPA *fpa = pmFPAAlloc(camera); // The FPA to fill out + + bool mdok = true; // Status from MD lookups + psMetadata *components = psMetadataLookupMD(&mdok, camera, "FPA"); // FPA components + psMetadataIterator *componentsIter = psMetadataIteratorAlloc(components, PS_LIST_HEAD, NULL); + psMetadataItem *componentsItem = NULL; // Item from components + while ((componentsItem = psMetadataGetAndIncrement(componentsIter))) { + const char *chipName = componentsItem->name; // Name of the chip + if (componentsItem->type != PS_DATA_STRING) { + psLogMsg(__func__, PS_LOG_WARN, "Element %s in FPA within the camera configuration is not of " + "type STR (type=%x) --- ignored.\n", chipName, componentsItem->type); + continue; + } + + pmChip *chip = pmChipAlloc(fpa, chipName); // The chip + psList *cellNames = psStringSplit(componentsItem->data.V, " ,;", true); // List of cell names + psListIterator *cellNamesIter = psListIteratorAlloc(cellNames, PS_LIST_HEAD, false); // Iterator + + psString cellName = NULL; // Name of cell + while ((cellName = psListGetAndIncrement(cellNamesIter))) { + pmCell *cell = pmCellAlloc(chip, cellName); // New cell + psFree(cell); // Drop reference + } + psFree(chip); // Drop reference + psFree(cellNamesIter); + psFree(cellNames); + } + psFree(componentsIter); + + return fpa; +} + +bool pmFPAAddSourceFromView(pmFPA *fpa, // The FPA + const pmFPAview *phuView, // The view, corresponding to the PHU + psMetadata *format // Format of file + ) +{ + assert(fpa); + assert(phuView); + assert(format); + + // Where does the PHU go? + bool mdok = true; // Status of MD lookup + psMetadata *fileInfo = psMetadataLookupMD(&mdok, format, "FILE"); // File information from the format + if (!mdok || !fileInfo) { + psError(PS_ERR_IO, false, "Unable to find FILE information in the camera format configuration.\n"); + return false; + } + const char *phuType = psMetadataLookupStr(&mdok, fileInfo, "PHU"); // What is the PHU? + if (!mdok || strlen(phuType) == 0) { + psError(PS_ERR_IO, false, "Unable to find PHU in the FILE information of the camera format.\n"); + return false; + } + + // Generate the PHU + pmHDU *phdu = pmHDUAlloc("PHU"); // The primary header data unit + phdu->format = psMemIncrRefCounter(format); + + // Put in the PHU + pmChip *chip = NULL; // The chip that corresponds to the PHU + pmCell *cell = NULL; // The cell that corresponds to the PHU + if (strcasecmp(phuType, "FPA") == 0) { + addHDUtoFPA(fpa, phdu); + psFree(phdu); + } else { + psArray *chips = fpa->chips; // Array of chips + if (phuView->chip < 0 || phuView->chip > chips->n) { + psError(PS_ERR_IO, true, "PHU specified by the camera format requires specification of a " + "particular chip, which cannot be determined from the view (%d).\n", phuView->chip); + psFree(phdu); + return false; + } + chip = chips->data[phuView->chip]; + if (strcasecmp(phuType, "CHIP") == 0) { + addHDUtoChip(chip, phdu); + psFree(phdu); + } else if (strcasecmp(phuType, "CELL") == 0) { + psArray *cells = chip->cells; // Array of cells + if (phuView->cell < 0 || phuView->cell > cells->n) { + psError(PS_ERR_IO, true, "PHU specified by the camera format requires specification of a " + "particular cell, which cannot be determined from the view (%d).\n", phuView->cell); + psFree(phdu); + return false; + } + addHDUtoCell(cell, phdu); + psFree(phdu); + } else { + psError(PS_ERR_IO, true, "PHU in the camera configuration format is not FPA, CHIP or CELL.\n"); + psFree(phdu); + return false; + } + } + + // Put in the extensions + pmFPALevel level = hduLevel(format);// The level for the HDUs to go + if (level == PM_FPA_LEVEL_NONE) { + // No extensions --- we're done + return true; + } + psMetadata *contents = psMetadataLookupMD(&mdok, format, "CONTENTS"); // The contents of the FITS file + if (!mdok || !contents) { + psError(PS_ERR_IO, false, "Unable to find CONTENTS in the camera format configuration.\n"); + return false; + } + + psMetadataIterator *contentsIter = psMetadataIteratorAlloc(contents, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *contentsItem = NULL;// Item from the contents iteration + while ((contentsItem = psMetadataGetAndIncrement(contentsIter))) { + if (contentsItem->type != PS_DATA_STRING) { + psLogMsg(__func__, PS_LOG_WARN, "Content item %s is not of type STR --- ignored.\n", + contentsItem->name); + continue; + } + pmHDU *hdu = pmHDUAlloc(contentsItem->name); // HDU to add + hdu->format = psMemIncrRefCounter(format); + const char *content = contentsItem->data.V; // The content data + processContents(fpa, chip, cell, hdu, level, content, format); + psFree(hdu); + } + psFree(contentsIter); + + return true; +} + + + +// Add an input file to the FPA +pmFPAview *pmFPAAddSourceFromHeader(pmFPA *fpa, // The FPA + psMetadata *phu, // Primary header of file + psMetadata *format // Format of file + ) +{ + assert(fpa); + assert(phu); + assert(format); + + bool mdok = true; // Status from metadata lookups + psMetadata *fileInfo = psMetadataLookupMD(&mdok, format, "FILE"); // The file information + if (!mdok || !fileInfo) { + psError(PS_ERR_IO, false, "Unable to find FILE in the camera format configuration.\n"); + return NULL; + } + + // Check the name of the FPA + psString newFPAname = phuNameFromHeader("FPA.NAME", fileInfo, phu); // New name for the FPA + const char *currentFPAname = psMetadataLookupStr(&mdok, fpa->concepts, "FPA.NAME"); // Current name + if (mdok && currentFPAname && strlen(currentFPAname) > 0 && strcmp(currentFPAname, newFPAname) != 0) { + psLogMsg(__func__, PS_LOG_WARN, "FPA.NAME for new source (%s) doesn't match FPA.NAME for current " + "fpa (%s).\n", newFPAname, currentFPAname); + } + psMetadataAddStr(fpa->concepts, PS_LIST_HEAD, "FPA.NAME", PS_META_REPLACE, "Name of FPA", newFPAname); + psFree(newFPAname); // Drop reference + + // Where does the PHU go? + const char *phuType = psMetadataLookupStr(&mdok, fileInfo, "PHU"); // What is the PHU? + if (!mdok || strlen(phuType) == 0) { + psError(PS_ERR_IO, false, "Unable to find PHU in the format specification.\n"); + return NULL; + } + pmHDU *phdu = pmHDUAlloc("PHU"); // The primary header data unit + phdu->header = psMemIncrRefCounter(phu); + phdu->format = psMemIncrRefCounter(format); + // Generate the view + pmFPAview *view = pmFPAviewAlloc(0); // The FPA view corresponding to the PHU, to be returned + view->chip = -1; + view->cell = -1; + view->readout = -1; + + // And what are the individual extensions? + const char *extType = psMetadataLookupStr(&mdok, fileInfo, "EXTENSIONS"); // What's in the extns? + if (!mdok || strlen(extType) == 0) { + psError(PS_ERR_IO, false, "Unable to find EXTENSIONS in the format specification.\n"); + psFree(view); + return NULL; + } + + + // This is a special case: PHU=FPA and EXTENSIONS=NONE. The only reason to do this is in the case of + // single CCD imagers, where the entire FPA is in the PHU image. In this case, the CONTENTS is of type + // STR (rather than METADATA), and has the usual chip:cell. + if (strcasecmp(phuType, "FPA") == 0 && strcasecmp(extType, "NONE") == 0) { + const char *contents = psMetadataLookupStr(&mdok, format, "CONTENTS"); // The contents of the file + if (!mdok || !contents || strlen(contents) == 0) { + psError(PS_ERR_IO, true, "Unable to find CONTENTS in the camera format configuration.\n"); + psFree(view); + return NULL; + } + + processContents(fpa, NULL, NULL, phdu, PM_FPA_LEVEL_FPA, contents, format); + psFree(phdu); + return view; + } + + // In all other cases, the CONTENTS is of type METADATA, and is either a menu if EXTENSIONS=NONE, or + // a list of extensions otherwise. + psMetadata *contents = psMetadataLookupMD(&mdok, format, "CONTENTS"); // The contents of the FITS file + if (!mdok || !contents) { + psError(PS_ERR_IO, false, "Unable to find CONTENTS in the camera format configuration.\n"); + psFree(view); + return NULL; + } + + // No extensions --- only have the PHU. In this case, the CONTENTS is a menu, with CONTENT indicating + // header keywords that provides a key to the CONTENTS menu. + if (strcasecmp(extType, "NONE") == 0) { + // We have already dealt with the case PHU=FPA, in a special case, above. + const char *content = getContent(fileInfo, contents, phu); // The content: string of chip:cell pairs + + // Need to look up what chip we have. + psString chipName = phuNameFromHeader("CHIP.NAME", fileInfo, phu); + psTrace(__func__, 5, "This is chip %s\n", chipName); + int chipNum = pmFPAFindChip(fpa, chipName); // Chip number + if (chipNum == -1) { + psError(PS_ERR_IO, true, "Unable to find chip %s in FPA.\n", chipName); + psFree(view); + return NULL; + } + pmChip *chip = fpa->chips->data[chipNum]; // Chip of interest + view->chip = chipNum; + + pmCell *cell = NULL; // Cell of interest + + pmFPALevel level = PM_FPA_LEVEL_NONE; // Level for HDU to be added + if (strcasecmp(phuType, "CHIP") == 0) { + level = PM_FPA_LEVEL_CHIP; + } else if (strcasecmp(phuType, "CELL") == 0) { + level = PM_FPA_LEVEL_CELL; + // Need to look up what cell we have. + psString cellName = phuNameFromHeader("CELL.NAME", fileInfo, phu); + int cellNum = pmChipFindCell(chip, cellName); // Cell number + if (cellNum == -1) { + psError(PS_ERR_IO, true, "Unable to find cell %s in chip %s.\n", cellName, chipName); + psFree(view); + return NULL; + } + cell = chip->cells->data[cellNum]; + view->cell = cellNum; + psFree(cellName); + } else { + psError(PS_ERR_IO, true, "PHU is not FPA, CHIP or CELL.\n"); + psFree(view); + return NULL; + } + psFree(chipName); + + processContents(fpa, chip, cell, phdu, level, content, format); + psFree(phdu); + return view; + } + + + // From here on, we have extensions that we iterate through. The CONTENTS is a list of extensions. + pmChip *chip = NULL; // The chip of interest + pmCell *cell = NULL; // The cell of interest + + // First, put in the PHU + if (strcasecmp(phuType, "FPA") == 0) { + addHDUtoFPA(fpa, phdu); + } else { + // Get the chip + psString chipName = phuNameFromHeader("CHIP.NAME", fileInfo, phu); // Name of the chip + int chipNum = pmFPAFindChip(fpa, chipName); // Chip number + if (chipNum == -1) { + psError(PS_ERR_IO, true, "Unable to find chip %s in FPA.\n", chipName); + psFree(view); + return NULL; + } + chip = fpa->chips->data[chipNum]; // The specified chip + view->chip = chipNum; + + if (strcasecmp(phuType, "CHIP") == 0) { + addHDUtoChip(chip, phdu); + } else if (strcasecmp(phuType, "CELL") == 0) { + psString cellName = phuNameFromHeader("CELL.NAME", fileInfo, phu); // Name of the cell + int cellNum = pmChipFindCell(chip, cellName); // Cell number + if (cellNum == -1) { + psError(PS_ERR_IO, true, "Unable to find cell %s in chip %s.\n", cellName, chipName); + psFree(view); + return NULL; + } + cell = chip->cells->data[cellNum]; // The specified cell + view->cell = cellNum; + psFree(cellName); + + addHDUtoCell(cell, phdu); + } else { + psError(PS_ERR_IO, true, "The format of the PHU (%s) is not FPA, CHIP or CELL.\n", phuType); + psFree(view); + return NULL; + } + psFree(chipName); + } + psFree(phdu); + + pmFPALevel level = hduLevel(format);// The level at which to plug in HDU + + // Now go through the contents + psMetadataIterator *contentsIter = psMetadataIteratorAlloc(contents, PS_LIST_HEAD, NULL); + psMetadataItem *contentsItem = NULL; // Item from contents + while ((contentsItem = psMetadataGetAndIncrement(contentsIter))) { + const char *extName = contentsItem->name; + if (contentsItem->type != PS_DATA_STRING) { + psLogMsg(__func__, PS_LOG_WARN, "CONTENTS item %s is not of type STR --- ignored.\n", extName); + continue; + } + + pmHDU *hdu = pmHDUAlloc(extName); // The extension + hdu->format = psMemIncrRefCounter(format); + + processContents(fpa, chip, cell, hdu, level, contentsItem->data.V, format); + psFree(hdu); + } + psFree(contentsIter); + + pmConceptsReadFPA(fpa, PM_CONCEPT_SOURCE_DEFAULTS, NULL); + + return view; +} + + +// Print out the focal plane structure +void pmFPAPrint(pmFPA *fpa, // FPA to print + bool header, // Print headers? + bool concepts // Print concepts? + ) +{ + 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 read in for extension %s\n", fpa->hdu->extname); + } + if (header) { + if (fpa->hdu->header) { + psTrace(__func__, 2, "---> Header:\n"); + psMetadataPrint(fpa->hdu->header, 8); + } else { + psTrace(__func__, 2, "---> NO HEADER read in for extension %s\n", fpa->hdu->extname); + } + } + } + if (concepts) { + 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 (header) { + if (chip->hdu->header) { + psTrace(__func__, 4, "---> Header:\n"); + psMetadataPrint(chip->hdu->header, 8); + } else { + psTrace(__func__, 4, "---> NO HEADER read in for extension %s\n", chip->hdu->extname); + } + } + if (! chip->hdu->images) { + psTrace(__func__, 4, "---> NO PIXELS read in for extension %s\n", chip->hdu->extname); + } + } + if (concepts) { + 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 (header) { + if (cell->hdu->header) { + psTrace(__func__, 6, "---> Header:\n"); + psMetadataPrint(cell->hdu->header, 8); + } else { + psTrace(__func__, 6, "---> NO HEADER read in for extension %s\n", cell->hdu->extname); + } + } + if (! cell->hdu->images) { + psTrace(__func__, 6, "---> NO PIXELS read in for extension %s\n", cell->hdu->extname); + } + } + if (concepts) { + 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 + psTrace(__func__, 8, "row0: %d\n", readout->row0); + psImage *image = readout->image; // The image + psList *bias = readout->bias; // The list of bias images + if (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); + } + if (bias) { + psListIterator *biasIter = psListIteratorAlloc(bias, PS_LIST_HEAD, false); // Iterator + psImage *biasImage = NULL; // Bias image from iteration + while ((biasImage = psListGetAndIncrement(biasIter))) { + psTrace(__func__, 8, "Bias: [%d:%d,%d:%d] (%dx%d)\n", biasImage->col0, + biasImage->col0 + biasImage->numCols, biasImage->row0, + biasImage->row0 + biasImage->numRows, biasImage->numCols, biasImage->numRows); + } + psFree(biasIter); + } + } // Iterating over cell + } // Iterating over chip + } // Iterating over FPA + +} Index: trunk/psModules/src/camera/pmFPAConstruct.h =================================================================== --- trunk/psModules/src/camera/pmFPAConstruct.h (revision 7017) +++ trunk/psModules/src/camera/pmFPAConstruct.h (revision 7017) @@ -0,0 +1,28 @@ +#ifndef PM_FPA_CONSTRUCT_H +#define PM_FPA_CONSTRUCT_H + +#include "pslib.h" +#include "pmFPA.h" +#include "pmFPAview.h" + +// Construct an FPA instance on the basis of a camera configuration +pmFPA *pmFPAConstruct(const psMetadata *camera // The camera configuration + ); + +bool pmFPAAddSourceFromView(pmFPA *fpa, // The FPA + const pmFPAview *phuView, // The view, corresponding to the PHU + psMetadata *format // Format of file + ); + +pmFPAview *pmFPAAddSourceFromHeader(pmFPA *fpa, // The FPA + psMetadata *phu, // Primary header of file + psMetadata *format // Format of file + ); + +// Print out the FPA +void pmFPAPrint(pmFPA *fpa, // FPA to print + bool header, // Print headers? + bool concepts // Print concepts? + ); + +#endif Index: trunk/psModules/src/camera/pmFPACopy.c =================================================================== --- trunk/psModules/src/camera/pmFPACopy.c (revision 7017) +++ trunk/psModules/src/camera/pmFPACopy.c (revision 7017) @@ -0,0 +1,693 @@ +#include +#include + +#include "pslib.h" +#include "psImageFlip.h" +#include "psRegionIsBad.h" + +#include "pmFPA.h" +#include "pmFPAUtils.h" +#include "pmHDU.h" +#include "pmHDUUtils.h" +#include "pmFPACopy.h" + +#define MAX(x,y) ((x) > (y) ? (x) : (y)) + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// File-static functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Add cells in a chip to a list +static bool addCellsFromChip(psList *list, // List of cells + const pmChip *chip // The chip from which to add cells + ) +{ + assert(list); + assert(chip); + + psArray *cells = chip->cells; // Array of cells + bool result = true; // Result of adding cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // A cell + result |= psListAdd(list, PS_LIST_TAIL, cell); + } + + return result; +} + +// Add cells in an FPA to a list +static bool addCellsFromFPA(psList *list, // List of cells + const pmFPA *fpa // The FPA from which to add cells + ) +{ + assert(list); + assert(fpa); + + psArray *chips = fpa->chips; // Array of chips + bool result = true; // Result of adding cells + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // A chip + result |= addCellsFromChip(list, chip); + } + + return result; +} + +// Get a list of cells that share the HDU for the target cell +static bool cellList(psList *targets, // The list of target cells + psList *sources, // The list of source cells + pmCell *targetCell, // The target cell + pmCell *sourceCell // The source cell + ) +{ + assert(targetCell); + assert(sourceCell); + + if (targetCell->hdu) { + if (targets) { + psListAdd(targets, PS_LIST_TAIL, targetCell); + } + if (sources) { + psListAdd(sources, PS_LIST_TAIL, sourceCell); + } + } else { + pmChip *targetChip = targetCell->parent; // The target parent chip + pmChip *sourceChip = sourceCell->parent; // The source parent chip + if (targetChip->hdu) { + if (targets) { + addCellsFromChip(targets, targetChip); + } + if (sources) { + addCellsFromChip(sources, sourceChip); + } + } else { + pmFPA *targetFPA = targetChip->parent; // The target parent FPA + pmFPA *sourceFPA = sourceChip->parent; // The source parent FPA + if (targetFPA->hdu) { + if (targets) { + addCellsFromFPA(targets, targetFPA); + } + if (sources) { + addCellsFromFPA(sources, sourceFPA); + } + } else { + psError(PS_ERR_IO, true, "Unable to find HDU for cell to generate list!\n"); + return false; + } + } + } + + return true; +} + +// Get the maximum extent of the HDU from the trimsec and biassecs +static bool sizeHDU(int *xSize, int *ySize, // Size of HDU + psList *cells // List of cells + ) +{ + psListIterator *cellsIter = psListIteratorAlloc(cells, PS_LIST_HEAD, false); // Iterator for cells + pmCell *cell = NULL; // The cell from iteration + bool mdok = true; // Status of MD lookup + *xSize = 0; + *ySize = 0; + while ((cell = psListGetAndIncrement(cellsIter))) { + psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section + if (mdok && trimsec && !psRegionIsBad(*trimsec)) { + *xSize = MAX(trimsec->x1, *xSize); + *ySize = MAX(trimsec->y1, *ySize); + } else { + psFree(cellsIter); + return false; + } + psList *biassecs = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.BIASSEC"); // Bias sections + if (mdok && biassecs) { + psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); // Iterator + psRegion *biassec = NULL; // The bias section + while ((biassec = psListGetAndIncrement(biassecsIter))) { + if (!psRegionIsBad(*trimsec)) { + *xSize = MAX(biassec->x1, *xSize); + *ySize = MAX(biassec->y1, *ySize); + } else { + psFree(biassecsIter); + psFree(cellsIter); + return false; + } + } + psFree(biassecsIter); + } + } + psFree(cellsIter); + + return (*xSize != 0 && *ySize != 0); +} + + +static psRegion sectionForImage(int *position, // Position on the output image, updated + const psImage *image, // Image containing the sizes and offsets + int readdir // Read direction, 1=rows, 2=cols + ) +{ + psRegion region; + switch (readdir) { + case 1: // Read direction is rows + region = psRegionSet(*position, *position + image->numCols, image->row0, + image->row0 + image->numRows); + *position += image->numCols; + break; + case 2: // Read direction is columns + region = psRegionSet(image->col0, image->col0 + image->numCols, *position, + *position + image->numRows); + *position += image->numRows; + break; + default: + psAbort(__func__, "Shouldn't ever get here!\n"); + } + + return region; +} + +static bool doBiasSections(int *position, // Position on the output image, updated + pmCell *target, // Target cell + const pmCell *source // Source cell + ) +{ + psMetadataItem *biassecItem = psMetadataLookup(target->concepts, "CELL.BIASSEC"); // Bias sections + if (!biassecItem) { + psLogMsg(__func__, PS_LOG_WARN, "CELL.BIASSEC has not been initialised in target cell --- " + "ignored.\n"); + return false; + } + psFree(biassecItem->data.V); // Blow away the old list + psList *biassecs = psListAlloc(NULL); + biassecItem->data.V = biassecs; + + bool mdok = true; // Status of MD lookup + int readdir = psMetadataLookupS32(&mdok, source->concepts, "CELL.READDIR"); // Read direction + if (!mdok || (readdir != 1 && readdir != 2)) { + // Probably unnecessary, but just in case.... + psLogMsg(__func__, PS_LOG_WARN, "CELL.READDIR is not set in source cell --- ignored.\n"); + return false; + } + + pmReadout *readout = source->readouts->data[0]; // The first source readout, as representative + psList *biases = readout->bias; // The bias images from the source readout + + psListIterator *biasIter = psListIteratorAlloc(biases, PS_LIST_HEAD, true); // Iterator for biases + psImage *bias = NULL; // Bias image from iteration + while ((bias = psListGetAndIncrement(biasIter))) { + // Construct a region + psRegion *biassec = psAlloc(sizeof(psRegion)); // The new region; need a psMemBlock + *biassec = sectionForImage(position, bias, readdir); + psListAdd(biassecs, PS_LIST_TAIL, biassec); + psFree(biassec); // Drop reference + } + psFree(biasIter); + + return true; +} + +// Generate CELL.TRIMSEC and CELL.BIASSEC for the target cells +static bool generateTrimBias(psList *targets, // List of target cells + psList *sources // List of source cells + ) +{ + pmCell *target = NULL, *source = NULL; // Cells from iteration + int numCells = targets->n; // Number of cells + int cellNum = 0; // The cell number + int position = 0; // Position on the image + bool mdok = true; // Status of MD lookup + + // First run through to do the LHS biases + psListIterator *targetsIter = psListIteratorAlloc(targets, PS_LIST_HEAD, false); // Iterator for targets + psListIterator *sourcesIter = psListIteratorAlloc(sources, PS_LIST_HEAD, false); // Iterator for sources + bool done = false; // Done with iteration? + while ((target = psListGetAndIncrement(targetsIter)) && (source = psListGetAndIncrement(sourcesIter)) && + !done) { + if (cellNum <= numCells/2 - 1) { + doBiasSections(&position, target, source); + cellNum++; + } else { + done = true; + } + } + + // Second run through to do the trim sections + psListIteratorSet(targetsIter, PS_LIST_HEAD); + psListIteratorSet(sourcesIter, PS_LIST_HEAD); + while ((target = psListGetAndIncrement(targetsIter)) && (source = psListGetAndIncrement(sourcesIter))) { + psRegion *trimsec = psMetadataLookupPtr(&mdok, target->concepts, "CELL.TRIMSEC"); // Trim section + if (!mdok || !trimsec) { + psLogMsg(__func__, PS_LOG_WARN, "CELL.TRIMSEC has not been initialised in target cell --- " + "ignored.\n"); + continue; + } + + int readdir = psMetadataLookupS32(&mdok, source->concepts, "CELL.READDIR"); // Read direction + if (!mdok || (readdir != 1 && readdir != 2)) { + // Probably unnecessary, but just in case.... + psLogMsg(__func__, PS_LOG_WARN, "CELL.READDIR is not set in source cell --- ignored.\n"); + continue; + } + + pmReadout *readout = source->readouts->data[0]; // The first source readout, as representative + psImage *image = readout->image;// The proper image + *trimsec = sectionForImage(&position, image, readdir); + } + + // A final run through to do the RHS biases + psListIteratorSet(targetsIter, cellNum); + psListIteratorSet(sourcesIter, cellNum); + while ((target = psListGetAndIncrement(targetsIter)) && (source = psListGetAndIncrement(sourcesIter))) { + doBiasSections(&position, target, source); + } + + // Clean up + psFree(targetsIter); + psFree(sourcesIter); + + return (position > 0); +} + + +// Generate an HDU with the pixels +static bool generateHDU(pmCell *target, // The target cell + pmCell *source, // The source cell + int xBin, int yBin // Binning in x and y + ) +{ + // Get the HDU and a list of cells below it + pmHDU *hdu = pmHDUFromCell(target); // The HDU in the target cell + psList *targetCells = psListAlloc(NULL); // List of target cells below the target HDU + psList *sourceCells = psListAlloc(NULL); // List of source cells below the target HDU + if (! cellList(targetCells, sourceCells, target, source)) { + psError(PS_ERR_IO, true, "Unable to find cells to generate HDU!\n"); + return false; + } + + // Check the number of readouts + int numReadouts = -1; // Number of readouts + { + psListIterator *iter = psListIteratorAlloc(sourceCells, PS_LIST_HEAD, false); // Iterator for cells + pmCell *cell = NULL; // The cell from iteration + while ((cell = psListGetAndIncrement(iter))) + { + psArray *readouts = cell->readouts; + if (numReadouts == -1) { + numReadouts = readouts->n; + } else if (readouts->n != numReadouts) { + psError(PS_ERR_IO, true, "Number of readouts doesn't match: %d vs %d\n", readouts->n, + numReadouts); + return false; + } + + } + psFree(iter); + } + + // Get the size of the HDU, either from existing trimsec and biassec, or generate these and try again + int xSize = 0, ySize = 0; // Size of HDU + if (!sizeHDU(&xSize, &ySize, targetCells) && !(generateTrimBias(targetCells, sourceCells) && + sizeHDU(&xSize, &ySize, targetCells))) { + psError(PS_ERR_IO, true, "Unable to determine size of HDU!\n"); + return false; + } + psFree(targetCells); + psFree(sourceCells); + + xSize = (int)ceilf((float)xSize/(float)xBin); + ySize = (int)ceilf((float)ySize/(float)yBin); + + hdu->images = psArrayAlloc(numReadouts); + hdu->images->n = numReadouts; + for (int i = 0; i < numReadouts; i++) { + psImage *image = psImageAlloc(xSize, ySize, PS_TYPE_F32); + psImageInit(image, 0.0); + hdu->images->data[i] = image; + } + + return true; +} + +// Copy pixels from a target image to a source image, with flips +static bool copyPixels(psImage *target, // Target image (HDU pixels) + psImage *source, // Source image (from source cell) + psRegion region, // Region for pasting + bool xFlip, // Flip in x? + bool yFlip // Flip in y? + ) +{ + psImage *overlay = psMemIncrRefCounter(source); + if (xFlip) { + psImage *temp = psImageFlipX(overlay); + psFree(overlay); + overlay = temp; + } + if (yFlip) { + psImage *temp = psImageFlipY(overlay); + psFree(overlay); + overlay = temp; + } + int numPix = psImageOverlaySection(target, overlay, region.x0, region.y0, "="); + psFree(overlay); + return (numPix > 0); +} + +// Bin a region down by specified factors in x and y +static void binRegion(psRegion *region, // Region to be binned + int xBin, int yBin// Binning in x and y + ) +{ + // Want to include the lower bound: 1 binned by 4 --> 0; 3 binned by 4 --> 0; 4 binned by 4 --> 1 + region->x0 = (int)(region->x0 / xBin); + region->y0 = (int)(region->y0 / yBin); + // Want to exclude the upper bound: 4 binned by 4 --> 1; 5 binned by 4 --> 2; 7 binned by 4 --> 2 + region->x1 = (int)((region->x1 + xBin - 1) / xBin); + region->y1 = (int)((region->y1 + yBin - 1) / yBin); +} + + +static pmHDU *findPHU(const pmCell *cell// The cell for which to find the PHU + ) +{ + if (cell->hdu && cell->hdu->phu) { + return cell->hdu; + } + pmChip *chip = cell->parent; // The parent chip + if (chip->hdu && chip->hdu->phu) { + return chip->hdu; + } + pmFPA *fpa = chip->parent; // The parent FPA + if (fpa->hdu && fpa->hdu->phu) { + return fpa->hdu; + } + + psError(PS_ERR_IO, true, "Unable to find PHU for target cell.\n"); + return NULL; +} + + +static bool copyPHU(pmCell *targetCell, // The target cell + pmCell *sourceCell // The source cell + ) +{ + // Find the respective PHUs + pmHDU *targetPHU = findPHU(targetCell); // The target PHU + if (targetPHU->header) { + return false; // No work required + } + // Copy the header over + pmHDU *sourcePHU = findPHU(sourceCell); // The source PHU + targetPHU->header = psMetadataCopy(NULL, sourcePHU->header); + return true; +} + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// File-static engine functions --- these do all the work. Actually, cellCopy does all the work; the others +// merely iterate on the higher-level components. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +static int cellCopy(pmCell *target, // The target cell + pmCell *source, // The source cell, to be copied + bool pixels, // Copy the pixels? + int xBin, int yBin // (Relative) binning factors in x and y + ) +{ + assert(target); + assert(source); + + psArray *sourceReadouts = source->readouts; // The source readouts + int numReadouts = sourceReadouts->n; // Number of readouts copied + + // Copy any headers + if (target->hdu && !target->hdu->phu) { + pmHDU *sourceHDU = pmHDUFromCell(source); + target->hdu->header = psMetadataCopy(target->hdu->header, sourceHDU->header); + } + + pmHDU *hdu = pmHDUFromCell(target); // The target HDU; we need to fix this up + if (!hdu->images) { + generateHDU(target, source, xBin, yBin); + } + if (!hdu->header) { + hdu->header = psMetadataAlloc(); + } + // Copy the PHU over as well, if required + copyPHU(target, source); + + // Need to check/change CELL.XPARITY and CELL.YPARITY + bool xFlip = (psMetadataLookupS32(NULL, target->concepts, "CELL.XPARITY") != + psMetadataLookupS32(NULL, source->concepts, "CELL.XPARITY")); // Switch parity in x? + bool yFlip = (psMetadataLookupS32(NULL, target->concepts, "CELL.YPARITY") != + psMetadataLookupS32(NULL, source->concepts, "CELL.YPARITY")); // Switch parity in y? + psTrace(__func__, 3, "xFlip: %d; yFlip: %d\n", xFlip, yFlip); + + bool mdok = true; // Status of MD lookup + psRegion *trimsec = psMetadataLookupPtr(&mdok, target->concepts, "CELL.TRIMSEC"); // The trim section + if (!mdok || !trimsec || psRegionIsBad(*trimsec)) { + psError(PS_ERR_IO, true, "CELL.TRIMSEC isn't set!\n"); + return 0; + } + psList *biassecs = psMetadataLookupPtr(&mdok, target->concepts, "CELL.BIASSEC"); // The bias sections + if (!mdok || !biassecs) { + psError(PS_ERR_IO, true, "CELL.BIASSEC isn't set!\n"); + return 0; + } + + for (int i = 0; i < numReadouts; i++) { + pmReadout *sourceReadout = sourceReadouts->data[i]; // The source readout + psImage *sourceImage = sourceReadout->image; // The source image + pmReadout *targetReadout = pmReadoutAlloc(target); // The target readout; this adds it to the cell + if (sourceImage->numCols != trimsec->x1 - trimsec->x0 || + sourceImage->numRows != trimsec->y1 - trimsec->y0) { + psString trimsecString = psRegionToString(*trimsec); // String with the trim section + psLogMsg(__func__, PS_LOG_WARN, "Source image size (%dx%d) for readout %d doesn't match " + "CELL.TRIMSEC for target (%s) -- ignored.\n", sourceImage->numCols, sourceImage->numRows, + i, trimsecString); + psFree(trimsecString); + } else { + binRegion(trimsec, xBin, yBin); + if (pixels) { + copyPixels(hdu->images->data[i], sourceImage, *trimsec, xFlip, yFlip); + } + targetReadout->image = psImageSubset(hdu->images->data[i], *trimsec); + } + + psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); // Iterator + psRegion *biassec = NULL; // Bias section from iteration + psListIterator *biasIter = psListIteratorAlloc(sourceReadout->bias, PS_LIST_HEAD, false); // Iterator + psImage *bias = NULL; // Bias image from iteration + while ((biassec = psListGetAndIncrement(biassecsIter)) && (bias = psListGetAndIncrement(biasIter))) { + if (psRegionIsBad(*biassec)) { + psString biassecString = psRegionToString(*biassec); // String for bias section + psLogMsg(__func__, PS_LOG_WARN, "Bias section (%s) isn't set --- ignored.\n", biassecString); + psFree(biassecString); + continue; + } + if (bias->numCols != biassec->x1 - biassec->x0 || + bias->numRows != biassec->y1 - biassec->y0) { + psString biassecString = psRegionToString(*biassec); // String with the bias section + psLogMsg(__func__, PS_LOG_WARN, "Source image size (%dx%d) for readout %d doesn't match " + "CELL.BIASSEC for target (%s) -- ignored.\n", bias->numCols, bias->numRows, i, + biassecString); + psFree(biassecString); + } else { + binRegion(biassec, xBin, yBin); + if (pixels) { + copyPixels(hdu->images->data[i], bias, *biassec, xFlip, yFlip); + } + psImage *newBias = psImageSubset(hdu->images->data[i], *biassec); + psListAdd(targetReadout->bias, PS_LIST_TAIL, newBias); + psFree(newBias); // Drop reference + } + } + psFree(targetReadout); // Drop reference + psFree(biassecsIter); + psFree(biasIter); + } + + // Copy the remaining "concepts" over + psMetadataIterator *conceptsIter = psMetadataIteratorAlloc(source->concepts, PS_LIST_HEAD, NULL); + psMetadataItem *conceptItem = NULL; // Item from iteration + while ((conceptItem = psMetadataGetAndIncrement(conceptsIter))) { + psString name = conceptItem->name; // Name of concept + if (strcmp(name, "CELL.TRIMSEC") != 0 && strcmp(name, "CELL.BIASSEC") != 0 && + strcmp(name, "CELL.XPARITY") != 0 && strcmp(name, "CELL.YPARITY") != 0 && + strcmp(name, "CELL.X0") != 0 && strcmp(name, "CELL.Y0") != 0) { + psMetadataAddItem(target->concepts, conceptItem, PS_LIST_TAIL, PS_META_REPLACE); + } + } + psFree(conceptsIter); + + // Need to update CELL.X0 and CELL.Y0 if we flipped + if (xFlip) { + int xZero = psMetadataLookupS32(NULL, source->concepts, "CELL.X0"); // CELL.X0 from source + int xParity = psMetadataLookupS32(NULL, target->concepts, "CELL.XPARITY"); // Parity in x + int xBin = psMetadataLookupS32(NULL, source->concepts, "CELL.XBIN"); // CELL.XBIN from source + pmReadout *readout = source->readouts->data[0]; // A representative readout + psTrace(__func__, 3, "CELL.X0: Before: %d After: %d\n", xZero, + xZero - (readout->image->numCols - 1) * xParity * xBin); + psTrace(__func__, 9, "(xParity: %d xBin: %d numCols: %d)\n", xParity, xBin, readout->image->numCols); + xZero -= (readout->image->numCols - 1) * xParity * xBin; // Change the parity on the X0 position + psMetadataItem *newItem = psMetadataLookup(target->concepts, "CELL.X0"); // CELL.X0 from target + newItem->data.S32 = xZero; + } + if (yFlip) { + int yZero = psMetadataLookupS32(NULL, source->concepts, "CELL.Y0"); // CELL.Y0 from source + int yParity = psMetadataLookupS32(NULL, target->concepts, "CELL.YPARITY"); // Parity in y + int yBin = psMetadataLookupS32(NULL, source->concepts, "CELL.YBIN"); // Parity in y + pmReadout *readout = source->readouts->data[0]; // A representative readout + psTrace(__func__, 3, "CELL.Y0: Before: %d After: %d\n", yZero, + yZero - (readout->image->numRows - 1) * yParity * yBin); + psTrace(__func__, 9, "(yParity: %d yBin: %d numRows: %d)\n", yParity, yBin, readout->image->numRows); + yZero -= (readout->image->numRows - 1) * yParity * yBin; // Change the parity on the Y0 position + psMetadataItem *newItem = psMetadataLookup(target->concepts, "CELL.Y0"); // CELL.Y0 from target + newItem->data.S32 = yZero; + } + + // Update the binning concepts + psMetadataItem *binItem = psMetadataLookup(target->concepts, "CELL.XBIN"); + binItem->data.S32 *= xBin; + binItem = psMetadataLookup(target->concepts, "CELL.YBIN"); + binItem->data.S32 *= yBin; + + return numReadouts; +} + +static int chipCopy(pmChip *target, // The target chip + pmChip *source, // The source chip, to be copied + bool pixels, // Copy the pixels? + int xBin, int yBin // (Relative) binning factors in x and y + ) +{ + assert(target); + assert(source); + + psArray *targetCells = target->cells; // The target cells + psArray *sourceCells = source->cells; // The source cells + if (targetCells->n != sourceCells->n) { + psError(PS_ERR_IO, true, "Number of source cells (%d) differs from the number of target cells (%d)\n", + sourceCells->n, targetCells->n); + return false; + } + + // Copy any headers + if (target->hdu && !target->hdu->phu) { + pmHDU *sourceHDU = pmHDUFromChip(source); + target->hdu->header = psMetadataCopy(target->hdu->header, sourceHDU->header); + } + + int numCells = 0; // Number of cells copied + for (int i = 0; i < targetCells->n; i++) { + pmCell *targetCell = targetCells->data[i]; // The target cell + const char *cellName = psMetadataLookupStr(NULL, targetCell->concepts, "CELL.NAME"); // Name of cell + int cellNum = pmChipFindCell(source, cellName); // Number of cell with that name + if (cellNum >= 0) { + pmCell *sourceCell = sourceCells->data[cellNum]; // The source cell + int numReadouts = cellCopy(targetCell, sourceCell, pixels, xBin, yBin); // Number of readouts + // copied + psTrace(__func__, 5, "Copied %d readouts for cell %s\n", numReadouts, cellName); + numCells++; + } + } + + // Update the concepts + psMetadataCopy(target->concepts, source->concepts); + + return numCells; + +} + +static int fpaCopy(pmFPA *target, // The target FPA + pmFPA *source, // The source FPA, to be copied + bool pixels, // Copy the pixels? + int xBin, int yBin // (Relative) binning factors in x and y + ) +{ + assert(target); + assert(source); + + psArray *targetChips = target->chips; // The target chips + psArray *sourceChips = source->chips; // The source chips + if (targetChips->n != sourceChips->n) { + psError(PS_ERR_IO, true, "Number of source chips (%d) differs from the number of target chips (%d)\n", + sourceChips->n, targetChips->n); + return false; + } + + // Copy any headers + if (target->hdu && !target->hdu->phu) { + pmHDU *sourceHDU = pmHDUFromFPA(source); + target->hdu->header = psMetadataCopy(target->hdu->header, sourceHDU->header); + } + + int numChips = 0; // Number of chips copied + for (int i = 0; i < targetChips->n; i++) { + pmChip *targetChip = targetChips->data[i]; // The target chip + const char *chipName = psMetadataLookupStr(NULL, targetChip->concepts, "CHIP.NAME"); // Name of chip + int chipNum = pmFPAFindChip(source, chipName); // Number of chip with that name + if (chipNum >= 0) { + pmChip *sourceChip = sourceChips->data[chipNum]; // The source chip + int numCells = chipCopy(targetChip, sourceChip, pixels, xBin, yBin); // Number of cells copied + psTrace(__func__, 5, "Copied %d cells for chip %s\n", numCells, chipName); + numChips++; + } + } + + // Update the concepts + psMetadataCopy(target->concepts, source->concepts); + + return numChips; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + +int pmFPACopy(pmFPA *target, // The target FPA + pmFPA *source // The source FPA, to be copied + ) +{ + return fpaCopy(target, source, true, 1, 1); +} + +int pmChipCopy(pmChip *target, // The target chip + pmChip *source // The source chip, to be copied + ) +{ + return chipCopy(target, source, true, 1, 1); +} + +int pmCellCopy(pmCell *target, // The target cell + pmCell *source // The source cell, to be copied + ) +{ + return cellCopy(target, source, true, 1, 1); +} + + +int pmFPACopyStructure(pmFPA *target, // The target FPA + pmFPA *source, // The source FPA, to be copied + int xBin, int yBin // Binning factors in x and y + ) +{ + return fpaCopy(target, source, false, xBin, yBin); +} + +int pmChipCopyStructure(pmChip *target, // The target chip + pmChip *source, // The source chip, to be copied + int xBin, int yBin // Binning factors in x and y + ) +{ + return chipCopy(target, source, false, xBin, yBin); +} + +int pmCellCopyStructure(pmCell *target, // The target cell + pmCell *source, // The source cell, to be copied + int xBin, int yBin // Binning factors in x and y + ) +{ + return cellCopy(target, source, false, xBin, yBin); +} Index: trunk/psModules/src/camera/pmFPACopy.h =================================================================== --- trunk/psModules/src/camera/pmFPACopy.h (revision 7017) +++ trunk/psModules/src/camera/pmFPACopy.h (revision 7017) @@ -0,0 +1,30 @@ +#ifndef PM_FPA_COPY_H +#define PM_FPA_COPY_H + +// Copy the FPA components, including the pixels +int pmFPACopy(pmFPA *target, // The target FPA + pmFPA *source // The source FPA, to be copied + ); +int pmChipCopy(pmChip *target, // The target chip + pmChip *source // The source chip, to be copied + ); +int pmCellCopy(pmCell *target, // The target cell + pmCell *source // The source cell, to be copied + ); + +// Versions that copy the structure and not the pixels; they also allow binning +int pmFPACopyStructure(pmFPA *target, // The target FPA + pmFPA *source, // The source FPA, to be copied + int xBin, int yBin // Binning factors in x and y + ); +int pmChipCopyStructure(pmChip *target, // The target chip + pmChip *source, // The source chip, to be copied + int xBin, int yBin // Binning factors in x and y + ); +int pmCellCopyStructure(pmCell *target, // The target cell + pmCell *source, // The source cell, to be copied + int xBin, int yBin // Binning factors in x and y + ); + + +#endif Index: trunk/psModules/src/camera/pmFPAHeader.c =================================================================== --- trunk/psModules/src/camera/pmFPAHeader.c (revision 7017) +++ trunk/psModules/src/camera/pmFPAHeader.c (revision 7017) @@ -0,0 +1,99 @@ +#include +#include +#include "pslib.h" +#include "pmFPA.h" +#include "pmHDU.h" +#include "pmConcepts.h" +#include "pmFPAHeader.h" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// File-static (private) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Read concepts for all cells in a chip +static bool chipConcepts(pmChip *chip // The chip for which to read the concepts + ) +{ + bool status = true; // Status of concept reading + status |= pmConceptsReadChip(chip, PM_CONCEPT_SOURCE_HEADER, false, NULL); + psArray *cells = chip->cells; // The cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // The cell of interest + if (!cell) { + continue; + } + status |= pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_HEADER, false, NULL); + } + + return status; +} + +// Read concepts for all chips in an FPA +static bool fpaConcepts(pmFPA *fpa // The FPA for which to read the concepts + ) +{ + bool status = true; // Status of concept reading + status |= pmConceptsReadFPA(fpa, PM_CONCEPT_SOURCE_HEADER, NULL); + psArray *chips = fpa->chips; // The chips + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // The chip of interest + if (!chip) { + continue; + } + status |= chipConcepts(chip); + } + + return status; +} + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +bool pmCellReadHeader(pmCell *cell, // Cell for which to read header + psFits *fits // FITS file handle + ) +{ + if (!cell->hdu) { + return pmChipReadHeader(cell->parent, fits); + } + if (!pmHDUReadHeader(cell->hdu, fits)) { + psError(PS_ERR_IO, false, "Unable to read header for cell.\n"); + return false; + } + + return pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_HEADER, false, NULL); +} + + +bool pmChipReadHeader(pmChip *chip, // Chip for which to read header + psFits *fits // FITS file handle + ) +{ + if (!chip->hdu) { + return pmFPAReadHeader(chip->parent, fits); + } + if (!pmHDUReadHeader(chip->hdu, fits)) { + psError(PS_ERR_IO, false, "Unable to read header for cell.\n"); + return false; + } + + return chipConcepts(chip); +} + + +bool pmFPAReadHeader(pmFPA *fpa, // FPA for which to read header + psFits *fits // FITS file handle + ) +{ + if (!fpa->hdu) { + return false; + } + if (!pmHDUReadHeader(fpa->hdu, fits)) { + psError(PS_ERR_IO, false, "Unable to read header for cell.\n"); + return false; + } + + return fpaConcepts(fpa); +} Index: trunk/psModules/src/camera/pmFPAHeader.h =================================================================== --- trunk/psModules/src/camera/pmFPAHeader.h (revision 7017) +++ trunk/psModules/src/camera/pmFPAHeader.h (revision 7017) @@ -0,0 +1,14 @@ +#ifndef PM_FPA_HEADER_H +#define PM_FPA_HEADER_H + +bool pmCellReadHeader(pmCell *cell, // Cell for which to read header + psFits *fits // FITS file handle + ); +bool pmChipReadHeader(pmChip *chip, // Chip for which to read header + psFits *fits // FITS file handle + ); +bool pmFPAReadHeader(pmFPA *fpa, // FPA for which to read header + psFits *fits // FITS file handle + ); + +#endif Index: trunk/psModules/src/camera/pmFPAMaskWeight.c =================================================================== --- trunk/psModules/src/camera/pmFPAMaskWeight.c (revision 7017) +++ trunk/psModules/src/camera/pmFPAMaskWeight.c (revision 7017) @@ -0,0 +1,101 @@ +#include +#include "pslib.h" +#include "pmFPA.h" +#include "pmFPAMaskWeight.h" + +bool pmReadoutSetMask(pmReadout *readout // Readout for which to set mask + ) +{ + PS_ASSERT_PTR_NON_NULL(readout, false); + + pmCell *cell = readout->parent; // The parent cell + bool mdok = true; // Status of MD lookup + + // Get the "concepts" of interest + float saturation = psMetadataLookupF32(&mdok, cell->concepts, "CELL.SATURATION"); // Saturation level + if (!mdok || isnan(saturation)) { + psError(PS_ERR_IO, true, "CELL.SATURATION is not set --- unable to set mask.\n"); + return false; + } + float bad = psMetadataLookupF32(NULL, cell->concepts, "CELL.BAD"); // Bad level + if (!mdok || isnan(bad)) { + psError(PS_ERR_IO, true, "CELL.BAD is not set --- unable to set mask.\n"); + return false; + } + + psImage *image = readout->image; // The image pixels + psImage *mask = readout->mask; // The mask pixels + if (!mask) { + mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8); + readout->mask = mask; + psImageInit(mask, 0); + } + + // Set up the mask + for (int i = 0; i < image->numRows; i++) { + for (int j = 0; j < image->numCols; j++) { + if (image->data.F32[i][j] > saturation) { + mask->data.U8[i][j] |= PM_MASK_SAT; + } + if (image->data.F32[i][j] < bad) { + mask->data.U8[i][j] |= PM_MASK_BAD; + } + } + } + + return true; +} + +bool pmReadoutSetWeight(pmReadout *readout // Readout for which to set weight + ) +{ + PS_ASSERT_PTR_NON_NULL(readout, false); + + pmCell *cell = readout->parent; // The parent cell + bool mdok = true; // Status of MD lookup + + // Get the "concepts" of interest + float gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); // Cell gain + if (!mdok || isnan(gain)) { + psError(PS_ERR_IO, true, "CELL.GAIN is not set --- unable to set weight.\n"); + return false; + } + float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise + if (!mdok || isnan(readnoise)) { + psError(PS_ERR_IO, true, "CELL.READNOISE is not set --- unable to set weight.\n"); + return false; + } + + psImage *image = readout->image; // The image pixels + psImage *weight = readout->weight; // The weight pixels + if (!weight) { + weight = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); + readout->weight = weight; + psImageInit(weight, 0.0); + } + + // Set weight image to the variance in ADU = f/g + rn^2 + psBinaryOp(readout->weight, image, "/", psScalarAlloc(gain, PS_TYPE_F32)); + psBinaryOp(readout->weight, readout->weight, "+", + psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32)); + + return true; +} + + + +bool pmCellSetMaskWeight(pmCell *cell // Cell for which to set weights + ) +{ + PS_ASSERT_PTR_NON_NULL(cell, false); + + bool success = true; // Was everything successful? + psArray *readouts = cell->readouts; // Array of readouts + for (int i = 0; i < readouts->n; i++) { + pmReadout *readout = readouts->data[i]; // The readout + success |= pmReadoutSetMask(readout); + success |= pmReadoutSetWeight(readout); + } + + return success; +} Index: trunk/psModules/src/camera/pmFPAMaskWeight.h =================================================================== --- trunk/psModules/src/camera/pmFPAMaskWeight.h (revision 7017) +++ trunk/psModules/src/camera/pmFPAMaskWeight.h (revision 7017) @@ -0,0 +1,23 @@ +#ifndef PM_READOUT_MASK_WEIGHT_H +#define PM_READOUT_MASK_WEIGHT_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_BAD = 0x0008, ///< The pixel is low + PM_MASK_FLAT = 0x0010 ///< The pixel is non-positive in the flat-field. +} pmMaskValue; + +bool pmReadoutSetMask(pmReadout *readout // Readout for which to set mask + ); +bool pmReadoutSetWeight(pmReadout *readout // Readout for which to set weight + ); +bool pmCellSetMaskWeight(pmCell *cell // Cell for which to set weights + ); + + +#endif Index: trunk/psModules/src/camera/pmFPARead.c =================================================================== --- trunk/psModules/src/camera/pmFPARead.c (revision 7017) +++ trunk/psModules/src/camera/pmFPARead.c (revision 7017) @@ -0,0 +1,295 @@ +#include +#include +#include +#include "pslib.h" + +#include "pmFPA.h" +#include "pmFPARead.h" +#include "pmHDU.h" +#include "pmHDUUtils.h" +#include "pmConcepts.h" +#include "psRegionIsBad.h" +#include "pmFPAHeader.h" + +#define MAX(x,y) ((x) > (y) ? (x) : (y)) +#define MIN(x,y) ((x) < (y) ? (x) : (y)) + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// File-static functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Carve a readout from the image pixels +static bool readoutCarve(pmReadout *readout, // Readout to be carved up + psImage *image, // Image that will be carved + const psRegion *trimsec, // Trim section + const psList *biassecs // Bias sections + ) +{ + // The image corresponding to the trim region + if (psRegionIsBad(*trimsec)) { + psString regionString = psRegionToString(*trimsec); + psError(PS_ERR_IO, true, "Invalid trim section: %s\n", regionString); + psFree(regionString); + psFree(readout); + return NULL; + } + psRegion region = psRegionSet(MAX(trimsec->x0 - readout->col0, 0), // x0 + MIN(trimsec->x1 - readout->col0, image->numCols), // x1 + MAX(trimsec->y0 - readout->row0, 0), // y0 + MIN(trimsec->y1 - readout->row0, image->numRows) // y1 + ); + if (readout->image) { + psFree(readout->image); // Make way! + } + readout->image = psMemIncrRefCounter(psImageSubset(image, region)); + + // Get the list of overscans + if (readout->bias->n != 0) { + // Make way! + psFree(readout->bias); + readout->bias = psListAlloc(NULL); + } + psListIterator *iter = psListIteratorAlloc((psList*)biassecs, PS_LIST_HEAD, false); // Iterator + psRegion *biassec = NULL; // A BIASSEC region from the list + while ((biassec = psListGetAndIncrement(iter))) { + if (psRegionIsBad(*biassec)) { + psString regionString = psRegionToString(*biassec); + psError(PS_ERR_IO, true, "Invalid bias section: %s\n", regionString); + psFree(regionString); + psFree(readout); + return NULL; + } + psRegion region = psRegionSet(MAX(biassec->x0 - readout->col0, 0), // x0 + MIN(biassec->x1 - readout->col0, image->numCols), // x1 + MAX(biassec->y0 - readout->row0, 0), // y0 + MIN(biassec->y1 - readout->row0, image->numRows) // y1 + ); + psImage *overscan = psMemIncrRefCounter(psImageSubset(image, region)); + psListAdd(readout->bias, PS_LIST_TAIL, overscan); + psFree(overscan); + } + psFree(iter); + + return readout; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Read the next readout; return true if we read pixels in. +bool pmReadoutReadNext(pmReadout *readout, // Readout into which to read + psFits *fits, // FITS file from which to read + int z, // Readout number/plane; zero-offset indexing + int numScans // The number of scans to read + ) +{ + // Get the HDU and read the header + pmCell *cell = readout->parent; // The parent cell + if (!pmCellReadHeader(cell, fits)) { + psError(PS_ERR_IO, false, "Unable to read header for cell!\n"); + return false; + } + + pmHDU *hdu = pmHDUFromCell(cell); // The HDU + if (!hdu) { + return false; + } + + // Make sure we have the information we need + pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_HEADER, false, NULL); + + // Get the trim and bias sections + bool mdok = true; // Status of MD lookup + psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim sections + if (!mdok || !trimsec || psRegionIsBad(*trimsec)) { + psError(PS_ERR_IO, true, "CELL.TRIMSEC is not set.\n"); + return false; + } + psList *biassecs = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.BIASSEC"); // Bias sections + if (!mdok || !biassecs) { + psError(PS_ERR_IO, true, "CELL.BIASSEC is not set.\n"); + return false; + } + int readdir = psMetadataLookupS32(&mdok, cell->concepts, "CELL.READDIR"); // Read direction + if (!mdok || readdir == 0 || (readdir != -1 && readdir != 1)) { + psError(PS_ERR_IO, true, "CELL.READDIR is not set to -1 or +1.\n"); + return false; + } + + // Check the third dimension + int naxis = psMetadataLookupS32(&mdok, hdu->header, "NAXIS"); // The number of axes + if (!mdok) { + psError(PS_ERR_IO, true, "Unable to find NAXIS in header for extension %s\n", hdu->extname); + return false; + } + if (naxis == 0) { + // No pixels to read, as for a PHU. + return false; + } + if (naxis < 2 || naxis > 3) { + psError(PS_ERR_IO, true, "NAXIS in header of extension %s (= %d) is not valid.\n", + hdu->extname, naxis); + return false; + } + int naxis3 = 1; // The number of image planes + if (naxis == 3) { + naxis3 = psMetadataLookupS32(&mdok, hdu->header, "NAXIS3"); + if (!mdok) { + psError(PS_ERR_IO, true, "Unable to find NAXIS3 in header for extension %s\n", hdu->extname); + return false; + } + } + if (z >= naxis3) { + // Nothing to see here. Move along. + return false; + } + + // Get the size of the image plane + int naxis1 = psMetadataLookupS32(&mdok, hdu->header, "NAXIS1"); // The number of columns + if (!mdok) { + psError(PS_ERR_IO, true, "Unable to find NAXIS1 in header for extension %s\n", hdu->extname); + return false; + } + int naxis2 = psMetadataLookupS32(&mdok, hdu->header, "NAXIS2"); // The number of columns + if (!mdok) { + psError(PS_ERR_IO, true, "Unable to find NAXIS2 in header for extension %s\n", hdu->extname); + return false; + } + + // Read the FITS image + int offset = readout->row0; // The row number to read + if (readout->image) { + if (readdir > 0) { + // Reading rows + offset += readout->image->numRows; + } else { + // Reading columns + offset += readout->image->numCols; + } + } + if ((readdir > 0 && offset >= naxis2) || (readdir < 0 && offset > naxis1)) { + // We've read everything there is + return false; + } + int upper = offset + numScans; // The upper limit for the pixel read + if (readdir > 0 && upper > naxis2) { + upper = naxis2; + } else if (readdir < 0 && upper > naxis1) { + upper = naxis1; + } + psRegion region = {0, 0, 0, 0}; // Region to be read + if (readdir > 0) { + region = psRegionSet(0, naxis1, offset, upper); // Read rows + } else { + region = psRegionSet(offset, upper, 0, naxis2); // Read columns + } + psImage *image = psFitsReadImage(fits, region, z); // The image + + // Stick the image into the HDU and the readout + if (!hdu->images) { + hdu->images = psArrayAlloc(naxis3); + hdu->images->n = naxis3; + } + if (hdu->images->nalloc != naxis3) { + hdu->images = psArrayRealloc(hdu->images, naxis3); + hdu->images->n = naxis3; + } + if (hdu->images->data[z]) { + psFree(hdu->images->data[z]); + } + hdu->images->data[z] = image; + readout->row0 = region.y0; + readout->col0 = region.x0; + readoutCarve(readout, image, trimsec, biassecs); + + return true; +} + + +// Read in the cell, and allocate the readouts +bool pmCellRead(pmCell *cell, // Cell to read into + psFits *fits, // FITS file from which to read + psDB *db // Database handle, for "concepts" ingest + ) +{ + pmHDU *hdu = pmHDUFromCell(cell); // The HDU + if (!hdu) { + return false; // Nothing to see here; move along + } + if (!hdu->images && !pmHDURead(hdu, fits)) { + psError(PS_ERR_IO, false, "Unable to read HDU for cell.\n"); + return false; + } + + pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_HEADER, false, NULL); + + // Having read the cell, we now have to cut it up + psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); + psList *biassecs = psMetadataLookupPtr(NULL, cell->concepts, "CELL.BIASSEC"); + + // Iterate over each of the image planes + for (int i = 0; i < hdu->images->n; i++) { + psImage *image = hdu->images->data[i]; // The i-th plane + + // XXX: Type conversion here to support the modules, which don't have multiple type support yet + if (image->type.type != PS_TYPE_F32) { + psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32); // Temporary image + psFree(hdu->images->data[i]); + hdu->images->data[i] = temp; + image = temp; + } + + pmReadout *readout = pmReadoutAlloc(cell); + readoutCarve(readout, image, trimsec, biassecs); + readout->mask = NULL; + readout->weight = NULL; + psFree(readout); // Drop reference + } + + pmCellSetDataStatus (cell, true); + return true; +} + +// Read in the component cells +bool pmChipRead(pmChip *chip, // Chip to read into + psFits *fits, // FITS file from which to read + psDB *db // Database handle, for "concepts" ingest + ) +{ + bool success = false; // Were we able to read at least one HDU? + psArray *cells = chip->cells; // Array of cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // The cell of interest + success |= pmCellRead(cell, fits, db); + } + if (success) { + pmConceptsReadChip(chip, PM_CONCEPT_SOURCE_HEADER, false, NULL); + // XXX probably could just use chip->data_exists + pmChipSetDataStatus (chip, true); + } + + return success; +} + +// Read in the component chips +bool pmFPARead(pmFPA *fpa, // FPA to read into + psFits *fits, // FITS file from which to read + psDB *db // Database handle, for "concepts" ingest + ) +{ + bool success = false; // Were we able to read at least one HDU? + psArray *chips = fpa->chips; // Array of chips + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // The cell of interest + success |= pmChipRead(chip, fits, db); + } + if (success) { + pmConceptsReadFPA(fpa, PM_CONCEPT_SOURCE_HEADER, NULL); + } else { + psLogMsg(__func__, PS_LOG_WARN, "Unable to read any chips in FPA.\n"); + } + + return success; +} + Index: trunk/psModules/src/camera/pmFPARead.h =================================================================== --- trunk/psModules/src/camera/pmFPARead.h (revision 7017) +++ trunk/psModules/src/camera/pmFPARead.h (revision 7017) @@ -0,0 +1,66 @@ +#ifndef PM_FPA_READ_H +#define PM_FPA_READ_H + +#include "pslib.h" +#include "pmFPA.h" + +bool pmReadoutReadNext(pmReadout *readout, // Readout into which to read + psFits *fits, // FITS file from which to read + int z, // Readout number/plane; zero-offset indexing + int numRows // The number of rows to read + ); + +bool pmCellRead(pmCell *cell, // Cell to read into + psFits *fits, // FITS file from which to read + psDB *db // Database handle, for "concepts" ingest + ); + +bool pmChipRead(pmChip *chip, // Chip to read into + psFits *fits, // FITS file from which to read + psDB *db // Database handle, for "concepts" ingest + ); + +bool pmFPARead(pmFPA *fpa, // FPA to read into + psFits *fits, // FITS file from which to read + psDB *db // Database handle, for "concepts" ingest + ); + +bool pmCellReadPHU(pmCell *cell, // Cell to read into + psFits *fits // FITS file from which to read + ); + +bool pmChipReadPHU(pmChip *chip, // Chip to read into + psFits *fits // FITS file from which to read + ); + +bool pmFPAReadPHU(pmFPA *fpa, // FPA to read into + psFits *fits // FITS file from which to read + ); + + +#if 0 +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 + +#endif Index: trunk/psModules/src/camera/pmFPAUtils.c =================================================================== --- trunk/psModules/src/camera/pmFPAUtils.c (revision 7017) +++ trunk/psModules/src/camera/pmFPAUtils.c (revision 7017) @@ -0,0 +1,41 @@ +#include +#include "pslib.h" +#include "pmFPA.h" +#include "pmFPAUtils.h" + +// Find a chip by name; return the index +int pmFPAFindChip(const pmFPA *fpa, // FPA in which to find the chip + const char *name // Name of the chip + ) +{ + psArray *chips = fpa->chips; // Array of chips + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // The chip of interest + psString testName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of this chip + if (strcmp(name, testName) == 0) { + return i; + } + } + + psError(PS_ERR_IO, true, "Unable to find chip %s\n", name); + return -1; +} + +// Find a cell by name; return the index +int pmChipFindCell(const pmChip *chip, // Chip in which to find the cell + const char *name // Name of the cell + ) +{ + psArray *cells = chip->cells; // Array of cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // The cell of interest + psString testName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of this cell + if (strcmp(name, testName) == 0) { + return i; + } + } + + psError(PS_ERR_IO, true, "Unable to find cell %s\n", name); + return -1; +} + Index: trunk/psModules/src/camera/pmFPAUtils.h =================================================================== --- trunk/psModules/src/camera/pmFPAUtils.h (revision 7017) +++ trunk/psModules/src/camera/pmFPAUtils.h (revision 7017) @@ -0,0 +1,17 @@ +#ifndef PM_FPA_UTILS_H +#define PM_FPA_UTILS_H + +#include "pmFPA.h" + +// Find a chip by name; return the index +int pmFPAFindChip(const pmFPA *fpa, // FPA in which to find the chip + const char *name // Name of the chip + ); + +// Find a cell by name; return the index +int pmChipFindCell(const pmChip *chip, // Chip in which to find the cell + const char *name // Name of the cell + ); + + +#endif Index: trunk/psModules/src/camera/pmFPAWrite.c =================================================================== --- trunk/psModules/src/camera/pmFPAWrite.c (revision 7017) +++ trunk/psModules/src/camera/pmFPAWrite.c (revision 7017) @@ -0,0 +1,153 @@ +#include +#include +#include "pslib.h" + +#include "pmFPA.h" +#include "pmHDU.h" +#include "pmHDUUtils.h" +#include "pmConcepts.h" + + +bool pmReadoutWriteNext(pmReadout *readout, // Readout to write + psFits *fits, // FITS file to which to write + int z // Image plane to write + ) +{ + pmHDU *hdu = pmHDUFromReadout(readout); // The HDU to which to write + psMetadata *header = hdu->header; // The FITS header + if (!header) { + psError(PS_ERR_IO, true, "No FITS header available in the HDU.\n"); + return false; + } + + // We have to rely to a great extent on the FITS header, because otherwise we simply don't know how many + // image planes there are (NAXIS3) or the size of the original image (NAXIS1, NAXIS2). + bool mdok = true; // Status of MD lookup + int naxis1 = psMetadataLookupS32(&mdok, header, "NAXIS1"); // Number of columns + if (!mdok || naxis1 <= 0) { + psError(PS_ERR_IO, true, "Can't find NAXIS1 in header.\n"); + return false; + } + int naxis2 = psMetadataLookupS32(&mdok, header, "NAXIS2"); // Number of rows + if (!mdok || naxis2 <= 0) { + psError(PS_ERR_IO, true, "Can't find NAXIS2 in header.\n"); + return false; + } + int naxis3 = psMetadataLookupS32(&mdok, header, "NAXIS3"); // Number of image planes + if (!mdok || naxis3 <= 0) { + naxis3 = 1; + } + if (z >= naxis3) { + psError(PS_ERR_IO, true, "Specified a plane number (%d) greater than NAXIS3 allows.\n", z); + return false; + } + + psImage *image = hdu->images->data[z]; // The image from the HDU to write + if (readout->row0 == 0 && readout->col0 == 0 && z == 0) { + // Then we can assume that nothing has been written to the FITS file for now + if (naxis1 == image->numCols && naxis2 == image->numRows) { + // We can write the whole lot at once + return psFitsWriteImage(fits, header, image, z, hdu->extname); + } + // Create a dummy image so we can write something larger than we actually have + psImage *dummy = psImageAlloc(naxis1, naxis2, image->type.type); // Dummy image + psImageInit(dummy, 0); + psImageOverlaySection(dummy, image, 0, 0, "="); + bool result = psFitsWriteImage(fits, header, dummy, z, hdu->extname); + psFree(dummy); + return result; + } + + // We can simply update an existing HDU + if (!psFitsMoveExtName(fits, hdu->extname)) { + psError(PS_ERR_IO, false, "Unable to move to extension %s\n", hdu->extname); + return false; + } + return psFitsUpdateImage(fits, image, readout->col0, readout->row0, z); +} + + + + +bool pmCellWrite(pmCell *cell, // Cell to write + psFits *fits, // FITS file to which to write + psDB *db, // Database handle for "concepts" update + bool pixels // Write the pixels? + ) +{ + pmHDU *hdu = cell->hdu; // The HDU + if (hdu && ((!pixels && hdu->phu) || pixels)) { + bool status = pmConceptsWriteCell(cell, PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CAMERA | + PM_CONCEPT_SOURCE_DEFAULTS, false, NULL); + status &= pmHDUWrite(hdu, fits); + if (!status) { + psError(PS_ERR_IO, false, "Unable to write HDU for Chip.\n"); + return false; + } + } + + return true; +} + + +bool pmChipWrite(pmChip *chip, // Chip to write + psFits *fits, // FITS file to which to write + psDB *db, // Database handle for "concepts" update + bool pixels // Write the pixels? + ) +{ + pmHDU *hdu = chip->hdu; // The HDU + if (hdu && ((!pixels && hdu->phu) || pixels)) { + bool status = pmConceptsWriteChip(chip, PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CAMERA | + PM_CONCEPT_SOURCE_DEFAULTS, false, NULL); + status &= pmHDUWrite(hdu, fits); + if (!status) { + psError(PS_ERR_IO, false, "Unable to write HDU for Chip.\n"); + return false; + } + } + + bool success = true; // Success of writing + if (pixels) { + psArray *cells = chip->cells; // Array of cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // The cell of interest + success |= pmCellWrite(cell, fits, db, true); + } + } + + return success; +} + + + + +bool pmFPAWrite(pmFPA *fpa, // FPA to write + psFits *fits, // FITS file to which to write + psDB *db, // Database handle for "concepts" update + bool pixels // Write the pixels? + ) +{ + pmHDU *hdu = fpa->hdu; // The HDU + if (hdu && ((!pixels && hdu->phu) || pixels)) { + bool status = pmConceptsWriteFPA(fpa, PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CAMERA | + PM_CONCEPT_SOURCE_DEFAULTS, NULL); + status &= pmHDUWrite(hdu, fits); + if (!status) { + psError(PS_ERR_IO, false, "Unable to write HDU for FPA.\n"); + return false; + } + } + + bool success = true; // Success of writing + if (pixels) { + psArray *chips = fpa->chips; // Array of chips + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // The chip of interest + success |= pmChipWrite(chip, fits, db, true); + } + } + + return success; +} + Index: trunk/psModules/src/camera/pmFPAWrite.h =================================================================== --- trunk/psModules/src/camera/pmFPAWrite.h (revision 7017) +++ trunk/psModules/src/camera/pmFPAWrite.h (revision 7017) @@ -0,0 +1,28 @@ +#ifndef PM_FPA_WRITE_H +#define PM_FPA_WRITE_H + +#include "pslib.h" +#include "pmFPA.h" + +bool pmReadoutWriteNext(pmReadout *readout, // Readout to write + psFits *fits, // FITS file to which to write + int z // Image plane to write + ); + +bool pmCellWrite(pmCell *cell, // Cell to write + psFits *fits, // FITS file to which to write + psDB *db, // Database handle for "concepts" update + bool pixels // Write the pixels? + ); +bool pmChipWrite(pmChip *chip, // Chip to write + psFits *fits, // FITS file to which to write + psDB *db, // Database handle for "concepts" update + bool pixels // Write the pixels? + ); +bool pmFPAWrite(pmFPA *fpa, // FPA to write + psFits *fits, // FITS file to which to write + psDB *db, // Database handle for "concepts" update + bool pixels // Write the pixels? + ); + +#endif Index: trunk/psModules/src/camera/pmFPA_JPEG.c =================================================================== --- trunk/psModules/src/camera/pmFPA_JPEG.c (revision 7017) +++ trunk/psModules/src/camera/pmFPA_JPEG.c (revision 7017) @@ -0,0 +1,163 @@ +/** @file pmFPA_JPEG.c + * + * This file contains functions to write JPEG images. + * + * @author EAM, IfA + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-05-01 01:55:43 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +/*****************************************************************************/ +/* INCLUDE FILES */ +/*****************************************************************************/ + +#include +# include "psImageJpeg.h" +# include "psImageUnbin.h" +# include "psAdditionals.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAfile.h" +#include "pmFPAview.h" +#include "pmFPA_JPEG.h" + + +bool pmFPAviewWriteJPEG (const pmFPAview *view, pmFPAfile *file) +{ + + pmFPA *fpa = file->fpa; + + if (view->chip == -1) { + pmFPAWriteJPEG (fpa, view, file); + return true; + } + + if (view->chip >= fpa->chips->n) { + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + pmChipWriteJPEG (chip, view, file); + return true; + } + + if (view->cell >= chip->cells->n) { + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + pmCellWriteJPEG (cell, view, file); + return true; + } + + if (view->readout >= cell->readouts->n) { + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + pmReadoutWriteJPEG (readout, view, file); + return true; +} + +// read in all chip-level JPEG files for this FPA +bool pmFPAWriteJPEG (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < fpa->chips->n; i++) { + + pmChip *chip = fpa->chips->data[i]; + pmChipWriteJPEG (chip, view, file); + } + return true; +} + +// read in all cell-level JPEG files for this chip +bool pmChipWriteJPEG (pmChip *chip, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < chip->cells->n; i++) { + + pmCell *cell = chip->cells->data[i]; + pmCellWriteJPEG (cell, view, file); + } + return true; +} + +// read in all readout-level JPEG files for this cell +bool pmCellWriteJPEG (pmCell *cell, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < cell->readouts->n; i++) { + + pmReadout *readout = cell->readouts->data[i]; + pmReadoutWriteJPEG (readout, view, file); + } + return true; +} + +// read in all readout-level Objects files for this cell +bool pmReadoutWriteJPEG (pmReadout *readout, const pmFPAview *view, pmFPAfile *file) +{ + char *name, *mode, *word, *mapname; + psArray *range; + psImageJpegColormap *map; + + switch (file->type) { + case PM_FPA_FILE_JPEG: + + name = pmFPAfileNameFromRule (file->filerule, file, view); + + mapname = pmFPAfileNameFromRule (file->filextra, file, view); + map = psImageJpegColormapSet (NULL, mapname); + if (!map) { + map = psImageJpegColormapSet (NULL, "-greyscale"); + } + + mode = pmFPAfileNameFromRule (file->extrule, file, view); + word = pmFPAfileNameFromRule (file->extxtra, file, view); + range = psStringSplitArray (word, ":", false); + + // XXX we need to decide where the scale will come from in the long term + psImageClippedStatsInit (10000); + + psStats *stats = psImageClippedStats (readout->image, NULL, 0, 0.25, 0.75); + float delta = stats->robustUQ - stats->robustLQ; + float min = stats->robustMedian - 3*delta; + float max = stats->robustMedian - 6*delta; + if (!strcasecmp (mode, "RANGE") && range) { + float fmin = atof(range->data[0]); + float fmax = atof(range->data[1]); + min = stats->robustMedian + fmin*delta; + max = stats->robustMedian + fmax*delta; + } + if (!strcasecmp (mode, "FRACTION")) { + float fmin = atof(range->data[0]); + float fmax = atof(range->data[1]); + min = fmin*stats->robustMedian; + max = fmax*stats->robustMedian; + } + psImageJpeg (map, readout->image, name, min, max); + + psFree (name); + psFree (mode); + psFree (word); + psFree (mapname); + psFree (map); + psFree (stats); + psImageClippedStatsCleanup (); + + return true; + + default: + fprintf (stderr, "warning: type mismatch\n"); + break; + } + return false; +} Index: trunk/psModules/src/camera/pmFPA_JPEG.h =================================================================== --- trunk/psModules/src/camera/pmFPA_JPEG.h (revision 7017) +++ trunk/psModules/src/camera/pmFPA_JPEG.h (revision 7017) @@ -0,0 +1,27 @@ +/** @file pmFPAview.h +* +* @brief Tools to manipulate the FPA structure elements. +* +* @ingroup AstroImage +* +* @author EAM, IfA +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-05-01 01:55:43 $ +* +* Copyright 2004-2005 Institute for Astronomy, University of Hawaii +*/ + +#ifndef PM_FPA_JPEG_H +#define PM_FPA_JPEG_H + +/// @addtogroup AstroImage +/// @{ + +bool pmFPAviewWriteJPEG (const pmFPAview *view, pmFPAfile *file); +bool pmFPAWriteJPEG (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file); +bool pmChipWriteJPEG (pmChip *chip, const pmFPAview *view, pmFPAfile *file); +bool pmCellWriteJPEG (pmCell *cell, const pmFPAview *view, pmFPAfile *file); +bool pmReadoutWriteJPEG (pmReadout *readout, const pmFPAview *view, pmFPAfile *file); + +# endif Index: trunk/psModules/src/camera/pmFPAfile.c =================================================================== --- trunk/psModules/src/camera/pmFPAfile.c (revision 7017) +++ trunk/psModules/src/camera/pmFPAfile.c (revision 7017) @@ -0,0 +1,1144 @@ +#include +#include "pslib.h" +#include "psAdditionals.h" +#include "pmConfig.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmMaskBadPixels.h" +#include "pmFPAConstruct.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" +#include "pmFPACopy.h" +#include "pmFPARead.h" +#include "pmFPAWrite.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceIO.h" +#include "pmGrowthCurve.h" +#include "pmPSF.h" +#include "pmPSF_IO.h" +#include "pmFPA_JPEG.h" + +static void pmFPAfileFree (pmFPAfile *file) +{ + + if (file == NULL) + return; + + psFree (file->fpa); + psFree (file->readout); + psFree (file->names); + + psFree (file->format); + psFree (file->name); + + if (file->fits != NULL) { + psFitsClose (file->fits); + } + + psFree (file->filerule); + psFree (file->filextra); + psFree (file->extrule); + psFree (file->extxtra); + + psFree (file->filename); + psFree (file->extname); + + // these are just views ?? + // psFree (file->phu); + // psFree (file->header); + + return; +} + +pmFPAfile *pmFPAfileAlloc () +{ + + pmFPAfile *file = psAlloc (sizeof(pmFPAfile)); + psMemSetDeallocator (file, (psFreeFunc) pmFPAfileFree); + + file->phu = NULL; + file->readout = NULL; + file->header = NULL; + + file->fpa = NULL; + file->fits = NULL; + file->names = psMetadataAlloc(); + + file->format = NULL; + file->name = NULL; + + file->filerule = NULL; + file->filextra = NULL; + file->extrule = NULL; + file->extxtra = NULL; + + file->filename = NULL; + file->extname = NULL; + + file->type = PM_FPA_FILE_NONE; + file->mode = PM_FPA_MODE_NONE; + file->state = PM_FPA_STATE_CLOSED; + + return (file); +} + +pmFPAfile *pmFPAfileDefine (psMetadata *files, psMetadata *camera, pmFPA *fpa, char *name) +{ + + bool status; + char *depth; + char *type; + + // select the FILERULES from the camera config + psMetadata *filerules = psMetadataLookupPtr (&status, camera, "FILERULES"); + if (filerules == NULL) { + psErrorStackPrint(stderr, "Can't find FILERULES in the CAMERA configuration!\n"); + return NULL; + } + + // select the name from the FILERULES + // check for alias name (type == STR, name is aliased name) + char *realname = psMetadataLookupStr (&status, filerules, name); + if (realname == NULL) + realname = name; + + psMetadata *data = psMetadataLookupPtr (&status, filerules, realname); + if (data == NULL) { + psErrorStackPrint(stderr, "Can't find file concept %s!\n", name); + return NULL; + } + + pmFPAfile *file = pmFPAfileAlloc (); + + // save the name of this pmFPAfile + file->name = psStringCopy (name); + + file->filerule = psMemIncrRefCounter(psMetadataLookupStr (&status, data, "FILENAME.RULE")); + file->filextra = psMemIncrRefCounter(psMetadataLookupStr (&status, data, "FILENAME.XTRA")); + file->extrule = psMemIncrRefCounter(psMetadataLookupStr (&status, data, "EXTNAME.RULE")); + file->extxtra = psMemIncrRefCounter(psMetadataLookupStr (&status, data, "EXTNAME.XTRA")); + + file->fileDepth = PM_FPA_DEPTH_NONE; + depth = psMetadataLookupStr (&status, data, "FILE.DEPTH"); + if (depth != NULL) { + if (!strcasecmp (depth, "FPA")) { + file->fileDepth = PM_FPA_DEPTH_FPA; + } + if (!strcasecmp (depth, "CHIP")) { + file->fileDepth = PM_FPA_DEPTH_CHIP; + } + if (!strcasecmp (depth, "CELL")) { + file->fileDepth = PM_FPA_DEPTH_CELL; + } + if (!strcasecmp (depth, "READOUT")) { + file->fileDepth = PM_FPA_DEPTH_READOUT; + } + } + if (file->fileDepth == PM_FPA_DEPTH_NONE) { + psLogMsg (__func__, 3, "warning: FILE.DEPTH is not set for %s\n", name); + } + + file->dataDepth = PM_FPA_DEPTH_NONE; + depth = psMetadataLookupStr (&status, data, "DATA.DEPTH"); + if (depth != NULL) { + if (!strcasecmp (depth, "FPA")) { + file->dataDepth = PM_FPA_DEPTH_FPA; + } + if (!strcasecmp (depth, "CHIP")) { + file->dataDepth = PM_FPA_DEPTH_CHIP; + } + if (!strcasecmp (depth, "CELL")) { + file->dataDepth = PM_FPA_DEPTH_CELL; + } + if (!strcasecmp (depth, "READOUT")) { + file->dataDepth = PM_FPA_DEPTH_READOUT; + } + } + if (file->dataDepth == PM_FPA_DEPTH_NONE) { + psLogMsg (__func__, 3, "warning: DATA.DEPTH is not set for %s\n", name); + } + + file->type = PM_FPA_FILE_NONE; + type = psMetadataLookupStr (&status, data, "FILE.TYPE"); + if (type != NULL) { + if (!strcasecmp (type, "SX")) { + file->type = PM_FPA_FILE_SX; + } + if (!strcasecmp (type, "OBJ")) { + file->type = PM_FPA_FILE_OBJ; + } + if (!strcasecmp (type, "CMP")) { + file->type = PM_FPA_FILE_CMP; + } + if (!strcasecmp (type, "CMF")) { + file->type = PM_FPA_FILE_CMF; + } + if (!strcasecmp (type, "RAW")) { + file->type = PM_FPA_FILE_RAW; + } + if (!strcasecmp (type, "IMAGE")) { + file->type = PM_FPA_FILE_IMAGE; + } + if (!strcasecmp (type, "PSF")) { + file->type = PM_FPA_FILE_PSF; + } + if (!strcasecmp (type, "JPEG")) { + file->type = PM_FPA_FILE_JPEG; + } + } + if (file->type == PM_FPA_FILE_NONE) { + psLogMsg (__func__, 3, "warning: FILE.TYPE is not defined for %s\n", name); + } + + file->mode = PM_FPA_MODE_NONE; + type = psMetadataLookupStr (&status, data, "FILE.MODE"); + if (type != NULL) { + if (!strcasecmp (type, "READ")) { + file->mode = PM_FPA_MODE_READ; + } + if (!strcasecmp (type, "WRITE")) { + file->mode = PM_FPA_MODE_WRITE; + } + } + if (file->mode == PM_FPA_MODE_NONE) { + psLogMsg (__func__, 3, "warning: FILE.MODE is not defined for %s\n", name); + } + + if (fpa != NULL) { + file->fpa = psMemIncrRefCounter(fpa); + } + + // for WRITE type of data, the output format needs to be determined + // this is only needed if the output file is not identical in structure to the input + char *formatName = psMetadataLookupStr (&status, data, "FILE.FORMAT"); + if (formatName && strcasecmp (formatName, "NONE")) { + psMetadata *formats = psMetadataLookupMD(&status, camera, "FORMATS"); // List of formats + char *formatFile = psMetadataLookupStr (&status, formats, formatName); + readConfig (&file->format, formatFile, formatName); + } + + psMetadataAddPtr (files, PS_LIST_TAIL, name, PS_DATA_UNKNOWN, "", file); + + // we free this copy of file, but 'files' still has a copy + psFree (file); + + // the returned value is a view into the version on 'files' + return (file); +} + +/* +pmFPAfile *pmFPAfileConstruct (psMetadata *files, psMetadata *format, psMetadata *camera, char *name) +{ + pmFPA *fpa = pmFPAConstruct (camera); + pmFPAfile *file = pmFPAfileDefine (files, format, fpa, name); + psFree (fpa); + return file; +} +*/ + +// open file (if not already opened) +bool pmFPAfileOpen (pmFPAfile *file, const pmFPAview *view) +{ + + bool status; + char *extra; + char *mode = NULL; + char *readMode = "r"; + char *writeMode = "w"; + + if (file->state & PM_FPA_STATE_INACTIVE) { + return false; + } + + if (file->state == PM_FPA_STATE_OPEN) { + return true; + } + + if (file->mode == PM_FPA_MODE_NONE) { + return false; + } + if (file->mode == PM_FPA_MODE_INTERNAL) { + return false; + } + if (file->mode == PM_FPA_MODE_READ) { + mode = readMode; + } + if (file->mode == PM_FPA_MODE_WRITE) { + mode = writeMode; + } + + // determine the file name + file->filename = pmFPAfileNameFromRule (file->filerule, file, view); + if (file->filename == NULL) + return false; + + // indirect filenames + if (!strcasecmp (file->filename, "@FILES")) { + psFree (file->filename); + extra = pmFPAfileNameFromRule (file->filextra, file, view); + file->filename = psMetadataLookupStr (&status, file->names, extra); + psFree (extra); + if (file->filename == NULL) + return false; + // psMetadataLookupStr just returns a view, file->filename must be protected + psMemIncrRefCounter (file->filename); + } + if (!strcasecmp (file->filename, "@DETDB")) { + psFree (file->filename); + extra = pmFPAfileNameFromRule (file->filextra, file, view); + // file->filename = pmDetrendSelect (extra); + psFree (extra); + if (file->filename == NULL) + return false; + psMemIncrRefCounter (file->filename); + } + + switch (file->type) { + // open the FITS types: + case PM_FPA_FILE_IMAGE: + case PM_FPA_FILE_CMF: + psTrace ("pmFPAfile", 5, "opening %s (type: %d)\n", file->filename, file->type); + file->fits = psFitsOpen (file->filename, mode); + if (file->fits == NULL) + psAbort (__func__, "error opening file %s\n", file->filename); + file->state = PM_FPA_STATE_OPEN; + break; + + // defer opening TEXT types: + case PM_FPA_FILE_SX: + case PM_FPA_FILE_OBJ: + case PM_FPA_FILE_CMP: + case PM_FPA_FILE_RAW: + case PM_FPA_FILE_PSF: + case PM_FPA_FILE_JPEG: + psTrace ("pmFPAfile", 5, "not opening %s (type: %d)\n", file->filename, file->type); + break; + + default: + fprintf (stderr, "warning: type mismatch for %s : %d\n", file->filename, file->type); + return false; + } + return true; +} + +bool pmFPAfileRead (pmFPAfile *file, const pmFPAview *view) +{ + if (file->state & PM_FPA_STATE_INACTIVE) + return false; + + if (file->mode != PM_FPA_MODE_READ) + return false; + + // get the current depth + pmFPAdepth depth = pmFPAviewDepth (view); + + // do we need to open this file? + if (depth == file->fileDepth) { + pmFPAfileOpen (file, view); + } + + // do we need to read this file? + if (depth != file->dataDepth) + return false; + + switch (file->type) { + case PM_FPA_FILE_IMAGE: + if (pmFPAviewReadFitsImage (view, file)) { + psTrace ("pmFPAfile", 5, "reading %s (type: %d)\n", file->filename, file->type); + } else { + psTrace ("pmFPAfile", 5, "skipping %s (type: %d)\n", file->filename, file->type); + } + break; + + case PM_FPA_FILE_SX: + case PM_FPA_FILE_RAW: + case PM_FPA_FILE_OBJ: + case PM_FPA_FILE_CMP: + case PM_FPA_FILE_CMF: + pmFPAviewReadObjects (view, file); + psTrace ("pmFPAfile", 5, "reading %s (type: %d)\n", file->filename, file->type); + break; + + case PM_FPA_FILE_PSF: + pmFPAviewReadPSFmodel (view, file); + psTrace ("pmFPAfile", 5, "reading %s (type: %d)\n", file->filename, file->type); + break; + + case PM_FPA_FILE_JPEG: + break; + + default: + fprintf (stderr, "warning: type mismatch\n"); + return false; + } + return true; +} + +bool pmFPAfileWrite (pmFPAfile *file, const pmFPAview *view) +{ + + if (file->state & PM_FPA_STATE_INACTIVE) { + return false; + } + + if (file->mode != PM_FPA_MODE_WRITE) { + return false; + } + + // get the current depth + pmFPAdepth depth = pmFPAviewDepth (view); + + // do we need to write this file? + if (depth != file->dataDepth) + return false; + + // do we need to open this file? + if (depth >= file->fileDepth) { + pmFPAfileOpen (file, view); + } + + switch (file->type) { + case PM_FPA_FILE_IMAGE: + pmFPAviewWriteFitsImage (view, file); + psTrace ("pmFPAfile", 5, "wrote image %s (fpa: %p)\n", file->filename, file->fpa); + break; + + case PM_FPA_FILE_SX: + case PM_FPA_FILE_RAW: + case PM_FPA_FILE_OBJ: + case PM_FPA_FILE_CMP: + case PM_FPA_FILE_CMF: + pmFPAviewWriteObjects (view, file); + psTrace ("pmFPAfile", 5, "wrote object %s (fpa: %p)\n", file->filename, file->fpa); + break; + + case PM_FPA_FILE_PSF: + pmFPAviewWritePSFmodel (view, file); + psTrace ("pmFPAfile", 5, "wrote PSF %s (fpa: %p)\n", file->filename, file->fpa); + break; + + case PM_FPA_FILE_JPEG: + pmFPAviewWriteJPEG (view, file); + psTrace ("pmFPAfile", 5, "wrote PSF %s (fpa: %p)\n", file->filename, file->fpa); + break; + + default: + fprintf (stderr, "warning: type mismatch\n"); + return false; + } + return true; +} + +// create the data elements (headers, images) appropriate for this view +bool pmFPAfileCreate (pmFPAfile *file, const pmFPAview *view) +{ + if (file->state & PM_FPA_STATE_INACTIVE) { + return false; + } + if (file->mode != PM_FPA_MODE_WRITE) + return false; + + // get the current depth + pmFPAdepth depth = pmFPAviewDepth (view); + + // do we need to write this file? + if (depth != file->dataDepth) + return false; + + // XXX is this a sufficient check to avoid creating elements? + if (file->src == NULL) + return false; + + switch (file->type) { + case PM_FPA_FILE_IMAGE: + /* create a PHU for thie file, if it does not exist */ + pmFPAfileCopyStructureView (file->fpa, file->src, file->format, file->xBin, file->yBin, view); + psTrace ("pmFPAfile", 5, "created fpa data elements for %s (fpa: %p)\n", file->filename, file->fpa); + break; + + case PM_FPA_FILE_SX: + case PM_FPA_FILE_RAW: + case PM_FPA_FILE_OBJ: + case PM_FPA_FILE_CMP: + case PM_FPA_FILE_CMF: + case PM_FPA_FILE_PSF: + case PM_FPA_FILE_JPEG: + break; + + default: + fprintf (stderr, "warning: type mismatch\n"); + return false; + } + return true; +} + +bool pmFPAfileClose (pmFPAfile *file, const pmFPAview *view) +{ + if (file->state & PM_FPA_STATE_INACTIVE) { + return false; + } + if (file->state == PM_FPA_STATE_CLOSED) { + return true; + } + + // is current depth == open depth? + pmFPAdepth depth = pmFPAviewDepth (view); + if (file->fileDepth != depth) { + return false; + } + + // check if we are actually open + switch (file->type) { + // check the FITS types + case PM_FPA_FILE_IMAGE: + case PM_FPA_FILE_CMF: + psFitsClose (file->fits); + psTrace ("pmFPAfile", 5, "closing %s (type: %d)\n", file->filename, file->type); + file->fits = NULL; + file->phu = NULL; + file->header = NULL; + file->state = PM_FPA_STATE_CLOSED; + break; + + // ignore the TEXT types + case PM_FPA_FILE_SX: + case PM_FPA_FILE_RAW: + case PM_FPA_FILE_OBJ: + case PM_FPA_FILE_CMP: + case PM_FPA_FILE_PSF: + case PM_FPA_FILE_JPEG: + break; + + default: + fprintf (stderr, "warning: type mismatch\n"); + return false; + } + return true; +} + +// set the state of the specified pmFPAfile to active (state == true) or inactive +// if name is NULL, set the state for all pmFPAfiles +bool pmFPAfileActivate (psMetadata *files, bool state, char *name) +{ + if (name == NULL) { + psMetadataItem *item = NULL; + psMetadataIterator *iter = psMetadataIteratorAlloc (files, PS_LIST_HEAD, NULL); + while ((item = psMetadataGetAndIncrement (iter)) != NULL) { + pmFPAfile *file = item->data.V; + if (state) { + file->state &= NOT_U8(PM_FPA_STATE_INACTIVE); + } else { + file->state |= PM_FPA_STATE_INACTIVE; + } + } + psFree (iter); + return true; + } + + bool status = false; + pmFPAfile *file = psMetadataLookupPtr (&status, files, name); + if (!file) { + return false; + } + if (state) { + file->state &= NOT_U8(PM_FPA_STATE_INACTIVE); + } else { + file->state |= PM_FPA_STATE_INACTIVE; + } + return true; +} + +// attempt open, read, write, or close pmFPAfiles in files +bool pmFPAfileIOChecks (psMetadata *files, const pmFPAview *view, pmFPAfilePlace place) +{ + // recipe override values (command-line options): + psMetadataItem *item = NULL; + psMetadataIterator *iter = psMetadataIteratorAlloc (files, PS_LIST_HEAD, NULL); + while ((item = psMetadataGetAndIncrement (iter)) != NULL) { + pmFPAfile *file = item->data.V; + + if (place == PM_FPA_BEFORE) { + pmFPAfileRead (file, view); + pmFPAfileCreate (file, view); + } else { + pmFPAfileWrite (file, view); + pmFPAfileClose (file, view); + } + } + psFree (iter); + return true; +} + +// create a file with the given name, assign it type "INTERNAL", and supply it with an image +// of the requested dimensions. (image only, mask and weight are ignored) +pmReadout *pmFPAfileCreateInternal (psMetadata *files, char *name, int Nx, int Ny, int type) +{ + pmReadout *readout = pmReadoutAlloc (NULL); + readout->image = psImageAlloc (Nx, Ny, type); + + // I want an image from the + pmFPAfile *file = pmFPAfileAlloc(); + file->mode = PM_FPA_MODE_INTERNAL; + + file->readout = readout; + psMetadataAddPtr (files, PS_LIST_TAIL, name, PS_DATA_UNKNOWN, "", file); + psFree (file); + // we free this copy of file, but 'files' still has a copy + + return (readout); +} + +bool pmFPAfileDropInternal (psMetadata *files, char *name) +{ + bool status; + + pmFPAfile *file = psMetadataLookupPtr (&status, files, name); + if (file == NULL) + return false; + + if (file->mode != PM_FPA_MODE_INTERNAL) + return false; + + psMetadataRemoveKey (files, name); + return true; +} + +// select the readout from the named pmFPAfile; if the named file does not exist, +pmReadout *pmFPAfileThisReadout (psMetadata *files, const pmFPAview *view, const char *name) +{ + bool status; + + pmFPAfile *file = psMetadataLookupPtr (&status, files, name); + if (file == NULL) + return NULL; + + // internal files have the readout as a separate element: + if (file->mode == PM_FPA_MODE_INTERNAL) { + return file->readout; + } + + pmReadout *readout = pmFPAviewThisReadout (view, file->fpa); + return readout; +} + +// given an already-opened fits file, read the components corresponding +// to the specified view +bool pmFPAviewReadFitsImage (const pmFPAview *view, pmFPAfile *file) +{ + bool status; + pmFPA *fpa = file->fpa; + psFits *fits = file->fits; + + if (view->chip == -1) { + status = pmFPARead (fpa, fits, NULL); + return status; + } + + if (view->chip >= fpa->chips->n) { + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + status = pmChipRead (chip, fits, NULL); + return status; + } + + if (view->cell >= chip->cells->n) { + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + status = pmCellRead (cell, fits, NULL); + return status; + } + return false; + + // XXX pmReadoutRead, pmReadoutReadSegement disabled for now + # if (0) + + if (view->readout >= cell->readouts->n) { + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + if (view->nRows == 0) { + pmReadoutRead (readout, fits, NULL); + } else { + pmReadoutReadSegment (readout, fits, view->nRows, view->iRows, NULL, NULL); + } + return true; + # endif +} + +// given an already-opened fits file, write the components corresponding +// to the specified view +bool pmFPAviewWriteFitsImage (const pmFPAview *view, pmFPAfile *file) +{ + + pmFPA *fpa = file->fpa; + psFits *fits = file->fits; + + // pmFPAWrite takes care of all PHUs as needed + if (view->chip == -1) { + pmFPAWrite (fpa, fits, NULL, true); + return true; + } + + if (view->chip >= fpa->chips->n) { + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + // do we need to write out a PHU for this entry? + if (file->phu == NULL) { + pmHDU *hdu = pmFPAviewThisHDU (view, file->fpa); + pmHDU *phu = pmFPAviewThisPHU (view, file->fpa); + // if this hdu is the phu, the write function below will create the phu + if (hdu != phu) { + // we assume that the PHU is just a header + psMetadata *outhead = psMetadataCopy (NULL, phu->header); + psMetadataAdd (outhead, PS_LIST_TAIL, "EXTEND", PS_DATA_BOOL | PS_META_REPLACE, "this file has extensions", true); + psFitsWriteHeaderNotImage (file->fits, outhead); + file->phu = phu->header; + psTrace ("pmFPAfile", 5, "wrote phu %s (type: %d)\n", file->filename, file->type); + psFree (outhead); + } + } + + if (view->cell == -1) { + pmChipWrite (chip, fits, NULL, true); + return true; + } + + if (view->cell >= chip->cells->n) { + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + pmCellWrite (cell, fits, NULL, true); + return true; + } + return false; + + // XXX disable readout write for now + # if (0) + + if (view->readout >= cell->readouts->n) { + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + if (view->nRows == 0) { + pmReadoutWrite (readout, fits, NULL, NULL); + } else { + pmReadoutWriteSegment (readout, fits, view->nRows, view->iRows, NULL, NULL); + } + return true; + # endif +} + +// look for the given name on the argument list. +// returns the file (a view to the one saved on config->files) +pmFPAfile *pmFPAfileFromArgs (bool *found, pmConfig *config, char *filename, char *argname) +{ + bool status; + pmFPA *fpa = NULL; + psFits *fits = NULL; + pmFPAfile *file = NULL; + psMetadata *phu = NULL; + psMetadata *format = NULL; + + if (*found) + return NULL; + + // we search the argument data for the named fileset (argname) + psArray *infiles = psMetadataLookupPtr(&status, config->arguments, argname); + if (!status) + return NULL; + if (infiles->n < 1) + return NULL; + + // determine the current format from the header + // if no camera has been specified, use the first image as a template for the rest. + fits = psFitsOpen (infiles->data[0], "r"); + phu = psFitsReadHeader (NULL, fits); + format = pmConfigCameraFormatFromHeader (config, phu); + psFitsClose (fits); + + // build the template fpa, set up the basic view + fpa = pmFPAConstruct (config->camera); + + // load the given filerule (from config->camera) and associate it with the fpa + // the output file is just a view to the file on config->files + file = pmFPAfileDefine (config->files, config->camera, fpa, filename); + if (!file) { + psErrorStackPrint(stderr, "file %s not defined\n", filename); + psFree (fpa); + psFree (format); + return NULL; + } + + // this file is (by virtue of being supplied in the argument list) coming from the fileset + psFree (file->filerule); + psFree (file->filextra); + + // adjust the rules to identify these files in the file->names data + file->filerule = psStringCopy ("@FILES"); + // XXX this rule does not work in general + file->filextra = psStringCopy ("{CHIP.NAME}.{CELL.NAME}"); + + // examine the list of input files and validate their cameras + for (int i = 0; i < infiles->n; i++) { + if (i > 0) { + fits = psFitsOpen (infiles->data[i], "r"); + phu = psFitsReadHeader (NULL, fits); + pmConfigValidateCameraFormat (format, phu); + psFitsClose (fits); + } + + // set the view to the corresponding entry for this phu + pmFPAview *view = pmFPAAddSourceFromHeader (fpa, phu, format); + + // XXX is this the correct psMD to save the filename? + char *name = pmFPAfileNameFromRule (file->filextra, file, view); + psMetadataAddStr (file->names, PS_LIST_TAIL, name, 0, "", infiles->data[i]); + + psFree (view); + psFree (name); + psFree (phu); + } + psFree (fpa); + psFree (format); + *found = true; + + return file; +} + +// XXX this this function through, then finish +# if 0 +// look for the given name on the argument list. +// returns the file (a view to the one saved on config->files) +// in this case, each file should correspond to the same view +// (except at the readout level), of multiple FPAs +// XXX think this through a bit more: +// - do we create multiple input fpas, or assign the same files to the single fpa? +// - do we save the multiple pmFPAfiles in the config->files psMD? +// - do we save the multiple filenames in the file->names psMD? +// - if we have a single FPA and multiple names, we need a method to +// turn on a specific name for the I/O actions. probably true if we +// have multiple FPAs as well. +pmFPAfile *pmFPAfileSetFromArgs (bool *found, pmConfig *config, char *filename, char *argname) +{ + bool status; + pmFPA *fpa = NULL; + psFits *fits = NULL; + pmFPAfile *file = NULL; + psMetadata *phu = NULL; + psMetadata *format = NULL; + + if (*found) + return NULL; + + // we search the argument data for the named fileset (argname) + psArray *infiles = psMetadataLookupPtr(&status, config->arguments, argname); + if (!status) + return NULL; + if (infiles->n < 1) + return NULL; + + // determine the current format from the header + // if no camera has been specified, use the first image as a template for the rest. + fits = psFitsOpen (infiles->data[0], "r"); + phu = psFitsReadHeader (NULL, fits); + format = pmConfigCameraFormatFromHeader (config, phu); + psFitsClose (fits); + + // build the template fpa, set up the basic view + fpa = pmFPAConstruct (config->camera); + + // load the given filerule (from config->camera) and associate it with the fpa + // the output file is just a view to the file on config->files + // XXX the filenames placed on config->files should include the seq number + file = pmFPAfileDefine (config->files, config->camera, fpa, filename); + if (!file) { + psErrorStackPrint(stderr, "file %s not defined\n", filename); + psFree (fpa); + psFree (format); + return NULL; + } + + // this file is (by virtue of being supplied in the argument list) coming from the fileset + psFree (file->filerule); + psFree (file->filextra); + + // adjust the rules to identify these files in the file->names data + file->filerule = psStringCopy ("@FILES"); + // XXX this rule does not work in general + file->filextra = psStringCopy ("{CHIP.NAME}.{CELL.NAME}"); + + // examine the list of input files and validate their cameras + for (int i = 0; i < infiles->n; i++) { + if (i > 0) { + fits = psFitsOpen (infiles->data[i], "r"); + phu = psFitsReadHeader (NULL, fits); + pmConfigValidateCameraFormat (format, phu); + psFitsClose (fits); + } + + // set the view to the corresponding entry for this phu + pmFPAview *view = pmFPAAddSourceFromHeader (fpa, phu, format); + + // XXX is this the correct psMD to save the filename? + char *name = pmFPAfileNameFromRule (file->filextra, file, view); + psMetadataAddStr (file->names, PS_LIST_TAIL, name, 0, "", infiles->data[i]); + + psFree (view); + psFree (name); + psFree (phu); + } + psFree (fpa); + psFree (format); + *found = true; + + return file; +} +# endif + +// create a new output pmFPAfile based on an existing FPA +pmFPAfile *pmFPAfileFromFPA (pmConfig *config, pmFPA *src, int xBin, int yBin, char *filename) +{ + // XXX pmFPAConstruct has many leaks (6919) + pmFPA *fpa = pmFPAConstruct (config->camera); + pmFPAfile *file = pmFPAfileDefine (config->files, config->camera, fpa, filename); + file->src = src; // inherit output elements from this source pmFPA + file->xBin = xBin; + file->yBin = yBin; + + return file; +} + +// look for the given name on the argument list. +pmFPAfile *pmFPAfileFromConf (bool *found, pmConfig *config, char *filename, pmFPA *input) +{ + psFits *fits = NULL; + pmFPAfile *file = NULL; + psMetadata *phu = NULL; + psMetadata *format = NULL; + psArray *infiles = NULL; + + if (*found) + return NULL; + + // a camera config is needed (as source of file rule) + if (config->camera == NULL) { + psErrorStackPrint (stderr, "camera is not defined\n"); + return NULL; + } + + // expect @DETDB in the config filerule + file = pmFPAfileDefine (config->files, config->camera, NULL, filename); + if (!file) { + psErrorStackPrint(stderr, "file %s not defined\n", filename); + return NULL; + } + + // image names come from the file->name list? + if (!strcasecmp (file->filerule, "@FILES")) + psAbort ("pmFPAfileFromConfig", "programming error"); + + // image needs to come from the detrend database + if (!strcasecmp (file->filerule, "@DETDB")) { + // char *extra = pmFPAfileNameFromRule (file->filextra, file, view); + // psArray *infiles = pmDetrendSelect (extra, input); + psAbort ("pmFPAfileFromConfig", "@DETDB not yet defined"); + } else { + infiles = psArrayAlloc(1); + infiles->n = 1; + infiles->data[0] = psStringCopy (file->filerule); + } + if (infiles == NULL) + return NULL; + if (infiles->n < 1) { + psFree (infiles); + return NULL; + } + + // determine the current format from the header + // if no camera has been specified, use the first image as a template for the rest. + fits = psFitsOpen (infiles->data[0], "r"); + phu = psFitsReadHeader (NULL, fits); + format = pmConfigCameraFormatFromHeader (config, phu); + psFitsClose (fits); + + // build the template fpa, set up the basic view + file->fpa = pmFPAConstruct (config->camera); + + // this file is (by virtue of being supplied in the argument list) coming from the fileset + psFree (file->filerule); + psFree (file->filextra); + + // adjust the rules to identify these files in the file->names data + file->filerule = psStringCopy ("@FILES"); + // XXX this rule does not work in general + file->filextra = psStringCopy ("{CHIP.NAME}.{CELL.NAME}"); + + // examine the list of input files and validate their cameras + for (int i = 0; i < infiles->n; i++) { + if (i > 0) { + fits = psFitsOpen (infiles->data[i], "r"); + phu = psFitsReadHeader (NULL, fits); + pmConfigValidateCameraFormat (format, phu); + psFitsClose (fits); + } + + // set the view to the corresponding entry for this phu + pmFPAview *view = pmFPAAddSourceFromHeader (file->fpa, phu, format); + + // XXX is this the correct psMD to save the filename? + char *name = pmFPAfileNameFromRule (file->filextra, file, view); + psMetadataAddStr (file->names, PS_LIST_TAIL, name, 0, "", infiles->data[i]); + + psFree (view); + psFree (name); + psFree (phu); + } + psFree (format); + psFree (infiles); + *found = true; + return file; +} + +bool pmFPAfileAddFileNames (psMetadata *files, char *name, char *value, int mode) +{ + + // add the output names to the output-type files + psMetadataItem *item = NULL; + psMetadataIterator *iter = psMetadataIteratorAlloc (files, PS_LIST_HEAD, NULL); + while ((item = psMetadataGetAndIncrement (iter)) != NULL) { + pmFPAfile *file = item->data.V; + + if (file->mode == mode) { + psMetadataAddStr (file->names, PS_LIST_TAIL, name, 0, "", value); + } + } + psFree (iter); + return true; +} + +// select the rule from the camera configuration, perform substitutions as needed +char *pmFPAfileNameFromRule (char *rule, pmFPAfile *file, const pmFPAview *view) +{ + + char *newName = NULL; // destination for resulting name + + newName = psStringCopy (rule); + + if (strstr (newName, "{CHIP.NAME}") != NULL) { + pmChip *chip = pmFPAviewThisChip (view, file->fpa); + if (chip != NULL) { + char *name = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); + if (name != NULL) { + newName = psStringSubstitute (newName, name, "{CHIP.NAME}"); + } + } + } + if (strstr (newName, "{CELL.NAME}") != NULL) { + pmCell *cell = pmFPAviewThisCell (view, file->fpa); + if (cell != NULL) { + char *name = psMetadataLookupStr (NULL, cell->concepts, "CELL.NAME"); + if (name != NULL) { + newName = psStringSubstitute (newName, name, "{CELL.NAME}"); + } + } + } + if (strstr (newName, "{EXTNAME}") != NULL) { + pmHDU *hdu = pmFPAviewThisHDU (view, file->fpa); + if (hdu->extname != NULL) { + newName = psStringSubstitute (newName, hdu->extname, "{EXTNAME}"); + } + } + if (strstr (newName, "{OUTPUT}") != NULL) { + char *name = psMetadataLookupStr (NULL, file->names, "OUTPUT"); + if (name != NULL) { + newName = psStringSubstitute (newName, name, "{OUTPUT}"); + } + } + return newName; +} + +// given an already-opened fits file, write the components corresponding +// to the specified view +bool pmFPAfileCopyView (pmFPA *out, pmFPA *in, const pmFPAview *view) +{ + // pmFPAWrite takes care of all PHUs as needed + if (view->chip == -1) { + pmFPACopy (out, in); + return true; + } + if (view->chip >= in->chips->n) { + return false; + } + pmChip *inChip = in->chips->data[view->chip]; + pmChip *outChip = out->chips->data[view->chip]; + + if (view->cell == -1) { + pmChipCopy (outChip, inChip); + return true; + } + if (view->cell >= inChip->cells->n) { + return false; + } + pmCell *inCell = inChip->cells->data[view->cell]; + pmCell *outCell = outChip->cells->data[view->cell]; + + if (view->readout == -1) { + pmCellCopy (outCell, inCell); + return true; + } + return false; + + // XXX add readout / segment equivalents +} + +// given an already-opened fits file, write the components corresponding +// to the specified view +bool pmFPAfileCopyStructureView (pmFPA *out, pmFPA *in, psMetadata *format, int xBin, int yBin, const pmFPAview *view) +{ + // pmFPAWrite takes care of all PHUs as needed + if (view->chip == -1) { + pmFPAAddSourceFromView (out, view, format); + pmFPACopyStructure (out, in, xBin, yBin); + return true; + } + if (view->chip >= in->chips->n) { + return false; + } + pmChip *inChip = in->chips->data[view->chip]; + pmChip *outChip = out->chips->data[view->chip]; + + if (view->cell == -1) { + pmFPAAddSourceFromView (out, view, format); + pmChipCopyStructure (outChip, inChip, xBin, yBin); + return true; + } + if (view->cell >= inChip->cells->n) { + return false; + } + pmCell *inCell = inChip->cells->data[view->cell]; + pmCell *outCell = outChip->cells->data[view->cell]; + + if (view->readout == -1) { + pmFPAAddSourceFromView (out, view, format); + pmCellCopyStructure (outCell, inCell, xBin, yBin); + return true; + } + return false; + + // XXX add readout / segment equivalents +} Index: trunk/psModules/src/camera/pmFPAfile.h =================================================================== --- trunk/psModules/src/camera/pmFPAfile.h (revision 7017) +++ trunk/psModules/src/camera/pmFPAfile.h (revision 7017) @@ -0,0 +1,160 @@ +/** @file pmFPAview.h +* +* @brief Tools to manipulate the FPA structure elements. +* +* @ingroup AstroImage +* +* @author EAM, IfA +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-05-01 01:55:43 $ +* +* Copyright 2004-2005 Institute for Astronomy, University of Hawaii +*/ + +#ifndef PM_FPA_FILE_H +#define PM_FPA_FILE_H + +#include "pslib.h" +#include "pmConfig.h" +#include "pmFPA.h" +#include "pmFPAview.h" + +typedef enum { + PM_FPA_BEFORE, + PM_FPA_AFTER, +} pmFPAfilePlace; + +typedef enum { + PM_FPA_FILE_NONE, + PM_FPA_FILE_SX, + PM_FPA_FILE_OBJ, + PM_FPA_FILE_CMP, + PM_FPA_FILE_CMF, + PM_FPA_FILE_RAW, + PM_FPA_FILE_IMAGE, + PM_FPA_FILE_PSF, + PM_FPA_FILE_JPEG, +} pmFPAfileType; + +typedef enum { + PM_FPA_MODE_NONE, + PM_FPA_MODE_READ, + PM_FPA_MODE_WRITE, + PM_FPA_MODE_INTERNAL, +} pmFPAfileMode; + +typedef enum { + PM_FPA_STATE_OPEN = 0x01, + PM_FPA_STATE_CLOSED = 0x02, + PM_FPA_STATE_INACTIVE = 0x04, +} pmFPAfileState; + +typedef struct +{ + pmFPAfileMode mode; // is this file read, written, or only used internally? + pmFPAfileType type; // what type of data is read from / written to disk? + pmFPAfileState state; // have we opened the file, etc? + + pmFPAdepth fileDepth; // what depth in the FPA hierarchy represents a unique file? + pmFPAdepth dataDepth; // at what depth do we read/write the data segment? + + pmFPA *fpa; // for I/O files, we carry a pointer to the complete fpa + psFits *fits; // for I/O files of fits type (IMAGE, CMP, CMF), we carry a file handle + + psMetadata *phu; // pointer (view) to the current phu header + psMetadata *header; // pointer (view) to the current hdu header + + pmReadout *readout; // for internal files, we only carry a single readout + + psMetadata *names; // filenames supplied by the cmdline or detdb are saved here + + char *filerule; // rule for constructing a filename when needed + char *filextra; // additional information used to define filenames (context dependent) + char *extrule; // rule for constructing an extension name when needed + char *extxtra; // additional information used to define extension names (context dependent) + + char *name; // the name of the rule (useful for debugging / tracing) + char *filename; // the current name of an active file + char *extname; // the current name of an active file extension + + // the following elements are used for WRITE-mode IMAGE-type pmFPAfiles to inform + // the creation of a new image based on an existing image + pmFPA *src; // if an output FPA, inherit from this FPA + int xBin; // desired binning in x direction + int yBin; // desired binning in y direction + psMetadata *format; +} +pmFPAfile; + +// allocate an empty pmFPAfile structure +pmFPAfile *pmFPAfileAlloc (); + +// load the pmFPAfile information from the camera configuration data +pmFPAfile *pmFPAfileDefine (psMetadata *files, psMetadata *camera, pmFPA *fpa, char *name); + +// load the pmFPAfile information from the camera configuration data, constructing the needed fpa structure +// XXX deprecate this function? +// pmFPAfile *pmFPAfileConstruct (psMetadata *files, psMetadata *format, psMetadata *camera, char *name); + +// open the real file corresponding to the given pmFPAfile appropriate to the current view +bool pmFPAfileOpen (pmFPAfile *file, const pmFPAview *view); + +// read from the real file corresponding to the given pmFPAfile for the current view +bool pmFPAfileRead (pmFPAfile *file, const pmFPAview *view); + +bool pmFPAfileCreate (pmFPAfile *file, const pmFPAview *view); + +// write to the real file corresponding to the given pmFPAfile for the current view +bool pmFPAfileWrite (pmFPAfile *file, const pmFPAview *view); + +// close the real file corresponding to the given pmFPAfile appropriate to the current view +bool pmFPAfileClose (pmFPAfile *file, const pmFPAview *view); + +// set the state of the specified pmFPAfile to active (state == true) or inactive +// if name is NULL, set the state for all pmFPAfiles +bool pmFPAfileActivate (psMetadata *files, bool state, char *name); + +// examine all pmFPAfiles listed in the files and perform the needed I/O operations (open,read,write,close) +bool pmFPAfileIOChecks (psMetadata *files, const pmFPAview *view, pmFPAfilePlace place); + +// return an image corresponding to the current readout, from the specified file. if the pmFPAfile does not +// exist, construct the image using the given size and type (save it in a pmFPAfile??) +// psImage *pmFPAfileReadoutImage (psMetadata *files, const pmFPAview *view, char *name, int Nx, int Ny, int type); + +// select the readout from the named pmFPAfile; if the named file does not exist, +pmReadout *pmFPAfileThisReadout (psMetadata *files, const pmFPAview *view, const char *name); + +// create a file with the given name, assign it type "INTERNAL", and supply it with an image +// of the requested dimensions. (image only, mask and weight are ignored) +pmReadout *pmFPAfileCreateInternal (psMetadata *files, char *name, int Nx, int Ny, int type); + +// delete the INTERNAL file of the given name (if it exists) +bool pmFPAfileDropInternal (psMetadata *files, char *name); + +// read an image into the current view +bool pmFPAviewReadFitsImage (const pmFPAview *view, pmFPAfile *file); + +// write the components for the specified view +bool pmFPAviewWriteFitsImage (const pmFPAview *view, pmFPAfile *file); + +// look for the given argname on the argument list. find the give filename from the file rules +pmFPAfile *pmFPAfileFromArgs (bool *found, pmConfig *config, char *filename, char *argname); + +// look for the given argname on the argument list. find the give filename from the file rules +pmFPAfile *pmFPAfileFromConf (bool *found, pmConfig *config, char *filename, pmFPA *input); + +// create a new output pmFPAfile based on an existing FPA +pmFPAfile *pmFPAfileFromFPA (pmConfig *config, pmFPA *src, int xBin, int yBin, char *filename); + +// add the specified filename info (value) to the files of the given mode using the given reference name +bool pmFPAfileAddFileNames (psMetadata *files, char *name, char *value, int mode); + +// convert the rule to a name based on the current view +char *pmFPAfileNameFromRule (char *rule, pmFPAfile *file, const pmFPAview *view); + +bool pmFPAfileCopyView (pmFPA *out, pmFPA *in, const pmFPAview *view); + +bool pmFPAfileCopyStructureView (pmFPA *out, pmFPA *in, psMetadata *format, int xBin, int yBin, const pmFPAview *view); + +# endif Index: trunk/psModules/src/camera/pmFPAview.c =================================================================== --- trunk/psModules/src/camera/pmFPAview.c (revision 7017) +++ trunk/psModules/src/camera/pmFPAview.c (revision 7017) @@ -0,0 +1,267 @@ +/** @file pmFPAview.c + * + * @author EAM, IfA + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-05-01 01:55:43 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmHDUUtils.h" +#include "pmFPA.h" +#include "pmFPAview.h" + +static void pmFPAviewFree (pmFPAview *view) +{ + if (view == NULL) + return; + return; +} + +pmFPAview *pmFPAviewAlloc (int nRows) +{ + pmFPAview *view = psAlloc (sizeof(pmFPAview)); + psMemSetDeallocator (view, (psFreeFunc) pmFPAviewFree); + + view->chip = -1; + view->cell = -1; + view->readout = -1; + view->iRows = 0; + view->nRows = nRows; + return (view); +} + +pmFPAdepth pmFPAviewDepth (const pmFPAview *view) +{ + + if (view->chip < 0) { + return PM_FPA_DEPTH_FPA; + } + if (view->cell < 0) { + return PM_FPA_DEPTH_CHIP; + } + if (view->readout < 0) { + return PM_FPA_DEPTH_CELL; + } + return PM_FPA_DEPTH_READOUT; +} + +pmChip *pmFPAviewThisChip (const pmFPAview *view, pmFPA *fpa) +{ + + if (view->chip < 0) { + return NULL; + } + + if (view->chip >= fpa->chips->n) { + return NULL; + } + + pmChip *chip = fpa->chips->data[view->chip]; + return chip; +} + +pmChip *pmFPAviewNextChip (pmFPAview *view, pmFPA *fpa, int nStep) +{ + view->cell = -1; + view->readout = -1; + view->iRows = 0; + + // if there are no available chips, return NULL + if (fpa->chips->n <= 0) { + view->chip = -1; + return NULL; + } + + // clean up < -1 values + if (view->chip < -1) { + view->chip = -1; + } + + // increment to the next chip + view->chip += nStep; + + // if we are at the end of the stack, return NULL + if (view->chip >= fpa->chips->n) { + view->chip = -1; + return NULL; + } + + // get the correct chip pointer + pmChip *chip = fpa->chips->data[view->chip]; + return (chip); +} + +pmCell *pmFPAviewThisCell (const pmFPAview *view, pmFPA *fpa) +{ + + if (view->cell < 0) { + return NULL; + } + + pmChip *chip = pmFPAviewThisChip (view, fpa); + if (chip == NULL) { + return NULL; + } + + if (view->cell >= chip->cells->n) { + return NULL; + } + + pmCell *cell = chip->cells->data[view->cell]; + return cell; +} + +pmCell *pmFPAviewNextCell (pmFPAview *view, pmFPA *fpa, int nStep) +{ + + pmChip *chip = pmFPAviewThisChip (view, fpa); + if (chip == NULL) { + return NULL; + } + + view->readout = -1; + view->iRows = 0; + + // if there are no available cells, return NULL + if (chip->cells->n <= 0) { + view->cell = -1; + return NULL; + } + + // clean up < -1 values + if (view->cell < -1) { + view->cell = -1; + } + + // increment to the next cell + view->cell += nStep; + + // if we are at the end of the stack, return NULL + if (view->cell >= chip->cells->n) { + view->cell = -1; + return NULL; + } + + // get the correct cell pointer + pmCell *cell = chip->cells->data[view->cell]; + return (cell); +} + +pmReadout *pmFPAviewThisReadout (const pmFPAview *view, pmFPA *fpa) +{ + + if (view->readout < 0) { + return NULL; + } + + pmCell *cell = pmFPAviewThisCell (view, fpa); + if (cell == NULL) { + return NULL; + } + + if (view->readout >= cell->readouts->n) { + return NULL; + } + + pmReadout *readout = cell->readouts->data[view->readout]; + return readout; +} + +pmReadout *pmFPAviewNextReadout (pmFPAview *view, pmFPA *fpa, int nStep) +{ + + pmCell *cell = pmFPAviewThisCell (view, fpa); + if (cell == NULL) { + return NULL; + } + + view->iRows = 0; + + // if there are no available cells, return NULL + if (cell->readouts->n <= 0) { + view->readout = -1; + return NULL; + } + + // clean up < -1 values + if (view->readout < -1) { + view->readout = -1; + } + + // increment to the next cell + view->readout += nStep; + + // if we are at the end of the stack, return NULL + if (view->readout >= cell->readouts->n) { + view->readout = -1; + return NULL; + } + + // get the correct cell pointer + pmReadout *readout = cell->readouts->data[view->readout]; + return (readout); +} + +pmHDU *pmFPAviewThisHDU (const pmFPAview *view, pmFPA *fpa) +{ + // the HDU is attached to a cell, chip or fpa + // if this view has a -1 for the level which contains the hdu, + // there is no unambiguous HDU + + if (view->chip < 0) { + return pmHDUFromFPA (fpa); + } + if (view->cell < 0) { + return pmHDUFromChip (pmFPAviewThisChip (view, fpa)); + } + if (view->readout < 0) { + return pmHDUFromCell (pmFPAviewThisCell (view, fpa)); + } + return pmHDUFromReadout (pmFPAviewThisReadout (view, fpa)); +} + +pmHDU *pmFPAviewThisPHU (const pmFPAview *view, pmFPA *fpa) +{ + // select the HDU which corresponds to the PHU containing this view + + pmHDU *hdu; + pmFPAview new; + pmChip *chip; + pmCell *cell; + + new = *view; + + if (view->chip < 0) { + hdu = pmHDUFromFPA (fpa); + if (hdu->phu) + return hdu; + return NULL; + } + if (view->cell < 0) { + chip = pmFPAviewThisChip (view, fpa); + hdu = pmHDUFromChip (chip); + if (hdu->phu) + return hdu; + new.chip = -1; + hdu = pmFPAviewThisPHU (&new, fpa); + return hdu; + } + if (view->readout < 0) { + cell = pmFPAviewThisCell (view, fpa); + hdu = pmHDUFromCell (cell); + if (hdu->phu) + return hdu; + new.cell = -1; + hdu = pmFPAviewThisPHU (&new, fpa); + return hdu; + } + return NULL; +} Index: trunk/psModules/src/camera/pmFPAview.h =================================================================== --- trunk/psModules/src/camera/pmFPAview.h (revision 7017) +++ trunk/psModules/src/camera/pmFPAview.h (revision 7017) @@ -0,0 +1,70 @@ +/** @file pmFPAview.h +* +* @brief Tools to manipulate the FPA structure elements. +* +* @ingroup AstroImage +* +* @author EAM, IfA +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-05-01 01:55:43 $ +* +* Copyright 2004-2005 Institute for Astronomy, University of Hawaii +*/ + +#ifndef PM_FPA_VIEW_H +#define PM_FPA_VIEW_H + +#include "pmFPA.h" + +/// @addtogroup AstroImage +/// @{ + +// depth of interest +typedef enum { + PM_FPA_DEPTH_NONE, + PM_FPA_DEPTH_FPA, + PM_FPA_DEPTH_CHIP, + PM_FPA_DEPTH_CELL, + PM_FPA_DEPTH_READOUT, +} pmFPAdepth; + +typedef struct +{ + int chip; // Number of the chip, or -1 for all + int cell; // Number of the cell, or -1 for all + int readout; // Number of the readout, or -1 for all + int nRows; // Maximum number of rows per readout segment read, or 0 for all + int iRows; // Starting point for this read +} +pmFPAview; + +// allocate a pmFPAview structure for this fpa and camera +pmFPAview *pmFPAviewAlloc (int nRows); + +// determine the current view depth +pmFPAdepth pmFPAviewDepth (const pmFPAview *view); + +// return the currently selected chip for this view +pmChip *pmFPAviewThisChip (const pmFPAview *view, pmFPA *fpa); + +// advance view to the next chip +pmChip *pmFPAviewNextChip (pmFPAview *view, pmFPA *fpa, int nStep); + +// return the currently selected cell for this view +pmCell *pmFPAviewThisCell (const pmFPAview *view, pmFPA *fpa); + +// advance view to the next cell +pmCell *pmFPAviewNextCell (pmFPAview *view, pmFPA *fpa, int nStep); + +// return the currently selected readout for this view +pmReadout *pmFPAviewThisReadout (const pmFPAview *view, pmFPA *fpa); + +// advance view to the next readout +pmReadout *pmFPAviewNextReadout (pmFPAview *view, pmFPA *fpa, int nStep); + +// return the HDU corresponding to the current view +pmHDU *pmFPAviewThisHDU (const pmFPAview *view, pmFPA *fpa); +pmHDU *pmFPAviewThisPHU (const pmFPAview *view, pmFPA *fpa); + +# endif Index: trunk/psModules/src/camera/pmHDU.c =================================================================== --- trunk/psModules/src/camera/pmHDU.c (revision 7017) +++ trunk/psModules/src/camera/pmHDU.c (revision 7017) @@ -0,0 +1,225 @@ +#include +#include + +#include "pslib.h" +#include "pmFPA.h" +#include "psAdditionals.h" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// File-static (private) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Move to the appropriate extension in FITS file for HDU +static bool hduMove(pmHDU *hdu, // HDU with extname + psFits *fits // FITS file in which to move + ) +{ + // Deal with the PHU case + if (strcasecmp(hdu->extname, "PHU") == 0 || hdu->phu) { + if (! psFitsMoveExtNum(fits, 0, false)) { + psError(PS_ERR_IO, false, "Unable to move to primary header!\n"); + return false; + } + hdu->phu = true; + return true; + } + + if (! psFitsMoveExtName(fits, hdu->extname)) { + psError(PS_ERR_IO, false, "Unable to move to extension %s\n", hdu->extname); + return false; + } + // Now, just in case for some reason the PHU has an extension name that we've moved to.... + if (psFitsGetExtNum(fits) == 0) { + hdu->phu = true; + } else { + hdu->phu = false; + } + + return true; +} + +static void hduFree(pmHDU *hdu) +{ + psFree(hdu->extname); + psFree(hdu->format); + psFree(hdu->header); + psFree(hdu->images); + psFree(hdu->table); + psFree(hdu->weights); + psFree(hdu->masks); +} + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +pmHDU *pmHDUAlloc(const char *extname // Extension name + ) +{ + pmHDU *hdu = psAlloc(sizeof(pmHDU)); + psMemSetDeallocator(hdu, (psFreeFunc)hduFree); + + if (! extname || strlen(extname) == 0) { + hdu->phu = true; + hdu->extname = psStringCopy("PHU"); + } else { + if (strcasecmp(extname, "PHU") == 0) { + hdu->phu = true; + } else { + hdu->phu = false; + } + hdu->extname = psStringCopy(extname); + } + hdu->format = NULL; + hdu->header = NULL; + hdu->images = NULL; + hdu->table = NULL; + hdu->weights = NULL; + hdu->masks = NULL; + + return hdu; +} + +bool pmHDUReadHeader(pmHDU *hdu, // HDU for which to read header + psFits *fits // FITS file from which to read + ) +{ + assert(hdu); + assert(fits); + + // Move to the appropriate extension + psTrace(__func__, 5, "Moving to extension %s...\n", hdu->extname); + if (!hduMove(hdu, fits)) { + return false; + } + + if (!hdu->header) { + psTrace(__func__, 5, "Reading the header...\n"); + hdu->header = psFitsReadHeader(NULL, fits); + if (! hdu->header) { + psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", hdu->extname); + return false; + } + } + + return true; +} + +// Read the HDU +// XXX: Add a region specifier? +bool pmHDURead(pmHDU *hdu, // HDU to read + psFits *fits // FITS file to read from + ) +{ + assert(hdu); + assert(fits); + + // Read the header; includes the move + if (!pmHDUReadHeader(hdu, fits)) { + return false; + } + + #ifdef FITS_TABLES + // What type is it? + if (psFitsIsImage(hdu->header)) { + #endif + if (hdu->images) { + psLogMsg(__func__, PS_LOG_WARN, "HDU %s has already been read --- overwriting.\n", hdu->extname); + psFree(hdu->images); // Blow away anything existing + } + psTrace(__func__, 5, "Reading the pixels...\n"); + hdu->images = psFitsReadImageCube(fits, psRegionSet(0,0,0,0)); + if (! hdu->images) { + psError(PS_ERR_IO, false, "Unable to read pixels for extension %s\n", hdu->extname); + return false; + } + return true; + #ifdef FITS_TABLES + + } + if (psFitsIsTable(hdu->table)) { + // Read the table + if (hdu->table) { + psFree(hdu->table); // Blow away anything existing + } + hdu->table = psFitsReadTable(fits); + return true; + } + + psError(PS_ERR_UNKNOWN, true, "No idea what this HDU consists of!\n"); + return false; + #endif +} + +// Write the HDU +// XXX: Add a region specifier? +bool pmHDUWrite(pmHDU *hdu, // HDU to write + psFits *fits // FITS file to write to + ) +{ + assert(hdu); + assert(fits); + + psTrace(__func__, 7, "Writing HDU %s\n", hdu->extname); + + if (hdu->images && hdu->table && !hdu->header) { + psError(PS_ERR_IO, true, "Both image and table data provided in HDU, but no header --- " + "don't know which to write!\n"); + return false; + } + + if (!hdu->images && !hdu->table && !hdu->header) { + psLogMsg(__func__, PS_LOG_WARN, "Nothing to write for HDU %s\n", hdu->extname); + return false; + } + + // Preserve the extension name, if it's the PHU + char *extname = hdu->extname; // The name of the extension + if (strcasecmp(extname, "PHU") == 0 && hdu->header) { + bool mdok = true; // Status of MD lookup + extname = psMetadataLookupStr(&mdok, hdu->header, "EXTNAME"); + if (!mdok || !extname || strlen(extname) == 0) { + extname = ""; + } + } + + // Only a header + if (!hdu->images && !hdu->table) { + // Tell CFITSIO there's nothing there + psMetadataItem *naxis = psMetadataLookup(hdu->header, "NAXIS"); + if (naxis) { + naxis->data.S32 = 0; + } + + if (!psFitsWriteHeader(hdu->header, fits)) { + psError(PS_ERR_IO, false, "Unable to write header for extension %s\n", extname); + } + } + + #ifdef FITS_TABLES + if (hdu->images && (!hdu->table || psFitsIsImage(hdu->header))) { + psFitsWriteImageCube(fits, hdu->header, hdu->images, extname); + return true; + } + + if (hdu->table && (!hdu->images || psFitsIsTable(hdu->header))) { + bool psFitsWriteTable(psFits* fits, const psMetadata *header, const psArray* table); + return true; + } + + psError(PS_ERR_IO, true, "No idea what this HDU consists of!\n"); + return false; + #else + + if (hdu->images) { + psTrace(__func__, 9, "Writing pixels for %s\n", hdu->extname); + if (!psFitsWriteImageCube(fits, hdu->header, hdu->images, extname)) { + psError(PS_ERR_IO, false, "Unable to write image to extension %s\n", hdu->extname); + return false; + } + } + return true; + #endif +} + Index: trunk/psModules/src/camera/pmHDU.h =================================================================== --- trunk/psModules/src/camera/pmHDU.h (revision 7017) +++ trunk/psModules/src/camera/pmHDU.h (revision 7017) @@ -0,0 +1,37 @@ +#ifndef PM_HDU_H +#define PM_HDU_H + +// An instance of the FITS Header Data Unit +typedef struct +{ + psString extname; // The extension name + bool phu; // Is this the FITS Primary Header Unit + psMetadata *format; // The camera format + psMetadata *header; // The FITS header, or NULL if primary for FITS; or section info + psArray *images; // The pixel data + psArray *weights; // The pixel data + psArray *masks; // The pixel data + psArray *table; // The table data +} +pmHDU; + + +// Allocators +pmHDU *pmHDUAlloc(const char *extname); + +// Read the header +bool pmHDUReadHeader(pmHDU *hdu, // HDU for which to read header + psFits *fits // FITS file from which to read + ); + +// Read the HDU +bool pmHDURead(pmHDU *hdu, // HDU to read + psFits *fits // FITS file to read from + ); + +// Write the HDU +bool pmHDUWrite(pmHDU *hdu, // HDU to write + psFits *fits // FITS file to write to + ); + +#endif Index: trunk/psModules/src/camera/pmHDUUtils.c =================================================================== --- trunk/psModules/src/camera/pmHDUUtils.c (revision 7017) +++ trunk/psModules/src/camera/pmHDUUtils.c (revision 7017) @@ -0,0 +1,60 @@ +#include +#include "pmFPA.h" +#include "pmHDU.h" +#include "pmHDUUtils.h" + + +pmHDU *pmHDUFromFPA(pmFPA *fpa // FPA for which to find HDU + ) +{ + return fpa->hdu; +} + +pmHDU *pmHDUFromChip(pmChip *chip // Chip for which to find HDU + ) +{ + pmHDU *hdu = chip->hdu; // The HDU information + if (!hdu) { + hdu = pmHDUFromFPA(chip->parent); // Grab HDU info from the FPA + } + + return hdu; +} + +pmHDU *pmHDUFromCell(pmCell *cell // Cell for which to find HDU + ) +{ + pmHDU *hdu = cell->hdu; // The HDU information + if (!hdu) { + hdu = pmHDUFromChip(cell->parent); // Grab HDU info from the chip + } + + return hdu; +} + +pmHDU *pmHDUFromReadout (pmReadout *readout) +{ + + pmCell *cell = readout->parent; // cell containing this readout; + pmHDU *hdu = pmHDUFromCell (cell); + return hdu; +} + +// Get the lowest HDU +pmHDU *pmHDUGetLowest(pmFPA *fpa, // The FPA + pmChip *chip, // The chip, or NULL + pmCell *cell // The cell, or NULL + ) +{ + pmHDU *hdu = NULL; // The HDU that's at the lowest level + if (cell) { + hdu = pmHDUFromCell(cell); + } else if (chip) { + hdu = pmHDUFromChip(chip); + } else if (fpa) { + hdu = pmHDUFromFPA(fpa); + } + + return hdu; +} + Index: trunk/psModules/src/camera/pmHDUUtils.h =================================================================== --- trunk/psModules/src/camera/pmHDUUtils.h (revision 7017) +++ trunk/psModules/src/camera/pmHDUUtils.h (revision 7017) @@ -0,0 +1,24 @@ +#ifndef PM_HDU_UTILS_H +#define PM_HDU_UTILS_H + +#include "pmFPA.h" +#include "pmHDU.h" + +// Get the lowest HDU in the hierarchy, as supplied +pmHDU *pmHDUGetLowest(pmFPA *fpa, // The FPA + pmChip *chip, // The chip, or NULL + pmCell *cell // The cell, or NULL + ); + +// Find the HDU in the FPA hierarchy +pmHDU *pmHDUFromFPA(pmFPA *fpa // FPA for which to find HDU + ); +pmHDU *pmHDUFromChip(pmChip *chip // Chip for which to find HDU + ); +pmHDU *pmHDUFromCell(pmCell *cell // Cell for which to find HDU + ); + +pmHDU *pmHDUFromReadout (pmReadout *readout // Readout for which to find HDU + ); + +#endif Index: trunk/psModules/src/camera/pmReadout.c =================================================================== --- trunk/psModules/src/camera/pmReadout.c (revision 7017) +++ trunk/psModules/src/camera/pmReadout.c (revision 7017) @@ -0,0 +1,111 @@ +#include +#include "pslib.h" +#include "pmFPA.h" +#include "pmMaskBadPixels.h" +#include "pmHDUUtils.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))) { + psImage *bias = psMemIncrRefCounter(psImageSubset(image, *region)); // Image from bias section + psListAdd(images, PS_LIST_TAIL, bias); + psFree(bias); + } + psFree(sectionsIter); + + return images; +} + +bool pmReadoutSetWeights(pmReadout *readout) +{ + + pmCell *cell = readout->parent; + + float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // Cell gain + float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE"); // Cell read noise + psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); // Trim section + float saturation = psMetadataLookupF32(NULL, cell->concepts, "CELL.SATURATION"); // Saturation level + float bad = psMetadataLookupF32(NULL, cell->concepts, "CELL.BAD"); // Bad level + + pmHDU *hdu = pmHDUFromCell (cell); + + psArray *pixels = hdu->images; // Array of images + psArray *weights = hdu->weights; // Array of weight images + psArray *masks = hdu->masks; // Array of mask images + // Generate the weights and masks if required + if (! weights) { + weights = psArrayAlloc(pixels->n); + weights->n = pixels->n; + for (int i = 0; i < pixels->n; i++) { + psImage *image = pixels->data[i]; + weights->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); + psImageInit(weights->data[i], 0.0); + } + hdu->weights = weights; + } + if (! masks) { + masks = psArrayAlloc(pixels->n); + masks->n = pixels->n; + for (int i = 0; i < pixels->n; i++) { + psImage *image = pixels->data[i]; + masks->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8); + psImageInit(masks->data[i], 0); + } + hdu->masks = masks; + } + + // Set the pixels + psArray *readouts = cell->readouts; // Array of readouts + for (int i = 0; i < readouts->n; i++) { + pmReadout *readout = readouts->data[i]; // The readout of interest + + if (! readout->weight) { + readout->weight = psMemIncrRefCounter(psImageSubset(weights->data[i], *trimsec)); + } + if (! readout->mask) { + readout->mask = psMemIncrRefCounter(psImageSubset(masks->data[i], *trimsec)); + } + + // Set up the mask + psImage *image = readout->image;// Pixels + psImage *mask = readout->mask; // Mask image + for (int i = 0; i < image->numRows; i++) { + for (int j = 0; j < image->numCols; j++) { + if (image->data.F32[i][j] >= saturation) { + mask->data.U8[i][j] = PM_MASK_SAT; + } + if (image->data.F32[i][j] <= bad) { + mask->data.U8[i][j] = PM_MASK_BAD; + } + } + } + + psImage *weight = readout->weight; // Mask image + float rnoise = PS_SQR(readnoise/gain); + for (int i = 0; i < image->numRows; i++) { + for (int j = 0; j < image->numCols; j++) { + if (!mask->data.U8[i][j]) { + weight->data.F32[i][j] = image->data.F32[i][j]/ gain + rnoise; + } + } + } + + // Set weight image to the variance = g*f + rn^2 + // psBinaryOp(readout->weight, image, "/", psScalarAlloc(gain, PS_TYPE_F32)); + // psBinaryOp(readout->weight, readout->weight, "+", + // psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32)); + } + + return true; +} + Index: trunk/psModules/src/camera/pmReadout.h =================================================================== --- trunk/psModules/src/camera/pmReadout.h (revision 7017) +++ trunk/psModules/src/camera/pmReadout.h (revision 7017) @@ -0,0 +1,13 @@ +#ifndef PM_READOUT_H +#define PM_READOUT_H + +#include "pslib.h" +#include "pmFPA.h" + +// Get the bias sections for a specific readout +psList *pmReadoutGetBias(pmReadout *readout // Readout for which to get the bias sections + ); + +bool pmReadoutSetWeights(pmReadout *readout); + +#endif