mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 67 Wiki Activity
Files
43627e22dc88fd250ce200389a24245a035dfa5a
Jungfraujoch/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp
T
leonarski_fandClaude Opus 5 43627e22dc docs: a rule for licences and academic credit, and apply it 
Several methods adopted recently came from other crystallographic packages - the screw-absence test
from POINTLESS, MINPK and the profile-fit reweighting from XDS/Otwinowski, the CC1/2 cutoff and merge
outlier rejection from DIALS, the per-frame indexing gate from CrystFEL - and nothing in the
repository said where such a debt is recorded. The licence side was already worked out (licences beside
the vendored code, verbatim texts in licenses/ collected by COLLECT.sh, a row in THIRD_PARTY_NOTICES.md,
all installed under share/doc/jfjoch); the credit side was ad hoc.

Write the rule into CLAUDE.md. It states the distinction that matters: vendoring or linking someone's
CODE creates a LICENCE obligation, discharged in licenses/ and THIRD_PARTY_NOTICES.md; reimplementing
an algorithm from a PAPER creates none of that but creates an obligation of academic CREDIT, discharged
in docs/ACKNOWLEDGEMENT.md and in a comment at the algorithm. Neither substitutes for the other, and
taking both source and paper incurs both. It also fixes the citation form (authors, title, year,
journal, volume, pages, verified DOI), and says in-source credit goes at the algorithm, not the file
header, in the one-line style the code already uses.

Then bring the repository into compliance for the works concerned: docs/ACKNOWLEDGEMENT.md gains a
section acknowledging XDS, DIALS, POINTLESS/CCP4, MOSFLM, CrystFEL, GEMMI, the Kabsch/Otwinowski
profile fit, the Diederichs & Karplus statistics and the IUCr nomenclature reports, each with a DOI
checked against Crossref; docs/CPU_DATA_ANALYSIS.md's reference list gains the ones it was missing;
and four algorithms gain a line naming their source where no adjacent comment carried one.

No licence change. licenses/ and THIRD_PARTY_NOTICES.md are untouched.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-08-11 16:34:48 +02:00

