Changeset 27839 for branches/simmosaic_branches/psModules/src/camera/pmFPAMaskWeight.c

Timestamp:

May 3, 2010, 8:45:22 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simmosaic_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/camera/pmFPAMaskWeight.c (modified) (10 diffs)

branches/simmosaic_branches/psModules/src/camera/pmFPAMaskWeight.c

@@ -111 +111 @@
         // psError(PS_ERR_IO, true, "CELL.SATURATION is not set --- unable to set mask.\n");
         // return false;
-        psWarning("CELL.SATURATION is not set --- completely masking cell.\n");
-        saturation = NAN;
+        psWarning("CELL.SATURATION is not set --- completely masking cell.\n");
+        saturation = NAN;
     }
     float bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); // Bad level
@@ -118 +118 @@
         // psError(PS_ERR_IO, true, "CELL.BAD is not set --- unable to set mask.\n");
         // return false;
-        psWarning("CELL.BAD is not set --- completely masking cell.\n");
-        bad = NAN;
+        psWarning("CELL.BAD is not set --- completely masking cell.\n");
+        bad = NAN;
     }
     psTrace("psModules.camera", 5, "Saturation: %f, bad: %f\n", saturation, bad);
 
-    // if CELL.GAIN or CELL.READNOISE are not set, then the variance will be set to NAN; 
+    // if CELL.GAIN or CELL.READNOISE are not set, then the variance will be set to NAN;
     // in this case, we have to set the mask as well
     float gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); // Cell gain
@@ -140 +140 @@
     // completely mask if SATURATION or BAD are invalid
     if (isnan(saturation) || isnan(bad) || isnan(gain) || isnan(readnoise)) {
-        psImageInit(mask, badMask);
-        return true;
+        psImageInit(mask, badMask);
+        return true;
     }
 
@@ -230 +230 @@
         // return false;
         psWarning("CELL.GAIN is not set --- setting variance to NAN\n");
-        gain = NAN;
+        gain = NAN;
     }
     float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise
@@ -237 +237 @@
         // return false;
         psWarning("CELL.READNOISE is not set --- setting variance to NAN\n");
-        readnoise = NAN;
+        readnoise = NAN;
     }
     // if we have a non-NAN readnoise, then we need to ensure it has been updated (not necessary if NAN)
@@ -248 +248 @@
     if (isnan(gain) || isnan(readnoise)) {
         if (!readout->variance) {
-            // generate the image if needed
+            // generate the image if needed
             readout->variance = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
         }
-        // XXX need to set the mask, if defined
+        // XXX need to set the mask, if defined
         psImageInit(readout->variance, NAN);
-        return true;
+        return true;
     }
 
@@ -262 +262 @@
 
         // a negative variance is non-sensical. if the image value drops below 1, the variance must be 1.
-        // XXX this calculation is wrong: limit is 1 e-, but this is in DN
+        // XXX this calculation is wrong: limit is 1 e-, but this is in DN
         readout->variance = (psImage*)psUnaryOp(readout->variance, readout->variance, "abs");
         readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "max",
@@ -276 +276 @@
     // apply a supplied readnoise map (NOTE: in DN, not electrons):
     if (noiseMap) {
-        psImage *rdVar = (psImage*)psBinaryOp(NULL, (const psPtr) noiseMap, "*", (const psPtr) noiseMap);
-        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", rdVar);
-        psFree (rdVar);
+        psImage *rdVar = (psImage*)psBinaryOp(NULL, (const psPtr) noiseMap, "*", (const psPtr) noiseMap);
+        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", rdVar);
+        psFree (rdVar);
     } else {
-        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32));
+        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32));
     }
 
