leonarski_fandClaude Opus 5
6f7b136ec2
Nothing kept a neighbour's flux out of a reflection's own signal disk. The union mask keeps neighbour cores out of the BACKGROUND ring, but the r1 disk was read whole, so on a dense pattern a crowded reflection measures part of its neighbour as its own. Ownership is decided once per image into a per-pixel (quantised distance, reflection) key written with an atomic minimum, so the nearest predicted centre wins whatever order the writes arrive in and the lowest index breaks a tie. `--overlap exclude`, now the default, drops the pixels a nearer neighbour owns from the profile fit. A profile fit is the amplitude of a normalised profile, so leaving pixels out renormalises the estimator by construction and the reflection stays unbiased rather than being discarded; the summation-fallback guard is scaled back to the disk the box-sum seed actually read, so it still compares like with like. `--overlap reject` is the XDS MINPK alternative - drop the reflection when less than `--overlap-minpk` of its expected profile is cleanly its own. A box sum has no profile to renormalise with, so `exclude` is a no-op there and only `reject` acts on it. Widening the split - keeping a pixel only where no other centre is within its distance PLUS a margin - was built and measured, and it is worse monotonically: the residual bias of the pixels that were kept grows from +0.072 to +0.209 in ln intensity at 0 to 3 px of margin. What the margin removes is the reflection's own profile, not the neighbour's tail, so the plain nearest-centre split is the rule. Measured on the full 38-crystal rotation battery against the same binary with the treatment off: ISa better 15 / worse 8, summed shortfall against XDS 39.7 -> 28.1. Three of the losses are the two-pass loop taking its other branch - their median mosaicity moves between the two known attractors - rather than the change under test; excluding those it is better 15 / worse 5 and the shortfall goes 31.3 -> 14.4. The two crowded crystals gain 38% and 52% of their ISa, one of them passing XDS. High-shell CC1/2 over the 35 crystals that neither flipped branch nor carry a collapsed error model is better 7 / worse 7. Space groups unchanged at 35/38. The owner map is built only when a treatment is asked for and costs 1.1% of the battery's wall clock - 23% on a genuinely crowded crystal, nothing where no two predictions touch. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
573 lines
31 KiB
C++
573 lines
31 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;
|
|
}
|
|
}
|
|
|
|
// --- 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). ---
|
|
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");
|
|
}
|
|
}
|