Index: /trunk/psphot/src/psphotSourceSize.c
===================================================================
--- /trunk/psphot/src/psphotSourceSize.c	(revision 18834)
+++ /trunk/psphot/src/psphotSourceSize.c	(revision 18835)
@@ -10,9 +10,16 @@
 // deviation from the psf model at the r = FWHM/2 position
 
-bool psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe, long first) {
+bool psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first) {
 
     bool status;
 
     psTimerStart ("psphot");
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    assert (maskVal);
+
+    // bit to mask the cosmic-ray pixels
+    psMaskType crMask  = pmConfigMaskGet("CR", config); // Mask value for cosmic rays
 
     float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT");
@@ -129,6 +136,70 @@
 	}
 
+	// this source is thought to be a cosmic ray.  flag the detection and mask the pixels
 	if (source->crNsigma > CR_NSIGMA_LIMIT) {
-	  source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+	    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+	    pmPeak *peak = source->peak;
+	    psAssert (peak, "NULL peak");
+
+	    // replace the source flux
+	    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+	    source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
+
+	    psImage *mask   = source->maskView;
+	    psImage *pixels = source->pixels;
+	    psImage *weight = source->weight;
+
+	    # define SN_LIMIT 5.0
+
+	    int xo = peak->x - pixels->col0;
+	    int yo = peak->y - pixels->row0;
+
+	    // mark the pixels in this row to the left, then the right
+	    for (int ix = xo; ix >= 0; ix--) {
+		float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]);
+		if (SN > SN_LIMIT) {
+		    mask->data.U8[yo][ix] |= crMask;
+		}
+	    }
+	    for (int ix = xo + 1; ix < pixels->numCols; ix++) {
+		float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]);
+		if (SN > SN_LIMIT) {
+		    mask->data.U8[yo][ix] |= crMask;
+		}
+	    }
+
+	    // for each of the neighboring rows, mark the high pixels if they have a marked neighbor 
+	    // first go up:
+	    for (int iy = yo; iy >= 0; iy--) {
+		// mark the pixels in this row to the left, then the right
+		for (int ix = 0; ix < pixels->numCols; ix++) {
+		    float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
+		    if (SN < SN_LIMIT) continue;
+		    
+		    bool valid = false;
+		    valid |= (mask->data.U8[iy+1][ix] & crMask);
+		    valid |= (ix > 0) ? (mask->data.U8[iy+1][ix-1] & crMask) : 0;
+		    valid |= (ix <= mask->numCols) ? (mask->data.U8[iy+1][ix+1] & crMask) : 0;
+
+		    if (!valid) continue;
+		    mask->data.U8[iy][ix] |= crMask;
+		}
+	    }
+	    // next go down:
+	    for (int iy = yo+1; iy < pixels->numRows; iy++) {
+		// mark the pixels in this row to the left, then the right
+		for (int ix = 0; ix < pixels->numCols; ix++) {
+		    float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
+		    if (SN < SN_LIMIT) continue;
+		    
+		    bool valid = false;
+		    valid |= (mask->data.U8[iy-1][ix] & crMask);
+		    valid |= (ix > 0) ? (mask->data.U8[iy-1][ix-1] & crMask) : 0;
+		    valid |= (ix <= mask->numCols) ? (mask->data.U8[iy-1][ix+1] & crMask) : 0;
+
+		    if (!valid) continue;
+		    mask->data.U8[iy][ix] |= crMask;
+		}
+	    }
 	}
     }