Index: trunk/psModules/src/astrom/pmChipMosaic.c =================================================================== --- trunk/psModules/src/astrom/pmChipMosaic.c (revision 6872) +++ trunk/psModules/src/astrom/pmChipMosaic.c (revision 6985) @@ -1,10 +1,20 @@ #include #include - #include "pslib.h" #include "pmFPA.h" -#include "pmChipMosaic.h" - -#define MEM_LEAKS 0 +#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 @@ -17,12 +27,204 @@ } +// 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 -psImage *p_pmImageMosaic(const psArray *source, // Images to splice in - const psVector *xFlip, const psVector *yFlip, // Need to flip x and y? - const psVector *xBinSource, const psVector *yBinSource, // Binning in x and y of - // source images - int xBinTarget, int yBinTarget, // Binning in x and y of target images - const psVector *x0, const psVector *y0 // Offsets for source images on target - ) +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); @@ -33,4 +235,14 @@ 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 @@ -40,4 +252,5 @@ 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 @@ -45,4 +258,5 @@ continue; } + numImages++; // Only implemented for F32 and U8 images so far. @@ -66,4 +280,7 @@ 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; } @@ -137,5 +354,6 @@ static bool cellConcepts(pmCell *target,// Target cell psArray *sources, // Source cells - int xBin, int yBin // Binning + int xBin, int yBin, // Binning + psRegion *trimsec // The trim section ) { @@ -143,6 +361,6 @@ float gain = 0.0; // Gain float readnoise = 0.0; // Read noise - float saturation = 0.0; // Saturation level - float bad = 0.0; // Bad level + float saturation = INFINITY; // Saturation level + float bad = -INFINITY; // Bad level float exposure = 0.0; // Exposure time float darktime = 0.0; // Dark time @@ -160,9 +378,8 @@ gain += psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); readnoise += psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE"); - saturation += psMetadataLookupF32(NULL, cell->concepts, "CELL.SATURATION"); - bad += psMetadataLookupF32(NULL, cell->concepts, "CELL.BAD"); exposure += psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); darktime += psMetadataLookupF32(NULL, cell->concepts, "CELL.DARKTIME"); - time += psTimeToMJD(psMetadataLookupPtr(NULL, cell->concepts, "CELL.TIME")); + psTime *cellTime = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TIME"); + time += psTimeToMJD(cellTime); if (i == 0) { timeSys = psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS"); @@ -172,35 +389,48 @@ success = false; } - } - gain /= (float)nCells; - readnoise /= (float)nCells; - saturation /= (float)nCells; - bad /= (float)nCells; - exposure /= (float)nCells; - darktime /= (float)nCells; - psTime *timePtr = psTimeFromMJD(time/(double)nCells); - timePtr = psTimeConvert(timePtr, timeSys); - - // XXX *REALLY* need a generic "concept update" function that handles the type and comments transparently. - psMetadataAddF32(target->concepts, PS_LIST_TAIL, "CELL.GAIN", PS_META_REPLACE, "Gain (e/ADU)", gain); - psMetadataAddF32(target->concepts, PS_LIST_TAIL, "CELL.READNOISE", PS_META_REPLACE, "Read noise (e)", readnoise); - psMetadataAddF32(target->concepts, PS_LIST_TAIL, "CELL.SATURATION", PS_META_REPLACE, "Saturation level (ADU)", saturation); - psMetadataAddF32(target->concepts, PS_LIST_TAIL, "CELL.BAD", PS_META_REPLACE, "Bad level (ADU)", bad); - psMetadataAddF32(target->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)", exposure); - psMetadataAddF32(target->concepts, PS_LIST_TAIL, "CELL.DARKTIME", PS_META_REPLACE, "Time since last CCD flush (sec)", darktime); - psMetadataAddPtr(target->concepts, PS_LIST_TAIL, "CELL.TIME", PS_DATA_TIME | PS_META_REPLACE, "Time of observation", timePtr); - psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.TIMESYS", PS_META_REPLACE, "Time system", timeSys); - psFree(timePtr); - - // Now fill in the ones I know by other means - psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.X0", PS_META_REPLACE, "Position of (0,0) on the chip", 0); - psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.Y0", PS_META_REPLACE, "Position of (0,0) on the chip", 0); - psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.XPARITY", PS_META_REPLACE, "Orientation in x compared to the rest of the FPA", 1); - psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.YPARITY", PS_META_REPLACE, "Orientation in x compared to the rest of the FPA", 1); - psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_META_REPLACE, "Binning in x", xBin); - psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_META_REPLACE, "Binning in x", yBin); - psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.READDIR", PS_META_REPLACE, "Read direction (faked)", 1); - psRegion *trimsec = psMetadataLookupPtr(NULL, target->concepts, "CELL.TRIMSEC"); - trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0.0; + + 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; @@ -208,52 +438,84 @@ - -// Mosaic a chip together into a single image -int pmChipMosaic(pmChip *chip,// Chip to mosaic - int xBinChip, int yBinChip // Binning of mosaic image in x and y - ) +// 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 - psArray *images = psArrayAlloc(cells->n); // Array of images that will be mosaicked - psArray *weights = psArrayAlloc(cells->n); // Array of weight images to be mosaicked - psArray *masks = psArrayAlloc(cells->n); // Array of mask images to be mosaicked - psVector *x0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin x coordinates - psVector *y0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin y coordinates - psVector *xBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in x - psVector *yBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in y - psVector *xFlip = psVectorAlloc(cells->n, PS_TYPE_U8); // Flip in x? - psVector *yFlip = psVectorAlloc(cells->n, PS_TYPE_U8); // Flip in y? + 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", cells->n); - for (int i = 0; i < cells->n; i++) { + 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; } - x0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); - y0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + 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(NULL, cell->concepts, "CELL.XBIN"); - yBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); - int xParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); - int yParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); - if (xParity == 1) { - xFlip->data.U8[i] = 0; - } else if (xParity == -1) { + 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 { - psLogMsg(__func__, PS_LOG_WARN, "The x parity of cell %d is not +/- 1 (it's %d) --- " - "assuming +1.\n", i, xParity); xFlip->data.U8[i] = 0; } - if (yParity == 1) { - yFlip->data.U8[i] = 0; - } else if (yParity == -1) { + if (yParity == -1) { yFlip->data.U8[i] = 1; } else { - psLogMsg(__func__, PS_LOG_WARN, "The y parity of cell %d is not +/- 1 (it's %d) --- " - "assuming +1.\n", i, yParity); yFlip->data.U8[i] = 0; } @@ -271,97 +533,102 @@ masks->data[i] = psMemIncrRefCounter(readout->mask); } + // Mosaic the images together and we're done - psImage *image = p_pmImageMosaic(images, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0); - psImage *weight = p_pmImageMosaic(weights, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0); - psImage *mask = p_pmImageMosaic(masks, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0); + 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(x0); - psFree(y0); - psFree(xBin); - psFree(yBin); - psFree(xFlip); - psFree(yFlip); psFree(images); psFree(weights); psFree(masks); - int nCells = cells->n; - - // Fix up the HDU - if (chip->parent->hdu) { - psLogMsg(__func__, PS_LOG_WARN, "The original format has the entire FPA in a single extension. " - "The FPA hierarchy may be invalid following the pmChipMosaic.\n"); + 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 { - if (! chip->hdu) { - psString chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); - chip->hdu = p_pmHDUAlloc(chipName); - } - p_pmHDU *hdu = chip->hdu; - psArrayElementsFree(hdu->images); - psArrayElementsFree(hdu->weights); - psArrayElementsFree(hdu->masks); - hdu->images = psArrayRealloc(hdu->images,1); - hdu->weights = psArrayRealloc(hdu->weights, 1); - hdu->masks = psArrayRealloc(hdu->masks, 1); - hdu->images->data[0] = image; - hdu->weights->data[0] = weight; - hdu->masks->data[0] = mask; - psMetadataAddS32(hdu->header, PS_LIST_TAIL, "NAXIS1", PS_META_REPLACE, "Number of columns", image->numCols); - psMetadataAddS32(hdu->header, PS_LIST_TAIL, "NAXIS2", PS_META_REPLACE, "Number of rows", image->numRows); - } - - // Chop off all the component cells, and construct a new one - pmCell *newCell = pmCellAlloc(NULL, NULL, __func__); // New cell - cellConcepts(newCell, cells, xBinChip, yBinChip); - pmChipFreeCells(chip); - // Have to put in the new cell manually, since we didn't want to put it in before blowing the cells away. - newCell->parent = chip; - psArrayAdd(chip->cells, 1, newCell); - newCell->exists = true; - newCell->process = true; + // 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 - // Want the readouts to contain a subimage, but that subimage is the whole image. - // This preserves the relationship there was before, where freeing the parent frees the child. - psRegion entire = {0.0, 0.0, 0.0, 0.0}; - newReadout->image = psMemIncrRefCounter(psImageSubset(image, entire)); - newReadout->weight = psMemIncrRefCounter(psImageSubset(weight, entire)); - newReadout->mask = psMemIncrRefCounter(psImageSubset(mask, entire)); - // Drop references - psFree(newReadout); - psFree(newCell); - - // Well, we've stuffed around with the camera configuration, so it's no longer valid... - #if 0 - - psFree(chip->parent->camera); - chip->parent->camera = NULL; - #endif - - return nCells; -} - - -int pmFPAMosaicCells(pmFPA *fpa, // FPA - int xBinChip, int yBinChip // Binning of mosaic image in x and y - ) -{ - assert(fpa); - - int numChips = 0; - psArray *chips = fpa->chips; // Component chips - for (int i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // The chip of interest - if (! chip || ! chip->exists || ! chip->process) { - continue; - } - - if (pmChipMosaic(chip, xBinChip, yBinChip) > 0) { - numChips++; - } - } - - return numChips; - -} + newReadout->image = mosaicImage; + newReadout->mask = mosaicMask; + newReadout->weight = mosaicWeights; + psFree(newReadout); // Drop reference + + return true; +}