Index: /branches/rel10_ifa/psModules/src/astrom/pmFPACopy.c =================================================================== --- /branches/rel10_ifa/psModules/src/astrom/pmFPACopy.c (revision 6682) +++ /branches/rel10_ifa/psModules/src/astrom/pmFPACopy.c (revision 6682) @@ -0,0 +1,522 @@ +#include +#include + +#include "pslib.h" +#include "pmFPA.h" +#include "pmFPAUtils.h" +#include "pmHDU.h" + +#define MAX(x,y) ((x) > (y) ? (x) : (y)) + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// File-static functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Add cells in a chip to a list +static bool addCellsFromChip(psList *list, // List of cells + const pmChip *chip // The chip from which to add cells + ) +{ + assert(list); + assert(chip); + + psArray *cells = chip->cells; // Array of cells + bool result = true; // Result of adding cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // A cell + result |= psListAdd(list, PS_LIST_TAIL, cell); + } + + return result; +} + +// Add cells in an FPA to a list +static bool addCellsFromFPA(psList *list, // List of cells + const pmFPA *fpa // The FPA from which to add cells + ) +{ + assert(list); + assert(fpa); + + psArray *chips = fpa->chips; // Array of chips + bool result = true; // Result of adding cells + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // A chip + result |= addCellsFromChip(list, chip); + } + + return result; +} + +// Get a list of cells that share the HDU for the target cell +static bool cellList(psList *targets, // The list of target cells + psList *sources, // The list of source cells + const pmCell *targetCell, // The target cell + const pmCell *sourceCell // The source cell + ) +{ + assert(targetCell); + assert(sourceCell); + + if (targetCell->hdu) { + if (targets) { + psListAdd(targets, PS_LIST_TAIL, targetCell); + } + if (sources) { + psListAdd(sources, PS_LIST_TAIL, sourceCell); + } + } else { + pmChip *targetChip = targetCell->parent; // The target parent chip + pmChip *sourceChip = sourceCell->parent; // The source parent chip + if (targetChip->hdu) { + if (targets) { + addCellsFromChip(targets, targetChip); + } + if (sources) { + addCellsFromChip(sources, sourceChip); + } + } else { + pmFPA *targetFPA = targetChip->parent; // The target parent FPA + pmFPA *sourceFPA = sourceChip->parent; // The source parent FPA + if (targetFPA->hdu) { + if (targets) { + addCellsFromFPA(targets, targetFPA); + } + if (sources) { + addCellsFromFPA(sources, sourceFPA); + } + } else { + psError(PS_ERR_IO, true, "Unable to find HDU for cell to generate list!\n"); + return false; + } + } + } + + return list; +} + +// Get the maximum extent of the HDU from the trimsec and biassecs +static bool sizeHDU(int *xSize, int *ySize, // Size of HDU + psList *cells // List of cells + ) +{ + psListIterator *cellsIter = psListIteratorAlloc(cells, PS_LIST_HEAD, false); // Iterator for cells + pmCell *cell = NULL; // The cell from iteration + bool mdok = true; // Status of MD lookup + *xSize = 0; + *ySize = 0; + while ((cell = psListGetAndIncrement(cellsIter))) { + psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section + if (mdok && trimsec && !psRegionIsBad(*trimsec)) { + *xSize = MAX(trimsec->x1, *xSize); + *ySize = MAX(trimsec->y1, *ySize); + } + psList *biassecs = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.BIASSEC"); // Bias sections + if (mdok && biassecs) { + psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); // Iterator + psRegion *biassec = NULL; // The bias section + while ((biassec = psListGetAndIncrement(biassecsIter))) { + if (!psRegionIsBad(*trimsec)) { + *xSize = MAX(biassec->x1, *xSize); + *ySize = MAX(biassec->y1, *ySize); + } + } + psFree(biassecsIter); + } + } + psFree(cellsIter); + + return (*xSize == 0 || *ySize == 0); +} + + +static psRegion sectionForImage(int *position, // Position on the output image, updated + const psImage *image, // Image containing the sizes and offsets + int readdir // Read direction, 1=rows, 2=cols + ) +{ + psRegion region; + switch (readdir) { + case 1: // Read direction is rows + region = psRegionSet(*position, image->row0, *position + image->numCols, + image->row0 + image->numRows); + *position += image->numCols; + break; + case 2: // Read direction is columns + region = psRegionSet(image->col0, *position, image->col0 + image->numCols, + *position + image->numRows); + *position += image->numRows; + break; + default: + psAbort(__func__, "Shouldn't ever get here!\n"); + } + + return region; +} + +static bool doBiasSections(int *position, // Position on the output image, updated + pmCell *target, // Target cell + const pmCell *source // Source cell + ) +{ + psMetadataItem *biassecItem = psMetadataLookup(target->concepts, "CELL.BIASSEC"); // Bias sections + if (!biassecItem) { + psLogMsg(__func__, PS_LOG_WARN, "CELL.BIASSEC has not been initialised in target cell --- " + "ignored.\n"); + return false; + } + psFree(biassecItem->data.V); // Blow away the old list + psList *biassecs = psListAlloc(NULL); + biassecItem->data.V = biassecs; + + bool mdok = true; // Status of MD lookup + int readdir = psMetadataLookupS32(&mdok, source->concepts, "CELL.READDIR"); // Read direction + if (!mdok || (readdir != 1 && readdir != 2)) { + // Probably unnecessary, but just in case.... + psLogMsg(__func__, PS_LOG_WARN, "CELL.READDIR is not set in source cell --- ignored.\n"); + return false; + } + + pmReadout *readout = source->readouts->data[0]; // The first source readout, as representative + psList *biases = readout->bias; // The bias images from the source readout + + psListIterator *biasIter = psListIteratorAlloc(biases, PS_LIST_HEAD, true); // Iterator for biases + psImage *bias = NULL; // Bias image from iteration + while ((bias = psListGetAndIncrement(biasIter))) { + // Construct a region + psRegion *biassec = psAlloc(sizeof(psRegion)); // The new region; need a psMemBlock + *biassec = sectionForImage(position, bias, readdir); + psListAdd(biassecs, PS_LIST_TAIL, biassec); + psFree(biassec); // Drop reference + } + psFree(biasIter); + + return true; +} + +// Generate CELL.TRIMSEC and CELL.BIASSEC for the target cells +static bool generateTrimBias(psList *targets, // List of target cells + psList *sources // List of source cells + ) +{ + pmCell *target = NULL, *source = NULL; // Cells from iteration + int numCells = targets->n; // Number of cells + int cellNum = 0; // The cell number + int position; // Position on the image + bool mdok = true; // Status of MD lookup + + // First run through to do the LHS biases + psListIterator *targetsIter = psListIteratorAlloc(targets, PS_LIST_HEAD, false); // Iterator for targets + psListIterator *sourcesIter = psListIteratorAlloc(sources, PS_LIST_HEAD, false); // Iterator for sources + bool done = false; // Done with iteration? + while ((target = psListGetAndIncrement(targetsIter)) && (source = psListGetAndIncrement(sourcesIter)) && + !done) { + if (cellNum < numCells/2) { + doBiasSections(&position, target, source); + } else { + done = true; + } + } + + // Second run through to do the trim sections + psListIteratorSet(targetsIter, PS_LIST_HEAD); + psListIteratorSet(sourcesIter, PS_LIST_HEAD); + while ((target = psListGetAndIncrement(targetsIter)) && (source = psListGetAndIncrement(sourcesIter))) { + psRegion *trimsec = psMetadataLookupPtr(&mdok, target->concepts, "CELL.TRIMSEC"); // Trim section + if (!mdok || !trimsec) { + psLogMsg(PS_ERR_IO, true, "CELL.TRIMSEC has not been initialised in target cell --- ignored.\n"); + continue; + } + + int readdir = psMetadataLookupS32(&mdok, source->concepts, "CELL.READDIR"); // Read direction + if (!mdok || (readdir != 1 && readdir != 2)) { + // Probably unnecessary, but just in case.... + psLogMsg(__func__, PS_LOG_WARN, "CELL.READDIR is not set in source cell --- ignored.\n"); + continue; + } + + pmReadout *readout = source->readouts->data[0]; // The first source readout, as representative + psImage *image = readout->image;// The proper image + *trimsec = sectionForImage(position, bias, readdir); + } + + // A final run through to do the RHS biases + psListIteratorSet(targetsIter, numCells/2); + psListIteratorSet(sourcesIter, numCells/2); + while ((target = psListGetAndIncrement(targetsIter)) && (source = psListGetAndIncrement(sourcesIter))) { + doBiasSections(&position, target, source); + } + + // Clean up + psFree(targetsIter); + psFree(sourcesIter); + + return (position > 0); +} + + +// Generate an HDU with the pixels +static bool generateHDU(pmCell *target // The target cell + const pmCell *source // The source cell + ) +{ + // Get the HDU and a list of cells below it + pmHDU *hdu = pmHDUFromCell(target); // The HDU in the target cell + psList *targetCells = psListAlloc(NULL); // List of target cells below the target HDU + psList *sourceCells = psListAlloc(NULL); // List of source cells below the target HDU + if (! cellList(targetCells, sources, target, source)) { + psError(PS_ERR_IO, true, "Unable to find cells to generate HDU!\n"); + return false; + } + + // Check the number of readouts + int numReadouts = -1; // Number of readouts + { + psListIterator *iter = psListIteratorAlloc(sourceCells, PS_LIST_HEAD, false); // Iterator for cells + pmCell *cell = NULL; // The cell from iteration + while ((cell = psListGetAndIncrement(iter))) + { + psArray *readouts = cell->readouts; + if (numReadouts == -1) { + numReadouts = readouts->n; + } else if (readouts->n != numReadouts) { + psError(PS_ERR_IO, true, "Number of readouts doesn't match: %d vs %d\n", readouts->n, + numReadouts); + return false; + } + + } + psFree(cellsIter); + } + + // Get the size of the HDU, either from existing trimsec and biassec, or generate these and try again + int xSize = 0, ySize = 0; // Size of HDU + if (!sizeHDU(&xSize, &ySize, targetCells) || !(generateTrimBias(targetCells, sourceCells) && + sizeHDU(&xSize, &ySize, targetCells))) { + psError(PS_ERR_IO, true, "Unable to determine size of HDU!\n"); + return false; + } + + hdu->images = psArrayAlloc(numReadouts); + for (int i = 0; i < numReadouts; i++) { + psImage *image = psImageAlloc(xSize, ySize, PS_TYPE_F32); + psImageInit(image, 0.0); + hdu->images->data[i] = image; + } + + return true; +} + + +#define UPDATE_CASE(TYPENAME,TYPE) \ +case TYPENAME: targetItem->data.TYPE = sourceItem->data.TYPE; \ +break; + +// Update a particular "concept" for the target cell from the source cell +static bool updateConcept(pmCell *target, // Target cell + const pmCell *source, // Source cell + const char *concept // Concept name + ) +{ + psMetadataItem *targetItem = psMetadataLookup(target->concepts, concept); // Concept from the target + psMetadataItem *sourceItem = psMetadataLookup(source->concepts, concept); // Concept from the source + + switch (targetItem->type) { + UPDATE_CASE(PS_TYPE_S32, S32); + UPDATE_CASE(PS_TYPE_F32, F32); + UPDATE_CASE(PS_TYPE_F64, F64); + UPDATE_CASE(PS_DATA_STR, V); + default: + psLogMsg(__func__, PS_LOG_WARN, "Unsupported type (%x) for concept %s --- ignored.\n", + targetItem->type, concept); + return false; + } + return true; +} + + +// Copy pixels from a target image to a source image, with flips +static bool copyPixels(psImage *target, // Target image (HDU pixels) + psImage *source, // Source image (from source cell) + psRegion region, // Region for pasting + bool xFlip, // Flip in x? + bool yFlip // Flip in y? + ) +{ + psImage *overlay = psMemIncrRefCounter(source); + if (xFlip) { + psImage *temp = psImageFlipX(overlay); + psFree(overlay); + overlay = temp; + } + if (yFlip) { + psImage *temp = psImageFlipY(overlay); + psFree(overlay); + overlay = temp; + } + int numPix = psImageOverlaySection(target, overlay, region.x0, region.y0, "+"); + psFree(overlay); + return (numPix > 0); +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + +int pmFPACopy(pmFPA *target, // The target FPA + pmFPA *source // The source FPA, to be copied + ) +{ + assert(target); + assert(source); + + psArray *targetChips = target->chips; // The target chips + psArray *sourceChips = source->chips; // The source chips + if (targetChips->n != sourceChips->n) { + psError(PS_ERR_IO, true, "Number of source chips (%d) differs from the number of target chips (%d)\n", + sourceChips->n, targetChips->n); + return false; + } + + int numChips = 0; // Number of chips copied + for (int i = 0; i < targetChips->n; i++) { + pmChip *targetChip = targetChips->data[i]; // The target chip + const char *chipName = psMetadataLookupStr(NULL, targetChip->concepts, "CHIP.NAME"); // Name of chip + int chipNum = pmFPAFindChip(source, chipName); // Number of chip with that name + if (chipNum >= 0) { + pmChip *sourceChip = sourceChips->data[chipNum]; // The source chip + int numCells = pmChipCopy(targetChip, sourceChip); // Number of cells copied + psTrace(__func__, 5, "Copied %d cells for chip %s\n", numCells, chipName); + numChips++; + } + } + + return numChips; +} + +int pmChipCopy(pmChip *target, // The target chip + pmChip *source // The source chip, to be copied + ) +{ + assert(target); + assert(source); + + psArray *targetCells = target->cells; // The target cells + psArray *sourceCells = source->cells; // The source cells + if (targetCells->n != sourceCells->n) { + psError(PS_ERR_IO, true, "Number of source cells (%d) differs from the number of target cells (%d)\n", + sourceCells->n, targetCells->n); + return false; + } + + int numCells = 0; // Number of cells copied + for (int i = 0; i < targetCells->n; i++) { + pmCell *targetCell = targetCells->data[i]; // The target cell + const char *cellName = psMetadataLookupStr(NULL, targetCell->concepts, "CELL.NAME"); // Name of cell + int cellNum = pmFPAFindCell(source, cellName); // Number of cell with that name + if (cellNum >= 0) { + pmCell *sourceCell = sourceCells->data[cellNum]; // The source cell + int numReadouts = pmCellCopy(targetCell, sourceCell); // Number of readouts copied + psTrace(__func__, 5, "Copied %d readouts for cell %s\n", numReadouts, cellName); + numCells++; + } + } + + return numCells; + +} + +int pmCellCopy(pmCell *target, // The target cell + pmCell *source // The source cell, to be copied + ) +{ + assert(target); + assert(source); + + psArray *targetReadouts = target->readouts; // The target readouts + psArray *sourceReadouts = source->readouts; // The source readouts + if (targetReadouts->n != sourceReadouts->n) { + psError(PS_ERR_IO, true, "Number of source readouts (%d) differs from the number of target readouts " + "(%d)\n", sourceReadouts->n, targetReadouts->n); + return false; + } + int numReadouts = targetReadouts->n; // Number of readouts copied + + pmHDU *hdu = pmHDUFromCell(target); // The target HDU; we need to fix this up + if (! hdu->images) { + generateHDU(cell); + } + + // Need to check CELL.X0, CELL.Y0, CELL.XBIN, CELL.YBIN should be unchanged; we can check this by + // setting the values and leaving it to the pmConceptWrite functions to do the checking against the + // camera configuration when it comes time to write them out. + updateConcept(target, source, "CELL.X0"); + updateConcept(target, source, "CELL.Y0"); + updateConcept(target, source, "CELL.XBIN"); + updateConcept(target, source, "CELL.YBIN"); + + // Need to check/change CELL.XPARITY and CELL.YPARITY + bool xFlip = (psMetadataLookupS32(NULL, target->concepts, "CELL.XPARITY") != + psMetadataLookupS32(NULL, target->concepts, "CELL.XPARITY")); // Switch parity in x? + bool yFlip = (psMetadataLookupS32(NULL, target->concepts, "CELL.YPARITY") != + psMetadataLookupS32(NULL, target->concepts, "CELL.YPARITY")); // Switch parity in y? + + bool mdok = true; // Status of MD lookup + psRegion *trimsec = psMetadataLookupPtr(&mdok, target->concepts, "CELL.TRIMSEC"); // The trim section + if (!mdok || !trimsec || psRegionIsBad(*trimsec)) { + psError(PS_ERR_IO, true, "CELL.TRIMSEC isn't set!\n"); + return 0; + } + psList *biassecs = psMetadataLookupPtr(&mdok, target->concepts, "CELL.BIASSEC"); // The bias sections + if (!mdok || !biassecs) { + psError(PS_ERR_IO, true, "CELL.BIASSEC isn't set!\n"); + return 0; + } + + for (int i = 0; i < targetReadouts->n; i++) { + pmReadout *sourceReadout = sourceReadouts->data[i]; // The source readout + psImage *sourceImage = sourceReadout->image; // The source image + pmReadout *targetReadout = targetReadouts->data[i]; // The target readout + if (sourceImage->numCols != trimsec->x1 - trimsec->x0 || + sourceImage->numRows != trimsec->y1 - trimsec->y0) { + psString trimsecString = psRegionToString(*trimsec); // String with the trim section + psLogMsg(__func__, PS_LOG_WARN, "Source image size (%dx%d) for readout %d doesn't match " + "CELL.TRIMSEC for target (%s) -- ignored.\n", sourceImage->numCols, sourceImage->numRows, + i, trimsecString); + psFree(trimsecString); + } else { + copyPixels(hdu->images->data[i], sourceImage, *trimsec, xFlip, yFlip); + targetReadout->image = psMemIncrRefCounter(psImageSubset(hdu->images->data[i], *trimsec)); + } + + psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); // Iterator + psRegion *biassec = NULL; // Bias section from iteration + while ((biassec = psListGetAndIncrement(biassecsIter))) { + if (psRegionIsBad(*biassec)) { + psString *biassecString = psRegionToString(biassec); // String for bias section + psLogMsg(__func__, PS_LOG_WARN, "Bias section (%s) isn't set --- ignored.\n", biassecString); + psFree(biassecString); + continue; + } + if (sourceImage->numCols != biassec->x1 - biassec->x0 || + sourceImage->numRows != biassec->y1 - biassec->y0) { + psString biassecString = psRegionToString(*biassec); // String with the bias section + psLogMsg(__func__, PS_LOG_WARN, "Source image size (%dx%d) for readout %d doesn't match " + "CELL.BIASSEC for target (%s) -- ignored.\n", sourceImage->numCols, + sourceImage->numRows, i, biassecString); + psFree(biassecString); + } else { + copyPixels(hdu->images->data[i], sourceImage, *biassec, xFlip, yFlip); + psImage *bias = psMemIncrRefCounter(psImageSubset(hdu->images->data[i], *biassec)); + psListAdd(targetReadout->bias, PS_LIST_TAIL, bias); + psFree(bias); // Drop reference + } + } + psFree(biassecsIter); + + } + + return numReadouts; +}