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:
1 parent
f199f23c54
commit
fac2a17428
11 files changed
+225
-51
No files matched your search
@@ -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
|
||||
|
||||
Reference in new issue
Block a user