Changeset 30322 for branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.c

Timestamp:

Jan 20, 2011, 11:53:03 AM (15 years ago)

Author:

eugene

Message:

moved subtractionGetStamps to pmSubtractionStamps and renamed as pmSubtractionStampsSelect; added return value to pmSubtractionStampsGetWindow to specify a size retry; if desired normalization window is too large, retry with a larger stamp size; move around structures to prevent annoying build collisions

File:

• branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.c (modified) (7 diffs)

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.c

@@ -27 +27 @@
 #include "pmSource.h"
 
+#include "pmSubtractionTypes.h"
 #include "pmSubtraction.h"
 #include "pmSubtractionStamps.h"
+#include "pmSubtractionVisual.h"
 
 #define STAMP_LIST_BUFFER 20            // Number of stamps to add to list at a time
@@ -367 +369 @@
 }
 
-
-pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image1,
-                                                const psImage *image2, const psImage *subMask,
-                                                const psRegion *region, float thresh1, float thresh2,
-                                                int size, int footprint, float spacing, float normFrac,
-                                                float sysErr, float skyErr, pmSubtractionMode mode)
+bool pmSubtractionStampsSelect(pmSubtractionStampList **stamps, // Stamps to read
+                               const pmReadout *ro1, // Readout 1
+                               const pmReadout *ro2, // Readout 2
+                               const psImage *subMask, // Mask for subtraction, or NULL
+                               psImage *variance,  // Variance map
+                               const psRegion *region, // Region of interest
+                               float thresh1,  // Threshold for stamp finding on readout 1
+                               float thresh2,  // Threshold for stamp finding on readout 2
+                               float stampSpacing, // Spacing between stamps
+                               float normFrac,     // Fraction of flux in window for normalisation window
+                               float sysError,     // Relative systematic error in images
+                               float skyError,     // Relative systematic error in images
+                               int size,         // Kernel half-size
+                               int footprint,     // Convolution footprint for stamps
+                               pmSubtractionMode mode // Mode for subtraction
+    )
+{
+    PS_ASSERT_PTR_NON_NULL(stamps, false);
+    PM_ASSERT_READOUT_NON_NULL(ro1, false);
+    PM_ASSERT_READOUT_NON_NULL(ro2, false);
+    PS_ASSERT_IMAGE_NON_NULL(subMask, false);
+    PS_ASSERT_IMAGE_NON_NULL(variance, false);
+    PS_ASSERT_PTR_NON_NULL(region, false);
+
+    psTrace("psModules.imcombine", 3, "Finding stamps...\n");
+
+    psImage *image1 = ro1 ? ro1->image : NULL, *image2 = ro2 ? ro2->image : NULL; // Images of interest
+
+    *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
+                                      size, footprint, stampSpacing, normFrac, sysError, skyError, mode);
+    if (!*stamps) {
+        psError(psErrorCodeLast(), false, "Unable to find stamps.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
+    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) {
+        psError(psErrorCodeLast(), false, "Unable to extract stamps.");
+        return false;
+    }
+
+    pmSubtractionVisualPlotStamps(*stamps, (pmReadout *) ro1);
+    return true;
+}
+
+pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, 
+                                                const psImage *image1,
+                                                const psImage *image2, 
+                                                const psImage *subMask,
+                                                const psRegion *region, 
+                                                float thresh1, 
+                                                float thresh2,
+                                                int size, 
+                                                int footprint, 
+                                                float spacing, 
+                                                float normFrac,
+                                                float sysErr, 
+                                                float skyErr, 
+                                                pmSubtractionMode mode)
 {
     if (!image1 && !image2) {
@@ -687 +742 @@
 }
 
-
-bool pmSubtractionStampsGetWindow(pmSubtractionStampList *stamps, int kernelSize)
+// we are essentially using aperture photometry to determine the photometric match between the
+// images.  we need to choose an appropriate-sized aperture for this analysis.  If it is too
+// large, the measurement will be noisy (and possibly biased) due to the sky noise.  If it is
+// too small, or inconsistent, the measurement will be biased.  We use Kron-mag like aperture
+// scaled by the first radial moment.
+bool pmSubtractionStampsGetWindow(bool *tryAgain, pmSubtractionStampList *stamps, int kernelSize)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
     PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
 
+    // if we succeed, or fail with an unrecoverable error, do not try again
+    if (tryAgain) {
+        *tryAgain = false;
+    }
+
     int size = stamps->footprint; // Size of postage stamps
 
+    // window for moments calculations downstream
+    psFree (stamps->window);
+    stamps->window = psKernelAlloc(-size, size, -size, size);
+    psImageInit(stamps->window->image, 0.0);
+
+    // window1 and window2 are mean stamp images used here to measure the 
+    // first radial moment, and thus the normalization window
     psFree (stamps->window1);
     stamps->window1 = psKernelAlloc(-size, size, -size, size);
@@ -703 +774 @@
     psImageInit(stamps->window2->image, 0.0);
 
-    // Generate a weighting window based on the fwhms (20% larger than the largest)
-    { 
-        float fwhm1, fwhm2;
-
-        // XXX this is annoyingly hack-ish
-        pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
-
-        float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35;
-
-        psFree (stamps->window);
-        stamps->window = psKernelAlloc(-size, size, -size, size);
-        psImageInit(stamps->window->image, 0.0);
-
-        for (int y = -size; y <= size; y++) {
-            for (int x = -size; x <= size; x++) {
-                stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma));
-            }
-        }
+    // Generate an initial weighting window based on the fwhms (50% larger than the largest)
+    float fwhm1, fwhm2;
+
+    // XXX this is annoyingly hack-ish
+    pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
+    
+    float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35;
+    
+    for (int y = -size; y <= size; y++) {
+        for (int x = -size; x <= size; x++) {
+            stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma));
+        }
     }
 