575 lines
31 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "BraggIntegrationEngineCPU.h"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <limits>
#include <vector>
#include "../../common/CompressedImage.h"
#include "../../common/JFJochException.h"
using namespace bragg_engine;
namespace {
// The engine reads pixels in the INT32_MIN(masked)/INT32_MAX(saturated) convention.
inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MAX; }
// Identity sampler over the preprocessed int32 buffer (already in that convention).
struct BufferSampler {
const int32_t *p;
int32_t operator[](size_t i) const { return p[i]; }
};
// Sampler over a raw detector image of pixel type T: masked pixels carry the type minimum, saturated
// the type maximum (the FPGA image has no lossy-codec +/-1 band). Only the pixels actually read - the
// reflection disks - are converted, so there is no whole-image pass.
template <class T>
struct ImageSampler {
const T *p;
int64_t special_value;
int64_t saturation;
int32_t operator[](size_t i) const {
const int64_t v = p[i];
if (v == special_value) return INT32_MIN;
if (v == saturation) return INT32_MAX;
return static_cast<int32_t>(v);
}
};
} // namespace
BraggIntegrationEngineCPU::BraggIntegrationEngineCPU(const DiffractionExperiment &experiment)
: BraggIntegrationEngine(experiment) {}
template <class Sampler>
std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
const std::vector<Reflection> &predicted,
size_t npredicted, int64_t image_number) {
std::vector<BraggFitResult> results(npredicted);
if (npredicted == 0)
return Finalize(predicted, npredicted, results, image_number);
const int W = static_cast<int>(xpixel), H = static_cast<int>(ypixel);
const bool do_clip = bkg_clip_nsigma > 0.0f && mode != IntegratorMode::BoxSum;
// Symmetric trimmed-mean background fraction (BraggIntegrationSettings, rugnux --background-trim):
// replaces the r2..r3 ring MEAN with an f-trimmed mean, robust to the high-side contamination
// (neighbour wings, tails) that biases the plain mean up and makes it over-subtract weak high-angle
// reflections. The base ctor has already forced it to 0 whenever the clip below is in force, which
// is the default - the two are alternatives.
const double bkg_trim_frac = bkg_trim;
auto grid_idx = [this](int dx, int dy) { return (dy + R) * G + (dx + R); };
// --- Reflection mask: mark the r2 signal region of every predicted reflection so a neighbour's
// signal is excluded from this reflection's r2..r3 background ring. The region is the INNER
// stencil ellipse, taken in each neighbour's OWN frame - an elongated neighbour whose streak
// is still marked as a disk would leak its tails into this reflection's ring. ---
std::vector<uint8_t> refl_mask(npixel, 0);
for (size_t i = 0; i < npredicted; ++i) {
const auto &r = predicted[i];
const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil);
const int x0 = std::max(0, static_cast<int>(std::floor(r.predicted_x - st.ex_in - 1.0f)));
const int x1 = std::min(W - 1, static_cast<int>(std::ceil(r.predicted_x + st.ex_in + 1.0f)));
const int y0 = std::max(0, static_cast<int>(std::floor(r.predicted_y - st.ey_in - 1.0f)));
const int y1 = std::min(H - 1, static_cast<int>(std::ceil(r.predicted_y + st.ey_in + 1.0f)));
for (int y = y0; y <= y1; ++y)
for (int x = x0; x <= x1; ++x) {
const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y);
if (d.inner < r2_sq) refl_mask[y * W + x] = 1;
}
}
// --- Signal-region ownership: a pixel inside two reflections' signal regions belongs to the
// NEARER predicted centre. The union mask above cannot answer that - it also marks a
// reflection's own core - so ownership gets its own map, one (distance, reflection) key per
// pixel (BraggStencil.h). Built only when an overlap treatment is asked for. ---
// A box sum has no profile to renormalise a disk it has taken pixels out of, so it never excludes.
const bool exclude = overlap == OverlapMode::Exclude && mode != IntegratorMode::BoxSum;
std::vector<uint32_t> owner;
if (overlap != OverlapMode::Off) {
owner.assign(npixel, BRAGG_OWNER_NONE);
const float claim_sq = claim * claim;
for (size_t i = 0; i < npredicted; ++i) {
const auto &r = predicted[i];
const int x0 = std::max(0, static_cast<int>(std::floor(r.predicted_x - claim)));
const int x1 = std::min(W - 1, static_cast<int>(std::ceil(r.predicted_x + claim)));
const int y0 = std::max(0, static_cast<int>(std::floor(r.predicted_y - claim)));
const int y1 = std::min(H - 1, static_cast<int>(std::ceil(r.predicted_y + claim)));
for (int y = y0; y <= y1; ++y)
for (int x = x0; x <= x1; ++x) {
const float dx = x - r.predicted_x, dy = y - r.predicted_y;
const float d2 = dx * dx + dy * dy;
if (d2 >= claim_sq) continue;
uint32_t &o = owner[y * W + x];
o = std::min(o, BraggOwnerKey(std::sqrt(d2), inv_claim, static_cast<int>(i)));
}
}
}
auto clean = [&](int x, int y, size_t i) {
return owner.empty()
|| BraggOwnedBy(owner[static_cast<size_t>(y) * W + x], static_cast<int>(i));
};
// --- Pass A: box-sum every reflection (rough I, background, centroid, strong flag). ---
struct Rough {
double I = 0.0, sigma = NAN, bkg = 0.0, obs_x = 0.0, obs_y = 0.0;
double bkg_var = 0.0; // variance of the background ESTIMATE itself, bkg / n_bkg
double var_bkg = 0.0; // total non-signal variance carried to the merge
int64_t I_sum = 0; // kept so I can be rebuilt after the radial background correction
int n_inner = 0;
int n_disk = 0, n_own = 0; // signal-disk pixels, and how many of them are this reflection's
int r_bin = 0; // rounded distance from the beam centre, indexes the radial curve
int k_bin = 0; // which radial-background kernel this reflection's stencil needs
int cx = 0, cy = 0, shell = -1;
bool ok = false, strong = false, has_obs = false;
};
std::vector<Rough> rough(npredicted);
double inv_d2_min = std::numeric_limits<double>::max(), inv_d2_max = 0.0;
std::vector<int32_t> bkg_vals; // reused per reflection for the trimmed-mean background (idea 1)
// Radial background curve, accumulated from the annulus pixels this pass already reads. A pixel's
// radius is the reflection's radius plus the pixel's projection on the beam->reflection direction,
// so no per-pixel sqrt is needed.
const int n_rad = static_cast<int>(std::ceil(std::hypot(std::max<double>(beam_x, W - beam_x),
std::max<double>(beam_y, H - beam_y)))) + 2;
std::vector<double> rad_sum;
std::vector<int> rad_cnt;
if (bkg_radial) {
rad_sum.assign(n_rad, 0.0);
rad_cnt.assign(n_rad, 0);
}
for (size_t i = 0; i < npredicted; ++i) {
const auto &r = predicted[i];
Rough out;
// This reflection's stencil: the r1 signal disk, and the r2..r3 ring elongated radially by
// the analytic smear. The bounding box spans the OUTER ellipse, tightly - taking the largest
// semi-axis in both directions instead would read up to 60% more pixels for nothing.
const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil);
const int x0 = std::max(0, static_cast<int>(std::floor(r.predicted_x - st.ex_out - 1.0f)));
const int x1 = std::min(W - 1, static_cast<int>(std::ceil(r.predicted_x + st.ex_out + 1.0f)));
const int y0 = std::max(0, static_cast<int>(std::floor(r.predicted_y - st.ey_out - 1.0f)));
const int y1 = std::min(H - 1, static_cast<int>(std::ceil(r.predicted_y + st.ey_out + 1.0f)));
// Both from the stencil's own radius, so nothing downstream is derived from a second one.
out.r_bin = std::clamp(static_cast<int>(std::lround(st.r0)), 0, n_rad - 1);
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;
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;
bkg_vals.clear();
for (int y = y0; y <= y1; ++y)
for (int x = x0; x <= x1; ++x) {
const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y);
const int32_t px = img[y * W + x];
if (d.signal < r1_sq) {
// A pixel a nearer neighbour owns carries that neighbour's flux, so in Exclude
// mode it leaves the disk entirely - the sum, the pixel count the background is
// subtracted with, and the all-or-nothing validity gate alike. The box sum only
// counts how many were lost, which is all it can act on.
++n_disk;
if (overlap != OverlapMode::Off) {
if (clean(x, y, i)) ++n_own;
else if (exclude) continue;
}
++n_inner;
if (!valid(px)) continue;
I_sum += px;
I_sum_x += static_cast<int64_t>(x) * px;
I_sum_y += static_cast<int64_t>(y) * px;
++n_inner_valid;
} else if (d.inner >= r2_sq && d.outer < r3_sq) {
if (refl_mask[y * W + x]) continue;
if (!valid(px)) continue;
bkg_sum += static_cast<double>(px);
if (bkg_trim_frac > 0.0) bkg_vals.push_back(px);
++n_bkg;
}
}
int n_bkg_used = n_bkg; // pixels behind the FINAL background value (trim/clip shrink it)
if (n_inner_valid == n_inner && n_bkg > 5) {
out.bkg = bkg_sum / n_bkg;
if (bkg_trim_frac > 0.0 && bkg_vals.size() > 5
&& bkg_vals.size() <= static_cast<size_t>(bragg_engine::BKG_TRIM_MAX)) {
// Symmetric trimmed mean over the background ring (idea 1): drop the lowest and highest
// bkg_trim_frac of the pixels, average the rest. Robust to the high-side contamination
// that biases the plain ring mean and makes it over-subtract at high resolution.
std::sort(bkg_vals.begin(), bkg_vals.end());
const size_t lo = static_cast<size_t>(bkg_vals.size() * bkg_trim_frac);
const size_t hi = bkg_vals.size() - lo;
if (hi > lo) {
double s = 0.0;
for (size_t t = lo; t < hi; ++t) s += bkg_vals[t];
out.bkg = s / static_cast<double>(hi - lo);
n_bkg_used = static_cast<int>(hi - lo);
}
} else if (do_clip) {
// One high-outlier sigma-clip pass on the background ring: reject pixels above
// mean + n*sqrt(mean) to strip a neighbour core or a zinger that biases the mean.
const double thr = out.bkg + bkg_clip_nsigma * std::sqrt(std::max(out.bkg, 1.0));
double s = 0.0;
int n = 0;
for (int y = y0; y <= y1; ++y)
for (int x = x0; x <= x1; ++x) {
const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y);
if (!(d.inner >= r2_sq && d.outer < r3_sq)) continue;
if (refl_mask[y * W + x]) continue;
const int32_t px = img[y * W + x];
if (!valid(px)) continue;
if (static_cast<double>(px) <= thr) {
s += px;
++n;
if (bkg_radial) {
// The radial curve is binned on the TRUE detector radius, so the
// offset here is the unshrunk radial projection.
const int b = std::clamp(static_cast<int>(std::lround(st.r0 + d.rad)), 0, n_rad - 1);
rad_sum[b] += static_cast<double>(px);
++rad_cnt[b];
}
}
}
if (n > 5) { out.bkg = s / n; n_bkg_used = n; }
}
out.I = static_cast<double>(I_sum) - static_cast<double>(n_inner) * out.bkg;
// I = I_sum - n_inner*bkg, and bkg is itself estimated from n_bkg_used pixels, so its
// error enters n_inner times over: var(I) = I_sum + n_inner^2 * bkg/n_bkg_used. Leaving
// the second term out understates sigma by sqrt(1 + n_inner/n_bkg) - 1.109x at the
// default r1=4/r2=6/r3=10 stencil, on every reflection of every dataset.
out.bkg_var = out.bkg / n_bkg_used;
out.var_bkg = static_cast<double>(n_inner) * out.bkg
+ static_cast<double>(n_inner) * n_inner * out.bkg_var;
out.I_sum = I_sum;
out.n_inner = static_cast<int>(n_inner);
out.n_disk = n_disk;
out.n_own = n_own;
const double var_bkg_term = static_cast<double>(n_inner) * n_inner * out.bkg_var;
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);
out.has_obs = true;
}
out.cx = static_cast<int>(std::lround(r.predicted_x));
out.cy = static_cast<int>(std::lround(r.predicted_y));
out.ok = true;
out.strong = out.sigma > 0.0 && out.I / out.sigma >= STRONG_I_OVER_SIGMA;
if (r.d > 0.0f) {
const double inv_d2 = 1.0 / (static_cast<double>(r.d) * r.d);
inv_d2_min = std::min(inv_d2_min, inv_d2);
inv_d2_max = std::max(inv_d2_max, inv_d2);
}
}
rough[i] = out;
}
// --- Radial background curvature correction. The annulus mean is blind to the curvature of the
// radial background (a linear background cancels between the concentric disk and annulus), so
// correct it by mean_annulus(B) - mean_disk(B) taken from the curve just accumulated. Reads no
// pixels: one short dot product per reflection. ---
if (bkg_radial) {
for (size_t i = 0; i < npredicted; ++i) {
auto &rh = rough[i];
if (!rh.ok) continue;
// An empty bin contributes the reflection's own background, so a fully empty
// neighbourhood gives corr == 0 exactly (the kernel weights sum to zero).
const float *kern = k_diff.data() + static_cast<size_t>(rh.k_bin) * k_len;
double corr = 0.0;
for (int k = 0; k < k_len; ++k) {
const int b = std::clamp(rh.r_bin + k - k_off, 0, n_rad - 1);
const double v = rad_cnt[b] > 0 ? rad_sum[b] / rad_cnt[b] : rh.bkg;
corr += static_cast<double>(kern[k]) * v;
}
rh.bkg -= corr; // annulus mean -> mean over the signal disk
rh.I = static_cast<double>(rh.I_sum) - static_cast<double>(rh.n_inner) * rh.bkg;
}
}
// --- BoxSum mode is BraggIntegrate2D: emit the rough result directly. ---
if (mode == IntegratorMode::BoxSum) {
for (size_t i = 0; i < npredicted; ++i) {
const auto &rh = rough[i];
if (!rh.ok) continue;
// A box sum measures what is in the disk with no model of what should be there, so it can
// neither renormalise nor tell a neighbour's photon from its own: dropping the reflection
// is the only treatment it has, and it is applied only where that is what was ASKED for.
// Exclude means "leave the shared pixels out of the fit", and a box sum has no fit, so it
// is a no-op here rather than a rejection the caller never requested. The fraction is by
// AREA, not by profile mass, so the same threshold cuts harder here than in the profile
// modes.
if (overlap == OverlapMode::Reject && rh.n_own < overlap_min_peak * rh.n_disk) continue;
results[i] = {static_cast<float>(rh.I), static_cast<float>(rh.sigma), static_cast<float>(rh.bkg),
static_cast<float>(rh.obs_x), static_cast<float>(rh.obs_y),
static_cast<float>(rh.var_bkg), true, rh.has_obs};
}
return Finalize(predicted, npredicted, results, image_number);
}
auto shell_of = [&](float d) {
if (!(d > 0.0f) || inv_d2_max <= inv_d2_min) return 0;
const double t = (1.0 / (static_cast<double>(d) * d) - inv_d2_min) / (inv_d2_max - inv_d2_min);
return std::clamp(static_cast<int>(t * N_SHELL), 0, N_SHELL - 1);
};
for (size_t i = 0; i < npredicted; ++i)
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;
if (x < 0 || y < 0 || x >= W || y >= H) continue;
const int32_t px = img[y * W + x];
if (!valid(px)) continue;
if (exclude && !clean(x, y, i)) continue;
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;
}
// 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 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);
double sum = 0.0;
for (int k = 0; k < GG; ++k) {
P[k] = std::max(0.0, grid[k]);
sum += P[k];
}
if (sum > 0.0)
for (double &p : P) p /= sum;
return P;
};
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) 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). The
// reweighting is the Kabsch/Otwinowski iteration: Kabsch, Acta Cryst D66, 133-144 (2010);
// Otwinowski & Minor, Methods Enzymol 276, 307-326 (1997). ---
std::vector<double> Pbuf;
for (size_t i = 0; i < npredicted; ++i) {
const auto &rh = rough[i];
if (!rh.ok) continue;
const int sh = rh.shell < 0 ? 0 : rh.shell;
int Rf = R;
const std::vector<double> *Pvec = &shell_P[sh];
if (!empirical) {
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].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 = 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;
elong = true;
}
}
// Build the Gaussian per reflection, centred on the sub-pixel predicted position and (when
// needed) radially elongated, on a grid grown to hold the streak.
const double fx = predicted[i].predicted_x - rh.cx, fy = predicted[i].predicted_y - rh.cy;
Rf = elong ? std::min(3 * R, static_cast<int>(std::ceil(r2 + 2.0 * std::sqrt(s2r)))) : R;
const int Gf = 2 * Rf + 1;
Pbuf.assign(static_cast<size_t>(Gf) * Gf, 0.0);
double gs = 0.0;
for (int dy = -Rf; dy <= Rf; ++dy)
for (int dx = -Rf; dx <= Rf; ++dx) {
const double ex = dx - fx, ey = dy - fy;
const double rad = ex * ux + ey * uy, tn = -ex * uy + ey * ux;
const double g = std::exp(-rad * rad / (2.0 * s2r) - tn * tn / (2.0 * s2t));
Pbuf[(dy + Rf) * Gf + (dx + Rf)] = g;
gs += g;
}
for (double &p : Pbuf) p /= gs;
Pvec = &Pbuf;
}
const int Gf = 2 * Rf + 1;
// --- How much of the expected profile is cleanly this reflection's own. p_own is the whole
// grid's clean mass, i.e. XDS's MINPK quantity, and it is what Reject cuts on. m_own /
// m_all is the same fraction over the r1 disk alone, which is what the summation seed the
// runaway guard below compares against actually saw; with nothing excluded it is 1 and
// the guard is untouched. ---
double p_own = 1.0, m_all = 0.0, m_own = 0.0;
if (overlap != OverlapMode::Off) {
p_own = 0.0;
for (int dy = -Rf; dy <= Rf; ++dy)
for (int dx = -Rf; dx <= Rf; ++dx) {
const double Pp = (*Pvec)[(dy + Rf) * Gf + (dx + Rf)];
if (Pp <= 0.0) continue;
const int x = rh.cx + dx, y = rh.cy + dy;
if (x < 0 || y < 0 || x >= W || y >= H) continue;
if (!valid(img[y * W + x])) continue;
const bool own = clean(x, y, i);
if (own) p_own += Pp;
if (dx * dx + dy * dy < r1_sq) {
m_all += Pp;
if (own) m_own += Pp;
}
}
}
if (overlap == OverlapMode::Reject && p_own < overlap_min_peak) continue;
const double B = std::max(rh.bkg, PIXEL_VARIANCE_FLOOR);
double I = rh.I, den = 0.0, wsum = 0.0;
for (int iter = 0; iter < 4; ++iter) {
double num = 0.0;
den = 0.0;
wsum = 0.0;
for (int dy = -Rf; dy <= Rf; ++dy)
for (int dx = -Rf; dx <= Rf; ++dx) {
const double Pp = (*Pvec)[(dy + Rf) * Gf + (dx + Rf)];
if (Pp <= 0.0) continue;
const int x = rh.cx + dx, y = rh.cy + dy;
if (x < 0 || y < 0 || x >= W || y >= H) continue;
const int32_t px = img[y * W + x];
if (!valid(px)) continue;
if (exclude && !clean(x, y, i)) continue;
const double v = std::max(B + I * Pp, WEIGHT_VARIANCE_MIN_FRACTION * B);
num += Pp * (static_cast<double>(px) - rh.bkg) / v;
den += Pp * Pp / v;
wsum += Pp / v;
}
if (den > 0.0) I = num / den; else break;
}
if (!(den > 0.0)) continue;
// Guard against profile-fit runaways: on a weak / near-zero reflection the reweighted Kabsch
// iteration has no real peak to lock onto and can manufacture intensity the box sum never sees.
// Keep the profile intensity only if it agrees with the summation seed within the margin;
// otherwise fall back to the summation, which is robust there.
// 1/den is the profile-fit variance with the background taken as exact. The fit is
// I = sum(P*(px-bkg)/v) / sum(P^2/v), so dI/dbkg = -wsum/den and the background estimate's
// own error adds (wsum/den)^2 * var(bkg) - the same term the box sum was missing.
double sigma = std::sqrt(1.0 / den + (wsum / den) * (wsum / den) * rh.bkg_var);
double var_bkg = std::max(0.0, 1.0 / den - std::max(0.0, I)
+ (wsum / den) * (wsum / den) * rh.bkg_var);
// The seed is a sum over the disk the box sum actually read, so when Exclude has taken pixels
// out of both, the fit's full-profile intensity has to be scaled down to that same disk
// before the two are comparable. Nothing excluded gives exactly 1.
const double guard_scale = exclude && m_all > 0.0 ? m_own / m_all : 1.0;
if (std::abs(I * guard_scale - rh.I) > PROFILE_SUMMATION_MAX_NSIGMA * rh.sigma) {
I = rh.I;
sigma = rh.sigma;
var_bkg = rh.var_bkg;
}
// Carry the Pass-A box-sum intensity-weighted centroid (observed spot position) through the
// profile path too - post-refinement uses it as the observed position (beam-centre / distance).
results[i] = {static_cast<float>(I), static_cast<float>(sigma),
static_cast<float>(rh.bkg),
static_cast<float>(rh.obs_x), static_cast<float>(rh.obs_y),
static_cast<float>(var_bkg), true, rh.has_obs};
}
return Finalize(predicted, npredicted, results, image_number);
}
std::vector<Reflection> BraggIntegrationEngineCPU::Run(const ImagePreprocessorBuffer &image,
const std::vector<Reflection> &predicted,
size_t npredicted, int64_t image_number) {
if (image.size() != npixel)
return Finalize(predicted, npredicted, std::vector<BraggFitResult>(npredicted), image_number);
return RunImpl(BufferSampler{image.data()}, predicted, npredicted, image_number);
}
std::vector<Reflection> BraggIntegrationEngineCPU::Run(const CompressedImage &image,
const std::vector<Reflection> &predicted,
size_t npredicted, int64_t image_number) {
if (image.GetWidth() * image.GetHeight() != npixel)
return Finalize(predicted, npredicted, std::vector<BraggFitResult>(npredicted), image_number);
std::vector<uint8_t> scratch;
const auto *ptr = image.GetUncompressedPtr(scratch);
switch (image.GetMode()) {
case CompressedImageMode::Int8:
return RunImpl(ImageSampler<int8_t>{reinterpret_cast<const int8_t *>(ptr), INT8_MIN, INT8_MAX},
predicted, npredicted, image_number);
case CompressedImageMode::Int16:
return RunImpl(ImageSampler<int16_t>{reinterpret_cast<const int16_t *>(ptr), INT16_MIN, INT16_MAX},
predicted, npredicted, image_number);
case CompressedImageMode::Int32:
return RunImpl(ImageSampler<int32_t>{reinterpret_cast<const int32_t *>(ptr), INT32_MIN, INT32_MAX},
predicted, npredicted, image_number);
case CompressedImageMode::Uint8:
return RunImpl(ImageSampler<uint8_t>{reinterpret_cast<const uint8_t *>(ptr), UINT8_MAX, UINT8_MAX},
predicted, npredicted, image_number);
case CompressedImageMode::Uint16:
return RunImpl(ImageSampler<uint16_t>{reinterpret_cast<const uint16_t *>(ptr), UINT16_MAX, UINT16_MAX},
predicted, npredicted, image_number);
case CompressedImageMode::Uint32:
return RunImpl(ImageSampler<uint32_t>{reinterpret_cast<const uint32_t *>(ptr), UINT32_MAX, UINT32_MAX},
predicted, npredicted, image_number);
default:
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Image mode not supported");
}
}
Reference in New Issue View Git Blame Copy Permalink