leonarski_fandClaude Opus 5
4e6eb18d93
48008e144 widens the signal radius on crystals with wide spots. On one
battery crystal - simultaneously the widest-spot and among the highest in
mosaicity - the wider radius left its neighbours' background rings with too
few clean pixels and cost 28.5% of its observations. Pass 1 now measures how
often that happens and, above a bound, hands pass 2 the settings from before
the pre-scan widened them.
The obvious quantity does not work. On the total rate of reflections dropped
for a starved ring, the losing crystal reads 4.08% and the rule's four
biggest winners read 1.89-2.04% - and a crystal at the shipped radius reads
2.37%, above all of them. Re-running the winners at r1 = 4 shows why: they
read 2.20-2.32% there too, and widening moves them down. That floor is module
gaps, the beam stop and the resolution mask, which are properties of the
detector and do not move with the radius.
So the counter separates the two. A ring is neighbour-starved when it would
have kept more than five pixels but for the pixels a neighbouring
reflection's signal region occupies. That is exact rather than estimated: the
reflection mask marks the disk inside r2 and the ring is everything outside
it, so a masked ring pixel always belongs to some other reflection's core.
The separation goes from a factor of 2 to a factor of 13 - over the twelve
crystals the radius moves, the rate is 0.000 five times, 0.001 three times,
then 0.004, 0.235, 0.315 and 4.082 - and the bound is the log-space midpoint
of that one gap, 0.0113, a factor 3.6 clear of the nearest measurement on
either side.
Predicted reflection spacing does not separate them at all: the losing
crystal is 19th of 38, a winner sits at 21.9 px, and the loosest pattern in
the battery starves 1.93% of its rings.
Battery: the space group is identical on all 38 and the merged .hkl is
byte-identical on 37, so it is inert wherever it does not fire. On the one
crystal it fires on, <I/sigma> is up 18.3%, R_meas down 29.4%, observations
up 6.3%, CC1/2 0.944 to 0.974, and its two empty top shells come back as
numbers. Its indexing rate, refined distance, beam centre and cell are
bit-identical between the two arms, so this is the guard and not the two-pass
gate.
The counters are a shared channel through both engines, summed across
workers and logged once per pass; on the GPU it is one atomic add per dropped
reflection. The profile-fit runaway guard reports on the same channel, which
is the first measurement of its trip rate.
This does not recover that crystal fully. With the adaptive radius on, pass 1
reaches a different lattice and pass 2 indexes 21% fewer frames - which
happens before the measurement this guard reads exists, and is unaffected by
it. At matched indexing rate the guard recovers 96% of the baseline's
observations against 90.5% without it.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P
627 lines
35 KiB
C++
627 lines
35 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);
|
|
counts.predicted += npredicted;
|
|
|
|
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;
|
|
bool full = false; // every pixel of the signal disk was readable
|
|
};
|
|
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;
|
|
// Ring pixels a NEIGHBOUR's signal region occupies. refl_mask marks d.inner < r2_sq and the
|
|
// ring is d.inner >= r2_sq, so the two are complementary and a masked ring pixel is always
|
|
// someone else's - never this reflection's own core. That makes this an exact count of what
|
|
// the pattern's density took, separable from what the detector took.
|
|
int n_bkg_neighbour = 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]) { ++n_bkg_neighbour; 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)
|
|
// A masked, untrusted, gapped or overloaded pixel inside the signal disk used to discard the
|
|
// reflection outright. A profile fit does not need it to: the fit is the amplitude of a
|
|
// NORMALISED profile, so leaving pixels out renormalises the estimator by construction and
|
|
// widens sigma by the information they carried - the same argument Exclude makes for a
|
|
// neighbour's pixels. Pass B keeps the reflection only while enough of the expected profile
|
|
// survived to constrain the amplitude (XDS's MINPK, dials' valid_foreground_threshold). A box
|
|
// sum has no profile to renormalise with, so there it stays all or nothing.
|
|
const bool full = n_inner_valid == n_inner;
|
|
// A ring left with five or fewer clean pixels cannot estimate a background, so the reflection
|
|
// is dropped whole - the one thing the stencil geometry does to the DATA rather than to a
|
|
// measurement. Counted here because it is the only direct evidence of a radius that has
|
|
// outgrown the pattern it is integrating (BraggIntegrationCounts).
|
|
const bool keep_partial = full || mode != IntegratorMode::BoxSum;
|
|
if (keep_partial && n_bkg <= 5) {
|
|
++counts.bkg_starved;
|
|
// Would the ring have been enough without the neighbours? Then it is the pattern, not the
|
|
// detector, that took it.
|
|
if (n_bkg + n_bkg_neighbour > 5) ++counts.bkg_starved_by_neighbour;
|
|
}
|
|
if (keep_partial && 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; }
|
|
}
|
|
// The sum is over the pixels actually READ, so that is the count the background is
|
|
// subtracted with; with nothing missing it is the whole disk, exactly as before.
|
|
out.I = static_cast<double>(I_sum) - static_cast<double>(n_inner_valid) * 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_valid) * out.bkg
|
|
+ static_cast<double>(n_inner_valid) * n_inner_valid * out.bkg_var;
|
|
out.I_sum = I_sum;
|
|
out.n_inner = static_cast<int>(n_inner_valid);
|
|
out.n_disk = n_disk;
|
|
out.n_own = n_own;
|
|
const double var_bkg_term = static_cast<double>(n_inner_valid) * n_inner_valid * 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);
|
|
// A disk with a hole in it gives a centroid pulled away from the hole, and the hole
|
|
// sits at a fixed place on the detector - post-refinement would read that as geometry.
|
|
out.has_obs = full;
|
|
}
|
|
out.cx = static_cast<int>(std::lround(r.predicted_x));
|
|
out.cy = static_cast<int>(std::lround(r.predicted_y));
|
|
out.ok = true;
|
|
out.full = full;
|
|
// The profile, its resolution shells and their widths are learned from COMPLETE
|
|
// reflections: a partial one has no intensity to normalise its grid by and a hole to
|
|
// stack into it. So the learning sees exactly what it saw before this rescue existed.
|
|
out.strong = full && out.sigma > 0.0 && out.I / out.sigma >= STRONG_I_OVER_SIGMA;
|
|
if (full && 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 the fit can actually see. p_valid is the readable
|
|
// fraction of it - XDS's MINPK, dials' valid_foreground_threshold - measured against the
|
|
// mass that falls on the detector at all, so a reflection is judged on the pixels that
|
|
// exist and not on where the sensor ends. p_own is the same quantity over neighbour-owned
|
|
// pixels, what Reject cuts on. m_read / m_all is the fraction of the r1 disk the
|
|
// summation seed the runaway guard compares against actually summed; with nothing missing
|
|
// and nothing excluded it is 1 and the guard is untouched. ---
|
|
double p_grid = 0.0, p_valid = 0.0, p_own = 0.0, m_all = 0.0, m_read = 0.0;
|
|
double p_peak = 0.0, p_lost_peak = 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 bool in_disk = dx * dx + dy * dy < r1_sq;
|
|
p_grid += Pp;
|
|
p_peak = std::max(p_peak, Pp);
|
|
if (in_disk) m_all += Pp;
|
|
if (!valid(img[y * W + x])) {
|
|
p_lost_peak = std::max(p_lost_peak, Pp);
|
|
continue;
|
|
}
|
|
p_valid += Pp;
|
|
const bool own = clean(x, y, i);
|
|
if (own) p_own += Pp;
|
|
if (in_disk && (own || !exclude)) m_read += Pp;
|
|
}
|
|
if (p_valid < overlap_min_peak * p_grid) continue;
|
|
// A hole in the profile's PEAK is a different defect from a hole in its wings, and the mass
|
|
// fraction above cannot tell them apart - the peak of a broad spot is a few percent of the
|
|
// mass, so MINPK passes a reflection that has lost the one part of the profile its amplitude
|
|
// is determined by. See the header.
|
|
if (p_lost_peak > MINPK_MAX_MISSING_PEAK * p_peak) continue;
|
|
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);
|
|
// var_bkg is the NON-SIGNAL part of that variance, and 1/den is the fit's variance taken at the
|
|
// fitted intensity, so the signal part to remove is I itself - not max(0, I). Clamping it leaves
|
|
// a down-fluctuated reflection's own (deflated) variance standing as its background variance,
|
|
// which is 2-3x too small; the merge then weights exactly the down-fluctuated observations up.
|
|
// The whole point of a separate var_bkg is that it does not move with the observation's own
|
|
// fluctuation, and 1/den - I is what does not (1/den grows with I one for one).
|
|
double var_bkg = std::max(0.0, 1.0 / den - I
|
|
+ (wsum / den) * (wsum / den) * rh.bkg_var);
|
|
// The seed is a sum over the disk the box sum actually read, so when pixels are missing from
|
|
// both - excluded to a neighbour, or unreadable - the fit's full-profile intensity has to be
|
|
// scaled down to that same disk before the two are comparable. Nothing dropped gives 1.
|
|
const double guard_scale = m_all > 0.0 ? m_read / 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;
|
|
++counts.profile_fallback;
|
|
}
|
|
// 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");
|
|
}
|
|
}
|