Index: trunk/psModules/src/camera/pmFPA.c =================================================================== --- trunk/psModules/src/camera/pmFPA.c (revision 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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 21352) +++ 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; } } Index: trunk/psModules/src/config/pmConfigMask.c =================================================================== --- trunk/psModules/src/config/pmConfigMask.c (revision 21352) +++ trunk/psModules/src/config/pmConfigMask.c (revision 21363) @@ -10,15 +10,15 @@ static pmConfigMaskInfo masks[] = { - { "DETECTOR", NULL, 0x00, true }, // Something is wrong with the detector - { "DARK", "DETECTOR", 0x00, true }, // Pixel doesn't dark-subtract properly - { "FLAT", "DETECTOR", 0x01, true }, // Pixel doesn't flat-field properly - { "BLANK", "DETECTOR", 0x01, true }, // Pixel doesn't contain valid data - { "RANGE", NULL, 0x00, true }, // Pixel is out-of-range of linearity - { "SAT", "RANGE", 0x01, true }, // Pixel is saturated - { "BAD", "RANGE", 0x01, true }, // Pixel is low - { "BAD.WARP", NULL, 0x01, true }, // Pixel is bad after convolution with a bad pixel - { "CR", NULL, 0x00, true }, // Pixel contains a cosmic ray - { "GHOST", NULL, 0x00, true }, // Pixel contains an optical ghost - { "POOR.WARP", NULL, 0x00, false }, // Pixel is poor after convolution with a bad pixel + { "DETECTOR", NULL, 0x00, true }, // Something is wrong with the detector + { "DARK", "DETECTOR", 0x00, true }, // Pixel doesn't dark-subtract properly + { "FLAT", "DETECTOR", 0x01, true }, // Pixel doesn't flat-field properly + { "BLANK", "DETECTOR", 0x01, true }, // Pixel doesn't contain valid data + { "RANGE", NULL, 0x00, true }, // Pixel is out-of-range of linearity + { "SAT", "RANGE", 0x01, true }, // Pixel is saturated + { "BAD", "RANGE", 0x01, true }, // Pixel is low + { "BAD.WARP", NULL, 0x01, true }, // Pixel is bad after convolution with a bad pixel + { "CR", NULL, 0x00, true }, // Pixel contains a cosmic ray + { "GHOST", NULL, 0x00, true }, // Pixel contains an optical ghost + { "POOR.WARP", NULL, 0x00, false }, // Pixel is poor after convolution with a bad pixel // "LOW" Pixel is low // "CONV" Pixel is bad after convolution with a bad pixel @@ -31,5 +31,5 @@ // bits in 8-bits of space). -// XXX this file does not have psError vs psWarning worked out correctly. some of the +// XXX this file does not have psError vs psWarning worked out correctly. some of the // failure modes should result in errors, not just warnings. @@ -40,13 +40,13 @@ bool pmConfigMaskSetInMetadata(psImageMaskType *outMaskValue, // Value of MASK.VALUE, returned - psImageMaskType *outMarkValue, // Value of MARK.VALUE, returned - psMetadata *source // Source of mask bits + psImageMaskType *outMarkValue, // Value of MARK.VALUE, returned + psMetadata *source // Source of mask bits ) { PS_ASSERT_METADATA_NON_NULL(source, false); - + // Ensure all the bad mask names exist, and set the value to catch all bad pixels psImageMaskType maskValue = 0; // Value to mask to catch all the bad pixels - psImageMaskType allMasks = 0; // Value to mask to catch all masked bits (to set MARK) + psImageMaskType allMasks = 0; // Value to mask to catch all masked bits (to set MARK) int nMasks = sizeof (masks) / sizeof (pmConfigMaskInfo); @@ -55,21 +55,21 @@ bool mdok; // Status of MD lookup psImageMaskType value = psMetadataLookupImageMaskFromGeneric(&mdok, source, masks[i].badMaskName); // Value of mask - if (!mdok) { - psWarning ("problem with mask value %s\n", masks[i].badMaskName); - } + if (!mdok) { + psWarning ("problem with mask value %s\n", masks[i].badMaskName); + } if (!value) { - if (masks[i].fallbackName) { - value = psMetadataLookupImageMaskFromGeneric(&mdok, source, masks[i].fallbackName); - } - if (!value) { - value = masks[i].defaultMaskValue; - } - psMetadataAddImageMask(source, PS_LIST_TAIL, masks[i].badMaskName, PS_META_REPLACE, NULL, value); - } - if (masks[i].isBad) { - maskValue |= value; - } - allMasks |= value; + if (masks[i].fallbackName) { + value = psMetadataLookupImageMaskFromGeneric(&mdok, source, masks[i].fallbackName); + } + if (!value) { + value = masks[i].defaultMaskValue; + } + psMetadataAddImageMask(source, PS_LIST_TAIL, masks[i].badMaskName, PS_META_REPLACE, NULL, value); + } + if (masks[i].isBad) { + maskValue |= value; + } + allMasks |= value; } @@ -107,5 +107,5 @@ // Get a mask value by name(s) psImageMaskType pmConfigMaskGetFromMetadata(psMetadata *source, // Source of masks - const char *masks // Mask values to get + const char *masks // Mask values to get ) { @@ -139,66 +139,91 @@ } -// lookup an image mask value by name from a psMetadata, without requiring the entry to +// lookup an image mask value by name from a psMetadata, without requiring the entry to // be of type psImageMaskType, but verifying that it will fit in psImageMaskType -psImageMaskType psMetadataLookupImageMaskFromGeneric (bool *status, const psMetadata *md, const char *name) { - +psImageMaskType psMetadataLookupImageMaskFromGeneric(bool *status, const psMetadata *md, const char *name) +{ + if (!md) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Metadata is NULL."); + if (status) { + *status = false; + } + return 0; + } + if (!name) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Keyword is NULL."); + if (status) { + *status = false; + } + return 0; + } *status = true; - // select the mask bit name from the header - psMetadataItem *item = psMetadataLookup (md, name); - if (!item) { - psWarning("Unable to find header keyword %s when parsing mask", name); - *status = false; - return 0; - } - - // the value may be any of the U8, U16, U32, U64 types : accept the value regardless of type size - psU64 fullValue = 0; - switch (item->type) { - case PS_DATA_U8: - fullValue = item->data.U8; - break; - case PS_DATA_U16: - fullValue = item->data.U16; - break; - case PS_DATA_U32: - fullValue = item->data.U32; - break; - case PS_DATA_U64: - fullValue = item->data.U64; - break; - case PS_DATA_S8: - fullValue = item->data.S8; - break; - case PS_DATA_S16: - fullValue = item->data.S16; - break; - case PS_DATA_S32: - fullValue = item->data.S32; - break; - case PS_DATA_S64: - fullValue = item->data.S64; - break; - default: - psWarning("Mask entry %s in metadata is not of a mask type", name); - *status = false; - return 0; - } - - // will the incoming value fit within the current image mask type? - if (fullValue > PS_MAX_IMAGE_MASK_TYPE) { - psWarning("Mask entry %s in metadata is larger than allowed by the psImageMaskType", name); - *status = false; - return 0; - } - psImageMaskType value = fullValue; - // XXX validate that value is a 2^n value? - - return value; + // select the mask bit name from the header + psMetadataItem *item = psMetadataLookup(md, name); + if (!item) { + if (status) { + *status = false; + } else { + psError(PS_ERR_BAD_PARAMETER_VALUE, "Unable to find keyword %s when parsing mask", name); + } + return 0; + } + + // the value may be any of the U8, U16, U32, U64 types : accept the value regardless of type size + psU64 fullValue = 0; + switch (item->type) { + case PS_DATA_U8: + fullValue = item->data.U8; + break; + case PS_DATA_U16: + fullValue = item->data.U16; + break; + case PS_DATA_U32: + fullValue = item->data.U32; + break; + case PS_DATA_U64: + fullValue = item->data.U64; + break; + case PS_DATA_S8: + fullValue = item->data.S8; + break; + case PS_DATA_S16: + fullValue = item->data.S16; + break; + case PS_DATA_S32: + fullValue = item->data.S32; + break; + case PS_DATA_S64: + fullValue = item->data.S64; + break; + default: + if (status) { + *status = false; + } else { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + "Mask entry %s in metadata is not of a mask type", name); + } + return 0; + } + + // will the incoming value fit within the current image mask type? + if (fullValue > PS_MAX_IMAGE_MASK_TYPE) { + if (status) { + *status = false; + } else { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Mask entry %s in metadata is larger than allowed by the psImageMaskType", name); + } + return 0; + } + psImageMaskType value = fullValue; + // XXX validate that value is a 2^n value? + + return value; } // Remove from the header keywords starting with the provided string int pmConfigMaskRemoveHeaderKeywords(psMetadata *header, // Header from which to remove keywords - const char *start // Remove keywords that start with this string + const char *start // Remove keywords that start with this string ) { @@ -229,5 +254,5 @@ return 0; } - + psImageMaskType mask = pmConfigMaskGetFromMetadata (recipe, masks); return mask; @@ -274,5 +299,5 @@ } - // How many mask values do we need to read? We raise an error if this is not found, + // How many mask values do we need to read? We raise an error if this is not found, // unless the MASK.FORCE is set to true in the camera config int nMask = psMetadataLookupS32(&status, header, "MSKNUM"); @@ -300,17 +325,17 @@ } - psImageMaskType headerValue = psMetadataLookupImageMaskFromGeneric (&status, header, valuekey); - if (!status) { + psImageMaskType headerValue = psMetadataLookupImageMaskFromGeneric (&status, header, valuekey); + if (!status) { psWarning("Failed to get mask value %s from header, skipping", valuekey); - continue; - } - - // since we may read multiple mask files, we need to warn (or error?) if any of the - // header mask values conflict with other header mask values; However, the original - // mask values from the recipe do not need to match the header values. - - // when we add a header mask value, we will also add the NAME.ALREADY entry; check for - // the NAME.ALREADY entry to see if we have previously added this mask value from a - // header. + continue; + } + + // since we may read multiple mask files, we need to warn (or error?) if any of the + // header mask values conflict with other header mask values; However, the original + // mask values from the recipe do not need to match the header values. + + // when we add a header mask value, we will also add the NAME.ALREADY entry; check for + // the NAME.ALREADY entry to see if we have previously added this mask value from a + // header. psString nameAlready = NULL; // Name of key with ".ALREADY" added @@ -318,16 +343,16 @@ bool already = psMetadataLookupBool(&status, recipe, nameAlready); // Already read this one? - bool inRecipe = false; - psImageMaskType recipeValue = psMetadataLookupImageMaskFromGeneric (&inRecipe, recipe, name); - if (!inRecipe) { + bool inRecipe = false; + psImageMaskType recipeValue = psMetadataLookupImageMaskFromGeneric (&inRecipe, recipe, name); + if (!inRecipe) { psWarning("Mask value %s is not defined in the recipe", name); - } + } if (already) { - assert (inRecipe); // XXX makes no sense for NAME.ALREADY to be in without NAME + assert (inRecipe); // XXX makes no sense for NAME.ALREADY to be in without NAME if (recipeValue != headerValue) { psWarning("New mask header value does not match previously loaded entry: %x vs %x", headerValue, recipeValue); - psMetadataAddImageMask(recipe, PS_LIST_TAIL, name, PS_META_REPLACE, "Bitmask bit value", headerValue); - // XXX alternatively, error here + psMetadataAddImageMask(recipe, PS_LIST_TAIL, name, PS_META_REPLACE, "Bitmask bit value", headerValue); + // XXX alternatively, error here } } else { @@ -369,27 +394,27 @@ while ((item = psMetadataGetAndIncrement(iter))) { - // XXX this would give a false positive for mask which include '.ALREADY' in their names - char *ptr = strstr (item->name, ".ALREADY"); - if (ptr) continue; - - psU64 fullValue = 0; + // XXX this would give a false positive for mask which include '.ALREADY' in their names + char *ptr = strstr (item->name, ".ALREADY"); + if (ptr) continue; + + psU64 fullValue = 0; switch (item->type) { - case PS_DATA_U8: - fullValue = item->data.U8; - break; - case PS_DATA_U16: - fullValue = item->data.U16; - break; - case PS_DATA_U32: - fullValue = item->data.U32; - break; - case PS_DATA_U64: - fullValue = item->data.U64; - break; - default: + case PS_DATA_U8: + fullValue = item->data.U8; + break; + case PS_DATA_U16: + fullValue = item->data.U16; + break; + case PS_DATA_U32: + fullValue = item->data.U32; + break; + case PS_DATA_U64: + fullValue = item->data.U64; + break; + default: psWarning("mask recipe entry %s is not a bit value\n", item->name); continue; } - assert (fullValue <= PS_MAX_IMAGE_MASK_TYPE); // this should have been asserted on read... + assert (fullValue <= PS_MAX_IMAGE_MASK_TYPE); // this should have been asserted on read... snprintf(namekey, 64, "MSKNAM%02d", nMask); Index: trunk/psModules/src/detrend/pmShutterCorrection.c =================================================================== --- trunk/psModules/src/detrend/pmShutterCorrection.c (revision 21352) +++ trunk/psModules/src/detrend/pmShutterCorrection.c (revision 21363) @@ -383,5 +383,5 @@ psArray *images = psArrayAlloc(num);// Array of images psArray *masks = NULL; // Array of masks - psArray *weights = NULL; // Array of weights + psArray *variances = NULL; // Array of variances psVector *exptimes = psVectorAlloc(num, PS_TYPE_F32); // Vector of exposure times @@ -392,7 +392,7 @@ masks = psArrayAlloc(num); } - if (readout->weight) + if (readout->variance) { - weights = psArrayAlloc(num); + variances = psArrayAlloc(num); } } @@ -472,24 +472,24 @@ } - psImage *weight = readout->weight; // Weight map of interest - if (weight) { - if (!weights) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Not all readouts have weights.\n"); + psImage *variance = readout->variance; // Variance map of interest + if (variance) { + if (!variances) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Not all readouts have variances.\n"); goto MEASURE_ERROR; } - weights->data[i] = psMemIncrRefCounter(weight); - - if (weight->type.type != PS_TYPE_F32) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Bad type for weights: %x\n", weight->type.type); + variances->data[i] = psMemIncrRefCounter(variance); + + if (variance->type.type != PS_TYPE_F32) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Bad type for variances: %x\n", variance->type.type); goto MEASURE_ERROR; } - if (weight->numRows != numRows || weight->numCols != numCols) { + if (variance->numRows != numRows || variance->numCols != numCols) { psError(PS_ERR_BAD_PARAMETER_SIZE, true, - "Weight sizes don't match: %dx%d vs %dx%d\n", weight->numCols, weight->numRows, + "Variance sizes don't match: %dx%d vs %dx%d\n", variance->numCols, variance->numRows, numCols, numRows); goto MEASURE_ERROR; } - } else if (weights) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Not all readouts have weights.\n"); + } else if (variances) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Not all readouts have variances.\n"); goto MEASURE_ERROR; } @@ -566,12 +566,12 @@ psTrace("psModules.detrend", 3, "Mean reference value: %f\n", meanRef); - // Check the weights - if (weights && nIter > 1) { - for (int i = 0; i < weights->n && nIter > 1; i++) { - psImage *weight = weights->data[i]; // Weight image - if (!weight) { - // We don't have weights, so no realistic errors: turn off iteration + // Check the variances + if (variances && nIter > 1) { + for (int i = 0; i < variances->n && nIter > 1; i++) { + psImage *variance = variances->data[i]; // Variance image + if (!variance) { + // We don't have variances, so no realistic errors: turn off iteration if (nIter > 0) { - psWarning("Not all images have weights --- turning iteration off.\n"); + psWarning("Not all images have variances --- turning iteration off.\n"); } nIter = 1; @@ -594,7 +594,7 @@ mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = (maskImage->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal); } - psImage *weight; // Weight image - if (weights && (weight = weights->data[i])) { - errors->data.F32[i] = sqrtf(weight->data.F32[y][x]) * refs->data.F32[i]; + psImage *variance; // Variance image + if (variances && (variance = variances->data[i])) { + errors->data.F32[i] = sqrtf(variance->data.F32[y][x]) * refs->data.F32[i]; } else { errors->data.F32[i] = sqrtf(image->data.F32[y][x]) * refs->data.F32[i]; @@ -645,5 +645,5 @@ psFree(images); psFree(masks); - psFree(weights); + psFree(variances); psFree(refs); psFree(regions); @@ -878,8 +878,8 @@ PS_ASSERT_IMAGE_SIZE(readout->mask, data->numCols, data->numRows, NULL); } - if (readout->weight) { - PS_ASSERT_IMAGE_NON_NULL(readout->weight, NULL); - PS_ASSERT_IMAGE_TYPE(readout->weight, PS_TYPE_F32, NULL); - PS_ASSERT_IMAGE_SIZE(readout->weight, data->numCols, data->numRows, NULL); + if (readout->variance) { + PS_ASSERT_IMAGE_NON_NULL(readout->variance, NULL); + PS_ASSERT_IMAGE_TYPE(readout->variance, PS_TYPE_F32, NULL); + PS_ASSERT_IMAGE_SIZE(readout->variance, data->numCols, data->numRows, NULL); } @@ -1133,6 +1133,6 @@ mask->data.PS_TYPE_VECTOR_MASK_DATA[r] = (readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskVal); } - if (readout->weight) { - errors->data.F32[r] = sqrtf(readout->weight->data.F32[yIn][xIn]) * ref; + if (readout->variance) { + errors->data.F32[r] = sqrtf(readout->variance->data.F32[yIn][xIn]) * ref; } else { // XXX guess that the input data is Poisson distributed; if we go negative, force high Index: trunk/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 21352) +++ trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 21363) @@ -36,6 +36,6 @@ #define PEAK_FLUX 1.0e4 // Peak flux for each source #define SKY_VALUE 0.0e0 // Sky value for fake image -#define WEIGHT_VAL 3.0 // Weighting for image -#define WEIGHT_FACTOR 10.0 // Factor to multiply image by to get weighting +#define VARIANCE_VAL 3.0 // Variance for image +#define VARIANCE_FACTOR 10.0 // Factor to multiply image by to get variance #define PSF_STATS PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV // Statistics options for measuring PSF #define SOURCE_FIT_ITERATIONS 100 // Number of iterations for source fitting @@ -149,5 +149,5 @@ pmModel *model = pmModelFromPSFforXY(psf, x, y, PEAK_FLUX); // Model for source psAssert (model, "failed to generate model: should this be an error or not?"); - float srcRadius = model->modelRadius(model->params, PS_SQR(WEIGHT_VAL)); // Radius for source + float srcRadius = model->modelRadius(model->params, PS_SQR(VARIANCE_VAL)); // Radius for source if (srcRadius > maxRadius) { maxRadius = srcRadius; @@ -216,11 +216,12 @@ psFree(envelope); - // XXX Setting the weight seems to be an art + // XXX Setting the variance seems to be an art // Can't set it too high so that pixels are rejected as insignificant // Can't set it too low so that it's hard to get to the minimum // Have also tried: - // *** readout->weight = (psImage*)psBinaryOp(NULL, readout->image, "*", readout->image); - // *** readout->weight = (psImage*)psBinaryOp(NULL, readout->image, "*", psScalarAlloc(WEIGHT_FACTOR, PS_TYPE_F32)); - readout->weight = (psImage*)psBinaryOp(NULL, readout->image, "+", psScalarAlloc(WEIGHT_VAL, PS_TYPE_F32)); + // *** readout->variance = (psImage*)psBinaryOp(NULL, readout->image, "*", readout->image); + // *** readout->variance = (psImage*)psBinaryOp(NULL, readout->image, "*", psScalarAlloc(VARIANCE_FACTOR, PS_TYPE_F32)); + readout->variance = (psImage*)psBinaryOp(NULL, readout->image, "+", + psScalarAlloc(VARIANCE_VAL, PS_TYPE_F32)); readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); psImageInit(readout->mask, 0); @@ -242,9 +243,9 @@ psFree(source->pixels); - psFree(source->weight); + psFree(source->variance); psFree(source->maskView); psFree(source->maskObj); source->pixels = NULL; - source->weight = NULL; + source->variance = NULL; source->maskView = NULL; source->maskObj = NULL; @@ -280,5 +281,5 @@ options->psfFieldYo = 0; - pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, WEIGHT_VAL, true); + pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, VARIANCE_VAL, true); pmPSFtry *try = pmPSFtryModel(fakes, modelName, options, 0, 0xff); Index: trunk/psModules/src/imcombine/pmReadoutCombine.c =================================================================== --- trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 21352) +++ trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 21363) @@ -37,5 +37,5 @@ params->iter = 1; params->rej = INFINITY; - params->weights = false; + params->variances = false; return params; @@ -75,5 +75,5 @@ psStringAppend(&comment, "Combining using statistic: %x", params->combine); if (!hdu->header) { - hdu->header = psMetadataAlloc(); + hdu->header = psMetadataAlloc(); } psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); @@ -82,14 +82,14 @@ // note the clipping parameters, if used if (params->combine == PS_STAT_CLIPPED_MEAN) { - psString comment = NULL; // Comment to add to header - psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma", params->iter, params->rej); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - } - - // note the use of weights - if (params->weights) { + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma", params->iter, params->rej); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + } + + // note the use of variances + if (params->variances) { psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, - "Using input weights to combine images", ""); + "Using input variances to combine images", ""); } @@ -121,5 +121,5 @@ // generate the required output images based on the specified sizes - pmReadoutStackDefineOutput(output, col0, row0, numCols, numRows, true, params->weights, params->blank); + pmReadoutStackDefineOutput(output, col0, row0, numCols, numRows, true, params->variances, params->blank); psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0); @@ -177,7 +177,7 @@ return false; } - if (params->weights && !readout->weight) { + if (params->variances && !readout->variance) { psError(PS_ERR_UNEXPECTED_NULL, true, - "Rejection based on weights requested, but no weights supplied for image %d.\n", i); + "Rejection based on variances requested, but no variances supplied for image %d.\n", i); return false; } @@ -190,10 +190,12 @@ } - // pthread_t id = pthread_self(); - // char name[64]; - // sprintf (name, "%x", (unsigned int) id); - // psTimerStart (name); - - psStatsOptions combineStdev = 0; // Statistics option for weights +#if 0 + pthread_t id = pthread_self(); + char name[64]; + sprintf(name, "%x", (unsigned int)id); + psTimerStart(name); +#endif + + psStatsOptions combineStdev = 0; // Statistics option for variances switch (params->combine) { case PS_STAT_SAMPLE_MEAN: @@ -250,10 +252,10 @@ psVectorMaskType *maskData = mask->data.PS_TYPE_VECTOR_MASK_DATA; // Dereference mask - psVector *weights = NULL; // Stack of weights - psVector *errors = NULL; // Stack of errors (sqrt of variance/weights), for psVectorStats - psF32 *weightsData = NULL; // Dereference weights - if (params->weights) { - weights = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of weights - weightsData = weights->data.F32; + psVector *variances = NULL; // Stack of variances + psVector *errors = NULL; // Stack of errors (sqrt of variance), for psVectorStats + psF32 *variancesData = NULL; // Dereference variances + if (params->variances) { + variances = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of variances + variancesData = variances->data.F32; } psVector *index = NULL; // The indices to sort the pixels @@ -279,7 +281,7 @@ psF32 **outputImage = output->image->data.F32; // Output image psImageMaskType **outputMask = output->mask->data.PS_TYPE_IMAGE_MASK_DATA; // Output mask - psF32 **outputWeight = NULL; // Output weight map - if (output->weight) { - outputWeight = output->weight->data.F32; + psF32 **outputVariance = NULL; // Output variance map + if (output->variance) { + outputVariance = output->variance->data.F32; } @@ -323,6 +325,6 @@ } - if (params->weights) { - weightsData[r] = readout->weight->data.F32[yIn][xIn]; + if (params->variances) { + variancesData[r] = readout->variance->data.F32[yIn][xIn]; } @@ -332,6 +334,6 @@ if (scale) { pixelsData[r] *= invScale->data.F32[r]; - if (params->weights) { - weightsData[r] *= invScale->data.F32[r] * invScale->data.F32[r]; + if (params->variances) { + variancesData[r] *= invScale->data.F32[r] * invScale->data.F32[r]; } } @@ -376,6 +378,6 @@ // XXXXX this step probably is very expensive : convert errors to variance everywhere? - if (params->weights) { - errors = (psVector*)psUnaryOp(errors, weights, "sqrt"); + if (params->variances) { + errors = (psVector*)psUnaryOp(errors, variances, "sqrt"); } @@ -386,6 +388,6 @@ outputImage[yOut][xOut] = NAN; outputMask[yOut][xOut] = params->blank; - if (params->weights) { - outputWeight[yOut][xOut] = NAN; + if (params->variances) { + outputVariance[yOut][xOut] = NAN; } sigma->data.F32[yOut][xOut] = NAN; @@ -393,7 +395,7 @@ outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine); outputMask[yOut][xOut] = isfinite(outputImage[yOut][xOut]) ? 0 : params->blank; - if (params->weights) { + if (params->variances) { float stdev = psStatsGetValue(stats, combineStdev); - outputWeight[yOut][xOut] = PS_SQR(stdev); // Variance + outputVariance[yOut][xOut] = PS_SQR(stdev); // Variance // XXXX this is not the correct formal error. // also, the weighted mean is not obviously the correct thing here @@ -412,5 +414,5 @@ psFree(pixels); psFree(mask); - psFree(weights); + psFree(variances); psFree(errors); psFree(stats); Index: trunk/psModules/src/imcombine/pmReadoutCombine.h =================================================================== --- trunk/psModules/src/imcombine/pmReadoutCombine.h (revision 21352) +++ trunk/psModules/src/imcombine/pmReadoutCombine.h (revision 21363) @@ -5,6 +5,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * Copyright 2004-2006 Institute for Astronomy, University of Hawaii */ @@ -22,6 +22,6 @@ typedef struct { psStatsOptions combine; ///< Statistic to use when performing the combination - psImageMaskType maskVal; ///< Mask value - psImageMaskType blank; ///< Mask value to give blank (i.e., no data) pixels + psImageMaskType maskVal; ///< Mask value + psImageMaskType blank; ///< Mask value to give blank (i.e., no data) pixels int nKeep; ///< Mimimum number of pixels to keep float fracHigh; ///< Fraction of high pixels to immediately throw @@ -29,5 +29,5 @@ int iter; ///< Number of iterations for clipping (for CLIPPED_MEAN only) float rej; ///< Rejection threshould for clipping (for CLIPPED_MEAN only) - bool weights; ///< Use the supplied weights (instead of calculated stdev)? + bool variances; ///< Use the supplied variances (instead of calculated stdev)? } pmCombineParams; @@ -41,5 +41,5 @@ /// Combine multiple readouts, applying zero and scale, with optional minmax clipping bool pmReadoutCombine(pmReadout *output,///< Output readout; altered and returned - const psArray *inputs, ///< Array of input readouts (F32 image and weight, U8 mask) + const psArray *inputs, ///< Array of input readouts const psVector *zero, ///< Zero corrections to subtract from input, or NULL const psVector *scale, ///< Scale corrections to divide into input, or NULL Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 21352) +++ trunk/psModules/src/imcombine/pmStack.c (revision 21363) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.46 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.47 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii * @@ -30,5 +30,6 @@ #define PIXEL_LIST_BUFFER 100 // Number of entries to add to pixel list at a time #define PIXEL_MAP_BUFFER 2 // Number of entries to add to pixel map at a time -//#define VARIANCE_FACTORS // Use variance factors when calculating the variances? +#define VARIANCE_FACTORS // Use variance factors when calculating the variances? +//#define ADD_VARIANCE // Allow additional variance (besides variance factor)? #define NUM_DIRECT_STDEV 5 // For less than this number of values, measure stdev directly @@ -263,5 +264,5 @@ psImage *image = data->readout->image; // Image of interest - psImage *variance = data->readout->weight; // Variance ("weight") map of interest + psImage *variance = data->readout->variance; // Variance map of interest pixelData->data.F32[num] = image->data.F32[yIn][xIn]; if (variance) { @@ -442,9 +443,9 @@ *numCols = data->readout->image->numCols; *numRows = data->readout->image->numRows; - if (data->readout->weight) { + if (data->readout->variance) { *haveVariances = true; - PS_ASSERT_IMAGE_NON_NULL(data->readout->weight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->weight, false); - PS_ASSERT_IMAGE_TYPE(data->readout->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(data->readout->variance, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->variance, false); + PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false); } *haveRejects = (data->reject != NULL); @@ -473,7 +474,7 @@ PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->mask, false); if (*haveVariances) { - PS_ASSERT_IMAGE_NON_NULL(data->readout->weight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->weight, false); - PS_ASSERT_IMAGE_TYPE(data->readout->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(data->readout->variance, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->variance, false); + PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false); } } @@ -606,14 +607,11 @@ weights->data.F32[i] = data->weight; stack->data[i] = psMemIncrRefCounter(data->readout); - // Variance factor - float vf = psMetadataLookupF32(NULL, data->readout->parent->concepts, "CELL.VARFACTOR"); // Var factor - if (!isfinite(vf)) { - psWarning("Non-finite CELL.VARFACTOR for image %d --- setting to unity.", i); - vf = 1.0; - } - varFactors->data.F32[i] = vf; + varFactors->data.F32[i] = data->readout->covariance ? + psImageCovarianceFactor(data->readout->covariance) : 1.0; +#if 0 if (isfinite(data->addVariance)) { varFactors->data.F32[i] *= data->addVariance; } +#endif if (!haveRejects && !data->inspect) { data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); @@ -649,5 +647,5 @@ psImage *combinedImage = combined->image; // Combined image psImage *combinedMask = combined->mask; // Combined mask - psImage *combinedVariance = combined->weight; // Combined variance map + psImage *combinedVariance = combined->variance; // Combined variance map psArray *pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols, @@ -702,8 +700,8 @@ } - psImage *combinedVariance = combined->weight; // Combined variance map + psImage *combinedVariance = combined->variance; // Combined variance map if (haveVariances && !combinedVariance) { - combined->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); - combinedVariance = combined->weight; + combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + combinedVariance = combined->variance; } Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 21352) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 21363) @@ -60,5 +60,5 @@ // Normalise so that the sum of the variance kernel is the square of the sum of the normal kernel - // This is required to keep the relative scaling between the image and the weight map + // This is required to keep the relative scaling between the image and the variance map psBinaryOp(out->image, out->image, "*", psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32)); @@ -287,6 +287,6 @@ // Convolve an image using FFT -static void convolveWeightFFT(psImage *target,// Place the result in here - psImage *weight, // Weight map to convolve +static void convolveVarianceFFT(psImage *target,// Place the result in here + psImage *variance, // Variance map to convolve psImage *sys, // Systematic error image psImage *mask, // Mask image @@ -302,13 +302,13 @@ bool threaded = pmSubtractionThreaded(); // Are we running threaded? - psImage *subWeight = convolveSubsetAlloc(weight, border, threaded); // Weight map + psImage *subVariance = convolveSubsetAlloc(variance, border, threaded); // Variance map psImage *subSys = convolveSubsetAlloc(sys, border, threaded); // Systematic error image psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Mask // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once - psImage *convWeight = psImageConvolveFFT(NULL, subWeight, subMask, maskVal, kernel); // Convolved weight + psImage *convVariance = psImageConvolveFFT(NULL, subVariance, subMask, maskVal, kernel); // Convolved variance psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys - convolveSubsetFree(weight, subWeight, threaded); + convolveSubsetFree(variance, subVariance, threaded); convolveSubsetFree(sys, subSys, threaded); convolveSubsetFree(mask, subMask, threaded); @@ -319,5 +319,5 @@ for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) { for (int xTarget = xMin, xSource = size; xTarget < xMax; xTarget++, xSource++) { - target->data.F32[yTarget][xTarget] = convWeight->data.F32[ySource][xSource] + + target->data.F32[yTarget][xTarget] = convVariance->data.F32[ySource][xSource] + convSys->data.F32[ySource][xSource]; } @@ -326,9 +326,9 @@ int numBytes = (xMax - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) { - memcpy(&target->data.F32[yTarget][xMin], &convWeight->data.F32[ySource][size], numBytes); - } - } - - psFree(convWeight); + memcpy(&target->data.F32[yTarget][xMin], &convVariance->data.F32[ySource][size], numBytes); + } + } + + psFree(convVariance); psFree(convSys); @@ -361,10 +361,10 @@ // Convolve a region of an image static inline void convolveRegion(psImage *convImage, // Convolved image (output) - psImage *convWeight, // Convolved weight map (output), or NULL + psImage *convVariance, // Convolved variance map (output), or NULL psImage *convMask, // Convolve mask (output), or NULL psKernel **kernelImage, // Convolution kernel for the image - psKernel **kernelWeight, // Convolution kernel for the weight map, or NULL + psKernel **kernelVariance, // Convolution kernel for the variance map, or NULL psImage *image, // Image to convolve - psImage *weight, // Weight map to convolve, or NULL + psImage *variance, // Variance map to convolve, or NULL psImage *sys, // Systematic error image, or NULL psImage *subMask, // Subtraction mask @@ -381,6 +381,6 @@ { *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual); - if (weight || subMask) { - *kernelWeight = varianceKernel(*kernelWeight, *kernelImage); + if (variance || subMask) { + *kernelVariance = varianceKernel(*kernelVariance, *kernelImage); } @@ -398,17 +398,17 @@ } - // Convolve the image and weight + // Convolve the image and variance if (useFFT) { // Use Fast Fourier Transform to do the convolution // This provides a big speed-up for large kernels convolveFFT(convImage, image, subMask, subBad, *kernelImage, region, background, kernels->size); - if (weight) { - convolveWeightFFT(convWeight, weight, sys, subMask, subBad, *kernelWeight, region, kernels->size); + if (variance) { + convolveVarianceFFT(convVariance, variance, sys, subMask, subBad, *kernelVariance, region, kernels->size); } } else { // XXX Direct convolution doesn't account for bad pixels yet convolveDirect(convImage, image, *kernelImage, region, background, kernels->size); - if (weight) { - convolveDirect(convWeight, weight, *kernelWeight, region, 0.0, kernels->size); + if (variance) { + convolveDirect(convVariance, variance, *kernelVariance, region, 0.0, kernels->size); } } @@ -885,6 +885,6 @@ psFree(stamp->image1); psFree(stamp->image2); - psFree(stamp->weight); - stamp->image1 = stamp->image2 = stamp->weight = NULL; + psFree(stamp->variance); + stamp->image1 = stamp->image2 = stamp->variance = NULL; psFree(stamp->matrix1); psFree(stamp->matrix2); @@ -915,5 +915,5 @@ } -psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, float x, float y, bool wantDual) +psKernel *pmSubtractionKernel(const pmSubtractionKernels *kernels, float x, float y, bool wantDual) { PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL); @@ -922,12 +922,19 @@ PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); - // Precalulate polynomial values - psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); - - // The appropriate kernel - psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual); - + psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); // Solved polynomial + psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual); // The appropriate kernel psFree(polyValues); + return kernel; +} + +psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, float x, float y, bool wantDual) +{ + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL); + PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL); + PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); + + psKernel *kernel = pmSubtractionKernel(kernels, x, y, wantDual); // Convolution kernel psImage *image = psMemIncrRefCounter(kernel->image); // Image of the kernel psFree(kernel); @@ -989,5 +996,5 @@ -// XXX Put kernelImage, kernelWeight and polyValues on thread-dependent data +// XXX Put kernelImage, kernelVariance and polyValues on thread-dependent data static bool subtractionConvolvePatch(int numCols, int numRows, // Size of image int x0, int y0, // Offsets for image @@ -1010,5 +1017,5 @@ psKernel *kernelImage = NULL; // Kernel for the images - psKernel *kernelWeight = NULL; // Kernel for the weight maps + psKernel *kernelVariance = NULL; // Kernel for the variance maps // Only generate polynomial values every kernel footprint, since we have already assumed @@ -1020,16 +1027,16 @@ if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - convolveRegion(out1->image, out1->weight, convMask, &kernelImage, &kernelWeight, - ro1->image, ro1->weight, sys1, subMask, kernels, polyValues, background, *region, + convolveRegion(out1->image, out1->variance, convMask, &kernelImage, &kernelVariance, + ro1->image, ro1->variance, sys1, subMask, kernels, polyValues, background, *region, maskBad, maskPoor, poorFrac, useFFT, false); } if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - convolveRegion(out2->image, out2->weight, convMask, &kernelImage, &kernelWeight, - ro2->image, ro2->weight, sys2, subMask, kernels, polyValues, background, *region, + convolveRegion(out2->image, out2->variance, convMask, &kernelImage, &kernelVariance, + ro2->image, ro2->variance, sys2, subMask, kernels, polyValues, background, *region, maskBad, maskPoor, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL); } psFree(kernelImage); - psFree(kernelWeight); + psFree(kernelVariance); psFree(polyValues); @@ -1148,12 +1155,12 @@ // XXX } } - if (ro1->weight) { - if (!out1->weight) { - out1->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + if (ro1->variance) { + if (!out1->variance) { + out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); // XXX if (threaded) { - // XXX psMutexInit(out1->weight); + // XXX psMutexInit(out1->variance); // XXX } } - psImageInit(out1->weight, 0.0); + psImageInit(out1->variance, 0.0); } } @@ -1165,12 +1172,12 @@ // XXX } } - if (ro2->weight) { - if (!out2->weight) { - out2->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + if (ro2->variance) { + if (!out2->variance) { + out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); // XXX if (threaded) { - // XXX psMutexInit(out2->weight); + // XXX psMutexInit(out2->variance); // XXX } } - psImageInit(out2->weight, 0.0); + psImageInit(out2->variance, 0.0); } } @@ -1232,5 +1239,5 @@ psImage *polyValues = NULL; // Pre-calculated polynomial values psKernel *kernelImage = NULL; // Kernel for the images - psKernel *kernelWeight = NULL; // Kernel for the weight maps + psKernel *kernelVariance = NULL; // Kernel for the variance maps #endif @@ -1299,7 +1306,7 @@ psFree(job); } else { - subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, - sys1, sys2, subMask, maskBad, maskPoor, poorFrac, subRegion, - kernels, doBG, useFFT); + subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2, + subMask, maskBad, maskPoor, poorFrac, subRegion, kernels, doBG, + useFFT); } psFree(subRegion); @@ -1334,9 +1341,23 @@ // XXX } } - psImageConvolveSetThreads(oldThreads); psFree(sys1); psFree(sys2); + + // Calculate covariances + // This can take a while, so we only do it for a single instance + // XXX psImageCovarianceCalculate could be multithreaded + if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel + out1->covariance = psImageCovarianceCalculate(kernel, ro1->covariance); + psFree(kernel); + } + if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, + kernels->mode == PM_SUBTRACTION_MODE_DUAL); // Conv. kernel + out2->covariance = psImageCovarianceCalculate(kernel, ro2->covariance); + psFree(kernel); + } // Copy anything that wasn't convolved @@ -1345,6 +1366,7 @@ if (out2) { out2->image = psMemIncrRefCounter(ro2->image); - out2->weight = psMemIncrRefCounter(ro2->weight); + out2->variance = psMemIncrRefCounter(ro2->variance); out2->mask = psMemIncrRefCounter(ro2->mask); + out2->covariance = psMemIncrRefCounter(ro2->covariance); } break; @@ -1352,6 +1374,7 @@ if (out1) { out1->image = psMemIncrRefCounter(ro1->image); - out1->weight = psMemIncrRefCounter(ro1->weight); + out1->variance = psMemIncrRefCounter(ro1->variance); out1->mask = psMemIncrRefCounter(ro1->mask); + out1->covariance = psMemIncrRefCounter(ro1->covariance); } break; Index: trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.h (revision 21352) +++ trunk/psModules/src/imcombine/pmSubtraction.h (revision 21363) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * Copyright 2004-207 Institute for Astronomy, University of Hawaii */ @@ -70,4 +70,10 @@ float sigmaRej, ///< Number of RMS deviations above zero at which to reject int footprint ///< Half-size of stamp + ); + +/// Generate the convolution kernel +psKernel *pmSubtractionKernel(const pmSubtractionKernels *kernels, ///< Kernel parameters + float x, float y, ///< Normalised position [-1,1] for which to generate image + bool wantDual ///< Calculate for the dual kernel? ); Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 21352) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 21363) @@ -24,5 +24,5 @@ static inline double calculateSumProduct(const psKernel *image1, // First image in multiplication const psKernel *image2, // Second image in multiplication - const psKernel *weight, // Weight image + const psKernel *variance, // Variance image int footprint // (Half-)Size of stamp ) @@ -31,5 +31,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // weight->kernel[y][x]; + sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // variance->kernel[y][x]; } } @@ -42,5 +42,5 @@ const psKernel *image1, // First image in multiplication const psKernel *image2, // Second image in multiplication - const psKernel *weight, // Weight image + const psKernel *variance, // Variance image const psImage *polyValues, // Spatial polynomial values int numKernels, // Number of kernel basis functions @@ -50,5 +50,5 @@ ) { - double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products + double sum = calculateSumProduct(image1, image2, variance, footprint); // Sum of the image products if (!isfinite(sum)) { return false; @@ -79,5 +79,5 @@ const psKernel *image1, // First image in multiplication const psKernel *image2, // Second image in multiplication - const psKernel *weight, // Weight image + const psKernel *variance, // Variance image const psImage *polyValues, // Spatial polynomial values int numKernels, // Number of kernel basis functions @@ -87,5 +87,5 @@ ) { - double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products + double sum = calculateSumProduct(image1, image2, variance, footprint); // Sum of the image products if (!isfinite(sum)) { return false; @@ -120,5 +120,5 @@ const psArray *convolutions1, // Convolutions for element 1 const psArray *convolutions2, // Convolutions for element 2 - const psKernel *weight, // Weight image + const psKernel *variance, // Variance image const psImage *polyValues, // Polynomial values int numKernels, // Number of kernel basis functions @@ -135,5 +135,5 @@ psKernel *jConv = convolutions2->data[j]; // Convolution for j-th element - if (!calculateMatrixElement2(matrix, i, j, iConv, jConv, weight, polyValues, numKernels, + if (!calculateMatrixElement2(matrix, i, j, iConv, jConv, variance, polyValues, numKernels, footprint, spatialOrder, symmetric)) { psTrace("psModules.imcombine", 2, "Bad sumCC at %d, %d", i, j); @@ -151,5 +151,5 @@ const psArray *convolutions, // Convolutions of source with kernels const psKernel *input, // Input stamp, or NULL - const psKernel *weight, // Weight stamp + const psKernel *variance, // Variance stamp const psImage *polyValues, // Spatial polynomial values int footprint, // (Half-)Size of stamp @@ -171,5 +171,5 @@ // Square part of the matrix (convolution-convolution products) - if (!calculateMatrixSquare(matrix, convolutions, convolutions, weight, polyValues, numKernels, + if (!calculateMatrixSquare(matrix, convolutions, convolutions, variance, polyValues, numKernels, spatialOrder, footprint)) { return false; @@ -185,5 +185,5 @@ // Normalisation-convolution terms - if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, weight, polyValues, numKernels, + if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, variance, polyValues, numKernels, footprint, spatialOrder, true)) { psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); @@ -195,5 +195,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x]; + sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; } } @@ -218,5 +218,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - double invNoise2 = 1.0 / 1.0; // weight->kernel[y][x]; + double invNoise2 = 1.0 / 1.0; // variance->kernel[y][x]; double value = input->kernel[y][x] * invNoise2; sumI += value; @@ -253,5 +253,5 @@ const psKernel *input, // Input stamp, or NULL if !normAndBG const psKernel *target, // Target stamp - const psKernel *weight, // Weight stamp + const psKernel *variance, // Variance stamp const psImage *polyValues, // Spatial polynomial values int footprint, // (Half-)Size of stamp @@ -277,5 +277,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // weight->kernel[y][x]; + sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; } } @@ -297,5 +297,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - float value = target->kernel[y][x] / 1.0; // weight->kernel[y][x]; + float value = target->kernel[y][x] / 1.0; // variance->kernel[y][x]; sumIT += value * input->kernel[y][x]; sumT += value; @@ -329,5 +329,5 @@ const psArray *convolutions2, // Convolutions of image 2 const psKernel *image1, // Image 1 stamp - const psKernel *weight, // Weight stamp + const psKernel *variance, // Variance stamp const psImage *polyValues, // Spatial polynomial values int footprint // (Half-)Size of stamp @@ -348,5 +348,5 @@ PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels); - if (!calculateMatrixSquare(matrix, convolutions1, convolutions2, weight, polyValues, numKernels, + if (!calculateMatrixSquare(matrix, convolutions1, convolutions2, variance, polyValues, numKernels, spatialOrder, footprint)) { return false; @@ -356,5 +356,5 @@ // Normalisation psKernel *conv = convolutions2->data[i]; // Convolution - if (!calculateMatrixElement1(matrix, i, normIndex, conv, image1, weight, polyValues, numKernels, + if (!calculateMatrixElement1(matrix, i, normIndex, conv, image1, variance, polyValues, numKernels, footprint, spatialOrder, false)) { psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); @@ -366,5 +366,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x]; + sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; } } @@ -559,13 +559,13 @@ case PM_SUBTRACTION_MODE_1: status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, - stamp->weight, polyValues, footprint, true); + stamp->variance, polyValues, footprint, true); status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, - stamp->image2, stamp->weight, polyValues, footprint, true); + stamp->image2, stamp->variance, polyValues, footprint, true); break; case PM_SUBTRACTION_MODE_2: status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2, - stamp->weight, polyValues, footprint, true); + stamp->variance, polyValues, footprint, true); status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2, - stamp->image1, stamp->weight, polyValues, footprint, true); + stamp->image1, stamp->variance, polyValues, footprint, true); break; case PM_SUBTRACTION_MODE_DUAL: @@ -581,14 +581,14 @@ #endif status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, - stamp->weight, polyValues, footprint, true); + stamp->variance, polyValues, footprint, true); status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL, - stamp->weight, polyValues, footprint, false); + stamp->variance, polyValues, footprint, false); status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1, - stamp->convolutions2, stamp->image1, stamp->weight, polyValues, + stamp->convolutions2, stamp->image1, stamp->variance, polyValues, footprint); status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, - stamp->image2, stamp->weight, polyValues, footprint, true); + stamp->image2, stamp->variance, polyValues, footprint, true); status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL, - stamp->image2, stamp->weight, polyValues, footprint, false); + stamp->image2, stamp->variance, polyValues, footprint, false); break; default: @@ -1033,5 +1033,5 @@ // Calculate residuals - psKernel *weight = stamp->weight; // Weight postage stamp + psKernel *variance = stamp->variance; // Variance postage stamp psImageInit(residual->image, 0.0); if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { @@ -1098,5 +1098,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - double dev = PS_SQR(residual->kernel[y][x]) / weight->kernel[y][x]; + double dev = PS_SQR(residual->kernel[y][x]) / variance->kernel[y][x]; deviation += dev; #ifdef TESTING @@ -1141,10 +1141,10 @@ psFitsClose(fits); } - if (stamp->weight) { + if (stamp->variance) { psString filename = NULL; - psStringAppend(&filename, "stamp_weight_%03d.fits", i); + psStringAppend(&filename, "stamp_variance_%03d.fits", i); psFits *fits = psFitsOpen(filename, "w"); psFree(filename); - psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL); + psFitsWriteImage(fits, NULL, stamp->variance->image, 0, NULL); psFitsClose(fits); } Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 21352) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 21363) @@ -60,5 +60,5 @@ const pmReadout *ro2, // Readout 2 const psImage *subMask, // Mask for subtraction, or NULL - psImage *weight, // Weight map + psImage *variance, // Variance map const psRegion *region, // Region of interest, or NULL float threshold, // Threshold for stamp finding @@ -80,5 +80,5 @@ psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); - if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, size)) { + if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, variance, size)) { psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps."); return false; @@ -111,7 +111,7 @@ conv1->mask = NULL; } - if (conv1->weight) { - psFree(conv1->weight); - conv1->weight = NULL; + if (conv1->variance) { + psFree(conv1->variance); + conv1->variance = NULL; } } @@ -126,7 +126,7 @@ conv2->mask = NULL; } - if (conv2->weight) { - psFree(conv2->weight); - conv2->weight = NULL; + if (conv2->variance) { + psFree(conv2->variance); + conv2->variance = NULL; } } @@ -191,15 +191,15 @@ } - // Where does our weight map come from? - // Getting the weight exactly right is not necessary --- it's just used for weighting. - psImage *weight = NULL; // Weight image to use - if (ro1->weight && ro2->weight) { - weight = (psImage*)psBinaryOp(NULL, ro1->weight, "+", ro2->weight); - } else if (ro1->weight) { - weight = psMemIncrRefCounter(ro1->weight); - } else if (ro2->weight) { - weight = psMemIncrRefCounter(ro2->weight); + // Where does our variance map come from? + // Getting the variance exactly right is not necessary --- it's just used for weighting. + psImage *variance = NULL; // Variance image to use + if (ro1->variance && ro2->variance) { + variance = (psImage*)psBinaryOp(NULL, ro1->variance, "+", ro2->variance); + } else if (ro1->variance) { + variance = psMemIncrRefCounter(ro1->variance); + } else if (ro2->variance) { + variance = psMemIncrRefCounter(ro2->variance); } else { - weight = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image); + variance = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image); } @@ -274,5 +274,5 @@ // We get the stamps here; we will also attempt to get stamps at the first iteration, but it // doesn't matter. - if (!getStamps(&stamps, ro1, ro2, subMask, weight, NULL, threshold, stampSpacing, + if (!getStamps(&stamps, ro1, ro2, subMask, variance, NULL, threshold, stampSpacing, size, footprint, subMode)) { goto MATCH_ERROR; @@ -336,5 +336,5 @@ psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); - if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing, + if (!getStamps(&stamps, ro1, ro2, subMask, variance, region, threshold, stampSpacing, size, footprint, subMode)) { goto MATCH_ERROR; @@ -439,8 +439,8 @@ psFree(subMask); subMask = NULL; - psFree(weight); - weight = NULL; - - if (!pmSubtractionBorder(conv1->image, conv1->weight, conv1->mask, size, maskBad)) { + psFree(variance); + variance = NULL; + + if (!pmSubtractionBorder(conv1->image, conv1->variance, conv1->mask, size, maskBad)) { psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image."); goto MATCH_ERROR; @@ -482,5 +482,5 @@ psFree(kernels); psFree(stamps); - psFree(weight); + psFree(variance); psFree(rng); return false; Index: trunk/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 21352) +++ trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 21363) @@ -71,5 +71,5 @@ double *sumII, // Sum of I(x)^2/sigma(x)^2 double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2 - const pmSubtractionStamp *stamp, // Stamp with weight + const pmSubtractionStamp *stamp, // Stamp with variance const psKernel *target, // Target stamp int kernelIndex, // Index for kernel component @@ -78,10 +78,10 @@ ) { - psKernel *weight = stamp->weight; // Weight, sigma(x)^2 + psKernel *variance = stamp->variance; // Variance, sigma(x)^2 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { psF32 *in = &target->kernel[y][-footprint]; // Dereference input - psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { @@ -98,5 +98,5 @@ static void accumulateConvolutions(double *sumC, // Sum of conv(x)/sigma(x)^2 double *sumCC, // Sum of conv(x)^2/sigma(x)^2 - const pmSubtractionStamp *stamp, // Stamp with input and weight + const pmSubtractionStamp *stamp, // Stamp with input and variance int kernelIndex, // Index for kernel component int footprint, // Size of region of interest @@ -104,9 +104,9 @@ ) { - psKernel *weight = stamp->weight; // Weight, sigma(x)^2 + psKernel *variance = stamp->variance; // Variance, sigma(x)^2 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { - psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution for (int x = -footprint; x <= footprint; x++, wt++, conv++) { @@ -120,5 +120,5 @@ static double accumulateChi2(const psKernel *target, // Target stamp - pmSubtractionStamp *stamp, // Stamp with weight + pmSubtractionStamp *stamp, // Stamp with variance int kernelIndex, // Index for kernel component double coeff, // Coefficient of convolution @@ -129,10 +129,10 @@ { double chi2 = 0.0; - psKernel *weight = stamp->weight; // Weight, sigma(x)^2 + psKernel *variance = stamp->variance; // Variance, sigma(x)^2 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { psF32 *in = &target->kernel[y][-footprint]; // Dereference input - psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { @@ -146,10 +146,10 @@ // Return the initial value of chi^2 static double initialChi2(const psKernel *target, // Target stamp - const pmSubtractionStamp *stamp, // Stamp with weight + const pmSubtractionStamp *stamp, // Stamp with variance int footprint, // Size of convolution pmSubtractionMode mode // Mode of subtraction ) { - psKernel *weight = stamp->weight; // Weight map + psKernel *variance = stamp->variance; // Variance map psKernel *source; // Source stamp switch (mode) { @@ -167,5 +167,5 @@ for (int y = -footprint; y <= footprint; y++) { psF32 *in = &target->kernel[y][-footprint]; // Dereference input - psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance psF32 *ref = &source->kernel[y][-footprint]; // Derference reference for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) { @@ -180,5 +180,5 @@ // Subtract a convolution from the input static void subtractConvolution(psKernel *target, // Target stamp - const pmSubtractionStamp *stamp, // Stamp with weight + const pmSubtractionStamp *stamp, // Stamp with variance int kernelIndex, // Index for kernel component float coeff, // Coefficient of subtraction @@ -288,7 +288,7 @@ // This sum is invariant to the kernel - psKernel *weight = stamp->weight; // Weight map for stamp + psKernel *variance = stamp->variance; // Variance map for stamp for (int v = -footprint; v <= footprint; v++) { - psF32 *wt = &weight->kernel[v][-footprint]; // Dereference weight map + psF32 *wt = &variance->kernel[v][-footprint]; // Dereference variance map for (int u = -footprint; u <= footprint; u++, wt++) { sum1 += 1.0 / *wt; Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 21352) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 21363) @@ -54,5 +54,5 @@ psFree(stamp->image1); psFree(stamp->image2); - psFree(stamp->weight); + psFree(stamp->variance); psFree(stamp->convolutions1); psFree(stamp->convolutions2); @@ -211,5 +211,5 @@ stamp->image1 = NULL; stamp->image2 = NULL; - stamp->weight = NULL; + stamp->variance = NULL; stamp->convolutions1 = NULL; stamp->convolutions2 = NULL; @@ -333,8 +333,8 @@ psFree(stamp->image1); psFree(stamp->image2); - psFree(stamp->weight); + psFree(stamp->variance); psFree(stamp->convolutions1); psFree(stamp->convolutions2); - stamp->image1 = stamp->image2 = stamp->weight = NULL; + stamp->image1 = stamp->image2 = stamp->variance = NULL; stamp->convolutions1 = stamp->convolutions2 = NULL; @@ -480,5 +480,5 @@ bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2, - psImage *weight, int kernelSize) + psImage *variance, int kernelSize) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -490,7 +490,7 @@ PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false); } - PS_ASSERT_IMAGE_NON_NULL(weight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image1, false); - PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(variance, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image1, false); + PS_ASSERT_IMAGE_TYPE(variance, PS_TYPE_F32, false); PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); @@ -520,5 +520,5 @@ assert(stamp->image1 == NULL); assert(stamp->image2 == NULL); - assert(stamp->weight == NULL); + assert(stamp->variance == NULL); psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest @@ -534,6 +534,6 @@ } - psImage *wtSub = psImageSubset(weight, region); // Subimage with stamp - stamp->weight = psKernelAllocFromImage(wtSub, size, size); + psImage *wtSub = psImageSubset(variance, region); // Subimage with stamp + stamp->variance = psKernelAllocFromImage(wtSub, size, size); psFree(wtSub); // Drop reference Index: trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 21352) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 21363) @@ -60,5 +60,5 @@ psKernel *image1; ///< Reference image postage stamp psKernel *image2; ///< Input image postage stamp - psKernel *weight; ///< Weight image postage stamp, or NULL + psKernel *variance; ///< Variance image postage stamp, or NULL psArray *convolutions1; ///< Convolutions of image 1 for each kernel component, or NULL psArray *convolutions2; ///< Convolutions of image 2 for each kernel component, or NULL @@ -120,5 +120,5 @@ psImage *image1, ///< Reference image psImage *image2, ///< Input image (or NULL) - psImage *weight, ///< Weight (variance) map + psImage *variance, ///< Variance map int kernelSize ///< Kernel half-size ); Index: trunk/psModules/src/imcombine/pmSubtractionThreads.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionThreads.c (revision 21352) +++ trunk/psModules/src/imcombine/pmSubtractionThreads.c (revision 21363) @@ -27,6 +27,6 @@ // XXX psMutexInit(ro->image); // XXX } - // XXX if (ro->weight) { - // XXX psMutexInit(ro->weight); + // XXX if (ro->variance) { + // XXX psMutexInit(ro->variance); // XXX } @@ -43,6 +43,6 @@ // XXX psMutexDestroy(ro->image); // XXX } - // XXX if (ro->weight) { - // XXX psMutexDestroy(ro->weight); + // XXX if (ro->variance) { + // XXX psMutexDestroy(ro->variance); // XXX } Index: trunk/psModules/src/objects/pmResiduals.c =================================================================== --- trunk/psModules/src/objects/pmResiduals.c (revision 21352) +++ trunk/psModules/src/objects/pmResiduals.c (revision 21363) @@ -4,6 +4,6 @@ * * @author EAM, IfA - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * Copyright 2004 IfA, University of Hawaii */ @@ -23,5 +23,5 @@ psFree (resid->Rx); psFree (resid->Ry); - psFree (resid->weight); + psFree (resid->variance); psFree (resid->mask); return; @@ -42,5 +42,5 @@ resid->Rx = psImageAlloc (nX, nY, PS_TYPE_F32); resid->Ry = psImageAlloc (nX, nY, PS_TYPE_F32); - resid->weight = psImageAlloc (nX, nY, PS_TYPE_F32); + resid->variance = psImageAlloc (nX, nY, PS_TYPE_F32); resid->mask = psImageAlloc (nX, nY, PM_TYPE_RESID_MASK); @@ -49,5 +49,5 @@ resid->xBin = xBin; resid->yBin = yBin; - resid->xCenter = 0.5*(nX - 1); + resid->xCenter = 0.5*(nX - 1); resid->yCenter = 0.5*(nY - 1); return resid; Index: trunk/psModules/src/objects/pmResiduals.h =================================================================== --- trunk/psModules/src/objects/pmResiduals.h (revision 21352) +++ trunk/psModules/src/objects/pmResiduals.h (revision 21363) @@ -4,6 +4,6 @@ * * @author EAM, IfA - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * Copyright 2004 IfA, University of Hawaii */ @@ -14,5 +14,5 @@ /// @{ -/** residual tables for sources +/** residual tables for sources */ typedef struct { @@ -20,5 +20,5 @@ psImage *Rx; psImage *Ry; - psImage *weight; + psImage *variance; psImage *mask; int xBin; @@ -31,5 +31,5 @@ bool psMemCheckResiduals(psPtr ptr); -// macros to abstract the resid mask type : these values must be consistent +// macros to abstract the resid mask type : these values must be consistent #define PM_TYPE_RESID_MASK PS_TYPE_U8 /**< the psElemType to use for mask image */ #define PM_TYPE_RESID_MASK_DATA U8 /**< the data member to use for mask image */ Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 21352) +++ trunk/psModules/src/objects/pmSource.c (revision 21363) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.68 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-06 01:09:50 $ + * @version $Revision: 1.69 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -44,5 +44,5 @@ psFree(tmp->peak); psFree(tmp->pixels); - psFree(tmp->weight); + psFree(tmp->variance); psFree(tmp->maskObj); psFree(tmp->maskView); @@ -66,5 +66,5 @@ psFree (source->pixels); - psFree (source->weight); + psFree (source->variance); psFree (source->maskObj); psFree (source->maskView); @@ -73,5 +73,5 @@ source->pixels = NULL; - source->weight = NULL; + source->variance = NULL; source->maskObj = NULL; source->maskView = NULL; @@ -101,5 +101,5 @@ source->peak = NULL; source->pixels = NULL; - source->weight = NULL; + source->variance = NULL; source->maskObj = NULL; source->maskView = NULL; @@ -144,5 +144,5 @@ } // this copy is used to allow multiple fit attempts on the - // same object. the pixels, weight, and mask arrays all point to + // same object. the pixels, variance, and mask arrays all point to // the same original subarrays. the peak and moments point at // the original values. @@ -167,5 +167,5 @@ // We want a new view, but pointing at the same pixels. source->pixels = psImageCopyView(NULL, in->pixels); - source->weight = psImageCopyView(NULL, in->weight); + source->variance = psImageCopyView(NULL, in->variance); source->maskView = in->maskView ? psImageCopyView(NULL, in->maskView) : NULL; @@ -199,6 +199,6 @@ // these images are subset images of the equivalent parents mySource->pixels = psImageSubset(readout->image, srcRegion); - if (readout->weight) { - mySource->weight = psImageSubset(readout->weight, srcRegion); + if (readout->variance) { + mySource->variance = psImageSubset(readout->variance, srcRegion); } if (readout->mask) { @@ -237,5 +237,5 @@ extend |= (mySource->pixels == NULL); - extend |= (mySource->weight == NULL); + extend |= (mySource->variance == NULL); extend |= (mySource->maskObj == NULL); extend |= (mySource->maskView == NULL); @@ -245,9 +245,9 @@ // re-create the subimage psFree (mySource->pixels); - psFree (mySource->weight); + psFree (mySource->variance); psFree (mySource->maskView); mySource->pixels = psImageSubset(readout->image, newRegion); - mySource->weight = psImageSubset(readout->weight, newRegion); + mySource->variance = psImageSubset(readout->variance, newRegion); mySource->maskView = psImageSubset(readout->mask, newRegion); mySource->region = newRegion; @@ -670,5 +670,5 @@ pmSource->peak pmSource->pixels - pmSource->weight + pmSource->variance pmSource->mask @@ -737,5 +737,5 @@ psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->weight->data.F32[row]; + psF32 *vWgt = source->variance->data.F32[row]; psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; @@ -924,5 +924,5 @@ if (mode & PM_MODEL_OP_NOISE) { // XXX need to scale by the gain... - target = source->weight->data.F32; + target = source->variance->data.F32; } @@ -946,5 +946,5 @@ psImage *target = source->pixels; if (mode & PM_MODEL_OP_NOISE) { - target = source->weight; + target = source->variance; } Index: trunk/psModules/src/objects/pmSource.h =================================================================== --- trunk/psModules/src/objects/pmSource.h (revision 21352) +++ trunk/psModules/src/objects/pmSource.h (revision 21363) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.27 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.28 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -72,5 +72,5 @@ pmPeak *peak; ///< Description of peak pixel. psImage *pixels; ///< Rectangular region including object pixels. - psImage *weight; ///< Image variance. + psImage *variance; ///< Image variance. psImage *maskObj; ///< unique mask for this object which marks included pixels associated with objects. psImage *maskView; ///< view into global image mask for this object region @@ -210,5 +210,5 @@ psMetadata *metadata, ///< Contains classification parameters pmPSFClump clump, ///< Statistics about the PSF clump - psImageMaskType maskSat ///< Mask value for saturated pixels + psImageMaskType maskSat ///< Mask value for saturated pixels ); Index: trunk/psModules/src/objects/pmSourceExtendedPars.c =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 21352) +++ trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 21363) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-12-08 02:51:14 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -71,5 +71,5 @@ psFree(profile->radius); psFree(profile->flux); - psFree(profile->weight); + psFree(profile->variance); return; } @@ -82,5 +82,5 @@ profile->radius = NULL; profile->flux = NULL; - profile->weight = NULL; + profile->variance = NULL; return profile; Index: trunk/psModules/src/objects/pmSourceExtendedPars.h =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 21352) +++ trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 21363) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-01-02 20:39:04 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -17,5 +17,5 @@ psVector *radius; psVector *flux; - psVector *weight; + psVector *variance; } pmSourceRadialProfile; Index: trunk/psModules/src/objects/pmSourceFitModel.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.c (revision 21352) +++ trunk/psModules/src/objects/pmSourceFitModel.c (revision 21363) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.30 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -62,5 +62,5 @@ PS_ASSERT_PTR_NON_NULL(source->pixels, false); PS_ASSERT_PTR_NON_NULL(source->maskObj, false); - PS_ASSERT_PTR_NON_NULL(source->weight, false); + PS_ASSERT_PTR_NON_NULL(source->variance, false); psBool fitStatus = true; @@ -84,6 +84,6 @@ continue; } - // skip zero-weight points - if (source->weight->data.F32[i][j] == 0) { + // skip zero-variance points + if (source->variance->data.F32[i][j] == 0) { continue; } @@ -102,8 +102,8 @@ // psMinimizeLMChi2 takes wt = 1/dY^2. suggestion from RHL is to use the local sky - // as weight to avoid the bias from systematic errors here we would just use the + // as variance to avoid the bias from systematic errors here we would just use the // source sky variance if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) { - yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; + yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j]; } else { yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT; Index: trunk/psModules/src/objects/pmSourceFitSet.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitSet.c (revision 21352) +++ trunk/psModules/src/objects/pmSourceFitSet.c (revision 21363) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii * @@ -55,12 +55,12 @@ if (!fitSets) { - fitSets = psArrayAlloc (PS_MAX (1, nThreads)); - } - - // the allocated elements should be NULL on psArrayAlloc, - // and a previously allocated array of fitSets should have been cleared + fitSets = psArrayAlloc (PS_MAX (1, nThreads)); + } + + // the allocated elements should be NULL on psArrayAlloc, + // and a previously allocated array of fitSets should have been cleared // before pmSourceFitSetInit is called for (int i = 0; i < fitSets->n; i++) { - psAssert (fitSets->data[i] == NULL, "failure to init or clear fitSets?"); + psAssert (fitSets->data[i] == NULL, "failure to init or clear fitSets?"); } return true; @@ -125,10 +125,10 @@ // we do not need to lock to do this.... for (int i = 0; i < fitSets->n; i++) { - thisSet = fitSets->data[i]; - if (!thisSet) continue; - if (thisSet->thread == id) { - break; - } - thisSet = NULL; + thisSet = fitSets->data[i]; + if (!thisSet) continue; + if (thisSet->thread == id) { + break; + } + thisSet = NULL; } psAssert (thisSet == NULL, "pmSourceFitSetDataSet() called but previous entry not cleared"); @@ -137,11 +137,11 @@ pthread_mutex_lock (&fitSetInitMutex); - + // find an entry that is NULL and place it there for (int i = 0; i < fitSets->n; i++) { - if (fitSets->data[i]) continue; - fitSets->data[i] = thisSet; - pthread_mutex_unlock (&fitSetInitMutex); - return thisSet; + if (fitSets->data[i]) continue; + fitSets->data[i] = thisSet; + pthread_mutex_unlock (&fitSetInitMutex); + return thisSet; } pthread_mutex_unlock (&fitSetInitMutex); @@ -160,10 +160,10 @@ pmSourceFitSetData *thisSet = NULL; for (int i = 0; i < fitSets->n; i++) { - thisSet = fitSets->data[i]; - if (!thisSet) continue; - if (thisSet->thread == id) { - break; - } - thisSet = NULL; + thisSet = fitSets->data[i]; + if (!thisSet) continue; + if (thisSet->thread == id) { + break; + } + thisSet = NULL; } psAssert (thisSet != NULL, "pmSourceFitSetDataGet() called, but no entry found"); @@ -184,10 +184,10 @@ pmSourceFitSetData *thisSet = NULL; for (i = 0; i < fitSets->n; i++) { - thisSet = fitSets->data[i]; - if (!thisSet) continue; - if (thisSet->thread == id) { - break; - } - thisSet = NULL; + thisSet = fitSets->data[i]; + if (!thisSet) continue; + if (thisSet->thread == id) { + break; + } + thisSet = NULL; } psAssert (thisSet != NULL, "pmSourceFitSetDataClear() called, but no entry found"); @@ -256,14 +256,14 @@ psVector *derivOne = set->derivSet->data[i]; - // one or the other (param or deriv) must be set - assert ((deriv != NULL) || (param != NULL)); - - // if we are setting derive, derivOne and paramOne must be same length + // one or the other (param or deriv) must be set + assert ((deriv != NULL) || (param != NULL)); + + // if we are setting derive, derivOne and paramOne must be same length assert ((deriv == NULL) || (paramOne->n == derivOne->n)); - + for (int j = 0; j < paramOne->n; j++, n++) { - if (param) { - param->data.F32[n] = paramOne->data.F32[j]; - } + if (param) { + param->data.F32[n] = paramOne->data.F32[j]; + } if (deriv) { deriv->data.F32[n] = derivOne->data.F32[j]; @@ -275,5 +275,5 @@ // distribute parameters from single param and deriv vectors into FitSet models -bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param) +bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param) { PS_ASSERT_VECTOR_NON_NULL(param, false); @@ -301,5 +301,5 @@ // set the model parameters for this fit set -bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, +bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, const psVector *param, pmSource *source, psMinimization *myMin, int nPix, bool fitStatus) @@ -420,9 +420,9 @@ for (int i = 0; i < modelSet->n; i++) { - pmModel *model = modelSet->data[i]; - int nParams = pmModelClassParameterCount (model->type); - for (int j = 0; j < nParams; j++) { - fprintf (file, "%d %d : %f %f\n", i, j, model->params->data.F32[j], model->dparams->data.F32[j]); - } + pmModel *model = modelSet->data[i]; + int nParams = pmModelClassParameterCount (model->type); + for (int j = 0; j < nParams; j++) { + fprintf (file, "%d %d : %f %f\n", i, j, model->params->data.F32[j], model->dparams->data.F32[j]); + } } return true; @@ -435,5 +435,5 @@ psVector *derivOne = set->derivSet->data[i]; for (int j = 0; j < paramOne->n; j++) { - fprintf (file, "%d %d : %f %f\n", i, j, paramOne->data.F32[j], derivOne->data.F32[j]); + fprintf (file, "%d %d : %f %f\n", i, j, paramOne->data.F32[j], derivOne->data.F32[j]); } } @@ -450,5 +450,5 @@ PS_ASSERT_PTR_NON_NULL(source->pixels, false); PS_ASSERT_PTR_NON_NULL(source->maskObj, false); - PS_ASSERT_PTR_NON_NULL(source->weight, false); + PS_ASSERT_PTR_NON_NULL(source->variance, false); bool fitStatus = true; @@ -471,6 +471,6 @@ continue; } - // skip zero-weight points - if (source->weight->data.F32[i][j] == 0) { + // skip zero-variance points + if (source->variance->data.F32[i][j] == 0) { continue; } @@ -489,8 +489,8 @@ // psMinimizeLMChi2 takes wt = 1/dY^2. suggestion from RHL is to use the local sky - // as weight to avoid the bias from systematic errors here we would just use the + // as variance to avoid the bias from systematic errors here we would just use the // source sky variance if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) { - yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; + yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j]; } else { yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT; @@ -543,5 +543,5 @@ psFree (params); psFree(constraint); - pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets + pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets return(false); } Index: trunk/psModules/src/objects/pmSourceMoments.c =================================================================== --- trunk/psModules/src/objects/pmSourceMoments.c (revision 21352) +++ trunk/psModules/src/objects/pmSourceMoments.c (revision 21363) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -47,5 +47,5 @@ pmSource->peak pmSource->pixels - pmSource->weight (optional) + pmSource->variance (optional) pmSource->mask (optional) @@ -99,5 +99,5 @@ psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->weight->data.F32[row]; + psF32 *vWgt = source->variance->data.F32[row]; psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; @@ -110,5 +110,5 @@ vMsk++; } - if (isnan(*vPix)) continue; + if (isnan(*vPix)) continue; psF32 xDiff = col - xPeak; @@ -189,5 +189,5 @@ psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->weight->data.F32[row]; + psF32 *vWgt = source->variance->data.F32[row]; psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; @@ -200,5 +200,5 @@ vMsk++; } - if (isnan(*vPix)) continue; + if (isnan(*vPix)) continue; psF32 xDiff = col - xCM; @@ -206,6 +206,6 @@ // radius is just a function of (xDiff, yDiff) - psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); - psF32 r = sqrt(r2); + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + psF32 r = sqrt(r2); if (r > radius) continue; @@ -278,5 +278,5 @@ psTrace ("psModules.objects", 4, "Mxx: %f Mxy: %f Myy: %f Mxxx: %f Mxxy: %f Mxyy: %f Myyy: %f Mxxxx: %f Mxxxy: %f Mxxyy: %f Mxyyy: %f Mxyyy: %f\n", - source->moments->Mxx, source->moments->Mxy, source->moments->Myy, + source->moments->Mxx, source->moments->Mxy, source->moments->Myy, source->moments->Mxxx, source->moments->Mxxy, source->moments->Mxyy, source->moments->Myyy, source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy); @@ -286,8 +286,8 @@ if (source->moments->Mxx < 0) { - fprintf (stderr, "error: neg second moment??\n"); + fprintf (stderr, "error: neg second moment??\n"); } if (source->moments->Myy < 0) { - fprintf (stderr, "error: neg second moment??\n"); + fprintf (stderr, "error: neg second moment??\n"); } @@ -341,5 +341,5 @@ psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->weight->data.F32[row]; + psF32 *vWgt = source->variance->data.F32[row]; psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; @@ -352,5 +352,5 @@ vMsk++; } - if (isnan(*vPix)) continue; + if (isnan(*vPix)) continue; psF32 xDiff = col - xPeak; @@ -407,5 +407,5 @@ } -# if (PS_TRACE_ON) +# if (PS_TRACE_ON) float Sxx = PS_MAX(X2/Sum - PS_SQR(x), 0); float Sxy = XY / Sum; Index: trunk/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.c (revision 21352) +++ trunk/psModules/src/objects/pmSourcePhotometry.c (revision 21363) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -446,5 +446,5 @@ assert (Pj != NULL); - const psImage *Wi = Mi->weight; + const psImage *Wi = Mi->variance; if (!unweighted_sum) { assert (Wi != NULL); @@ -511,5 +511,5 @@ assert (Pj != NULL); - const psImage *Wi = Mi->weight; + const psImage *Wi = Mi->variance; if (!unweighted_sum) { assert (Wi != NULL); @@ -567,5 +567,5 @@ const psImage *Pi = Mi->pixels; assert (Pi != NULL); - const psImage *Wi = Mi->weight; + const psImage *Wi = Mi->variance; if (!unweighted_sum) { assert (Wi != NULL); @@ -612,5 +612,5 @@ # endif -bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, +bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal) { @@ -618,5 +618,5 @@ PS_ASSERT_PTR_NON_NULL(image, false); PS_ASSERT_PTR_NON_NULL(mask, false); - PS_ASSERT_PTR_NON_NULL(weight, false); + PS_ASSERT_PTR_NON_NULL(variance, false); double dC = 0.0; @@ -626,7 +626,7 @@ if (mask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] & maskVal) continue; - if (weight->data.F32[j][i] <= 0) - continue; - dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i]; + if (variance->data.F32[j][i] <= 0) + continue; + dC += PS_SQR (image->data.F32[j][i]) / variance->data.F32[j][i]; Npix ++; } @@ -648,5 +648,5 @@ const psImage *Pi = Mi->modelFlux; assert (Pi != NULL); - const psImage *Wi = Mi->weight; + const psImage *Wi = Mi->variance; if (!unweighted_sum) { assert (Wi != NULL); @@ -706,5 +706,5 @@ assert (Pj != NULL); - const psImage *Wi = Mi->weight; + const psImage *Wi = Mi->variance; if (!unweighted_sum) { assert (Wi != NULL); @@ -770,5 +770,5 @@ assert (Pj != NULL); - const psImage *Wi = Mi->weight; + const psImage *Wi = Mi->variance; if (!unweighted_sum) { assert (Wi != NULL); Index: trunk/psModules/src/objects/pmSourceSky.c =================================================================== --- trunk/psModules/src/objects/pmSourceSky.c (revision 21352) +++ trunk/psModules/src/objects/pmSourceSky.c (revision 21363) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -113,5 +113,5 @@ psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__); PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_IMAGE_NON_NULL(source->weight, false); + PS_ASSERT_IMAGE_NON_NULL(source->variance, false); PS_ASSERT_IMAGE_NON_NULL(source->maskObj, false); PS_ASSERT_PTR_NON_NULL(source->peak, false); @@ -127,5 +127,5 @@ } - psImage *image = source->weight; + psImage *image = source->variance; psImage *mask = source->maskObj; pmPeak *peak = source->peak;