Index: branches/eam_branches/ipp-20100823/psModules/src/objects/pmSource.h
===================================================================
--- branches/eam_branches/ipp-20100823/psModules/src/objects/pmSource.h	(revision 29472)
+++ branches/eam_branches/ipp-20100823/psModules/src/objects/pmSource.h	(revision 29473)
@@ -249,5 +249,6 @@
   psF32 sigma, 
   psF32 minSN, 
-  psImageMaskType maskVal);
+  psImageMaskType maskVal, 
+  float xGuess, float yGuess);
 
 pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source);
Index: branches/eam_branches/ipp-20100823/psModules/src/objects/pmSourceMoments.c
===================================================================
--- branches/eam_branches/ipp-20100823/psModules/src/objects/pmSourceMoments.c	(revision 29472)
+++ branches/eam_branches/ipp-20100823/psModules/src/objects/pmSourceMoments.c	(revision 29473)
@@ -64,6 +64,4 @@
 void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; }
 
-bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal);
-
 // if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid)
 
@@ -96,5 +94,11 @@
     // (int) so they can be used in the image index below.
 
-    if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal)) {
+    // 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;
     }
@@ -314,5 +318,5 @@
 }
 
-bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal) { 
+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)
@@ -333,8 +337,6 @@
     psF32 rsigma2 = 0.5 / PS_SQR(sigma);
 
-    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
+    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?)
@@ -344,7 +346,12 @@
     // 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 - yPeak;
+	psF32 yDiff = row + 0.5 - yPeak;
 	if (fabs(yDiff) > radius) continue;
 
@@ -363,5 +370,5 @@
 	    if (isnan(*vPix)) continue;
 
-	    psF32 xDiff = col - xPeak;
+	    psF32 xDiff = col + 0.5 - xPeak;
 	    if (fabs(xDiff) > radius) continue;
 
@@ -439,6 +446,6 @@
 	source->moments->My = source->peak->yf;
     } else {
-	source->moments->Mx = Mx + xOff + 0.5;
-	source->moments->My = My + yOff + 0.5;
+	source->moments->Mx = Mx + xGuess;
+	source->moments->My = My + yGuess;
     }