Index: /trunk/ppSub/src/Makefile.am
===================================================================
--- /trunk/ppSub/src/Makefile.am	(revision 29002)
+++ /trunk/ppSub/src/Makefile.am	(revision 29003)
@@ -39,4 +39,5 @@
 	ppSubDefineOutput.c		\
 	ppSubExtras.c			\
+	ppSubFlagNeighbors.c		\
 	ppSubMakePSF.c			\
 	ppSubMatchPSFs.c		\
Index: /trunk/ppSub/src/ppSub.h
===================================================================
--- /trunk/ppSub/src/ppSub.h	(revision 29002)
+++ /trunk/ppSub/src/ppSub.h	(revision 29003)
@@ -175,4 +175,8 @@
     );
 
+bool ppSubFlagNeighbors(pmConfig *config, pmFPAview *view, psArray *sources, bool matchRef);
+bool ppSubMatchSources (psArray *objects, psArray *sources, float RADIUS, float MIN_SN);
+bool ppSubSetSourceImageIDs (psArray *sources, int imageID);
+
 ///@}
 #endif
Index: /trunk/ppSub/src/ppSubBackground.c
===================================================================
--- /trunk/ppSub/src/ppSubBackground.c	(revision 29002)
+++ /trunk/ppSub/src/ppSubBackground.c	(revision 29003)
@@ -67,15 +67,17 @@
     for (int y = 0; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
+	    // special case 1: NAN the masked pixels
             if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskBad) {
-                image->data.F32[y][x] = 0.0;
-            } else {
-                float value = psImageUnbinPixel(x + 0.5, y + 0.5, modelImage, binning); // Background value
-                if (!isfinite(value)) {
-                    image->data.F32[y][x] = NAN;
-                    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskBad;
-                } else {
-                    image->data.F32[y][x] -= value;
-                }
-            }
+                image->data.F32[y][x] = NAN;
+		continue;
+            } 
+	    // special case 1: NAN & mask pixels without a valid background model
+	    float value = psImageUnbinPixel(x + 0.5, y + 0.5, modelImage, binning); // Background value
+	    if (!isfinite(value)) {
+	      image->data.F32[y][x] = NAN;
+	      mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskBad;
+	      continue;
+	    } 
+	    image->data.F32[y][x] -= value;
         }
     }
