Changeset 20486 for trunk/psModules/src/camera/pmFPAMaskWeight.c

Timestamp:

Oct 31, 2008, 11:50:16 AM (18 years ago)

Author:

Paul Price

Message:

Adding functions to renormalise the weight maps based on pixels or on (weighted aperture) photometry.

File:

• trunk/psModules/src/camera/pmFPAMaskWeight.c (modified) (2 diffs)

trunk/psModules/src/camera/pmFPAMaskWeight.c

@@ -14 +14 @@
 
 #define PIXELS_BUFFER 100               // Size of buffer for allocating pixel lists
+#define RENORM_NUM_SIGMA 3.0            // Number of standard deviations for Gaussian kernel
+#define RENORM_PEAK 4.0                 // Number of standard deviations of noise for fake source peak flux
+
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -315 +318 @@
 
     return success;
+}
+
+
+bool pmReadoutWeightRenormPixels(const pmReadout *readout, psMaskType 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
+
+    if (!psMemIncrRefCounter(rng)) {
+        rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
+    }
+
+    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);
+        psFree(rng);
+        return false;
+    }
+    float meanVariance = weightStats->robustMedian; // Mean variance
+    psFree(weightStats);
+
+    psStats *imageStats = psStatsAlloc(stdevStat);// Statistics for mean
+    if (!psImageBackground(imageStats, NULL, image, mask, maskVal, rng)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure stdev of image");
+        psFree(imageStats);
+        psFree(rng);
+        return false;
+    }
+    float stdevImage = imageStats->robustStdev; // Standard deviation of image
+    psFree(imageStats);
+    psFree(rng);
+
+    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));
+
+    return true;
+}
+
+
+bool pmReadoutWeightRenormPhot(const pmReadout *readout, psMaskType 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);
+
+    if (!psMemIncrRefCounter(rng)) {
+        rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
+    }
+
+    psImage *image = readout->image, *mask = readout->mask, *weight = readout->weight; // Readout images
+
+    // Measure background
+    psStats *bgStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);// Statistics for background
+    if (!psImageBackground(bgStats, NULL, image, mask, maskVal, rng)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image");
+        psFree(bgStats);
+        psFree(rng);
+        return false;
+    }
+    float bgMean = bgStats->robustMedian; // Background level
+    float bgNoise = bgStats->robustStdev; // Background standard deviation
+    psFree(bgStats);
+    psTrace("psModules.camera", 5, "Background is %f +/- %f\n", bgMean, bgNoise);
+
+
+    // Construct kernels for flux measurement
+    // We use N(0,width) and N(0,width/sqrt(2)) kernels, following psphotSignificanceImage.
+    float sigFactor = 4.0 * M_PI * PS_SQR(width); // Factor for conversion from im/wt ratio to significance
+    int size = RENORM_NUM_SIGMA, fullSize = 2 * size + 1; // Half-size and full size of Gaussian
+    psVector *gauss = psVectorAlloc(fullSize, PS_TYPE_F32); // Gaussian for weighting
+    psVector *gauss2 = psVectorAlloc(fullSize, PS_TYPE_F32); // Gaussian squared
+    for (int i = 0, x = -size; i < fullSize; i++, x++) {
+        gauss->data.F32[i] = expf(-0.5 * PS_SQR(x) / PS_SQR(width));
+        gauss2->data.F32[i] = expf(-PS_SQR(x) / PS_SQR(width));
+    }
+
+    // Size of image
+    int numCols = image->numCols, numRows = image->numRows; // Size of images
+    int xSize = numCols - fullSize, ySize = numRows - fullSize; // Size of consideration
+    int xOffset = size, yOffset = size;       // Offset to region of consideration
+
+    // Measure fluxes
+    float peakFlux = RENORM_PEAK * bgNoise;     // Peak flux for fake sources
+    psVector *noise = psVectorAlloc(num, PS_TYPE_F32); // Measurements of the noise
+    psVector *source = psVectorAlloc(num, PS_TYPE_F32); // Measurements of fake sources
+    psVector *guess = psVectorAlloc(num, PS_TYPE_F32); // Guess at significance
+    psVector *photMask = psVectorAlloc(num, PS_TYPE_MASK); // Mask for fluxes
+    for (int i = 0; i < num; i++) {
+        // Coordinates of interest
+        int xPix = psRandomUniform(rng) * xSize + xOffset;
+        int yPix = psRandomUniform(rng) * ySize + yOffset;
+        psAssert(xPix - size >= 0 && xPix + size < numCols &&
+                 yPix - size >= 0 && yPix + size > numRows,
+                 "Bad pixel position");
+
+        // Weighted aperture photometry
+        // This has the same effect as smoothing the image by the window function
+        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 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 xSumGauss = 0.0, xSumGauss2 = 0.0; // Sums of Gaussian kernels in x
+            float yGauss = peakFlux * gauss->data.F32[v]; // Value of Gaussian in y
+            for (int u = 0, x = xPix - size; u < fullSize; u++, x++) {
+                if (mask && mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) {
+                    continue;
+                }
+                float value = image->data.F32[y][x] - bgMean; // Value of image
+                float xGauss = gauss->data.F32[u]; // Value of Gaussian in x
+                float xGauss2 = gauss2->data.F32[u]; // Value of Gaussian^2 in x
+                xSumNoise += value * xGauss;
+                xSumSource += (value + xGauss * yGauss) * xGauss;
+                xSumWeight += weight->data.F32[y][x] * xGauss2;
+                xSumGauss += xGauss;
+                xSumGauss2 += xGauss2;
+            }
+            float yGauss2 = gauss2->data.F32[v]; // Value of Gaussian^2 in y
+            sumNoise += xSumNoise * yGauss;
+            sumSource += xSumSource * yGauss;
+            sumWeight += xSumWeight * yGauss2;
+            sumGauss += xSumGauss * yGauss;
+            sumGauss2 += xSumGauss2 * yGauss2;
+        }
+
+        float weightValue = sumWeight / sumGauss2; // Value of convolved weight
+        photMask->data.PS_TYPE_MASK_DATA[i] = ((isfinite(sumNoise) && isfinite(sumSource) &&
+                                                isfinite(weightValue) && sumGauss > 0 && sumGauss2 > 0) ?
+                                               0 : 0xFF);
+
+        noise->data.F32[i] = sumNoise / sumGauss;
+        source->data.F32[i] = sumSource / sumGauss;
+        guess->data.F32[i] = (sumSource > 0) ? sigFactor * PS_SQR(sumSource) / sumWeight * sumGauss2 : 0.0;
+        psTrace("psModules.camera", 10, "Flux %d (%d,%d): %f, %f, %f\n",
+                i, xPix, yPix, source->data.F32[i], noise->data.F32[i], guess->data.F32[i]);
+    }
+    psFree(gauss);
+    psFree(gauss2);
+    psFree(rng);
+
+    // Standard deviation of fluxes gives us the real significance
+    psStats *stdevStats = psStatsAlloc(stdevStat); // Statistics
+    if (!psVectorStats(stdevStats, noise, NULL, photMask, 0xFF)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure standard deviation of fluxes");
+        psFree(stdevStats);
+        psFree(noise);
+        psFree(source);
+        psFree(guess);
+        psFree(photMask);
+        return false;
+    }
+    float stdev = psStatsGetValue(stdevStats, stdevStat); // Stadard deviation of fluxes
+    psFree(stdevStats);
+    psFree(noise);
+    psTrace("psModules.camera", 5, "Standard deviation of fluxes is %f\n", stdev);
+
+    // Ratio of measured significance to guessed significance
+    psVector *ratio = psVectorAlloc(num, PS_TYPE_F32); // Ratio of measured to guess
+    for (int i = 0; i < num; i++) {
+        float measuredSig = PS_SQR(source->data.F32[i] / stdev); // Measured significance
+        if (source->data.F32[i] <= 0.0) {
+            photMask->data.PS_TYPE_MASK_DATA[i] = 0xFF;
+        }
+        ratio->data.F32[i] = measuredSig / guess->data.F32[i];
+    }
+    psFree(source);
+    psFree(guess);
+
+    psStats *meanStats = psStatsAlloc(meanStat | stdevStat); // Statistics
+    if (!psVectorStats(meanStats, ratio, NULL, photMask, 0xFF)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure mean ratio");
+        psFree(meanStats);
+        psFree(ratio);
+        psFree(photMask);
+        return false;
+    }
+    float meanRatio = psStatsGetValue(meanStats, meanStat); // Mean ratio
+    psTrace("psModules.camera", 5, "Mean significance ratio is %f +/- %f\n",
+            meanRatio, psStatsGetValue(meanStats, stdevStat));
+    psFree(meanStats);
+    psFree(ratio);
+    psFree(photMask);
+
+    psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising weight map by %f", meanRatio);
+    psBinaryOp(weight, weight, "*", psScalarAlloc(meanRatio, PS_TYPE_F32));
+
+    return true;
 }