Index: /branches/eam_branches/20100225/psphot/src/psphotSourceSize.c
===================================================================
--- /branches/eam_branches/20100225/psphot/src/psphotSourceSize.c	(revision 27269)
+++ /branches/eam_branches/20100225/psphot/src/psphotSourceSize.c	(revision 27270)
@@ -13,4 +13,5 @@
     float soft;
     int grow;
+    int xtest, ytest;
 } psphotSourceSizeOptions;
 
@@ -101,4 +102,9 @@
     assert (status);
 
+    // XXX recipe name is not great
+    options.xtest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.XTEST");
+    options.ytest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.YTEST");
+    assert (status);
+
     options.grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs
     if (!status || options.grow < 0) {
@@ -239,8 +245,11 @@
             continue;
         }
+	// psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
     }
 
     return true;
 }
+
+# define SIZE_SN_LIM 10
 
 bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
@@ -285,5 +294,5 @@
         }
 
-        // we are basically classifying by moments; use the default if not found
+        // we are basically classifying by moments
         psAssert (source->moments, "why is this source missing moments?");
         if (source->mode & noMoments) {
@@ -292,4 +301,5 @@
         }
 
+	// convert to Mmaj, Mmin:
         psF32 Mxx = source->moments->Mxx;
         psF32 Myy = source->moments->Myy;
@@ -315,18 +325,36 @@
         float apMag = -2.5*log10(source->moments->Sum);
         float dMag = source->psfMag - apMag;
-        float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
-
-        source->extNsigma = nSigma;
-        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+
+	// set nSigma to include both systematic and poisson error terms
+	// XXX the 'poisson error' contribution for size is probably wrong...
+        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
+	float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag);
+	float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag);
+
+	// partially-masked sources are more likely to be mis-measured PSFs
+	float sizeBias = 1.0;
+	if (source->pixWeight < 0.9) {
+	    sizeBias = 3.0;
+	} 
+
+	float minMxx = psfClump->X - sizeBias*options->nSigmaMoments*psfClump->dX;
+	float minMyy = psfClump->Y - sizeBias*options->nSigmaMoments*psfClump->dY;
+
+	// include MAG, MXX, and MYY?
+        source->extNsigma = nSigmaMAG;
+
+	// notes to clarify the source size classification rules:
+	// * a defect should be functionally equivalent to a cosmic ray
+	// * CR & defect should have a faintess limit (min S/N)
+	// * SAT stars should not be faint, but defects may?
 
         // Anything within this region is a probably PSF-like object. Saturated stars may land
         // in this region, but are detected elsewhere on the basis of their peak value.
-        bool isPSF = ((fabs(nSigma) < options->nSigmaApResid) &&
-		      (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) &&
-		      (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY));
+        bool isPSF = (fabs(nSigmaMAG) < options->nSigmaApResid) && (fabs(nSigmaMXX) < sizeBias*options->nSigmaMoments) && (fabs(nSigmaMYY) < sizeBias*options->nSigmaMoments);
         if (isPSF) {
-	  psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
-		  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
-		  options->nSigmaApResid,options->nSigmaMoments);
+	  psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
+		  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
+		  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
+	  source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
 	  Npsf ++;
 	  continue;
@@ -336,10 +364,13 @@
         // Defects may also be marked as SATSTAR -- XXX deactivate this flag?
         // XXX this rule is not great
-        if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
-
-	  psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
-		  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
-		  options->nSigmaApResid,options->nSigmaMoments);
+	// XXX only accept brightish detections as CRs
+	// (nSigmaMAG < -options->nSigmaApResid) || 
+	bool isCR = isCR = (source->errMag < 1.0 / SIZE_SN_LIM) && ((Mxx < minMxx) || (Myy < minMyy));
+	if (isCR) {
+	    psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g %f\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
+		    source->peak->xf,source->peak->yf,source->pixWeight,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
+		    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
             source->mode |= PM_SOURCE_MODE_DEFECT;
+	    source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
             Ncr ++;
             continue;
@@ -349,8 +380,8 @@
         // just large saturated regions.
         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
-	  psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
-		  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
-		  options->nSigmaApResid,options->nSigmaMoments);
-	  
+	    psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
+		    source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
+		    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
+	    source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
             Nsat ++;
             continue;
