Index: trunk/psModules/src/camera/pmFPAMaskWeight.c =================================================================== --- trunk/psModules/src/camera/pmFPAMaskWeight.c (revision 21183) +++ 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);