Index: trunk/psModules/src/detrend/Makefile.am =================================================================== --- trunk/psModules/src/detrend/Makefile.am (revision 25930) +++ trunk/psModules/src/detrend/Makefile.am (revision 26893) @@ -17,5 +17,6 @@ pmDark.c \ pmRemnance.c \ - pmPattern.c + pmPattern.c \ + pmPatternIO.c # pmSkySubtract.c @@ -35,5 +36,6 @@ pmDark.h \ pmRemnance.h \ - pmPattern.h + pmPattern.h \ + pmPatternIO.h # pmSkySubtract.h Index: trunk/psModules/src/detrend/pmPattern.c =================================================================== --- trunk/psModules/src/detrend/pmPattern.c (revision 25930) +++ trunk/psModules/src/detrend/pmPattern.c (revision 26893) @@ -1,3 +1,4 @@ #include +#include #include @@ -9,5 +10,5 @@ // Mask a row as bad static void patternMaskRow(pmReadout *ro, // Readout to mask - int y, // Row to mask + int y, // Row to mask psImageMaskType bad // Mask value to give ) @@ -17,5 +18,5 @@ psAssert(y < image->numRows, "Row not in image"); - int numCols = image->numCols; // Size of image + int numCols = image->numCols; // Number of columns for (int x = 0; x < numCols; x++) { image->data.F32[y][x] = NAN; @@ -29,4 +30,8 @@ return; } + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Measurement and application +////////////////////////////////////////////////////////////////////////////////////////////////////////////// bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh, @@ -63,4 +68,5 @@ float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data + float background = stats->robustMedian; psFree(stats); psFree(rng); @@ -79,4 +85,7 @@ psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, order); // Polynomial to fit psVector *data = psVectorAlloc(numCols, PS_TYPE_F32); // Data to fit + + psImage *corr = psImageAlloc(order + 1, numRows, PS_TYPE_F64); // Corrections applied + psImageInit(corr, NAN); for (int y = 0; y < numRows; y++) { @@ -104,4 +113,7 @@ continue; } + + poly->coeff[0] -= background; + memcpy(corr->data.F64[y], poly->coeff, (order + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector if (!solution) { @@ -117,4 +129,8 @@ psFree(solution); } + + psMetadataAddImage(ro->analysis, PS_LIST_TAIL, PM_PATTERN_ROW_CORRECTION, PS_META_REPLACE, + "Pattern row correction", corr); + psFree(corr); psFree(indices); @@ -126,2 +142,239 @@ return true; } + +bool pmPatternRowApply(pmReadout *ro, psImageMaskType maskBad) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_IMAGE(ro, false); + + bool mdok; // Status of MD lookup + psImage *corr = psMetadataLookupPtr(&mdok, ro->analysis, PM_PATTERN_ROW_CORRECTION); // Correction + if (!mdok) { + // No correction to apply + return true; + } + + psImage *image = ro->image; // Image of interest + int numCols = image->numCols, numRows = image->numRows; // Size of image + + if (corr->numRows != numRows) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "Number of rows of image (%d) does not match number of rows of pattern correction (%d)\n", + numRows, corr->numRows); + return false; + } + + int order = corr->numCols - 1; // Polynomial order + psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, order); // Polynomial to apply + psVector *indices = psVectorAlloc(numCols, PS_TYPE_F32); // Indices for polynomial + float norm = 2.0 / (float)numCols; // Normalisation for indices + for (int x = 0; x < numCols; x++) { + indices->data.F32[x] = x * norm - 1.0; + } + + for (int y = 0; y < numRows; y++) { + memcpy(poly->coeff, corr->data.F64[y], (order + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector + if (!solution) { + psWarning("Unable to evaluate polynomial for row %d", y); + psErrorClear(); + patternMaskRow(ro, y, maskBad); + continue; + } + + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] -= solution->data.F32[x]; + } + psFree(solution); + } + + psFree(poly); + psFree(indices); + + return true; +} + + +bool pmPatternCell(pmChip *chip, const psVector *tweak, psStatsOptions bgStat, psStatsOptions cellStat, + psImageMaskType maskVal, psImageMaskType maskBad) +{ + PS_ASSERT_PTR_NON_NULL(chip, false); + PS_ASSERT_VECTOR_NON_NULL(tweak, false); + PS_ASSERT_VECTOR_SIZE(tweak, chip->cells->n, false); + PS_ASSERT_VECTOR_TYPE(tweak, PS_TYPE_U8, false); + + int numCells = tweak->n; // Number of cells + + psVector *mean = psVectorAlloc(numCells, PS_TYPE_F32); // Mean for each cell + psVector *meanMask = psVectorAlloc(numCells, PS_TYPE_VECTOR_MASK); // Mask for means + psVectorInit(mean, NAN); + psVectorInit(meanMask, 0); + + // Mask bits + enum { + PM_PATTERN_IGNORE = 0x01, // Ignore this cell + PM_PATTERN_TWEAK = 0x02, // Tweak this cell + PM_PATTERN_ERROR = 0x04, // Error in calculating background + PM_PATTERN_ALL = 0xFF, // All causes + }; + + // Count number of cells to tweak + int numTweak = 0; // Number of cells to tweak + int numIgnore = 0; // Number of cells to ignore + for (int i = 0; i < numCells; i++) { + pmCell *cell = chip->cells->data[i]; // Cell of interest + if (!cell || !cell->data_exists || !cell->process || + cell->readouts->n == 0 || cell->readouts->n > 1 || !cell->readouts->data[0]) { + numIgnore++; + meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PM_PATTERN_IGNORE; + continue; + } + if (tweak->data.U8[i]) { + meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PM_PATTERN_TWEAK; + numTweak++; + } + } + if (numTweak == 0) { + // Nothing to do + psFree(mean); + psFree(meanMask); + return true; + } + if (numTweak >= numCells - numIgnore) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot pattern-correct all cells within a chip."); + psFree(mean); + psFree(meanMask); + return false; + } + + // Measure mean of each cell + // This is not really the perfect thing to do, which would be to take a common mean for the set of cells + // which aren't being tweaked (because some cells will be heavily masked, so shouldn't be weighted the + // same as pure cells), but it's simple and fast. + psStats *bgStats = psStatsAlloc(bgStat); // Statistics on background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + for (int i = 0; i < numCells; i++) { + if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) { + continue; + } + pmCell *cell = chip->cells->data[i]; // Cell of interest + pmReadout *ro = cell->readouts->data[0]; // Readout of interest + + psStatsInit(bgStats); +#if 1 + if (!psImageBackground(bgStats, NULL, ro->image, ro->mask, maskVal, rng)) { + psWarning("Unable to measure background for cell %d\n", i); + psErrorClear(); + meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PM_PATTERN_ERROR; + continue; + } +#else + if (!psImageStats(bgStats, ro->image, ro->mask, maskVal)) { + psWarning("Unable to measure background for cell %d\n", i); + psErrorClear(); + meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PM_PATTERN_ERROR; + continue; + } +#endif + mean->data.F32[i] = psStatsGetValue(bgStats, bgStat); + if (!isfinite(mean->data.F32[i])) { + psWarning("Non-finite background for cell %d\n", i); + psErrorClear(); + meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PM_PATTERN_ERROR; + continue; + } + } + psFree(bgStats); + psFree(rng); + + psStats *cellStats = psStatsAlloc(cellStat); // Statistics on cells + if (!psVectorStats(cellStats, mean, NULL, meanMask, PM_PATTERN_ALL)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate mean cell background."); + psFree(mean); + psFree(meanMask); + psFree(cellStats); + return false; + } + + float background = psStatsGetValue(cellStats, cellStat); // Background value for chip + psFree(cellStats); + if (!isfinite(background)) { + psError(PS_ERR_UNKNOWN, false, "Non-finite mean cell background."); + psFree(mean); + psFree(meanMask); + return false; + } + + psLogMsg("psModules.detrend", PS_LOG_DETAIL, "Mean chip background is %f", background); + + for (int i = 0; i < numCells; i++) { + if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) { + continue; + } + if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_TWEAK)) { + continue; + } + pmCell *cell = chip->cells->data[i]; // Cell of interest + pmReadout *ro = cell->readouts->data[0]; // Readout of interest + if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_ERROR) { + psImageInit(ro->image, NAN); + psBinaryOp(ro->mask, ro->mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK)); + psMetadataAddF32(ro->analysis, PS_LIST_TAIL, PM_PATTERN_CELL_CORRECTION, PS_META_REPLACE, + "Pattern cell correction solution", NAN); + continue; + } + float correction = background - mean->data.F32[i]; // Correction to apply + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell + psLogMsg("psModules.detrend", PS_LOG_DETAIL, "Correcting background of cell %s by %f", + cellName, correction); + psBinaryOp(ro->image, ro->image, "+", psScalarAlloc(correction, PS_TYPE_F32)); + psMetadataAddF32(ro->analysis, PS_LIST_TAIL, PM_PATTERN_CELL_CORRECTION, PS_META_REPLACE, + "Pattern cell correction solution", correction); + } + + psFree(mean); + psFree(meanMask); + + return true; +} + +bool pmPatternCellApply(pmReadout *ro, psImageMaskType maskBad) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_IMAGE(ro, false); + + bool mdok; // Status of MD lookup + float corr = psMetadataLookupF32(&mdok, ro->analysis, PM_PATTERN_CELL_CORRECTION); // Correction to apply + if (!mdok) { + // No correction to apply + return true; + } + + psImage *image = ro->image, *mask = ro->mask; // Image and mask of interest + int numCols = image->numCols, numRows = image->numRows; // Size of image + + if (!isfinite(corr)) { + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] = NAN; + } + } + if (mask) { + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskBad; + } + } + } + } else { + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += corr; + } + } + } + + return true; +} + + Index: trunk/psModules/src/detrend/pmPattern.h =================================================================== --- trunk/psModules/src/detrend/pmPattern.h (revision 25930) +++ trunk/psModules/src/detrend/pmPattern.h (revision 26893) @@ -18,4 +18,7 @@ /// @{ +#define PM_PATTERN_ROW_CORRECTION "PATTERN.ROW.CORRECTION" // Pattern row correction on analysis metadata +#define PM_PATTERN_CELL_CORRECTION "PATTERN.CELL.CORRECTION" // Pattern cell correction on analysis metadata + /// Fit and remove pattern noise over rows bool pmPatternRow( @@ -31,4 +34,25 @@ ); +/// Apply previously measured row pattern correction +bool pmPatternRowApply(pmReadout *ro, ///< Readout to correct + psImageMaskType maskBad ///< Mask value to give bad pixels + ); + +/// Fix the background on cells known to be troublesome +bool pmPatternCell( + pmChip *chip, ///< Chip to correct + const psVector *tweak, ///< U8 vector indicating whether to tweak the corresponding cell + psStatsOptions bgStat, ///< Statistic to use for background measurement + psStatsOptions cellStat, ///< Statistic to use for combination of cell background measurements + psImageMaskType maskVal, ///< Mask value to use + psImageMaskType maskBad ///< Mask value to give bad pixels + ); + +/// Apply previously measured cell pattern correction +bool pmPatternCellApply(pmReadout *ro, ///< Readout to correct + psImageMaskType maskBad ///< Mask value to give bad pixels + ); + + /// @} #endif Index: trunk/psModules/src/detrend/pmPatternIO.c =================================================================== --- trunk/psModules/src/detrend/pmPatternIO.c (revision 26893) +++ trunk/psModules/src/detrend/pmPatternIO.c (revision 26893) @@ -0,0 +1,460 @@ +#include +#include +#include + +#include "pmHDU.h" +#include "pmHDUUtils.h" +#include "pmFPA.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" +#include "pmFPAfileFitsIO.h" +#include "pmFPAHeader.h" +#include "pmConceptsRead.h" +#include "pmConceptsWrite.h" + +#include "pmPattern.h" +#include "pmPatternIO.h" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +bool pmReadoutWritePattern(pmReadout *ro, psFits *fits) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PS_ASSERT_FITS_NON_NULL(fits, false); + + bool gotRow, gotCell; + psImage *rowCorr = psMetadataLookupPtr(&gotRow, ro->analysis, PM_PATTERN_ROW_CORRECTION); // Row correction + float cellCorr = psMetadataLookupF32(&gotCell, ro->analysis, PM_PATTERN_CELL_CORRECTION); // Cell corr. + + pmCell *cell = ro->parent; // Cell of interest + if (!cell) { + psError(PS_ERR_UNEXPECTED_NULL, true, "No cell associated with readout."); + return false; + } + pmChip *chip = cell->parent; // Chip of interest + pmFPA *fpa = chip->parent; // FPA of interest + pmHDU *hdu = pmHDUGetLowest(fpa, chip, cell); // HDU for readout + if (!hdu) { + psError(PS_ERR_UNEXPECTED_NULL, true, "No HDU associated with readout."); + return false; + } + + psMetadata *header = psMetadataCopy(NULL, hdu->header); // Header for output + + if (gotCell) { + psMetadataAddF32(header, PS_LIST_TAIL, PM_PATTERN_CELL_CORRECTION, PS_META_REPLACE, + "Pattern cell correction value", cellCorr); + } + + if (gotRow) { + if (!psFitsWriteImage(fits, header, rowCorr, 0, hdu->extname)) { + psError(PS_ERR_IO, false, "Unable to write pattern row correction."); + psFree(header); + return false; + } + } else { + if (!psFitsWriteBlank(fits, header, hdu->extname)) { + psError(PS_ERR_IO, false, "Unable to write pattern cell correction."); + psFree(header); + return false; + } + } + + psFree(header); + return true; +} + +static bool pmCellWritePattern(pmCell *cell, const pmFPAview *view, + pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(cell, false); + PS_ASSERT_PTR_NON_NULL(cell->readouts, false); + PS_ASSERT_PTR_NON_NULL(view, false); + + if (!pmConceptsWriteCell(cell, true, config)) { + psError(PS_ERR_IO, false, "Unable to write concepts for cell."); + return false; + } + + pmFPAview *thisView = pmFPAviewAlloc(view->nRows); // Copy of input view + *thisView = *view; + + for (int i = 0; i < cell->readouts->n; i++) { + pmReadout *readout = cell->readouts->data[i]; + thisView->readout = i; + if (!pmReadoutWritePattern(readout, file->fits)) { + psError(PS_ERR_IO, false, "Failed to write %dth readout", i); + psFree(thisView); + return false; + } + } + psFree(thisView); + return true; +} + +static bool pmChipWritePattern(pmChip *chip, const pmFPAview *view, + pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(chip, false); + PS_ASSERT_PTR_NON_NULL(chip->cells, false); + PS_ASSERT_PTR_NON_NULL(view, false); + + if (!pmConceptsWriteChip(chip, true, true, config)) { + psError(PS_ERR_IO, false, "Unable to write concepts for chip.\n"); + return false; + } + + pmFPAview *thisView = pmFPAviewAlloc(view->nRows); // Copy of input view + *thisView = *view; + + for (int i = 0; i < chip->cells->n; i++) { + pmCell *cell = chip->cells->data[i]; + thisView->cell = i; + if (!pmCellWritePattern(cell, thisView, file, config)) { + psError(PS_ERR_IO, false, "Failed to write %dth cell", i); + psFree(thisView); + return false; + } + } + psFree(thisView); + return true; +} + +static bool pmFPAWritePattern(pmFPA *fpa, const pmFPAview *view, + pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(fpa, false); + PS_ASSERT_PTR_NON_NULL(fpa->chips, false); + + if (!pmConceptsWriteFPA(fpa, true, config)) { + psError(PS_ERR_IO, false, "Unable to write concepts for FPA.\n"); + return false; + } + + pmFPAview *thisView = pmFPAviewAlloc(view->nRows); // Copy of input view + *thisView = *view; + + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + thisView->chip = i; + if (!pmChipWritePattern(chip, thisView, file, config)) { + psError(PS_ERR_IO, false, "Failed to write %dth chip", i); + psFree(thisView); + return false; + } + } + psFree(thisView); + return true; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + +bool pmReadoutReadPattern(pmReadout *ro, psFits *fits) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PS_ASSERT_FITS_NON_NULL(fits, false); + + pmCell *cell = ro->parent; // Cell of interest + if (!cell) { + psError(PS_ERR_UNEXPECTED_NULL, true, "No cell associated with readout."); + return false; + } + pmChip *chip = cell->parent; // Chip of interest + pmFPA *fpa = chip->parent; // FPA of interest + pmHDU *hdu = pmHDUGetLowest(fpa, chip, cell); // HDU for readout + if (!hdu) { + psError(PS_ERR_UNEXPECTED_NULL, true, "No HDU associated with readout."); + return false; + } + + if (!psFitsMoveExtName(fits, hdu->extname)) { + psError(PS_ERR_IO, false, "Unable to move to pattern correction."); + return false; + } + + psMetadata *header = psFitsReadHeader(NULL, fits); // Header + if (!header) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to read header for pattern correction."); + return false; + } + + bool data = true; // Did we find any data? + + psMetadataItem *cellCorr = psMetadataLookup(header, PM_PATTERN_CELL_CORRECTION); // Cell pattern correction + if (cellCorr) { + psMetadataAddItem(ro->analysis, cellCorr, PS_LIST_TAIL, PS_META_REPLACE); + data = true; + } + + int naxis = psMetadataLookupS32(NULL, header, "NAXIS"); // Number of axes + if (naxis > 0) { + psImage *rowCorr = psFitsReadImage(fits, psRegionSet(0, 0, 0, 0), 0); // Row pattern correction + psMetadataAddImage(ro->analysis, PS_LIST_TAIL, PM_PATTERN_ROW_CORRECTION, PS_META_REPLACE, + "Pattern row correction", rowCorr); + psFree(rowCorr); + data = true; + } + + if (data) { + ro->data_exists = true; + ro->parent->data_exists = true; + ro->parent->parent->data_exists = true; + } + + psFree(header); + + return true; +} + +static bool pmCellReadPattern(pmCell *cell, const pmFPAview *view, + pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(cell, false); + PS_ASSERT_PTR_NON_NULL(cell->readouts, false); + + pmFPAview *thisView = pmFPAviewAlloc(view->nRows); // Copy of input view + *thisView = *view; + + // Create a readout if none exists + if (!cell->readouts || cell->readouts->n == 0) { + pmReadout *readout = pmReadoutAlloc(cell); // New readout + psFree(readout); // Drop reference + } + + cell->data_exists = false; + for (int i = 0; i < cell->readouts->n; i++) { + pmReadout *readout = cell->readouts->data[i]; + thisView->readout = i; + if (!pmReadoutReadPattern(readout, file->fits)) { + psError(PS_ERR_IO, false, "Unable to read pattern correction."); + return false; + } + } + psFree(thisView); + + if (!pmCellReadHeader(cell, file->fits, config)) { + psError(PS_ERR_IO, false, "Unable to read header for cell."); + return false; + } + // load in the concept information for this cell + if (!pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_HEADER, true, NULL)) { + psErrorClear(); + psWarning("Difficulty reading concepts for cell; attempting to proceed."); + } + + return true; +} + +static bool pmChipReadPattern(pmChip *chip, const pmFPAview *view, + pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(chip, false); + PS_ASSERT_PTR_NON_NULL(chip->cells, false); + + pmFPAview *thisView = pmFPAviewAlloc(view->nRows); // Copy of input view + *thisView = *view; + + chip->data_exists = false; + for (int i = 0; i < chip->cells->n; i++) { + pmCell *cell = chip->cells->data[i]; + thisView->cell = i; + pmCellReadPattern(cell, thisView, file, config); + if (!cell->data_exists) { + continue; + } + chip->data_exists = true; + } + psFree(thisView); + + if (!pmChipReadHeader(chip, file->fits, config)) { + psError(PS_ERR_IO, false, "Unable to read header for cell."); + return false; + } + if (!pmConceptsReadChip(chip, PM_CONCEPT_SOURCE_HEADER, true, true, NULL)) { + psError(PS_ERR_IO, false, "Failed to read concepts for chip.\n"); + return false; + } + + return true; +} + +static bool pmFPAReadPattern(pmFPA *fpa, const pmFPAview *view, + pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + PS_ASSERT_PTR_NON_NULL(file->fpa, false); + PS_ASSERT_PTR_NON_NULL(file->fpa->chips, false); + + pmFPAview *thisView = pmFPAviewAlloc(view->nRows); // Copy of input view + *thisView = *view; + + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + thisView->chip = i; + pmChipReadPattern(chip, thisView, file, config); + } + psFree(thisView); + + if (!pmFPAReadHeader(fpa, file->fits, config)) { + psError(PS_ERR_IO, false, "Unable to read header for cell."); + return false; + } + if (!pmConceptsReadFPA(fpa, PM_CONCEPT_SOURCE_HEADER, true, NULL)) { + psError(PS_ERR_IO, false, "Failed to read concepts for fpa.\n"); + return false; + } + + return true; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +bool pmPatternWrite(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + PS_ASSERT_PTR_NON_NULL(file->fpa, false); + + pmFPA *fpa = pmFPAfileSuitableFPA(file, view, config, false); // Suitable FPA for writing + + if (view->chip == -1) { + if (!pmFPAWritePattern(fpa, view, file, config)) { + psError(PS_ERR_IO, false, "Failed to write pattern correction from fpa"); + psFree(fpa); + return false; + } + psFree(fpa); + return true; + } + + if (view->chip >= fpa->chips->n) { + psError(PS_ERR_UNKNOWN, false, "Writing chip == %d (>= chips->n == %ld)", view->chip, fpa->chips->n); + psFree(fpa); + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + if (!pmChipWritePattern(chip, view, file, config)) { + psError(PS_ERR_IO, false, "Failed to write pattern correction from chip"); + psFree(fpa); + return false; + } + psFree(fpa); + return true; + } + + if (view->cell >= chip->cells->n) { + psError(PS_ERR_UNKNOWN, false, "Writing cell == %d (>= cells->n == %ld)", + view->cell, chip->cells->n); + psFree(fpa); + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + if (!pmCellWritePattern(cell, view, file, config)) { + psError(PS_ERR_IO, false, "Failed to write pattern correction from cell"); + psFree(fpa); + return false; + } + psFree(fpa); + return true; + } + + if (view->readout >= cell->readouts->n) { + psError(PS_ERR_UNKNOWN, false, "Writing readout == %d (>= readouts->n == %ld)", + view->readout, cell->readouts->n); + psFree(fpa); + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + if (!pmReadoutWritePattern(readout, file->fits)) { + psError(PS_ERR_IO, false, "Failed to write pattern correction from readout %d", view->readout); + psFree(fpa); + return false; + } + + psFree(fpa); + return true; +} + +bool pmPatternWritePHU(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + + if (file->wrote_phu) { + return true; + } + + // find the FPA phu + pmFPA *fpa = pmFPAfileSuitableFPA(file, view, config, false); // Suitable FPA for writing + pmHDU *phu = psMemIncrRefCounter(pmFPAviewThisPHU(view, fpa)); + psFree(fpa); + + // if there is no PHU, this is a single header+image (extension-less) file. This could be the case for an + // input SPLIT set of files being written out as a MEF. if there is a PHU, write it out as a 'blank' + psMetadata *outhead = psMetadataAlloc(); + if (phu) { + psMetadataCopy (outhead, phu->header); + } + psFree(phu); + + pmConfigConformHeader(outhead, file->format); + + psFitsWriteBlank(file->fits, outhead, ""); + file->wrote_phu = true; + + psTrace("pmFPAfile", 5, "wrote phu %s (type: %d)\n", file->filename, file->type); + psFree(outhead); + + return true; +} + +bool pmPatternRead(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + PS_ASSERT_PTR_NON_NULL(file->fpa, false); + + pmFPA *fpa = file->fpa; + + if (view->chip == -1) { + pmFPAReadPattern(fpa, view, file, config); + return true; + } + + if (view->chip >= fpa->chips->n) { + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + pmChipReadPattern(chip, view, file, config); + return true; + } + + if (view->cell >= chip->cells->n) { + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + pmCellReadPattern(cell, view, file, config); + return true; + } + + if (view->readout >= cell->readouts->n) { + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + return pmReadoutReadPattern(readout, file->fits); +} Index: trunk/psModules/src/detrend/pmPatternIO.h =================================================================== --- trunk/psModules/src/detrend/pmPatternIO.h (revision 26893) +++ trunk/psModules/src/detrend/pmPatternIO.h (revision 26893) @@ -0,0 +1,37 @@ +#ifndef PM_PATTERN_IO_H +#define PM_PATTERN_IO_H + +#include + +#include +#include + +/// Write pattern correction within a readout to a FITS file +bool pmReadoutWritePattern( + pmReadout *readout, ///< Readout for which to write pattern correction (in analysis MD) + psFits *fits ///< FITS file to which to write + ); + + +bool pmReadoutReadPattern( + pmReadout *readout, ///< Readout for which to read pattern correction (into analysis MD) + psFits *fits ///< FITS file to which to write + ); + + +bool pmPatternRead(const pmFPAview *view, ///< View into which to read + pmFPAfile *file, ///< File from which to read + pmConfig *config ///< Configuration + ); + +bool pmPatternWrite(const pmFPAview *view, ///< View from which to write + pmFPAfile *file, ///< File to which to write + pmConfig *config ///< Configuration + ); + +bool pmPatternWritePHU(const pmFPAview *view, // View to PHU + pmFPAfile *file, ///< File to which to write + pmConfig *config ///< Configuration + ); + +#endif