@@ -358,19 +389,20 @@
 
         // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
-        bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y));
+        bool isEXT = (nSigmaMAG > options->nSigmaApResid) || (Mxx > minMxx) || (Myy > minMyy);
         if (isEXT) {
-	  psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
-		  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
-		  options->nSigmaApResid,options->nSigmaMoments);
+	  psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
+		  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
+		  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
 	  
             source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
+	    source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
             Next ++;
             continue;
         }
-	psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
-		source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
-		options->nSigmaApResid,options->nSigmaMoments);
+	psTrace("psphotSourceClassRegion.MISS",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
+		source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
+		options->nSigmaApResid,sizeBias*options->nSigmaMoments);
 	
-        psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma);
+        psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigmaMAG);
         Nmiss ++;
     }
@@ -385,4 +417,7 @@
 bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
 
+    psTimerStart ("psphot.cr");
+
+    int nMasked = 0;
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
@@ -393,4 +428,14 @@
             continue;
         }
+
+        // only check candidates marked above
+        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) {
+            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
+            continue;
+        }
+
+        // skip unless this source is thought to be a cosmic ray.  flag the detection and mask the pixels
+	// XXX this may be degenerate with the above test
+	if (!(source->mode & PM_SOURCE_MODE_DEFECT)) continue;
 
         // Integer position of peak
@@ -402,17 +447,23 @@
             yPeak < 1 || yPeak > source->pixels->numRows - 2) {
             psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
-	    //	    psTrace("psphot.czw", 2, "Not calculating crNsigma due to edge\n");
             continue;
         }
 	
-        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
-	if (source->mode & PM_SOURCE_MODE_DEFECT) {
-	    // XXX this is running slowly and is too agressive, but it more-or-less works
-	    // psphotMaskCosmicRayCZW(readout, source, options->crMask);
-	    
-	    // XXX these are the old versions which we are not using
-	    // psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
-	    // psphotMaskCosmicRay_Old (source, maskVal, crMask);
-        }
+	// XXX for testing, only CRMASK a single source:
+	if (options->xtest && (fabs(source->peak->xf - options->xtest) > 5)) continue;
+	if (options->ytest && (fabs(source->peak->yf - options->ytest) > 5)) continue;
+
+        // replace object in image
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
+        }
+
+	// XXX this is running slowly and is too agressive, but it more-or-less works
+	psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf);
+	psphotMaskCosmicRay(readout, source, options->crMask);
+	nMasked ++;
+
+        // re-subtract the object, leave local sky
+        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
     }
     
@@ -430,10 +481,16 @@
     }
 
+    psLogMsg ("psphot.cr", PS_LOG_INFO, "mask CR: %d masked in %f sec\n", nMasked, psTimerMark ("psphot.cr"));
+
     // XXX test : save the mask image
-    if (0) {
+    if (1) {
 	psphotSaveImage (NULL, readout->mask,   "mask.fits");
     }
+
     return true;
 }
+
+# define LIMIT_XRANGE(X, IMAGE) { X = PS_MIN(PS_MAX(0, X), IMAGE->numCols); }
+# define LIMIT_YRANGE(Y, IMAGE) { Y = PS_MIN(PS_MAX(0, Y), IMAGE->numRows); }
 
 // Comments by CZW 20091209 : Mechanics of how to identify CR pixels taken from "Cosmic-Ray
@@ -448,27 +505,23 @@
     pmFootprint *footprint = peak->footprint;
 
