Changeset 29449 for branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c

Timestamp:

Oct 17, 2010, 9:10:35 AM (16 years ago)

Author:

eugene

Message:

generate a weight window based on the fwhms; pass around fwhm to ensure the window can be calculated; plug some leaks; better test outputs

File:

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

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

@@ -51 +51 @@
     psFree(list->y);
     psFree(list->flux);
+    psFree(list->window);
     psFree(list->window1);
     psFree(list->window2);
@@ -96 +97 @@
 
     // fprintf (stderr, "xMin, xMax: %d, %d -> ", xMin, xMax);
-    float xRaw = xMin;
-    float yRaw = yMin;
+    // float xRaw = xMin;
+    // float yRaw = yMin;
 
     // Ensure we're not going to go outside the bounds of the image
@@ -126 +127 @@
             if (subMask && subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] &
                 (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) {
-              fprintf (stderr, "%f,%f : masked\n", xRaw, yRaw);
+                // fprintf (stderr, "%f,%f : masked\n", xRaw, yRaw);
                 return false;
             }
@@ -142 +143 @@
     }
 
-    if (!found) {
-      fprintf (stderr, "%f,%f : fails flux test\n", xRaw, yRaw);
-    }
+    // if (!found) {
+    //  fprintf (stderr, "%f,%f : fails flux test\n", xRaw, yRaw);
+    // }
+
     return found;
 }
@@ -245 +247 @@
     list->normFrac = normFrac;
     list->normValue = NAN;
+    list->window = NULL;
     list->window1 = NULL;
     list->window2 = NULL;
@@ -269 +272 @@
     out->y = NULL;
     out->flux = NULL;
+    out->window = psMemIncrRefCounter(in->window);
     out->window1 = psMemIncrRefCounter(in->window1);
     out->window2 = psMemIncrRefCounter(in->window2);
@@ -487 +491 @@
                                             yList->data.F32[j], yList->data.F32[j], numCols, numRows, border);
 #endif
-                    if (0 && goodStamp) {
-                        fprintf (stderr, "find: %6.1f < %6.1f < %6.1f, %6.1f < %6.1f %6.1f\n",
-                                 region->x0 + size, xStamp, region->x1 - size,
-                                 region->y0 + size, yStamp, region->y1 - size);
-                    }
                 }
             } else {
@@ -684 +683 @@
     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 normalizations for each stamp
     psVector *norm1 = psVectorAlloc(stamps->num, PS_TYPE_F32);
@@ -759 +778 @@
     {
         psFits *fits = NULL;
-        fits = psFitsOpen ("window1.fits", "w");
+        fits = psFitsOpen ("window1.raw.fits", "w");
         psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL);
         psFitsClose (fits);
-        fits = psFitsOpen ("window2.fits", "w");
+        fits = psFitsOpen ("window2.raw.fits", "w");
         psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL);
         psFitsClose (fits);
@@ -818 +837 @@
             // interpolate to the radius at which delta2 is normFrac:
             stamps->normWindow1 = radius - (stamps->normFrac - delta1) / (delta1o - delta1);
-            fprintf (stderr, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1);
+            psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1);
             done1 = true;
         }
@@ -825 +844 @@
             // interpolate to the radius at which delta2 is normFrac:
             stamps->normWindow2 = radius - (stamps->normFrac - delta2) / (delta2o - delta2);
-            fprintf (stderr, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2);
+            psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2);
             done2 = true;
         }
@@ -873 +892 @@
     {
         psFits *fits = NULL;
-        fits = psFitsOpen ("window1.fits", "w");
+        fits = psFitsOpen ("window1.norm.fits", "w");
         psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL);
         psFitsClose (fits);
-        fits = psFitsOpen ("window2.fits", "w");
+        fits = psFitsOpen ("window2.norm.fits", "w");
         psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL);
         psFitsClose (fits);
@@ -962 +981 @@
     psAssert (isfinite(kernels->penalties2->data.F32[index]), "invalid penalty");
 
-    fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2);
+    // fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2);
 
     return true;
@@ -989 +1008 @@
     penalty *= 1.0 / sum2;
 
-    if (1) {
-        fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2);
+    if (0) {
+        // fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2);
         // psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty);
     }