Index: /trunk/psModules/src/camera/Makefile.am =================================================================== --- /trunk/psModules/src/camera/Makefile.am (revision 7486) +++ /trunk/psModules/src/camera/Makefile.am (revision 7487) @@ -9,4 +9,5 @@ pmFPAHeader.c \ pmFPAMaskWeight.c \ + pmFPAMosaic.c \ pmFPARead.c \ pmFPAUtils.c \ @@ -15,5 +16,4 @@ pmHDUUtils.c \ pmHDUGenerate.c \ - pmChipMosaic.c \ pmFPA_JPEG.c \ pmFPAview.c \ @@ -27,4 +27,5 @@ pmFPAHeader.h \ pmFPAMaskWeight.h \ + pmFPAMosaic.h \ pmFPARead.h \ pmFPAUtils.h \ @@ -33,5 +34,4 @@ pmHDUUtils.h \ pmHDUGenerate.h \ - pmChipMosaic.h \ pmFPA_JPEG.h \ pmFPAview.h \ Index: unk/psModules/src/camera/pmChipMosaic.c =================================================================== --- /trunk/psModules/src/camera/pmChipMosaic.c (revision 7486) +++ (revision ) @@ -1,1101 +1,0 @@ -#include -#include -#include -#include "pmFPA.h" -#include "pmHDU.h" -#include "pmHDUUtils.h" -#include "pmChipMosaic.h" - -#define CELL_LIST_BUFFER 50 // Buffer size for cell lists - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// 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; \ -} - - -// Get the bounds for an chip's pixels on the HDU -static bool chipBounds(psRegion *bounds, // The bounds for the chip - const pmChip *chip // The chip to examine for contiguity - ) -{ - assert(chip); - - 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 || psRegionIsNaN(*trimsec)) { - psError(PS_ERR_UNKNOWN, true, "CELL.TRIMSEC hasn't been set for cell %d.\n", i); - return false; - } - - 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; - } - } - - return true; -} - -// Make sure the TRIMSEC doesn't overlap with the established image bounds -static bool chipContiguousTrimsec(psRegion *bounds, // The bounds of the image, altered if primary==true - const pmChip *chip // The chip to examine for contiguity - ) -{ - assert(bounds); - assert(chip); - - 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 || psRegionIsNaN(*trimsec)) { - psError(PS_ERR_UNKNOWN, true, "CELL.TRIMSEC hasn't been set for cell %d.\n", i); - return false; - } - - if (REGIONS_OVERLAP(trimsec, bounds)) { - return false; - } - } - - return true; -} - -// Make sure the BIASSEC doesn't overlap with the established image bounds -static bool chipContiguousBiassec(psRegion *bounds, // The bounds of the image, altered if primary==true - const pmChip *chip // The chip to examine for contiguity - ) -{ - assert(bounds); - assert(chip); - - // Check that the biases don't get in the way - 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 - psList *biassecs = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.BIASSEC"); // Bias sections - if (!mdok || !biassecs) { - psError(PS_ERR_UNKNOWN, 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 (psRegionIsNaN(*biassec)) { - continue; - } - if (REGIONS_OVERLAP(biassec, bounds)) { - psFree(biassecsIter); - return false; - } - } - psFree(biassecsIter); - } - - // If we've gotten this far, everything is fine. - return true; -} - -// Are the pixels for the FPA contiguous on the HDU? -// Work this out by examining all the CELL.TRIMSEC and CELL.BIASSEC regions for the component cells -static bool fpaContiguous(psRegion *bounds, // The bounds of the image, returned - const pmFPA *fpa // The FPA to examine for contiguity - ) -{ - assert(bounds); - assert(fpa); - - *bounds = psRegionSet(INFINITY, 0, INFINITY, 0); - - // Get the size of the pixels on the HDU - psArray *chips = fpa->chips; // The array of chips - for (int i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // Chip of interest - if (!chipBounds(bounds, chip)) { - return false; - } - } - - // Make sure the bias regions don't get in the way of the HDU - for (int i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // Chip of interest - if (!chipContiguousBiassec(bounds, chip)) { - return false; - } - } - - // If we got through it all, they must all be contiguous - return true; -} - - - -// Check a cell for niceness in the parity and binning -static bool niceCellParityBinning(int *xBin, int *yBin, // Binning for cell, to be returned - const pmCell *cell // Cell to check for niceness - ) -{ - assert(xBin); - assert(yBin); - assert(cell); - - // A "nice" cell must have only a single readout - if (cell->readouts->n != 1) { - return false; - } - - // A "nice" cell must have parity == 1 - bool mdok = true; // Status of MD lookup - int xParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XPARITY"); // Parity in x - if (!mdok || xParity == 0) { - psError(PS_ERR_UNKNOWN, true, "CELL.XPARITY hasn't been set for cell.\n"); - return false; - } - if (xParity != 1) { - return false; - } - int yParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YPARITY"); // Parity in y - if (!mdok || yParity == 0) { - psError(PS_ERR_UNKNOWN, true, "CELL.YPARITY hasn't been set for cell.\n"); - return false; - } - if (yParity != 1) { - return false; - } - - // A "nice" cell must have consistent binning - int xBinCell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XBIN"); // Binning in x - if (!mdok || xBin <= 0) { - psError(PS_ERR_UNKNOWN, true, "CELL.XPARITY hasn't been set for cell.\n"); - return false; - } - int yBinCell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YBIN"); // Binning in y - if (!mdok || yBin <= 0) { - psError(PS_ERR_UNKNOWN, true, "CELL.YPARITY hasn't been set for cell.\n"); - return false; - } - if (*xBin == 0 || *yBin == 0) { - *xBin = xBinCell; - *yBin = yBinCell; - } else if (xBinCell != *xBin || yBinCell != *yBin) { - return false; - } - - return true; -} - - -// Check a cell for niceness in the boundaries -static bool niceCellBounds(const pmCell *cell, // Cell to check for niceness - const psRegion *imageBounds // Bounds of the image on the HDU - ) -{ - // A "nice" cell must have the (0,0) pixel at CELL.X0,CELL.Y0 - bool mdok = true; // Status of MD lookup - int x0 = psMetadataLookupS32(&mdok, cell->concepts, "CELL.X0"); // Position of (0,0) on chip - if (!mdok) { - psError(PS_ERR_UNKNOWN, true, "CELL.X0 hasn't been set for cell.\n"); - return false; - } - int y0 = psMetadataLookupS32(&mdok, cell->concepts, "CELL.Y0"); // Position of (0,0) on chip - if (!mdok) { - psError(PS_ERR_UNKNOWN, true, "CELL.Y0 hasn't been set for cell.\n"); - return false; - } - pmReadout *readout = cell->readouts->data[0]; // A representative readout - if (!readout) { - return false; // 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); - return false; - } - - 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 - const pmChip *chip // Chip to examine for "niceness". - ) -{ - assert(xBinChip); - assert(yBinChip); - assert(chip); - - // 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 - *xBinChip = 0; - *yBinChip = 0; - for (int i = 0; i < chip->cells->n; i++) { - pmCell *cell = chip->cells->data[i]; // The cell of interest - if (!niceCellParityBinning(xBinChip, yBinChip, cell)) { - return NULL; - } - } - - // Now check that the pixels are all contiguous - psRegion *imageBounds = psRegionAlloc(INFINITY, 0, INFINITY, 0); // Bound of image on HDU - if (!chipBounds(imageBounds, chip) || !chipContiguousBiassec(imageBounds, chip)) { - 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 - if (!niceCellBounds(cell, imageBounds)) { - 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 (!chipContiguousTrimsec(imageBounds, testChip) || - !chipContiguousBiassec(imageBounds, testChip)) { - psTrace(__func__, 5, "Image isn't contiguous.\n"); - psFree(imageBounds); - return NULL; - } - } - } - - return imageBounds; -} - -// Is the FPA "nice"? If so, return the region containing the FPA pixels -static psRegion *niceFPA(int *xBinFPA, int *yBinFPA, // Binning for FPA, to be returned - const pmFPA *fpa // FPA to examine for "niceness". - ) -{ - assert(xBinFPA); - assert(yBinFPA); - assert(fpa); - - // Check that we've got the HDU in the chip or the FPA - if (!fpa->hdu || !fpa->hdu->images) { - return NULL; - } - - // Check parity and binning for component cells - *xBinFPA = 0; - *yBinFPA = 0; - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; // The chip of interest - for (int j = 0; i < chip->cells->n; i++) { - pmCell *cell = chip->cells->data[j]; // The cell of interest - if (!niceCellParityBinning(xBinFPA, yBinFPA, cell)) { - return NULL; - } - } - } - - // Now check that the pixels are all contiguous - psRegion *imageBounds = psRegionAlloc(0, 0, 0, 0); // Bound of image on HDU - if (!fpaContiguous(imageBounds, fpa)) { - 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 < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; // The chip of interest - for (int j = 0; i < chip->cells->n; i++) { - pmCell *cell = chip->cells->data[j]; // The cell of interest - if (!niceCellBounds(cell, imageBounds)) { - 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 - psTrace(__func__, 3, "Mosaicking %ld cells.\n", source->n); - 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] - xMin, y0->data.S32[i] - yMin, "+"); - } 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 - yMin) / - 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 - xMin) / - 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 - psList *sources, // Source cells - int xBin, int yBin, // Binning - psRegion *trimsec // The trim section - ) -{ - assert(target); - assert(sources); - assert(xBin > 0 && yBin > 0); - assert(trimsec); - - 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 readdir = 0; // Cell read direction - - int nCells = 0; // Number of cells; - psListIterator *sourcesIter = psListIteratorAlloc(sources, PS_LIST_HEAD, false); // Iterator for sources - pmCell *cell = NULL; // Source cell from iteration - while ((cell = psListGetAndIncrement(sourcesIter))) { - 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 (nCells == 1) { - timeSys = psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS"); - readdir = psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR"); - } 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; - } - if (readdir != psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR")) { - psLogMsg(__func__, PS_LOG_ERROR, "Differing cell read directions in use: %d vs %d\n", readdir, - psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR")); - 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; - } - } - psFree(sourcesIter); - - 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); - MD_UPDATE(target->concepts, "CELL.READDIR", S32, readdir); - - // 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; -} - - -// Add a cell and its various properties to the arrays -static bool addCell(psArray *images, // Array of images - psArray *masks, // Array of masks - psArray *weights, // Array of weights - psVector *x0, // Array of X0 - psVector *y0, // Array of Y0 - psVector *xBin, // Array of XBIN - psVector *yBin, // Array of YBIN - psVector *xFlip, // Array indicating whether x axis should be flipped - psVector *yFlip, // Array indicating whether y axis should be flipped - const pmCell *cell, // Cell to add - int *xBinMin, // The minimum x binning, returned - int *yBinMin, // The minimum y binning, returned - bool chipStuff // Worry about chip stuff as well? - ) -{ - if (!cell) { - return false; - } - - // Expand the arrays and vectors to handle new data - long index = images->n; // The index to use - if (images->n == images->nalloc) { - images = psArrayRealloc(images, index + CELL_LIST_BUFFER); - masks = psArrayRealloc(masks, index + CELL_LIST_BUFFER); - weights = psArrayRealloc(weights, index + CELL_LIST_BUFFER); - x0 = psVectorRealloc(x0, index+ CELL_LIST_BUFFER); - y0 = psVectorRealloc(y0, index+ CELL_LIST_BUFFER); - xBin = psVectorRealloc(xBin, index+ CELL_LIST_BUFFER); - yBin = psVectorRealloc(yBin, index+ CELL_LIST_BUFFER); - xFlip = psVectorRealloc(xFlip, index+ CELL_LIST_BUFFER); - yFlip = psVectorRealloc(yFlip, index+ CELL_LIST_BUFFER); - } - - images->n = index + 1; - masks->n = index + 1; - weights->n = index + 1; - x0->n = index + 1; - y0->n = index + 1; - xBin->n = index + 1; - yBin->n = index + 1; - xFlip->n = index + 1; - yFlip->n = index + 1; - - bool mdok = true; // Status of MD lookup - bool good = true; // Is everything good? - - // Offset of the cell on the chip - int x0Cell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.X0"); - if (!mdok) { - psError(PS_ERR_UNKNOWN, true, "CELL.X0 for cell is not set.\n"); - good = false; - } - int y0Cell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.Y0"); - if (!mdok) { - psError(PS_ERR_UNKNOWN, true, "CELL.Y0 for cell is not set.\n"); - good = false; - } - psTrace(__func__, 5, "Cell %d: x0=%d y0=%d\n", index, x0Cell, y0Cell); - - // Offset of the chip on the FPA - int x0Chip = 0, y0Chip = 0; - if (chipStuff) { - pmChip *chip = cell->parent; // The parent chip - if (!chip) { - psError(PS_ERR_UNKNOWN, true, "Cell has no parent chip --- can't find CHIP.X0 and CHIP.Y0\n"); - good = false; - } - x0Chip = psMetadataLookupS32(&mdok, chip->concepts, "CHIP.X0"); - if (!mdok) { - psError(PS_ERR_UNKNOWN, true, "CHIP.X0 for chip is not set.\n"); - good = false; - } - y0Chip = psMetadataLookupS32(&mdok, chip->concepts, "CHIP.Y0"); - if (!mdok) { - psError(PS_ERR_UNKNOWN, true, "CHIP.Y0 for chip is not set.\n"); - good = false; - } - } - - // Binning - xBin->data.S32[index] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XBIN"); - if (!mdok || xBin->data.S32[index] == 0) { - psError(PS_ERR_UNKNOWN, true, "CELL.XBIN for cell is not set.\n"); - good = false; - } - if (xBin->data.S32[index] < *xBinMin) { - *xBinMin = xBin->data.S32[index]; - } - yBin->data.S32[index] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YBIN"); - if (!mdok || yBin->data.S32[index] == 0) { - psError(PS_ERR_UNKNOWN, true, "CELL.YBIN for cell is not set.\n"); - good = false; - } - if (yBin->data.S32[index] < *yBinMin) { - *yBinMin = yBin->data.S32[index]; - } - - // Do we need to flip? - int xParityCell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XPARITY"); - if (!mdok || (xParityCell != 1 && xParityCell != -1)) { - psError(PS_ERR_UNKNOWN, true, "CELL.XPARITY for cell is not set.\n"); - good = false; - } - int yParityCell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YPARITY"); - if (!mdok || (yParityCell != 1 && yParityCell != -1)) { - psError(PS_ERR_UNKNOWN, true, "CELL.YPARITY for cell is not set.\n"); - good = false; - } - - // Parity of the chip on the FPA - int xParityChip = 0, yParityChip = 0; - if (chipStuff) { - pmChip *chip = cell->parent; // The parent chip - xParityChip = psMetadataLookupS32(&mdok, chip->concepts, "CHIP.XPARITY"); - if (!mdok || (xParityChip != 1 && xParityChip != -1)) { - psError(PS_ERR_UNKNOWN, true, "CHIP.XPARITY for chip is not set.\n"); - good = false; - } - yParityChip = psMetadataLookupS32(&mdok, chip->concepts, "CHIP.YPARITY"); - if (!mdok || (yParityChip != 1 && yParityChip != -1)) { - psError(PS_ERR_UNKNOWN, true, "CHIP.YPARITY for chip is not set.\n"); - good = false; - } - } - - // Set the flips on the basis of the parity - if (xParityCell * xParityChip == -1) { - xFlip->data.U8[index] = 1; - } else { - xFlip->data.U8[index] = 0; - } - if (yParityCell * yParityChip == -1) { - yFlip->data.U8[index] = 1; - } else { - yFlip->data.U8[index] = 0; - } - - x0->data.S32[index] = x0Chip + x0Cell; - y0->data.S32[index] = y0Chip + y0Cell; - - // Add the readout to the array of images to be mosaicked - psArray *readouts = cell->readouts; // The array of readouts - if (readouts->n > 1) { - psLogMsg(__func__, PS_LOG_WARN, "Cell contains more than one readout (%ld) --- only the first will " - "be mosaicked.\n", readouts->n); - } - pmReadout *readout = readouts->data[0]; // The only readout we'll bother with - - // The images to put into the mosaic - images->data[index] = psMemIncrRefCounter(readout->image); - weights->data[index] = psMemIncrRefCounter(readout->weight); - masks->data[index] = psMemIncrRefCounter(readout->mask); - - psTrace(__func__, 9, "Added cell (%x) %ld: %d,%d; %d,%d, %d,%d.\n", cell, index, - x0->data.S32[index], y0->data.S32[index], xBin->data.S32[index], yBin->data.S32[index], - xFlip->data.U8[index], yFlip->data.U8[index]); - - return true; -} - - -// Mosaic together the cells in a chip -static 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 - const pmChip *chip // Chip to mosaic - ) -{ - assert(mosaicImage); - assert(mosaicMask); - assert(mosaicWeights); - assert(chip); - - psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked - psArray *weights = psArrayAlloc(0); // Array of weight images to be mosaicked - psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked - psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates - psVector *y0 = psVectorAlloc(0, PS_TYPE_S32); // Origin y coordinates - psVector *xBin = psVectorAlloc(0, PS_TYPE_S32); // Binning in x - psVector *yBin = psVectorAlloc(0, PS_TYPE_S32); // Binning in y - psVector *xFlip = psVectorAlloc(0, PS_TYPE_U8); // Flip in x? - psVector *yFlip = psVectorAlloc(0, PS_TYPE_U8); // Flip in y? - - // Binning for the mosaicked chip is the minimum binning allowed by the cells - *xBinChip = INT_MAX; - *yBinChip = INT_MAX; - - // Set up the required inputs - bool allGood = true; // Is everything good, well-behaved? - psArray *cells = chip->cells; // The array of cells - for (int i = 0; i < cells->n; i++) { - pmCell *cell = cells->data[i]; // The cell of interest - if (!cell) { - continue; - } - allGood |= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip, - cell, xBinChip, yBinChip, false); - } - - // Mosaic the images together and we're done - if (allGood) { - *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 allGood; -} - -// Mosaic together the cells in a FPA -static bool fpaMosaic(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 *xBinFPA, int *yBinFPA, // The binning in x and y, to be returned - const pmFPA *fpa // FPA to mosaic - ) -{ - assert(mosaicImage); - assert(mosaicMask); - assert(mosaicWeights); - assert(fpa); - - psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked - psArray *weights = psArrayAlloc(0); // Array of weight images to be mosaicked - psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked - psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates - psVector *y0 = psVectorAlloc(0, PS_TYPE_S32); // Origin y coordinates - psVector *xBin = psVectorAlloc(0, PS_TYPE_S32); // Binning in x - psVector *yBin = psVectorAlloc(0, PS_TYPE_S32); // Binning in y - psVector *xFlip = psVectorAlloc(0, PS_TYPE_U8); // Flip in x? - psVector *yFlip = psVectorAlloc(0, PS_TYPE_U8); // Flip in y? - - // Binning for the mosaicked chip is the minimum binning allowed by the cells - *xBinFPA = INT_MAX; - *yBinFPA = INT_MAX; - - // Set up the required inputs - bool allGood = true; // Is everything good, well-behaved? - psArray *chips = fpa->chips; // Array of chips - for (int i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // The chip of interest - if (!chip) { - continue; - } - psArray *cells = chip->cells; // The array of cells - for (int j = 0; j < cells->n; j++) { - pmCell *cell = cells->data[j]; // The cell of interest - if (!cell) { - continue; - } - allGood |= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip, - cell, xBinFPA, yBinFPA, true); - } - } - - // Mosaic the images together and we're done - if (allGood) { - *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0); - *mosaicWeights = imageMosaic(weights, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0); - *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0); - } - - // Clean up - psFree(images); - psFree(weights); - psFree(masks); - psFree(xFlip); - psFree(yFlip); - psFree(xBin); - psFree(yBin); - psFree(x0); - psFree(y0); - - return allGood; -} - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// 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 *target, // Target chip --- may contain only a single cell - const pmChip *source // Source chip whose cells will be mosaicked - ) -{ - // Target exists, and has only a single cell - PS_ASSERT_PTR_NON_NULL(target, false); - PS_ASSERT_PTR_NON_NULL(target->cells, false); - if (target->cells->n != 1) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Target chip for mosaicking must contain a single cell.\n"); - return false; - } - pmCell *targetCell = target->cells->data[0]; // The target cell - PS_ASSERT_PTR_NON_NULL(targetCell, false); - // Source exists - PS_ASSERT_PTR_NON_NULL(source, false); - - - 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, source))) { - // Case 1 --- we need only cut out the region - psTrace(__func__, 1, "Case 1 mosaicking: simple cut-out.\n"); - pmHDU *hdu = source->hdu; // The HDU that has the pixels - if (!hdu || !hdu->images) { - hdu = source->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, source)) { - psError(PS_ERR_UNKNOWN, 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 - } - - // Set the concepts for the target cell - psList *sourceCells = psArrayToList(source->cells); // List of cells - cellConcepts(targetCell, sourceCells, xBin, yBin, chipRegion); - psFree(sourceCells); - psFree(chipRegion); - - // Currently, there's nothing interesting in the chip concepts that needs to be updated. - - // Now make a new readout to go in the target cell - pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout - newReadout->image = mosaicImage; - newReadout->mask = mosaicMask; - newReadout->weight = mosaicWeights; - psFree(newReadout); // Drop reference - - // Update the headers - pmHDU *sourceHDU = pmHDUFromChip(source); // The HDU for the source - pmHDU *targetHDU = pmHDUFromChip(target); // The HDU for the target - targetHDU->header = psMetadataCopy(targetHDU->header, sourceHDU->header); - - return true; -} - - -// Same deal with the FPA -bool pmFPAMosaic(pmFPA *target, // Target FPA --- may contain only a single chip with a single cell - const pmFPA *source // FPA whose chips and cells will be mosaicked - ) -{ - // Target exists, and has only a single chip with single cell - PS_ASSERT_PTR_NON_NULL(target, false); - PS_ASSERT_PTR_NON_NULL(target->chips, false); - if (target->chips->n != 1) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Target FPA for mosaicking must contain a single chip.\n"); - return false; - } - pmChip *targetChip = target->chips->data[0]; // The target chip - PS_ASSERT_PTR_NON_NULL(targetChip, false); - PS_ASSERT_PTR_NON_NULL(targetChip->cells, false); - if (target->chips->n != 1) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Target FPA for mosaicking must contain a single cell.\n"); - return false; - } - pmCell *targetCell = targetChip->cells->data[0]; // The target cell - PS_ASSERT_PTR_NON_NULL(targetCell, false); - // Source exists - PS_ASSERT_PTR_NON_NULL(source, false); - - psImage *mosaicImage = NULL; // The mosaic image - psImage *mosaicMask = NULL; // The mosaic mask - psImage *mosaicWeights = NULL; // The mosaic weights - - // Find the HDU - psRegion *fpaRegion = NULL; // Region on the HDU that corresponds to the FPA - int xBin = 0, yBin = 0; // Binning for the FPA mosaic - if ((fpaRegion = niceFPA(&xBin, &yBin, source))) { - // Case 1 --- we need only cut out the region - psTrace(__func__, 1, "Case 1 mosaicking: simple cut-out.\n"); - pmHDU *hdu = source->hdu; // The HDU that has the pixels - mosaicImage = psMemIncrRefCounter(psImageSubset(hdu->images->data[0], *fpaRegion)); - if (hdu->masks) { - mosaicMask = psMemIncrRefCounter(psImageSubset(hdu->masks->data[0], *fpaRegion)); - } - if (hdu->weights) { - mosaicWeights = psMemIncrRefCounter(psImageSubset(hdu->weights->data[0], *fpaRegion)); - } - } else { - // Case 2 --- we need to mosaic by cut and paste - psTrace(__func__, 1, "Case 2 mosaicking: cut and paste.\n"); - if (!fpaMosaic(&mosaicImage, &mosaicMask, &mosaicWeights, &xBin, &yBin, source)) { - psError(PS_ERR_UNKNOWN, false, "Unable to mosaic chips.\n"); - return false; - } - fpaRegion = psRegionAlloc(NAN, NAN, NAN, NAN); // We've cut and paste, so there's no valid trimsec - } - - // Set the concepts for the target cell, and add the mosaic in - // First we need a list of cells - psList *sourceCells = psListAlloc(NULL); // List of source cells - psArray *chips = source->chips; // Array of chips - for (long i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // Chip of interest - if (!chip) { - continue; - } - psArray *cells = chip->cells; - for (long j = 0; j < cells->n; j++) { - pmCell *cell = cells->data[j]; // Cell of interest - if (!cell) { - continue; - } - psListAdd(sourceCells, PS_LIST_TAIL, cell); - } - } - cellConcepts(targetCell, sourceCells, xBin, yBin, fpaRegion); - psFree(sourceCells); - psFree(fpaRegion); - - // Currently, there's nothing interesting in the chip concepts that needs to be updated. - - // Copy the concepts for the target FPA - target->concepts = psMetadataCopy(target->concepts, source->concepts); - - // Now make a new readout to go in the new cell - pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout - newReadout->image = mosaicImage; - newReadout->mask = mosaicMask; - newReadout->weight = mosaicWeights; - psFree(newReadout); // Drop reference - - // Update the headers - pmHDU *sourceHDU = pmHDUFromFPA(source); // The HDU for the source - pmHDU *targetHDU = pmHDUFromFPA(target); // The HDU for the target - targetHDU->header = psMetadataCopy(targetHDU->header, sourceHDU->header); - - return true; -} - Index: unk/psModules/src/camera/pmChipMosaic.h =================================================================== --- /trunk/psModules/src/camera/pmChipMosaic.h (revision 7486) +++ (revision ) @@ -1,16 +1,0 @@ -#ifndef PM_CHIP_MOSAIC_H -#define PM_CHIP_MOSAIC_H - -#include "pmFPA.h" - -// Mosaic all cells within a chip -bool pmChipMosaic(pmChip *target, // Target chip --- may contain only a single cell - const pmChip *source // Source chip whose cells will be mosaicked - ); - -// Mosaic all cells within an FPA -bool pmFPAMosaic(pmFPA *target, // Target FPA --- may contain only a single chip with a single cell - const pmFPA *source // FPA whose chips and cells will be mosaicked - ); - -#endif Index: /trunk/psModules/src/camera/pmFPAMosaic.c =================================================================== --- /trunk/psModules/src/camera/pmFPAMosaic.c (revision 7487) +++ /trunk/psModules/src/camera/pmFPAMosaic.c (revision 7487) @@ -0,0 +1,1101 @@ +#include +#include +#include +#include "pmFPA.h" +#include "pmHDU.h" +#include "pmHDUUtils.h" +#include "pmChipMosaic.h" + +#define CELL_LIST_BUFFER 50 // Buffer size for cell lists + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// 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; \ +} + + +// Get the bounds for an chip's pixels on the HDU +static bool chipBounds(psRegion *bounds, // The bounds for the chip + const pmChip *chip // The chip to examine for contiguity + ) +{ + assert(chip); + + 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 || psRegionIsNaN(*trimsec)) { + psError(PS_ERR_UNKNOWN, true, "CELL.TRIMSEC hasn't been set for cell %d.\n", i); + return false; + } + + 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; + } + } + + return true; +} + +// Make sure the TRIMSEC doesn't overlap with the established image bounds +static bool chipContiguousTrimsec(psRegion *bounds, // The bounds of the image, altered if primary==true + const pmChip *chip // The chip to examine for contiguity + ) +{ + assert(bounds); + assert(chip); + + 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 || psRegionIsNaN(*trimsec)) { + psError(PS_ERR_UNKNOWN, true, "CELL.TRIMSEC hasn't been set for cell %d.\n", i); + return false; + } + + if (REGIONS_OVERLAP(trimsec, bounds)) { + return false; + } + } + + return true; +} + +// Make sure the BIASSEC doesn't overlap with the established image bounds +static bool chipContiguousBiassec(psRegion *bounds, // The bounds of the image, altered if primary==true + const pmChip *chip // The chip to examine for contiguity + ) +{ + assert(bounds); + assert(chip); + + // Check that the biases don't get in the way + 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 + psList *biassecs = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.BIASSEC"); // Bias sections + if (!mdok || !biassecs) { + psError(PS_ERR_UNKNOWN, 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 (psRegionIsNaN(*biassec)) { + continue; + } + if (REGIONS_OVERLAP(biassec, bounds)) { + psFree(biassecsIter); + return false; + } + } + psFree(biassecsIter); + } + + // If we've gotten this far, everything is fine. + return true; +} + +// Are the pixels for the FPA contiguous on the HDU? +// Work this out by examining all the CELL.TRIMSEC and CELL.BIASSEC regions for the component cells +static bool fpaContiguous(psRegion *bounds, // The bounds of the image, returned + const pmFPA *fpa // The FPA to examine for contiguity + ) +{ + assert(bounds); + assert(fpa); + + *bounds = psRegionSet(INFINITY, 0, INFINITY, 0); + + // Get the size of the pixels on the HDU + psArray *chips = fpa->chips; // The array of chips + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // Chip of interest + if (!chipBounds(bounds, chip)) { + return false; + } + } + + // Make sure the bias regions don't get in the way of the HDU + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // Chip of interest + if (!chipContiguousBiassec(bounds, chip)) { + return false; + } + } + + // If we got through it all, they must all be contiguous + return true; +} + + + +// Check a cell for niceness in the parity and binning +static bool niceCellParityBinning(int *xBin, int *yBin, // Binning for cell, to be returned + const pmCell *cell // Cell to check for niceness + ) +{ + assert(xBin); + assert(yBin); + assert(cell); + + // A "nice" cell must have only a single readout + if (cell->readouts->n != 1) { + return false; + } + + // A "nice" cell must have parity == 1 + bool mdok = true; // Status of MD lookup + int xParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XPARITY"); // Parity in x + if (!mdok || xParity == 0) { + psError(PS_ERR_UNKNOWN, true, "CELL.XPARITY hasn't been set for cell.\n"); + return false; + } + if (xParity != 1) { + return false; + } + int yParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YPARITY"); // Parity in y + if (!mdok || yParity == 0) { + psError(PS_ERR_UNKNOWN, true, "CELL.YPARITY hasn't been set for cell.\n"); + return false; + } + if (yParity != 1) { + return false; + } + + // A "nice" cell must have consistent binning + int xBinCell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XBIN"); // Binning in x + if (!mdok || xBin <= 0) { + psError(PS_ERR_UNKNOWN, true, "CELL.XPARITY hasn't been set for cell.\n"); + return false; + } + int yBinCell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YBIN"); // Binning in y + if (!mdok || yBin <= 0) { + psError(PS_ERR_UNKNOWN, true, "CELL.YPARITY hasn't been set for cell.\n"); + return false; + } + if (*xBin == 0 || *yBin == 0) { + *xBin = xBinCell; + *yBin = yBinCell; + } else if (xBinCell != *xBin || yBinCell != *yBin) { + return false; + } + + return true; +} + + +// Check a cell for niceness in the boundaries +static bool niceCellBounds(const pmCell *cell, // Cell to check for niceness + const psRegion *imageBounds // Bounds of the image on the HDU + ) +{ + // A "nice" cell must have the (0,0) pixel at CELL.X0,CELL.Y0 + bool mdok = true; // Status of MD lookup + int x0 = psMetadataLookupS32(&mdok, cell->concepts, "CELL.X0"); // Position of (0,0) on chip + if (!mdok) { + psError(PS_ERR_UNKNOWN, true, "CELL.X0 hasn't been set for cell.\n"); + return false; + } + int y0 = psMetadataLookupS32(&mdok, cell->concepts, "CELL.Y0"); // Position of (0,0) on chip + if (!mdok) { + psError(PS_ERR_UNKNOWN, true, "CELL.Y0 hasn't been set for cell.\n"); + return false; + } + pmReadout *readout = cell->readouts->data[0]; // A representative readout + if (!readout) { + return false; // 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); + return false; + } + + 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 + const pmChip *chip // Chip to examine for "niceness". + ) +{ + assert(xBinChip); + assert(yBinChip); + assert(chip); + + // 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 + *xBinChip = 0; + *yBinChip = 0; + for (int i = 0; i < chip->cells->n; i++) { + pmCell *cell = chip->cells->data[i]; // The cell of interest + if (!niceCellParityBinning(xBinChip, yBinChip, cell)) { + return NULL; + } + } + + // Now check that the pixels are all contiguous + psRegion *imageBounds = psRegionAlloc(INFINITY, 0, INFINITY, 0); // Bound of image on HDU + if (!chipBounds(imageBounds, chip) || !chipContiguousBiassec(imageBounds, chip)) { + 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 + if (!niceCellBounds(cell, imageBounds)) { + 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 (!chipContiguousTrimsec(imageBounds, testChip) || + !chipContiguousBiassec(imageBounds, testChip)) { + psTrace(__func__, 5, "Image isn't contiguous.\n"); + psFree(imageBounds); + return NULL; + } + } + } + + return imageBounds; +} + +// Is the FPA "nice"? If so, return the region containing the FPA pixels +static psRegion *niceFPA(int *xBinFPA, int *yBinFPA, // Binning for FPA, to be returned + const pmFPA *fpa // FPA to examine for "niceness". + ) +{ + assert(xBinFPA); + assert(yBinFPA); + assert(fpa); + + // Check that we've got the HDU in the chip or the FPA + if (!fpa->hdu || !fpa->hdu->images) { + return NULL; + } + + // Check parity and binning for component cells + *xBinFPA = 0; + *yBinFPA = 0; + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; // The chip of interest + for (int j = 0; i < chip->cells->n; i++) { + pmCell *cell = chip->cells->data[j]; // The cell of interest + if (!niceCellParityBinning(xBinFPA, yBinFPA, cell)) { + return NULL; + } + } + } + + // Now check that the pixels are all contiguous + psRegion *imageBounds = psRegionAlloc(0, 0, 0, 0); // Bound of image on HDU + if (!fpaContiguous(imageBounds, fpa)) { + 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 < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; // The chip of interest + for (int j = 0; i < chip->cells->n; i++) { + pmCell *cell = chip->cells->data[j]; // The cell of interest + if (!niceCellBounds(cell, imageBounds)) { + 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 + psTrace(__func__, 3, "Mosaicking %ld cells.\n", source->n); + 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] - xMin, y0->data.S32[i] - yMin, "+"); + } 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 - yMin) / + 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 - xMin) / + 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 + psList *sources, // Source cells + int xBin, int yBin, // Binning + psRegion *trimsec // The trim section + ) +{ + assert(target); + assert(sources); + assert(xBin > 0 && yBin > 0); + assert(trimsec); + + 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 readdir = 0; // Cell read direction + + int nCells = 0; // Number of cells; + psListIterator *sourcesIter = psListIteratorAlloc(sources, PS_LIST_HEAD, false); // Iterator for sources + pmCell *cell = NULL; // Source cell from iteration + while ((cell = psListGetAndIncrement(sourcesIter))) { + 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 (nCells == 1) { + timeSys = psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS"); + readdir = psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR"); + } 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; + } + if (readdir != psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR")) { + psLogMsg(__func__, PS_LOG_ERROR, "Differing cell read directions in use: %d vs %d\n", readdir, + psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR")); + 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; + } + } + psFree(sourcesIter); + + 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); + MD_UPDATE(target->concepts, "CELL.READDIR", S32, readdir); + + // 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; +} + + +// Add a cell and its various properties to the arrays +static bool addCell(psArray *images, // Array of images + psArray *masks, // Array of masks + psArray *weights, // Array of weights + psVector *x0, // Array of X0 + psVector *y0, // Array of Y0 + psVector *xBin, // Array of XBIN + psVector *yBin, // Array of YBIN + psVector *xFlip, // Array indicating whether x axis should be flipped + psVector *yFlip, // Array indicating whether y axis should be flipped + const pmCell *cell, // Cell to add + int *xBinMin, // The minimum x binning, returned + int *yBinMin, // The minimum y binning, returned + bool chipStuff // Worry about chip stuff as well? + ) +{ + if (!cell) { + return false; + } + + // Expand the arrays and vectors to handle new data + long index = images->n; // The index to use + if (images->n == images->nalloc) { + images = psArrayRealloc(images, index + CELL_LIST_BUFFER); + masks = psArrayRealloc(masks, index + CELL_LIST_BUFFER); + weights = psArrayRealloc(weights, index + CELL_LIST_BUFFER); + x0 = psVectorRealloc(x0, index+ CELL_LIST_BUFFER); + y0 = psVectorRealloc(y0, index+ CELL_LIST_BUFFER); + xBin = psVectorRealloc(xBin, index+ CELL_LIST_BUFFER); + yBin = psVectorRealloc(yBin, index+ CELL_LIST_BUFFER); + xFlip = psVectorRealloc(xFlip, index+ CELL_LIST_BUFFER); + yFlip = psVectorRealloc(yFlip, index+ CELL_LIST_BUFFER); + } + + images->n = index + 1; + masks->n = index + 1; + weights->n = index + 1; + x0->n = index + 1; + y0->n = index + 1; + xBin->n = index + 1; + yBin->n = index + 1; + xFlip->n = index + 1; + yFlip->n = index + 1; + + bool mdok = true; // Status of MD lookup + bool good = true; // Is everything good? + + // Offset of the cell on the chip + int x0Cell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.X0"); + if (!mdok) { + psError(PS_ERR_UNKNOWN, true, "CELL.X0 for cell is not set.\n"); + good = false; + } + int y0Cell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.Y0"); + if (!mdok) { + psError(PS_ERR_UNKNOWN, true, "CELL.Y0 for cell is not set.\n"); + good = false; + } + psTrace(__func__, 5, "Cell %d: x0=%d y0=%d\n", index, x0Cell, y0Cell); + + // Offset of the chip on the FPA + int x0Chip = 0, y0Chip = 0; + if (chipStuff) { + pmChip *chip = cell->parent; // The parent chip + if (!chip) { + psError(PS_ERR_UNKNOWN, true, "Cell has no parent chip --- can't find CHIP.X0 and CHIP.Y0\n"); + good = false; + } + x0Chip = psMetadataLookupS32(&mdok, chip->concepts, "CHIP.X0"); + if (!mdok) { + psError(PS_ERR_UNKNOWN, true, "CHIP.X0 for chip is not set.\n"); + good = false; + } + y0Chip = psMetadataLookupS32(&mdok, chip->concepts, "CHIP.Y0"); + if (!mdok) { + psError(PS_ERR_UNKNOWN, true, "CHIP.Y0 for chip is not set.\n"); + good = false; + } + } + + // Binning + xBin->data.S32[index] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XBIN"); + if (!mdok || xBin->data.S32[index] == 0) { + psError(PS_ERR_UNKNOWN, true, "CELL.XBIN for cell is not set.\n"); + good = false; + } + if (xBin->data.S32[index] < *xBinMin) { + *xBinMin = xBin->data.S32[index]; + } + yBin->data.S32[index] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YBIN"); + if (!mdok || yBin->data.S32[index] == 0) { + psError(PS_ERR_UNKNOWN, true, "CELL.YBIN for cell is not set.\n"); + good = false; + } + if (yBin->data.S32[index] < *yBinMin) { + *yBinMin = yBin->data.S32[index]; + } + + // Do we need to flip? + int xParityCell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XPARITY"); + if (!mdok || (xParityCell != 1 && xParityCell != -1)) { + psError(PS_ERR_UNKNOWN, true, "CELL.XPARITY for cell is not set.\n"); + good = false; + } + int yParityCell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YPARITY"); + if (!mdok || (yParityCell != 1 && yParityCell != -1)) { + psError(PS_ERR_UNKNOWN, true, "CELL.YPARITY for cell is not set.\n"); + good = false; + } + + // Parity of the chip on the FPA + int xParityChip = 0, yParityChip = 0; + if (chipStuff) { + pmChip *chip = cell->parent; // The parent chip + xParityChip = psMetadataLookupS32(&mdok, chip->concepts, "CHIP.XPARITY"); + if (!mdok || (xParityChip != 1 && xParityChip != -1)) { + psError(PS_ERR_UNKNOWN, true, "CHIP.XPARITY for chip is not set.\n"); + good = false; + } + yParityChip = psMetadataLookupS32(&mdok, chip->concepts, "CHIP.YPARITY"); + if (!mdok || (yParityChip != 1 && yParityChip != -1)) { + psError(PS_ERR_UNKNOWN, true, "CHIP.YPARITY for chip is not set.\n"); + good = false; + } + } + + // Set the flips on the basis of the parity + if (xParityCell * xParityChip == -1) { + xFlip->data.U8[index] = 1; + } else { + xFlip->data.U8[index] = 0; + } + if (yParityCell * yParityChip == -1) { + yFlip->data.U8[index] = 1; + } else { + yFlip->data.U8[index] = 0; + } + + x0->data.S32[index] = x0Chip + x0Cell; + y0->data.S32[index] = y0Chip + y0Cell; + + // Add the readout to the array of images to be mosaicked + psArray *readouts = cell->readouts; // The array of readouts + if (readouts->n > 1) { + psLogMsg(__func__, PS_LOG_WARN, "Cell contains more than one readout (%ld) --- only the first will " + "be mosaicked.\n", readouts->n); + } + pmReadout *readout = readouts->data[0]; // The only readout we'll bother with + + // The images to put into the mosaic + images->data[index] = psMemIncrRefCounter(readout->image); + weights->data[index] = psMemIncrRefCounter(readout->weight); + masks->data[index] = psMemIncrRefCounter(readout->mask); + + psTrace(__func__, 9, "Added cell (%x) %ld: %d,%d; %d,%d, %d,%d.\n", cell, index, + x0->data.S32[index], y0->data.S32[index], xBin->data.S32[index], yBin->data.S32[index], + xFlip->data.U8[index], yFlip->data.U8[index]); + + return true; +} + + +// Mosaic together the cells in a chip +static 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 + const pmChip *chip // Chip to mosaic + ) +{ + assert(mosaicImage); + assert(mosaicMask); + assert(mosaicWeights); + assert(chip); + + psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked + psArray *weights = psArrayAlloc(0); // Array of weight images to be mosaicked + psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked + psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates + psVector *y0 = psVectorAlloc(0, PS_TYPE_S32); // Origin y coordinates + psVector *xBin = psVectorAlloc(0, PS_TYPE_S32); // Binning in x + psVector *yBin = psVectorAlloc(0, PS_TYPE_S32); // Binning in y + psVector *xFlip = psVectorAlloc(0, PS_TYPE_U8); // Flip in x? + psVector *yFlip = psVectorAlloc(0, PS_TYPE_U8); // Flip in y? + + // Binning for the mosaicked chip is the minimum binning allowed by the cells + *xBinChip = INT_MAX; + *yBinChip = INT_MAX; + + // Set up the required inputs + bool allGood = true; // Is everything good, well-behaved? + psArray *cells = chip->cells; // The array of cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // The cell of interest + if (!cell) { + continue; + } + allGood |= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip, + cell, xBinChip, yBinChip, false); + } + + // Mosaic the images together and we're done + if (allGood) { + *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 allGood; +} + +// Mosaic together the cells in a FPA +static bool fpaMosaic(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 *xBinFPA, int *yBinFPA, // The binning in x and y, to be returned + const pmFPA *fpa // FPA to mosaic + ) +{ + assert(mosaicImage); + assert(mosaicMask); + assert(mosaicWeights); + assert(fpa); + + psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked + psArray *weights = psArrayAlloc(0); // Array of weight images to be mosaicked + psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked + psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates + psVector *y0 = psVectorAlloc(0, PS_TYPE_S32); // Origin y coordinates + psVector *xBin = psVectorAlloc(0, PS_TYPE_S32); // Binning in x + psVector *yBin = psVectorAlloc(0, PS_TYPE_S32); // Binning in y + psVector *xFlip = psVectorAlloc(0, PS_TYPE_U8); // Flip in x? + psVector *yFlip = psVectorAlloc(0, PS_TYPE_U8); // Flip in y? + + // Binning for the mosaicked chip is the minimum binning allowed by the cells + *xBinFPA = INT_MAX; + *yBinFPA = INT_MAX; + + // Set up the required inputs + bool allGood = true; // Is everything good, well-behaved? + psArray *chips = fpa->chips; // Array of chips + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // The chip of interest + if (!chip) { + continue; + } + psArray *cells = chip->cells; // The array of cells + for (int j = 0; j < cells->n; j++) { + pmCell *cell = cells->data[j]; // The cell of interest + if (!cell) { + continue; + } + allGood |= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip, + cell, xBinFPA, yBinFPA, true); + } + } + + // Mosaic the images together and we're done + if (allGood) { + *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0); + *mosaicWeights = imageMosaic(weights, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0); + *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0); + } + + // Clean up + psFree(images); + psFree(weights); + psFree(masks); + psFree(xFlip); + psFree(yFlip); + psFree(xBin); + psFree(yBin); + psFree(x0); + psFree(y0); + + return allGood; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// 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 *target, // Target chip --- may contain only a single cell + const pmChip *source // Source chip whose cells will be mosaicked + ) +{ + // Target exists, and has only a single cell + PS_ASSERT_PTR_NON_NULL(target, false); + PS_ASSERT_PTR_NON_NULL(target->cells, false); + if (target->cells->n != 1) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Target chip for mosaicking must contain a single cell.\n"); + return false; + } + pmCell *targetCell = target->cells->data[0]; // The target cell + PS_ASSERT_PTR_NON_NULL(targetCell, false); + // Source exists + PS_ASSERT_PTR_NON_NULL(source, false); + + + 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, source))) { + // Case 1 --- we need only cut out the region + psTrace(__func__, 1, "Case 1 mosaicking: simple cut-out.\n"); + pmHDU *hdu = source->hdu; // The HDU that has the pixels + if (!hdu || !hdu->images) { + hdu = source->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, source)) { + psError(PS_ERR_UNKNOWN, 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 + } + + // Set the concepts for the target cell + psList *sourceCells = psArrayToList(source->cells); // List of cells + cellConcepts(targetCell, sourceCells, xBin, yBin, chipRegion); + psFree(sourceCells); + psFree(chipRegion); + + // Currently, there's nothing interesting in the chip concepts that needs to be updated. + + // Now make a new readout to go in the target cell + pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout + newReadout->image = mosaicImage; + newReadout->mask = mosaicMask; + newReadout->weight = mosaicWeights; + psFree(newReadout); // Drop reference + + // Update the headers + pmHDU *sourceHDU = pmHDUFromChip(source); // The HDU for the source + pmHDU *targetHDU = pmHDUFromChip(target); // The HDU for the target + targetHDU->header = psMetadataCopy(targetHDU->header, sourceHDU->header); + + return true; +} + + +// Same deal with the FPA +bool pmFPAMosaic(pmFPA *target, // Target FPA --- may contain only a single chip with a single cell + const pmFPA *source // FPA whose chips and cells will be mosaicked + ) +{ + // Target exists, and has only a single chip with single cell + PS_ASSERT_PTR_NON_NULL(target, false); + PS_ASSERT_PTR_NON_NULL(target->chips, false); + if (target->chips->n != 1) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Target FPA for mosaicking must contain a single chip.\n"); + return false; + } + pmChip *targetChip = target->chips->data[0]; // The target chip + PS_ASSERT_PTR_NON_NULL(targetChip, false); + PS_ASSERT_PTR_NON_NULL(targetChip->cells, false); + if (target->chips->n != 1) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Target FPA for mosaicking must contain a single cell.\n"); + return false; + } + pmCell *targetCell = targetChip->cells->data[0]; // The target cell + PS_ASSERT_PTR_NON_NULL(targetCell, false); + // Source exists + PS_ASSERT_PTR_NON_NULL(source, false); + + psImage *mosaicImage = NULL; // The mosaic image + psImage *mosaicMask = NULL; // The mosaic mask + psImage *mosaicWeights = NULL; // The mosaic weights + + // Find the HDU + psRegion *fpaRegion = NULL; // Region on the HDU that corresponds to the FPA + int xBin = 0, yBin = 0; // Binning for the FPA mosaic + if ((fpaRegion = niceFPA(&xBin, &yBin, source))) { + // Case 1 --- we need only cut out the region + psTrace(__func__, 1, "Case 1 mosaicking: simple cut-out.\n"); + pmHDU *hdu = source->hdu; // The HDU that has the pixels + mosaicImage = psMemIncrRefCounter(psImageSubset(hdu->images->data[0], *fpaRegion)); + if (hdu->masks) { + mosaicMask = psMemIncrRefCounter(psImageSubset(hdu->masks->data[0], *fpaRegion)); + } + if (hdu->weights) { + mosaicWeights = psMemIncrRefCounter(psImageSubset(hdu->weights->data[0], *fpaRegion)); + } + } else { + // Case 2 --- we need to mosaic by cut and paste + psTrace(__func__, 1, "Case 2 mosaicking: cut and paste.\n"); + if (!fpaMosaic(&mosaicImage, &mosaicMask, &mosaicWeights, &xBin, &yBin, source)) { + psError(PS_ERR_UNKNOWN, false, "Unable to mosaic chips.\n"); + return false; + } + fpaRegion = psRegionAlloc(NAN, NAN, NAN, NAN); // We've cut and paste, so there's no valid trimsec + } + + // Set the concepts for the target cell, and add the mosaic in + // First we need a list of cells + psList *sourceCells = psListAlloc(NULL); // List of source cells + psArray *chips = source->chips; // Array of chips + for (long i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // Chip of interest + if (!chip) { + continue; + } + psArray *cells = chip->cells; + for (long j = 0; j < cells->n; j++) { + pmCell *cell = cells->data[j]; // Cell of interest + if (!cell) { + continue; + } + psListAdd(sourceCells, PS_LIST_TAIL, cell); + } + } + cellConcepts(targetCell, sourceCells, xBin, yBin, fpaRegion); + psFree(sourceCells); + psFree(fpaRegion); + + // Currently, there's nothing interesting in the chip concepts that needs to be updated. + + // Copy the concepts for the target FPA + target->concepts = psMetadataCopy(target->concepts, source->concepts); + + // Now make a new readout to go in the new cell + pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout + newReadout->image = mosaicImage; + newReadout->mask = mosaicMask; + newReadout->weight = mosaicWeights; + psFree(newReadout); // Drop reference + + // Update the headers + pmHDU *sourceHDU = pmHDUFromFPA(source); // The HDU for the source + pmHDU *targetHDU = pmHDUFromFPA(target); // The HDU for the target + targetHDU->header = psMetadataCopy(targetHDU->header, sourceHDU->header); + + return true; +} + Index: /trunk/psModules/src/camera/pmFPAMosaic.h =================================================================== --- /trunk/psModules/src/camera/pmFPAMosaic.h (revision 7487) +++ /trunk/psModules/src/camera/pmFPAMosaic.h (revision 7487) @@ -0,0 +1,16 @@ +#ifndef PM_CHIP_MOSAIC_H +#define PM_CHIP_MOSAIC_H + +#include "pmFPA.h" + +// Mosaic all cells within a chip +bool pmChipMosaic(pmChip *target, // Target chip --- may contain only a single cell + const pmChip *source // Source chip whose cells will be mosaicked + ); + +// Mosaic all cells within an FPA +bool pmFPAMosaic(pmFPA *target, // Target FPA --- may contain only a single chip with a single cell + const pmFPA *source // FPA whose chips and cells will be mosaicked + ); + +#endif Index: /trunk/psModules/src/psmodules.h =================================================================== --- /trunk/psModules/src/psmodules.h (revision 7486) +++ /trunk/psModules/src/psmodules.h (revision 7487) @@ -25,8 +25,8 @@ #include #include +#include #include #include #include -#include // the following headers are from psModule:astrom