@@ -810 +875 @@
     psTrace("psModules.imcombine", 3, "Window total (2): %f, threshold: %f\n", sum2, (1.0 - stamps->normFrac) * sum2);
 
-# if (1)
-    // this block attempts to calculate the radius based on the first radial moment
+    // attempt to calculate the normalization window based on the first radial moment
     double Sr1 = 0.0;
     double Sr2 = 0.0;
@@ -833 +897 @@
     psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2);
 
+    // if the calculated normWindows are too large, we will fall off the stamps.  In this case, we need to try again.
+    if ((stamps->normWindow1 > size) || (stamps->normWindow2 > size)) { 
+        if (tryAgain) {
+            *tryAgain = true;
+        }
+        psFree (stats);
+        psFree (flux1);
+        psFree (flux2);
+        psFree (norm1);
+        psFree (norm2);
+        return false; 
+    }
+
+    // this is an unrecoverable error : something really bogus in the data
+    if (stamps->normWindow1 == 0) {
+        psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (1).");
+        psFree (stats);
+        psFree (flux1);
+        psFree (flux2);
+        psFree (norm1);
+        psFree (norm2);
+        return false;
+    }
+    if (stamps->normWindow2 == 0) {
+        psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (2).");
+        psFree (stats);
+        psFree (flux1);
+        psFree (flux2);
+        psFree (norm1);
+        psFree (norm2);
+        return false;
+    }
+
     // Generate a weighting window based on the kron radii
-    { 
-        float radius = 2.0 * PS_MAX(R1, R2);
-
-        psImageInit(stamps->window->image, 0.0);
-
-        for (int y = -size; y <= size; y++) {
-            for (int x = -size; x <= size; x++) {
-                if (hypot(x,y) > radius) continue;
-                stamps->window->kernel[y][x] = 1.0;
-            }
-        }
-    }
-
-    // complain if normWindow1 or normWindow2 are larger than size
-    psAssert(stamps->normWindow1 < size, "norm Window 1 too large");
-    psAssert(stamps->normWindow2 < size, "norm Window 2 too large");
-
-# else
-    // XXX : this block attempts to calculate the radius by looking at the curve of growth (or something vaguely equivalent).
-    // It did not do very well (though a true curve-of-growth analysis might be better...)
-    bool done1 = false;
-    bool done2 = false;
-    double prior1 = 0.0;
-    double prior2 = 0.0;
-    double delta1o = 1.0;
-    double delta2o = 1.0;
-    for (int radius = 1; radius <= size && !(done1 && done2); radius++) {
-        double within1 = 0.0;
-        double within2 = 0.0;
-        for (int y = -radius; y <= radius; y++) {
-            for (int x = -radius; x <= radius; x++) {
-                if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
-                    within1 += stamps->window1->kernel[y][x];
-                }
-                if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
-                    within2 += stamps->window2->kernel[y][x];
-                }
-            }
-        }
-        double delta1 = (within1 - prior1) / within1;
-        if (!done1 && (fabs(delta1) < stamps->normFrac)) {
-            // interpolate to the radius at which delta2 is normFrac:
-            stamps->normWindow1 = radius - (stamps->normFrac - delta1) / (delta1o - delta1);
-            psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1);
-            done1 = true;
-        }
-        double delta2 = (within2 - prior2) / within2;
-        if (!done2 && (fabs(delta2) < stamps->normFrac)) {
-            // interpolate to the radius at which delta2 is normFrac:
-            stamps->normWindow2 = radius - (stamps->normFrac - delta2) / (delta2o - delta2);
-            psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2);
-            done2 = true;
-        }
-        psTrace("psModules.imcombine", 5, "Radius %d: %f (%f) and %f (%f)\n", radius, within1, delta1, within2, delta2);
-
-        prior1 = within1;
-        prior2 = within2;
-        delta1o = delta1;
-        delta2o = delta2;
-
-        // if (!done1 && (within1 > (1.0 - stamps->normFrac) * sum1)) {
-        //     stamps->normWindow1 = radius;
-        //     done1 = true;
-        // }
-        // if (!done2 && (within2 > (1.0 - stamps->normFrac) * sum2)) {
-        //     stamps->normWindow2 = radius;
-        //     done2 = true;
-        // }
-
-    }
-# endif
-
-    psTrace("psModules.imcombine", 3, "Normalisation window radii set to %f and %f\n", stamps->normWindow1, stamps->normWindow2);
-    if (stamps->normWindow1 == 0 || stamps->normWindow1 >= size) {
-        psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (1).");
-        return false;
-    }
-    if (stamps->normWindow2 == 0 || stamps->normWindow2 >= size) {
-        psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (2).");
-        return false;
+    float radius = 2.0 * PS_MAX(R1, R2);
+    psImageInit(stamps->window->image, 0.0);
+
+    // we use a top-hat window for the moments analysis
+    for (int y = -size; y <= size; y++) {
+        for (int x = -size; x <= size; x++) {
+            if (hypot(x,y) > radius) continue;
+            stamps->window->kernel[y][x] = 1.0;
+        }
     }
 
@@ -928 +954 @@
         }
     }
-
-#if 0
-    {
-        psFits *fits = NULL;
-        fits = psFitsOpen ("window1.norm.fits", "w");
-        psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL);
-        psFitsClose (fits);
-        fits = psFitsOpen ("window2.norm.fits", "w");
-        psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL);
-        psFitsClose (fits);
-    }
-#endif
 
     psFree (stats);