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 860d0b6f1..a8a5a6bdc 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; @@ -851,10 +861,90 @@ 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); ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) { - for (int i = lo; i < hi; i++) - dim[i] = valid[i] && !region[i] && ratio[i] < PENUMBRA_RATIO && deficit[i] > MIN_DEFICIT_SIGMA; + 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); diff --git a/tests/ShadowFinderTest.cpp b/tests/ShadowFinderTest.cpp index 9d8617c5d..99473260e 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