Index: trunk/psModules/src/camera/pmFPAWrite.c =================================================================== --- trunk/psModules/src/camera/pmFPAWrite.c (revision 21279) +++ trunk/psModules/src/camera/pmFPAWrite.c (revision 21363) @@ -9,4 +9,5 @@ #include "pmConfig.h" +#include "pmConfigMask.h" #include "pmHDU.h" #include "pmFPA.h" @@ -29,5 +30,5 @@ FPA_WRITE_TYPE_IMAGE, // Write image FPA_WRITE_TYPE_MASK, // Write mask - FPA_WRITE_TYPE_WEIGHT // Write weight map + FPA_WRITE_TYPE_VARIANCE // Write variance map } fpaWriteType; @@ -46,6 +47,6 @@ case FPA_WRITE_TYPE_MASK: return &hdu->masks; - case FPA_WRITE_TYPE_WEIGHT: - return &hdu->weights; + case FPA_WRITE_TYPE_VARIANCE: + return &hdu->variances; default: psAbort("Unknown write type: %x\n", type); @@ -66,6 +67,6 @@ case FPA_WRITE_TYPE_MASK: return pmHDUWriteMask(hdu, fits, config); - case FPA_WRITE_TYPE_WEIGHT: - return pmHDUWriteWeight(hdu, fits, config); + case FPA_WRITE_TYPE_VARIANCE: + return pmHDUWriteVariance(hdu, fits, config); default: psAbort("Unknown write type: %x\n", type); @@ -74,5 +75,70 @@ } -// Write a cell image/mask/weight +// Indicate whether a covariance matrix is defined +static bool readoutSearchCovariances(pmReadout *ro) +{ + return ro->covariance ? true : false; +} + +// Search for a covariance matrix +#define SEARCH_COVARIANCES(NAME, PARENT, CHILD, CHILDREN, TESTFUNC) \ +static bool NAME(PARENT *parent) \ +{ \ + if (!parent || !parent->CHILDREN) { \ + return false; \ + } \ + psArray *children = parent->CHILDREN; /* Array of children */ \ + for (int i = 0; i < children->n; i++) { \ + CHILD *child = children->data[i]; /* Child of interest */ \ + if (child && TESTFUNC(child)) { \ + return true; \ + } \ + } \ + return false; \ +} + +SEARCH_COVARIANCES(cellSearchCovariances, pmCell, pmReadout, readouts, readoutSearchCovariances); +SEARCH_COVARIANCES(chipSearchCovariances, pmChip, pmCell, cells, cellSearchCovariances); +SEARCH_COVARIANCES(fpaSearchCovariances, pmFPA, pmChip, chips, chipSearchCovariances); + +// Some type-specific additions to the header +static bool writeUpdateHeader(pmFPA *fpa, // FPA of interest + pmChip *chip, // Chip of interest, or NULL + pmCell *cell, // Cell of interest, or NULL + fpaWriteType type, // Type to write + pmConfig *config // Configuration + ) +{ + switch (type) { + case FPA_WRITE_TYPE_MASK: { + pmHDU *phu = pmHDUGetHighest(fpa, chip, cell); // Primary header + if (!pmConfigMaskWriteHeader(config, phu->header)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set the mask names in the PHU header"); + return false; + } + break; + } + case FPA_WRITE_TYPE_VARIANCE: { + bool covar = false; // Are covariances present? + if ((cell && cellSearchCovariances(cell)) || + (!cell && ((chip && chipSearchCovariances(chip)) || + (!chip && fpa && fpaSearchCovariances(fpa))))) { + covar = true; + } + + pmHDU *hdu = pmHDUGetLowest(fpa, chip, cell); // Header being written + psMetadataAddBool(hdu->header, PS_LIST_TAIL, PM_HDU_COVARIANCE_KEYWORD, PS_META_REPLACE, + "Is a covariance matrix present?", covar); + break; + } + default: + break; + } + + return true; +} + + +// Write a cell image/mask/variance static bool cellWrite(pmCell *cell, // Cell to write psFits *fits, // FITS file to which to write @@ -97,5 +163,5 @@ // generate the HDU, but only copies the structure. if (!blank && !hdu->blankPHU && !*imageArray && (!pmHDUGenerateForCell(cell) || !*imageArray)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate HDU for cell --- likely programming error.\n"); + psError(PS_ERR_UNKNOWN, false, "Unable to generate HDU for cell --- likely programming error."); return false; } @@ -106,9 +172,12 @@ if (writeBlank || writeImage) { - pmConceptSource source = PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CELLS | PM_CONCEPT_SOURCE_DEFAULTS | PM_CONCEPT_SOURCE_DATABASE; if (!pmConceptsWriteCell(cell, source, true, config)) { - psError(PS_ERR_IO, false, "Unable to write concepts for cell.\n"); + psError(PS_ERR_IO, false, "Unable to write concepts for cell."); + return false; + } + if (!writeUpdateHeader(NULL, NULL, cell, type, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to update header for writing"); return false; } @@ -123,5 +192,5 @@ } -// Write a chip image/mask/weight +// Write a chip image/mask/variance static bool chipWrite(pmChip *chip, // Chip to write psFits *fits, // FITS file to which to write @@ -162,4 +231,10 @@ return false; } + + if (!writeUpdateHeader(NULL, chip, NULL, type, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to update header for writing"); + return false; + } + if (!appropriateWriteFunc(hdu, fits, config, type)) { psError(PS_ERR_IO, false, "Unable to write HDU for chip.\n"); @@ -186,5 +261,5 @@ -// Write an FPA image/mask/weight +// Write an FPA image/mask/variance static bool fpaWrite(pmFPA *fpa, // FPA to write psFits *fits, // FITS file to which to write @@ -223,4 +298,8 @@ if (!pmConceptsWriteFPA(fpa, source, true, config)) { psError(PS_ERR_IO, false, "Unable to write concepts for FPA.\n"); + return false; + } + if (!writeUpdateHeader(fpa, NULL, NULL, type, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to update header for writing"); return false; } @@ -433,23 +512,41 @@ -bool pmCellWriteWeight(pmCell *cell, psFits *fits, pmConfig *config, bool blank) +bool pmCellWriteVariance(pmCell *cell, psFits *fits, pmConfig *config, bool blank) { PS_ASSERT_PTR_NON_NULL(cell, false); PS_ASSERT_PTR_NON_NULL(fits, false); - return cellWrite(cell, fits, config, blank, FPA_WRITE_TYPE_WEIGHT); -} - -bool pmChipWriteWeight(pmChip *chip, psFits *fits, pmConfig *config, bool blank, bool recurse) + if (!cellWrite(cell, fits, config, blank, FPA_WRITE_TYPE_VARIANCE)) { + return false; + } + if (!pmCellWriteCovariance(fits, cell)) { + return false; + } + return true; +} + +bool pmChipWriteVariance(pmChip *chip, psFits *fits, pmConfig *config, bool blank, bool recurse) { PS_ASSERT_PTR_NON_NULL(chip, false); PS_ASSERT_PTR_NON_NULL(fits, false); - return chipWrite(chip, fits, config, blank, recurse, FPA_WRITE_TYPE_WEIGHT); -} - -bool pmFPAWriteWeight(pmFPA *fpa, psFits *fits, pmConfig *config, bool blank, bool recurse) + if (!chipWrite(chip, fits, config, blank, recurse, FPA_WRITE_TYPE_VARIANCE)) { + return false; + } + if (!pmChipWriteCovariance(fits, chip)) { + return false; + } + return true; +} + +bool pmFPAWriteVariance(pmFPA *fpa, psFits *fits, pmConfig *config, bool blank, bool recurse) { PS_ASSERT_PTR_NON_NULL(fpa, false); PS_ASSERT_PTR_NON_NULL(fits, false); - return fpaWrite(fpa, fits, config, blank, recurse, FPA_WRITE_TYPE_WEIGHT); + if (!fpaWrite(fpa, fits, config, blank, recurse, FPA_WRITE_TYPE_VARIANCE)) { + return false; + } + if (!pmFPAWriteCovariance(fits, fpa)) { + return false; + } + return true; } @@ -526,2 +623,120 @@ return numWrite; } + +bool pmCellWriteCovariance(psFits *fits, const pmCell *cell) +{ + PS_ASSERT_PTR_NON_NULL(cell, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + + int numCovar = 0; + psArray *readouts = cell->readouts; // Array of readouts + for (int i = 0; i < readouts->n; i++) { + pmReadout *readout = readouts->data[i]; // The readout of interest + if (readout && readout->covariance) { + numCovar++; + } + } + if (numCovar == 0) { + return true; + } + if (numCovar != readouts->n) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Number of covariances (%d) doesn't match number of readouts (%ld)", + numCovar, readouts->n); + return false; + } + + // Check size of covariances + int xMinCovar = INT_MAX, xMaxCovar = INT_MIN, yMinCovar = INT_MAX, yMaxCovar = INT_MIN; // Size + for (int i = 0; i < readouts->n; i++) { + pmReadout *readout = readouts->data[i]; // The readout of interest + psAssert(readout, "Should be defined."); + psKernel *covar = readout->covariance; // Covariance matrix + psAssert(covar, "Should be defined."); + xMinCovar = PS_MIN(xMinCovar, covar->xMin); + xMaxCovar = PS_MAX(xMaxCovar, covar->xMax); + yMinCovar = PS_MIN(yMinCovar, covar->yMin); + yMaxCovar = PS_MAX(yMaxCovar, covar->yMax); + } + + // Correct covariances to common size + psArray *images = psArrayAlloc(numCovar); // Array of images + for (int i = 0; i < readouts->n; i++) { + pmReadout *readout = readouts->data[i]; // The readout of interest + psAssert(readout, "Should be defined."); + psKernel *covar = readout->covariance; // Covariance matrix + psAssert(covar, "Should be defined."); + int xMin = covar->xMin, xMax = covar->xMax, yMin = covar->yMin, yMax = covar->yMax;// Size + if (xMin == xMinCovar && xMax == xMaxCovar && yMin == yMinCovar && yMax == yMaxCovar) { + images->data[i] = psMemIncrRefCounter(covar->image); + } else { + psImage *new = psImageAlloc(xMaxCovar - xMinCovar + 1, yMaxCovar - yMinCovar + 1, PS_TYPE_F32); + psImageInit(new, 0); + psImageOverlaySection(new, covar->image, xMinCovar - xMin, yMinCovar - yMin, "="); + images->data[i] = new; + } + } + + // Determine extension name + const char *chipName = psMetadataLookupStr(NULL, cell->parent->concepts, "CHIP.NAME"); // Name of chip + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell + psString extname = NULL; // Extension name + psStringAppend(&extname, "COVAR_%s_%s", chipName, cellName); + + // Generate header + pmHDU *hdu = pmHDUFromCell(cell); // HDU for cell + psMetadata *header = psMetadataCopy(NULL, hdu->header); // Header to write + psMetadataAddS32(header, PS_LIST_TAIL, "COVARIANCE.CENTRE.X", PS_META_REPLACE, + "Centre of covariance matrix in x", -xMinCovar); + psMetadataAddS32(header, PS_LIST_TAIL, "COVARIANCE.CENTRE.Y", PS_META_REPLACE, + "Centre of covariance matrix in y", -yMinCovar); + + // Write images + if (!psFitsWriteImageCube(fits, header, images, extname)) { + psError(PS_ERR_UNKNOWN, false, "Unable to write covariances from chip %s, cell %s to extension %s", + chipName, cellName, extname); + psFree(extname); + psFree(header); + psFree(images); + return 0; + } + psFree(extname); + psFree(header); + psFree(images); + + return true; +} + + +bool pmChipWriteCovariance(psFits *fits, const pmChip *chip) +{ + PS_ASSERT_PTR_NON_NULL(chip, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + + psArray *cells = chip->cells; // Array of cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // Cell of interest + if (!pmCellWriteCovariance(fits, cell)) { + return false; + } + } + + return true; +} + + +bool pmFPAWriteCovariance(psFits *fits, const pmFPA *fpa) +{ + PS_ASSERT_PTR_NON_NULL(fpa, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + + psArray *chips = fpa->chips; // Array of chips + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // Chip of interest + if (!pmChipWriteCovariance(fits, chip)) { + return false; + } + } + + return true; +}