Changeset 29543 for trunk/psModules/src/imcombine/pmSubtractionStamps.c

Timestamp:

Oct 25, 2010, 2:45:41 PM (16 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-20100823

File:

• trunk/psModules/src/imcombine/pmSubtractionStamps.c (modified) (23 diffs)

trunk/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;
 
     // Ensure we're not going to go outside the bounds of the image
@@ -103 +106 @@
     yMax = PS_MIN(numRows - border - 1, yMax);
 
-    if (xMax < xMin) return false;
-    if (yMax < yMin) return false;
+    if (xMax < xMin) {
+        // fprintf (stderr, "%f,%f : x-border\n", xRaw, yRaw);
+        return false;
+    }
+    if (yMax < yMin) {
+        // fprintf (stderr, "%f,%f : y-border\n", xRaw, yRaw);
+        return false;
+    }
 
     psAssert (xMin <= xMax, "x mismatch?");
@@ -118 +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);
                 return false;
             }
@@ -133 +143 @@
     }
 
+    // if (!found) {
+    //  fprintf (stderr, "%f,%f : fails flux test\n", xRaw, yRaw);
+    // }
+
     return found;
 }
@@ -232 +246 @@
     list->flux = NULL;
     list->normFrac = normFrac;
+    list->normValue = NAN;
+    list->window = NULL;
     list->window1 = NULL;
     list->window2 = NULL;
@@ -256 +272 @@
     out->y = NULL;
     out->flux = NULL;
+    out->window = psMemIncrRefCounter(in->window);
     out->window1 = psMemIncrRefCounter(in->window1);
     out->window2 = psMemIncrRefCounter(in->window2);
@@ -331 +348 @@
     stamp->vector = NULL;
     stamp->norm = NAN;
+    stamp->normI1 = NAN;
+    stamp->normI2 = NAN;
+    stamp->normSquare1 = NAN;
+    stamp->normSquare2 = NAN;
+
+    stamp->MxxI1 = NULL;
+    stamp->MyyI1 = NULL;
+    stamp->MxxI2 = NULL;
+    stamp->MyyI2 = NULL;
+
+    stamp->MxxI1raw = NAN;
+    stamp->MyyI1raw = NAN;
+    stamp->MxxI2raw = NAN;
+    stamp->MyyI2raw = NAN;
 
     return stamp;
@@ -431 +462 @@
                 // Take stamps off the top of the (sorted) list
                 for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) {
+                    // fprintf (stderr, "%d : xList: %ld elements\n", i, xList->n);
                     // Chop off the top of the list
                     xList->n = j;
@@ -459 +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 {
@@ -656 +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);
@@ -678 +725 @@
 
     // storage vector for flux data
-    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
     psVector *flux1 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
     psVector *flux2 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
@@ -706 +753 @@
                 psAbort ("failed to generate stats");
             }
-            float f1 = stats->robustMedian;
+            float f1 = stats->sampleMedian;
 
             psStatsInit (stats);
@@ -712 +759 @@
                 psAbort ("failed to generate stats");
             }
-            float f2 = stats->robustMedian;
+            float f2 = stats->sampleMedian;
 
             stamps->window1->kernel[y][x] = f1;
@@ -728 +775 @@
     }
 
+#if 0
+    {
+        psFits *fits = NULL;
+        fits = psFitsOpen ("window1.raw.fits", "w");
+        psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL);
+        psFitsClose (fits);
+        fits = psFitsOpen ("window2.raw.fits", "w");
+        psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL);
+        psFitsClose (fits);
+    }
+#endif
+
     psTrace("psModules.imcombine", 3, "Window total (1): %f, threshold: %f\n", sum1, (1.0 - stamps->normFrac) * sum1);
     psTrace("psModules.imcombine", 3, "Window total (2): %f, threshold: %f\n", sum2, (1.0 - stamps->normFrac) * sum2);
 
-# if (0)
+# if (1)
     // this block attempts to calculate the radius based on the first radial moment
-    bool done1 = false;
-    bool done2 = false;
-    double prior1 = 0.0;
-    double prior2 = 0.0;
+    double Sr1 = 0.0;
+    double Sr2 = 0.0;
+    double Sf1 = 0.0;
+    double Sf2 = 0.0;
     for (int y = -size; y <= size; y++) {
         for (int x = -size; x <= size; x++) {
@@ -750 +809 @@
     float R2 = Sr2 / Sf2;
 
-    stamps->normWindow1 = 2.5*R1;
-    stamps->normWindow1 = 2.5*R2;
+    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).
@@ -778 +839 @@
             // 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;
         }
@@ -785 +846 @@
             // 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;
         }
@@ -833 +894 @@
     {
         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);
@@ -906 +967 @@
         float M2 = pmSubtractionKernelPenaltySingle(kernel->kernel, zeroNull);
 
-        penalty1 = M2 + PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment 
-        penalty2 = M2 + PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+        if (1) {
+            penalty1 = M2 * PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+            penalty2 = M2 * PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+            // penalty1 = M2 + PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+            // penalty2 = M2 + PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+        } else {
+            penalty1 = PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+            penalty2 = PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+        }
     }
     kernels->penalties1->data.F32[index] = kernels->penalty * penalty1;
@@ -915 +983 @@
     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;
@@ -942 +1010 @@
     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);
     }
@@ -983 +1051 @@
 
         if (x < bounds.x0 + size || x > bounds.x1 - size || y < bounds.y0 + size || y > bounds.y1 - size) {
-            psError(PM_ERR_PROG, false, "Stamp %d (%d,%d) is within the region border.\n", i, x, y);
-            return false;
+          psLogMsg("psModules.imcombine", 3, "Stamp %d (%d,%d) is within the region border.\n", i, x, y);
+          stamp->status = PM_SUBTRACTION_STAMP_NONE;
+          continue;
         }