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

Bragg integration: separate the three things a bandwidth used to switch

Browse Source 
Setting a bandwidth flipped three unrelated switches at once: it changed the profile's
radial capture term, it moved the width measurement from the signal disk to the whole
fit grid, and it silently overrode the background clip and trim, so --background-clip
under --bandwidth was ignored - the two runs were bit-identical.

The width measurement was the damaging one. The fit grid is an azimuthally averaged
stack, so its second moment is sigma_r^2 + sigma_t^2 and the radial smear of a
bandwidth leaked into the tangential model - a tangential width of 3.04 px against a
1.06 px truth, inflating the effective background pixel count where the weak signal is.
The result was a step rather than a slope: on genuinely monochromatic data, declaring a
0.2% bandwidth cost ISa 28.4 -> 22.2.

Measure the two widths separately, accumulated in each spot's own radial/tangential
frame over the signal disk, from the signed profile cells - away from the peak a
learned cell is background noise centred on zero, so the signed sum is unbiased, while
clamping it at zero turns that noise into a pedestal the r^2 weight reads as width. The
radial term is then the measured excess or the analytic floor, whichever is larger.

With the two widths separated there is nothing left for the broadband switch to select,
so it is gone - which is the proof the three were independent. The background clip and
trim now come from the settings in every case; the tuned 3-sigma broadband default
moves to the rugnux front end, which is the only place that knows whether the user gave
a value.

