From 4d3434e2a59ef48cc7c8f0d959c45a51d2db9418 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 9 Aug 2026 05:31:57 +0200 Subject: [PATCH] Beam stop: compare each pixel only against its own ring The background belongs to the beam and the shadow to the stop, and the two are not concentric - fitting the stop edge per azimuth gives offsets of 13.4 px on an 85.8 px disk, 22.2 px on 67.3 px and 6.9 px on 23.7 px, 8 to 33 per cent of the stop radius on every crystal measured. The finder bridged that gap with a radial envelope, the largest ring background over an outward window, used as the reference for an individual pixel. That quantity exceeds the local background wherever the background rises outward, so sound pixels near the stop scored below the penumbra threshold and were masked. Measured against the fitted edge on a long-distance disk stop, the mask was displaced rather than mis-sized: short by up to 20 px on one side, over-reaching by up to 45 px on the other, with eight of twenty-four azimuth sectors falling short. The ring median is already the right reference wherever a ring still has unshadowed pixels to measure, which is every ring except those lying wholly inside the disk - and it needs no assumption about where the stop sits. So the envelope is gone from the per-pixel test, and the rings it existed to cover are handled directly: walking outward, a ring whose background is a fraction of the background further out is shadow in its entirety. That comparison is only ever asked whether a whole ring is inside the stop, never to judge a pixel, which is where its failure mode lives. Blockage is deliberately not a counting test - on a bright dataset the shadow interior is still well counted. Detection is now one channel instead of two, and 113 lines shorter. Measured: no azimuth sector falls short by more than 3.4 px, over-reach drops on all three fitted crystals, and mask area moves by at most 0.04 per cent of the detector on six crystals, so this corrects the shape rather than resizing. Battery: space-group agreement with XDS unchanged at 34/37, median change in R_meas and in the lowest shell 0.000 pp. The crystal that suffered worst when masking was introduced recovers to its unmasked quality - R_meas 25.1 -> 17.2 per cent, ISa 4.45 -> 10.04 - which is what removing the over-masking should do. Co-Authored-By: Claude Opus 5 (1M context) --- image_analysis/beam_stop/ShadowFinder.cpp | 70 +++++++++++++---------- image_analysis/beam_stop/ShadowFinder.h | 13 +++-- 2 files changed, 49 insertions(+), 34 deletions(-) diff --git a/image_analysis/beam_stop/ShadowFinder.cpp b/image_analysis/beam_stop/ShadowFinder.cpp index 4687ce0c..93928281 100644 --- a/image_analysis/beam_stop/ShadowFinder.cpp +++ b/image_analysis/beam_stop/ShadowFinder.cpp @@ -32,10 +32,13 @@ constexpr int64_t MIN_REFLECTION = 25; // and testing anyway masks whole detectors on low-background data. constexpr double MIN_EXPECTED_COUNTS = 60; -// Side of the box the background is pooled over before testing, and how far out the radial -// comparison looks for unshadowed background. -constexpr int POOL_PX = 5; -constexpr float ENVELOPE_MM = 6.0f; +// Side of the box the background is pooled over before testing. Its area is how many pixels back +// a ring's countability test, which decides where an azimuthal comparison is possible at all. +constexpr int POOL_PX = 5; +constexpr double MEAN_POOLED_PIXELS = POOL_PX * POOL_PX; + +// A ring with fewer valid pixels than this says nothing about whether it was counted. +constexpr int MIN_RING_PIXELS = 32; // Binary-image helpers on a width*height frame stored row-major as char (0/1). All run once, // at GetMask() time; the BFS forms keep them O(pixels) rather than O(pixels * radius). @@ -162,28 +165,15 @@ std::vector ring_median(const std::vector &values, const std::vect for (size_t i = 0; i < values.size(); i++) if (use[i]) bins[radius[i]].push_back(values[i]); - std::vector median(max_radius + 1, 0.0f); + std::vector out(max_radius + 1, 0.0f); for (int r = 0; r <= max_radius; r++) { auto &b = bins[r]; if (!b.empty()) { const size_t k = b.size() / 2; std::nth_element(b.begin(), b.begin() + k, b.end()); - median[r] = b[k]; + out[r] = b[k]; } } - return median; -} - -// Largest baseline over [r, r + win] - the background just outside radius r. -std::vector outer_envelope(const std::vector &baseline, int win) { - const int n = static_cast(baseline.size()); - std::vector out(n, 0.0f); - for (int r = 0; r < n; r++) { - float v = baseline[r]; - for (int k = 1; k <= win; k++) - v = std::max(v, baseline[std::min(r + k, n - 1)]); - out[r] = v; - } return out; } @@ -194,7 +184,6 @@ ShadowFinder::ShadowFinder(const DiffractionExperiment &experiment, const PixelM height(static_cast(experiment.GetYPixelsNumConv())), beam_x(experiment.GetBeamX_pxl()), beam_y(experiment.GetBeamY_pxl()), - envelope_px(std::max(4, static_cast(std::lround(ENVELOPE_MM / experiment.GetPixelSize_mm())))), pixel_mask(mask.GetMask(experiment)), max_value(static_cast(width) * height, 0), sum_value(static_cast(width) * height, 0), @@ -314,22 +303,45 @@ std::vector ShadowFinder::GetMask() const { excluded[i] = valid[i] && ratio[i] < SHADOW_RATIO; } - // Radial comparison: the background just outside this radius. The disk blocks its rings - // completely, so their median is the shadow itself and only this comparison sees it. - const auto envelope = outer_envelope(baseline, envelope_px); - std::vector ratio_radial(n_pixels, 1.0f); + // A ring whose background was never counted carries no information to test a pixel against. + // Walking outward, every ring before the first countable one lies wholly inside the stop - a + // ring fully within the disk has no unshadowed pixel for the median to find, which is exactly + // where an azimuthal comparison must fail. Those rings are shadow in their entirety. + // Innermost rings hold only a handful of pixels, too few to judge, so they are stepped over + // rather than allowed to end the walk. + std::vector ring_pixels(max_radius + 1, 0); for (int i = 0; i < n_pixels; i++) if (valid[i]) - ratio_radial[i] = pooled[i] / std::max(envelope[radius[i]], 1e-6f); + ring_pixels[radius[i]]++; + + // A ring lies inside the stop when its background is a fraction of the background further out. + // Counting statistics cannot decide this: on a bright dataset the shadow is still well counted. + // The comparison is only ever used to answer "is this whole ring inside the stop", never to + // judge an individual pixel, so taking the largest background over an outward window is safe + // here in a way it would not be per pixel. + std::vector outward_max(max_radius + 2, 0.0f); + for (int rad = max_radius; rad >= 0; rad--) + outward_max[rad] = std::max(baseline[rad], outward_max[rad + 1]); + + int blocked_out_to = -1; + for (int rad = 0; rad <= max_radius; rad++) { + if (ring_pixels[rad] < MIN_RING_PIXELS) + continue; + if (baseline[rad] >= SHADOW_RATIO * outward_max[rad]) + break; + blocked_out_to = rad; + } - // A dip counts only where the background it is compared against was actually counted. std::vector low(n_pixels, 0); for (int i = 0; i < n_pixels; i++) { if (!valid[i]) continue; + if (radius[i] <= blocked_out_to) { + low[i] = 1; + continue; + } const double counted = frames * pooled_count[i]; - low[i] = (ratio[i] < SHADOW_RATIO && baseline[radius[i]] * counted >= MIN_EXPECTED_COUNTS) - || (ratio_radial[i] < SHADOW_RATIO && envelope[radius[i]] * counted >= MIN_EXPECTED_COUNTS); + low[i] = ratio[i] < SHADOW_RATIO && baseline[radius[i]] * counted >= MIN_EXPECTED_COUNTS; } // The shadow is the low region connected to the beam centre, bridging the gaps it crosses. @@ -365,7 +377,7 @@ std::vector ShadowFinder::GetMask() const { // Grow the soft boundary, round it and fill the disk interior. const std::vector penumbra = dilate(region, W, H, PENUMBRA_MAX_PX); for (int i = 0; i < n_pixels; i++) - if (penumbra[i] && valid[i] && std::min(ratio[i], ratio_radial[i]) < PENUMBRA_RATIO) + if (penumbra[i] && valid[i] && ratio[i] < PENUMBRA_RATIO) region[i] = 1; region = erode(dilate(region, W, H, 2), W, H, 2); diff --git a/image_analysis/beam_stop/ShadowFinder.h b/image_analysis/beam_stop/ShadowFinder.h index ca7aea34..ef56fa7d 100644 --- a/image_analysis/beam_stop/ShadowFinder.h +++ b/image_analysis/beam_stop/ShadowFinder.h @@ -18,10 +18,14 @@ // and is 1 where the beam stop shadows the detector. // // The shadow is a place where the background is missing, so it is found by comparing each -// pixel's mean against the typical background at the same radius. Two comparisons are -// needed: an azimuthal one (median over the ring) finds the arm, which is a minority of -// its ring, and a radial one (the background just outside) finds the disk, which is not - -// inside a fully blocked ring the ring median is itself the shadow. +// pixel's mean against the typical background at the same radius - the median over its ring, +// taken over the pixels not already known to be shadowed. That comparison holds wherever the +// ring still has unshadowed pixels to measure. Where it does not - a ring lying wholly inside +// the stop - there is nothing to compare against, and such a ring is shadow in its entirety. +// +// The background belongs to the beam and the shadow to the stop, and the two are not concentric: +// the stop sits off the beam by a sizeable fraction of its own radius. Only the per-ring +// comparison is used, so nothing here assumes they share a centre. // // Frames are chosen by the caller; the detection needs enough of them that the background // is counted rather than guessed (see MIN_EXPECTED_COUNTS in the .cpp). @@ -33,7 +37,6 @@ class ShadowFinder { const int height; const float beam_x; const float beam_y; - const int envelope_px; std::vector pixel_mask; // pixels already masked carry no background to test