Changeset 25338

Timestamp:

Sep 10, 2009, 7:11:55 PM (17 years ago)

Author:

Paul Price

Message:

Only convolving the pixels of interest, rather than the entire image. On the plus side, this runs a LOT faster. On the minus side, the code has to be kept in sync with psphotSignificanceImage.c.

File:

• branches/pap/psphot/src/psphotEfficiency.c (modified) (10 diffs)

branches/pap/psphot/src/psphotEfficiency.c

@@ -15 +15 @@
                      float *minFlux,          // Minimum flux for fake sources, to return
                      float *norm,             // Normalisation of PSF (conversion: peak --> integrated flux)
+                     float *covarFactor,// Covariance factor
                      const pmReadout *ro,     // Readout of interest
                      const pmPSF *psf,        // Point-spread function
@@ -29 +30 @@
     assert(minFlux);
     assert(norm);
+    assert(covarFactor);
 
     psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing
     psKernel *newCovar = psImageCovarianceCalculate(kernel, ro->covariance); // Covariance matrix
     psFree(kernel);
-    float factor = psImageCovarianceFactor(newCovar); // Variance factor
+    *covarFactor = psImageCovarianceFactor(newCovar); // Variance factor
     psFree(newCovar);
 
@@ -68 +70 @@
     // with itself has a normalisation of 1/2pi(w^2+w^2) = 1/4pi.w^2.  Therefore:
     // I_smooth = Flux / 2*norm.
-    float peakLim = thresh * sqrtf(meanVar * factor); // Limiting peak value in smoothed image
+    float peakLim = thresh * sqrtf(meanVar * *covarFactor); // Limiting peak value in smoothed image
     float fluxLim = 2.0 * *norm * peakLim; // Limiting flux in original
     *magLim = -2.5 * log10f(fluxLim);
@@ -106 +108 @@
     *ySrc = psImageRecycle(*ySrc, numSources, numBins, PS_TYPE_F32);
 
-    // Master list, for image creation
     psVector *xAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32);
     psVector *yAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32);
@@ -132 +133 @@
         psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image");
         psFree(fakeRO);
-        psFree(xAll);
-        psFree(yAll);
         psFree(magAll);
         return false;
     }
-    psFree(xAll);
-    psFree(yAll);
     psFree(magAll);
 
@@ -187 +184 @@
         return NULL;
     }
+    float minGauss = psMetadataLookupF32(NULL, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel
+    if (!isfinite(minGauss)) {
+        psWarning("PEAKS_MIN_GAUSS is not set in recipe; using default value");
+        minGauss = 0.5;
+    }
 
     float smoothSigma = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrtf(2.0*log(2.0))); // Gaussian smoothing sigma
     int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image
+    int numBins = magOffsets->n;                          // Number of bins
 
     psImageMaskType maskVal = psMetadataLookupImageMask(NULL, recipe, "MASK.PSPHOT"); // Value to mask
@@ -220 +223 @@
     float minFlux;                      // Minimum flux for fake sources
     float norm;                         // Normalisation of PSF
