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