Index: /trunk/psModules/src/imcombine/pmSubtractionStamps.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtractionStamps.c	(revision 19200)
+++ /trunk/psModules/src/imcombine/pmSubtractionStamps.c	(revision 19201)
@@ -24,4 +24,7 @@
 
 #define STAMP_LIST_BUFFER 20            // Number of stamps to add to list at a time
+
+#define SOURCE_MASK (PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_SATURATED | PM_SOURCE_MODE_DEFECT | \
+                     PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_FAIL) // Mask for bad sources
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -582,20 +585,25 @@
     // Let pmSubtractionStampsSet take care of the rest of the assertions
 
-    int numSources = sources->n;          // Number of stars
-
-    psVector *x = psVectorAlloc(numSources, PS_TYPE_F32); // x coordinates
-    psVector *y = psVectorAlloc(numSources, PS_TYPE_F32); // y coordinates
-    psVector *flux = psVectorAlloc(numSources, PS_TYPE_F32); // Fluxes
-
+    int numIn = sources->n;             // Number of sources in input list
+
+    psVector *x = psVectorAlloc(numIn, PS_TYPE_F32); // x coordinates
+    psVector *y = psVectorAlloc(numIn, PS_TYPE_F32); // y coordinates
+    psVector *flux = psVectorAlloc(numIn, PS_TYPE_F32); // Fluxes
+
+    int numOut = 0;                     // Number of sources in output list
     for (int i = 0; i < numSources; i++) {
         pmSource *source = sources->data[i]; // Source of interest
+        if (source->mode & SOURCE_MASK) {
+            continue;
+        }
         if (source->modelPSF) {
-            x->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_XPOS];
-            y->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+            x->data.F32[numOut] = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+            y->data.F32[numOut] = source->modelPSF->params->data.F32[PM_PAR_YPOS];
         } else {
-            x->data.F32[i] = source->peak->xf;
-            y->data.F32[i] = source->peak->yf;
-        }
-        flux->data.F32[i] = powf(10.0, -0.4 * source->psfMag);
+            x->data.F32[numOut] = source->peak->xf;
+            y->data.F32[numOut] = source->peak->yf;
+        }
+        flux->data.F32[numOut] = powf(10.0, -0.4 * source->psfMag);
+        numOut++;
     }