Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionStamps.c	(revision 29019)
+++ trunk/psModules/src/imcombine/pmSubtractionStamps.c	(revision 29543)
@@ -51,4 +51,5 @@
     psFree(list->y);
     psFree(list->flux);
+    psFree(list->window);
     psFree(list->window1);
     psFree(list->window2);
@@ -96,4 +97,6 @@
 
     // 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,6 +106,12 @@
     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,4 +127,5 @@
             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,4 +143,8 @@
     }
 
+    // if (!found) {
+    // 	fprintf (stderr, "%f,%f : fails flux test\n", xRaw, yRaw);
+    // }
+
     return found;
 }
@@ -232,4 +246,6 @@
     list->flux = NULL;
     list->normFrac = normFrac;
+    list->normValue = NAN;
+    list->window = NULL;
     list->window1 = NULL;
     list->window2 = NULL;
@@ -256,4 +272,5 @@
     out->y = NULL;
     out->flux = NULL;
+    out->window = psMemIncrRefCounter(in->window);
     out->window1 = psMemIncrRefCounter(in->window1);
     out->window2 = psMemIncrRefCounter(in->window2);
@@ -331,4 +348,18 @@
     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,4 +462,5 @@
                 // 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,9 +491,4 @@
                                             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,4 +683,24 @@
     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,5 +725,5 @@
 
     // 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,5 +753,5 @@
                 psAbort ("failed to generate stats");
             }
-            float f1 = stats->robustMedian;
+            float f1 = stats->sampleMedian;
 
             psStatsInit (stats);
@@ -712,5 +759,5 @@
                 psAbort ("failed to generate stats");
             }
-            float f2 = stats->robustMedian;
+            float f2 = stats->sampleMedian;
 
             stamps->window1->kernel[y][x] = f1;
@@ -728,13 +775,25 @@
     }
 
+#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,6 +809,8 @@
     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,5 +839,5 @@
 	    // 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,5 +846,5 @@
 	    // 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,8 +894,8 @@
     {
 	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,6 +967,13 @@
 	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,5 +983,5 @@
     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,6 +1010,6 @@
     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,6 +1051,7 @@
 
         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;
         }