Index: /trunk/ppSub/src/ppSubFlagNeighbors.c
===================================================================
--- /trunk/ppSub/src/ppSubFlagNeighbors.c	(revision 29003)
+++ /trunk/ppSub/src/ppSubFlagNeighbors.c	(revision 29003)
@@ -0,0 +1,256 @@
+/** @file ppSubFlagNeighbors.c
+ *
+ *  @ingroup ppSub
+ *  @author EAM, IfA
+ *  @date $Date: 2009-02-18 00:31:20 $
+ *  Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+#include <psphot.h>
+
+#include "ppSub.h"
+
+// match the diff sources to the detections from the positive or negative input images
+// matchRef : false if this is the standard diff image, true if this is the inverse (in-ref vs ref-in)
+// this is needed to choose which value is positive and which is negative relative to the difference
+bool ppSubFlagNeighbors(pmConfig *config, pmFPAview *view, psArray *sources, bool matchRef) {
+
+    bool status;
+
+    psAssert(config, "Require config");
+    psAssert(view, "Require config");
+    psAssert (sources, "missing sources?");
+
+    // Look up recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSUB_RECIPE); // should this be psphot?
+    psAssert(recipe, "We checked this earlier, so it should be here.");
+
+    float radius = psMetadataLookupF32(&status, recipe, "INPUT.MATCH.RADIUS");
+    psAssert (status, "cannot find INPUT.MATCH.RADIUS");
+
+    float minSN  = psMetadataLookupF32(&status, recipe, "INPUT.MATCH.MIN.SN");
+    psAssert (status, "cannot find INPUT.MATCH.MIN.SN");
+
+    // Input sources
+    pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES");
+    pmDetections *inDetections = inSourceRO ? psMetadataLookupPtr(NULL, inSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Input sources
+    if (!inDetections) {
+        psWarning("Unable to filter detections based on image A.");
+    }
+
+    pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES");
+    pmDetections *refDetections = refSourceRO ? psMetadataLookupPtr(NULL, refSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Ref sources
+    if (!refDetections) {
+        psWarning("Unable to filter detections based on image B.");
+    }
+    if (!inDetections && !refDetections) {
+	return true;
+    }
+
+    psArray *inSources = inDetections->allSources;
+    psAssert (inSources, "missing inSources?");
+
+    psArray *refSources = refDetections->allSources;
+    psAssert (refSources, "missing refSources?");
+
+    ppSubSetSourceImageIDs(sources, 0);
+    ppSubSetSourceImageIDs(inSources, 1);
+    ppSubSetSourceImageIDs(refSources, 2);
+
+    // define the objects from the diff sources (1-to-1)
+    psArray *objects = psArrayAllocEmpty(sources->n);
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *src = sources->data[i]; 
+	pmPhotObj *obj = pmPhotObjAlloc();
+	pmPhotObjAddSource(obj, src);
+	psArrayAdd (objects, 100, obj);
+	psFree (obj);
+    }
+
+    // Match the sources to both inSources and refSources.  We want to find the closest source
+    // in the positive source lists to each source in the diff source list.
+    // XXX we need to either (a) apply S/N and other filters before we pass the in/ref sources
+    // to this function, or (b) add those filters in the function
+    ppSubMatchSources (objects, inSources, radius, minSN);
+    ppSubMatchSources (objects, refSources, radius, minSN);
+    
+    // now mark or flag the sources with matches that meet some criterion
+    for (int i = 0; i < objects->n; i++) {
+        pmPhotObj *obj = objects->data[i]; 
+
+	// the first entry should be from image 0 (the diff image).
+	// we should have only 1 of 4 possible results:
+	// a) only one entry -- the diff source (no positive matches)
+	// b) 2 matches -- diff and a source from image 1 
+	// c) 2 matches -- diff and a source from image 2
+	// d) 3 matches -- diff and a source from image 1 & 2
+
+	switch (obj->sources->n) {
+	  case 1: {
+	      pmSource *source = obj->sources->data[0];
+	      psAssert (source->imageID == 0, "error in pmPhotObj construction: first entry is not image ID 0");
+	      break;
+	  }
+	  case 2: {
+	      pmSource *source = obj->sources->data[0];
+	      psAssert (source->imageID == 0, "error in pmPhotObj construction: first entry not 0 (case 2)");
+
+	      pmSource *sourceM1 = obj->sources->data[1];
+	      psAssert ((sourceM1->imageID == 1) || (sourceM1->imageID == 2), "error in pmPhotObj construction (second entry is not from 1 or 2");
+
+	      // only one match.  set a flag?
+	      source->mode2 |= PM_SOURCE_MODE2_DIFF_WITH_SINGLE;
+	      if (source->diffStats) {
+		  // which is match need an xor here...
+		  bool positive = !matchRef && (sourceM1->imageID == 1);
+		  positive |= matchRef && (sourceM1->imageID == 2);
+		  float SN1 = isfinite(sourceM1->errMag) ? 1.0 / sourceM1->errMag : NAN;
+		  if (positive) {
+		      source->diffStats->SNp = SN1;
+		      source->diffStats->Rp  = hypot(source->peak->xf - sourceM1->peak->xf, source->peak->yf - sourceM1->peak->yf);
+		  } else {
+		      source->diffStats->SNm = SN1;
+		      source->diffStats->Rm  = hypot(source->peak->xf - sourceM1->peak->xf, source->peak->yf - sourceM1->peak->yf);
+		  }
+	      }
+	      break;
+	  }
+	  case 3: {
+	      pmSource *source = obj->sources->data[0];
+	      psAssert (source->imageID == 0, "error in pmPhotObj construction: first entry not 0 (case 3)");
+
+	      pmSource *sourceM1 = obj->sources->data[1];
+	      pmSource *sourceM2 = obj->sources->data[2];
+	      psAssert ((sourceM1->imageID == 1) || (sourceM1->imageID == 2), "error in pmPhotObj construction (case 3.1)");
+	      psAssert ((sourceM2->imageID == 1) || (sourceM2->imageID == 2), "error in pmPhotObj construction (case 3.2)");
+	      psAssert (sourceM1->imageID != sourceM2->imageID, "error in pmPhotObj construction (case 3.3)");
+
+	      source->mode2 |= PM_SOURCE_MODE2_DIFF_WITH_DOUBLE;
+	      if (source->diffStats) {
+		  // which is match need an xor here...
+		  bool positive = !matchRef && (sourceM1->imageID == 1);
+		  positive |= matchRef && (sourceM1->imageID == 2);
+		  float SN1 = isfinite(sourceM1->errMag) ? 1.0 / sourceM1->errMag : NAN;
+		  float SN2 = isfinite(sourceM2->errMag) ? 1.0 / sourceM2->errMag : NAN;
+		  if (positive) {
+		      source->diffStats->SNp = SN1;
+		      source->diffStats->Rp  = hypot(source->peak->xf - sourceM1->peak->xf, source->peak->yf - sourceM1->peak->yf);
+		      source->diffStats->SNm = SN2;
+		      source->diffStats->Rm  = hypot(source->peak->xf - sourceM2->peak->xf, source->peak->yf - sourceM2->peak->yf);
+		  } else {
+		      source->diffStats->SNp = SN2;
+		      source->diffStats->Rp  = hypot(source->peak->xf - sourceM2->peak->xf, source->peak->yf - sourceM2->peak->yf);
+		      source->diffStats->SNm = SN1;
+		      source->diffStats->Rm  = hypot(source->peak->xf - sourceM1->peak->xf, source->peak->yf - sourceM1->peak->yf);
+		  }
+	      }
+	      break;
+	  }
+	  default:
+	    psAbort ("error in pmPhotObj construction (unexpected matches)");
+	}
+    }
+    psFree (objects);
+    return true;
+}
+
+// XXX based on psphotMatchSourcesToObjects
+# define NEXT1 { i++; continue; } 
+# define NEXT2 { j++; continue; } 
+bool ppSubMatchSources (psArray *objects, psArray *sources, float RADIUS, float MIN_SN) { 
+ 
+    float dx, dy; 
+ 
+    float RADIUS2 = RADIUS*RADIUS;
+
+    // sort the source list by X 
+    sources = psArraySort (sources, pmSourceSortByX); 
+    objects = psArraySort (objects, pmPhotObjSortByX); 
+ 
+    // match sources to existing objects
+
+    psLogMsg ("psphot", PS_LOG_DETAIL, "attempt to match sources (%ld vs %ld) (SN > %f)", sources->n, objects->n, MIN_SN);
+
+    int nMatch = 0;
+    int i = 0;
+    int j = 0;
+    for (i = j = 0; (i < sources->n) && (j < objects->n); ) { 
+ 
+        pmSource  *src = sources->data[i]; 
+        pmPhotObj *obj = objects->data[j]; 
+ 
+        if (!src) NEXT1; 
+        if (!src->peak) NEXT1; 
+        if (!finite(src->peak->xf)) NEXT1; 
+        if (!finite(src->peak->yf)) NEXT1; 
+        if (!finite(src->peak->SN)) NEXT1; 
+        if (src->peak->SN < MIN_SN) NEXT1; 
+ 
+        if (!obj) NEXT2; 
+        if (!finite(obj->x)) NEXT2; 
+        if (!finite(obj->y)) NEXT2; 
+ 
+        dx = src->peak->xf - obj->x; 
+        if (dx < -1.02*RADIUS) NEXT1; 
+        if (dx > +1.02*RADIUS) NEXT2; 
+ 
+        // we are within match range, look for closest match to object j
+	int Imin = -1;
+	float Rmin = RADIUS2;
+        for (int I = i; (dx < +1.02*RADIUS) && (I < sources->n); I++) { 
+ 
+	    src = sources->data[I]; 
+	    
+	    if (!src) NEXT1; 
+	    if (!src->peak) NEXT1; 
+	    if (!finite(src->peak->xf)) NEXT1; 
+	    if (!finite(src->peak->yf)) NEXT1; 
+	    if (!finite(src->peak->SN)) NEXT1; 
+	    if (src->peak->SN < MIN_SN) NEXT1; 
+
+	    dx = src->peak->xf - obj->x; 
+            dy = src->peak->yf - obj->y; 
+ 
+            float dr = dx*dx + dy*dy; 
+            if (dr > RADIUS2) continue; 
+	    if (dr > Rmin) continue;
+	    Rmin = dr;
+	    Imin  = I;
+	    // fprintf (stderr, "j: %d, I: %d, Imin: %d, dr: %f, dx: %f, dy: %f\n", j, I, Imin, dr, dx, dy);
+	}
+
+	// no match, try next object
+	if (Imin == -1) {
+	    // fprintf (stderr, "*** missed j: %d, Imin: %d, obj x,y: %f, %f\n", j, Imin, obj->x, obj->y);
+	    j++;
+	    continue;
+	}
+	src = sources->data[Imin]; 
+
+	fprintf (stderr, "j: %d, Imin: %d, obj x,y: %f, %f  src x,y: %f, %f, SN: %f, ID: %d\n", j, Imin, obj->x, obj->y, src->peak->xf, src->peak->yf, src->peak->SN, src->id);
+
+	// add source to object
+	pmPhotObjAddSource (obj, src);
+	nMatch ++;
+        j++; 
+    } 
+
+    psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%d matches)", nMatch);
+    return true;
+} 
+
+bool ppSubSetSourceImageIDs (psArray *sources, int imageID) {
+
+    for (int i = 0; i < sources->n; i++) {
+	pmSource *source = sources->data[i];
+	source->imageID = imageID;
+    }
+    return true;
+}
Index: /trunk/ppSub/src/ppSubLoop.c
===================================================================
--- /trunk/ppSub/src/ppSubLoop.c	(revision 29002)
+++ /trunk/ppSub/src/ppSubLoop.c	(revision 29003)
@@ -17,6 +17,16 @@
 #include <pslib.h>
 #include <psmodules.h>
+#include <psphot.h>
 
 #include "ppSub.h"
+
+bool dumpout(pmConfig *config, char *name) 
+{
+  pmFPAview *view = ppSubViewReadout(); // View to readout
+  pmReadout *out = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT");
+  psphotSaveImage (NULL, out->image, name);
+  return true;
+}
+
 
 bool ppSubLoop(ppSubData *data)
@@ -105,4 +115,5 @@
         return false;
     }
+    // dumpout(config, "diff.1.fits");
 
     // Close convolved files
@@ -111,4 +122,5 @@
         return false;
     }
+    // dumpout(config, "diff.2a.fits");
 
     // Higher order background subtraction using psphot
@@ -117,4 +129,5 @@
         return false;
     }
+    // dumpout(config, "diff.2b.fits");
 
     // Perform Variance correction (rescale within a modest range)
@@ -123,4 +136,5 @@
         return false;
     }
+    // dumpout(config, "diff.2c.fits");
 
     if (data->quality) {
@@ -138,4 +152,5 @@
         return false;
     }
+    // dumpout(config, "diff.3.fits");
 
     if (!ppSubFilesIterateUp(config, PPSUB_FILES_PHOT_SUB)) {
@@ -157,4 +172,5 @@
       }
     }
+    // dumpout(config, "diff.4.fits");
     
     // generate the binned image used to write the jpeg
@@ -204,4 +220,5 @@
     } else {
         // Close subtraction files
+	// dumpout(config, "diff.5.fits");
         if (!ppSubFilesIterateUp(config, PPSUB_FILES_SUB)) {
             psError(PPSUB_ERR_IO, false, "Unable to close subtraction files.");
Index: /trunk/ppSub/src/ppSubReadoutPhotometry.c
===================================================================
--- /trunk/ppSub/src/ppSubReadoutPhotometry.c	(revision 29002)
+++ /trunk/ppSub/src/ppSubReadoutPhotometry.c	(revision 29003)
@@ -97,4 +97,10 @@
     psAssert (sources, "missing sources?");
 
+    // a likely source of false positives is poor subtractions.  this results in
+    // detections in the wings (or cores) of bright(er) stars found in both images.
+    // flag detections based on their distance from the bright(er) input sources.
+    bool matchRef = !strcasecmp(name, "PPSUB.INVERSE");
+    ppSubFlagNeighbors (config, view, sources, matchRef);
+
     if (data->stats) {
         bool mdok;
Index: /trunk/ppSub/src/ppSubReadoutSubtract.c
===================================================================
--- /trunk/ppSub/src/ppSubReadoutSubtract.c	(revision 29002)
+++ /trunk/ppSub/src/ppSubReadoutSubtract.c	(revision 29003)
@@ -51,9 +51,16 @@
     outRO->variance = (psImage*)psBinaryOp(outRO->variance, minuend->variance, "+", subtrahend->variance);
 
+    // NAN the masked pixels in the diff image (pixels masked in A are not yet NAN'ed in B)
+    psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config) | pmConfigMaskGet("BLANK", config); // Bits to mask in inputs
+    for (int iy = 0; iy < outRO->image->numRows; iy++) {
+	for (int ix = 0; ix < outRO->image->numCols; ix++) {
+	    if ((outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) == 0) continue;
+	    outRO->image->data.F32[iy][ix] = NAN;
+	}
+    }
+
     // Measure the variance scales
     psStats *varStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance images
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);           // Random number generator
-    psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config) |
-        pmConfigMaskGet("BLANK", config); // Bits to mask in inputs
     psImageBackground(varStats, NULL, minuend->variance, minuend->mask, maskVal, rng);
     float minuendVar = varStats->robustMedian; // Mean variance for minuend