Index: trunk/psModules/src/objects/pmSourceMoments.c
===================================================================
--- trunk/psModules/src/objects/pmSourceMoments.c	(revision 31670)
+++ trunk/psModules/src/objects/pmSourceMoments.c	(revision 32347)
@@ -96,11 +96,18 @@
     // (int) so they can be used in the image index below.
 
-    // do 2 passes : the first pass should use a somewhat smaller radius and no sigma window to 
-    // get an unbiased (but probably noisy) centroid
-    if (!pmSourceMomentsGetCentroid (source, 0.75*radius, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
-	return false;
-    }
-    // second pass applies the Gaussian window and uses the centroid from the first pass
-    if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal, source->moments->Mx, source->moments->My)) {
+    // XXX // do 2 passes : the first pass should use a somewhat smaller radius and no sigma window to 
+    // XXX // get an unbiased (but probably noisy) centroid
+    // XXX if (!pmSourceMomentsGetCentroid (source, 0.75*radius, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
+    // XXX 	return false;
+    // XXX }
+    // XXX // second pass applies the Gaussian window and uses the centroid from the first pass
+    // XXX if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal, source->moments->Mx, source->moments->My)) {
+    // XXX 	return false;
+    // XXX }
+
+    // If we use a large radius for the centroid, it will be biased by any neighbors.  The flux
+    // of any object drops pretty quickly outside 1-2 sigmas.  (The exception is bright
+    // saturated stars, for which we need to use a very large radius here)
+    if (!pmSourceMomentsGetCentroid (source, 1.5*sigma, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
 	return false;
     }
@@ -108,11 +115,4 @@
     // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
     // Xn  = SUM (x - xc)^n * (z - sky)
-
-    float RFW = 0.0;
-    float RHW = 0.0;
-
-    float RF = 0.0;
-    float RH = 0.0;
-    float RS = 0.0;
     float XX = 0.0;
     float XY = 0.0;
@@ -143,4 +143,5 @@
     float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
 
+    // calculate the higher-order moments using Xo,Yo
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
@@ -194,12 +195,5 @@
 	    Sum += pDiff;
 
-	    // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
 	    float r = sqrt(r2);
-	    float rf = r * fDiff;
-	    float rh = sqrt(r) * fDiff;
-	    float rs = fDiff;
-
-	    float rfw = r * pDiff;
-	    float rhw = sqrt(r) * pDiff;
 
 	    float x = xDiff * pDiff;
@@ -221,11 +215,4 @@
 	    float yyyy = yDiff * yyy / r2;
 
-	    RF  += rf;
-	    RH  += rh;
-	    RS  += rs;
-
-	    RFW  += rfw;
-	    RHW  += rhw;
-
 	    XX  += xx;
 	    XY  += xy;
@@ -242,10 +229,22 @@
 	    XYYY  += xyyy;
 	    YYYY  += yyyy;
+
+	    // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
+	    // XXX float r = sqrt(r2);
+	    // XXX float rf = r * fDiff;
+	    // XXX float rh = sqrt(r) * fDiff;
+	    // XXX float rs = fDiff;
+	    // XXX 
+	    // XXX float rfw = r * pDiff;
+	    // XXX float rhw = sqrt(r) * pDiff;
+	    // XXX 
+	    // XXX RF  += rf;
+	    // XXX RH  += rh;
+	    // XXX RS  += rs;
+	    // XXX 
+	    // XXX RFW  += rfw;
+	    // XXX RHW  += rhw;
 	}
     }
-
-    source->moments->Mrf = RF/RS;
-    source->moments->Mrh = RH/RS;
-
     source->moments->Mxx = XX/Sum;
     source->moments->Mxy = XY/Sum;
@@ -262,4 +261,64 @@
     source->moments->Mxyyy = XYYY/Sum;
     source->moments->Myyyy = YYYY/Sum;
