From 7b55a31a3a3f0b6e111949cb6166efad7fa29ae4 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 27 Sep 2026 09:59:17 +0200 Subject: [PATCH] ShadowFinder: judge the arm search's dim pixels against the ring's azimuthal harmonic The transmitting-piece step added in a1816e905 called a pixel dim when it sat below 0.75 of its ring's median. When the polarization factor the background is divided by is not the beam's, every outer ring is left with a cos 2phi modulation: on one open-arm marCCD sweep (0.75 A, 90 mm, default p = 0.99) the corrected background runs 1.5-1.8x the ring median along x and 0.6-0.75x along y at r = 1600-2200 px. Joined, the y lobes made two pieces of 1.01 M and 0.39 M dim pixels holding only 0.17 % deep pixels (noise and a real ~20 % partial shadow) - past the absolute gate of 200 - and the mask grew from 16,584 to 1,428,328 pixels, covering the top and bottom of the detector. The step now fits, per 64 px radial band, the 24 sector medians with m + p cos 2phi + q sin 2phi (sectors below 0.75 of the fit left out, three times over) and judges dim against min(1, model). Capped at 1 the model only removes dim pixels, so a sweep with no piece keeps its mask bit for bit and a piece can only shrink or go. Shadow pixels (b90b377b8 -> this), same binary otherwise: that sweep 1,428,328 -> 29,821 (the stop's stick extension kept); arm sets 7atg 724,032 -> 723,708, 6oel 223,087 -> 216,126, 7rji 311,951 -> 309,816, 9fhc 173,597 -> 171,372, 5ebi 793,217 -> 785,477, 6zqr 289,324 -> 251,764, 5j23 530,541 -> 526,115, metformin 341,368 -> 231,123 (two speckled pieces in the high-angle lobe dropped, the arm below the beam kept). Merges: 6oel, 7rji, 6zqr, metformin unchanged in space group and resolution, R_meas/ISa within 0.1 %/0.2; 7atg fails either way on this base (the beam-centre defect fixed on atg-arm-centre). The sweep above: same space group, resolution, CC1/2 and ISa (6.72 -> 6.73); R_meas 21.7 -> 22.7 %, because the spurious mask also hid part of a real ~20 % partial shadow at the top right that neither version masks (it is above penumbra_ratio). New test ShadowFinder_AWrongPolarizationFactorIsNotAnArm (unpolarized background, 0.99 stated, small absorbing speck in the dim lobe) masks 4547 px before and 0 after; all *Shadow* tests pass. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- image_analysis/beam_stop/SHADOW_FINDER.md | 8 +- image_analysis/beam_stop/ShadowFinder.cpp | 96 ++++++++++++++++++++++- tests/ShadowFinderTest.cpp | 41 ++++++++++ 3 files changed, 142 insertions(+), 3 deletions(-) diff --git a/image_analysis/beam_stop/SHADOW_FINDER.md b/image_analysis/beam_stop/SHADOW_FINDER.md index 3e446ef10..f40c5f129 100644 --- a/image_analysis/beam_stop/SHADOW_FINDER.md +++ b/image_analysis/beam_stop/SHADOW_FINDER.md @@ -105,7 +105,13 @@ own radius**, not to a global threshold. recorded attenuated. The step only ever adds such pieces - the dim edge of an ordinary stop is left as steps 6-10 drew it, so a sweep without transmitting hardware keeps its mask bit for bit (an earlier version redrew every stop's edge from the dim pixels and moved every mask in the - corpus by 0.01-0.7 % of the detector). + corpus by 0.01-0.7 % of the detector). "Dim" here is against the ring **as it varies in + azimuth**: each 64 px radial band fits its sector medians (24 sectors, hardware sectors left out) + with a second harmonic, cos 2φ, and a pixel is dim below `penumbra_ratio` times that model where + the model is below 1. A polarization factor that is not the beam's leaves exactly such a harmonic + after the division, dim lobes at high angle over a quarter of the detector, which the median alone + read as one enormous transmitting piece. The model is capped at 1, so it only ever removes dim + pixels and a sweep without such a piece still keeps its mask bit for bit. Three findings are baked in: the **mean** ratio (not the max) does the detecting — the max projection is a great *picture* but a noisier detector; **size and counting statistics**, not a diff --git a/image_analysis/beam_stop/ShadowFinder.cpp b/image_analysis/beam_stop/ShadowFinder.cpp index 944e965e4..81845db38 100644 --- a/image_analysis/beam_stop/ShadowFinder.cpp +++ b/image_analysis/beam_stop/ShadowFinder.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include @@ -55,6 +56,15 @@ constexpr int MIN_SHADOW_PIXELS = 2000; // and deep nowhere. constexpr int MIN_CORE_PIXELS = 200; +// The arm search compares a pixel with its ring as the ring actually varies around the beam. A ring +// of background is not flat once divided by the polarization factor when that factor is not the +// beam's: the remainder is a second harmonic in azimuth, cos 2phi, which reaches tens of percent at +// high angle. It is measured over radial bands of HARMONIC_BAND_PX, from the median of each of +// HARMONIC_SECTORS sectors that holds at least MIN_SECTOR_PIXELS pixels. +constexpr int HARMONIC_BAND_PX = 64; +constexpr int HARMONIC_SECTORS = 24; +constexpr int MIN_SECTOR_PIXELS = 200; + // 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; @@ -846,9 +856,91 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { // It only adds whole pieces. The dim pixels along the edge of an ordinary stop are left as the // region above drew them: each such edge feeds decisions downstream that can sit at a margin, so // a sweep with no transmitting hardware keeps its mask bit for bit. + // + // Dim is judged against the ring's own azimuthal variation, not against its median alone. A + // polarization factor that is not the beam's leaves every ring brighter along one axis and dimmer + // along the other, smoothly and symmetrically about the beam; at high angle the dim lobes sit + // below PENUMBRA_RATIO over a quarter of the detector, deep in places from counting noise, and + // joined they made one piece that covered it. Hardware is not a second harmonic of the ring: it + // stands in a few sectors, and those sectors are left out of the fit. The harmonic is only ever + // allowed to explain a dim sector away - where it would ask for more than the median, the median + // stands - so the step can only drop what it found before, never find something new. + const int n_bands = max_radius / HARMONIC_BAND_PX + 1; + std::vector> sector_values(static_cast(n_bands) * HARMONIC_SECTORS); + for (int i = 0; i < n_pixels; i++) { + if (!valid[i] || region[i]) + continue; + const float dx = static_cast(i % W) - beam_x, dy = static_cast(i / W) - beam_y; + const double phi = std::atan2(dy, dx) + std::numbers::pi; + const int sector = std::min(HARMONIC_SECTORS - 1, static_cast(phi / (2.0 * std::numbers::pi) * HARMONIC_SECTORS)); + sector_values[static_cast(radius[i] / HARMONIC_BAND_PX) * HARMONIC_SECTORS + sector].push_back(ratio[i]); + } + std::vector harm_c(n_bands, 0.0f), harm_s(n_bands, 0.0f); + ParallelFor(n_bands, nthreads, [&](int band) { + std::vector med(HARMONIC_SECTORS, -1.0), c(HARMONIC_SECTORS), s(HARMONIC_SECTORS); + for (int k = 0; k < HARMONIC_SECTORS; k++) { + auto &v = sector_values[static_cast(band) * HARMONIC_SECTORS + k]; + if (v.size() >= MIN_SECTOR_PIXELS) { + std::nth_element(v.begin(), v.begin() + v.size() / 2, v.end()); + med[k] = v[v.size() / 2]; + } + const double phi = (k + 0.5) * 2.0 * std::numbers::pi / HARMONIC_SECTORS - std::numbers::pi; + c[k] = std::cos(2 * phi); + s[k] = std::sin(2 * phi); + } + // Least squares of med = m + p cos 2phi + q sin 2phi over the sectors that are not themselves + // dim against the fit, three times over as the baseline is; the model is relative, 1 + (p cos + q sin)/m. + std::vector use(HARMONIC_SECTORS); + for (int k = 0; k < HARMONIC_SECTORS; k++) + use[k] = med[k] >= 0; + for (int iter = 0; iter < 3; iter++) { + double n = 0, sc = 0, ss = 0, scc = 0, sss = 0, scs = 0, y = 0, yc = 0, ys = 0; + for (int k = 0; k < HARMONIC_SECTORS; k++) { + if (!use[k]) continue; + n++; sc += c[k]; ss += s[k]; scc += c[k] * c[k]; sss += s[k] * s[k]; scs += c[k] * s[k]; + y += med[k]; yc += med[k] * c[k]; ys += med[k] * s[k]; + } + // Sectors crowded into a narrow arc cannot tell a harmonic from a level; the determinant of + // the normal matrix, per sector cubed, is 1/4 on a full ring. + const double det = n * (scc * sss - scs * scs) - sc * (sc * sss - scs * ss) + ss * (sc * scs - scc * ss); + if (n < 6 || det < 0.01 * n * n * n) { + harm_c[band] = harm_s[band] = 0.0f; + return; + } + const double m = (y * (scc * sss - scs * scs) - sc * (yc * sss - scs * ys) + ss * (yc * scs - scc * ys)) / det; + const double p = (n * (yc * sss - ys * scs) - y * (sc * sss - scs * ss) + ss * (sc * ys - yc * ss)) / det; + const double q = (n * (scc * ys - scs * yc) - sc * (sc * ys - yc * ss) + y * (sc * scs - scc * ss)) / det; + if (m <= 0) { + harm_c[band] = harm_s[band] = 0.0f; + return; + } + harm_c[band] = static_cast(p / m); + harm_s[band] = static_cast(q / m); + for (int k = 0; k < HARMONIC_SECTORS; k++) + use[k] = med[k] >= 0 && med[k] >= PENUMBRA_RATIO * (1.0 + harm_c[band] * c[k] + harm_s[band] * s[k]); + } + }); + std::vector dim(n_pixels, 0); - for (int i = 0; i < n_pixels; i++) - dim[i] = valid[i] && !region[i] && ratio[i] < PENUMBRA_RATIO && deficit[i] > MIN_DEFICIT_SIGMA; + ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) { + for (int i = lo; i < hi; i++) { + if (!valid[i] || region[i] || ratio[i] >= PENUMBRA_RATIO) + continue; + // cos 2phi and sin 2phi from the offset to the beam. + const float dx = static_cast(i % W) - beam_x, dy = static_cast(i / W) - beam_y; + const float r2 = std::max(dx * dx + dy * dy, 1e-6f); + const int band = radius[i] / HARMONIC_BAND_PX; + const float model = std::min(1.0f, 1.0f + harm_c[band] * (dx * dx - dy * dy) / r2 + + harm_s[band] * 2.0f * dx * dy / r2); + if (model >= 1.0f) { + dim[i] = deficit[i] > MIN_DEFICIT_SIGMA; + continue; + } + const double counted = static_cast(frames) * pooled_count[i] * pol[i]; + dim[i] = ratio[i] < PENUMBRA_RATIO * model + && poisson_deficit_sigma(pooled[i] * counted, baseline[radius[i]] * model * counted) > MIN_DEFICIT_SIGMA; + } + }); const std::vector joined = bridge_gaps(dilate(dim, W, H, BRIDGE_PX, nthreads), valid, W, H); std::vector seen(n_pixels, 0); std::vector component; diff --git a/tests/ShadowFinderTest.cpp b/tests/ShadowFinderTest.cpp index 0cae4a368..eb7596719 100644 --- a/tests/ShadowFinderTest.cpp +++ b/tests/ShadowFinderTest.cpp @@ -383,6 +383,47 @@ TEST_CASE("ShadowFinder_APolarizedBackgroundIsNotAShadow", "[ShadowFinder]") { CHECK(mask_for(std::nullopt) > static_cast(W) * H / 20); } +// The other side of that correction: a polarization factor that is not the beam's. Here the +// background is unpolarized and the experiment states 0.99, so after the division every outer ring is +// bright along x and dim along y - a smooth lobe below PENUMBRA_RATIO across the top and bottom of the +// detector with nothing in the beam. On its own the lobe is never deep, but a small absorbing speck +// inside it is, and the transmitting-hardware step joined the two into one piece "dim everywhere and +// deep in places" and masked the whole lobe. The step compares a pixel with its ring's own second +// harmonic, so the lobe is explained away and the speck, far below a shadow's size, stays unmasked. +TEST_CASE("ShadowFinder_AWrongPolarizationFactorIsNotAnArm", "[ShadowFinder]") { + constexpr float DISTANCE_MM = 6.0f; + constexpr int LEVEL = 24; + constexpr int SPECK_X = C, SPECK_Y = 20, SPECK_R = 11; + + DiffractionExperiment x = TestExperiment(); + x.DetectorDistance_mm(DISTANCE_MM).PolarizationFactor(0.99f); + const float pixel_mm = x.GetPixelSize_mm(); + + std::vector frame(static_cast(W) * H); + for (int y = 0; y < H; y++) + for (int xi = 0; xi < W; xi++) { + const double rho = std::hypot(xi - C, y - C) * pixel_mm; + const double cos_2theta_sq = DISTANCE_MM * DISTANCE_MM / (DISTANCE_MM * DISTANCE_MM + rho * rho); + frame[I(xi, y)] = static_cast(std::lround(LEVEL * 0.5 * (1.0 + cos_2theta_sq))); + const int sx = xi - SPECK_X, sy = y - SPECK_Y; + if (sx * sx + sy * sy <= SPECK_R * SPECK_R) + frame[I(xi, y)] /= 10; + } + + const PixelMask pixel_mask(x); + ShadowFinder finder(x, pixel_mask); + std::vector buffer; + for (int f = 0; f < NFRAMES; f++) { + DataMessage msg{}; + msg.image = CompressedImage(frame, W, H); + finder.AddImage(msg, buffer); + } + const auto mask = finder.GetMask(); + CHECK(mask[I(C, 5)] == 0); // the dim lobe along y + CHECK(mask[I(SPECK_X, SPECK_Y)] == 0); // the speck, too small to be a shadow + CHECK(std::count(mask.begin(), mask.end(), 1u) == 0); +} + // Inside the stop a whole ring is blocked, so its own median is blocked too and the per-pixel // comparison has nothing to work with; the walk that finds those rings is what covers the disk. // It decides by comparing a ring against what this detector's background typically is - and NOT