Index: /branches/eam_branches/20100225/psphot/doc/notes.20100131.txt
===================================================================
--- /branches/eam_branches/20100225/psphot/doc/notes.20100131.txt	(revision 27439)
+++ /branches/eam_branches/20100225/psphot/doc/notes.20100131.txt	(revision 27440)
@@ -1,2 +1,16 @@
+
+2010.03.24
+
+ I have finished another iteration on source size analysis.  I feel
+ fairly confident of the extended source and CR identifications.  They
+ are not perfect, but they are in pretty good shape.  The CR IDs use
+ the moments to find things smaller in one dimension than a PSF, while
+ the EXT is set for things larger than a PSF in one dimention. 
+
+ There still seems to be some trouble with SATSTARS : these are being
+ set for sources which are fainter than the SAT limit -- I'm not sure
+ if this is due to the quality of the PSF fit or if they really are
+ SAT, but have integrated fluxes lower than expected (underfitted in
+ PSF).
 
 20100131
Index: /branches/eam_branches/20100225/psphot/src/psphotSourceSize.c
===================================================================
--- /branches/eam_branches/20100225/psphot/src/psphotSourceSize.c	(revision 27439)
+++ /branches/eam_branches/20100225/psphot/src/psphotSourceSize.c	(revision 27440)
@@ -23,4 +23,5 @@
 bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal);
 bool psphotMaskCosmicRayFootprintCheck (psArray *sources);
+int  psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh);
 
 // we need to call this function after sources have been fitted to the PSF model and
@@ -252,5 +253,4 @@
 
 # define SIZE_SN_LIM 10
-
 bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
 
@@ -484,5 +484,5 @@
 
     // XXX test : save the mask image
-    if (1) {
+    if (0) {
 	psphotSaveImage (NULL, readout->mask,   "mask.fits");
     }
@@ -491,4 +491,5 @@
 }
 
+# define DUMPPICS 0
 # 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); }
@@ -524,6 +525,5 @@
     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.
+    int max_iter = 1; // XXX with isophot masking, we only want to do a single pass
 
     // Temporary images.
@@ -591,23 +591,17 @@
 	}
 
-	// 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 ++;
-		}
-	    }
-	}
-
-# if 1
+	// coordinate of peak in subimage pixels:
+	int xPeak = peak->x - xs;
+	int yPeak = peak->y - ys;
+
+	// Modify my mask if we're above the significance threshold, but only for connected pixels
+	nCRpix = psphotMaskCosmicRayConnected (xPeak, yPeak, mymask, myvar, edges, binning, sigma_thresh);
+	
+# if DUMPPICS
 	psphotSaveImage (NULL, mypix,   "crmask.pix.fits");
 # endif
-
+	
+// XXX do not repair the pixels in isophot version
+# if 0
 	// "Repair" Masked pixels for the next round.
 	for (int y = 1; y < dy - 1; y++) {
@@ -629,6 +623,7 @@
 	    }
 	}
-
-# if 1
+# endif
+
+# if DUMPPICS
 	fprintf (stderr, "CRMASK %d %d %d %d %d\n", xs, ys, dx, dy, iteration);
 	psphotSaveImage (NULL, mypix,   "crmask.fix.fits");
@@ -642,4 +637,5 @@
     }
 
+# if 0
     // A solitary masked pixel is likely a lie. Remove those
     // XXX can't we use nCRpix == 1 to test for these?
@@ -662,4 +658,5 @@
 	}
     }
+# endif
 
     // transfer temporary mask to real mask & count masked pixels
@@ -677,5 +674,7 @@
     if (nCRpix > 1) {
 	source->mode |= PM_SOURCE_MODE_CR_LIMIT;
-    }
+	source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+    }
+    // fprintf (stderr, "CRMASK %d %d %d %d %d\n", peak->x, peak->y, dx, dy, nCRpix);
 
     psFree(mypix);
@@ -753,4 +752,147 @@
 }
 