@@ -362 +362 @@
 
 
-bool pmReadoutVarianceRenormPixels(const pmReadout *readout, psImageMaskType maskVal,
-                                 psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng)
+bool pmReadoutVarianceRenormalise(const pmReadout *readout, psImageMaskType maskVal,
+                                  int sample, float minValid, float maxValid)
 {
     PM_ASSERT_READOUT_NON_NULL(readout, false);
@@ -370 +370 @@
 
     psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout parts
-
-    if (!psMemIncrRefCounter(rng)) {
-        rng = psRandomAlloc(PS_RANDOM_TAUS);
-    }
-
-    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 = varianceStats->robustMedian; // Mean variance
-    psFree(varianceStats);
-
-    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 variance map by %f", correction);
-    psBinaryOp(variance, variance, "*", psScalarAlloc(correction, PS_TYPE_F32));
-
-    return true;
-}
-
-
-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_VARIANCE(readout, false);
-
-    if (!psMemIncrRefCounter(rng)) {
-        rng = psRandomAlloc(PS_RANDOM_TAUS);
-    }
-
-    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // 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_VECTOR_MASK); // Mask for fluxes
-    for (int i = 0; i < num; i++) {
-        // Coordinates of interest
-        int xPix = psRandomUniform(rng) * xSize + xOffset + 0.5;
-        int yPix = psRandomUniform(rng) * ySize + yOffset + 0.5;
-        psAssert(xPix - size >= 0 && xPix + size < numCols &&
-                 yPix - size >= 0 && yPix + size < numRows,
-                 "Bad pixel position: %d,%d", xPix, yPix);
-
-        // 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 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 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
-            for (int u = 0, x = xPix - size; u < fullSize; u++, x++) {
-                if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) {
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+
+    int xMin, xMax, yMin, yMax;         // Bounds of image
+    if (mask) {
+        xMin = numCols;
+        xMax = 0;
+        yMin = numRows;
+        yMax = 0;
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                if (mask->data.PS_TYPE_IMAGE_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 + peakFlux * xGauss * yGauss) * xGauss;
-                xSumVariance += variance->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;
-            sumVariance += xSumVariance * yGauss2;
-            sumGauss += xSumGauss * yGauss;
-            sumGauss2 += xSumGauss2 * yGauss2;
-        }
-
-        photMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = ((isfinite(sumNoise) && isfinite(sumSource) &&
-                                                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 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) / smoothVariance : 0.0;
-        psTrace("psModules.camera", 10, "Flux %d (%d,%d): %f, %f, %f\n",
-                i, xPix, yPix, smoothImageNoise, smoothImageSource, smoothVariance);
-    }
-    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);
+                xMin = PS_MIN(xMin, x);
+                xMax = PS_MAX(xMax, x);
+                yMin = PS_MIN(yMin, y);
+                yMax = PS_MAX(yMax, y);
+            }
+        }
+    } else {
+        xMin = 0;
+        xMax = numCols;
+        yMin = 0;
+        yMax = numRows;
+    }
+
+    int xNum = xMax - xMin, yNum = yMax - yMin; // Number of pixels
+
+    int numPix = xNum * yNum;                                  // Number of pixels
+    int num = PS_MIN(sample, numPix);                          // Number we care about
+    psVector *signoise = psVectorAllocEmpty(num, PS_TYPE_F32);   // Signal-to-noise values
+
+    if (num >= numPix) {
+        // We have an image smaller than Nsubset, so just loop over the image pixels
+        int index = 0;                  // Index for vector
+        for (int y = yMin; y < yMax; y++) {
+            for (int x = xMin; x < xMax; x++) {
+                if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
+                    !isfinite(image->data.F32[y][x]) || !isfinite(variance->data.F32[y][x])) {
+                    continue;
+                }
+
+                signoise->data.F32[index] = image->data.F32[y][x] / sqrtf(variance->data.F32[y][x]);
+                index++;
+            }
+        }
+        signoise->n = index;
+    } else {
+        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+        int index = 0;                  // Index for vector
+        for (long i = 0; i < num; i++) {
+            // Pixel coordinates
+            int pixel = numPix * psRandomUniform(rng);
+            int x = xMin + pixel % xNum;
+            int y = yMin + pixel / xNum;
+
+            if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
+                !isfinite(image->data.F32[y][x]) || !isfinite(variance->data.F32[y][x])) {
+                continue;
+            }
+
+            signoise->data.F32[index] = image->data.F32[y][x] / sqrtf(variance->data.F32[y][x]);
+            index++;
+        }
+        signoise->n = index;
+        psFree(rng);
+    }
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics
+
+    if (!psVectorStats(stats, signoise, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics on S/N image");
+        psFree(signoise);
         return false;
     }
-    float stdev = psStatsGetValue(stdevStats, stdevStat); // Standard 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
-        ratio->data.F32[i] = measuredSig / guess->data.F32[i];
-        if (guess->data.F32[i] <= 0.0 || source->data.F32[i] <= 0.0 || !isfinite(ratio->data.F32[i])) {
-            photMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xFF;
-        }
-        psTrace("psModules.camera", 9, "Ratio %d: %f, %f, %f\n",
-                i, guess->data.F32[i], measuredSig, ratio->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);
+    psFree(signoise);
+
+    float covar = sqrtf(psImageCovarianceFactor(readout->covariance)); // Covariance factor
+    float correction = stats->robustStdev / covar; // Correction factor
+    psFree(stats);
+    psLogMsg("psModules.camera", PS_LOG_DETAIL, "Variance renormalisation factor is %f", correction);
+
+    // Check valid range of correction factor
+    if ((isfinite(minValid) && correction < minValid) || (isfinite(maxValid) && correction > maxValid)) {
+        psError(PS_ERR_UNKNOWN, true, "Variance renormalisation is outside valid range: %f vs %f:%f --- no correction made", correction, minValid, maxValid);
+        psMetadataAddF32(readout->analysis, PS_LIST_TAIL, PM_READOUT_ANALYSIS_RENORM, 0, "Renormalisation of variance", PS_SQR(correction));
         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 variance map by %f", meanRatio);
-    psBinaryOp(variance, variance, "*", psScalarAlloc(meanRatio, PS_TYPE_F32));
-
-    return true;
-}
-
-
-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_VARIANCE(readout, false);
-    PS_ASSERT_INT_POSITIVE(width, false);
-
-    if (!psMemIncrRefCounter(rng)) {
-        rng = psRandomAlloc(PS_RANDOM_TAUS);
-    }
-
-    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
-    float xSize = numCols / (float)xNum, ySize = numRows / (float)yNum; // Size of renormalisation regions
-
-    psStats *meanStats = psStatsAlloc(meanStat), *stdevStats = psStatsAlloc(stdevStat); // Statistics
-    psVector *buffer = NULL;
-
-    for (int j = 0; j < yNum; j++) {
-        // Bounds in y
-        int yMin = j * ySize;
-        int yMax = (j + 1) * ySize;
-        for (int i = 0; i < xNum; i++) {
-            // Bounds in x
-            int xMin = i * xSize;
-            int xMax = (i + 1) * xSize;
-
-            psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest
-            psImage *subImage = psImageSubset(image, region); // Sub-image of the image 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, 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
-                psWarning("Unable to measure statistics over %s", regionStr);
-                psFree(regionStr);
-                psErrorClear();
-                psImageInit(subVariance, NAN);
-                if (subMask) {
-                    psImageInit(subMask, maskVal);
-                }
-            } else {
-                double meanVar = psStatsGetValue(meanStats, meanStat); // Mean of variance map
-                double stdev = psStatsGetValue(stdevStats, stdevStat); // Standard deviation of image
-                psTrace("psModules.camera", 3,
-                        "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(subVariance, subVariance, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32));
-            }
-
-            psFree(subImage);
-            psFree(subVariance);
-            psFree(subMask);
-        }
-    }
-    psFree(meanStats);
-    psFree(stdevStats);
-    psFree(rng);
-    psFree(buffer);
-
-    return true;
-}
+
+    psImage *subImage = psImageSubset(variance, psRegionSet(xMin, xMax, yMin, yMax)); // Smaller image
+    psBinaryOp(subImage, subImage, "*", psScalarAlloc(PS_SQR(correction), PS_TYPE_F32));
+    psFree(subImage);
+
+    pmHDU *hdu = pmHDUFromReadout(readout); // HDU for readout
+    if (hdu)  {
+        psString history = NULL;
+        psStringAppend(&history, "Rescaled variance by %6.4f (stdev by %6.4f)",
+                       PS_SQR(correction), correction);
+        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history);
+        psFree(history);
+    }
+
+    return psMetadataAddF32(readout->analysis, PS_LIST_TAIL, PM_READOUT_ANALYSIS_RENORM, 0,
+                            "Renormalisation of variance", PS_SQR(correction));
+}
+