Index: trunk/psModules/src/camera/pmFPA.c =================================================================== --- trunk/psModules/src/camera/pmFPA.c (revision 21314) +++ trunk/psModules/src/camera/pmFPA.c (revision 21363) @@ -32,5 +32,6 @@ psFree(readout->image); psFree(readout->mask); - psFree(readout->weight); + psFree(readout->variance); + psFree(readout->covariance); psFree(readout->analysis); psFree(readout->bias); @@ -163,5 +164,6 @@ psFree(readout->image); psFree(readout->mask); - psFree(readout->weight); + psFree(readout->variance); + psFree(readout->covariance); psFree(readout->bias); @@ -169,5 +171,6 @@ readout->image = NULL; - readout->weight = NULL; + readout->variance = NULL; + readout->covariance = NULL; readout->mask = NULL; @@ -179,9 +182,9 @@ readout->thisImageScan = 0; readout->thisMaskScan = 0; - readout->thisWeightScan = 0; + readout->thisVarianceScan = 0; readout->lastImageScan = 0; readout->lastMaskScan = 0; - readout->lastWeightScan = 0; + readout->lastVarianceScan = 0; } @@ -197,10 +200,10 @@ if (cell->hdu) { psFree(cell->hdu->images); - psFree(cell->hdu->weights); + psFree(cell->hdu->variances); psFree(cell->hdu->masks); // psFree(cell->hdu->header); cell->hdu->images = NULL; - cell->hdu->weights = NULL; + cell->hdu->variances = NULL; cell->hdu->masks = NULL; // cell->hdu->header = NULL; @@ -219,10 +222,10 @@ if (chip->hdu) { psFree(chip->hdu->images); - psFree(chip->hdu->weights); + psFree(chip->hdu->variances); psFree(chip->hdu->masks); // psFree(chip->hdu->header); chip->hdu->images = NULL; - chip->hdu->weights = NULL; + chip->hdu->variances = NULL; chip->hdu->masks = NULL; // chip->hdu->header = NULL; @@ -241,10 +244,10 @@ if (fpa->hdu) { psFree(fpa->hdu->images); - psFree(fpa->hdu->weights); + psFree(fpa->hdu->variances); psFree(fpa->hdu->masks); // psFree(fpa->hdu->header); fpa->hdu->images = NULL; - fpa->hdu->weights = NULL; + fpa->hdu->variances = NULL; fpa->hdu->masks = NULL; // fpa->hdu->header = NULL; @@ -259,5 +262,6 @@ tmpReadout->image = NULL; tmpReadout->mask = NULL; - tmpReadout->weight = NULL; + tmpReadout->variance = NULL; + tmpReadout->covariance = NULL; tmpReadout->bias = psListAlloc(NULL); tmpReadout->analysis = psMetadataAlloc(); @@ -276,9 +280,9 @@ tmpReadout->thisImageScan = 0; tmpReadout->thisMaskScan = 0; - tmpReadout->thisWeightScan = 0; + tmpReadout->thisVarianceScan = 0; tmpReadout->lastImageScan = 0; tmpReadout->lastMaskScan = 0; - tmpReadout->lastWeightScan = 0; + tmpReadout->lastVarianceScan = 0; tmpReadout->forceScan = false; Index: trunk/psModules/src/camera/pmFPA.h =================================================================== --- trunk/psModules/src/camera/pmFPA.h (revision 21314) +++ trunk/psModules/src/camera/pmFPA.h (revision 21363) @@ -6,6 +6,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-08-12 02:51:20 $ + * @version $Revision: 1.26 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:24 $ * Copyright 2005-2006 Institute for Astronomy, University of Hawaii */ @@ -113,5 +113,5 @@ /// /// A readout corresponds to an individual read of a cell (e.g., a single image as part of a video sequence, -/// or one of multiple coadds). It contains the actual pixels used in analysis (along with mask and weight +/// or one of multiple coadds). It contains the actual pixels used in analysis (along with mask and variance /// maps). When reading from a FITS file, the images are subimages (from CELL.TRIMSEC) of the pixels read /// from the appropriate HDU (at the FPA, chip or cell level). The readout also contains a list of bias @@ -123,5 +123,6 @@ psImage *image; ///< Imaging area of readout (corresponds to CELL.TRIMSEC region) psImage *mask; ///< Mask of input image (corresponds to CELL.TRIMSEC region) - psImage *weight; ///< Weight of input image (corresponds to CELL.TRIMSEC region) + psImage *variance; ///< Variance of input image (corresponds to CELL.TRIMSEC region) + psKernel *covariance; ///< Covariance pseudo-matrix (covariance factors for single pixel) psList *bias; ///< List of bias (prescan/overscan) images psMetadata *analysis; ///< Readout-level analysis metadata @@ -130,10 +131,10 @@ bool file_exists; ///< Does the file for this readout exist (read case only)? bool data_exists; ///< Does the data for this readout exist (read case only)? - int thisImageScan; ///< start scan for next/current read - int lastImageScan; ///< start scan of the last read - int thisMaskScan; ///< start scan for next/current read - int lastMaskScan; ///< start scan of the last read - int thisWeightScan; ///< start scan for next/current read - int lastWeightScan; ///< start scan of the last read + int thisImageScan; ///< start scan for next/current read of image + int lastImageScan; ///< start scan of the last read of image + int thisMaskScan; ///< start scan for next/current read of mask + int lastMaskScan; ///< start scan of the last read of mask + int thisVarianceScan; ///< start scan for next/current read of variance + int lastVarianceScan; ///< start scan of the last read of variance bool forceScan; ///< Force pmFPARead to obey the above commanded this and last scans. } pmReadout; @@ -217,10 +218,12 @@ } \ } \ - psImage *weight = (READOUT)->weight; /* Weight map pixels */ \ - if (weight) { \ - PS_ASSERT_IMAGE_NON_NULL((READOUT)->weight, RETVAL); \ - if ((numCols != 0 || numRows != 0) && (weight->numCols != numCols || weight->numRows != numRows)) { \ - psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Weight in readout %s has wrong size (%dx%d vs %dx%d)", \ - #READOUT, weight->numCols, weight->numRows, numCols, numRows); \ + psImage *variance = (READOUT)->variance; /* Variance map pixels */ \ + if (variance) { \ + PS_ASSERT_IMAGE_NON_NULL((READOUT)->variance, RETVAL); \ + if ((numCols != 0 || numRows != 0) && \ + (variance->numCols != numCols || variance->numRows != numRows)) { \ + psError(PS_ERR_BAD_PARAMETER_SIZE, true, \ + "Variance in readout %s has wrong size (%dx%d vs %dx%d)", \ + #READOUT, variance->numCols, variance->numRows, numCols, numRows); \ return RETVAL; \ } \ @@ -236,7 +239,11 @@ PS_ASSERT_IMAGE_NON_NULL((READOUT)->mask, RETVAL); -/// Assert that a readout contains a weight map -#define PM_ASSERT_READOUT_WEIGHT(READOUT, RETVAL) \ - PS_ASSERT_IMAGE_NON_NULL((READOUT)->weight, RETVAL); +/// Assert that a readout contains a variance map +#define PM_ASSERT_READOUT_VARIANCE(READOUT, RETVAL) \ + PS_ASSERT_IMAGE_NON_NULL((READOUT)->variance, RETVAL); + +/// Assert that a readout contains a covariance matrix +#define PM_ASSERT_READOUT_COVARIANCE(READOUT, RETVAL) \ + PS_ASSERT_KERNEL_NON_NULL((READOUT)->covariance, RETVAL); /// @} Index: trunk/psModules/src/camera/pmFPAConstruct.c =================================================================== --- trunk/psModules/src/camera/pmFPAConstruct.c (revision 21314) +++ trunk/psModules/src/camera/pmFPAConstruct.c (revision 21363) @@ -1499,11 +1499,12 @@ psImage *image = readout->image; // The image psImage *mask = readout->mask; // The mask - psImage *weight = readout->weight; // The weight + psImage *variance = readout->variance; // The variance psList *bias = readout->bias; // The list of bias images if (image) { INDENT(fd, 7); - fprintf(fd, "Image: [%d:%d,%d:%d] (%dx%d)\n", image->col0, image->col0 + - image->numCols, image->row0, image->row0 + image->numRows, image->numCols, - image->numRows); + fprintf(fd, "Image: [%d:%d,%d:%d] (%dx%d)\n", + image->col0, image->col0 + image->numCols, + image->row0, image->row0 + image->numRows, + image->numCols, image->numRows); } if (bias) { @@ -1512,7 +1513,8 @@ while ((biasImage = psListGetAndIncrement(biasIter))) { INDENT(fd, 7); - fprintf(fd, "Bias: [%d:%d,%d:%d] (%dx%d)\n", biasImage->col0, - biasImage->col0 + biasImage->numCols, biasImage->row0, - biasImage->row0 + biasImage->numRows, biasImage->numCols, biasImage->numRows); + fprintf(fd, "Bias: [%d:%d,%d:%d] (%dx%d)\n", + biasImage->col0, biasImage->col0 + biasImage->numCols, + biasImage->row0, biasImage->row0 + biasImage->numRows, + biasImage->numCols, biasImage->numRows); } psFree(biasIter); @@ -1520,13 +1522,15 @@ if (mask) { INDENT(fd, 7); - fprintf(fd, "Mask: [%d:%d,%d:%d] (%dx%d)\n", mask->col0, mask->col0 + - mask->numCols, mask->row0, mask->row0 + mask->numRows, mask->numCols, - mask->numRows); + fprintf(fd, "Mask: [%d:%d,%d:%d] (%dx%d)\n", + mask->col0, mask->col0 + mask->numCols, + mask->row0, mask->row0 + mask->numRows, + mask->numCols, mask->numRows); } - if (weight) { + if (variance) { INDENT(fd, 7); - fprintf(fd, "Weight: [%d:%d,%d:%d] (%dx%d)\n", weight->col0, weight->col0 + - weight->numCols, weight->row0, weight->row0 + weight->numRows, weight->numCols, - weight->numRows); + fprintf(fd, "Variance: [%d:%d,%d:%d] (%dx%d)\n", + variance->col0, variance->col0 + variance->numCols, + variance->row0, variance->row0 + variance->numRows, + variance->numCols, variance->numRows); } } // Iterating over cell Index: trunk/psModules/src/camera/pmFPACopy.c =================================================================== --- trunk/psModules/src/camera/pmFPACopy.c (revision 21314) +++ trunk/psModules/src/camera/pmFPACopy.c (revision 21363) @@ -224,8 +224,9 @@ targetReadout->data_exists = sourceReadout->data_exists; - // Copy all three image components (image, mask, weight) - readoutCopyComponent(&targetReadout->image, sourceReadout->image, binning, xFlip, yFlip, pixels); - readoutCopyComponent(&targetReadout->mask, sourceReadout->mask, binning, xFlip, yFlip, pixels); - readoutCopyComponent(&targetReadout->weight, sourceReadout->weight, binning, xFlip, yFlip, pixels); + // Copy all three image components (image, mask, variance) + readoutCopyComponent(&targetReadout->image, sourceReadout->image, binning, xFlip, yFlip, pixels); + readoutCopyComponent(&targetReadout->mask, sourceReadout->mask, binning, xFlip, yFlip, pixels); + readoutCopyComponent(&targetReadout->variance, sourceReadout->variance, binning, xFlip, yFlip, + pixels); // Copy bias @@ -328,5 +329,5 @@ int xSize = psMetadataLookupS32(NULL, source->concepts, "CELL.XSIZE"); // CELL.XSIZE of source - psAssert (abs(xParity) == 1, "CELL.XPARITY not set for source"); + psAssert (abs(xParity) == 1, "CELL.XPARITY not set for source"); if (sourceBin == 0) { @@ -478,8 +479,8 @@ pmCell *targetCell = pmCellAlloc (targetChip, cellName); int xParityTarget = psMetadataLookupS32(&status, sourceCell->concepts, "CELL.XPARITY"); // Target x parity - psAssert (abs(xParityTarget) == 1, "CELL.XPARITY not set for target"); + psAssert (abs(xParityTarget) == 1, "CELL.XPARITY not set for target"); psMetadataAddS32 (targetCell->concepts, PS_LIST_TAIL, "CELL.XPARITY", PS_META_REPLACE, "", xParityTarget); int yParityTarget = psMetadataLookupS32(&status, sourceCell->concepts, "CELL.YPARITY"); // Target y parity - psAssert (abs(xParityTarget) == 1, "CELL.YPARITY not set for target"); + psAssert (abs(xParityTarget) == 1, "CELL.YPARITY not set for target"); psMetadataAddS32 (targetCell->concepts, PS_LIST_TAIL, "CELL.YPARITY", PS_META_REPLACE, "", yParityTarget); if (!cellCopy(targetCell, sourceCell, true, 1, 1)) { Index: trunk/psModules/src/camera/pmFPAExtent.c =================================================================== --- trunk/psModules/src/camera/pmFPAExtent.c (revision 21314) +++ trunk/psModules/src/camera/pmFPAExtent.c (revision 21363) @@ -16,5 +16,5 @@ } if (!image) { - image = readout->weight; + image = readout->variance; } @@ -23,22 +23,22 @@ if (!image) { - pmHDU *hdu = pmHDUFromReadout (readout); - if (hdu && hdu->header) { - bool status; - xSize = psMetadataLookupS32(&status, hdu->header, "NAXIS1"); - if (!status) { - xSize = psMetadataLookupS32(&status, hdu->header, "IMNAXIS1"); - if (!status) return NULL; - } - ySize = psMetadataLookupS32(&status, hdu->header, "NAXIS2"); - if (!status) { - ySize = psMetadataLookupS32(&status, hdu->header, "IMNAXIS2"); - if (!status) return NULL; - } - } else { + pmHDU *hdu = pmHDUFromReadout (readout); + if (hdu && hdu->header) { + bool status; + xSize = psMetadataLookupS32(&status, hdu->header, "NAXIS1"); + if (!status) { + xSize = psMetadataLookupS32(&status, hdu->header, "IMNAXIS1"); + if (!status) return NULL; + } + ySize = psMetadataLookupS32(&status, hdu->header, "NAXIS2"); + if (!status) { + ySize = psMetadataLookupS32(&status, hdu->header, "IMNAXIS2"); + if (!status) return NULL; + } + } else { // Don't have anything to base the true extent on, so have to give the hardwired value (largest possible extent) - xSize = psMetadataLookupS32(NULL, readout->parent->concepts, "CELL.XSIZE"); - ySize = psMetadataLookupS32(NULL, readout->parent->concepts, "CELL.YSIZE"); - } + xSize = psMetadataLookupS32(NULL, readout->parent->concepts, "CELL.XSIZE"); + ySize = psMetadataLookupS32(NULL, readout->parent->concepts, "CELL.YSIZE"); + } return psRegionAlloc(0, xSize, 0, ySize); } @@ -70,10 +70,10 @@ // Don't have anything to base the true extent on, so have to give the hardwired value (largest possible extent) if (readouts->n == 0) { - int xSize = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE"); - int ySize = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE"); - cellExtent->x0 = 0; - cellExtent->x1 = xSize; - cellExtent->y0 = 0; - cellExtent->y1 = ySize; + int xSize = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE"); + int ySize = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE"); + cellExtent->x0 = 0; + cellExtent->x1 = xSize; + cellExtent->y0 = 0; + cellExtent->y1 = ySize; } @@ -97,20 +97,20 @@ // CELL.X0,Y0 are the coordinate of the amp on the chip, subtract size if parity flipped if (xParityCell > 0) { - cellExtent->x0 += x0; - cellExtent->x1 += x0; + cellExtent->x0 += x0; + cellExtent->x1 += x0; } else { - float x0Cell = x0 - cellExtent->x1; - float x1Cell = x0 - cellExtent->x0; - cellExtent->x0 = x0Cell; - cellExtent->x1 = x1Cell; + float x0Cell = x0 - cellExtent->x1; + float x1Cell = x0 - cellExtent->x0; + cellExtent->x0 = x0Cell; + cellExtent->x1 = x1Cell; } if (yParityCell > 0) { - cellExtent->y0 += y0; - cellExtent->y1 += y0; + cellExtent->y0 += y0; + cellExtent->y1 += y0; } else { - float y0Cell = y0 - cellExtent->y1; - float y1Cell = y0 - cellExtent->y0; - cellExtent->y0 = y0Cell; - cellExtent->y1 = y1Cell; + float y0Cell = y0 - cellExtent->y1; + float y1Cell = y0 - cellExtent->y0; + cellExtent->y0 = y0Cell; + cellExtent->y1 = y1Cell; } Index: trunk/psModules/src/camera/pmFPAMaskWeight.c =================================================================== --- trunk/psModules/src/camera/pmFPAMaskWeight.c (revision 21314) +++ trunk/psModules/src/camera/pmFPAMaskWeight.c (revision 21363) @@ -49,6 +49,6 @@ } -// Create the parent weight images that reside in the HDU -static void createParentWeights(pmHDU *hdu // The HDU for which to create +// Create the parent variance images that reside in the HDU +static void createParentVariances(pmHDU *hdu // The HDU for which to create ) { @@ -58,16 +58,16 @@ // Generate the parent mask images psArray *images = hdu->images; // Array of images - psArray *weights = hdu->weights; // Array of weight images - if (!weights) { - weights = psArrayAlloc(images->n); - hdu->weights = weights; + psArray *variances = hdu->variances; // Array of variance images + if (!variances) { + variances = psArrayAlloc(images->n); + hdu->variances = variances; } for (long i = 0; i < images->n; i++) { psImage *image = images->data[i]; // The image for this readout - if (!image || weights->data[i]) { + if (!image || variances->data[i]) { continue; } - weights->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); + variances->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); } @@ -199,5 +199,5 @@ } -bool pmReadoutSetWeight(pmReadout *readout, bool poisson) +bool pmReadoutSetVariance(pmReadout *readout, bool poisson) { PS_ASSERT_PTR_NON_NULL(readout, false); @@ -209,10 +209,10 @@ float gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); // Cell gain if (!mdok || isnan(gain)) { - psError(PS_ERR_IO, true, "CELL.GAIN is not set --- unable to set weight.\n"); + psError(PS_ERR_IO, true, "CELL.GAIN is not set --- unable to set variance.\n"); return false; } float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise if (!mdok || isnan(readnoise)) { - psError(PS_ERR_IO, true, "CELL.READNOISE is not set --- unable to set weight.\n"); + psError(PS_ERR_IO, true, "CELL.READNOISE is not set --- unable to set variance.\n"); return false; } @@ -223,21 +223,21 @@ if (poisson) { - // Set weight image to the variance in ADU = f/g + rn^2 + // Set variance image to the variance in ADU = f/g + rn^2 psImage *image = readout->image; // The image pixels - readout->weight = (psImage*)psBinaryOp(readout->weight, image, "/", psScalarAlloc(gain, PS_TYPE_F32)); - - // a negative weight is non-sensical. if the image value drops below 1, the weight must be 1. - readout->weight = (psImage*)psUnaryOp(readout->weight, readout->weight, "abs"); - readout->weight = (psImage*)psBinaryOp(readout->weight, readout->weight, "max", + readout->variance = (psImage*)psBinaryOp(readout->variance, image, "/", psScalarAlloc(gain, PS_TYPE_F32)); + + // a negative variance is non-sensical. if the image value drops below 1, the variance must be 1. + readout->variance = (psImage*)psUnaryOp(readout->variance, readout->variance, "abs"); + readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "max", psScalarAlloc(1, PS_TYPE_F32)); } else { // Just use the read noise - if (!readout->weight) { - readout->weight = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32); - } - psImageInit(readout->weight, 0.0); - } - - readout->weight = (psImage*)psBinaryOp(readout->weight, readout->weight, "+", + if (!readout->variance) { + readout->variance = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32); + } + psImageInit(readout->variance, 0.0); + } + + readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32)); @@ -245,7 +245,7 @@ } -// this function creates the weight pixels, or uses the existing weight pixels. it will set -// the noise pixel values only if the weight image is not supplied -bool pmReadoutGenerateWeight(pmReadout *readout, bool poisson) +// this function creates the variance pixels, or uses the existing variance pixels. it will set +// the noise pixel values only if the variance image is not supplied +bool pmReadoutGenerateVariance(pmReadout *readout, bool poisson) { PS_ASSERT_PTR_NON_NULL(readout, false); @@ -254,11 +254,11 @@ bool mdok = true; // Status of MD lookup - // Create the weight image if required - if (readout->weight) + // Create the variance image if required + if (readout->variance) return true; psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section if (!mdok || psRegionIsNaN(*trimsec)) { - psError(PS_ERR_IO, true, "CELL.TRIMSEC is not set --- unable to set weight.\n"); + psError(PS_ERR_IO, true, "CELL.TRIMSEC is not set --- unable to set variance.\n"); return false; } @@ -271,5 +271,5 @@ } - createParentWeights(hdu); + createParentVariances(hdu); // Need to identify which readout we're working with.... @@ -280,19 +280,19 @@ } - psImage *weight = psImageSubset(hdu->weights->data[index], *trimsec); // The weight pixels - if (!weight) { + psImage *variance = psImageSubset(hdu->variances->data[index], *trimsec); // The variance pixels + if (!variance) { psString trimsecString = psRegionToString(*trimsec); - psError(PS_ERR_UNKNOWN, false, "Unable to set weight from HDU with trimsec: %s.\n", + psError(PS_ERR_UNKNOWN, false, "Unable to set variance from HDU with trimsec: %s.\n", trimsecString); psFree(trimsecString); return false; } - psImageInit(weight, 0); - readout->weight = weight; - - return pmReadoutSetWeight(readout, poisson); -} - -bool pmReadoutGenerateMaskWeight(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson) + psImageInit(variance, 0); + readout->variance = variance; + + return pmReadoutSetVariance(readout, poisson); +} + +bool pmReadoutGenerateMaskVariance(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson) { PS_ASSERT_PTR_NON_NULL(readout, false); @@ -301,10 +301,10 @@ success &= pmReadoutGenerateMask(readout, satMask, badMask); - success &= pmReadoutGenerateWeight(readout, poisson); + success &= pmReadoutGenerateVariance(readout, poisson); return success; } -bool pmCellGenerateMaskWeight(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson) +bool pmCellGenerateMaskVariance(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson) { PS_ASSERT_PTR_NON_NULL(cell, false); @@ -314,5 +314,5 @@ for (int i = 0; i < readouts->n; i++) { pmReadout *readout = readouts->data[i]; // The readout - success &= pmReadoutGenerateMaskWeight(readout, poisson, satMask, badMask); + success &= pmReadoutGenerateMaskVariance(readout, poisson, satMask, badMask); } @@ -321,12 +321,12 @@ -bool pmReadoutWeightRenormPixels(const pmReadout *readout, psImageMaskType maskVal, +bool pmReadoutVarianceRenormPixels(const pmReadout *readout, psImageMaskType maskVal, psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng) { PM_ASSERT_READOUT_NON_NULL(readout, false); PM_ASSERT_READOUT_IMAGE(readout, false); - PM_ASSERT_READOUT_WEIGHT(readout, false); - - psImage *image = readout->image, *mask = readout->mask, *weight = readout->weight; // Readout parts + PM_ASSERT_READOUT_VARIANCE(readout, false); + + psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout parts if (!psMemIncrRefCounter(rng)) { @@ -334,13 +334,13 @@ } - psStats *weightStats = psStatsAlloc(meanStat);// Statistics for mean - if (!psImageBackground(weightStats, NULL, weight, mask, maskVal, rng)) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure mean weight for image"); - psFree(weightStats); + psStats *varianceStats = psStatsAlloc(meanStat);// Statistics for mean + if (!psImageBackground(varianceStats, NULL, variance, mask, maskVal, rng)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for image"); + psFree(varianceStats); psFree(rng); return false; } - float meanVariance = weightStats->robustMedian; // Mean variance - psFree(weightStats); + float meanVariance = varianceStats->robustMedian; // Mean variance + psFree(varianceStats); psStats *imageStats = psStatsAlloc(stdevStat);// Statistics for mean @@ -356,6 +356,6 @@ float correction = PS_SQR(stdevImage) / meanVariance; // Correction to take variance to what it should be - psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising weight map by %f", correction); - psBinaryOp(weight, weight, "*", psScalarAlloc(correction, PS_TYPE_F32)); + psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", correction); + psBinaryOp(variance, variance, "*", psScalarAlloc(correction, PS_TYPE_F32)); return true; @@ -363,10 +363,10 @@ -bool pmReadoutWeightRenormPhot(const pmReadout *readout, psImageMaskType maskVal, int num, float width, +bool pmReadoutVarianceRenormPhot(const pmReadout *readout, psImageMaskType maskVal, int num, float width, psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng) { PM_ASSERT_READOUT_NON_NULL(readout, false); PM_ASSERT_READOUT_IMAGE(readout, false); - PM_ASSERT_READOUT_WEIGHT(readout, false); + PM_ASSERT_READOUT_VARIANCE(readout, false); if (!psMemIncrRefCounter(rng)) { @@ -374,5 +374,5 @@ } - psImage *image = readout->image, *mask = readout->mask, *weight = readout->weight; // Readout images + psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images // Measure background @@ -424,10 +424,10 @@ float sumNoise = 0.0; // Sum for noise measurement float sumSource = 0.0; // Sum for source measurement - float sumWeight = 0.0; // Sum for weight measurement + float sumVariance = 0.0; // Sum for variance measurement float sumGauss = 0.0, sumGauss2 = 0.0; // Sums of Gaussian kernels for (int v = 0, y = yPix - size; v < fullSize; v++, y++) { float xSumNoise = 0.0; // Sum for noise measurement in x float xSumSource = 0.0; // Sum for source measurement in x - float xSumWeight = 0.0; // Sum for weight measurement in x + float xSumVariance = 0.0; // Sum for variance measurement in x float xSumGauss = 0.0, xSumGauss2 = 0.0; // Sums of Gaussian kernels in x float yGauss = gauss->data.F32[v]; // Value of Gaussian in y @@ -441,5 +441,5 @@ xSumNoise += value * xGauss; xSumSource += (value + peakFlux * xGauss * yGauss) * xGauss; - xSumWeight += weight->data.F32[y][x] * xGauss2; + xSumVariance += variance->data.F32[y][x] * xGauss2; xSumGauss += xGauss; xSumGauss2 += xGauss2; @@ -448,5 +448,5 @@ sumNoise += xSumNoise * yGauss; sumSource += xSumSource * yGauss; - sumWeight += xSumWeight * yGauss2; + sumVariance += xSumVariance * yGauss2; sumGauss += xSumGauss * yGauss; sumGauss2 += xSumGauss2 * yGauss2; @@ -454,16 +454,16 @@ photMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = ((isfinite(sumNoise) && isfinite(sumSource) && - isfinite(sumWeight) && sumGauss > 0 && sumGauss2 > 0) ? + isfinite(sumVariance) && sumGauss > 0 && sumGauss2 > 0) ? 0 : 0xFF); float smoothImageNoise = sumNoise / sumGauss; // Value of smoothed image pixel for noise float smoothImageSource = sumSource / sumGauss; // Value of smoothed image pixel for source - float smoothWeight = sumWeight / sumGauss2; // Value of smoothed weight pixel + float smoothVariance = sumVariance / sumGauss2; // Value of smoothed variance pixel noise->data.F32[i] = smoothImageNoise; source->data.F32[i] = smoothImageSource; - guess->data.F32[i] = (sumSource > 0) ? sigFactor * PS_SQR(smoothImageSource) / smoothWeight : 0.0; + guess->data.F32[i] = (sumSource > 0) ? sigFactor * PS_SQR(smoothImageSource) / smoothVariance : 0.0; psTrace("psModules.camera", 10, "Flux %d (%d,%d): %f, %f, %f\n", - i, xPix, yPix, smoothImageNoise, smoothImageSource, smoothWeight); + i, xPix, yPix, smoothImageNoise, smoothImageSource, smoothVariance); } psFree(gauss); @@ -516,6 +516,6 @@ psFree(photMask); - psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising weight map by %f", meanRatio); - psBinaryOp(weight, weight, "*", psScalarAlloc(meanRatio, PS_TYPE_F32)); + psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", meanRatio); + psBinaryOp(variance, variance, "*", psScalarAlloc(meanRatio, PS_TYPE_F32)); return true; @@ -523,10 +523,10 @@ -bool pmReadoutWeightRenorm(const pmReadout *readout, psImageMaskType maskVal, psStatsOptions meanStat, +bool pmReadoutVarianceRenorm(const pmReadout *readout, psImageMaskType maskVal, psStatsOptions meanStat, psStatsOptions stdevStat, int width, psRandom *rng) { PM_ASSERT_READOUT_NON_NULL(readout, false); PM_ASSERT_READOUT_IMAGE(readout, false); - PM_ASSERT_READOUT_WEIGHT(readout, false); + PM_ASSERT_READOUT_VARIANCE(readout, false); PS_ASSERT_INT_POSITIVE(width, false); @@ -535,5 +535,5 @@ } - psImage *image = readout->image, *mask = readout->mask, *weight = readout->weight; // Readout images + psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images int numCols = image->numCols, numRows = image->numRows; // Size of images int xNum = numCols / width + 1, yNum = numRows / width + 1; // Number of renormalisation regions @@ -554,9 +554,9 @@ psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest psImage *subImage = psImageSubset(image, region); // Sub-image of the image pixels - psImage *subWeight = psImageSubset(weight, region); // Sub image of the weight pixels + psImage *subVariance = psImageSubset(variance, region); // Sub image of the variance pixels psImage *subMask = mask ? psImageSubset(mask, region) : NULL; // Sub-image of the mask pixels if (!psImageBackground(stdevStats, &buffer, subImage, subMask, maskVal, rng) || - !psImageBackground(meanStats, &buffer, subWeight, subMask, maskVal, rng)) { + !psImageBackground(meanStats, &buffer, subVariance, subMask, maskVal, rng)) { // Nothing we can do about it, but don't want to keel over and die, so do our best to flag it. psString regionStr = psRegionToString(region); // String with region @@ -564,5 +564,5 @@ psFree(regionStr); psErrorClear(); - psImageInit(subWeight, NAN); + psImageInit(subVariance, NAN); if (subMask) { psImageInit(subMask, maskVal); @@ -574,9 +574,9 @@ "Region [%d:%d,%d:%d] has variance %lf, but mean of variance map is %lf\n", xMin, xMax, yMin, yMax, PS_SQR(stdev), meanVar); - psBinaryOp(subWeight, subWeight, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32)); + psBinaryOp(subVariance, subVariance, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32)); } psFree(subImage); - psFree(subWeight); + psFree(subVariance); psFree(subMask); } @@ -597,5 +597,5 @@ psImage *image = readout->image; // Readout's image - psImage *weight = readout->weight; // Readout's weight + psImage *variance = readout->variance; // Readout's variance int numCols = image->numCols, numRows = image->numRows; // Size of image @@ -607,5 +607,5 @@ for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { - if (!isfinite(image->data.F32[y][x]) || (weight && !isfinite(weight->data.F32[y][x]))) { + if (!isfinite(image->data.F32[y][x]) || (variance && !isfinite(variance->data.F32[y][x]))) { mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal; } @@ -627,5 +627,5 @@ psImageMaskType **maskData = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA; // Dereference mask psF32 **imageData = readout->image->data.F32;// Dereference image - psF32 **weightData = readout->weight ? readout->weight->data.F32 : NULL; // Dereference weight map + psF32 **varianceData = readout->variance ? readout->variance->data.F32 : NULL; // Dereference variance map for (int y = 0; y < numRows; y++) { @@ -633,6 +633,6 @@ if (maskData[y][x] & maskVal) { imageData[y][x] = NAN; - if (weightData) { - weightData[y][x] = NAN; + if (varianceData) { + varianceData[y][x] = NAN; } } @@ -656,7 +656,7 @@ psImage *image = readout->image; // Image psImage *mask = readout->mask; // Mask - psImage *weight = readout->weight; // Weight map - - psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, weight, mask, maskVal, + psImage *variance = readout->variance; // Variance map + + psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, variance, mask, maskVal, NAN, NAN, maskBad, maskPoor, poorFrac, 0); interp->shifting = false; // Turn off "exact shifts" so we get proper interpolation @@ -666,5 +666,5 @@ psPixels *pixels = psPixelsAllocEmpty(PIXELS_BUFFER); // Pixels that have been interpolated psVector *imagePix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for image - psVector *weightPix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for weight + psVector *variancePix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for variance psVector *maskPix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_IMAGE_MASK); // Corresponding values for mask // NOTE: maskPix carries the actual image mask values -- do NOT use @@ -675,7 +675,7 @@ for (int x = 0; x < numCols; x++) { if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) { - double imageValue, weightValue; // Image and weight value from interpolation + double imageValue, varianceValue; // Image and variance value from interpolation psImageMaskType maskValue = 0; // Mask value from interpolation - psImageInterpolateStatus status = psImageInterpolate(&imageValue, &weightValue, &maskValue, x, y, interp); + psImageInterpolateStatus status = psImageInterpolate(&imageValue, &varianceValue, &maskValue, x, y, interp); if (status == PS_INTERPOLATE_STATUS_ERROR || status == PS_INTERPOLATE_STATUS_OFF) { psError(PS_ERR_UNKNOWN, false, "Unable to interpolate readout at %d,%d", x, y); @@ -683,5 +683,5 @@ psFree(pixels); psFree(imagePix); - psFree(weightPix); + psFree(variancePix); psFree(maskPix); return false; @@ -694,8 +694,8 @@ pixels = psPixelsAdd(pixels, PIXELS_BUFFER, x, y); imagePix = psVectorExtend(imagePix, PIXELS_BUFFER, 1); - weightPix = psVectorExtend(weightPix, PIXELS_BUFFER, 1); + variancePix = psVectorExtend(variancePix, PIXELS_BUFFER, 1); maskPix = psVectorExtend(maskPix, PIXELS_BUFFER, 1); imagePix->data.F32[numBad] = imageValue; - weightPix->data.F32[numBad] = weightValue; + variancePix->data.F32[numBad] = varianceValue; maskPix->data.PS_TYPE_IMAGE_MASK_DATA[numBad] = maskValue; numBad++; @@ -709,5 +709,5 @@ int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of pixel image->data.F32[y][x] = imagePix->data.F32[i]; - weight->data.F32[y][x] = weightPix->data.F32[i]; + variance->data.F32[y][x] = variancePix->data.F32[i]; mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskPix->data.PS_TYPE_IMAGE_MASK_DATA[i]; } @@ -715,5 +715,5 @@ psFree(pixels); psFree(imagePix); - psFree(weightPix); + psFree(variancePix); psFree(maskPix); Index: trunk/psModules/src/camera/pmFPAMaskWeight.h =================================================================== --- trunk/psModules/src/camera/pmFPAMaskWeight.h (revision 21314) +++ trunk/psModules/src/camera/pmFPAMaskWeight.h (revision 21363) @@ -5,6 +5,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:24 $ * Copyright 2005-2006 Institute for Astronomy, University of Hawaii */ @@ -47,13 +47,13 @@ -/// Set a temporary readout weight map using CELL.GAIN and CELL.READNOISE +/// Set a temporary readout variance map using CELL.GAIN and CELL.READNOISE /// -/// Calculates weights (actually variances) for each pixel using photon statistics and the cell gain -/// (CELL.GAIN) and read noise (CELL.READNOISE). The weight map that is produced within the readout is -/// temporary --- it is not added to the HDU. This is intended for when the user is iterating using -/// pmReadoutReadNext, in which case the HDU can't be generated. -bool pmReadoutSetWeight(pmReadout *readout, ///< Readout for which to set weight - bool poisson ///< Use poisson weights (in addition to read noise)? - ); +/// Calculates variances for each pixel using photon statistics and the cell gain (CELL.GAIN) and read noise +/// (CELL.READNOISE). The weight map that is produced within the readout is temporary --- it is not added to +/// the HDU. This is intended for when the user is iterating using pmReadoutReadNext, in which case the HDU +/// can't be generated. +bool pmReadoutSetVariance(pmReadout *readout, ///< Readout for which to set variance + bool poisson ///< Include poisson variance (in addition to read noise)? + ); /// Generate a readout mask (suitable for output) using CELL.SATURATION and CELL.BAD @@ -66,35 +66,35 @@ ); -/// Generate a weight map (suitable for output) using CELL.GAIN and CELL.READNOISE +/// Generate a variance map (suitable for output) using CELL.GAIN and CELL.READNOISE /// -/// Calculates weights (actually variances) for each pixel using photon statistics and the cell gain -/// (CELL.GAIN) and read noise (CELL.READNOISE). The weight map that is produced within the readout is -/// suitable for output (complete with HDU entry). This is intended for most operations. -bool pmReadoutGenerateWeight(pmReadout *readout, ///< Readout for which to generate weight - bool poisson ///< Use poisson weights (in addition to read noise)? - ); +/// Calculates variances for each pixel using photon statistics and the cell gain (CELL.GAIN) and read noise +/// (CELL.READNOISE). The variance map that is produced within the readout is suitable for output (complete +/// with HDU entry). This is intended for most operations. +bool pmReadoutGenerateVariance(pmReadout *readout, ///< Readout for which to generate variance + bool poisson ///< Include poisson variance (in addition to read noise)? + ); -/// Generate mask and weight map for a readout +/// Generate mask and variance map for a readout /// -/// Calls pmReadoutGenerateMask and pmReadoutGenerateWeight for the readout -bool pmReadoutGenerateMaskWeight(pmReadout *readout, ///< Readout for which to generate mask and weights - psImageMaskType sat, ///< Mask value to give saturated pixels - psImageMaskType bad, ///< Mask value to give bad (low) pixels - bool poisson ///< Use poisson weights (in addition to read noise)? - ); +/// Calls pmReadoutGenerateMask and pmReadoutGenerateVariance for the readout +bool pmReadoutGenerateMaskVariance(pmReadout *readout, ///< Readout for which to generate mask and variance + psImageMaskType sat, ///< Mask value to give saturated pixels + psImageMaskType bad, ///< Mask value to give bad (low) pixels + bool poisson ///< Include poisson variance (in addition to read noise)? + ); -/// Generate mask and weight maps for all readouts within a cell +/// Generate mask and variance maps for all readouts within a cell /// -/// Calls pmReadoutGenerateMaskWeight for each readout within the cell. -bool pmCellGenerateMaskWeight(pmCell *cell, ///< Cell for which to generate mask and weights - psImageMaskType sat, ///< Mask value to give saturated pixels - psImageMaskType bad, ///< Mask value to give bad (low) pixels - bool poisson ///< Use poisson weights (in addition to read noise)? - ); +/// Calls pmReadoutGenerateMaskVariance for each readout within the cell. +bool pmCellGenerateMaskVariance(pmCell *cell, ///< Cell for which to generate mask and variance + psImageMaskType sat, ///< Mask value to give saturated pixels + psImageMaskType bad, ///< Mask value to give bad (low) pixels + bool poisson ///< Include poisson variance (in addition to read noise)? + ); -/// Renormalise the weight map to match the actual pixel variance +/// Renormalise the variance map to match the actual pixel variance /// -/// The weight (variance) map is adjusted so that the mean matches the actual pixel variance in the image -bool pmReadoutWeightRenormPixels( +/// The variance map is adjusted so that the mean matches the actual pixel variance in the image +bool pmReadoutVarianceRenormPixels( const pmReadout *readout, ///< Readout to normalise psImageMaskType maskVal, ///< Value to mask @@ -104,11 +104,11 @@ ); -/// Renormalise the weight map to match the actual photometry variance +/// Renormalise the variance map to match the actual photometry variance /// -/// The weight (variance) map is adjusted so that the actual significance of fake sources matches the +/// The variance map is adjusted so that the actual significance of fake sources matches the /// guestimated significance -bool pmReadoutWeightRenormPhot( +bool pmReadoutVarianceRenormPhot( const pmReadout *readout, ///< Readout to normalise - psImageMaskType maskVal, ///< Value to mask + psImageMaskType maskVal, ///< Value to mask int num, ///< Number of instances to measure over the image float width, ///< Photometry width @@ -118,14 +118,14 @@ ); -/// Renormalise the weight map to match the actual variance +/// Renormalise the variance map to match the actual variance /// /// The variance in the image is measured in patches, and the variance map is adjusted so that the mean for /// that patch corresponds. -bool pmReadoutWeightRenorm(const pmReadout *readout, // Readout to normalise - psImageMaskType maskVal, // Value to mask - psStatsOptions meanStat, // Statistic to measure the mean (of the variance map) - psStatsOptions stdevStat, // Statistic to measure the stdev (of the image) - int width, // Width of patch (pixels) - psRandom *rng // Random number generator (for sub-sampling images) +bool pmReadoutVarianceRenorm(const pmReadout *readout, // Readout to normalise + psImageMaskType maskVal, // Value to mask + psStatsOptions meanStat, // Statistic to measure the mean (of the variance map) + psStatsOptions stdevStat, // Statistic to measure the stdev (of the image) + int width, // Width of patch (pixels) + psRandom *rng // Random number generator (for sub-sampling images) ); @@ -133,5 +133,5 @@ /// /// Since unmasked non-finite pixels can occur (e.g., by out-of-range in quantisation), it is sometimes -/// necessary to mask them explicitly. Non-finite pixels in the image or weight have their mask OR-ed with +/// necessary to mask them explicitly. Non-finite pixels in the image or variance have their mask OR-ed with /// the provided value. bool pmReadoutMaskNonfinite(pmReadout *readout, ///< Readout to mask @@ -139,9 +139,9 @@ ); -/// Apply a mask to the image and weight map +/// Apply a mask to the image and variance map /// /// Unfortunately, image subtraction may result in a bi-modal image in masked areas, which can upset image /// statistics (very important for quantising images so that a product can be written out!). This function -/// sets masked areas to NAN in the image and weight. +/// sets masked areas to NAN in the image and variance. bool pmReadoutMaskApply(pmReadout *readout, ///< Readout to mask psImageMaskType maskVal ///< Mask value for which to apply mask Index: trunk/psModules/src/camera/pmFPAMosaic.c =================================================================== --- trunk/psModules/src/camera/pmFPAMosaic.c (revision 21314) +++ trunk/psModules/src/camera/pmFPAMosaic.c (revision 21363) @@ -554,14 +554,14 @@ // We have to do all of the hard work ourselves switch (type) { - FILL_IN(U8); - FILL_IN(U16); - FILL_IN(U32); - FILL_IN(U64); - FILL_IN(S8); - FILL_IN(S16); - FILL_IN(S32); - FILL_IN(S64); - FILL_IN(F32); - FILL_IN(F64); + FILL_IN(U8); + FILL_IN(U16); + FILL_IN(U32); + FILL_IN(U64); + FILL_IN(S8); + FILL_IN(S16); + FILL_IN(S32); + FILL_IN(S64); + FILL_IN(F32); + FILL_IN(F64); default: psAbort("Should never get here.\n"); @@ -575,5 +575,5 @@ static bool addCell(psArray *images, // Array of images psArray *masks, // Array of masks - psArray *weights, // Array of weights + psArray *variances, // Array of variances psVector *x0, // Array of X0 psVector *y0, // Array of Y0 @@ -604,5 +604,5 @@ images = psArrayRealloc(images, index + CELL_LIST_BUFFER); masks = psArrayRealloc(masks, index + CELL_LIST_BUFFER); - weights = psArrayRealloc(weights, index + CELL_LIST_BUFFER); + variances = psArrayRealloc(variances, index + CELL_LIST_BUFFER); x0 = psVectorRealloc(x0, index+ CELL_LIST_BUFFER); y0 = psVectorRealloc(y0, index+ CELL_LIST_BUFFER); @@ -615,5 +615,5 @@ images->n = index + 1; masks->n = index + 1; - weights->n = index + 1; + variances->n = index + 1; x0->n = index + 1; y0->n = index + 1; @@ -726,5 +726,5 @@ // The images to put into the mosaic images->data[index] = psMemIncrRefCounter(readout->image); - weights->data[index] = psMemIncrRefCounter(readout->weight); + variances->data[index] = psMemIncrRefCounter(readout->variance); masks->data[index] = psMemIncrRefCounter(readout->mask); @@ -740,5 +740,5 @@ static bool chipMosaic(psImage **mosaicImage, // The mosaic image, to be returned psImage **mosaicMask, // The mosaic mask, to be returned - psImage **mosaicWeights, // The mosaic weights, to be returned + psImage **mosaicVariances, // The mosaic variances, to be returned int *xBinChip, int *yBinChip, // The binning in x and y, to be returned const pmChip *chip, // Chip to mosaic @@ -749,5 +749,5 @@ assert(mosaicImage); assert(mosaicMask); - assert(mosaicWeights); + assert(mosaicVariances); assert(xBinChip); assert(yBinChip); @@ -756,5 +756,5 @@ psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked - psArray *weights = psArrayAlloc(0); // Array of weight images to be mosaicked + psArray *variances = psArrayAlloc(0); // Array of variance images to be mosaicked psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates @@ -802,5 +802,5 @@ continue; } - allGood &= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip, + allGood &= addCell(images, masks, variances, x0, y0, xBin, yBin, xFlip, yFlip, cell, xBinChip, yBinChip, false, x0Target, y0Target, xParityTarget, yParityTarget); @@ -826,5 +826,5 @@ if (allGood) { *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE); - *mosaicWeights = imageMosaic(weights, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE); + *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE); *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, blank); } @@ -832,5 +832,5 @@ // Clean up psFree(images); - psFree(weights); + psFree(variances); psFree(masks); psFree(xFlip); @@ -847,5 +847,5 @@ static bool fpaMosaic(psImage **mosaicImage, // The mosaic image, to be returned psImage **mosaicMask, // The mosaic mask, to be returned - psImage **mosaicWeights, // The mosaic weights, to be returned + psImage **mosaicVariances, // The mosaic variances, to be returned int *xBinFPA, int *yBinFPA, // The binning in x and y, to be returned const pmFPA *fpa, // FPA to mosaic @@ -857,5 +857,5 @@ assert(mosaicImage); assert(mosaicMask); - assert(mosaicWeights); + assert(mosaicVariances); assert(xBinFPA); assert(yBinFPA); @@ -865,5 +865,5 @@ psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked - psArray *weights = psArrayAlloc(0); // Array of weight images to be mosaicked + psArray *variances = psArrayAlloc(0); // Array of variance images to be mosaicked psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates @@ -941,5 +941,5 @@ continue; } - allGood |= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip, + allGood |= addCell(images, masks, variances, x0, y0, xBin, yBin, xFlip, yFlip, cell, xBinFPA, yBinFPA, true, x0Target, y0Target, xParityTarget, yParityTarget); @@ -960,5 +960,5 @@ if (allGood) { *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE); - *mosaicWeights = imageMosaic(weights, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE); + *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE); *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, blank); } @@ -966,5 +966,5 @@ // Clean up psFree(images); - psFree(weights); + psFree(variances); psFree(masks); psFree(xFlip); @@ -1025,5 +1025,5 @@ psImage *mosaicImage = NULL; // The mosaic image psImage *mosaicMask = NULL; // The mosaic mask - psImage *mosaicWeights = NULL; // The mosaic weights + psImage *mosaicVariances = NULL; // The mosaic variances // Find the HDU @@ -1051,8 +1051,8 @@ } } - if (hdu->weights) { - mosaicWeights = psImageSubset(hdu->weights->data[0], bounds); - if (!mosaicWeights) { - psError(PS_ERR_UNKNOWN, false, "Unable to select weight pixels.\n"); + if (hdu->variances) { + mosaicVariances = psImageSubset(hdu->variances->data[0], bounds); + if (!mosaicVariances) { + psError(PS_ERR_UNKNOWN, false, "Unable to select variance pixels.\n"); return false; } @@ -1061,5 +1061,5 @@ // Case 2 --- we need to mosaic by cut and paste psTrace("psModules.camera", 1, "Case 2 mosaicking: cut and paste.\n"); - if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicWeights, &xBin, &yBin, source, targetCell, blank)) { + if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicVariances, &xBin, &yBin, source, targetCell, blank)) { psError(PS_ERR_UNKNOWN, false, "Unable to mosaic cells.\n"); return false; @@ -1094,5 +1094,5 @@ newReadout->image = mosaicImage; newReadout->mask = mosaicMask; - newReadout->weight = mosaicWeights; + newReadout->variance = mosaicVariances; psFree(newReadout); // Drop reference @@ -1262,5 +1262,5 @@ psImage *mosaicImage = NULL; // The mosaic image psImage *mosaicMask = NULL; // The mosaic mask - psImage *mosaicWeights = NULL; // The mosaic weights + psImage *mosaicVariances = NULL; // The mosaic variances // Find the HDU @@ -1275,11 +1275,11 @@ mosaicMask = psImageSubset(hdu->masks->data[0], *fpaRegion); } - if (hdu->weights) { - mosaicWeights = psImageSubset(hdu->weights->data[0], *fpaRegion); + if (hdu->variances) { + mosaicVariances = psImageSubset(hdu->variances->data[0], *fpaRegion); } } else { // Case 2 --- we need to mosaic by cut and paste psTrace("psModules.camera", 1, "Case 2 mosaicking: cut and paste.\n"); - if (!fpaMosaic(&mosaicImage, &mosaicMask, &mosaicWeights, &xBin, &yBin, source, + if (!fpaMosaic(&mosaicImage, &mosaicMask, &mosaicVariances, &xBin, &yBin, source, targetChip, targetCell, blank)) { psError(PS_ERR_UNKNOWN, false, "Unable to mosaic chips.\n"); @@ -1338,5 +1338,5 @@ newReadout->image = mosaicImage; newReadout->mask = mosaicMask; - newReadout->weight = mosaicWeights; + newReadout->variance = mosaicVariances; psFree(newReadout); // Drop reference Index: trunk/psModules/src/camera/pmFPARead.c =================================================================== --- trunk/psModules/src/camera/pmFPARead.c (revision 21314) +++ trunk/psModules/src/camera/pmFPARead.c (revision 21363) @@ -9,4 +9,6 @@ #include +#include "pmConfig.h" +#include "pmConfigMask.h" #include "pmHDU.h" #include "pmFPA.h" @@ -28,5 +30,5 @@ FPA_READ_TYPE_IMAGE, // Read image FPA_READ_TYPE_MASK, // Read mask - FPA_READ_TYPE_WEIGHT, // Read weight map + FPA_READ_TYPE_VARIANCE, // Read variance map FPA_READ_TYPE_HEADER // Read header } fpaReadType; @@ -54,6 +56,6 @@ case FPA_READ_TYPE_MASK: return readout->thisMaskScan; - case FPA_READ_TYPE_WEIGHT: - return readout->thisWeightScan; + case FPA_READ_TYPE_VARIANCE: + return readout->thisVarianceScan; default: psAbort("Unknown read type: %x\n", type); @@ -74,7 +76,7 @@ readout->thisMaskScan = thisScan; return readout->lastMaskScan; - case FPA_READ_TYPE_WEIGHT: - readout->thisWeightScan = thisScan; - return readout->lastWeightScan; + case FPA_READ_TYPE_VARIANCE: + readout->thisVarianceScan = thisScan; + return readout->lastVarianceScan; default: psAbort("Unknown read type: %x\n", type); @@ -93,6 +95,6 @@ case FPA_READ_TYPE_MASK: return readout->lastMaskScan; - case FPA_READ_TYPE_WEIGHT: - return readout->lastWeightScan; + case FPA_READ_TYPE_VARIANCE: + return readout->lastVarianceScan; default: psAbort("Unknown read type: %x\n", type); @@ -113,7 +115,7 @@ readout->lastMaskScan = lastScan; return readout->lastMaskScan; - case FPA_READ_TYPE_WEIGHT: - readout->lastWeightScan = lastScan; - return readout->lastWeightScan; + case FPA_READ_TYPE_VARIANCE: + readout->lastVarianceScan = lastScan; + return readout->lastVarianceScan; default: psAbort("Unknown read type: %x\n", type); @@ -132,6 +134,6 @@ case FPA_READ_TYPE_MASK: return &readout->mask; - case FPA_READ_TYPE_WEIGHT: - return &readout->weight; + case FPA_READ_TYPE_VARIANCE: + return &readout->variance; default: psAbort("Unknown read type: %x\n", type); @@ -193,4 +195,25 @@ return naxis3; +} + +// Determine whether a FITS file contains covariance matrices +static bool hduCovariance(pmHDU *hdu, // Header data unit + psFits *fits // FITS file + ) +{ + if (hdu->extname && !psFitsMoveExtName(fits, hdu->extname)) { + psError(PS_ERR_IO, false, "Unable to move to extension %s", hdu->extname); + return false; + } + // Need to explicitly read the header, since the HDU may not contain the variance header + psMetadata *header = psFitsReadHeader(NULL, fits); // Header + if (!header) { + psError(PS_ERR_IO, false, "Unable to read variance header."); + return false; + } + bool mdok; // Status of MD lookup + bool covar = psMetadataLookupBool(&mdok, header, PM_HDU_COVARIANCE_KEYWORD); // Got covariance? + psFree(header); + return covar; } @@ -350,5 +373,5 @@ *target = psImageSubset(image, region); - // Get the list of overscans: only for IMAGE types (no overscan for MASK and WEIGHT) + // Get the list of overscans: only for IMAGE types (no overscan for MASK and VARIANCE) if (type == FPA_READ_TYPE_IMAGE) { if (readout->bias->n != 0) { @@ -515,11 +538,29 @@ } - // 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: { + // Need to explicitly read the header, since what's in the pmHDU may not correspond to the mask + psMetadata *header = psFitsReadHeader(NULL, fits); + if (!header) { + psError(PS_ERR_IO, false, "Unable to read mask header."); + return false; + } + if (!pmConfigMaskReadHeader(config, header)) { + psError(PS_ERR_IO, false, "Unable to determine mask bits"); + psFree(header); + return false; + } + psFree(header); + 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); } @@ -582,5 +623,5 @@ *image = readoutReadComponent(*image, fits, trimsec, readdir, thisScan, lastScan, z, bad, pixelTypes[type]); - // Read overscans only for "image" type --- weights and masks shouldn't record overscans + // Read overscans only for "image" type --- variances and masks shouldn't record overscans if (type == FPA_READ_TYPE_IMAGE) { // Blow away existing data @@ -608,5 +649,5 @@ } -// Read into an cell; this is the engine for pmCellRead, pmCellReadMask, pmCellReadWeight +// Read into an cell; this is the engine for pmCellRead, pmCellReadMask, pmCellReadVariance // Does most of the work for the reading --- reads the HDU, and portions the HDU into readouts. static bool cellRead(pmCell *cell, // Cell into which to read @@ -640,7 +681,7 @@ dataPointer = hdu->masks; break; - case FPA_READ_TYPE_WEIGHT: - hduReadFunc = pmHDUReadWeight; - dataPointer = hdu->weights; + case FPA_READ_TYPE_VARIANCE: + hduReadFunc = pmHDUReadVariance; + dataPointer = hdu->variances; break; default: @@ -680,6 +721,6 @@ imageArray = hdu->masks; break; - case FPA_READ_TYPE_WEIGHT: - imageArray = hdu->weights; + case FPA_READ_TYPE_VARIANCE: + imageArray = hdu->variances; break; default: @@ -729,5 +770,5 @@ -// Read into an chip; this is the engine for pmChipRead, pmChipReadMask, pmChipReadWeight +// Read into an chip; this is the engine for pmChipRead, pmChipReadMask, pmChipReadVariance // Iterates over component cells, reading each static bool chipRead(pmChip *chip, // Chip into which to read @@ -760,5 +801,5 @@ -// Read into an FPA; this is the engine for pmFPARead, pmFPAReadMask, pmFPAReadWeight +// Read into an FPA; this is the engine for pmFPARead, pmFPAReadMask, pmFPAReadVariance // Iterates over component chips, reading each static bool fpaRead(pmFPA *fpa, // FPA into which to read @@ -852,5 +893,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; } @@ -1091,16 +1132,16 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Reading the weight map -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -bool pmReadoutMoreWeight(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config) +// Reading the variance map +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +bool pmReadoutMoreVariance(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config) { PS_ASSERT_PTR_NON_NULL(readout, false); PS_ASSERT_FITS_NON_NULL(fits, false); - return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_WEIGHT, config); -} - -bool pmReadoutReadChunkWeight(pmReadout *readout, psFits *fits, int z, int numScans, int overlap, + return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_VARIANCE, config); +} + +bool pmReadoutReadChunkVariance(pmReadout *readout, psFits *fits, int z, int numScans, int overlap, pmConfig *config) { @@ -1110,37 +1151,58 @@ PS_ASSERT_INT_NONNEGATIVE(numScans, false); - return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_WEIGHT, config); -} - -bool pmReadoutReadWeight(pmReadout *readout, psFits *fits, int z, pmConfig *config) + return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_VARIANCE, config); +} + +bool pmReadoutReadVariance(pmReadout *readout, psFits *fits, int z, pmConfig *config) { PS_ASSERT_PTR_NON_NULL(readout, false); PS_ASSERT_FITS_NON_NULL(fits, false); - return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_WEIGHT, config); -} - -bool pmCellReadWeight(pmCell *cell, psFits *fits, pmConfig *config) + return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_VARIANCE, config); +} + +bool pmCellReadVariance(pmCell *cell, psFits *fits, pmConfig *config) { PS_ASSERT_PTR_NON_NULL(cell, false); PS_ASSERT_FITS_NON_NULL(fits, false); - return cellRead(cell, fits, config, FPA_READ_TYPE_WEIGHT); -} - -bool pmChipReadWeight(pmChip *chip, psFits *fits, pmConfig *config) + if (!cellRead(cell, fits, config, FPA_READ_TYPE_VARIANCE)) { + return false; + } + pmHDU *hdu = pmHDUFromCell(cell); // Header data unit + if (hduCovariance(hdu, fits)) { + return pmCellReadCovariance(cell, fits); + } + return true; +} + +bool pmChipReadVariance(pmChip *chip, psFits *fits, pmConfig *config) { PS_ASSERT_PTR_NON_NULL(chip, false); PS_ASSERT_FITS_NON_NULL(fits, false); - return chipRead(chip, fits, config, FPA_READ_TYPE_WEIGHT); -} - -bool pmFPAReadWeight(pmFPA *fpa, psFits *fits, pmConfig *config) + if (!chipRead(chip, fits, config, FPA_READ_TYPE_VARIANCE)) { + return false; + } + pmHDU *hdu = pmHDUFromChip(chip); // Header data unit + if (hduCovariance(hdu, fits)) { + return pmChipReadCovariance(chip, fits); + } + return true; +} + +bool pmFPAReadVariance(pmFPA *fpa, psFits *fits, pmConfig *config) { PS_ASSERT_PTR_NON_NULL(fpa, false); PS_ASSERT_FITS_NON_NULL(fits, false); - return fpaRead(fpa, fits, config, FPA_READ_TYPE_WEIGHT); + if (!fpaRead(fpa, fits, config, FPA_READ_TYPE_VARIANCE)) { + return false; + } + pmHDU *hdu = pmHDUFromFPA(fpa); // Header data unit + if (hduCovariance(hdu, fits)) { + return pmFPAReadCovariance(fpa, fits); + } + return true; } @@ -1268,2 +1330,111 @@ return numRead; } + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Reading covariance matrices +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +bool pmCellReadCovariance(pmCell *cell, psFits *fits) +{ + PS_ASSERT_PTR_NON_NULL(cell, false); + PS_ASSERT_FITS_NON_NULL(fits, false); + + 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); + + if (!psFitsMoveExtName(fits, extname)) { + psError(PS_ERR_IO, false, "Unable to move to extension %s\n", extname); + psFree(extname); + return false; + } + psFree(extname); + + psMetadata *header = psFitsReadHeader(NULL, fits); // The FITS header + if (!header) { + psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", extname); + psFree(header); + return false; + } + + bool mdok; // Status of MD lookup + int x0 = psMetadataLookupS32(&mdok, header, "COVARIANCE.CENTRE.X"); // Centre of matrix in x + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to read covariance centre"); + psFree(header); + return false; + } + int y0 = psMetadataLookupS32(&mdok, header, "COVARIANCE.CENTRE.Y"); // Centre of matrix in y + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to read covariance centre"); + psFree(header); + return false; + } + psFree(header); + + psArray *images = psFitsReadImageCube(fits, psRegionSet(0,0,0,0)); // Covariance matrices + if (!images) { + psError(PS_ERR_IO, false, "Unable to read covariance matrices for chip %s, cell %s", + chipName, cellName); + return false; + } + + psArray *readouts = cell->readouts; // Readouts of cell + if (images->n != readouts->n) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "Number of covariance matrices (%ld) doesn't match number of readouts (%ld)", + images->n, readouts->n); + psFree(images); + return false; + } + + for (int i = 0; i < readouts->n; i++) { + pmReadout *ro = readouts->data[i]; // Readout of interest + psImage *image = images->data[i]; // Image of interest + if (ro->covariance) { + psWarning("Clobbering extant covariance matrix in chip %s, cell %s, readout %d", + chipName, cellName, i); + psFree(ro->covariance); + } + ro->covariance = psKernelAllocFromImage(image, x0, y0); + } + psFree(images); + + return true; +} + + +bool pmChipReadCovariance(pmChip *chip, psFits *fits) +{ + PS_ASSERT_PTR_NON_NULL(chip, false); + PS_ASSERT_FITS_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 (!pmCellReadCovariance(cell, fits)) { + return false; + } + } + + return true; +} + + +bool pmFPAReadCovariance(pmFPA *fpa, psFits *fits) +{ + PS_ASSERT_PTR_NON_NULL(fpa, false); + PS_ASSERT_FITS_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 (!pmChipReadCovariance(chip, fits)) { + return false; + } + } + + return true; +} Index: trunk/psModules/src/camera/pmFPARead.h =================================================================== --- trunk/psModules/src/camera/pmFPARead.h (revision 21314) +++ trunk/psModules/src/camera/pmFPARead.h (revision 21363) @@ -4,6 +4,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-06-17 22:16:38 $ + * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:24 $ * Copyright 2005-2006 Institute for Astronomy, University of Hawaii */ @@ -60,9 +60,9 @@ int numRows, // The number of rows to read pmConfig *config - ); + ); /// Return the number of readouts within a cell /// -/// This function is type-independent (doesn't matter if you are interested in the image/mask/weight). +/// This function is type-independent (doesn't matter if you are interested in the image/mask/variance). int pmCellNumReadouts(pmCell *cell, psFits *fits, pmConfig *config); @@ -74,5 +74,5 @@ psFits *fits, // FITS file from which to read pmConfig *config // Configuration - ); + ); /// Read a chip @@ -82,5 +82,5 @@ psFits *fits, // FITS file from which to read pmConfig *config // Configuration - ); + ); /// Read an FPA @@ -90,5 +90,5 @@ psFits *fits, // FITS file from which to read pmConfig *config // Configuration - ); + ); // Mask functions follow @@ -126,5 +126,5 @@ psFits *fits, // FITS file from which to read pmConfig *config // Configuration - ); + ); /// Read an entire chip into the mask @@ -134,5 +134,5 @@ psFits *fits, // FITS file from which to read pmConfig *config // Configuration - ); + ); /// Read an entire FPA into the mask @@ -142,57 +142,57 @@ psFits *fits, // FITS file from which to read pmConfig *config // Configuration - ); - -// Weight functions follow - -/// Check to see if there is more to read when reading a readout incrementally into the weight -bool pmReadoutMoreWeight(pmReadout *readout, ///< Readout of interest - psFits *fits, ///< FITS file from which to read - int z, ///< Readout number/plane; zero-offset indexing - int numScans, ///< Number of scans (rows/cols) to read - pmConfig *config ///< Configuration - ); - -/// Read a chunk of a readout into the weight + ); + +// Variance functions follow + +/// Check to see if there is more to read when reading a readout incrementally into the variance +bool pmReadoutMoreVariance(pmReadout *readout, ///< Readout of interest + psFits *fits, ///< FITS file from which to read + int z, ///< Readout number/plane; zero-offset indexing + int numScans, ///< Number of scans (rows/cols) to read + pmConfig *config ///< Configuration + ); + +/// Read a chunk of a readout into the variance /// /// Allows reading the readout incrementally -bool pmReadoutReadChunkWeight(pmReadout *readout, ///< Readout of interest - psFits *fits, ///< FITS file from which to read - int z, ///< Readout number/plane; zero-offset indexing - int numScans, ///< Number of scans (rows/cols) to read - int overlap, ///< Overlap between consecutive reads - pmConfig *config ///< Configuration - ); - -/// Read the entire readout into the weight -bool pmReadoutReadWeight(pmReadout *readout, ///< Readout of interest - psFits *fits, ///< FITS file from which to read - int z, ///< Readout number/plane; zero-offset indexing - pmConfig *config ///< Configuration - ); - -/// Read an entire cell into the weight -/// -/// Same as pmCellRead, but reads into the weight element of the readouts. -bool pmCellReadWeight(pmCell *cell, // Cell to read into - psFits *fits, // FITS file from which to read - pmConfig *config // Configuration - ); - -/// Read an entire chip into the weight -/// -/// Same as pmChipRead, but reads into the weight element of the readouts. -bool pmChipReadWeight(pmChip *chip, // Chip to read into - psFits *fits, // FITS file from which to read - pmConfig *config // Configuration - ); - -/// Read an entire FPA into the weight -/// -/// Same as pmFPARead, but reads into the weight element of the readouts. -bool pmFPAReadWeight(pmFPA *fpa, // FPA to read into - psFits *fits, // FITS file from which to read - pmConfig *config // Configuration - ); +bool pmReadoutReadChunkVariance(pmReadout *readout, ///< Readout of interest + psFits *fits, ///< FITS file from which to read + int z, ///< Readout number/plane; zero-offset indexing + int numScans, ///< Number of scans (rows/cols) to read + int overlap, ///< Overlap between consecutive reads + pmConfig *config ///< Configuration + ); + +/// Read the entire readout into the variance +bool pmReadoutReadVariance(pmReadout *readout, ///< Readout of interest + psFits *fits, ///< FITS file from which to read + int z, ///< Readout number/plane; zero-offset indexing + pmConfig *config ///< Configuration + ); + +/// Read an entire cell into the variance +/// +/// Same as pmCellRead, but reads into the variance element of the readouts. +bool pmCellReadVariance(pmCell *cell, // Cell to read into + psFits *fits, // FITS file from which to read + pmConfig *config // Configuration + ); + +/// Read an entire chip into the variance +/// +/// Same as pmChipRead, but reads into the variance element of the readouts. +bool pmChipReadVariance(pmChip *chip, // Chip to read into + psFits *fits, // FITS file from which to read + pmConfig *config // Configuration + ); + +/// Read an entire FPA into the variance +/// +/// Same as pmFPARead, but reads into the variance element of the readouts. +bool pmFPAReadVariance(pmFPA *fpa, // FPA to read into + psFits *fits, // FITS file from which to read + pmConfig *config // Configuration + ); /// Read cell headers @@ -226,5 +226,5 @@ /// also read and included in the cell analysis metadata under "name.HEADER". int pmCellReadTable(pmCell *cell, ///< Cell for which to read table - psFits *fits, ///< FITS file from which the table + psFits *fits, ///< FITS file from which to read the table const char *name ///< Specifies the extension name, and target in the analysis metadata ); @@ -233,6 +233,6 @@ /// /// Iterates over component cells, calling pmCellReadTable. -int pmChipReadTable(pmChip *chip, ///< Cell for which to read table - psFits *fits, ///< FITS file from which the table +int pmChipReadTable(pmChip *chip, ///< Chip for which to read table + psFits *fits, ///< FITS file from which to read the table const char *name ///< Specifies the extension name, and target in the analysis metadata ); @@ -241,8 +241,26 @@ /// /// Iterates over component chips, calling pmChipReadTable. -int pmFPAReadTable(pmFPA *fpa, ///< Cell for which to read table - psFits *fits, ///< FITS file from which the table +int pmFPAReadTable(pmFPA *fpa, ///< FPA for which to read table + psFits *fits, ///< FITS file from which to read the table const char *name ///< Specifies the extension name, and target in the analysis metadata ); + +/// Read covariance matrices for a cell +bool pmCellReadCovariance(pmCell *cell, ///< Cell for which to read covariance matrices + psFits *fits ///< FITS file from which to read + ); + +/// Read covariance matrices for a chip +bool pmChipReadCovariance(pmChip *chip, ///< Chip for which to read covariance matrices + psFits *fits ///< FITS file from which to read + ); + +/// Read covariance matrices for a cell +bool pmFPAReadCovariance(pmFPA *fpa, ///< FPA for which to read covariance matrices + psFits *fits ///< FITS file from which to read + ); + + + /// @} #endif Index: trunk/psModules/src/camera/pmFPAWrite.c =================================================================== --- trunk/psModules/src/camera/pmFPAWrite.c (revision 21314) +++ 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; +} Index: trunk/psModules/src/camera/pmFPAWrite.h =================================================================== --- trunk/psModules/src/camera/pmFPAWrite.h (revision 21314) +++ trunk/psModules/src/camera/pmFPAWrite.h (revision 21363) @@ -4,6 +4,7 @@ * @author Paul Price, IfA * - * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-04 02:39:36 $ + * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:24 $ + * * Copyright 2005-2006 Institute for Astronomy, University of Hawaii */ @@ -24,5 +25,5 @@ psFits *fits, ///< FITS file to which to write int z ///< Image plane to write - ); + ); /// Write a cell to a FITS file @@ -36,5 +37,5 @@ pmConfig *config, ///< Configuration bool blank ///< Write a blank PHU? - ); + ); /// Write a chip to a FITS file @@ -49,5 +50,5 @@ bool blank, ///< Write a blank PHU? bool recurse ///< Recurse to lower levels? - ); + ); /// Write an FPA to a FITS file @@ -62,5 +63,5 @@ bool blank, ///< Write a blank PHU? bool recurse ///< Recurse to lower levels? - ); + ); /// Write a cell mask to a FITS file @@ -74,5 +75,5 @@ pmConfig *config, ///< Configuration bool blank ///< Write a blank PHU? - ); + ); /// Write a chip mask to a FITS file @@ -88,5 +89,5 @@ bool blank, ///< Write a blank PHU? bool recurse ///< Recurse to lower levels? - ); + ); /// Write an FPA mask to a FITS file @@ -102,45 +103,45 @@ bool blank, ///< Write a blank PHU? bool recurse ///< Recurse to lower levels? - ); + ); -/// Write a cell weight to a FITS file +/// Write a cell variance to a FITS file /// -/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required. A blank (i.e., image-less +/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required. A blank (i.e., image-less /// header) is written only if specifically requested. Writes the concepts to the various locations, and then -/// the HDU weight to the FITS file. This function should be called at the beginning of the output cell loop -/// with blank=true in order to produce the correct file structure. -bool pmCellWriteWeight(pmCell *cell, ///< Cell to write - psFits *fits, ///< FITS file to which to write - pmConfig *config, ///< Configuration - bool blank ///< Write a blank PHU? - ); +/// the HDU variance to the FITS file. This function should be called at the beginning of the output cell +/// loop with blank=true in order to produce the correct file structure. +bool pmCellWriteVariance(pmCell *cell, ///< Cell to write + psFits *fits, ///< FITS file to which to write + pmConfig *config, ///< Configuration + bool blank ///< Write a blank PHU? + ); -/// Write a chip weight to a FITS file +/// Write a chip variance to a FITS file /// -/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required. A blank (i.e., image-less +/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required. A blank (i.e., image-less /// header) is written only if specifically requested. Writes the concepts to the various locations, and then -/// the HDU weight to the FITS file, optionally recursing to lower levels. This function should be called at -/// the beginning of the output chip loop with blank=true and recurse=false in order to produce the correct +/// the HDU variance to the FITS file, optionally recursing to lower levels. This function should be called +/// at the beginning of the output chip loop with blank=true and recurse=false in order to produce the correct /// file structure. -bool pmChipWriteWeight(pmChip *chip, ///< Chip to write - psFits *fits, ///< FITS file to which to write - pmConfig *config, ///< Configuration - bool blank, ///< Write a blank PHU? - bool recurse ///< Recurse to lower levels? - ); +bool pmChipWriteVariance(pmChip *chip, ///< Chip to write + psFits *fits, ///< FITS file to which to write + pmConfig *config, ///< Configuration + bool blank, ///< Write a blank PHU? + bool recurse ///< Recurse to lower levels? + ); -/// Write an FPA weight to a FITS file +/// Write an FPA variance to a FITS file /// -/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required. A blank (i.e., image-less +/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required. A blank (i.e., image-less /// header) is written only if specifically requested. Writes the concepts to the various locations, and then -/// the HDU weight to the FITS file, optionally recursing to lower levels. This function should be called at -/// the beginning of the output FPA loop with blank=true and recurse=false in order to produce the correct +/// the HDU variance to the FITS file, optionally recursing to lower levels. This function should be called +/// at the beginning of the output FPA loop with blank=true and recurse=false in order to produce the correct /// file structure. -bool pmFPAWriteWeight(pmFPA *fpa, ///< FPA to write - psFits *fits, ///< FITS file to which to write - pmConfig *config, ///< Configuration - bool blank, ///< Write a blank PHU? - bool recurse ///< Recurse to lower levels? - ); +bool pmFPAWriteVariance(pmFPA *fpa, ///< FPA to write + psFits *fits, ///< FITS file to which to write + pmConfig *config, ///< Configuration + bool blank, ///< Write a blank PHU? + bool recurse ///< Recurse to lower levels? + ); @@ -153,15 +154,31 @@ const pmCell *cell, ///< Cell containing FITS table in the analysis metadata const char *name ///< Name for the table data, and the extension name - ); + ); int pmChipWriteTable(psFits *fits, ///< FITS file to which to write const pmChip *chip, ///< Chip containing cells with tables to write const char *name ///< Name for the table data, and the extension name - ); + ); int pmFPAWriteTable(psFits *fits, ///< FITS file to which to write const pmFPA *fpa, ///< FPA containing cells with tables to write const char *name ///< Name for the table data, and the extension name - ); + ); + +/// Write covariance matrix to a FITS file +/// +/// The covariance matrices for a cell are written to an independent extension, named after the chip and cell +/// name. +bool pmCellWriteCovariance(psFits *fits,///< FITS file to which to write + const pmCell *cell ///< Cell for which to write covariance + ); + +bool pmChipWriteCovariance(psFits *fits,///< FITS file to which to write + const pmChip *chip ///< Chip for which to write covariance + ); + +bool pmFPAWriteCovariance(psFits *fits,///< FITS file to which to write + const pmFPA *fpa ///< FPA for which to write covariance + ); // Update the header before writing to be consistent Index: trunk/psModules/src/camera/pmFPAfile.c =================================================================== --- trunk/psModules/src/camera/pmFPAfile.c (revision 21314) +++ trunk/psModules/src/camera/pmFPAfile.c (revision 21363) @@ -476,6 +476,6 @@ return PM_FPA_FILE_MASK; } - if (!strcasecmp(type, "WEIGHT")) { - return PM_FPA_FILE_WEIGHT; + if (!strcasecmp(type, "VARIANCE")) { + return PM_FPA_FILE_VARIANCE; } if (!strcasecmp(type, "FRINGE")) { @@ -530,6 +530,6 @@ case PM_FPA_FILE_MASK: return ("MASK"); - case PM_FPA_FILE_WEIGHT: - return ("WEIGHT"); + case PM_FPA_FILE_VARIANCE: + return ("VARIANCE"); case PM_FPA_FILE_FRINGE: return ("FRINGE"); Index: trunk/psModules/src/camera/pmFPAfile.h =================================================================== --- trunk/psModules/src/camera/pmFPAfile.h (revision 21314) +++ trunk/psModules/src/camera/pmFPAfile.h (revision 21363) @@ -4,6 +4,7 @@ * @author EAM, IfA * - * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-04 02:39:36 $ + * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:24 $ + * * Copyright 2004-2005 Institute for Astronomy, University of Hawaii */ @@ -37,5 +38,5 @@ PM_FPA_FILE_IMAGE, PM_FPA_FILE_MASK, - PM_FPA_FILE_WEIGHT, + PM_FPA_FILE_VARIANCE, PM_FPA_FILE_FRINGE, PM_FPA_FILE_DARK, Index: trunk/psModules/src/camera/pmFPAfileFitsIO.c =================================================================== --- trunk/psModules/src/camera/pmFPAfileFitsIO.c (revision 21314) +++ trunk/psModules/src/camera/pmFPAfileFitsIO.c (revision 21363) @@ -148,5 +148,5 @@ case PM_FPA_FILE_IMAGE: case PM_FPA_FILE_MASK: - case PM_FPA_FILE_WEIGHT: + case PM_FPA_FILE_VARIANCE: case PM_FPA_FILE_HEADER: case PM_FPA_FILE_FRINGE: @@ -329,9 +329,10 @@ } -bool pmFPAviewReadFitsWeight(const pmFPAview *view, pmFPAfile *file, pmConfig *config) -{ - PS_ASSERT_PTR_NON_NULL(view, false); - PS_ASSERT_PTR_NON_NULL(file, false); - return fpaViewReadFitsImage(view, file, config, pmFPAReadWeight, pmChipReadWeight, pmCellReadWeight); +bool pmFPAviewReadFitsVariance(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + return fpaViewReadFitsImage(view, file, config, pmFPAReadVariance, pmChipReadVariance, + pmCellReadVariance); } @@ -347,5 +348,6 @@ PS_ASSERT_PTR_NON_NULL(view, false); PS_ASSERT_PTR_NON_NULL(file, false); - return fpaViewReadFitsImage(view, file, config, pmFPAReadHeaderSet, pmChipReadHeaderSet, pmCellReadHeaderSet); + return fpaViewReadFitsImage(view, file, config, pmFPAReadHeaderSet, pmChipReadHeaderSet, + pmCellReadHeaderSet); } @@ -432,9 +434,10 @@ } -bool pmFPAviewWriteFitsWeight(const pmFPAview *view, pmFPAfile *file, pmConfig *config) -{ - PS_ASSERT_PTR_NON_NULL(view, false); - PS_ASSERT_PTR_NON_NULL(file, false); - return fpaViewWriteFitsImage(view, file, config, pmFPAWriteWeight, pmChipWriteWeight, pmCellWriteWeight); +bool pmFPAviewWriteFitsVariance(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + return fpaViewWriteFitsImage(view, file, config, pmFPAWriteVariance, pmChipWriteVariance, + pmCellWriteVariance); } Index: trunk/psModules/src/camera/pmFPAfileFitsIO.h =================================================================== --- trunk/psModules/src/camera/pmFPAfileFitsIO.h (revision 21314) +++ trunk/psModules/src/camera/pmFPAfileFitsIO.h (revision 21363) @@ -5,6 +5,6 @@ * @author PAP, IfA * - * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-07-17 22:38:15 $ + * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:24 $ * Copyright 2004-2005 Institute for Astronomy, University of Hawaii */ @@ -27,6 +27,6 @@ pmConfig *config ); -/// Read a weight map into the current view -bool pmFPAviewReadFitsWeight(const pmFPAview *view, ///< View specifying level of interest +/// Read a variance map into the current view +bool pmFPAviewReadFitsVariance(const pmFPAview *view, ///< View specifying level of interest pmFPAfile *file, ///< FPA file into which to read pmConfig *config @@ -57,6 +57,6 @@ ); -/// Write the weight map for the specified view -bool pmFPAviewWriteFitsWeight(const pmFPAview *view, ///< View specifying level of interest +/// Write the variance map for the specified view +bool pmFPAviewWriteFitsVariance(const pmFPAview *view, ///< View specifying level of interest pmFPAfile *file, ///< FPA file to write pmConfig *config ///< Configuration Index: trunk/psModules/src/camera/pmFPAfileIO.c =================================================================== --- trunk/psModules/src/camera/pmFPAfileIO.c (revision 21314) +++ trunk/psModules/src/camera/pmFPAfileIO.c (revision 21363) @@ -173,6 +173,6 @@ status = pmFPAviewReadFitsMask(view, file, config); break; - case PM_FPA_FILE_WEIGHT: - status = pmFPAviewReadFitsWeight(view, file, config); + case PM_FPA_FILE_VARIANCE: + status = pmFPAviewReadFitsVariance(view, file, config); break; case PM_FPA_FILE_HEADER: @@ -262,5 +262,5 @@ case PM_FPA_FILE_IMAGE: case PM_FPA_FILE_MASK: - case PM_FPA_FILE_WEIGHT: + case PM_FPA_FILE_VARIANCE: case PM_FPA_FILE_FRINGE: case PM_FPA_FILE_DARK: { @@ -427,6 +427,6 @@ status = pmFPAviewWriteFitsMask(view, file, config); break; - case PM_FPA_FILE_WEIGHT: - status = pmFPAviewWriteFitsWeight(view, file, config); + case PM_FPA_FILE_VARIANCE: + status = pmFPAviewWriteFitsVariance(view, file, config); break; case PM_FPA_FILE_HEADER: @@ -524,5 +524,5 @@ case PM_FPA_FILE_IMAGE: case PM_FPA_FILE_MASK: - case PM_FPA_FILE_WEIGHT: + case PM_FPA_FILE_VARIANCE: case PM_FPA_FILE_HEADER: case PM_FPA_FILE_FRINGE: @@ -589,5 +589,5 @@ case PM_FPA_FILE_IMAGE: case PM_FPA_FILE_MASK: - case PM_FPA_FILE_WEIGHT: + case PM_FPA_FILE_VARIANCE: case PM_FPA_FILE_HEADER: case PM_FPA_FILE_FRINGE: @@ -730,6 +730,6 @@ psString tmpName = pmConfigConvertFilename (file->filename, config, create, false); if (!tmpName) { - psError(PS_ERR_IO, false, "failed to determine real name from template for %s\n", file->filename); - return false; + psError(PS_ERR_IO, false, "failed to determine real name from template for %s\n", file->filename); + return false; } psFree (file->filename); @@ -740,5 +740,5 @@ case PM_FPA_FILE_IMAGE: case PM_FPA_FILE_MASK: - case PM_FPA_FILE_WEIGHT: + case PM_FPA_FILE_VARIANCE: case PM_FPA_FILE_HEADER: case PM_FPA_FILE_FRINGE: @@ -771,20 +771,11 @@ } - // 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 // XXX can I actually reach this with camera not specified?? - psMetadata *camera = NULL; - psString formatName = NULL; - psString cameraName = NULL; - file->format = pmConfigCameraFormatFromHeader (&camera, &cameraName, &formatName, config, phu, true); + psMetadata *camera = NULL; + psString formatName = NULL; + psString cameraName = NULL; + 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); @@ -794,18 +785,18 @@ psFree(phu); - pmFPA *newFPA = pmFPAConstruct (camera, formatName); - if (!newFPA) { - psError(PS_ERR_IO, false, "Failed to construct FPA from %s for %s", file->filename, formatName); - psFree(camera); - psFree(formatName); - return NULL; - } - psFree(camera); - psFree(formatName); - psFree(cameraName); - - // XXX this is really dangerous... - psFree (file->fpa); - file->fpa = newFPA; + pmFPA *newFPA = pmFPAConstruct (camera, formatName); + if (!newFPA) { + psError(PS_ERR_IO, false, "Failed to construct FPA from %s for %s", file->filename, formatName); + psFree(camera); + psFree(formatName); + return NULL; + } + psFree(camera); + psFree(formatName); + psFree(cameraName); + + // XXX this is really dangerous... + psFree (file->fpa); + file->fpa = newFPA; } @@ -938,5 +929,5 @@ case PM_FPA_FILE_IMAGE: case PM_FPA_FILE_MASK: - case PM_FPA_FILE_WEIGHT: + case PM_FPA_FILE_VARIANCE: case PM_FPA_FILE_DARK: case PM_FPA_FILE_FRINGE: Index: trunk/psModules/src/camera/pmHDU.c =================================================================== --- trunk/psModules/src/camera/pmHDU.c (revision 21314) +++ trunk/psModules/src/camera/pmHDU.c (revision 21363) @@ -45,5 +45,5 @@ psFree(hdu->header); psFree(hdu->images); - psFree(hdu->weights); + psFree(hdu->variances); psFree(hdu->masks); } @@ -69,5 +69,5 @@ hdu->header = NULL; hdu->images = NULL; - hdu->weights = NULL; + hdu->variances = NULL; hdu->masks = NULL; @@ -155,10 +155,10 @@ } -bool pmHDUReadWeight(pmHDU *hdu, psFits *fits) -{ - PS_ASSERT_PTR_NON_NULL(hdu, false); - PS_ASSERT_PTR_NON_NULL(fits, false); - - return hduRead(hdu, &hdu->weights, fits); +bool pmHDUReadVariance(pmHDU *hdu, psFits *fits) +{ + PS_ASSERT_PTR_NON_NULL(hdu, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + + return hduRead(hdu, &hdu->variances, fits); } @@ -235,5 +235,5 @@ } -bool pmHDUWriteWeight(pmHDU *hdu, psFits *fits, const pmConfig *config) +bool pmHDUWriteVariance(pmHDU *hdu, psFits *fits, const pmConfig *config) { PS_ASSERT_PTR_NON_NULL(hdu, false); @@ -241,5 +241,5 @@ psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Value to mask - return hduWrite(hdu, hdu->weights, hdu->masks, maskVal, fits); + return hduWrite(hdu, hdu->variances, hdu->masks, maskVal, fits); } Index: trunk/psModules/src/camera/pmHDU.h =================================================================== --- trunk/psModules/src/camera/pmHDU.h (revision 21314) +++ trunk/psModules/src/camera/pmHDU.h (revision 21363) @@ -4,6 +4,7 @@ * @author Paul Price, IfA * - * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-04 02:39:36 $ + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:24 $ + * * Copyright 2005-2006 Institute for Astronomy, University of Hawaii */ @@ -18,19 +19,19 @@ /// @{ +#define PM_HDU_COVARIANCE_KEYWORD "PS_COVAR" // FITS keyword to indicate presence of a covariance matrix + /// An instance of the FITS Header Data Unit /// -/// Of course, it is not an exact replica of a FITS HDU --- they have no mask and weight data, but these are +/// 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 *weights; ///< 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 +} pmHDU; @@ -60,9 +61,9 @@ ); -/// Read the HDU header and weight map +/// Read the HDU header and variance map /// /// Moves to the appropriate extension -bool pmHDUReadWeight(pmHDU *hdu, ///< HDU to read - psFits *fits ///< FITS file to read from +bool pmHDUReadVariance(pmHDU *hdu, ///< HDU to read + psFits *fits ///< FITS file to read from ); @@ -79,8 +80,8 @@ ); -/// Write the HDU header and weight map -bool pmHDUWriteWeight(pmHDU *hdu, ///< HDU to write - psFits *fits, ///< FITS file to write to - const pmConfig *config ///< Configuration +/// Write the HDU header and variance map +bool pmHDUWriteVariance(pmHDU *hdu, ///< HDU to write + psFits *fits, ///< FITS file to write to + const pmConfig *config ///< Configuration ); Index: trunk/psModules/src/camera/pmHDUGenerate.c =================================================================== --- trunk/psModules/src/camera/pmHDUGenerate.c (revision 21314) +++ trunk/psModules/src/camera/pmHDUGenerate.c (revision 21363) @@ -274,5 +274,5 @@ pmReadout *readout = cell->readouts->data[0]; // The first readout, as representative // The proper image, used to get the size - psImage *image = readout->image ? readout->image : (readout->mask ? readout->mask : readout->weight); + psImage *image = readout->image ? readout->image : (readout->mask ? readout->mask : readout->variance); if (!image) { continue; @@ -283,8 +283,8 @@ image->numCols, image->numRows, readout->mask->numCols, readout->mask->numRows); } - if (readout->weight && - (readout->weight->numCols != image->numCols || readout->weight->numRows != image->numRows)) { - psWarning("Image and weight have different sizes (%dx%d vs %dx%d)!\n", - image->numCols, image->numRows, readout->weight->numCols, readout->weight->numRows); + if (readout->variance && + (readout->variance->numCols != image->numCols || readout->variance->numRows != image->numRows)) { + psWarning("Image and variance have different sizes (%dx%d vs %dx%d)!\n", + image->numCols, image->numRows, readout->variance->numCols, readout->variance->numRows); } // New reference @@ -355,5 +355,5 @@ psElemType imageType = 0; // Type of readout images psElemType maskType = 0; // Type of readout masks - psElemType weightType = 0; // Type of readout weights + psElemType varianceType = 0; // Type of readout variances { psListIterator *iter = psListIteratorAlloc(cells, PS_LIST_HEAD, false); // Iterator for cells @@ -381,6 +381,6 @@ maskType = checkTypes(maskType, readout->mask->type.type); } - if (!hdu->weights && readout->weight) { - weightType = checkTypes(weightType, readout->weight->type.type); + if (!hdu->variances && readout->variance) { + varianceType = checkTypes(varianceType, readout->variance->type.type); } } @@ -388,5 +388,5 @@ psFree(iter); } - if (numReadouts == 0 || (imageType == 0 && maskType == 0 && weightType == 0)) { + if (numReadouts == 0 || (imageType == 0 && maskType == 0 && varianceType == 0)) { // Nothing from which to create an HDU psFree(cells); @@ -421,10 +421,10 @@ } } - if (weightType) { - hdu->weights = psArrayAlloc(numReadouts); + if (varianceType) { + hdu->variances = psArrayAlloc(numReadouts); for (int i = 0; i < numReadouts; i++) { - psImage *weight = psImageAlloc(xSize, ySize, weightType); - psImageInit(weight, 0.0); - hdu->weights->data[i] = weight; + psImage *variance = psImageAlloc(xSize, ySize, varianceType); + psImageInit(variance, 0.0); + hdu->variances->data[i] = variance; } } @@ -448,7 +448,8 @@ psArray *readouts = cell->readouts; // Array of readouts + psArray *hduImages = hdu->images; // Array of images in the HDU psArray *hduMasks = hdu->masks; // Array of masks in the HDU - psArray *hduWeights = hdu->weights; // Array of weights in the HDU + psArray *hduVariances = hdu->variances; // Array of variances in the HDU for (int i = 0; i < readouts->n; i++) { pmReadout *readout = readouts->data[i]; // The readout of interest @@ -467,8 +468,8 @@ readout->mask = new; } - if (readout->weight) { - psImage *new = pasteImage(hduWeights->data[i], readout->weight, trimsec); - psFree(readout->weight); - readout->weight = new; + if (readout->variance) { + psImage *new = pasteImage(hduVariances->data[i], readout->variance, trimsec); + psFree(readout->variance); + readout->variance = new; } @@ -504,5 +505,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 +513,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 +582,5 @@ return false; } - if (hdu->images && hdu->masks && hdu->weights) { + if (hdu->images && hdu->masks && hdu->variances) { // It's already here! return true; @@ -631,5 +632,5 @@ return generateForCells(chip); } - if (hdu->images && hdu->masks && hdu->weights) { + if (hdu->images && hdu->masks && hdu->variances) { // It's already here! return true; @@ -680,5 +681,5 @@ return generateForChips(fpa); } - if (hdu->images && hdu->masks && hdu->weights) { + if (hdu->images && hdu->masks && hdu->variances) { // It's already here! return true; Index: trunk/psModules/src/camera/pmHDUUtils.c =================================================================== --- trunk/psModules/src/camera/pmHDUUtils.c (revision 21314) +++ trunk/psModules/src/camera/pmHDUUtils.c (revision 21363) @@ -18,13 +18,13 @@ for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - if (!chip) continue; - if (chip->hdu) return chip->hdu; - if (!chip->cells) continue; - for (int j = 0; j < chip->cells->n; j++) { - pmCell *cell = chip->cells->data[j]; - if (!cell) continue; - if (cell->hdu) return cell->hdu; - } + pmChip *chip = fpa->chips->data[i]; + if (!chip) continue; + if (chip->hdu) return chip->hdu; + if (!chip->cells) continue; + for (int j = 0; j < chip->cells->n; j++) { + pmCell *cell = chip->cells->data[j]; + if (!cell) continue; + if (cell->hdu) return cell->hdu; + } } return NULL; @@ -158,13 +158,13 @@ INDENT(fd, level + 1); - if (hdu->weights) { - fprintf(fd, "Weights:\n"); - for (long i = 0; i < hdu->weights->n; i++) { - psImage *weight = hdu->weights->data[i]; // Weight image of interest + if (hdu->variances) { + fprintf(fd, "Variances:\n"); + for (long i = 0; i < hdu->variances->n; i++) { + psImage *variance = hdu->variances->data[i]; // Variance image of interest INDENT(fd, level + 2); - fprintf(fd, "%ld: %dx%d\n", i, weight->numCols, weight->numRows); + fprintf(fd, "%ld: %dx%d\n", i, variance->numCols, variance->numRows); } } else { - fprintf(fd, "NO weights.\n"); + fprintf(fd, "NO variances.\n"); } Index: trunk/psModules/src/camera/pmReadoutStack.c =================================================================== --- trunk/psModules/src/camera/pmReadoutStack.c (revision 21314) +++ trunk/psModules/src/camera/pmReadoutStack.c (revision 21363) @@ -12,8 +12,8 @@ // XXX should it be an error for the image to exist? psImage *pmReadoutSetAnalysisImage(pmReadout *readout, // Readout containing image - const char *name, // Name of image in analysis metadata - int numCols, int numRows, // Expected size of image - psElemType type, // Expected type of image - double init // Initial value + const char *name, // Name of image in analysis metadata + int numCols, int numRows, // Expected size of image + psElemType type, // Expected type of image + double init // Initial value ) { @@ -24,5 +24,5 @@ if (!psMetadataAddImage(readout->analysis, PS_LIST_TAIL, name, 0, "Analysis image from " __FILE__, image)) { - psAbort ("analysis image already exists"); + psAbort ("analysis image already exists"); } psImageInit(image, init); @@ -35,5 +35,5 @@ // XXX not sure why this should call psMemIncrRefCounter psImage *pmReadoutGetAnalysisImage(pmReadout *readout, // Readout containing image - const char *name // Name of image in analysis metadata + const char *name // Name of image in analysis metadata ) { @@ -78,5 +78,5 @@ // XXX for the moment, use col0, row0, numCols, numRows supplied from the outside -bool pmReadoutStackDefineOutput(pmReadout *readout, int col0, int row0, int numCols, int numRows, bool mask, bool weight, psImageMaskType blank) +bool pmReadoutStackDefineOutput(pmReadout *readout, int col0, int row0, int numCols, int numRows, bool mask, bool variance, psImageMaskType blank) { PS_ASSERT_PTR_NON_NULL(readout, false); @@ -92,15 +92,15 @@ if (mask) { - // XXX is this an error? + // XXX is this an error? if (readout->mask) return false; - readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); - psImageInit(readout->mask, blank); - } - - if (weight) { - // XXX is this an error? - if (readout->weight) return false; - readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); - psImageInit(readout->weight, NAN); + readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + psImageInit(readout->mask, blank); + } + + if (variance) { + // XXX is this an error? + if (readout->variance) return false; + readout->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(readout->variance, NAN); } @@ -162,5 +162,5 @@ bool pmReadoutUpdateSize(pmReadout *readout, int minCols, int minRows, - int numCols, int numRows, bool mask, bool weight, + int numCols, int numRows, bool mask, bool variance, psImageMaskType blank) { @@ -203,15 +203,15 @@ } - if (weight) { - if (!readout->weight) { - readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); - psImageInit(readout->weight, NAN); - } - if (readout->weight->numCols < numCols || readout->weight->numRows < numRows) { - psImage *newWeight = psImageAlloc(numCols, numRows, PS_TYPE_F32); - psImageInit(newWeight, NAN); - psImageOverlaySection(newWeight, readout->weight, readout->col0, readout->row0, "="); - psFree(readout->weight); - readout->weight = newWeight; + if (variance) { + if (!readout->variance) { + readout->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(readout->variance, NAN); + } + if (readout->variance->numCols < numCols || readout->variance->numRows < numRows) { + psImage *newVariance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(newVariance, NAN); + psImageOverlaySection(newVariance, readout->variance, readout->col0, readout->row0, "="); + psFree(readout->variance); + readout->variance = newVariance; } }