-    int xm = footprint->bbox.x0;
-    int xM = footprint->bbox.x1;
-    int ym = footprint->bbox.y0;
-    int yM = footprint->bbox.y1;
-
-    if (xm < 0) { xm = 0; }
-    if (ym < 0) { ym = 0; }
-    if (xM > mask->numCols) { xM = mask->numCols; }
-    if (yM > mask->numRows) { yM = mask->numRows; }
-    int dx = xM - xm;
-    int dy = yM - ym;
-
     // Bounding boxes are inclusive of final pixel
-    // XXX ensure dx,dy do not become too large here
-    dx++;
-    dy++;
+    int xs = footprint->bbox.x0;
+    int xe = footprint->bbox.x1 + 1;
+    int ys = footprint->bbox.y0;
+    int ye = footprint->bbox.y1 + 1;
+
+    LIMIT_XRANGE(xs, mask);
+    LIMIT_XRANGE(xe, mask);
+    LIMIT_YRANGE(ys, mask);
+    LIMIT_YRANGE(ye, mask);
+
+    int dx = xe - xs;
+    int dy = ye - ys;
 
     psImage *image= readout->image;
     psImage *variance = readout->variance;
       
-    int binning = 1;
-    float sigma_thresh = 2.0;
-    int iteration = 0;
+    int binning = 2;
+    float sigma_thresh = 3.0;
     int max_iter = 5;
     float noise_factor = 5.0 / 4.0;  // Intrinsic to the Laplacian making noise spikes spikier.
@@ -476,4 +529,5 @@
     // Temporary images.
     psImage *mypix  = psImageAlloc(dx,dy,image->type.type);
+    psImage *myfix  = psImageAlloc(dx,dy,image->type.type);
     psImage *myvar  = psImageAlloc(dx,dy,image->type.type);
     psImage *binned = psImageAlloc(dx * binning,dy * binning,image->type.type);
@@ -482,10 +536,9 @@
     psImage *mymask = psImageAlloc(dx,dy,PS_TYPE_IMAGE_MASK);
       
-    int x,y;
     // Load my copy of things.
