mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests 1 Actions 3 Packages Projects Releases 72 Wiki Activity

integration: the observed centroid is the signal's, not the disk's

Browse Source 
The centroid was a first moment of the RAW counts over the signal disk, so it weighted signal plus
background. The background is flat over a disk centred on the PREDICTION, which makes its own centroid
the prediction exactly: it adds nothing to the displacement and everything to the denominator, and the
measured offset comes out shrunk by I/(I + n*bkg).

That factor is worst where the background dominates the signal, which is at high resolution - so the
one consumer of this quantity, the geometry post-refinement, fits the beam centre and the detector
distance on displacements that are systematically too small, by a factor that varies with resolution.
An estimator whose bias depends on the very coordinate it is correcting.

Subtracting a flat pedestal from a first moment is exact, and the background is not known until the
ring has been read, so the positions of the pixels behind the intensity sum are accumulated alongside
it and the correction is applied afterwards: sum(x*(px-bkg)) = sum(x*px) - bkg*sum(x). Both engines,
identically. Where nothing rises above background there is no signal centroid to compute and the raw
one is kept.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
This commit is contained in:
leonarski_f
2026-09-01 19:55:47 +02:00
co-authored by Claude Opus 5
parent 4026ebc6ab
commit ecb0571d8d
2 changed files with 41 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp
+22 -2
View File
@@ -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;
image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu
+19 -4
View File
@@ -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;
}