Radial background: with the frame's curve, the correction is taken over the pixels read, at their 2 theta
With the frame's curve in use (D5: stills, and rotation runs whose empty apertures ask for it) the correction was a precomputed annulus-minus-disk kernel over radial offset, applied to the curve at the reflection's distance from the PONI. The empty apertures that decide it compute the same quantity over the pixels they actually read, at each pixel's own 2 theta. The two differ exactly next to the beam stop: the ring's beam-side pixels are masked by the shadow, but the kernel still averaged the curve over them (the curve there holds the shadow's edge, a cliff of several hundred counts), and on a tilted detector the PONI is a few pixels from where the beam lands, where the curve is steepest. On a binned-Rayonix still run with the beam stop off the beam this sent low-resolution reflections next to the stop to large negative intensities (median -854 -> -4171 at 42-50 px), destabilised the per-image scales and collapsed the lowest shell's CC1/2 (0.906 -> 0.448, "not usable" warning). Now, where the frame's curve is in use: - the correction is the curve's mean over the ring pixels the background was taken from (after the clip) minus its mean over the readable disk pixels (BraggCurveCorrection), summed in pass A while the box is read anyway - in fixed point (1/256 count), so CPU and GPU sums are exact and order-independent; - the curve is read at the 2 theta of the reflection's position (CurveRadius, from the geometry), the radius it is tabulated against, and the clip offset follows it there; - BackgroundAperture takes the curve's mean over the ring pixels the clip kept, so the decision measures exactly what is applied. The kernel stays for the curve accumulated from the reflections' rings (rotation's per-image ice gate), unchanged. Metal unchanged (it hands such frames to the CPU engine). Validation. Rotation, CPU A/B vs 87684012f over 18 sets: .hkl byte-identical on the 12 sets where D5 is off (and on the 3 ice-gated D14 sets identical to D5); the three D5-on sets improve: a home-source CCD set R-free 0.210 -> 0.194 (D5 0.201), ISa 10.4 -> 10.8, anomalous 3.70 -> 3.94; the ice set 0.315 -> 0.312, ISa 22.3 -> 22.7 (D5 lost ISa there); the third neutral. Stills (0-cycle REFMAC, 9 sets): R-free within +0.0021/-0.0039 of 87684012f on every set, verdicts unchanged except one draw -> win, no "not usable" warning on any set (two base runs had one). The lowest shell's CC1/2 on the run above stays low (0.44): it is set by three to five reflections at 30-41 A in the penumbra of the stop, which integrate as garbage with every variant (the uncorrected run included); with the merge's low-resolution limit at 30 A it reads 0.79. Tests: BraggIntegrationEngineCPU_RadialCurveOverThePixelsRead (shadow edge on a tilted detector; fails with the kernel or with the PONI radius), the GPU parity test now with a masked ring. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu (cherry picked from commit d225b546c613ab2266c63ff2add02624b04d5e25)
This commit is contained in:
leonarski_f committed
1 parent
f199f23c54
commit
fac2a17428
11 files changed
+225
-51
No files matched your search
@@ -151,9 +151,9 @@ $
|
||||
$
|
||||
with $\kappa$ the annulus-minus-disk histogram of the stencil over radial offset, averaged over azimuth, and $\bar B(r)$ the image's own radial background curve. With the fixed circular stencil ($k=0$, §9.1) that single kernel serves every reflection. An elongated ring does not: its radial-offset histogram depends on how far that particular reflection's ring was grown, so $\kappa$ becomes a small table of kernels, indexed by the growth rounded to whole pixels. The azimuthal average survives the change unaltered, because the stencil is rebuilt in the reflection's own frame at each azimuth and so stays radially aligned: what is averaged over is the sub-pixel phase of the detector grid against the radius, which is what genuinely differs between reflections. Applying it costs one short dot product per reflection and no extra pixel reads; correcting the background *scalar* means the box sum, the profile fit and the non-signal variance all pick it up. The variance carried to the merge, $\mathrm{var}_\mathrm{bkg} = g\,(n_S\hat b + n_S^2\,\mathrm{var}(\hat b))$, moves with the corrected $\hat b$ in its first term — the Poisson share of the background *under the disk* — while $\mathrm{var}(\hat b)$, the noise of the ring estimate, and $\sigma$, taken from the counts in the disk, do not depend on a smooth curve and stay.
|
||||
|
||||
The curve $\bar B(r)$ is the adaptive spot finder's peak-excluded background per resolution ring of the **whole frame** (§3.2), one value per pixel of radius: every pixel of the frame, against the few dozen to few thousand reflections whose own rings would otherwise sample it. Where that curve is in use the ring's **high-side clip follows it**: a ring pixel at radial offset $u$ is clipped against the threshold plus the curve's rise $\bar B(r_0+u)-\bar B(r_0)$, in whole counts. On a background falling steeply away from the beam a flat threshold cuts the beam-side half of the ring as outliers and the background reads low — up to 150× too strong a low-resolution reflection on a binned CCD near the beam stop. Without a frame curve (fixed-threshold spot finding) the curve is accumulated from the clipped annulus pixels of the predicted reflections themselves, and the clip stays flat.
|
||||
The curve $\bar B(r)$ is the adaptive spot finder's peak-excluded background per resolution ring of the **whole frame** (§3.2), one value per pixel of radius: every pixel of the frame, against the few dozen to few thousand reflections whose own rings would otherwise sample it. Where that curve is in use the ring's **high-side clip follows it**: a ring pixel at radial offset $u$ is clipped against the threshold plus the curve's rise $\bar B(r_0+u)-\bar B(r_0)$, in whole counts. On a background falling steeply away from the beam a flat threshold cuts the beam-side half of the ring as outliers and the background reads low — up to 150× too strong a low-resolution reflection on a binned CCD near the beam stop. With the frame's curve the correction is not the kernel above but the curve itself averaged over the pixels the reflection actually reads: its mean over the ring pixels the clip kept minus its mean over the readable disk pixels, summed while the box is read anyway. The kernel is an average over a whole annulus, so it counts the curve wherever the ring has lost pixels — to the clip, to a neighbour's signal, to a masked or shadowed pixel — and next to the beam-stop shadow, where the curve falls steeply into the shadow and the ring's beam-side part is masked, that turned a correction of a few counts into tens and made low-resolution reflections near the stop read strongly negative. The curve is read at the $2\theta$ of each reflection's position — the radius $F\tan 2\theta$ it is tabulated against, read off the geometry as the curve's rings are — and not at its distance from the PONI, which on a tilted detector is a few pixels off exactly where the curve is steepest. Without a frame curve (fixed-threshold spot finding) the curve is accumulated from the clipped annulus pixels of the predicted reflections themselves, the kernel applies it, and the clip stays flat.
|
||||
|
||||
**When it is applied.** `--background-radial` takes `on`, `off` or `auto`; in Rugnux it is **`auto` by default** (the broker keeps it off), and auto is decided **per run, in the pre-scan, on empty apertures** — the same measurement that, on stills, also chooses how far the ring reaches (the $r_3$, up to the one the run starts from, whose ring predicts the empty disk best). On a grid of apertures away from every found spot, the ring's prediction of the empty disk is taken twice: plain, and with the correction applied exactly as the integrator applies it (the clip following the curve, the kernel's correction). With nothing in the aperture each difference is the error the ring makes in a reflection's background, noise plus bias. Rotation measures every stencil the spot width can give it (and the one it starts from, which the geometry pre-pass and the starvation guard return to) and reads the one it gets.
|
||||
**When it is applied.** `--background-radial` takes `on`, `off` or `auto`; in Rugnux it is **`auto` by default** (the broker keeps it off), and auto is decided **per run, in the pre-scan, on empty apertures** — the same measurement that, on stills, also chooses how far the ring reaches (the $r_3$, up to the one the run starts from, whose ring predicts the empty disk best). On a grid of apertures away from every found spot, the ring's prediction of the empty disk is taken twice: plain, and with the correction applied exactly as the integrator applies it (the clip following the curve, and the curve's mean over the ring pixels kept minus its mean over the disk). With nothing in the aperture each difference is the error the ring makes in a reflection's background, noise plus bias. Rotation measures every stencil the spot width can give it (and the one it starts from, which the geometry pre-pass and the starvation guard return to) and reads the one it gets.
|
||||
|
||||
- **Stills** have no other correction, so it is applied to every image of a run where it lowers the apertures' mean square error at all.
|
||||
- **Rotation** has one: under auto, each image whose peak-excluded ice score (§3.3) reaches `--ice-min-score` is corrected with the curve of its reflections' rings. The frame's curve on every image replaces that only where it takes **at least 10 %** of the apertures' mean square error away. On photon-counting sweeps the apertures gain 0.1–5 % (statistically significant, $z\approx3$–4), and correcting every image there was measured to cost ISa and anomalous signal; where it pays — the rings of home-source CCD frames, a smooth ice ring — they gain 13–31 %. A run below the margin keeps the per-image ice correction, bit for bit.
|
||||
|
||||
@@ -101,17 +101,17 @@ void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width
|
||||
// Every candidate ring as the integrator reads it: the mean, then one high-side clip pass. With the
|
||||
// radial correction, as the integrator applies it: the clip follows the curve (each pixel's
|
||||
// threshold raised by the curve's rise from the centre to it, in whole counts -
|
||||
// RadialClipOffset), and the correction is taken over the ring's geometry, as the kernel is.
|
||||
// RadialClipOffset), and the correction is the curve's mean over the ring pixels the clip kept
|
||||
// minus its mean over the disk (BraggCurveCorrection).
|
||||
std::vector<double> d_plain(r3.size()), d_radial(r3.size());
|
||||
bool measured = true;
|
||||
for (size_t k = 0; k < r3.size() && measured; ++k) {
|
||||
const float r3_sq = r3[k] * r3[k];
|
||||
double sum = 0.0, sum_b = 0.0;
|
||||
double sum = 0.0;
|
||||
int n = 0;
|
||||
for (const auto &p : ring)
|
||||
if (p.d2 < r3_sq) {
|
||||
sum += p.v;
|
||||
sum_b += p.b;
|
||||
++n;
|
||||
}
|
||||
if (n == 0) {
|
||||
@@ -120,7 +120,7 @@ void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width
|
||||
}
|
||||
const double mean = sum / n;
|
||||
const double thr = mean + clip_nsigma * std::sqrt(gain * std::max(mean, static_cast<double>(gain)));
|
||||
double kept = 0.0, kept_radial = 0.0;
|
||||
double kept = 0.0, kept_radial = 0.0, kept_radial_b = 0.0;
|
||||
int n_kept = 0, n_kept_radial = 0;
|
||||
for (const auto &p : ring) {
|
||||
if (p.d2 >= r3_sq)
|
||||
@@ -131,13 +131,14 @@ void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width
|
||||
}
|
||||
if (clip_nsigma <= 0.0f || p.v <= thr + static_cast<double>(std::lround(p.b - centre_b))) {
|
||||
kept_radial += p.v;
|
||||
kept_radial_b += p.b;
|
||||
++n_kept_radial;
|
||||
}
|
||||
}
|
||||
measured = n_kept > 0 && n_kept_radial > 0;
|
||||
if (measured) {
|
||||
d_plain[k] = disk_v - kept / n_kept;
|
||||
d_radial[k] = disk_v - kept_radial / n_kept_radial - (disk_b - sum_b / n);
|
||||
d_radial[k] = disk_v - kept_radial / n_kept_radial - (disk_b - kept_radial_b / n_kept_radial);
|
||||
}
|
||||
}
|
||||
if (!measured)
|
||||
|
||||
@@ -34,7 +34,8 @@
|
||||
// The same apertures also say whether the radial background-curvature correction helps: each difference
|
||||
// is taken a second time with the correction applied as the integrator applies it - the ring's clip
|
||||
// following the frame's own ring background (RingBackground), and that background averaged over the
|
||||
// ring's pixels minus over the disk's, which is what the integrator's kernel takes off a reflection.
|
||||
// ring pixels the clip kept minus over the disk's, which is what the integrator takes off a reflection
|
||||
// (BraggCurveCorrection).
|
||||
// Where the background really is a function of the radius the corrected ring predicts the disk better;
|
||||
// where it is not - a strongly polarised beam, panels of different response - it predicts it worse, and
|
||||
// the run measures which. On stills the ring is chosen on the plain differences first, so its reach does
|
||||
|
||||
@@ -302,18 +302,34 @@ public:
|
||||
// a few dozen reflections, against every pixel of the frame here. An empty profile goes back to
|
||||
// accumulating the curve from the reflections.
|
||||
void RadialBackgroundFromRings(const std::vector<float> &ring_bkg, int32_t q_bins, float low_q, float dq);
|
||||
// Where a reflection sits on that curve. The curve is tabulated against F_px tan(2 theta), so this is
|
||||
// the 2 theta of the reflection's position, read off the geometry as the curve's rings and the empty
|
||||
// apertures read every pixel's; a pixel of its stencil is at that radius plus its radial offset. Not
|
||||
// the distance from the PONI (BraggStencil::r0): on a tilted detector the beam lands elsewhere, and
|
||||
// near the beam stop, where the curve is steepest, those few pixels put a reflection on the wrong
|
||||
// part of it.
|
||||
[[nodiscard]] float CurveRadius(const BraggPanel &pn) const {
|
||||
const double s = geom.GetWavelength_A() / (2.0 * geom.PxlToRes(pn.plane_x, pn.plane_y));
|
||||
return static_cast<float>(F_px * std::tan(2.0 * std::asin(s)));
|
||||
}
|
||||
// Where that curve is in use (the correction on), the ring's high-side clip follows it as well: a ring
|
||||
// pixel at radial offset `rad` from a reflection at radius r0 is clipped against the threshold plus
|
||||
// the curve's rise from r0 to it, rounded to whole counts so that every engine keeps the same pixels.
|
||||
// On a background falling steeply away from the beam the flat threshold otherwise cuts the beam-side
|
||||
// half of the ring as outliers, and the background reads low. 0 where no curve is in use.
|
||||
[[nodiscard]] int32_t RadialClipOffset(float r0, float rad) const {
|
||||
// pixel at radial offset `rad` from a reflection at curve radius rc (CurveRadius) is clipped against
|
||||
// the threshold plus the curve's rise from rc to it, rounded to whole counts so that every engine
|
||||
// keeps the same pixels. On a background falling steeply away from the beam the flat threshold
|
||||
// otherwise cuts the beam-side half of the ring as outliers, and the background reads low. 0 where no
|
||||
// curve is in use.
|
||||
[[nodiscard]] int32_t RadialClipOffset(float rc, float rad) const {
|
||||
if (!bkg_radial || radial_curve.empty())
|
||||
return 0;
|
||||
return static_cast<int32_t>(std::lround(CurveAt(rc + rad) - CurveAt(rc)));
|
||||
}
|
||||
// The curve at a radius, to the nearest bin, and the same in whole 1/BRAGG_CURVE_FIXED counts.
|
||||
[[nodiscard]] float CurveAt(float r) const {
|
||||
const int n = static_cast<int>(radial_curve.size());
|
||||
const int b = std::clamp(static_cast<int>(std::lround(r0 + rad)), 0, n - 1);
|
||||
const int b0 = std::clamp(static_cast<int>(std::lround(r0)), 0, n - 1);
|
||||
return static_cast<int32_t>(std::lround(radial_curve[b] - radial_curve[b0]));
|
||||
return radial_curve[std::clamp(static_cast<int>(std::lround(r)), 0, n - 1)];
|
||||
}
|
||||
[[nodiscard]] int64_t CurveFixed(float r) const {
|
||||
return static_cast<int64_t>(CurveAt(r) * BRAGG_CURVE_FIXED);
|
||||
}
|
||||
[[nodiscard]] const std::vector<float> &GetRadialBackground() const { return radial_curve; }
|
||||
|
||||
|
||||
@@ -259,6 +259,12 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
|
||||
BraggBoxSums s{};
|
||||
bkg_vals.clear();
|
||||
clip_ring.clear();
|
||||
// Where the frame's curve is in use: the reflection's place on it, and the curve summed over the
|
||||
// readable disk pixels and over the ring pixels the background is taken from
|
||||
// (BraggCurveCorrection).
|
||||
const bool curve = bkg_radial && !radial_curve.empty();
|
||||
const float rc = curve ? CurveRadius(pn) : st.r0;
|
||||
int64_t curve_disk = 0, curve_ring = 0, curve_kept = 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);
|
||||
@@ -292,6 +298,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
|
||||
s.x_sum += x;
|
||||
s.y_sum += y;
|
||||
++s.n_inner_valid;
|
||||
if (curve) curve_disk += CurveFixed(rc + d.rad);
|
||||
} else if (d.inner >= r2_sq && d.outer < r3_sq) {
|
||||
if (!valid(px)) continue;
|
||||
++s.n_ring;
|
||||
@@ -301,6 +308,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
|
||||
s.bkg_sum += px;
|
||||
if (bkg_trim > 0.0f) bkg_vals.push_back(px);
|
||||
++s.n_bkg;
|
||||
if (curve) curve_ring += CurveFixed(rc + d.rad);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -310,9 +318,10 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
|
||||
if (ClipsRing(s)) {
|
||||
const int32_t thr = ClipThreshold(s);
|
||||
for (const auto &[px, rad] : clip_ring)
|
||||
if (static_cast<int64_t>(px) <= static_cast<int64_t>(thr) + RadialClipOffset(st.r0, rad)) {
|
||||
if (static_cast<int64_t>(px) <= static_cast<int64_t>(thr) + RadialClipOffset(rc, rad)) {
|
||||
clip.sum += px;
|
||||
++clip.n;
|
||||
if (curve) curve_kept += CurveFixed(rc + rad);
|
||||
if (bkg_radial && radial_curve.empty()) {
|
||||
// The radial curve is binned on the TRUE detector radius, so the
|
||||
// offset here is the unshrunk radial projection.
|
||||
@@ -323,6 +332,11 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
|
||||
}
|
||||
}
|
||||
FinishPassA(a, i, r, st, s, clip, bkg_vals);
|
||||
if (curve) {
|
||||
const bool kept = ClipsRing(s) && clip.n > 5; // the pixels FinishPassA took the background from
|
||||
a.rough[i].curve_corr = BraggCurveCorrection(kept ? curve_kept : curve_ring, kept ? clip.n : s.n_bkg,
|
||||
curve_disk, static_cast<int>(s.n_inner_valid));
|
||||
}
|
||||
}
|
||||
|
||||
return ProfilePasses(img, owner, predicted, npredicted, panels, a, image_number);
|
||||
@@ -609,24 +623,23 @@ std::vector<Reflection> BraggIntegrationEngineCPU::ProfilePasses(const Sampler &
|
||||
|
||||
// --- 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. ---
|
||||
// correct it by mean_ring(B) - mean_disk(B). With the frame's curve that was taken in pass A over
|
||||
// the pixels themselves (BraggCurveCorrection); with the curve just accumulated from the
|
||||
// reflections it is the kernel, 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, a.n_rad - 1);
|
||||
double v;
|
||||
if (!radial_curve.empty())
|
||||
v = radial_curve[b];
|
||||
else
|
||||
v = a.rad_cnt[b] > 0 ? static_cast<double>(a.rad_sum[b]) / a.rad_cnt[b] : rh.bkg;
|
||||
corr += static_cast<double>(kern[k]) * v;
|
||||
double corr = rh.curve_corr;
|
||||
if (radial_curve.empty()) {
|
||||
const float *kern = k_diff.data() + static_cast<size_t>(rh.k_bin) * k_len;
|
||||
for (int k = 0; k < k_len; ++k) {
|
||||
const int b = std::clamp(rh.r_bin + k - k_off, 0, a.n_rad - 1);
|
||||
const double v = a.rad_cnt[b] > 0 ? static_cast<double>(a.rad_sum[b]) / a.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;
|
||||
|
||||
@@ -105,6 +105,7 @@ protected:
|
||||
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
|
||||
double curve_corr = 0.0; // the correction with the frame's curve (BraggCurveCorrection)
|
||||
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
|
||||
|
||||
@@ -41,6 +41,12 @@ struct BraggGpuParams {
|
||||
|
||||
__device__ inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MAX; }
|
||||
|
||||
// The frame's curve at a radius, to the nearest bin, in whole 1/BRAGG_CURVE_FIXED counts
|
||||
// (BraggIntegrationEngine::CurveFixed).
|
||||
__device__ inline long long curve_fixed(const float *curve, int n, float r) {
|
||||
return (long long) (curve[min(max((int) lroundf(r), 0), n - 1)] * BRAGG_CURVE_FIXED);
|
||||
}
|
||||
|
||||
// Learned second moments, (sum v*rad^2, sum v*tan^2, sum v) per shell plus one global slot at N_SHELL.
|
||||
constexpr int MOM_STRIDE = 3;
|
||||
|
||||
@@ -169,7 +175,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
|
||||
uint8_t *hasobs_o, uint8_t *overloaded_o, int *saturated_o, unsigned long long *invd2mm,
|
||||
float *isum_o, int *ninner_o, int *rbin_o, int *kbin_o,
|
||||
unsigned long long *rad_sum, int *rad_cnt, int n_rad, int accumulate_rad,
|
||||
const float *rad_curve, unsigned long long *counts) {
|
||||
const float *rad_curve, const float *curve_r, float *curvecorr_o,
|
||||
unsigned long long *counts) {
|
||||
const int i = blockIdx.x;
|
||||
if (i >= n) return;
|
||||
|
||||
@@ -184,6 +191,11 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
|
||||
__shared__ long long s_thr_i;
|
||||
__shared__ int s_clipn;
|
||||
__shared__ float s_r0;
|
||||
// Where the frame's curve is in use: the reflection's place on it (CurveRadius), and the curve
|
||||
// summed over the readable disk pixels, over the ring pixels behind the background, and over those
|
||||
// the clip keeps - in whole 1/BRAGG_CURVE_FIXED counts, so the sums are exact (BraggCurveCorrection).
|
||||
__shared__ float s_rc;
|
||||
__shared__ unsigned long long s_cdisk, s_cring, s_ckept;
|
||||
// The stencil spans r0 +- the radial semi-axis of the outer ellipse, so the window has to be
|
||||
// sized from the widest aperture on the detector - which the host does, into p.rad_w. A fixed
|
||||
// 32-bin window silently dropped the outer bins as soon as anything was elongated.
|
||||
@@ -200,6 +212,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
|
||||
if (threadIdx.x == 0) {
|
||||
const float rx = pn.plane_x - p.beam_x, ry = pn.plane_y - p.beam_y;
|
||||
s_r0 = sqrtf(rx * rx + ry * ry);
|
||||
s_rc = rad_curve ? curve_r[i] : s_r0;
|
||||
s_cdisk = 0; s_cring = 0; s_ckept = 0;
|
||||
s_Isum = 0; s_Ix = 0; s_Iy = 0; s_x = 0; s_y = 0;
|
||||
s_ninner = 0; s_ninner_valid = 0; s_nbkg = 0; s_bkgsum = 0;
|
||||
s_ndisk = 0; s_nown = 0; s_nring = 0; s_nclean = 0; s_nsat = 0; s_nmax = 0;
|
||||
@@ -221,6 +235,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
|
||||
long long l_x = 0, l_y = 0; // positions behind l_Isum, for the background-free centroid
|
||||
int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0, l_nring = 0, l_nclean = 0, l_nsat = 0, l_nmax = 0;
|
||||
long long l_bkg = 0;
|
||||
long long l_cdisk = 0, l_cring = 0;
|
||||
for (int t = threadIdx.x; t < area; t += blockDim.x) {
|
||||
const int x = x0 + t % bw, y = y0 + t / bw;
|
||||
const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy);
|
||||
@@ -242,7 +257,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
|
||||
if (px == INT32_MAX || (px == INT32_MIN && on && !static_mask[y * p.W + x])) ++l_nsat;
|
||||
if (px == INT32_MAX) ++l_nmax;
|
||||
if (valid(px)) { l_Isum += px; l_Ix += (long long) x * px; l_Iy += (long long) y * px;
|
||||
l_x += x; l_y += y; ++l_niv; }
|
||||
l_x += x; l_y += y; ++l_niv;
|
||||
if (rad_curve) l_cdisk += curve_fixed(rad_curve, n_rad, s_rc + d.rad); }
|
||||
} else if (d.inner >= p.r2_sq && d.outer < p.r3_sq) {
|
||||
if (!BraggPanelHolds(pn, x, y)) continue;
|
||||
const int32_t px = img[y * p.W + x];
|
||||
@@ -252,6 +268,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
|
||||
if (m == 0) ++l_nclean;
|
||||
if (m == MASK_FLUX) continue;
|
||||
l_bkg += px; ++l_nb;
|
||||
if (rad_curve) l_cring += curve_fixed(rad_curve, n_rad, s_rc + d.rad);
|
||||
}
|
||||
}
|
||||
WARP_ATOMIC_ADD(s_Isum, (unsigned long long) l_Isum);
|
||||
@@ -260,6 +277,10 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
|
||||
WARP_ATOMIC_ADD(s_x, (unsigned long long) l_x);
|
||||
WARP_ATOMIC_ADD(s_y, (unsigned long long) l_y);
|
||||
WARP_ATOMIC_ADD(s_ninner, l_ni);
|
||||
if (rad_curve) {
|
||||
WARP_ATOMIC_ADD(s_cdisk, (unsigned long long) l_cdisk);
|
||||
WARP_ATOMIC_ADD(s_cring, (unsigned long long) l_cring);
|
||||
}
|
||||
WARP_ATOMIC_ADD(s_ninner_valid, l_niv);
|
||||
WARP_ATOMIC_ADD(s_nbkg, l_nb);
|
||||
WARP_ATOMIC_ADD(s_bkgsum, (unsigned long long) l_bkg);
|
||||
@@ -352,7 +373,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
|
||||
|
||||
// Second ring pass for the high-side sigma-clip (re-reads the annulus; avoids storing bkg values).
|
||||
if (s_accept && p.bkg_clip_nsigma > 0.0f && !do_trim) {
|
||||
long long c_l = 0; int cn_l = 0;
|
||||
long long c_l = 0, cc_l = 0; int cn_l = 0;
|
||||
for (int t = threadIdx.x; t < area; t += blockDim.x) {
|
||||
const int x = x0 + t % bw, y = y0 + t / bw;
|
||||
const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy);
|
||||
@@ -363,12 +384,13 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
|
||||
// Where the frame's radial curve is in use the clip follows it (RadialClipOffset).
|
||||
long long lim = s_thr_i;
|
||||
if (rad_curve) {
|
||||
const int b = min(max((int) lroundf(s_r0 + d.rad), 0), n_rad - 1);
|
||||
const int b0 = min(max((int) lroundf(s_r0), 0), n_rad - 1);
|
||||
const int b = min(max((int) lroundf(s_rc + d.rad), 0), n_rad - 1);
|
||||
const int b0 = min(max((int) lroundf(s_rc), 0), n_rad - 1);
|
||||
lim += lroundf(rad_curve[b] - rad_curve[b0]);
|
||||
}
|
||||
if ((long long) px <= lim) {
|
||||
c_l += px; ++cn_l;
|
||||
if (rad_curve) cc_l += curve_fixed(rad_curve, n_rad, s_rc + d.rad);
|
||||
if (accumulate_rad) {
|
||||
// The radial curve is binned on the TRUE detector radius, so this is the
|
||||
// unshrunk radial projection - no per-pixel sqrt.
|
||||
@@ -381,6 +403,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
|
||||
}
|
||||
}
|
||||
WARP_ATOMIC_ADD(s_clipsum, (unsigned long long) c_l); WARP_ATOMIC_ADD(s_clipn, cn_l);
|
||||
if (rad_curve) WARP_ATOMIC_ADD(s_ckept, (unsigned long long) cc_l);
|
||||
}
|
||||
__syncthreads();
|
||||
|
||||
@@ -445,6 +468,12 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
|
||||
isum_o[i] = (float) Isum;
|
||||
ninner_o[i] = s_ninner_valid;
|
||||
rbin_o[i] = min(max((int) lroundf(s_r0), 0), n_rad > 0 ? n_rad - 1 : 0);
|
||||
if (rad_curve) {
|
||||
// From the ring pixels the background above was taken from (see the CPU engine).
|
||||
const bool kept = p.bkg_clip_nsigma > 0.0f && s_clipn > 5;
|
||||
curvecorr_o[i] = (float) BraggCurveCorrection((long long) (kept ? s_ckept : s_cring), kept ? s_clipn : s_nbkg,
|
||||
(long long) s_cdisk, s_ninner_valid);
|
||||
}
|
||||
kbin_o[i] = BraggStencilKernelIndex(st, p.n_kern);
|
||||
obsx_o[i] = (float) ox; obsy_o[i] = (float) oy; hasobs_o[i] = hasobs;
|
||||
overloaded_o[i] = s_nsat > 0 ? 1 : 0;
|
||||
@@ -610,9 +639,11 @@ __global__ void build_profiles(const unsigned long long *shell_grid, const unsig
|
||||
// histogram over radial offset) turns that into one short dot product. An empty radial bin
|
||||
// contributes the reflection's own background, so an empty neighbourhood gives exactly zero
|
||||
// correction because the kernel weights sum to zero. The kernel depends on how far this
|
||||
// reflection's ring was elongated, so k_diff is a table and kbin_a picks the row.
|
||||
// reflection's ring was elongated, so k_diff is a table and kbin_a picks the row. With the
|
||||
// frame's curve handed over, boxsum already took the correction over the pixels themselves
|
||||
// (curve_corr, BraggCurveCorrection) and only its application is left here.
|
||||
// Mirrors BraggIntegrationEngineCPU. ---
|
||||
__global__ void radial_correct(const unsigned long long *rad_sum, const int *rad_cnt, const float *curve, int n_rad,
|
||||
__global__ void radial_correct(const unsigned long long *rad_sum, const int *rad_cnt, const float *curve_corr, int n_rad,
|
||||
const float *k_diff, int k_len, int k_off,
|
||||
const float *isum_a, const int *ninner_a, const int *rbin_a,
|
||||
const int *kbin_a, float pixel_gain,
|
||||
@@ -620,14 +651,18 @@ __global__ void radial_correct(const unsigned long long *rad_sum, const int *rad
|
||||
const int i = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
if (i >= n || !ok_a[i]) return;
|
||||
const float bkg = bkg_o[i];
|
||||
const float *kern = k_diff + (size_t) kbin_a[i] * k_len;
|
||||
// With the frame's curve the correction was taken over the pixels themselves in boxsum.
|
||||
float corr = 0.0f;
|
||||
for (int k = 0; k < k_len; ++k) {
|
||||
int b = rbin_a[i] + k - k_off;
|
||||
b = min(max(b, 0), n_rad - 1);
|
||||
const float v = curve ? curve[b]
|
||||
: rad_cnt[b] > 0 ? (float) (long long) rad_sum[b] / (float) rad_cnt[b] : bkg;
|
||||
corr += kern[k] * v;
|
||||
if (curve_corr)
|
||||
corr = curve_corr[i];
|
||||
else {
|
||||
const float *kern = k_diff + (size_t) kbin_a[i] * k_len;
|
||||
for (int k = 0; k < k_len; ++k) {
|
||||
int b = rbin_a[i] + k - k_off;
|
||||
b = min(max(b, 0), n_rad - 1);
|
||||
const float v = rad_cnt[b] > 0 ? (float) (long long) rad_sum[b] / (float) rad_cnt[b] : bkg;
|
||||
corr += kern[k] * v;
|
||||
}
|
||||
}
|
||||
const float bkg_new = bkg - corr;
|
||||
bkg_o[i] = bkg_new;
|
||||
@@ -952,6 +987,8 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) {
|
||||
d_isum = CudaDevicePtr<float>(new_capacity);
|
||||
d_ninner = CudaDevicePtr<int>(new_capacity);
|
||||
d_rbin = CudaDevicePtr<int>(new_capacity);
|
||||
d_curve_corr = CudaDevicePtr<float>(new_capacity);
|
||||
d_curve_r = CudaDevicePtr<float>(new_capacity);
|
||||
d_kbin = CudaDevicePtr<int>(new_capacity);
|
||||
d_obs_x = CudaDevicePtr<float>(new_capacity);
|
||||
d_obs_y = CudaDevicePtr<float>(new_capacity);
|
||||
@@ -964,6 +1001,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) {
|
||||
h_px_x = CudaHostPtr<float>(new_capacity);
|
||||
h_px_y = CudaHostPtr<float>(new_capacity);
|
||||
h_d = CudaHostPtr<float>(new_capacity);
|
||||
h_curve_r = CudaHostPtr<float>(new_capacity);
|
||||
h_panel = CudaHostPtr<BraggPanel>(new_capacity);
|
||||
h_mark = CudaHostPtr<uint8_t>(new_capacity);
|
||||
h_I = CudaHostPtr<float>(new_capacity);
|
||||
@@ -1007,6 +1045,12 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
|
||||
cuda_err(cudaMemcpyAsync(d_d, h_d.get(), sizeof(float) * npredicted, cudaMemcpyHostToDevice, *stream));
|
||||
cuda_err(cudaMemcpyAsync(d_panel, h_panel.get(), sizeof(BraggPanel) * npredicted, cudaMemcpyHostToDevice, *stream));
|
||||
cuda_err(cudaMemcpyAsync(d_mark, h_mark.get(), sizeof(uint8_t) * npredicted, cudaMemcpyHostToDevice, *stream));
|
||||
// Where the frame's curve is in use, where each reflection sits on it (CurveRadius).
|
||||
if (bkg_radial && !radial_curve.empty()) {
|
||||
for (size_t i = 0; i < npredicted; ++i)
|
||||
h_curve_r[i] = CurveRadius(panels[i]);
|
||||
cuda_err(cudaMemcpyAsync(d_curve_r, h_curve_r.get(), sizeof(float) * npredicted, cudaMemcpyHostToDevice, *stream));
|
||||
}
|
||||
|
||||
BraggGpuParams p{
|
||||
.W = static_cast<int>(xpixel), .H = static_cast<int>(ypixel),
|
||||
@@ -1074,14 +1118,14 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
|
||||
d_ok, d_strong, d_has_obs, d_overloaded, d_saturated, d_invd2,
|
||||
d_isum, d_ninner, d_rbin, d_kbin,
|
||||
d_rad_sum, d_rad_cnt, rad_n, rad_n > 0 && !radial_given,
|
||||
radial_given ? d_rad_curve.get() : nullptr, d_counts);
|
||||
radial_given ? d_rad_curve.get() : nullptr, d_curve_r, d_curve_corr, d_counts);
|
||||
cuda_err(cudaGetLastError());
|
||||
|
||||
// Correct the flat annulus background for the curvature of the radial background before anything
|
||||
// downstream (profile fit, variance) reads it.
|
||||
if (radial_on) {
|
||||
radial_correct<<<(n + threads - 1) / threads, threads, 0, *stream>>>(
|
||||
d_rad_sum, d_rad_cnt, radial_given ? d_rad_curve.get() : nullptr, n_rad, d_k_diff, k_len, k_off,
|
||||
d_rad_sum, d_rad_cnt, radial_given ? d_curve_corr.get() : nullptr, n_rad, d_k_diff, k_len, k_off,
|
||||
d_isum, d_ninner, d_rbin, d_kbin, static_cast<float>(pixel_gain), d_ok, d_bkg, d_I, d_var_bkg, n);
|
||||
cuda_err(cudaGetLastError());
|
||||
}
|
||||
|
||||
@@ -56,6 +56,8 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine {
|
||||
CudaDevicePtr<float> d_k_diff;
|
||||
CudaDevicePtr<int> d_rad_cnt;
|
||||
CudaDevicePtr<float> d_rad_curve; // the curve handed over (RadialBackgroundFromRings)
|
||||
CudaDevicePtr<float> d_curve_corr; // its correction per reflection (BraggCurveCorrection)
|
||||
CudaDevicePtr<float> d_curve_r; // where each reflection sits on it (CurveRadius)
|
||||
|
||||
// --- fixed-size device arrays ---
|
||||
// The learning/fit math is single precision: FP64 is heavily throttled on consumer GPUs and the
|
||||
@@ -81,6 +83,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine {
|
||||
// driver has staged it through a bounce buffer, so eight of them in a row are eight serialised
|
||||
// round-trips rather than eight queued transfers.
|
||||
CudaHostPtr<float> h_px_x, h_px_y, h_d;
|
||||
CudaHostPtr<float> h_curve_r; // where each reflection sits on the frame's curve (CurveRadius)
|
||||
CudaHostPtr<BraggPanel> h_panel;
|
||||
CudaHostPtr<uint8_t> h_mark;
|
||||
CudaHostPtr<float> h_I, h_sigma, h_bkg, h_var_bkg, h_obs_x, h_obs_y;
|
||||
|
||||
@@ -290,6 +290,25 @@ struct BraggStencilDist {
|
||||
float signal, inner, outer, rad;
|
||||
};
|
||||
|
||||
// --- The frame's radial background curve (BraggIntegrationEngine::RadialBackgroundFromRings) -------
|
||||
// The curve is summed over a reflection's pixels in whole 1/BRAGG_CURVE_FIXED counts: integer sums,
|
||||
// so they do not depend on the order the pixels are visited in, as the engines' other sums do not.
|
||||
BRAGG_STENCIL_CONSTANT constexpr float BRAGG_CURVE_FIXED = 256.0f;
|
||||
|
||||
#ifndef __METAL_VERSION__
|
||||
// The radial background correction with the frame's curve: the curve's mean over the ring pixels the
|
||||
// background was taken from, minus its mean over the disk pixels the intensity was summed over, from
|
||||
// their fixed-point sums. It is what the ring mean got wrong about the disk if the background is the
|
||||
// curve - taken over the pixels actually read, so it holds whatever the ring lost (the clip, a
|
||||
// neighbour's signal, a masked or shadowed pixel), where a kernel averaged over a whole annulus does
|
||||
// not. Nothing to correct on a disk with no readable pixel.
|
||||
BRAGG_STENCIL_HD double BraggCurveCorrection(long long ring_sum, int ring_n, long long disk_sum, int disk_n) {
|
||||
if (ring_n <= 0 || disk_n <= 0)
|
||||
return 0.0;
|
||||
return ((double) ring_sum / ring_n - (double) disk_sum / disk_n) / BRAGG_CURVE_FIXED;
|
||||
}
|
||||
#endif
|
||||
|
||||
// --- Ownership of a shared signal region ---------------------------------------------------------
|
||||
// Two reflections whose centres are closer than 2*r1 have signal disks that intersect, and a pixel in
|
||||
// the intersection carries both. It belongs to the NEARER centre. The union mask above cannot say
|
||||
|
||||
@@ -181,7 +181,7 @@ TEST_CASE("BraggIntegrationEngineCPU_RadialCurveFromRings", "[Integration][porta
|
||||
r.h = 1; r.k = 0; r.l = 0;
|
||||
r.predicted_x = bx + 300.0f;
|
||||
r.predicted_y = by;
|
||||
r.d = 2.0f;
|
||||
r.d = static_cast<float>(lambda / (2.0 * std::sin(0.5 * std::atan(300.0 / F_px)))); // where the curve has it
|
||||
r.prescaling_corr = 1.0f;
|
||||
r.partiality = 1.0f;
|
||||
const std::vector<Reflection> predicted{r};
|
||||
@@ -235,7 +235,7 @@ TEST_CASE("BraggIntegrationEngineCPU_ClipFollowsRadialCurve", "[Integration][por
|
||||
r.h = 1; r.k = 0; r.l = 0;
|
||||
r.predicted_x = bx + 200.0f;
|
||||
r.predicted_y = by;
|
||||
r.d = 2.0f;
|
||||
r.d = static_cast<float>(lambda / (2.0 * std::sin(0.5 * std::atan(200.0 / F_px))));
|
||||
r.prescaling_corr = 1.0f;
|
||||
r.partiality = 1.0f;
|
||||
const std::vector<Reflection> predicted{r};
|
||||
@@ -254,6 +254,77 @@ TEST_CASE("BraggIntegrationEngineCPU_ClipFollowsRadialCurve", "[Integration][por
|
||||
CHECK(std::fabs(background_read(true) - truth) < 20.0);
|
||||
}
|
||||
|
||||
// An empty aperture next to the edge of the beam-stop shadow: the background rises towards the beam until
|
||||
// the shadow cuts it off, and the shadowed pixels are masked, so the ring has lost its beam-side part. The
|
||||
// correction is taken over the ring pixels actually read (BraggCurveCorrection) and recovers the
|
||||
// background under the disk (51.7 against 51.5; the plain ring reads 54.3). An average of the curve over
|
||||
// the whole annulus counts the curve inside the shadow, where the ring has no pixels, and on the same
|
||||
// edge without the tilt moved the background the wrong way (57 against 50.4). The detector is tilted, so
|
||||
// the curve has to be read at the reflection's own 2 theta (CurveRadius): the PONI is 17 px from where
|
||||
// the beam lands, and read at the distance from it the reflection sits on the flat part of the curve and
|
||||
// nothing is corrected (54.3).
|
||||
TEST_CASE("BraggIntegrationEngineCPU_RadialCurveOverThePixelsRead", "[Integration][portable]") {
|
||||
DiffractionExperiment experiment(DetJF(2));
|
||||
experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f)
|
||||
.PoniRot1_rad(0.01f);
|
||||
const DiffractionGeometry geom = experiment.GetDiffractionGeometry();
|
||||
const double F_px = 100.0 / 0.075, lambda = 1.0;
|
||||
// The radius the curve is tabulated against, F tan(2 theta), of a pixel.
|
||||
const auto radius = [&](double x, double y) {
|
||||
return F_px * std::tan(2.0 * std::asin(lambda / (2.0 * geom.PxlToRes(static_cast<float>(x), static_cast<float>(y)))));
|
||||
};
|
||||
const double edge = 296.0; // the shadow's edge
|
||||
const auto background = [](double r) { return 50.0 + 4.0 * std::max(0.0, 302.0 - r); };
|
||||
const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum();
|
||||
ImagePreprocessorBuffer image(npixel);
|
||||
for (size_t i = 0; i < npixel; ++i) {
|
||||
const double r = radius(static_cast<double>(i % width), static_cast<double>(i / width));
|
||||
image[i] = r < edge ? INT32_MIN : static_cast<int32_t>(std::lround(background(r)));
|
||||
}
|
||||
// The rings inside the shadow hold what is left there, little; fine enough that the edge is sharp.
|
||||
const float low_q = 0.05f, dq = 0.0002f;
|
||||
std::vector<float> rings(10000);
|
||||
for (size_t i = 0; i < rings.size(); ++i) {
|
||||
const double q = low_q + (i + 0.5) * dq;
|
||||
const double r = F_px * std::tan(2.0 * std::asin(q * lambda / (4.0 * PI)));
|
||||
rings[i] = static_cast<float>(r < edge ? 10.0 : background(r));
|
||||
}
|
||||
// 304 px out from where the beam lands, on the side away from the PONI.
|
||||
const auto [dbx, dby] = geom.GetDirectBeam_pxl();
|
||||
REQUIRE(std::fabs(dbx - 400.0f) > 10.0f);
|
||||
Reflection r{};
|
||||
r.h = 1; r.k = 0; r.l = 0;
|
||||
r.predicted_x = std::round(dbx + (dbx > 400.0f ? 304.0f : -304.0f));
|
||||
r.predicted_y = std::round(dby);
|
||||
r.d = geom.PxlToRes(r.predicted_x, r.predicted_y);
|
||||
r.prescaling_corr = 1.0f;
|
||||
r.partiality = 1.0f;
|
||||
const std::vector<Reflection> predicted{r};
|
||||
auto integrate = [&](bool correct) {
|
||||
experiment.ImportBraggIntegrationSettings(
|
||||
BraggIntegrationSettings().Integrator(IntegratorMode::BoxSum).BackgroundRadialCorrection(correct));
|
||||
BraggIntegrationEngineCPU engine(experiment, PixelMask(experiment));
|
||||
engine.RadialBackgroundFromRings(rings, static_cast<int32_t>(rings.size()), low_q, dq);
|
||||
const auto out = engine.Run(image, predicted, predicted.size(), 0);
|
||||
REQUIRE(out.size() == 1);
|
||||
return out[0];
|
||||
};
|
||||
// The background under the disk, from the image itself.
|
||||
const float r1 = experiment.GetBraggIntegrationSettings().GetR1();
|
||||
double truth = 0.0;
|
||||
int n = 0;
|
||||
for (int dy = -8; dy <= 8; ++dy)
|
||||
for (int dx = -8; dx <= 8; ++dx)
|
||||
if (dx * dx + dy * dy < r1 * r1) {
|
||||
truth += image[static_cast<size_t>(r.predicted_y + dy) * width + static_cast<size_t>(r.predicted_x + dx)];
|
||||
++n;
|
||||
}
|
||||
truth /= n;
|
||||
const Reflection plain = integrate(false), corrected = integrate(true);
|
||||
CHECK(plain.bkg > truth + 2.0);
|
||||
CHECK(corrected.bkg == Catch::Approx(truth).margin(0.5));
|
||||
}
|
||||
|
||||
// A still has no rocking event to drop, so an overloaded observation is kept as a lower bound: the box sum
|
||||
// with every saturated pixel at the saturation limit. A spot with its core clipped at the limit, the
|
||||
// clipped pixels marked saturated as the preprocessor marks them: the still gives exactly the clipped
|
||||
|
||||
@@ -496,8 +496,8 @@ TEST_CASE("BraggIntegrationEngineGPU_RadialCurveFromRingsMatchesCPU", "[Integrat
|
||||
ImagePreprocessorBuffer cpu_image(npixel);
|
||||
ImagePreprocessorBufferGPU gpu_image(npixel);
|
||||
for (size_t i = 0; i < npixel; ++i) {
|
||||
cpu_image[i] = static_cast<int32_t>(std::lround(
|
||||
background(std::hypot(static_cast<double>(i % width) - bx, static_cast<double>(i / width) - by))));
|
||||
const double r = std::hypot(static_cast<double>(i % width) - bx, static_cast<double>(i / width) - by);
|
||||
cpu_image[i] = r < 296.0 ? INT32_MIN : static_cast<int32_t>(std::lround(background(r)));
|
||||
gpu_image[i] = cpu_image[i];
|
||||
}
|
||||
auto stream = std::make_shared<CudaStream>();
|
||||
@@ -512,9 +512,14 @@ TEST_CASE("BraggIntegrationEngineGPU_RadialCurveFromRingsMatchesCPU", "[Integrat
|
||||
rings[i] = static_cast<float>(background(F_px * std::tan(2.0 * std::asin(q * lambda / (4.0 * PI)))));
|
||||
}
|
||||
|
||||
// Half of them next to a masked disc around the beam - a beam-stop shadow - whose edge their rings
|
||||
// reach, so the correction is taken over a ring that lost part of itself (BraggCurveCorrection).
|
||||
std::vector<Reflection> predicted;
|
||||
for (int i = 0; i < 8; ++i)
|
||||
predicted.push_back(MakeReflection(bx + 300.0f * std::cos(0.7f * i), by + 300.0f * std::sin(0.7f * i), 2.0f, i + 1));
|
||||
for (int i = 0; i < 8; ++i) {
|
||||
const float rr = i % 2 ? 306.0f : 300.0f;
|
||||
const float d = static_cast<float>(lambda / (2.0 * std::sin(0.5 * std::atan(rr / F_px))));
|
||||
predicted.push_back(MakeReflection(bx + rr * std::cos(0.7f * i), by + rr * std::sin(0.7f * i), d, i + 1));
|
||||
}
|
||||
|
||||
for (const auto mode : {IntegratorMode::BoxSum, IntegratorMode::ProfileGaussian}) {
|
||||
experiment.ImportBraggIntegrationSettings(BraggIntegrationSettings().Integrator(mode).BackgroundRadialCorrection(true));
|
||||
|
||||
Reference in new issue
Block a user