diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp
index 015053f6b..339182e65 100644
--- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp
+++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp
@@ -164,6 +164,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
         out.k_bin = BraggStencilKernelIndex(st, n_kern);
 
         int64_t I_sum = 0, I_sum_x = 0, I_sum_y = 0, n_inner = 0, n_inner_valid = 0;
+        int64_t x_sum = 0, y_sum = 0;   // positions of the pixels behind I_sum, for the centroid
         int n_disk = 0, n_own = 0;   // pixels in the signal disk, and how many are this reflection's
         double bkg_sum = 0.0;
         int n_bkg = 0;
@@ -192,6 +193,10 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                     I_sum += px;
                     I_sum_x += static_cast<int64_t>(x) * px;
                     I_sum_y += static_cast<int64_t>(y) * px;
+                    // Position sums over the very same pixels, so the background can be taken out of
+                    // the centroid below. It is not known yet - the ring is read in this same loop.
+                    x_sum += x;
+                    y_sum += y;
                     ++n_inner_valid;
                 } else if (d.inner >= r2_sq && d.outer < r3_sq) {
                     if (refl_mask[y * W + x]) { ++n_bkg_neighbour; continue; }
@@ -283,8 +288,23 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
             out.sigma = 1.0;
             if (I_sum > 0) {
                 out.sigma = std::max(out.sigma, std::sqrt(static_cast<double>(I_sum) + var_bkg_term));
-                out.obs_x = static_cast<double>(I_sum_x) / static_cast<double>(I_sum);
-                out.obs_y = static_cast<double>(I_sum_y) / static_cast<double>(I_sum);
+                // The centroid of the SIGNAL, not of the disk. Weighting by the raw counts weights by
+                // signal plus background, and the background is flat over a disk centred on the
+                // PREDICTION - so its own centroid is the prediction exactly, and it pulls the answer
+                // there. The measured displacement comes out shrunk by I/(I + n*bkg), which is worst
+                // where the background dominates, i.e. at high resolution. Post-refinement fits the
+                // beam and the distance on these centroids, so it under-corrects by that factor and
+                // does so resolution-dependently. Subtracting a flat pedestal from a first moment is
+                // exact: sum(x*(px-bkg)) = I_sum_x - bkg*x_sum over the same pixels.
+                const double net = static_cast<double>(I_sum) - n_inner_valid * out.bkg;
+                if (net > 0.0) {
+                    out.obs_x = (static_cast<double>(I_sum_x) - out.bkg * x_sum) / net;
+                    out.obs_y = (static_cast<double>(I_sum_y) - out.bkg * y_sum) / net;
+                } else {
+                    // Nothing above background to take a centroid of; the raw one is all there is.
+                    out.obs_x = static_cast<double>(I_sum_x) / static_cast<double>(I_sum);
+                    out.obs_y = static_cast<double>(I_sum_y) / static_cast<double>(I_sum);
+                }
                 // A disk with a hole in it gives a centroid pulled away from the hole, and the hole
                 // sits at a fixed place on the detector - post-refinement would read that as geometry.
                 out.has_obs = full;
diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu
index 2cc96b2a2..c224c2870 100644
--- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu
+++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu
@@ -163,6 +163,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     if (i >= n) return;
 
     __shared__ unsigned long long s_Isum, s_Ix, s_Iy;
+    __shared__ unsigned long long s_x, s_y;   // positions behind s_Isum, to take the background out of the centroid
     __shared__ int s_ninner, s_ninner_valid, s_nbkg, s_ndisk, s_nown;
     __shared__ int s_nbkg_nb;   // ring pixels a NEIGHBOUR's signal region holds; see the CPU engine
     __shared__ unsigned long long s_bkgsum;
@@ -187,7 +188,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     if (threadIdx.x == 0) {
         const float rx = px_x[i] - p.beam_x, ry = px_y[i] - p.beam_y;
         s_r0 = sqrtf(rx * rx + ry * ry);
-        s_Isum = 0; s_Ix = 0; s_Iy = 0;
+        s_Isum = 0; s_Ix = 0; s_Iy = 0; s_x = 0; s_y = 0;
         s_ninner = 0; s_ninner_valid = 0; s_nbkg = 0; s_bkgsum = 0;
         s_ndisk = 0; s_nown = 0; s_nbkg_nb = 0;
     }
@@ -205,6 +206,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     const int area = (bw > 0 && bh > 0) ? bw * bh : 0;
 
     long long l_Isum = 0, l_Ix = 0, l_Iy = 0;
+    long long l_x = 0, l_y = 0;   // positions behind l_Isum, for the background-free centroid
     int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0, l_nbnb = 0;
     long long l_bkg = 0;
     for (int t = threadIdx.x; t < area; t += blockDim.x) {
@@ -222,7 +224,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
             }
             ++l_ni;
             const int32_t px = img[y * p.W + x];
-            if (valid(px)) { l_Isum += px; l_Ix += (long long) x * px; l_Iy += (long long) y * px; ++l_niv; }
+            if (valid(px)) { l_Isum += px; l_Ix += (long long) x * px; l_Iy += (long long) y * px;
+                             l_x += x; l_y += y; ++l_niv; }
         } else if (d.inner >= p.r2_sq && d.outer < p.r3_sq) {
             if (mask[y * p.W + x]) { ++l_nbnb; continue; }
             const int32_t px = img[y * p.W + x];
@@ -233,6 +236,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     WARP_ATOMIC_ADD(s_Isum, (unsigned long long) l_Isum);
     WARP_ATOMIC_ADD(s_Ix, (unsigned long long) l_Ix);
     WARP_ATOMIC_ADD(s_Iy, (unsigned long long) l_Iy);
+    WARP_ATOMIC_ADD(s_x, (unsigned long long) l_x);
+    WARP_ATOMIC_ADD(s_y, (unsigned long long) l_y);
     WARP_ATOMIC_ADD(s_ninner, l_ni);
     WARP_ATOMIC_ADD(s_ninner_valid, l_niv);
     WARP_ATOMIC_ADD(s_nbkg, l_nb);
@@ -380,8 +385,18 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     uint8_t hasobs = 0; double ox = 0.0, oy = 0.0;
     if (Isum > 0) {
         sigma = fmax(sigma, sqrt((double) Isum + var_bkg_term));
-        ox = (double) (long long) s_Ix / (double) Isum;
-        oy = (double) (long long) s_Iy / (double) Isum;
+        // The centroid of the SIGNAL. Weighting by raw counts weights by signal plus background,
+        // and the background is flat over a disk centred on the PREDICTION, so its centroid is the
+        // prediction and it pulls the answer there - the displacement comes out shrunk by
+        // I/(I + n*bkg), worst where the background dominates. See the CPU engine.
+        const double net = (double) Isum - (double) s_ninner_valid * bkg;
+        if (net > 0.0) {
+            ox = ((double) (long long) s_Ix - bkg * (double) (long long) s_x) / net;
+            oy = ((double) (long long) s_Iy - bkg * (double) (long long) s_y) / net;
+        } else {
+            ox = (double) (long long) s_Ix / (double) Isum;
+            oy = (double) (long long) s_Iy / (double) Isum;
+        }
         // A disk with a hole gives a centroid pulled away from it; see the CPU engine.
         hasobs = s_full ? 1 : 0;
     }