-    if (!effLimit(&magLim, &radius, &minFlux, &norm, readout,
-                   psf, thresh, smoothSigma, smoothNsigma, maskVal)) {
+    float covarFactor;                  // Covariance factor
+    if (!effLimit(&magLim, &radius, &minFlux, &norm, &covarFactor, readout,
+                  psf, thresh, smoothSigma, smoothNsigma, maskVal)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to determine limits for image");
         return false;
@@ -233 +237 @@
 
     psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row
-    if (!effGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, radius, minFlux)) {
+    if (!effGenerate(&xFake, &yFake, readout, psf, magOffsets,
+                     numSources, magLim, radius, minFlux)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources");
         psFree(xFake);
@@ -247 +252 @@
 
     // XXX Could speed this up significantly by only convolving the central pixels of each fake source
-    psImage *significance = psphotSignificanceImage(readout, recipe, 2, maskVal);
-    if (!significance) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate significance image");
-        psFree(xFake);
-        psFree(yFake);
-        return false;
-    }
-
-#ifdef TESTING
-    psphotSaveImage(NULL, significance, "fake_sig.fits");
-#endif
-
-    thresh *= thresh;                   // Significance image is actually significance-squared
-
-    int numBins = magOffsets->n;                          // Number of bins
+    psVector *significance = NULL;       // Significance image
+    {
+        int num = numSources * numBins; // Total number of sources
+        // Pixel coordinates of sources
+        psVector *xPix = psVectorAlloc(num, PS_TYPE_S32);
+        psVector *yPix = psVectorAlloc(num, PS_TYPE_S32);
+        for (int i = 0, index = 0; i < numBins; i++) {
+            for (int j = 0; j < numSources; j++, index++) {
+                float x = xFake->data.F32[i][j];
+                float y = yFake->data.F32[i][j];
+                xPix->data.S32[index] = PS_MIN(x + 0.5, numCols - 1);
+                yPix->data.S32[index] = PS_MIN(y + 0.5, numRows - 1);
+            }
+        }
+
+        psVector *convImage = psImageSmoothMaskPixels(readout->image, readout->mask, maskVal, xPix, yPix,
+                                                      smoothSigma, smoothNsigma, minGauss); // Convolved image
+        if (!convImage) {
+            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to smooth image pixels");
+            psFree(xPix);
+            psFree(yPix);
+            psFree(xFake);
+            psFree(yFake);
+            return false;
+        }
+        psVector *convVar = psImageSmoothMaskPixels(readout->variance, readout->mask, maskVal, xPix, yPix,
+                                                    smoothSigma, smoothNsigma, minGauss); // Convolved varianc
+        if (!convVar) {
+            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to smooth variance pixels");
+            psFree(xPix);
+            psFree(yPix);
+            psFree(xFake);
+            psFree(yFake);
+            return false;
+        }
+
+        float factor = 1.0 / covarFactor; // Correction for covariance
+
+        const psImage *mask = readout->mask;  // Mask for readout
+        for (int i = 0; i < num; i++) {
+            float imageVal = convImage->data.F32[i]; // Image value
+            float varVal = convVar->data.F32[i]; // Variance value
+            int x = xPix->data.S32[i], y = yPix->data.S32[i]; // Coordinates
+            if (imageVal < 0 || varVal <= 0 || (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) {
+                convImage->data.F32[i] = 0.0;
+            } else {
+                convImage->data.F32[i] = factor * PS_SQR(imageVal) / varVal;
+            }
+        }
+
+        psFree(xPix);
+        psFree(yPix);
+
+        significance = convImage;
+        psFree(convVar);
+    }
+
+    thresh *= thresh;                   // "Significance" is actually significance-squared
+
     psVector *count = psVectorAlloc(numBins, PS_TYPE_S32); // Number of sources found in each bin
     psArray *fakeSources = psArrayAlloc(numSources);            // Fake sources in each bin
     psArray *fakeSourcesAll = psArrayAllocEmpty(numSources * numBins); // All fake sources
-    for (int i = 0; i < numBins; i++) {
+    for (int i = 0, index = 0; i < numBins; i++) {
         int numFound = 0;               // Number found
 
@@ -281 +330 @@
 
         psArray *sources = psArrayAllocEmpty(numSources); // Sources in this bin
-        for (int j = 0; j < numSources; j++) {
+        for (int j = 0; j < numSources; j++, index++) {
             // Coordinates of interest
-            float x = xFake->data.F32[i][j];
-            float y = yFake->data.F32[i][j];
-            int xPix = PS_MIN(x + 0.5, numCols - 1), yPix = PS_MIN(y + 0.5, numRows - 1); // Pixel coordinates
-            float sig = significance->data.F32[yPix][xPix]; // Significance of pixel
+            float sig = significance->data.F32[index]; // Significance of pixel
             if (sig > thresh) {
                 pmSource *source = pmSourceAlloc(); // Fake source
+                float x = xFake->data.F32[i][j];
+                float y = yFake->data.F32[i][j];
                 source->peak = pmPeakAlloc(x, y, sig, PM_PEAK_LONE);
                 if (!pmSourceDefinePixels(source, readout, x, y, sourceRadius)) {