Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionAnalysis.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionAnalysis.c	(revision 26592)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionAnalysis.c	(revision 26593)
@@ -117,4 +117,8 @@
     }
 
+    // sample difference images
+    { 
+        psMetadataAddArray(analysis, PS_LIST_TAIL, "SUBTRACTION.SAMPLE.STAMP.SET", PS_META_DUPLICATE_OK, "Sample Difference Stamps", kernels->sampleStamps);
+    }
 
 #ifdef TESTING
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26592)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26593)
@@ -1331,4 +1331,44 @@
 }
 
+bool pmSubtractionResidualStats(psVector *fSigRes, psVector *fMaxRes, psVector *fMinRes, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) {
+
+    // XXX measure some useful stats on the residuals
+    float sum = 0.0;
+    float peak = 0.0;
+    for (int y = - footprint; y <= footprint; y++) {
+	for (int x = - footprint; x <= footprint; x++) {
+	    sum += 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
+	    peak = PS_MAX(peak, 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm));
+	}
+    }
+
+    // only count pixels with more than X% of the source flux
+    // calculate stdev(dflux)
+    float dflux1 = 0.0;
+    float dflux2 = 0.0;
+    int npix = 0;
+
+    float dmax = 0.0;
+    float dmin = 0.0;
+
+    for (int y = - footprint; y <= footprint; y++) {
+	for (int x = - footprint; x <= footprint; x++) {
+	    float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
+	    if (dflux < 0.02*sum) continue;
+	    dflux1 += residual->kernel[y][x];
+	    dflux2 += PS_SQR(residual->kernel[y][x]);
+	    dmax = PS_MAX(residual->kernel[y][x], dmax);
+	    dmin = PS_MIN(residual->kernel[y][x], dmin);
+	    npix ++;
+	}
+    }
+    float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
+    // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak);
+    psVectorAppend(fSigRes, sigma/sum);
+    psVectorAppend(fMaxRes, dmax/peak);
+    psVectorAppend(fMinRes, dmin/peak);
+    return true;
+}
+
 psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps,
                                            pmSubtractionKernels *kernels)
@@ -1353,4 +1393,35 @@
     psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+
+    // we want to save the residual images for the 9 brightest stamps.
+    // identify the 9 brightest stamps
+    psVector *keepStamps  = psVectorAlloc(stamps->num, PS_TYPE_S32);
+    psVectorInit (keepStamps, 0);
+    { 
+	psVector *flux  = psVectorAlloc(stamps->num, PS_TYPE_F32);
+	psVectorInit (flux, 0.0);
+	
+	for (int i = 0; i < stamps->num; i++) {
+	    pmSubtractionStamp *stamp = stamps->stamps->data[i];
+	    if (!isfinite(stamp->flux)) continue;
+	    flux->data.F32[i] = stamp->flux;
+	}
+
+	psVector *index = psVectorSortIndex(NULL, flux);
+	for (int i = 0; (i < stamps->num) && (i < 9); i++) {
+	    int n = stamps->num - i - 1;
+	    keepStamps->data.S32[index->data.S32[n]] = 1;
+	    fprintf (stderr, "keeping %d (%d of %d)\n", index->data.S32[n], n, 9);
+	}
+	psFree (flux);
+	psFree (index);
+
+	// this function is called multiple times in the iteration, but 
+	// we only know after the interation is done if we will try again.
+	// therefore we must save the sample each time, and blow away the old one
+	// if it exists.
+	psFree (kernels->sampleStamps);
+	kernels->sampleStamps = psArrayAllocEmpty(9);
+    }
 
     for (int i = 0; i < stamps->num; i++) {
@@ -1427,39 +1498,12 @@
             }
 
-            // XXX measure some useful stats on the residuals
-            float sum = 0.0;
-            float peak = 0.0;
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = - footprint; x <= footprint; x++) {
-                    sum += 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
-                    peak = PS_MAX(peak, 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm));
-                }
-            }
-
-            // only count pixels with more than X% of the source flux
-            // calculate stdev(dflux)
-            float dflux1 = 0.0;
-            float dflux2 = 0.0;
-            int npix = 0;
-
-            float dmax = 0.0;
-            float dmin = 0.0;
-
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = - footprint; x <= footprint; x++) {
-                    float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
-                    if (dflux < 0.02*sum) continue;
-                    dflux1 += residual->kernel[y][x];
-                    dflux2 += PS_SQR(residual->kernel[y][x]);
-                    dmax = PS_MAX(residual->kernel[y][x], dmax);
-                    dmin = PS_MIN(residual->kernel[y][x], dmin);
-                    npix ++;
-                }
-            }
-            float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
-            // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak);
-            psVectorAppend(fSigRes, sigma/sum);
-            psVectorAppend(fMaxRes, dmax/peak);
-            psVectorAppend(fMinRes, dmin/peak);
+	    if (keepStamps->data.S32[i]) {
+		psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32);
+		psArrayAdd (kernels->sampleStamps, 9, sample);
+		psFree (sample);
+	    }
+
+	    pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint);
+
         } else {
             // Dual convolution
@@ -1490,4 +1534,11 @@
                 }
             }
+	    if (keepStamps->data.S32[i]) {
+		psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32);
+		psArrayAdd (kernels->sampleStamps, 9, sample);
+		psFree (sample);
+	    }
+
+	    pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint);
         }
 
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c	(revision 26592)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c	(revision 26593)
@@ -29,4 +29,5 @@
     psFree(kernels->solution1);
     psFree(kernels->solution2);
+    psFree(kernels->sampleStamps);
 }
 
@@ -226,5 +227,5 @@
     if (!forceZeroNull && ((uOrder % 2) || (vOrder % 2))) {
         // Odd function
-        scale2D = 1.0 / sqrt(sum2);
+	scale2D = 1.0 / sqrt(sum2);
     }
 
@@ -573,4 +574,5 @@
     kernels->rms = NAN;
     kernels->numStamps = 0;
+    kernels->sampleStamps = NULL;
 
     kernels->fSigResMean  = NAN;
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h	(revision 26592)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h	(revision 26593)
@@ -55,4 +55,5 @@
     float fMinResMean;			///< mean fractional negative swing in residuals
     float fMinResStdev;			///< stdev of fractional negative swing in residuals
+    psArray *sampleStamps;		///< array of brightest set of stamps for output visualizations
 } pmSubtractionKernels;