+
+    // *** now calculate the 1st radial moment (for kron flux) -- symmetrical averaging
+
+    float **vPix = source->pixels->data.F32;
+    float **vWgt = source->variance->data.F32;
+    psImageMaskType  **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
+
+    float RF = 0.0;
+    float RH = 0.0;
+    float RS = 0.0;
+
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+	float yDiff = row - yCM;
+	if (fabs(yDiff) > radius) continue;
+
+	// coordinate of mirror pixel
+	int yFlip = yCM - yDiff;
+	if (yFlip < 0) continue;
+	if (yFlip >= source->pixels->numRows) continue;
+
+	for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+	    // check mask and value for this pixel
+	    if (vMsk && (vMsk[row][col] & maskVal)) continue;
+	    if (isnan(vPix[row][col])) continue;
+
+	    float xDiff = col - xCM;
+	    if (fabs(xDiff) > radius) continue;
+
+	    // coordinate of mirror pixel
+	    int xFlip = xCM - xDiff;
+	    if (xFlip < 0) continue;
+	    if (xFlip >= source->pixels->numCols) continue;
+
+	    // check mask and value for mirror pixel
+	    if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
+	    if (isnan(vPix[yFlip][xFlip])) continue;
+
+	    // radius is just a function of (xDiff, yDiff)
+	    float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+	    if (r2 > R2) continue;
+
+	    float fDiff1 = vPix[row][col] - sky;
+	    float fDiff2 = vPix[yFlip][xFlip] - sky;
+	    float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
+
+	    // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
+	    float r = sqrt(r2);
+	    float rf = r * pDiff;
+	    float rh = sqrt(r) * pDiff;
+	    float rs = 0.5 * (fDiff1 + fDiff2);
+
+	    RF  += rf;
+	    RH  += rh;
+	    RS  += rs;
+	}
+    }
+
+    source->moments->Mrf = RF/RS;
+    source->moments->Mrh = RH/RS;
 
     // if Mrf (first radial moment) is very small, we are getting into low-significance
@@ -270,4 +329,6 @@
 	kronRefRadius = MIN(radius, kronRefRadius);
     }
+
+    // *** now calculate the kron flux values using the 1st radial moment
 
     float radKinner = 1.0*kronRefRadius;
@@ -283,14 +344,130 @@
     float SumOuter = 0.0;
 
+    // calculate the Kron flux, and related fluxes (NO symmetrical averaging)
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
-
+	
 	float yDiff = row - yCM;
 	if (fabs(yDiff) > radKouter) continue;
