Index: trunk/psphot/src/psphotKronMasked.c
===================================================================
--- trunk/psphot/src/psphotKronMasked.c	(revision 31673)
+++ trunk/psphot/src/psphotKronMasked.c	(revision 32348)
@@ -1,4 +1,8 @@
 # include "psphotInternal.h"
-
+# ifndef ROUND
+# define ROUND(X) ((int) ((X) + 0.5*SIGN(X)))
+# endif
+
+bool psphotKronMask (pmSource *source, psImage *kronMask, psImageMaskType markVal);
 bool psphotKronMag (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal);
 
@@ -45,4 +49,8 @@
 }
 
+int psphotKapaChannel (int channel);
+bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel);
+bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
+
 bool psphotKronMaskedReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) {
 
@@ -72,4 +80,6 @@
     }
 
+    float EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS");
+
     // bit-masks to test for good/bad pixels
     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
@@ -83,9 +93,14 @@
     maskVal |= markVal;
 
+    // psphotSaveImage (NULL, readout->mask, "kron.unmasked.fits");
+
+    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY;
+
+    // XXX tmp visualization
+    // int kapa = psphotKapaChannel (1);
+
     // generate the mask image: increment counter for every source overlapping the pixel
-    psImage *kronMask = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_U8);
+    psImage *kronMask = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
     psImageInit (kronMask, 0);
-    int Nx = kronMask->numCols;
-    int Ny = kronMask->numRows;
     for (int i = 0; i < sources->n; i++) {
 
@@ -96,45 +111,59 @@
         if (!source->moments) continue;
 
-        // replace object in image
-        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
-            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-        }
-
-	psphotKronMag (source, RADIUS, MIN_KRON_RADIUS, maskVal);
-
-        // re-subtract the object, leave local sky
-        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-
-	continue;
-
-	// XXX skip this code
-	// generate the pixel masks
-	// int Xo = source->moments->Mx;
-	// int Yo = source->moments->My;
-	float dX = source->moments->Mx - source->peak->xf;
-	float dY = source->moments->My - source->peak->yf;
-	float dR = hypot(dX, dY);
-	
-	float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
-	float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
-	
-	int Kr = 2.5*source->moments->Mrf;
-	int Kr2 = Kr*Kr;
-	
-	for (int iy = Yo - Kr; iy < Yo + Kr + 1; iy++) {
-	    if (iy < 0) continue;
-	    if (iy >= Ny) continue;
-	    for (int ix = Xo - Kr; ix < Xo + Kr + 1; ix++) {
-		if (ix < 0) continue;
-		if (ix >= Nx) continue;
-
-		if (PS_SQR(ix - Xo) + PS_SQR(iy - Yo) > Kr2) continue;
-		if (kronMask->data.U8[iy][ix] < 0xff) {
-		    kronMask->data.U8[iy][ix] ++;
-		}
-	    }
-	}
-    }		
+	// replace object in image
+	if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+	    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+	}
+
+	// save the PSF-based values
+	// source->moments->KronFluxPSF    = source->moments->KronFlux;
+	// source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
+	// source->moments->KronRadiusPSF  = source->moments->Mrf;
+
+	// iterate to the window radius
+	float windowRadius = RADIUS;
+	for (int j = 0; j < 4; j++) {
+	    // XXX use some S/N criterion to limit?
+	    // if (source->moments->KronFlux / source->moments->KronFluxErr > 10.0) {
+	    // windowRadius = PS_MAX (RADIUS, 10*source->moments->Mrf);
+	    // }
+
+	    // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
+	    pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
+
+	    // mask the pixels not contained by the footprint
+	    // if (psphotMaskFootprint (readout, source, markVal)) {
+	    // 	// psphotVisualShowMask (kapa, source->maskObj, "mask", 2);
+	    // 	// psphotVisualRangeImage (kapa, source->pixels, "image", 1, -50, 300);
+	    // 	// fprintf (stderr, "masked\n");
+	    // }	    
+
+	    // mask the pixels associated with near neighbors
+	    if (psphotKronMask (source, kronMask, markVal)) {
+		// psphotVisualShowMask (kapa, source->maskObj, "mask", 2);
+		// psphotVisualRangeImage (kapa, source->pixels, "image", 1, -50, 300);
+	     	// fprintf (stderr, "masked\n");
+	    }	    
+	    // this function populates moments->Mrf,KronFlux,KronFluxErr
+	    // XXX what about KronInner, KronOuter, etc?
+	    psphotKronMag (source, windowRadius, MIN_KRON_RADIUS, maskVal);
+	    windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
+	    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+	}
+        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
+	float kmag = -2.5*log10(source->moments->KronFlux);
+	if (source->psfMag - kmag > 0.25) {
+	    // psphotVisualShowMask (kapa, source->maskObj, "mask", 1);
+	    // psphotVisualRangeImage (kapa, source->pixels, "image", 0, -50, 300);
+	    // fprintf (stderr, "psf: %f, kron: %f, dmag: %f\n", source->psfMag, kmag, source->psfMag - kmag);
+	    // fprintf (stderr, "continue\n");
+	}
+	// re-subtract the object, leave local sky
+	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+    }
+
+    // psphotSaveImage (NULL, readout->mask, "kron.masked.fits");
     // psphotSaveImage (NULL, kronMask, "kronmask.fits");
