Index: trunk/psModules/src/objects/pmSourceMoments.c
===================================================================
--- trunk/psModules/src/objects/pmSourceMoments.c	(revision 29044)
+++ trunk/psModules/src/objects/pmSourceMoments.c	(revision 29546)
@@ -64,4 +64,6 @@
 void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; }
 
+// if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid)
+
 bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal)
 {
@@ -71,30 +73,13 @@
     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
 
-    // use sky from moments if defined, 0.0 otherwise 
-
-    // XXX this value comes from the sky model at the source center, and tends to over-estimate
-    // the sky in the vicinity of bright sources.  we are better off assuming the model worked
-    // well:
+    // this function assumes the sky has been well-subtracted for the image
     psF32 sky = 0.0;
+
     if (source->moments == NULL) {
       source->moments = pmMomentsAlloc();
     }
-    // XXX if (source->moments == NULL) {
-    // XXX     source->moments = pmMomentsAlloc();
-    // XXX } else {
-    // XXX     sky = source->moments->Sky;
-    // XXX }
-
-    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
-    // Sum = SUM (z - sky)
-    // X1  = SUM (x - xc)*(z - sky)
-    // .. etc
-
-    psF32 peakPixel = -PS_MAX_F32;
-    psS32 numPixels = 0;
+
     psF32 Sum = 0.0;
     psF32 Var = 0.0;
-    psF32 X1 = 0.0;
-    psF32 Y1 = 0.0;
     psF32 R2 = PS_SQR(radius);
     psF32 minSN2 = PS_SQR(minSN);
@@ -109,107 +94,13 @@
     // (int) so they can be used in the image index below.
 
-    int xOff  = source->peak->x;
-    int yOff  = source->peak->y;
-    int xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage
-    int yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage
-
-    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
-    // psF32 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");
-
-    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
-
-        psF32 yDiff = row - yPeak;
-	if (fabs(yDiff) > radius) continue;
-
-        psF32 *vPix = source->pixels->data.F32[row];
-        psF32 *vWgt = source->variance->data.F32[row];
-        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->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;
-
-            psF32 xDiff = col - xPeak;
-	    if (fabs(xDiff) > radius) continue;
-
-            // radius is just a function of (xDiff, yDiff)
-            if (!VALID_RADIUS(xDiff, yDiff, R2)) continue;
-
-            psF32 pDiff = *vPix - sky;
-            psF32 wDiff = *vWgt;
-
-	    // skip pixels below specified significance level.  for a PSFs, this
-	    // over-weights the wings of bright stars compared to those of faint stars.
-	    // for the estimator used for extended source analysis (where the window
-	    // function is allowed to be arbitrarily large), we need to clip to avoid
-	    // negative second moments.
-            if (PS_SQR(pDiff) < minSN2*wDiff) continue; // 
-            if ((minSN > 0.0) && (pDiff < 0)) continue; // 
-
-	    // Apply a Gaussian window function.  Be careful with the window function.  S/N
-	    // weighting over weights the sky for faint sources
-	    if (sigma > 0.0) {
-		// XXX a lot of extra flops; can we pre-calculate?
-		psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
-		assert (z >= 0.0);
-		psF32 weight  = exp(-z);
-
-		wDiff *= weight;
-		pDiff *= weight;
-	    } 
-
-	    Var += wDiff;
-	    Sum += pDiff;
-
-	    psF32 xWght = xDiff * pDiff;
-	    psF32 yWght = yDiff * pDiff;
-
-	    X1  += xWght;
-	    Y1  += yWght;
-
-	    peakPixel = PS_MAX (*vPix, peakPixel);
-	    numPixels++;
-	}
-    }
-
-    // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
-    int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
-
-    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
-    if ((numPixels < minPixels) || (Sum <= 0)) {
-	psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
-	return (false);
-    }
-
-    // calculate the first moment.
-    float Mx = X1/Sum;
-    float My = Y1/Sum;
-    if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
-	psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
-	return (false);
-    }
-
-    psTrace ("psModules.objects", 5, "sky: %f  Mx: %f  My: %f  Sum: %f  X1: %f  Y1: %f  Npix: %d\n", sky, Mx, My, Sum, X1, Y1, numPixels);
-
-    // add back offset of peak in primary image
-    // also offset from pixel index to pixel coordinate
-    // (the calculation above uses pixel index instead of coordinate)
-    // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
-    source->moments->Mx = Mx + xOff + 0.5;
-    source->moments->My = My + yOff + 0.5;
-
-    source->moments->Sum = Sum;
-    source->moments->SN  = Sum / sqrt(Var);
-    source->moments->Peak = peakPixel;
-    source->moments->nPixels = numPixels;
+    // 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)) {
+	return false;
+    }
 
     // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