Monochromatic data: declaring a 0.2% bandwidth now costs ISa 28.4 -> 27.9 rather than
22.2, and forcing the old 3-sigma clip in the new build reproduces the good result, so
none of the step came from the clip. On large-bandwidth data CC1/2 improves in 8 of 10
shells. Across 12 monochromatic crystals the space groups are unchanged and CC1/2 moves
by at most 0.2 points.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-08-09 21:08:29 +02:00
co-authored by Claude Opus 5
parent 3d3fb0e58b
commit 1239c49731
9 changed files with 174 additions and 124 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp
+51 -46
View File
@@ -262,13 +262,25 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
if (rough[i].ok) rough[i].shell = shell_of(predicted[i].d);
// --- Learn the profile per shell (+ global) from the strong spots. ---
// Two things are learned. The empirical profile is the average grid in the DETECTOR frame, which is
// where it is applied. The width is a pair of second moments taken in each spot's OWN radial /
// tangential frame: a grid stacked in the detector frame is azimuthally averaged, so its <r^2> is
// sigma_r^2 + sigma_t^2 with no way back, and a radially smeared spot reads as a wide TANGENTIAL
// one. Rotating each contribution into the spot's frame keeps the two apart.
struct Moments { double rad = 0.0, tan = 0.0, w = 0.0; };
struct Sigma2 { double rad = 1.0, tan = 1.0; };
std::vector<std::vector<double>> shell_grid(N_SHELL, std::vector<double>(GG, 0.0));
std::vector<Moments> shell_mom(N_SHELL);
std::vector<int> shell_n(N_SHELL, 0);
std::vector<double> global_grid(GG, 0.0);
Moments global_mom;
int global_n = 0;
for (size_t i = 0; i < npredicted; ++i) {
const auto &rh = rough[i];
if (!rh.ok || !rh.strong || rh.I <= 0.0) continue;
const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y;
const double Rpx = std::hypot(rx, ry);
const double ux = Rpx > 1e-6 ? rx / Rpx : 1.0, uy = Rpx > 1e-6 ? ry / Rpx : 0.0;
for (int dy = -R; dy <= R; ++dy)
for (int dx = -R; dx <= R; ++dx) {
const int x = rh.cx + dx, y = rh.cy + dy;
@@ -278,65 +290,54 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
const double v = (static_cast<double>(px) - rh.bkg) / rh.I;
shell_grid[rh.shell][grid_idx(dx, dy)] += v;
global_grid[grid_idx(dx, dy)] += v;
if (dx * dx + dy * dy >= r1_sq) continue;
const double rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux;
shell_mom[rh.shell].rad += v * rad * rad;
shell_mom[rh.shell].tan += v * tn * tn;
shell_mom[rh.shell].w += v;
global_mom.rad += v * rad * rad;
global_mom.tan += v * tn * tn;
global_mom.w += v;
}
++shell_n[rh.shell];
++global_n;
}
// Isotropic width (2nd moment) of a learned grid: over the r1 disk (monochromatic) or the full
// grid (broadband); <r^2> = 2 sigma^2 in 2D. The cells are signed: away from the peak a learned
// cell is pure background noise centred on zero, and clamping it at zero turns that noise into a
// positive pedestal spread over the whole domain, which the r^2 weight then reads as extra width.
auto measure_sigma2 = [&](const std::vector<double> &grid) {
double m2 = 0.0, m2w = 0.0;
for (int dy = -R; dy <= R; ++dy)
for (int dx = -R; dx <= R; ++dx) {
if (!broadband && dx * dx + dy * dy >= r1_sq) continue;
const double g = grid[grid_idx(dx, dy)];
m2 += g * (dx * dx + dy * dy);
m2w += g;
}
return m2w > 0.0 ? std::max(0.25, (m2 / m2w) / 2.0) : 1.0;
// Radial and tangential variances from the moments. The domain is the r1 disk, which is
// azimuthally symmetric and so adds no anisotropy of its own. The cells are SIGNED: away from the
// peak a learned cell is background noise centred on zero, and clamping it at zero turns that
// noise into a positive pedestal that the rad^2 / tan^2 weights read as extra width.
auto widths = [](const Moments &m) {
Sigma2 s;
if (m.w > 0.0) {
s.rad = std::max(0.25, m.rad / m.w);
s.tan = std::max(0.25, m.tan / m.w);
}
return s;
};
// Normalised profile (sum = 1): empirical average grid, or an isotropic Gaussian of the measured
// 2nd moment (only used by ProfileEmpirical; ProfileGaussian rebuilds per reflection in Pass B).
auto build_profile = [&](const std::vector<double> &grid, int n) {
// Normalised empirical profile (sum = 1), the average grid over the strong spots of a shell.
// ProfileGaussian does not use it - it rebuilds a per-reflection Gaussian in Pass B.
auto build_profile = [&](const std::vector<double> &grid) {
std::vector<double> P(GG, 0.0);
if (n <= 0) return P;
double sum = 0.0;
for (int k = 0; k < GG; ++k) {
const double g = std::max(0.0, grid[k]);
sum += g;
if (empirical) P[k] = g;
P[k] = std::max(0.0, grid[k]);
sum += P[k];
}
if (sum <= 0.0) return P;
if (empirical) {
if (sum > 0.0)
for (double &p : P) p /= sum;
} else {
const double sigma2 = measure_sigma2(grid);
double gsum = 0.0;
for (int dy = -R; dy <= R; ++dy)
for (int dx = -R; dx <= R; ++dx) {
const double g = std::exp(-(dx * dx + dy * dy) / (2.0 * sigma2));
P[grid_idx(dx, dy)] = g;
gsum += g;
}
for (double &p : P) p /= gsum;
}
return P;
};
const std::vector<double> global_P = build_profile(global_grid, global_n);
const double global_sigma2 = global_n > 0 ? measure_sigma2(global_grid) : 1.0;
std::vector<std::vector<double>> shell_P(N_SHELL);
std::vector<double> shell_sigma2(N_SHELL, global_sigma2);
const std::vector<double> global_P = empirical && global_n > 0 ? build_profile(global_grid)
: std::vector<double>(GG, 0.0);
const Sigma2 global_sigma2 = widths(global_mom);
std::vector<std::vector<double>> shell_P(N_SHELL, global_P);
std::vector<Sigma2> shell_sigma2(N_SHELL, global_sigma2);
for (int s = 0; s < N_SHELL; ++s) {
if (shell_n[s] >= MIN_STRONG_PER_SHELL) {
shell_P[s] = build_profile(shell_grid[s], shell_n[s]);
shell_sigma2[s] = measure_sigma2(shell_grid[s]);
} else {
shell_P[s] = global_P;
}
if (shell_n[s] < MIN_STRONG_PER_SHELL) continue;
if (empirical) shell_P[s] = build_profile(shell_grid[s]);
shell_sigma2[s] = widths(shell_mom[s]);
}
// --- Pass B: profile-fit each reflection (Kabsch, de-biased variance v = B + I*P; iterate). ---
@@ -352,12 +353,16 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y;
const double Rpx = std::hypot(rx, ry);
const double tan2t = Rpx / F_px;
const double s2t = shell_sigma2[sh];
const double s2t = shell_sigma2[sh].tan;
double s2r = s2t, ux = 1.0, uy = 0.0;
bool elong = false;
if (use_ellipse) {
// Radial excess over the tangential width: measured where the peak is resolved inside
// the r1 disk, with the analytic bandwidth + parallax/capture term as the floor. The
// analytic term is what carries a streak the disk is too small to measure.
const double sbw = bw_sigma * Rpx;
const double radial_extra = sbw * sbw + c_radial * tan2t * tan2t;
const double radial_extra = std::max(shell_sigma2[sh].rad - s2t,
sbw * sbw + c_radial * tan2t * tan2t);
if (Rpx > 1e-6 && radial_extra > 0.25) {
ux = rx / Rpx; uy = ry / Rpx;
s2r = s2t + radial_extra;