Changeset 26747 for branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionStamps.c

Timestamp:

Jan 31, 2010, 5:00:42 PM (16 years ago)

Author:

eugene

Message:

updates relative to 20091201, fixes for all psphot variants

Location:

branches/eam_branches/psModules.stack.20100120

Files:

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

branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionStamps.c

@@ -176 +176 @@
 
 pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region,
-                                                    int footprint, float spacing, float sysErr, float skyErr)
+                                                    int footprint, float spacing, float normFrac,
+                                                    float sysErr, float skyErr)
 {
     pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return
@@ -217 +218 @@
     list->flux = NULL;
     list->window = NULL;
+    list->normFrac = normFrac;
+    list->normWindow = 0;
     list->footprint = footprint;
     list->sysErr = sysErr;
@@ -239 +242 @@
     out->window = psMemIncrRefCounter(in->window);
     out->footprint = in->footprint;
+    out->normWindow = in->normWindow;
 
     for (int i = 0; i < num; i++) {
@@ -308 +312 @@
     stamp->matrix = NULL;
     stamp->vector = NULL;
+    stamp->norm = NAN;
 
     return stamp;
@@ -316 +321 @@
                                                 const psImage *image2, const psImage *subMask,
                                                 const psRegion *region, float thresh1, float thresh2,
-                                                int size, int footprint, float spacing, float sysErr, float skyErr,
-                                                pmSubtractionMode mode)
+                                                int size, int footprint, float spacing, float normFrac,
+                                                float sysErr, float skyErr, pmSubtractionMode mode)
 {
     if (!image1 && !image2) {
@@ -373 +378 @@
 
     if (!stamps) {
-        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr, skyErr);
+        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing,
+                                             normFrac, sysErr, skyErr);
     }
 
@@ -407 +413 @@
                 // Take stamps off the top of the (sorted) list
                 for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) {
-                    int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre
-
                     // Chop off the top of the list
                     xList->n = j;
@@ -414 +418 @@
                     fluxList->n = j;
 
+#if 0
                     // Fish around a bit to see if we can find a pixel that isn't masked
+                    // This is not a good idea if we're using the window feature
                     psTrace("psModules.imcombine", 7, "Searching for stamp %d around %d,%d\n",
                             i, xCentre, yCentre);
 
                     // Search bounds
+                    int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre
                     int search = footprint - size; // Search radius
                     int xMin = PS_MAX(border, xCentre - search);
@@ -427 +434 @@
                     goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
                                             subMask, xMin, xMax, yMin, yMax, numCols, numRows, border);
-
-                    // XXX reset for a test:
-                    xStamp = xList->data.F32[j];
-                    yStamp = yList->data.F32[j];
+#else
+                    // Only search the exact centre pixel
+                    goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
+                                            subMask, xList->data.F32[j], xList->data.F32[j],
+                                            yList->data.F32[j], yList->data.F32[j], numCols, numRows, border);
+#endif
                     // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp);
                 }
@@ -486 +495 @@
                                                const psImage *image, const psImage *subMask,
                                                const psRegion *region, int size, int footprint,
-                                               float spacing, float sysErr, float skyErr, pmSubtractionMode mode)
+                                               float spacing, float normFrac, float sysErr, float skyErr,
+                                               pmSubtractionMode mode)
 
 {
@@ -507 +517 @@
     pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
                                                                  region, footprint, spacing,
-                                                                 sysErr, skyErr); // Stamp list
+                                                                 normFrac, sysErr, skyErr); // Stamp list
     int numStamps = stamps->num;        // Number of stamps
 
@@ -612 +622 @@
     PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
 
-    int size = kernelSize + stamps->footprint; // Size of postage stamps
+    int size = stamps->footprint; // Size of postage stamps
 
     psFree (stamps->window);
@@ -641 +651 @@
     // storage vector for flux data
     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    psVector *flux = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
-
-    double maxValue = 0.0;
+    psVector *flux1 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
+    psVector *flux2 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
 
     // generate the window pixels