+
     psLogMsg ("psphot.kron", PS_LOG_DETAIL, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
 
@@ -142,4 +171,6 @@
     return true;
 }
+
+# define WEIGHTED 0
 
 bool psphotKronMag (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal) {
@@ -150,7 +181,9 @@
     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
 
-    psF32 Sum = 0.0;
-    psF32 Var = 0.0;
     psF32 R2 = PS_SQR(radius);
+
+# if (WEIGHTED)
+    float rsigma2 = 0.5 / PS_SQR(radius / 4.0); // use a Gaussian window with sigma = R_window / 2
+# endif
 
     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
@@ -168,5 +201,5 @@
     psF32 RS = 0.0;
 
-# if (1)
+    // the peak position is less accurate but less subject to extreme deviations
     float dX = source->moments->Mx - source->peak->xf;
     float dY = source->moments->My - source->peak->yf;
@@ -174,8 +207,4 @@
     float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
     float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
-# else
-    float Xo = source->moments->Mx;
-    float Yo = source->moments->My;
-# endif
 
     // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
@@ -210,9 +239,14 @@
 	    if (r2 > R2) continue;
 
+# if (WEIGHTED)
+	    float z = r2 * rsigma2;
+	    assert (z >= 0.0);
+	    float weight  = exp(-z);
+	    psF32 pDiff = *vPix * weight;
+# else
 	    psF32 pDiff = *vPix;
-
-	    Sum += pDiff;
-
-	    // Kron Flux uses the 1st radial moment (NOT Gaussian windowed)
+# endif
+
+	    // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
 	    psF32 rf = pDiff * sqrt(r2);
 	    psF32 rs = pDiff;
@@ -234,5 +268,6 @@
 
     int nKronPix = 0;
-    Sum = Var = 0.0;
+    float Sum = 0.0;
+    float Var = 0.0;
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
@@ -263,6 +298,14 @@
 	    if (r2 > radKron2) continue;
 
+# if (WEIGHTED)
+	    float z = r2 * rsigma2;
+	    assert (z >= 0.0);
+	    float weight  = exp(-z);
+	    psF32 pDiff = *vPix * weight;
+	    psF32 wDiff = *vWgt * weight;
+# else
 	    psF32 pDiff = *vPix;
 	    psF32 wDiff = *vWgt;
+# endif
 
 	    Sum += pDiff;
@@ -272,13 +315,173 @@
     }
 
-    source->moments->Mrh = Mrf;
-
-    // XXX for a test, save the old values here:
-    source->moments->KronFouter    = source->moments->KronFlux;
-    source->moments->KronCoreErr = source->moments->KronFluxErr;
-
-    source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
-    source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
+    source->moments->Mrf = Mrf;
+    source->moments->KronFlux    = Sum;
+    source->moments->KronFluxErr = sqrt(Var);
+    // source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
+    // source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
 
     return true;
 }