@@ -234,8 +125,8 @@
     Sum = 0.0;  // the second pass may include slightly different pixels, re-determine Sum
 
-    // center of mass in subimage.  Note: the calculation below uses pixel index, so we do not
-    // correct xCM, yCM to pixel coordinate here.
-    psF32 xCM = Mx + xPeak; // coord of peak in subimage
-    psF32 yCM = My + yPeak; // coord of peak in subimage
+    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
+    // xCM, yCM from pixel coords to pixel index here.
+    psF32 xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
+    psF32 yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
@@ -263,6 +154,5 @@
 	    // radius is just a function of (xDiff, yDiff)
 	    psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-	    psF32 r  = sqrt(r2);
-	    if (r > radius) continue;
+	    if (r2 > R2) continue;
 
 	    psF32 fDiff = *vPix - sky;
@@ -274,13 +164,10 @@
 	    // stars.
 	    if (PS_SQR(pDiff) < minSN2*wDiff) continue;
-	    // if (pDiff < 0) continue;
+	    if ((minSN > 0.0) && (pDiff < 0)) continue; // 
 
 	    // Apply a Gaussian window function.  Be careful with the window function.  S/N
 	    // weighting over weights the sky for faint sources
 	    if (sigma > 0.0) {
-		// XXX a lot of extra flops; can we do pre-calculate?
-		// XXX we were re-calculating r2 (maybe the compiler caught this?)
-		// psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
-		psF32 z  = r2 * rsigma2;
+		psF32 z = r2 * rsigma2;
 		assert (z >= 0.0);
 		psF32 weight  = exp(-z);
@@ -292,16 +179,9 @@
 	    Sum += pDiff;
 
-# if (1)
-# if (0)
-	    if (fDiff < 0) continue;
-# endif
+	    // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
+	    psF32 r = sqrt(r2);
 	    psF32 rf = r * fDiff;
 	    psF32 rh = sqrt(r) * fDiff;
 	    psF32 rs = fDiff;
-# else
-	    psF32 rf = r * pDiff;
-	    psF32 rh = sqrt(r) * pDiff;
-	    psF32 rs = pDiff;
-# endif
 
 	    psF32 x = xDiff * pDiff;
@@ -363,5 +243,4 @@
 
     // Calculate the Kron magnitude (make this block optional?)
-    // float radKron = 2.5*source->moments->Mrf;
     float radKinner = 1.0*source->moments->Mrf;
     float radKron   = 2.5*source->moments->Mrf;
@@ -397,5 +276,4 @@
 	    // radKron is just a function of (xDiff, yDiff)
 	    psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-	    psF32 r  = sqrt(r2);
 
 	    psF32 pDiff = *vPix - sky;
@@ -407,4 +285,5 @@
 	    if (PS_SQR(pDiff) < minSN2*wDiff) continue;
 
+	    psF32 r  = sqrt(r2);
 	    if (r < radKron) {
 		Sum += pDiff;
@@ -434,6 +313,146 @@
 
     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->flux, source->peak->SN, source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, numPixels);
+	     source->peak->xf, source->peak->yf, source->peak->flux, source->peak->SN, 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, psF32 radius, psF32 sigma, psF32 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
+
+    psF32 sky = 0.0;
+
+    psF32 peakPixel = -PS_MAX_F32;
+    psS32 numPixels = 0;
+    psF32 Sum = 0.0;
+    psF32 Var = 0.0;
+    psF32 X1 = 0.0;
+    psF32 Y1 = 0.0;
+    psF32 R2 = PS_SQR(radius);
+    psF32 minSN2 = PS_SQR(minSN);
+    psF32 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?)
+    // psF32 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++) {
+
+	psF32 yDiff = row + 0.5 - yPeak;
+	if (fabs(yDiff) > radius) continue;
+
+	psF32 *vPix = source->pixels->data.F32[row];
+	psF32 *vWgt = source->variance->data.F32[row];
+	psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->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;
+
+	    psF32 xDiff = col + 0.5 - xPeak;
+	    if (fabs(xDiff) > radius) continue;
+
+	    // radius is just a function of (xDiff, yDiff)
+	    psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+	    if (r2 > R2) continue;
+
+	    psF32 pDiff = *vPix - sky;
+	    psF32 wDiff = *vWgt;
+
+	    // skip pixels below specified significance level.  for a PSFs, this
+	    // over-weights the wings of bright stars compared to those of faint stars.
+	    // for the estimator used for extended source analysis (where the window
+	    // function is allowed to be arbitrarily large), we need to clip to avoid
+	    // negative second moments.
+	    if (PS_SQR(pDiff) < minSN2*wDiff) continue; // 
+	    if ((minSN > 0.0) && (pDiff < 0)) continue; // 
+
+	    // Apply a Gaussian window function.  Be careful with the window function.  S/N
+	    // weighting over weights the sky for faint sources
+	    if (sigma > 0.0) {
+		psF32 z  = r2*rsigma2;
+		assert (z >= 0.0);
+		psF32 weight  = exp(-z);
+
+		wDiff *= weight;
+		pDiff *= weight;
+	    } 
+
+	    Var += wDiff;
+	    Sum += pDiff;
+
+	    psF32 xWght = xDiff * pDiff;
+	    psF32 yWght = yDiff * pDiff;
+
+	    X1  += xWght;
+	    Y1  += yWght;
+
+	    peakPixel = PS_MAX (*vPix, peakPixel);
+	    numPixels++;
+	}
+    }
+
+    // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
+    int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
+
+    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
+    if ((numPixels < minPixels) || (Sum <= 0)) {
+	psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
+	return (false);
+    }
+
+    // calculate the first moment.
+    float Mx = X1/Sum;
+    float My = Y1/Sum;
+    if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
+	psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
+	return (false);
+    }
+
+    psTrace ("psModules.objects", 5, "sky: %f  Mx: %f  My: %f  Sum: %f  X1: %f  Y1: %f  Npix: %d\n", sky, Mx, My, Sum, X1, Y1, numPixels);
+
+    // add back offset of peak in primary image
+    // also offset from pixel index to pixel coordinate
+    // (the calculation above uses pixel index instead of coordinate)
+    // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
+
+    // we only update the centroid if the position is not supplied from elsewhere
+    bool skipCentroid = false;
+    skipCentroid |= (source->mode  & PM_SOURCE_MODE_EXTERNAL); // skip externally supplied positions
+    skipCentroid |= (source->mode2 & PM_SOURCE_MODE2_MATCHED); // skip sources defined by other image positions
+
+    if (skipCentroid) {
+	source->moments->Mx = source->peak->xf;
+	source->moments->My = source->peak->yf;
+    } else {
+	source->moments->Mx = Mx + xGuess;
+	source->moments->My = My + yGuess;
+    }
+
+    source->moments->Sum = Sum;
+    source->moments->SN  = Sum / sqrt(Var);
+    source->moments->Peak = peakPixel;
+    source->moments->nPixels = numPixels;
+
+    return true;
+}