+# define VERBOSE 0
+int psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh) {
+    
+    int xLo, xRo;
+    int nCRpix = 0;
+
+    float noise_factor = 5.0 / 4.0;  // Intrinsic to the Laplacian making noise spikes spikier.
+    
+    // mark the pixels in this row to the left, then the right. stay within footprint
+    int xL = xPeak; // find the range of valid pixels in this row
+    int xR = xPeak;
+    for (int ix = xPeak; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix--) {
+	float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
+	float value = edges->data.F32[yPeak][ix] / noise;
+	if (value < sigma_thresh ) break;
+	mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
+	xL = ix;
+	nCRpix ++;
+	if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
+    }
+    for (int ix = xPeak; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix++) {
+	float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
+	float value = edges->data.F32[yPeak][ix] / noise;
+	if (value < sigma_thresh ) break;
+	mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
+	xR = ix;
+	nCRpix ++;
+	if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
+    }
+    // xL and xR mark the first and last valid pixel in the row
+
+    // for each of the neighboring rows, mark the high pixels if they touch the range xL to xR 
+    xLo = PS_MAX(xL - 1, 0);
+    xRo = PS_MIN(xR + 1, mymask->numCols);
+
+    // first go down:
+    for (int iy = yPeak - 1; iy >= 0; iy--) {
+
+	int xLn = -1;
+	int xRn = -1;
+	int newPix = 0;
+
+	// mark the pixels in the good range
+	for (int ix = xLo; ix < xRo; ix++) {
+	    if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
+ 	    float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+	    float value = edges->data.F32[iy][ix] / noise;
+	    if (value < sigma_thresh ) continue;
+	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+	    if (xLn == -1) xLn = ix; // first valid pixel in this row
+	    xRn = ix;		     // last valid pixel in this row
+	    nCRpix ++;
+	    newPix ++;
+	    if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+	}
+
+	// mark the pixels to the left of the good range
+	for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
+	    float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+	    float value = edges->data.F32[iy][ix] / noise;
+	    if (value < sigma_thresh ) break;
+	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+	    if (xRn == -1) xRn = ix; // last valid pixel in this row
+	    xLn = ix;
+	    nCRpix ++;
+	    newPix ++;
+	    if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+	}
+
+	// mark the pixels to the right of the good range
+	for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
+	    float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+	    float value = edges->data.F32[iy][ix] / noise;
+	    if (value < sigma_thresh ) break;
+	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+	    if (xLn == -1) xLn = ix; // first valid pixel in this row
+	    xRn = ix;
+	    nCRpix ++;
+	    newPix ++;
+	    if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+	}
+	if (newPix == 0) break;
+	xLo = PS_MAX(xLn - 1, 0);
+	xRo = PS_MIN(xRn + 1, mymask->numCols);
+    }
+
+    xLo = PS_MAX(xL - 1, 0);
+    xRo = PS_MIN(xR + 1, mymask->numCols);
+
+    // next go up:
+    for (int iy = yPeak + 1; iy < mymask->numRows; iy++) {
+
+	int xLn = -1;
+	int xRn = -1;
+	int newPix = 0;
+
+	// mark the pixels in the good range
+	for (int ix = xLo; ix < xRo; ix++) {
+	    if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
+ 	    float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+	    float value = edges->data.F32[iy][ix] / noise;
+	    if (value < sigma_thresh ) continue;
+	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+	    if (xLn == -1) xLn = ix; // first valid pixel in this row
+	    xRn = ix;		     // last valid pixel in this row
+	    nCRpix ++;
+	    newPix ++;
+	    if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+	}
+
+	// mark the pixels to the left of the good range
+	for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
+	    float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+	    float value = edges->data.F32[iy][ix] / noise;
+	    if (value < sigma_thresh ) break;
+	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+	    if (xRn == -1) xRn = ix; // last valid pixel in this row
+	    xLn = ix;
+	    nCRpix ++;
+	    newPix ++;
+	    if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+	}
+
+	// mark the pixels to the right of the good range
+	for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
+	    float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+	    float value = edges->data.F32[iy][ix] / noise;
+	    if (value < sigma_thresh ) break;
+	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+	    if (xLn == -1) xLn = ix; // first valid pixel in this row
+	    xRn = ix;
+	    nCRpix ++;
+	    newPix ++;
+	    if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+	}
+	if (newPix == 0) break;
+	xLo = PS_MAX(xLn - 1, 0);
+	xRo = PS_MIN(xRn + 1, mymask->numCols);
+    }
+
+    return nCRpix;
+}
+
 bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
 
Index: /branches/eam_branches/20100225/psphot/src/psphotVisual.c
===================================================================
--- /branches/eam_branches/20100225/psphot/src/psphotVisual.c	(revision 27439)
+++ /branches/eam_branches/20100225/psphot/src/psphotVisual.c	(revision 27440)
@@ -180,5 +180,5 @@
 bool psphotVisualShowImage (pmReadout *readout) {
 
-    // if (!pmVisualIsVisual()) return true;
+    if (!pmVisualIsVisual()) return true;
 
     int kapa = psphotKapaChannel (1);
@@ -1514,5 +1514,4 @@
         if (source == NULL) continue;
 
-        // if (source->type != type) continue;
         if (mode) {
             if (keep) {
@@ -1551,5 +1550,5 @@
 bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) {
 
-    // if (!pmVisualIsVisual()) return true;
+    if (!pmVisualIsVisual()) return true;
 
     int myKapa = psphotKapaChannel (1);
@@ -1578,5 +1577,5 @@
     KapaSection section;
 
-    // if (!pmVisualIsVisual()) return true;
+    if (!pmVisualIsVisual()) return true;
 
     int myKapa = psphotKapaChannel (2);