Index: trunk/archive/scripts/src/phase2/pmChipMosaic.c =================================================================== --- trunk/archive/scripts/src/phase2/pmChipMosaic.c (revision 5621) +++ trunk/archive/scripts/src/phase2/pmChipMosaic.c (revision 5786) @@ -1,3 +1,4 @@ #include +#include #include "pslib.h" @@ -14,16 +15,20 @@ } -psImage *pmChipMosaic(pmChip *chip, // Chip to mosaic - int xBinChip, int yBinChip // Binning of mosaic image in x and y +// Mosaic multiple images, with flips, binning and offsets +psImage *p_pmImageMosaic(const psArray *source, // Images to splice in + const psVector *xFlip, const psVector *yFlip, // Need to flip x and y? + const psVector *xBinSource, const psVector *yBinSource, // Binning in x and y of + // source images + int xBinTarget, int yBinTarget, // Binning in x and y of target images + const psVector *x0, const psVector *y0 // Offsets for source images on target ) { - psArray *cells = chip->cells; // The array of cells - psVector *x0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin x coordinates - psVector *y0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin y coordinates - psVector *xBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in x - psVector *yBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in y - psVector *xParity = psVectorAlloc(cells->n, PS_TYPE_S32); // Parity in x - psVector *yParity = psVectorAlloc(cells->n, PS_TYPE_S32); // Parity in y - psArray *trimsec = psArrayAlloc(cells->n); // Trim section + assert(source); + assert(xFlip && xFlip->type.type == PS_TYPE_U8); + assert(yFlip && yFlip->type.type == PS_TYPE_U8); + assert(xBinSource && xBinSource->type.type == PS_TYPE_S32); + assert(yBinSource && yBinSource->type.type == PS_TYPE_S32); + assert(x0 && x0->type.type == PS_TYPE_S32); + assert(y0 && y0->type.type == PS_TYPE_S32); // Get the maximum extent of the mosaic image @@ -32,93 +37,138 @@ int yMin = INT_MAX; int yMax = 0; - for (int i = 0; i < cells->n; i++) { - pmCell *cell = cells->data[i]; // The cell of interest - x0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); - y0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); - xBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); - yBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); - xParity->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); - yParity->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); - trimsec->data[i] = psMemIncrRefCounter(psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC")); - psTrace(__func__, 7, "Cell %d trimsec: [%.0f:%.0f,%.0f:%.0f]\n", i, ((psRegion*)trimsec->data[i])->x0, - ((psRegion*)trimsec->data[i])->x1, ((psRegion*)trimsec->data[i])->y0, - ((psRegion*)trimsec->data[i])->y1); - // Size of cell in x and y - int nx = (int)(((psRegion*)trimsec->data[i])->x1 - (int)((psRegion*)trimsec->data[i])->x0); - int ny = (int)(((psRegion*)trimsec->data[i])->y1 - (int)((psRegion*)trimsec->data[i])->y0); - psTrace(__func__, 5, "Extent of cell %d: %d -> %d , %d -> %d\n", i, x0->data.S32[i], - x0->data.S32[i] + xParity->data.S32[i] * xBin->data.S32[i] * nx, y0->data.S32[i], - y0->data.S32[i] + yParity->data.S32[i] * yBin->data.S32[i] * ny); + for (int i = 0; i < source->n; i++) { + psImage *image = source->data[i]; // The image of interest - 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->data.S32[i] * xBin->data.S32[i] * nx - xParity->data.S32[i], xMin, xMax); - COMPARE(y0->data.S32[i] + yParity->data.S32[i] * yBin->data.S32[i] * ny - yParity->data.S32[i], yMin, yMax); + assert(image->type.type == PS_TYPE_F32); // Only implemented for F32 images so far. + + // Size of cell in x and y + int xParity = xFlip->data.U8[i] ? -1 : 1; + int yParity = yFlip->data.U8[i] ? -1 : 1; + psTrace(__func__, 5, "Extent of cell %d: %d -> %d , %d -> %d\n", i, x0->data.S32[i], + x0->data.S32[i] + xParity * xBinSource->data.S32[i] * image->numCols, y0->data.S32[i], + y0->data.S32[i] + yParity * yBinSource->data.S32[i] * image->numRows); + + COMPARE(x0->data.S32[i], xMin, xMax); + COMPARE(y0->data.S32[i], yMin, yMax); + // Subtract the parity to get the inclusive limit (not exclusive) + COMPARE(x0->data.S32[i] + xParity * xBinSource->data.S32[i] * image->numCols - xParity, xMin, xMax); + COMPARE(y0->data.S32[i] + yParity * yBinSource->data.S32[i] * image->numRows - yParity, yMin, yMax); } // Set up the image // Since both upper and lower values are inclusive, we need to add one to the size - float xSize = (float)(xMax - xMin + 1) / (float)xBinChip; + float xSize = (float)(xMax - xMin + 1) / (float)xBinTarget; if (xSize - (int)xSize > 0) { - xSize += 1; + xSize += 1; } - float ySize = (float)(yMax - yMin + 1) / (float)yBinChip; + float ySize = (float)(yMax - yMin + 1) / (float)yBinTarget; if (ySize - (int)ySize > 0) { - ySize += 1; + ySize += 1; } - psTrace(__func__, 3, "Mosaicked chip will be %dx%d\n", (int)xSize, (int)ySize); + psTrace(__func__, 3, "Spliced image will be %dx%d\n", (int)xSize, (int)ySize); psImage *mosaic = psImageAlloc((int)xSize, (int)ySize, PS_TYPE_F32); // The mosaic image psImageInit(mosaic, 0.0); - // Set the image + // Next pass through the images to do the mosaicking + for (int i = 0; i < source->n; i++) { + psImage *image = source->data[i]; // The image of interest + if (xBinSource->data.S32[i] == xBinTarget && yBinSource->data.S32[i] == yBinTarget && + xFlip->data.U8[i] == 0 && yFlip->data.U8[i] == 0) { + // Let someone else do the hard work; useful to test psImageOverlaySection if no other reason + psImageOverlaySection(mosaic, image, x0->data.S32[i], y0->data.S32[i], "+"); + } else { + // We have to do the hard work ourself + for (int y = 0; y < image->numRows; y++) { + int yParity = yFlip->data.U8[i] ? -1 : 1; + float yTargetBase = (y0->data.S32[i] + yParity * yBinSource->data.S32[i] * y) / yBinTarget; + for (int x = 0; x < image->numCols; x++) { + int xParity = xFlip->data.U8[i] ? -1 : 1; + float xTargetBase = (x0->data.S32[i] + xParity * xBinSource->data.S32[i] * x) / + xBinTarget; + + // In case the original image is binned but the mosaic is not, we need to fill in the + // values in the mosaic. + for (int j = 0; j < yBinSource->data.S32[i]; j++) { + int yTarget = (int)(yTargetBase + yParity * (float)j / (float)yBinTarget); + for (int i = 0; i < xBinSource->data.S32[i]; i++) { + int xTarget = (int)(xTargetBase + xParity * (float)i / (float)xBinTarget); + + mosaic->data.F32[yTarget][xTarget] += image->data.F32[y][x]; + } + } // Iterating over mosaic image for binned input image + } + } // Iterating over input image + } + } + + return mosaic; +} + +psImage *pmChipMosaic(pmChip *chip, // Chip to mosaic + int xBinChip, int yBinChip // Binning of mosaic image in x and y + ) +{ + psArray *cells = chip->cells; // The array of cells + psArray *images = psArrayAlloc(cells->n); // Array of images that will be mosaicked + psVector *x0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin x coordinates + psVector *y0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin y coordinates + psVector *xBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in x + psVector *yBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in y + psVector *xFlip = psVectorAlloc(cells->n, PS_TYPE_U8); // Flip in x? + psVector *yFlip = psVectorAlloc(cells->n, PS_TYPE_U8); // Flip in y? + + // Set up the required inputs for (int i = 0; i < cells->n; i++) { - pmCell *cell = cells->data[i]; // The cell of interest - psArray *readouts = cell->readouts; // The array of readouts - if (readouts->n > 1) { - psLogMsg(__func__, PS_LOG_WARN, "Cell %d contains more than one readout --- only the first will " - "be mosaicked.\n", i); - } - psImage *image = ((pmReadout*)readouts->data[0])->image; // The image to put into the mosaic - psImage *trimmed = psImageSubset(image, *(psRegion*)(trimsec->data[i])); // Trimmed image (no overscan) + pmCell *cell = cells->data[i]; // The cell of interest + x0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); + y0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + xBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); + yBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); + int xParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); + int yParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); + if (xParity == 1) { + xFlip->data.U8[i] = 0; + } else if (xParity == -1) { + xFlip->data.U8[i] = 1; + } else { + psLogMsg(__func__, PS_LOG_WARN, "The x parity of cell %d is not +/- 1 (it's %d) --- " + "assuming +1.\n", i, xParity); + xFlip->data.U8[i] = 0; + } + if (yParity == 1) { + yFlip->data.U8[i] = 0; + } else if (yParity == -1) { + yFlip->data.U8[i] = 1; + } else { + psLogMsg(__func__, PS_LOG_WARN, "The y parity of cell %d is not +/- 1 (it's %d) --- " + "assuming +1.\n", i, yParity); + yFlip->data.U8[i] = 0; + } - if (xBin->data.S32[i] == xBinChip && yBin->data.S32[i] == yBinChip && xParity->data.S32[i] == 1 && - yParity->data.S32[i] == 1) { - // Let someone else do the hard work; useful to test psImageOverlaySection if no other reason - psImageOverlaySection(mosaic, trimmed, x0->data.S32[i], y0->data.S32[i], "+"); - } else { - // We have to do the hard work ourself - for (int y = 0; y < trimmed->numRows; y++) { - float yTargetBase = (y0->data.S32[i] + yParity->data.S32[i] * yBin->data.S32[i] * y) / - yBinChip; - for (int x = 0; x < trimmed->numCols; x++) { - float xTargetBase = (x0->data.S32[i] + xParity->data.S32[i] * xBin->data.S32[i] * x) / - xBinChip; + // Trim the image to get rid of the overscan + psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); + psTrace(__func__, 7, "Cell %d trimsec: [%.0f:%.0f,%.0f:%.0f]\n", i, trimsec->x0, trimsec->x1, + trimsec->y0, trimsec->y1); + psArray *readouts = cell->readouts; // The array of readouts + if (readouts->n > 1) { + psLogMsg(__func__, PS_LOG_WARN, "Cell %d contains more than one readout --- only the first will " + "be mosaicked.\n", i); + } + psImage *image = ((pmReadout*)readouts->data[0])->image; // The image to put into the mosaic + images->data[i] = psImageSubset(image, *trimsec); // Trimmed image + } - // In case the original image is binned but the mosaic is not, we need to fill in the - // values in the mosaic. - for (int j = 0; j < yBin->data.S32[i]; j++) { - int yTarget = (int)(yTargetBase + yParity->data.S32[i] * (float)j / (float)yBinChip); - for (int i = 0; i < xBin->data.S32[i]; i++) { - int xTarget = (int)(xTargetBase + - xParity->data.S32[i] * (float)i / (float)xBinChip); + // Mosaic the images together and we're done + psImage *mosaic = p_pmImageMosaic(images, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0); - mosaic->data.F32[yTarget][xTarget] += trimmed->data.F32[y][x]; - } - } // Iterating over mosaic image for binned input image - } - } // Iterating over input image - } - psFree(trimmed); - } // Iterating over cells - + // Clean up psFree(x0); psFree(y0); psFree(xBin); psFree(yBin); - psFree(xParity); - psFree(yParity); - psFree(trimsec); + psFree(xFlip); + psFree(yFlip); + psFree(images); return mosaic;