+    double sum = 0.0;                   // Sum inside the window
+    float maxValue = 0.0;               // Maximum value, for normalisation
     for (int y = -size; y <= size; y++) {
         for (int x = -size; x <= size; x++) {
 
-            flux->n = 0;
+            flux1->n = 0;
+            flux2->n = 0;
             for (int i = 0; i < stamps->num; i++) {
                 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
@@ -656 +668 @@
                 if (!stamp->image2) continue;
 
-                psVectorAppend (flux, stamp->image1->kernel[y][x] / norm1->data.F32[i]);
-                psVectorAppend (flux, stamp->image2->kernel[y][x] / norm2->data.F32[i]);
+                psVectorAppend(flux1, stamp->image1->kernel[y][x] / norm1->data.F32[i]);
+                psVectorAppend(flux2, stamp->image2->kernel[y][x] / norm2->data.F32[i]);
             }
 
             psStatsInit (stats);
-            if (!psVectorStats (stats, flux, NULL, NULL, 0)) {
+            if (!psVectorStats (stats, flux1, NULL, NULL, 0)) {
                 psAbort ("failed to generate stats");
             }
-            stamps->window->kernel[y][x] = stats->robustMedian;
-            if (maxValue < stats->robustMedian) {
-                maxValue = stats->robustMedian;
-            }
-        }
+            float f1 = stats->robustMedian;
+            psStatsInit (stats);
+            if (!psVectorStats (stats, flux2, NULL, NULL, 0)) {
+                psAbort ("failed to generate stats");
+            }
+            float f2 = stats->robustMedian;
+
+            stamps->window->kernel[y][x] = f1 + f2;
+            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) {
+                sum += stamps->window->kernel[y][x];
+            }
+            maxValue = PS_MAX(maxValue, stamps->window->kernel[y][x]);
+        }
+    }
+
+    psTrace("psModules.imcombine", 3, "Window total: %f, threshold: %f\n",
+            sum, (1.0 - stamps->normFrac) * sum);
+    bool done = false;
+    for (int radius = 1; radius <= size && !done; radius++) {
+        double within = 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)) {
+                    within += stamps->window->kernel[y][x];
+                }
+            }
+        }
+        psTrace("psModules.imcombine", 5, "Radius %d: %f\n", radius, within);
+        if (within > (1.0 - stamps->normFrac) * sum) {
+            stamps->normWindow = radius;
+            done = true;
+        }
+    }
+
+    psTrace("psModules.imcombine", 3, "Normalisation window radius set to %d\n", stamps->normWindow);
+    if (stamps->normWindow == 0 || stamps->normWindow >= size) {
+        psError(PS_ERR_UNKNOWN, true, "Unable to determine normalisation window size.");
+        return false;
     }
 
@@ -678 +723 @@
     }
 
-# ifdef TESTING
+#if 0
     {
         psFits *fits = psFitsOpen ("window.fits", "w");
@@ -684 +729 @@
         psFitsClose (fits);
     }
-# endif
+#endif
 
     psFree (stats);
-    psFree (flux);
+    psFree (flux1);
+    psFree(flux2);
     psFree (norm1);
     psFree (norm2);
@@ -694 +740 @@
 
 bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,
-                                psImage *variance, int kernelSize)
+                                psImage *variance, int kernelSize, psRegion bounds)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
@@ -711 +757 @@
     }
 
-    int numCols = image1->numCols, numRows = image1->numRows; // Size of images
     int size = kernelSize + stamps->footprint; // Size of postage stamps
 
@@ -720 +765 @@
         }
 
-        if (isnan(stamp->xNorm)) {
-            stamp->xNorm = 2.0 * (stamp->x - (float)numCols/2.0) / (float)numCols;
-        }
-        if (isnan(stamp->yNorm)) {
-            stamp->yNorm = 2.0 * (stamp->y - (float)numRows/2.0) / (float)numRows;
-        }
+        p_pmSubtractionPolynomialNormCoords(&stamp->xNorm, &stamp->yNorm, stamp->x, stamp->y,
+                                            bounds.x0, bounds.x1, bounds.y0, bounds.y1);
 
         int x = stamp->x + 0.5, y = stamp->y + 0.5; // Stamp coordinates
-        if (x < size || x > numCols - size || y < size || y > numRows - size) {
-            psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the image border.\n", i, x, y);
+        if (x < bounds.x0 + size || x > bounds.x1 - size || y < bounds.y0 + size || y > bounds.y1 - size) {
+            psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the region border.\n", i, x, y);
             return false;
         }
@@ -788 +829 @@
 pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image,
                                                           const psImage *subMask, const psRegion *region,
-                                                          int size, int footprint, float spacing, float sysErr, float skyErr,
+                                                          int size, int footprint, float spacing,
+                                                          float normFrac, float sysErr, float skyErr,
                                                           pmSubtractionMode mode)
 {
@@ -819 +861 @@
 
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size,
-                                                            footprint, spacing, sysErr, skyErr, mode); // Stamps
+                                                            footprint, spacing, normFrac,
+                                                            sysErr, skyErr, mode); // Stamps
     psFree(x);
     psFree(y);
@@ -833 +876 @@
 pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image,
                                                        const psImage *subMask, const psRegion *region,
-                                                       int size, int footprint, float spacing, float sysErr, float skyErr,
-                                                       pmSubtractionMode mode)
+                                                       int size, int footprint, float spacing, float normFrac,
+                                                       float sysErr, float skyErr, pmSubtractionMode mode)
 {
     PS_ASSERT_STRING_NON_EMPTY(filename, NULL);
@@ -851 +894 @@
 
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint,
-                                                            spacing, sysErr, skyErr, mode);
+                                                            spacing, normFrac, sysErr, skyErr, mode);
     psFree(data);