Index: /trunk/psModules/src/camera/pmChipMosaic.c =================================================================== --- /trunk/psModules/src/camera/pmChipMosaic.c (revision 7453) +++ /trunk/psModules/src/camera/pmChipMosaic.c (revision 7454) @@ -4,5 +4,7 @@ #include "pmFPA.h" #include "pmHDU.h" - +#include "pmChipMosaic.h" + +#define CELL_LIST_BUFFER 50 // Buffer size for cell lists ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -33,17 +35,12 @@ } -// 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? - ) -{ - assert(bounds); + +// 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); - if (primary) { - *bounds = psRegionSet(INFINITY, 0, INFINITY, 0); - } psArray *cells = chip->cells; // The array of cells bool mdok = true; // Status of MD lookup @@ -52,28 +49,67 @@ psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section if (!mdok || !trimsec || psRegionIsNaN(*trimsec)) { - psError(PS_ERR_IO, true, "CELL.TRIMSEC hasn't been set for cell %d.\n", i); + psError(PS_ERR_UNKNOWN, 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)) { + 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_IO, true, "CELL.BIASSEC hasn't been set for cell %d.\n", i); + psError(PS_ERR_UNKNOWN, true, "CELL.BIASSEC hasn't been set for cell %d.\n", i); return false; } @@ -100,8 +136,129 @@ } +// 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 - pmChip *chip // Chip to examine for "niceness". + const pmChip *chip // Chip to examine for "niceness". ) { @@ -116,55 +273,16 @@ // 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) { + if (!niceCellParityBinning(xBinChip, yBinChip, cell)) { 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)) { + 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); @@ -178,25 +296,5 @@ 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); + if (!niceCellBounds(cell, imageBounds)) { psFree(imageBounds); return NULL; @@ -214,6 +312,60 @@ continue; } - if (!chipContiguous(imageBounds, testChip, false)) { + 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; @@ -259,4 +411,5 @@ 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 @@ -315,13 +468,14 @@ xFlip->data.U8[i] == 0 && yFlip->data.U8[i] == 0) { // Let someone else do the hard work; useful to test psImageOverlaySection if no other reason - psImageOverlaySection(mosaic, image, x0->data.S32[i], y0->data.S32[i], "+"); + 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) / yBinTarget; + 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) / + float xTargetBase = (x0->data.S32[i] + xParity * xBinSource->data.S32[i] * x - xMin) / xBinTarget; @@ -359,5 +513,5 @@ // Set the concepts in the new cell, based on the values in the old one static bool cellConcepts(pmCell *target,// Target cell - psArray *sources, // Source cells + psList *sources, // Source cells int xBin, int yBin, // Binning psRegion *trimsec // The trim section @@ -378,8 +532,10 @@ double time = 0.0; // Time of observation psTimeType timeSys = 0; // Time system + int readdir = 0; // Cell read direction int nCells = 0; // Number of cells; - for (int i = 0; i < sources->n; i++) { - pmCell *cell = sources->data[i];// The cell of interest + 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; @@ -393,10 +549,18 @@ psTime *cellTime = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TIME"); time += psTimeToMJD(cellTime); - if (i == 0) { + if (nCells == 1) { timeSys = psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS"); - } else if (timeSys != psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS")) { - psLogMsg(__func__, PS_LOG_ERROR, "Differing time systems in use: %d vs %d\n", timeSys, - psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS")); - success = false; + 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; + } } @@ -410,4 +574,6 @@ } } + psFree(sourcesIter); + gain /= (float)nCells; readnoise /= (float)nCells; @@ -429,4 +595,5 @@ 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 @@ -448,11 +615,173 @@ +// 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 -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 - ) +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); @@ -461,17 +790,13 @@ assert(chip); - psArray *cells = chip->cells; // The array of cells - int numCells = cells->n; // Number of cells - psArray *images = psArrayAlloc(numCells); // Array of images that will be mosaicked - psArray *weights = psArrayAlloc(numCells); // Array of weight images to be mosaicked - psArray *masks = psArrayAlloc(numCells); // Array of mask images to be mosaicked - psVector *x0 = psVectorAlloc(numCells, PS_TYPE_S32); // Origin x coordinates - psVector *y0 = psVectorAlloc(numCells, PS_TYPE_S32); // Origin y coordinates - psVector *xBin = psVectorAlloc(numCells, PS_TYPE_S32); // Binning in x - psVector *yBin = psVectorAlloc(numCells, PS_TYPE_S32); // Binning in y - psVector *xFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Flip in x? - psVector *yFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Flip in y? - images->n = weights->n = masks->n = numCells; - x0->n = y0->n = xBin->n = yBin->n = xFlip->n = yFlip->n = numCells; + 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 @@ -480,75 +805,17 @@ // Set up the required inputs - psTrace(__func__, 1, "Mosaicking %d cells...\n", numCells); - bool good = true; // Is everything good, well-behaved? - for (int i = 0; i < numCells && good; i++) { + 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; } - bool mdok = true; // Status of MD lookup - x0->data.S32[i] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.X0"); - if (!mdok) { - psError(PS_ERR_IO, true, "CELL.X0 for cell %d is not set.\n", i); - good = false; - } - y0->data.S32[i] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.Y0"); - if (!mdok) { - psError(PS_ERR_IO, true, "CELL.Y0 for cell %d is not set.\n", i); - good = false; - } - psTrace(__func__, 5, "Cell %d: x0=%d y0=%d\n", i, x0->data.S32[i], y0->data.S32[i]); - xBin->data.S32[i] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XBIN"); - if (!mdok || xBin->data.S32[i] == 0) { - psError(PS_ERR_IO, true, "CELL.XBIN for cell %d is not set.\n", i); - good = false; - } - if (xBin->data.S32[i] < *xBinChip) { - *xBinChip = xBin->data.S32[i]; - } - yBin->data.S32[i] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YBIN"); - if (!mdok || yBin->data.S32[i] == 0) { - psError(PS_ERR_IO, true, "CELL.YBIN for cell %d is not set.\n", i); - good = false; - } - if (yBin->data.S32[i] < *yBinChip) { - *yBinChip = yBin->data.S32[i]; - } - int xParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XPARITY"); - if (!mdok || (xParity != 1 && xParity != -1)) { - psError(PS_ERR_IO, true, "CELL.XPARITY for cell %d is not set.\n", i); - good = false; - } - int yParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YPARITY"); - if (!mdok || (yParity != 1 && yParity != -1)) { - psError(PS_ERR_IO, true, "CELL.YPARITY for cell %d is not set.\n", i); - good = false; - } - if (xParity == -1) { - xFlip->data.U8[i] = 1; - } else { - xFlip->data.U8[i] = 0; - } - if (yParity == -1) { - yFlip->data.U8[i] = 1; - } else { - yFlip->data.U8[i] = 0; - } - - psArray *readouts = cell->readouts; // The array of readouts - if (readouts->n > 1) { - psLogMsg(__func__, PS_LOG_WARN, "Cell %d contains more than one readout --- only the first will " - "be mosaicked.\n", i); - } - pmReadout *readout = readouts->data[0]; // The only readout we'll bother with - - // The images to put into the mosaic - images->data[i] = psMemIncrRefCounter(readout->image); - weights->data[i] = psMemIncrRefCounter(readout->weight); - masks->data[i] = psMemIncrRefCounter(readout->mask); + allGood |= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip, + cell, xBinChip, yBinChip, false); } // Mosaic the images together and we're done - if (good) { + if (allGood) { *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0); *mosaicWeights = imageMosaic(weights, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0); @@ -567,5 +834,72 @@ psFree(y0); - return good; + 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; } @@ -599,8 +933,20 @@ // 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 +bool pmChipMosaic(pmChip *target, // Target chip --- may contain only a single cell + const pmChip *source // Source chip whose cells will be mosaicked ) { - PS_ASSERT_PTR_NON_NULL(chip, false); + // 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 @@ -611,14 +957,14 @@ 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))) { + 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 = chip->hdu; // The HDU that has the pixels + pmHDU *hdu = source->hdu; // The HDU that has the pixels if (!hdu || !hdu->images) { - hdu = chip->parent->hdu; - } - mosaicImage = psMemIncrRefCounter(psImageSubset(hdu->images->data[0], *chipRegion)); + hdu = source->parent->hdu; + } + mosaicImage = psMemIncrRefCounter(psImageSubset(hdu->images->data[0], *chipRegion)); if (hdu->masks) { - mosaicMask = psMemIncrRefCounter(psImageSubset(hdu->masks->data[0], *chipRegion)); + mosaicMask = psMemIncrRefCounter(psImageSubset(hdu->masks->data[0], *chipRegion)); } if (hdu->weights) { @@ -628,6 +974,6 @@ // 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"); + if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicWeights, &xBin, &yBin, source)) { + psError(PS_ERR_UNKNOWN, false, "Unable to mosaic cells.\n"); return false; } @@ -635,12 +981,12 @@ } - // 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); + // 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); - chip->mosaic = (struct pmCell*)newCell; - - // Now make a new readout to go in the new cell - pmReadout *newReadout = pmReadoutAlloc(newCell); // New readout + + // Now make a new readout to go in the target cell + pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout newReadout->image = mosaicImage; newReadout->mask = mosaicMask; @@ -650,2 +996,88 @@ 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); + + // 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 + + return true; +} + Index: /trunk/psModules/src/camera/pmChipMosaic.h =================================================================== --- /trunk/psModules/src/camera/pmChipMosaic.h (revision 7453) +++ /trunk/psModules/src/camera/pmChipMosaic.h (revision 7454) @@ -2,10 +2,15 @@ #define PM_CHIP_MOSAIC_H -#include "pslib.h" #include "pmFPA.h" // Mosaic all cells within a chip -bool pmChipMosaic(pmChip *chip // Chip whose cells will be mosaicked +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