Index: trunk/psphot/src/psphotSourceSize.c
===================================================================
--- trunk/psphot/src/psphotSourceSize.c	(revision 20878)
+++ trunk/psphot/src/psphotSourceSize.c	(revision 20938)
@@ -4,4 +4,7 @@
 static float psphotModelContour(const psImage *image, const psImage *weight, const psImage *mask,
                                 psMaskType maskVal, const pmModel *model, float Ro);
+
+bool psphotMaskCosmicRay_Old (pmSource *source, psMaskType maskVal, psMaskType crMask);
+bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psMaskType maskVal, psMaskType crMask);
 
 // we need to call this function after sources have been fitted to the PSF model and
@@ -157,69 +160,6 @@
         // 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;
-            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;
-
-            // XXX This should be a recipe variable
-            # 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 = PS_MIN(yo, mask->numRows-2); 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 = PS_MIN(yo+1, mask->numRows-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;
-                }
-            }
+	    // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask);
+	    psphotMaskCosmicRay_Old (source, maskVal, crMask);
         }
     }
@@ -242,5 +182,5 @@
 
     psphotVisualPlotSourceSize (sources);
-    psphotVisualShowSourceSize (sources);
+    psphotVisualShowSourceSize (readout, sources);
 
     return true;
@@ -319,2 +259,108 @@
     return nSigma;
 }
+
+bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psMaskType maskVal, psMaskType crMask) {
+
+    // replace the source flux
+    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+    source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
+
+    // flag this as a CR
+    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+    pmPeak *peak = source->peak;
+    psAssert (peak, "NULL peak");
+
+    // grab the matching footprint
+    pmFootprint *footprint = peak->footprint;
+    if (!footprint) {
+	// if we have not footprint, use the old code to mask by isophot
+	psphotMaskCosmicRay_Old (source, maskVal, crMask);
+	return true;
+    }
+
+    if (!footprint->spans) {
+	// if we have not footprint, use the old code to mask by isophot
+	psphotMaskCosmicRay_Old (source, maskVal, crMask);
+	return true;
+    }
+
+    // mask all of the pixels covered by the spans of the footprint
+    for (int j = 1; j < footprint->spans->n; j++) {
+	pmSpan *span1 = footprint->spans->data[j];
+
+	int iy = span1->y;
+	int xs = span1->x0;
+	int xe = span1->x1;
+	
+	for (int ix = xs; ix < xe; ix++) {
+	    mask->data.U8[iy][ix] |= crMask;
+	}
+    }
+    return true;
+}
+
+bool psphotMaskCosmicRay_Old (pmSource *source, psMaskType maskVal, psMaskType crMask) {
+
+    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+    pmPeak *peak = source->peak;
+    psAssert (peak, "NULL peak");
+
+    psImage *mask   = source->maskView;
+    psImage *pixels = source->pixels;
+    psImage *weight = source->weight;
+
+    // XXX This should be a recipe variable
+# 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 = PS_MIN(yo, mask->numRows-2); 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 = PS_MIN(yo+1, mask->numRows-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;
+	}
+    }
+    return true;
+}