+
+enum {Y_P, Y_M, X_P, X_M};
+
+bool psphotKronMask (pmSource *source, psImage *kronMask, psImageMaskType markVal) {
+
+    int X, Y, e, dXs, dYs;
+
+    if (!source->peak) return false;
+    if (!source->peak->saddlePoints) return false;
+    if (!source->peak->saddlePoints->n) return false;
+
+    // return true;
+
+    int xMax = source->maskObj->numCols;
+    int yMax = source->maskObj->numRows;
+    // int xOff = source->maskObj->col0;
+    // int yOff = source->maskObj->row0;
+
+    int Xo = ROUND(source->peak->xf) - source->maskObj->col0;
+    int Yo = ROUND(source->peak->yf) - source->maskObj->row0;
+
+    for (int i = 0; i < source->peak->saddlePoints->n; i++) {
+
+	psVector *saddlePoint = source->peak->saddlePoints->data[i];
+	int Xs = ROUND(saddlePoint->data.S32[0]) - source->maskObj->col0;
+	int Ys = ROUND(saddlePoint->data.S32[1]) - source->maskObj->row0;
+	
+	// fprintf (stderr, "mask %d,%d @ %d,%d (%d,%d)\n", Xo, Yo, Xs, Ys, saddlePoint->data.S32[0], saddlePoint->data.S32[1]);
+
+	// We want to mask the pixels between the edge of the image and the line perpendicular
+	// to the connecting line.  Call dX,dY = (Xs-Xo,Ys-Yo).  There are 4 cases: 
+	// +y : dY > 0, fabs(dY) > fabs(dX)
+	// -y : dY < 0, fabs(dY) > fabs(dX)
+	// +x : dX > 0, fabs(dX) > fabs(dY)
+	// -x : dX < 0, fabs(dX) > fabs(dY)
+
+	int dX = Xs - Xo;
+	int dY = Ys - Yo;
+
+	// +y : dY > 0, fabs(dY) > fabs(dX)
+	// -y : dY < 0, fabs(dY) > fabs(dX)
+	if (fabs(dY) > fabs(dX)) {
+	    // points to right of saddle
+	    Y = Ys;
+	    e = 0;
+	    dXs = (dY > 0) ? +dY : -dY;  // this forces dXs > 0
+	    dYs = (dY > 0) ? -dX : +dX;
+	    for (X = Xs; X < xMax; X++) {
+		if (dY > 0) {
+		    for (int Ym = Y; Ym < yMax; Ym++) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+			// kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+		    }
+		} else {
+		    for (int Ym = Y; Ym >= 0; Ym--) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+			// kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+		    }
+		}
+		e += dYs;
+		int e2 = 2 * e;
+		if (e2 > dXs) {
+		    Y ++;
+		    e -= dXs;
+		}
+		if (e2 < -dXs) {
+		    Y --;
+		    e += dXs;
+		}
+		if (Y >= yMax) break;
+		if (Y < 0) break;
+	    }
+	    // points to left of saddle
+	    Y = Ys;
+	    e = 0;
+	    for (X = Xs; X >= 0; X--) {
+		if (dY > 0) {
+		    for (int Ym = Y; Ym < yMax; Ym++) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+			// kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+		    }
+		} else {
+		    for (int Ym = Y; Ym >= 0; Ym--) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+			// kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+		    }
+		}
+		e -= dYs;
+		int e2 = 2 * e;
+		if (e2 > dXs) {
+		    Y ++;
+		    e -= dXs;
+		}
+		if (e2 < -dXs) {
+		    Y --;
+		    e += dXs;
+		}
+		if (Y >= yMax) break;
+		if (Y < 0) break;
+	    }
+	} else {
+	    // points to right of saddle
+	    X = Xs;
+	    e = 0;
+	    dXs = (dX > 0) ? -dY : +dY;
+	    dYs = (dX > 0) ? +dX : -dX;  // this forces dYs > 0
+	    for (Y = Ys; Y < yMax; Y++) {
+		if (dX > 0) {
+		    for (int Xm = X; Xm < xMax; Xm++) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+			// kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+		    }
+		} else {
+		    for (int Xm = X; Xm >= 0; Xm--) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+			// kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+		    }
+		}
+		e += dXs;
+		int e2 = 2 * e;
+		if (e2 > dYs) {
+		    X ++;
+		    e -= dYs;
+		}
+		if (e2 < -dYs) {
+		    X --;
+		    e += dYs;
+		}
+		if (X >= yMax) break;
+		if (X < 0) break;
+	    }
+	    // points to left of saddle
+	    X = Xs;
+	    e = 0;
+	    for (Y = Ys; Y >= 0; Y--) {
+		if (dX > 0) {
+		    for (int Xm = X; Xm < xMax; Xm++) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+			// kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+		    }
+		} else {
+		    for (int Xm = X; Xm >= 0; Xm--) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+			// kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+		    }
+		}
+		e -= dXs;
+		int e2 = 2 * e;
+		if (e2 > dYs) {
+		    X ++;
+		    e -= dYs;
+		}
+		if (e2 < -dYs) {
+		    X --;
+		    e += dYs;
+		}
+		if (X >= yMax) break;
+		if (X < 0) break;
+	    }
+	}
+    }
+    return true;
+}