leonarski_fandClaude Opus 5
61d24db59f
The signal disk and the r2..r3 background ring were fixed pixel circles, identical for every reflection at every resolution. A reflection is not round: a finite bandwidth streaks it radially by bw_sigma*Rpx, so at high resolution the ring sits within 1.3-2.2 sigma of the reflection's own profile and measures its tails as background. --integration-stencil <k> makes the RING an ellipse, elongated along the beam->reflection direction by k times that streak, capped at 2*r3. The tangential half-widths stay r2 and r3, and the r1 signal disk stays a circle: r1 drives the all-or-nothing n_inner_valid == n_inner gate, so growing it rejects any reflection carrying one bad pixel along a long streak, and the flux a circular r1 loses is a function of resolution alone, which the per-shell scale absorbs. The geometry lives in one shared header compiled by both the host compiler and nvcc, so the seven pixel-classification sites - the CPU mask/main/clip loops and the GPU mark_mask/main/trim/clip kernels - cannot drift apart. Rather than evaluate an ellipse, each pixel's squared distance has its radial part scaled down, d2 - q*rad^2 against r2^2/r3^2 with q = 1 - (r/(r+grow))^2, so grow = 0 gives q = 0 and both tests collapse onto d2 exactly in floating point. The width is the bandwidth streak alone, not the profile's full radial variance, which also carries the sensor parallax and weak-spot capture terms. Deriving the growth from those was implemented first and measured on the rotation battery: at k=1 it took Thau_9's high-shell CC1/2 from 75.8 to 27.9 and Benas_3's from 14.1 to 6.0, against cytC_10 +1.2 and lyso_ref flat. On a monochromatic beam they are the only terms there are, and C_CAPTURE is 64% of them. Keeping only the streak also makes the option exactly inert without a bandwidth, rather than merely small. Default 0. Measured on broadband rotation data with the bandwidth set to its spectroscopic value, matched resolution limits: high-shell CC1/2 30.6 -> 46.4 at k=4, and better in EVERY shell in both CC1/2 and R_meas (top shell R_meas 194.7% -> 138.7%), with completeness, multiplicity and space group unchanged and 28 of 98833 unique reflections lost. Anomalous peak height over 18 sites +0.107 +- 0.039 sigma (p = 0.013). The full 38-crystal rotation battery is unchanged to every reported digit, base against k=3. Two consequences of an elongated ring are handled rather than inherited. The neighbour exclusion marks the inner ELLIPSE in each neighbour's own frame, or an elongated neighbour leaks its tails into this reflection's ring. And the radial-background curvature kernel becomes a small table indexed by the growth, because its azimuthal average makes one kernel serve every reflection only while their stencils are identical; the GPU's radial window, previously a fixed 32 bins, is now sized on the host from the widest ring on the detector. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
489 lines
26 KiB
C++
489 lines
26 KiB
C++
// 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;
|
|
}
|
|
}
|
|
|
|
// --- 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 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;
|
|
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) {
|
|
++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);
|
|
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;
|
|
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;
|
|
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). ---
|
|
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;
|
|
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;
|
|
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);
|
|
if (std::abs(I - 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");
|
|
}
|
|
}
|