leonarski_fandClaude Opus 5
e4ce573075
MINPK asks how MUCH of the expected profile is readable. It does not ask WHERE, and the two are not the same question. The renormalisation argument the rescue rests on - a fit over a subset of a normalised profile is unbiased - needs the pixels to go missing for reasons unrelated to the reflection. A gap, a mask or the edge of the sensor is such a reason: the loss is set by the detector, and the fit renormalises over what is left. A pixel invalidated BY THE FLUX IT SAW is not: it goes missing because the reflection was bright, and it is the peak. Measured on the combined fulls, against the mean of the complete observations of the same reflection, in the innermost resolution shell of the high-multiplicity control and of a weaker crystal: a rescued reflection whose unreadable pixel sits within a pixel of the predicted centre reads |I - <I>|/I of 0.50 and 0.53, against 0.073 and 0.212 for a complete observation - 6.8x and 2.5x - and carries several times the mean intensity of its shell. On the control that is 0.21% of the shell's observations supplying 1.77% of the R_meas numerator; on the weaker crystal 0.52% supplying 6.82%. Rescues that lost only rim pixels are unremarkable by the same measure, 1.19x and 0.88x. Dropping the peak-losers alone takes the shell's R_meas from 7.440% back to 7.315% (unrescued: 7.307%) and from 22.03% to 21.06% (unrescued: 21.28%) - which is the whole of the low-resolution R_meas the rescue cost, and on the second crystal rather more. Raw frames say what they are. The pattern is a dead-centre invalid pixel with 5878, 9875 and 27583 counts around it: the detector's per-frame invalid marker on the brightest reflections. MINPK cannot catch them because it cuts on profile MASS, and the peak of a broad spot is a few percent of the mass. So a second condition, in the loop that already measures the readable fraction: no unreadable pixel may carry more than 0.9 of the profile's own peak value. A fraction of the peak rather than a radius in pixels because the peak is as wide as the spot - for a Gaussian the cut is at sqrt(-2 ln f) sigma, 0.46 sigma here, which is the peak pixel alone where sigma is 0.8 px and the crest of the ridge where the profile is a bandwidth streak. Swept against the alternatives on two crystals: a fixed radius needs 1.0-1.5 px to do the same work and costs 3-9x more observations for it, and 0.5 px does not reach the peak of a sub-pixel-offset prediction at all; tightening the fraction to 0.5 or 0.2 buys nothing beyond 0.9 and costs 7x more. Six crystals, three detectors, against the rescue as it stands: the rule keeps 99.86-99.96% of the recovered observations and returns R_meas to its unrescued value or below (4.6 -> 4.5%, 6.7 -> 6.6%, 25.1 -> 25.0%), R_meas in the innermost shell likewise (2.7 -> 2.6%, 5.9 -> 5.3% against 5.4% unrescued, 16.5 -> 16.4%), <I/sigma> up or level everywhere, and every unique reflection the rescue won is kept. Raising --overlap-minpk to 0.90 instead reaches the same place on two of them and short of it on the third, while discarding 0.8% of the recovered observations rather than 0.05%. An elongated pink-beam profile on a 9M detector and an EIGER2 16M dataset are both untouched at 99.9%, so the crest protection does not over-reject a streak. One crystal is not improved: a dataset whose error model rugnux declines to fit for want of strong reflections, whose <I/sigma> is <= 0 in eight of its ten shells and whose R_meas is undefined in as many. There the rule costs about 3% of <I/sigma> in the one shell that has signal, reproducibly, on top of the 9% the rescue itself costs there - while its overall R_meas moves 1.5 points on nothing but the thread count. The parity test gains four sections. Unreadable pixels were only ever punched into empty sky, so neither the rescue nor this rule had any CPU/GPU coverage at all; they now go into the signal disks - the peak of every fifth reflection, ~1.1 sigma out of every seventh, the disk edge of every eleventh - for both profile modes, a box sum and an elongated stencil, with a check that the clipping actually costs reflections so the coverage cannot go quietly vacuous. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Full 38-crystal rotation battery against the rescue without this rule, both on the same base: ISa better 19 / worse 4, +0.73 CC1/2 better 3 / worse 1, +1.3 R_meas_lo better 4 / worse 3, -0.3 space groups unchanged for 17 770 observations, 0.09 % of the run total and under 2 % of what the rescue had won. The two crystals whose peak-loss population was measured beforehand land on their predicted values: a tetragonal reference goes R_meas_lo 2.7 -> 2.6 % and ISa 27.11 -> 27.42, a cubic insulin 5.9 -> 5.3 % and 20.34 -> 20.65. One crystal pays: a cubic case with 2381 unique reflections goes R_meas 8.8 -> 9.6 % and ISa 4.08 -> 3.49. It is the crystal in the battery with the fewest uniques, so its rescued population is small and its shell statistics are coarse, but the loss is real and not noise in the R_meas. The R_meas sum over the battery reads +1.3, of which +3.2 is one crystal whose R_meas moves 1.5 points on thread count alone; without it the sum is negative. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
609 lines
34 KiB
C++
609 lines
34 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;
|
|
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;
|
|
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)
|
|
// 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;
|
|
if ((full || mode != IntegratorMode::BoxSum) && 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;
|
|
}
|
|
// 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");
|
|
}
|
|
}
|