Index: /branches/pap_branch_20090128/psModules/src/camera/pmFPARead.c =================================================================== --- /branches/pap_branch_20090128/psModules/src/camera/pmFPARead.c (revision 21263) +++ /branches/pap_branch_20090128/psModules/src/camera/pmFPARead.c (revision 21264) @@ -28,5 +28,5 @@ FPA_READ_TYPE_IMAGE, // Read image FPA_READ_TYPE_MASK, // Read mask - FPA_READ_TYPE_VARIANCE, // Read variance map + FPA_READ_TYPE_VARIANCE, // Read variance map FPA_READ_TYPE_HEADER // Read header } fpaReadType; @@ -515,11 +515,22 @@ } - // XXX for IMAGE, we need the CELL.BAD value, but for MASK, we need the BAD mask value - - float bad = 0; - if (type == FPA_READ_TYPE_MASK) { - bad = 1.0; - } else { - bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); // Bad level + // Need to set the invalid (unread) pixels appropriately, and to read the mask bits + float bad = 0; // Bad level + switch (type) { + case FPA_READ_TYPE_MASK: { + pmHDU *phu = pmHDUGetHighest(cell->parent->parent, cell->parent, cell); // Primary header + if (!pmConfigMaskReadHeader(config, phu->header)) { + psError(PS_ERR_IO, false, "Unable to determine mask bits"); + return false; + } + bad = pmConfigMaskGet("BAD", config); + break; + } + case FPA_READ_TYPE_IMAGE: + case FPA_READ_TYPE_VARIANCE: + bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); + break; + default: + psAbort("Unrecognised type: %x", type); } @@ -721,4 +732,10 @@ return NULL; } + + if (type == FPA_READ_TYPE_VARIANCE && hdu->covariances) { + psArray *covariances = hdu->covariances; // Covariances in HDU + for (int i = 0; i < covariances->n; i++) { + pmHDUCovariances *covar = covariances->data[i]; // Covariance information + psFree(readout); // Drop reference } @@ -852,5 +869,5 @@ float bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); // Bad level if (!mdok) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.BAD is not set --- assuming zero.\n"); + psWarning("CELL.BAD is not set --- assuming zero.\n"); bad = 0.0; } Index: /branches/pap_branch_20090128/psModules/src/camera/pmFPAWrite.c =================================================================== --- /branches/pap_branch_20090128/psModules/src/camera/pmFPAWrite.c (revision 21263) +++ /branches/pap_branch_20090128/psModules/src/camera/pmFPAWrite.c (revision 21264) @@ -29,5 +29,5 @@ FPA_WRITE_TYPE_IMAGE, // Write image FPA_WRITE_TYPE_MASK, // Write mask - FPA_WRITE_TYPE_VARIANCE // Write variance map + FPA_WRITE_TYPE_VARIANCE // Write variance map } fpaWriteType; @@ -74,4 +74,28 @@ } + + + +// 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 + ) +{ + 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: { + pmHDU *hdu = pmHDUGetLowest(fpa, chip, cell); // Header being written + + + // Write a cell image/mask/variance static bool cellWrite(pmCell *cell, // Cell to write @@ -97,5 +121,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,12 +130,13 @@ if (writeBlank || writeImage) { - pmFPAUpdateNames(cell->parent->parent, cell->parent, cell); 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 (!appropriateWriteFunc(hdu, fits, config, type)) { psError(PS_ERR_IO, false, "Unable to write HDU for cell.\n"); Index: /branches/pap_branch_20090128/psModules/src/camera/pmFPAfileIO.c =================================================================== --- /branches/pap_branch_20090128/psModules/src/camera/pmFPAfileIO.c (revision 21263) +++ /branches/pap_branch_20090128/psModules/src/camera/pmFPAfileIO.c (revision 21264) @@ -771,13 +771,4 @@ } - // XXX if we have a mask file, then we need to read the mask bit names - // defined for this file - if (file->type == PM_FPA_FILE_MASK) { - if (!pmConfigMaskReadHeader (config, phu)) { - psError(PS_ERR_IO, false, "error in mask bits"); - return false; - } - } - // determine the current format from the header // determine camera if not specified already @@ -786,5 +777,5 @@ psString formatName = NULL; psString cameraName = NULL; - file->format = pmConfigCameraFormatFromHeader (&camera, &cameraName, &formatName, config, phu, true); + file->format = pmConfigCameraFormatFromHeader(&camera, &cameraName, &formatName, config, phu, true); if (!file->format) { psError(PS_ERR_IO, false, "Failed to read CCD format from %s\n", file->filename); Index: /branches/pap_branch_20090128/psModules/src/camera/pmHDU.c =================================================================== --- /branches/pap_branch_20090128/psModules/src/camera/pmHDU.c (revision 21263) +++ /branches/pap_branch_20090128/psModules/src/camera/pmHDU.c (revision 21264) @@ -12,4 +12,6 @@ #include "pmHDU.h" #include "pmFPA.h" + +#define COVARIANCE_INDICATOR "PS_COVAR" // FITS keyword to indicate if covariance matrices are present ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -160,5 +162,9 @@ PS_ASSERT_PTR_NON_NULL(fits, false); - return hduRead(hdu, &hdu->variances, fits); + if (!hduRead(hdu, &hdu->variances, fits)) { + return false; + } + + return true; } @@ -240,5 +246,17 @@ PS_ASSERT_PTR_NON_NULL(fits, false); + pmHDUCovariancesWriteHeader(hdu); + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Value to mask - return hduWrite(hdu, hdu->variances, hdu->masks, maskVal, fits); -} + if (!hduWrite(hdu, hdu->variances, hdu->masks, maskVal, fits)) { + return false; + } + pmHDUCovarianceClearHeader(hdu); + + if (!pmHDUCovarianceWrite(hdu, fits)) { + psError(PS_ERR_UNKNOWN, false, "Unable to write covariance"); + return false; + } + + return true; +} Index: /branches/pap_branch_20090128/psModules/src/camera/pmHDU.h =================================================================== --- /branches/pap_branch_20090128/psModules/src/camera/pmHDU.h (revision 21263) +++ /branches/pap_branch_20090128/psModules/src/camera/pmHDU.h (revision 21264) @@ -4,6 +4,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.9.22.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-29 00:33:51 $ + * @version $Revision: 1.9.22.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-02 18:05:38 $ * Copyright 2005-2006 Institute for Astronomy, University of Hawaii */ @@ -22,15 +22,14 @@ /// Of course, it is not an exact replica of a FITS HDU --- they have no mask and variance data, but these are /// stored here for convenience --- it keeps all the relevant data about the image in one place. -typedef struct -{ +typedef struct { psString extname; ///< The extension name bool blankPHU; ///< Is this a blank FITS Primary Header Unit, i.e., no data? psMetadata *format; ///< The camera format psMetadata *header; ///< The FITS header, or NULL if primary for FITS; or section info - psArray *images; ///< The pixel data - psArray *variances; ///< The pixel data - psArray *masks; ///< The pixel data -} -pmHDU; + psArray *images; ///< Pixel data + psArray *variances; ///< Variance in the pixel data, or NULL + psArray *masks; ///< Mask for the pixel data, or NULL + psArray *covariances; ///< Covariance matrices (pmHDUCovariance), or NULL +} pmHDU; Index: /branches/pap_branch_20090128/psModules/src/camera/pmHDUCovariances.c =================================================================== --- /branches/pap_branch_20090128/psModules/src/camera/pmHDUCovariances.c (revision 21264) +++ /branches/pap_branch_20090128/psModules/src/camera/pmHDUCovariances.c (revision 21264) @@ -0,0 +1,167 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "pmHDU.h" +#include "pmFPA.h" + +#include "pmHDUCovariances.h" + +/// Storage for covariance matrices of a cell +/// +/// We want to store the covariance matrices applicable for a single HDU together. However, when we write an +/// HDU, we know nothing about chips and cells, so we require that it be populated for us. This structure is +/// more convenient than carrying embedded psMetadata. +typedef struct { + psString chipName; ///< Name of chip + psString cellName; ///< Name of cell + int xMin, xMax, yMin, yMax; ///< Size of covariance matrices + psArray *covariances; ///< Covariance matrices for cell (one per readout) +} pmHDUCovariance; + +// Deallocator +static void hduCovarianceFree(pmHDUCovariance *covar) +{ + psFree(covar->chipName); + psFree(covar->cellName); + psFree(covar->covariances); + return; +} + +/// Allocator +pmHDUCovariance *pmHDUCovarianceAlloc(void) +{ + pmHDUCovariance *covar = psAlloc(sizeof(pmHDUCovariance)); // Covariance information to return + psMemSetDeallocator(covar, (psFreeFunc)hduCovarianceFree); + + covar->chipName = NULL; + covar->cellName = NULL; + covar->xMin = 0; + covar->xMax = 0; + covar->yMin = 0; + covar->yMax = 0; + covar->covariances = NULL; + + return covar; +} + +bool pmHDUCovarianceWriteHeader(pmHDU *hdu) +{ + PM_ASSERT_HDU_NON_NULL(hdu, false); + + return psMetadataAddBool(hdu->header, PS_LIST_TAIL, COVARIANCE_INDICATOR, PS_META_REPLACE, + "Covariance matrices present?", hdu->covariances ? true : false); +} + +bool pmHDUCovarianceClearHeader(pmHDU *hdu) +{ + PM_ASSERT_HDU_NON_NULL(hdu, false); + + return psMetadataRemoveKey(hdu->header, COVARIANCE_INDICATOR); +} + +bool pmHDUCovarianceWrite(pmHDU *hdu, psFits *fits) +{ + PM_ASSERT_HDU_NON_NULL(hdu, false); + PS_ASSERT_FITS_NON_NULL(fits, false); + + psArray *covariances = hdu->covariances; // Covariance matrices + if (!covariances) { + return true; + } + + psArray *table = psArrayAlloc(covariances->n); // FITS table to write + for (int i = 0; i < covariances->n; i++) { + pmHDUCovariance *covar = covariances->data[i]; // Covariance information + psMetadata *row = psMetadataAlloc(); // Table row + psMetadataAddStr(row, PS_LIST_TAIL, "CHIP", 0, "Chip name", covar->chipName); + psMetadataAddStr(row, PS_LIST_TAIL, "CELL", 0, "Cell name", covar->cellName); + psMetadataAddS32(row, PS_LIST_TAIL, "READOUT", 0, "Readout number", covar->readoutNum); + psMetadataAddS32(row, PS_LIST_TAIL, "XMIN", 0, "Minimum x value", covar->xMin); + psMetadataAddS32(row, PS_LIST_TAIL, "XMAX", 0, "Maximum x value", covar->xMax); + psMetadataAddS32(row, PS_LIST_TAIL, "YMIN", 0, "Minimum y value", covar->yMin); + psMetadataAddS32(row, PS_LIST_TAIL, "YMAX", 0, "Maximum x value", covar->yMax); + int numCols = covar->xMax - covar->xMin + 1, numRows = covar->yMax - covar->yMin + 1; // Size + psVector *data = psVectorAlloc(numCols * numRows, PS_TYPE_F32); // Vector with data + for (int i = 0, j = 0; i < numRows; i++, j += numCols) { + memcpy(&data->data.F32[j], kernel->image[i], numCols * PS_ELEMTYPE_SIZEOF(PS_TYPE_F32)); + } + psMetadataAddVector(row, PS_LIST_TAIL, "COVAR", 0, "Covariance matrix, row-major order", data); + psFree(data); // Drop reference + table->data[i] = row; + } + psString *extname = NULL; // Extension name for covariance data + psStringAppend(&extname, "%s_COVAR", hdu->extname); + if (!psFitsWriteTable(fits, hdu->header, table, extname)) { + psError(PS_ERR_IO, false, "Unable to write covariance table."); + psFree(extname); + psFree(table); + return false; + } + psFree(extname); + psFree(table); + + return true; +} + + +/// Read covariance into an HDU +bool pmHDUCovarianceRead(pmHDU *hdu) +{ + PM_ASSERT_HDU_NON_NULL(hdu, false); + PS_ASSERT_METADATA_NON_NULL(hdu->header, false); + + bool mdok; // Status of MD lookup + if (!psMetadataLookupBool(&mdok, hdu->header, COVARIANCE_INDICATOR)) { + return true; + } + + psString extname = NULL; // Name of extension + psStringAppend(&extname, "%s_COVAR", hdu->extname); + if (!psFitsMoveExtName(fits, extname)) { + psError(PS_ERR_IO, false, "Unable to move to covariance extension"); + psFree(extname); + return false; + } + psFree(extname); + + psArray *table = psFitsReadTable(fits); // Table with covariances + if (!table) { + psError(PS_ERR_IO, false, "Unable to read covariances table"); + return false; + } + psArray *covariances = hdu->covariances = psArrayAlloc(table->n); // Covariance matrices + for (int i = 0; i < table->n; i++) { + psMetadata *row = table->data[i]; // Table row + pmHDUCovariance *covar = pmHDUCovarianceAlloc(); + covar->chipName = psMetadataLookupStr(NULL, row, "CHIP"); + covar->cellName = psMetadataLookupStr(NULL, row, "CELL"); + covar->xMin = psMetadataLookupStr(NULL, row, "XMIN"); + covar->xMax = psMetadataLookupStr(NULL, row, "XMAX"); + covar->yMin = psMetadataLookupStr(NULL, row, "YMIN"); + covar->yMax = psMetadataLookupStr(NULL, row, "YMAX"); + + covariances->data[i] = covar; + } +} + + +bool pmHDUCovarianceAdd(pmHDU *hdu, pmReadout *readout) +{ + PM_ASSERT_HDU_NON_NULL(hdu, false); + PM_ASSERT_READOUT_NON_NULL(readout, false); + PM_ASSERT_READOUT_PARENT_NON_NULL(readout, false); + +} + +/// Select the appropriate covariance from an HDU +bool pmHDUCovarianceSelect( + pmReadout *readout, ///< Readout for which to select covariance + pmHDU *hdu, ///< Header data unit with covariance + ); + + +#endif Index: /branches/pap_branch_20090128/psModules/src/camera/pmHDUCovariances.h =================================================================== --- /branches/pap_branch_20090128/psModules/src/camera/pmHDUCovariances.h (revision 21264) +++ /branches/pap_branch_20090128/psModules/src/camera/pmHDUCovariances.h (revision 21264) @@ -0,0 +1,41 @@ +#ifndef PM_HDU_COVARIANCES_H +#define PM_HDU_COVARIANCES_H + +#include +#include "pmHDU.h" +#include "pmFPA.h" + +/// Storage for covariance matrices of a cell +/// +/// We want to store the covariance matrices applicable for a single HDU together. However, when we write an +/// HDU, we know nothing about chips and cells, so we require that it be populated for us. This structure is +/// more convenient than carrying embedded psMetadata. +typedef struct { + psString chipName; ///< Name of chip + psString cellName; ///< Name of cell + int readoutNum; ///< Number of readout + int xMin, xMax, yMin, yMax; ///< Size of covariance matrices + psArray *covariances; ///< Covariance matrices for cell (one per readout) +} pmHDUCovariance; + +/// Allocator +pmHDUCovariance *pmHDUCovarianceAlloc(void); + +/// Read covariance into an HDU +bool pmHDUCovarianceRead( + pmHDU *hdu ///< Header data unit into which to read + ); + +/// Write covariance from an HDU +bool pmHDUCovarianceWrite( + pmHDU *hdu ///< Header data unit from which to write + ); + +/// Select the appropriate covariance from an HDU +bool pmHDUCovarianceSelect( + pmReadout *readout, ///< Readout for which to select covariance + pmHDU *hdu, ///< Header data unit with covariance + ); + + +#endif Index: /branches/pap_branch_20090128/psModules/src/camera/pmHDUGenerate.c =================================================================== --- /branches/pap_branch_20090128/psModules/src/camera/pmHDUGenerate.c (revision 21263) +++ /branches/pap_branch_20090128/psModules/src/camera/pmHDUGenerate.c (revision 21264) @@ -355,5 +355,6 @@ psElemType imageType = 0; // Type of readout images psElemType maskType = 0; // Type of readout masks - psElemType varianceType = 0; // Type of readout variances + psElemType varianceType = 0; // Type of readout variances + psElemType covarianceType = 0; // Type of readout covariances { psListIterator *iter = psListIteratorAlloc(cells, PS_LIST_HEAD, false); // Iterator for cells @@ -384,9 +385,12 @@ varianceType = checkTypes(varianceType, readout->variance->type.type); } + if (!hdu->covariances && readout->covariance) { + covarianceType = checkTypes(covarianceType, PS_TYPE_F32); + } } } psFree(iter); } - if (numReadouts == 0 || (imageType == 0 && maskType == 0 && varianceType == 0)) { + if (numReadouts == 0 || (imageType == 0 && maskType == 0 && varianceType == 0 && covarianceType == 0)) { // Nothing from which to create an HDU psFree(cells); @@ -428,4 +432,7 @@ hdu->variances->data[i] = variance; } + } + if (covarianceType) { + hdu->covariances = psHashAlloc(psListLength(cells)); } @@ -448,7 +455,25 @@ psArray *readouts = cell->readouts; // Array of readouts + + // Check size of covariances + int xMinCovar = INT_MAX, xMaxCovar = INT_MIN, yMinCovar = INT_MAX, yMaxCovar = INT_MIN; // Size + if (covarianceType) { + for (int i = 0; i < readouts->n; i++) { + pmReadout *readout = readouts->data[i]; // The readout of interest + if (!readout) { + continue; + } + psKernel *covar = readout->covariance; // Covariance matrix + xMinCovar = PS_MIN(xMinCovar, covar->xMin); + xMaxCovar = PS_MAX(xMaxCovar, covar->xMax); + yMinCovar = PS_MIN(yMinCovar, covar->yMin); + yMaxCovar = PS_MAX(yMaxCovar, covar->yMax); + } + } + psArray *hduImages = hdu->images; // Array of images in the HDU psArray *hduMasks = hdu->masks; // Array of masks in the HDU psArray *hduVariances = hdu->variances; // Array of variances in the HDU + psArray *covariances = (covarianceType ? psArrayAlloc(readouts->n) : NULL; // Covariance images for (int i = 0; i < readouts->n; i++) { pmReadout *readout = readouts->data[i]; // The readout of interest @@ -471,4 +496,17 @@ psFree(readout->variance); readout->variance = new; + } + if (readout->covariance) { + psKernel *covar = readout->covariance; // Covariance matrix + int xMin = covar->xMin, xMax = covar->xMax, yMin = covar->yMin, yMax = covar->yMax;// Size + if (xMin == xMinCovar && xMax == xMaxCovar && yMin == yMinCovar && yMax == yMaxCovar) { + variances->data[i] = psMemIncrRefCounter(readout->covariance); + } else { + psKernel *new = psKernelAlloc(xMinCovar, xMaxCovar, yMinCovar, yMaxCovar);// New covar + psImageInit(covar->image, 0); + psImageOverlaySection(new->image, covar->image, + xMinCovar - xMin, yMinCovar - yMin, "="); + variances->data[i] = new; + } } @@ -495,4 +533,21 @@ } psFree(biassecsIter); + + if (covariances) { + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell + pmChip *chip = cell->parent; // Parent chip + const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip + if (!cellName || !chipName) { + psError(PS_ERR_UNKNOWN, false, "Unable to find cell or chip name."); + return false; + } + psString name = NULL; // Name for covariance extension + psStringAppend(&name, "COVAR_%s_%s", chipName, cellName); + if (psHashLookup(hdu->covariances, name)) { + psHashRemove(hash, name); + } + psHashAdd(hdu->covariances, name, covariance); + psFree(name); + } } // Iterating over cells within the HDU psFree(iter); @@ -504,5 +559,5 @@ // Return the level that an extension applies to -static pmFPALevel extensionLevel(pmHDU *hdu // HDU to check +static pmFPALevel extensionLevel(const pmHDU *hdu // HDU to check ) { @@ -512,5 +567,5 @@ } bool mdok = true; // Status of MD lookup - psMetadata *file = psMetadataLookupMetadata(&mdok, hdu->format, "FILE"); // File information for camera format + psMetadata *file = psMetadataLookupMetadata(&mdok, hdu->format, "FILE"); // File info for camera format if (!mdok || !file) { psError(PS_ERR_UNEXPECTED_NULL, true, "Can't file FILE information for camera format " @@ -581,5 +636,5 @@ return false; } - if (hdu->images && hdu->masks && hdu->variances) { + if (hdu->images && hdu->masks && hdu->variances && hdu->covariances) { // It's already here! return true; @@ -631,5 +686,5 @@ return generateForCells(chip); } - if (hdu->images && hdu->masks && hdu->variances) { + if (hdu->images && hdu->masks && hdu->variances && hdu->covariances) { // It's already here! return true; @@ -680,5 +735,5 @@ return generateForChips(fpa); } - if (hdu->images && hdu->masks && hdu->variances) { + if (hdu->images && hdu->masks && hdu->variances && hdu->covariances) { // It's already here! return true;