+	
+	for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+	    // check mask and value for this pixel
+	    if (vMsk && (vMsk[row][col] & maskVal)) continue;
+	    if (isnan(vPix[row][col])) continue;
+	    
+	    float xDiff = col - xCM;
+	    if (fabs(xDiff) > radKouter) continue;
+	    
+	    // radKron is just a function of (xDiff, yDiff)
+	    float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+
+	    float fDiff1 = vPix[row][col] - sky;
+	    float pDiff = fDiff1;
+	    float wDiff = vWgt[row][col];
+				    
+	    // skip pixels below specified significance level.  this is allowed, but should be
+	    // avoided -- the over-weights the wings of bright stars compared to those of faint
+	    // stars.
+	    if (PS_SQR(pDiff) < minSN2*wDiff) continue;
+	    
+	    float r  = sqrt(r2);
+	    if (r < radKron) {
+		Sum += pDiff;
+		Var += wDiff;
+		nKronPix ++;
+		// if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
+	    }
+
+	    // use sigma (fixed by psf) not a radKron based value
+	    if (r < sigma) {
+		SumCore += pDiff;
+		VarCore += wDiff;
+		nCorePix ++;
+	    }
+
+	    if ((r > radKinner) && (r < radKron)) {
+		SumInner += pDiff;
+		nInner ++;
+	    }
+	    if ((r > radKron)  && (r < radKouter)) {
+		SumOuter += pDiff;
+		nOuter ++;
+	    }
+	}
+    }
+    // *** should I rescale these fluxes by pi R^2 / nNpix?
+    // XXX source->moments->KronCore    = SumCore       * M_PI * PS_SQR(sigma) / nCorePix;
+    // XXX source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix;
+    // XXX source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
+    // XXX source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
+    // XXX source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
+    // XXX source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
+
+    source->moments->KronCore    = SumCore;
+    source->moments->KronCoreErr = sqrt(VarCore);
+    source->moments->KronFlux    = Sum;
+    source->moments->KronFluxErr = sqrt(Var);
+    source->moments->KronFinner = SumInner;
+    source->moments->KronFouter = SumOuter;
+
+    // XXX not sure I should save this here...
+    source->moments->KronFluxPSF    = source->moments->KronFlux;
+    source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
+    source->moments->KronRadiusPSF  = source->moments->Mrf;
+
+    psTrace ("psModules.objects", 4, "Mrf: %f  KronFlux: %f  Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",
+	     source->moments->Mrf,   source->moments->KronFlux, 
+	     source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
+	     source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
+	     source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
+
+    psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f  My: %f  Sum: %f  Mxx: %f  Mxy: %f  Myy: %f  sky: %f  Npix: %d\n",
+	     source->peak->xf, source->peak->yf, source->peak->rawFlux, sqrt(source->peak->detValue), source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, source->moments->nPixels);
+
+    return(true);
+}
+
+bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { 
+
+    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
+    // Sum = SUM (z - sky)
+    // X1  = SUM (x - xc)*(z - sky)
+    // .. etc
+
+    float sky = 0.0;
+
+    float peakPixel = -PS_MAX_F32;
+    psS32 numPixels = 0;
+    float Sum = 0.0;
+    float Var = 0.0;
+    float X1 = 0.0;
+    float Y1 = 0.0;
+    float R2 = PS_SQR(radius);
+    float minSN2 = PS_SQR(minSN);
+    float rsigma2 = 0.5 / PS_SQR(sigma);
+
+    float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
+    float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
+
+    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
+    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
+    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
+    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
+    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
+
+    // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
+    // not depend on the fractional pixel location of the source.  However, the aperture
+    // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
+    // position of the expected centroid
+
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+	float yDiff = row + 0.5 - yPeak;
+	if (fabs(yDiff) > radius) continue;
 
 	float *vPix = source->pixels->data.F32[row];
 	float *vWgt = source->variance->data.F32[row];
 
-	psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
-	// psImageMaskType  *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
+	psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+	// psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
 
 	for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
@@ -304,117 +481,4 @@
 	    if (isnan(*vPix)) continue;
 
-	    float xDiff = col - xCM;
-	    if (fabs(xDiff) > radKouter) continue;
-
-	    // radKron is just a function of (xDiff, yDiff)
-	    float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-
-	    float pDiff = *vPix - sky;
-	    float wDiff = *vWgt;
-
-	    // skip pixels below specified significance level.  this is allowed, but should be
-	    // avoided -- the over-weights the wings of bright stars compared to those of faint
-	    // stars.
-	    if (PS_SQR(pDiff) < minSN2*wDiff) continue;
-
-	    float r  = sqrt(r2);
-	    if (r < radKron) {
-		Sum += pDiff;
-		Var += wDiff;
-		nKronPix ++;
-		// if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
-	    }
-
-	    // use sigma (fixed by psf) not a radKron based value
-	    if (r < sigma) {
-		SumCore += pDiff;
-		VarCore += wDiff;
-		nCorePix ++;
-	    }
-
-	    if ((r > radKinner) && (r < radKron)) {
-		SumInner += pDiff;
-		nInner ++;
-	    }
-	    if ((r > radKron)  && (r < radKouter)) {
-		SumOuter += pDiff;
-		nOuter ++;
-	    }
-	}
-    }
-    // *** should I rescale these fluxes by pi R^2 / nNpix?
-    source->moments->KronCore    = SumCore       * M_PI * PS_SQR(sigma) / nCorePix;
-    source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix;
-    source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
-    source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
-    source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
-    source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
-
-    psTrace ("psModules.objects", 4, "Mrf: %f  KronFlux: %f  Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",
-	     source->moments->Mrf,   source->moments->KronFlux, 
-	     source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
-	     source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
-	     source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
-
-    psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f  My: %f  Sum: %f  Mxx: %f  Mxy: %f  Myy: %f  sky: %f  Npix: %d\n",
-	     source->peak->xf, source->peak->yf, source->peak->rawFlux, sqrt(source->peak->detValue), source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, source->moments->nPixels);
-
-    return(true);
-}
-
-bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { 
-
-    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
-    // Sum = SUM (z - sky)
-    // X1  = SUM (x - xc)*(z - sky)
-    // .. etc
-
-    float sky = 0.0;
-
-    float peakPixel = -PS_MAX_F32;
-    psS32 numPixels = 0;
-    float Sum = 0.0;
-    float Var = 0.0;
-    float X1 = 0.0;
-    float Y1 = 0.0;
-    float R2 = PS_SQR(radius);
-    float minSN2 = PS_SQR(minSN);
-    float rsigma2 = 0.5 / PS_SQR(sigma);
-
-    float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
-    float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
-
-    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
-    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
-    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
-    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
-    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
-
-    // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
-    // not depend on the fractional pixel location of the source.  However, the aperture
-    // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
-    // position of the expected centroid
-
-    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
-
-	float yDiff = row + 0.5 - yPeak;
-	if (fabs(yDiff) > radius) continue;
-
-	float *vPix = source->pixels->data.F32[row];
-	float *vWgt = source->variance->data.F32[row];
-
-	psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
-	// psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
-
-	for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
-	    if (vMsk) {
-		if (*vMsk & maskVal) {
-		    vMsk++;
-		    continue;
-		}
-		vMsk++;
-	    }
-	    if (isnan(*vPix)) continue;
-
 	    float xDiff = col + 0.5 - xPeak;
 	    if (fabs(xDiff) > radius) continue;