-    for (x = 0; x < dx; x++) {
-	for (y = 0; y < dy; y++) {
-	    psImageSet(mypix,x,y,psImageGet(image,x+xm,y+ym));
-	    psImageSet(myvar,x,y,psImageGet(variance,x+xm,y+ym));
+    for (int y = 0; y < dy; y++) {
+	for (int x = 0; x < dx; x++) {
+	    mypix->data.F32[y][x] = image->data.F32[y+ys][x+xs];
+	    myvar->data.F32[y][x] = variance->data.F32[y+ys][x+xs];
 	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0x00;
 	}
@@ -495,111 +548,126 @@
 	pmSpan *sp = footprint->spans->data[i];
 	for (int j = sp->x0; j <= sp->x1; j++) {
-	    y = sp->y - ym;
-	    x = j - xm;
+	    int y = sp->y - ys;
+	    int x = j - xs;
 	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01;
 	}
     }
 
-    int CRpix_count = 0;      
-    do {
-	CRpix_count = 0;
-	// Zero out my temp images.
-	for (x = 0; x < binning * dx; x++) {
-	    for (y = 0; y < binning * dy; y++) {
-		psImageSet(binned,x,y,0.0);
-		psImageSet(conved,x,y,0.0);
-		psImageSet(edges,x/binning,y/ binning,0.0);
+    int nCRpix = 1; // force at least one pass...
+    for (int iteration = 0; (iteration < max_iter) && (nCRpix > 0); iteration++) {
+	nCRpix = 0;
+	psImageInit (binned, 0.0);
+	psImageInit (conved, 0.0);
+	psImageInit (edges, 0.0);
+
+	// Make subsampled image. Maybe this should be called "unbinned" or something
+	for (int y = 0; y < binning * dy; y++) {
+	    int yraw = y / binning;
+	    for (int x = 0; x < binning * dx; x++) {
+		int xraw = x / binning;
+		binned->data.F32[y][x] = mypix->data.F32[yraw][xraw];
 	    }
-	}      
-	// Make subsampled image. Maybe this should be called "unbinned" or something
-	for (x = 0; x < binning * dx; x++) {
-	    for (y = 0; y < binning * dy; y++) {
-		psImageSet(binned,x,y,psImageGet(mypix,x / binning,y / binning));
+	}
+
+	// Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
+	for (int y = 1; y < binning * dy - 1; y++) {
+	    for (int x = 1; x < binning * dx - 1; x++) {
+		float value = binned->data.F32[y][x] - 0.25 * 
+		    (binned->data.F32[y+0][x-1] + binned->data.F32[y+0][x+1] + 
+		     binned->data.F32[y-1][x+0] + binned->data.F32[y+1][x+0]);
+		value = PS_MAX(0.0, value);
+
+		conved->data.F32[y][x] = value;
 	    }
 	}
-	// Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
-	for (x = 1; x < dx - 1; x++) {
-	    for (y = 1; y < dy - 1; y++) {
-		psImageSet(conved,x,y,psImageGet(binned,x,y) - 0.25 *
-			   (psImageGet(binned,x-1,y) + psImageGet(binned,x+1,y) +
-			    psImageGet(binned,x,y-1) + psImageGet(binned,x,y+1)));
-		if (psImageGet(conved,x,y) < 0.0) {
-		    psImageSet(conved,x,y,0.0);
+
+	// Create an edge map by rebinning
+	for (int y = 0; y < binning * dy; y++) {
+	    int yraw = y / binning;
+	    for (int x = 0; x < binning * dx; x++) {
+		int xraw = x / binning;
+		edges->data.F32[yraw][xraw] += conved->data.F32[y][x];
+	    }
+	}
+
+	// Modify my mask if we're above the significance threshold
+	for (int y = 0; y < dy; y++) {
+	    for (int x = 0; x < dx; x++) {
+		if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x01)) continue;
+
+		float noise = binning * sqrt(noise_factor * myvar->data.F32[y][x]);
+		float value = edges->data.F32[y][x] / noise;
+
+		if (value > sigma_thresh ) {
+		    mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x40;
+		    nCRpix ++;
 		}
 	    }
 	}
-	// Create an edge map by rebinning
-	for (x = 0; x < binning * dx; x++) {
-	    for (y = 0; y < binning * dy; y++) {
-		psImageSet(edges,x / binning, y / binning,
-			   psImageGet(edges, x / binning, y / binning) +
-			   psImageGet(conved,x,y));
-	    }
-	}
-	// Modify my mask if we're above the significance threshold
-	for (x = 0; x < dx; x++) {
-	    for (y = 0; y < dy; y++) {
-		if ( psImageGet(edges,x,y) / (binning * sqrt(noise_factor * psImageGet(myvar,x,y))) > sigma_thresh ) {
-		    if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x01) {
-			mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x40;
-			CRpix_count++;
-		    }
-		}
-	    }
-	}
+
+# if 1
+	psphotSaveImage (NULL, mypix,   "crmask.pix.fits");
+# endif
 
 	// "Repair" Masked pixels for the next round.
-	for (x = 1; x < dx - 1; x++) {
-	    for (y = 1; y < dy - 1; y++) {
-		if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
-		    psImageSet(mypix,x,y,0.25 *
-			       (psImageGet(mypix,x-1,y) + psImageGet(mypix,x+1,y) +
-				psImageGet(mypix,x,y-1) + psImageGet(mypix,x,y+1)));
-		}
-	    }
-	}
-	
-# if 0
- 	if ((psTraceGetLevel("psphot.czw") >= 2)&&(abs(xm - 2770) < 5)&&(abs(ym - 2581) < 5)&&(iteration == 0)) {
- 	  psTrace("psphot.czw",2,"TRACEMOTRON %d %d %d %d %d\n",xm,ym,dx,dy,iteration);
- 	  psphotSaveImage (NULL, mypix,   "czw.pix.fits");
- 	  psphotSaveImage (NULL, myvar,   "czw.var.fits");
- 	  psphotSaveImage (NULL, binned,  "czw.binn.fits");
- 	  psphotSaveImage (NULL, conved,  "czw.conv.fits");
- 	  psphotSaveImage (NULL, edges,   "czw.edge.fits");
- 	  psphotSaveImage (NULL, mymask,  "czw.mask.fits");
- 	}
-# endif
-
-	psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration,CRpix_count);
-	iteration++;
-    } while ((iteration < max_iter)&&(CRpix_count > 0));
-
-    // A solitary masked pixel is likely a lie. Remove those.
-    for (x = 0; x < dx; x++) {
-	for (y = 0; y < dy; y++) {
-	    if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
-		if ((x-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
+	for (int y = 1; y < dy - 1; y++) {
+	    for (int x = 1; x < dx - 1; x++) {
+		if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) {
+		    myfix->data.F32[y][x] = mypix->data.F32[y][x];
 		    continue;
 		}
-		if ((y-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
-		    continue;
-		}
-		if ((x+1 < dx)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
-		    continue;
-		}
-		if ((y+1 < dy)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
-		    continue;
-		}
-		mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
+		myfix->data.F32[y][x] = 0.25 * 
+		    (mypix->data.F32[y+0][x-1] + mypix->data.F32[y+0][x+1] +
+		     mypix->data.F32[y-1][x+0] + mypix->data.F32[y+1][x+0]);
 	    }
 	}
-    }
-
-    // transfer temporary mask to real mask
-    for (x = 0; x < dx; x++) {
-	for (y = 0; y < dy; y++) {
+	
+	// "Repair" Masked pixels for the next round.
+	for (int y = 1; y < dy - 1; y++) {
+	    for (int x = 1; x < dx - 1; x++) {
+		mypix->data.F32[y][x] = myfix->data.F32[y][x];
+	    }
+	}
+
+# if 1
+	fprintf (stderr, "CRMASK %d %d %d %d %d\n", xs, ys, dx, dy, iteration);
+	psphotSaveImage (NULL, mypix,   "crmask.fix.fits");
+	psphotSaveImage (NULL, myvar,   "crmask.var.fits");
+	psphotSaveImage (NULL, binned,  "crmask.binn.fits");
+	psphotSaveImage (NULL, conved,  "crmask.conv.fits");
+	psphotSaveImage (NULL, edges,   "crmask.edge.fits");
+	psphotSaveImage (NULL, mymask,  "crmask.mask.fits");
+# endif
+	psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration, nCRpix);
+    }
+
+    // A solitary masked pixel is likely a lie. Remove those
+    // XXX can't we use nCRpix == 1 to test for these?
+    for (int x = 0; x < dx; x++) {
+	for (int y = 0; y < dy; y++) {
+	    if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) continue;
+	    if ((x-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
+		continue;
+	    }
+	    if ((y-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
+		continue;
+	    }
+	    if ((x+1 < dx) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
+		continue;
+	    }
+	    if ((y+1 < dy) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
+		continue;
+	    }
+	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
+	}
+    }
+
+    // transfer temporary mask to real mask & count masked pixels
+    nCRpix = 0;
+    for (int x = 0; x < dx; x++) {
+	for (int y = 0; y < dy; y++) {
 	    if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
-		mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ym+mask->row0][x+xm+mask->col0] |= maskVal;
+		mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ys+mask->row0][x+xs+mask->col0] |= maskVal;
+		nCRpix ++;
 	    }
 	}
@@ -607,7 +675,10 @@
       
     // XXX if we decide this REALLY is a cosmic ray, set the CR_LIMIT bit
-    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+    if (nCRpix > 1) {
+	source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+    }
 
     psFree(mypix);
+    psFree(myfix);
     psFree(myvar);
     psFree(binned);