Changeset 30631 for branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionStamps.c

Timestamp:

Feb 14, 2011, 1:05:28 PM (15 years ago)

Author:

watersc1

Message:

Merge of trunk back into branch.

Location:

branches/czw_branch/20101203

Files:

• . (modified) (1 prop)
• psModules/src/imcombine/pmSubtractionStamps.c (modified) (12 diffs)

branches/czw_branch/20101203/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
@@ -122 +124 @@
             if ((image1 && image1->data.F32[y][x] < thresh1) ||
                 (image2 && image2->data.F32[y][x] < thresh2)) {
+                // fprintf (stderr, "%f,%f : thresh\n", xRaw, yRaw);
                 continue;
             }
@@ -366 +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) {
@@ -429 +485 @@
         stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing,
                                              normFrac, sysErr, skyErr);
+    }
+
+    // XXX TEST : dump all stars in the stamps here
+    if (0) {
+        FILE *f = fopen ("stamp.dat", "w");
+        for (int i = 0; i < stamps->num; i++) {
+            psVector *xList = stamps->x->data[i];
+            psVector *yList = stamps->y->data[i]; // Coordinate lists
+            psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes
+
+            for (int j = 0; j < xList->n; j++) {
+                fprintf (f, "%d %d  %f %f  %f\n", i, j, xList->data.F32[j], yList->data.F32[j], fluxList->data.F32[j]);
+            }
+        }
+        fclose (f);
     }
 
@@ -636 +707 @@
     }
 
+    int nTotal = 0;
+
     // Sort the list by flux, with the brightest last
     for (int i = 0; i < numStamps; i++) {
@@ -662 +735 @@
         stamps->y->data[i] = ySorted;
         stamps->flux->data[i] = fluxSorted;
-    }
-
+        nTotal += num;
+    }
+    // fprintf (stderr, "nTotal %d\n", nTotal);
+    
     return stamps;
 }
 
-
-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);
@@ -683 +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));
+        }
     }
 
@@ -790 +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;
@@ -809 +893 @@
     float R2 = Sr2 / Sf2;
 
-    stamps->normWindow1 = 2.0*R1;
-    stamps->normWindow2 = 2.0*R2;
-    psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2);
-
-# 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) {
+    // Compare the Kron Radii (R1 & R2) to above to the FWHMs : if they are too discrepant, we will need to rescale
+    psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii vs FWHMs 1: fwhm: %f, kron %f\n", fwhm1, R1);
+    psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii vs FWHMs 2: fwhm: %f, kron %f\n", fwhm2, R2);
+
+    // XXX CAREFUL : in pmSubtractionMatch.c:703, we rely on this factor of 2.75..
+    stamps->normWindow1 = 2.75*R1;
+    stamps->normWindow2 = 2.75*R2;
+    psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Windows from 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 || stamps->normWindow2 >= size) {
+    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);
+
+    // 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;
+        }
     }
 
@@ -890 +959 @@
         }
     }
-
-#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);
@@ -1128 +1185 @@
             continue;
         }
+       
+        // XXX this is somewhat arbitrary...
+        if (source->errMag > 0.05) continue;
+        if (fabs(source->psfMag - source->apMag) > 0.5) continue;
+
         if (source->modelPSF) {
             x->data.F32[numOut] = source->modelPSF->params->data.F32[PM_PAR_XPOS];
@@ -1135 +1197 @@
             y->data.F32[numOut] = source->peak->yf;
         }
-        // fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]);
         numOut++;
     }