From 59ae6fc5c8bd80d8686b3a62fc068c9d3b624d9b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 23 Oct 2026 02:47:42 +0200 Subject: [PATCH] ShadowFinder: measure the background's azimuthal modulation instead of predicting it from the polarization The shadow finder divided the background by the Kahn factor computed from the experiment's polarization degree, plane and detector distance before comparing each pixel with its ring. On a stills run none of the three is trustworthy at the pre-scan: the degree is not identifiable from the data, the plane can lie along the image's y axis rather than lab x (a detector mounted turned), and the distance is still the file's placeholder, so 2 theta is overestimated. On such a run the division made the ring anisotropy worse instead of removing it: a residual cos 2phi of 0.8 at high angle, dim lobes over the detector corners sitting on the PENUMBRA_RATIO threshold, and the transmitting-hardware search turned them into false pieces of 17-21k pixels whose extent followed the third decimal of p (22973 / 25102 / 26190 pixels masked at p = 0.99 / 0.995 / 0.998, of which only 5248 were the actual stop). A polarized background is exactly a second harmonic in azimuth whatever the degree, plane and 2 theta, so it is now measured: the pooled background relative to its ring median, the medians of 24 sectors per 64 px band fitted with 1 + c cos 2phi + s sin 2phi (HarmonicFit, dim sectors left out), band values interpolated in radius, a band that cannot be fitted taking the nearest fitted one. The background is divided by that before the ring baseline is taken; the Poisson test still puts the factor back. The finder no longer reads the polarization at all, on the host and on the device. The arm search's own capped harmonic, which patched the old residual, is gone: the ratio it works on is flat now. On the stills set that showed it, the mask is the 5248-pixel stop at every p and at either plane, and the merge statistics are the same at p = 0.99, 0.995 and 0.998 except for continuous changes from the integration's own polarization correction. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu --- image_analysis/beam_stop/SHADOW_FINDER.md | 60 ++-- image_analysis/beam_stop/ShadowFinder.cpp | 294 ++++++++++-------- image_analysis/beam_stop/ShadowFinder.h | 16 +- .../beam_stop/ShadowFinderInternal.h | 34 +- image_analysis/beam_stop/ShadowMaskGPU.cu | 144 +++++---- image_analysis/beam_stop/ShadowMaskGPU.h | 11 +- tests/ShadowFinderTest.cpp | 90 ++++-- 7 files changed, 369 insertions(+), 280 deletions(-) diff --git a/image_analysis/beam_stop/SHADOW_FINDER.md b/image_analysis/beam_stop/SHADOW_FINDER.md index a85a4d26f..a40b000d4 100644 --- a/image_analysis/beam_stop/SHADOW_FINDER.md +++ b/image_analysis/beam_stop/SHADOW_FINDER.md @@ -44,20 +44,27 @@ own radius**, not to a global threshold. 1. **Projection.** Per pixel accumulate the **max** and the **mean** (sum + count) over the added frames, skipping the no-data sentinel. The stop is stationary, so it stays dark while spots and the air-scatter halo build up around it. -2. **Polarization.** The background is not flat around a ring even with nothing in the beam: a - polarized source suppresses it in its own plane by `1 - sin^2(2 theta)`, which is a factor of - three at `2 theta = 55 deg` and four at 70 deg — several times the dip being looked for, so on a - short-distance geometry the two horizontal lobes of every outer ring read as shadow. The mean is - therefore divided by the Kahn factor before the comparison, read off the experiment's own - geometry so that detector tilt, a quarter-turned image and an in-plane rotation all move the - polarization plane correctly, about the centre named in step 3. On a **tilted** detector the - solid angle of a pixel also varies along a ring (by tens of per cent at a 5 deg 2theta swing and - 50 mm), so the mean is divided by that too: the pixel's solid angle over the one an untilted - detector has at the same `2 theta`, `(distance / l_z)^3`, which is exactly 1 untilted. Detector and - air absorption are functions of `2 theta` alone, nearly enough, and the ring's own median absorbs - them. The Poisson test of step 5 puts the factors back, so the significance is still the - significance of the counts that were recorded. No polarization is corrected when the experiment - states none. +2. **Azimuthal modulation.** The background is not flat around a ring even with nothing in the + beam: a polarized source suppresses it in its own plane by `1 - sin^2(2 theta)`, which is a factor + of three at `2 theta = 55 deg` and four at 70 deg — several times the dip being looked for, so on a + short-distance geometry the two lobes of every outer ring read as shadow. Whatever the degree of + polarization, the plane and the `2 theta`, that modulation is exactly a second harmonic in + azimuth, `cos 2(phi - phi0)`, so it is **measured, not predicted**: the pooled mean is taken + relative to its ring median (step 3), each 64 px radial band fits its 24 sector medians with + `1 + c cos 2phi + s sin 2phi` (sectors dim against the fit, where hardware stands, left out), the + band harmonics are interpolated linearly in radius (a band that cannot be fitted - inside the stop, + or an arc at a detector's corner - takes the nearest fitted band's), and the background is divided + by the result. The finder reads no polarization from the experiment: a prediction needs the degree + of polarization, which the data cannot confirm, the plane, which differs between experiments on the + same detector, and the detector distance, which a stills pre-scan does not know yet - and where any + of them was wrong, the division made a ring less flat than it found it, and a change of p in the + third decimal moved thousands of pixels in and out of the mask. On a **tilted** detector the solid + angle of a pixel also varies along a ring (by tens of per cent at a 5 deg 2theta swing and 50 mm), so + the mean is divided by that too: the pixel's solid angle over the one an untilted detector has at the + same `2 theta`, `(distance / l_z)^3`, which is exactly 1 untilted. Detector and air absorption are + functions of `2 theta` alone, nearly enough, and the ring's own median absorbs them. The Poisson test + of step 5 puts both factors back, so the significance is still the significance of the counts that + were recorded. 3. **Robust radial baseline.** Take the **median of the mean-projection in each 1-pixel ring** about the beam. A ring is a cone of constant `2 theta`, indexed by `distance * tan(2 theta)` in pixels: a circle about the centre on an untilted detector, but on a tilted one (a swung 2theta @@ -109,15 +116,10 @@ own radius**, not to a global threshold. piece with at least `min_shadow_pixels` dim pixels of which `min_core_pixels` are core is added whole. Nothing is given back inside it: the reflections behind transmitting hardware are 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). "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. + left as steps 6-10 drew it, so a sweep without transmitting hardware keeps its mask bit for bit. + "Dim" is against the ring with its azimuthal modulation divided out (step 2): left in, the dim + lobes of a polarized background, at high angle over a quarter of the detector, joined into one + enormous transmitting piece. 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 @@ -197,13 +199,11 @@ beam-stop mask is an **offline** (Rugnux + viewer) feature and never goes to the fit moves it by a few pixels, three orders below what the ring comparison notices. What remains is the case where the fit declines to place a centre at all, where the file's is used as before. Another caller that reads the mask at the file's centre keeps the old exposure. -- **The polarization plane is taken from the experiment, not measured.** The correction follows - `DiffractionGeometry::CalcAzIntPolarizationCorr`, which fixes the plane on the lab x axis and - expresses a vertically polarized source as a negative coefficient — the same assumption azimuthal - integration already makes, so a beamline where it is wrong is wrong there first. Measured over the - corpus the plane sits within a degree of lab x on every dataset, and the amplitude within tens of - per cent of the prediction, which is far inside the margin the 0.50 threshold leaves. An - experiment that states no polarization is left uncorrected, as before. +- **The azimuthal modulation is measured per 64 px band.** A real polarization factor is a second + harmonic of exactly this kind, but its amplitude changes within a band; the linear interpolation + between band centres follows it to a few per cent at high angle, far inside the margin the 0.50 and + 0.75 thresholds leave. Anything else that is a second harmonic about the beam is divided out with + it, including part of the fall-off seen about a badly wrong centre. - **Unsigned saturation**: for unsigned images a detector-saturating reflection collides with the no-data sentinel and is skipped; harmless in practice (such reflections do not sit under the stop), and absent entirely on the signed/preprocessed path. diff --git a/image_analysis/beam_stop/ShadowFinder.cpp b/image_analysis/beam_stop/ShadowFinder.cpp index 3fa23b110..dfc85a220 100644 --- a/image_analysis/beam_stop/ShadowFinder.cpp +++ b/image_analysis/beam_stop/ShadowFinder.cpp @@ -36,16 +36,6 @@ namespace { template using Plane = std::vector>; -// A Plane that starts at `value` everywhere, for the ones that are filled piecemeal afterwards. -template -Plane filled_plane(int n, T value, size_t nthreads) { - Plane out(n); - ParallelChunks(n, nthreads, [&](int lo, int hi) { - std::fill(out.begin() + lo, out.begin() + hi, value); - }); - return out; -} - // Deficit of an observed count against its expectation, in standard deviations, from the Poisson // likelihood ratio. Zero when the observation is not below the expectation. double poisson_deficit_sigma(double observed, double expected) { @@ -399,6 +389,78 @@ RingValues bin_by_ring(const Plane &values, const Plane &valid, return rv; } +// The median of each ring, iterated so the shadow stays out of it: three times, each excluding the +// ring's values below SHADOW_RATIO of the median before. +// +// The iteration only ever excludes values below a cut, and dividing by a positive median is monotone - +// so the excluded values of a ring are exactly the lowest ones, and the values the next median is taken +// over are exactly the rest. Each iteration's median is therefore an order statistic of the ring's +// sorted values: each iteration picks a rank and counts a prefix. +std::vector ring_baseline(const RingValues &rings, int max_radius, size_t nthreads) { + std::vector baseline(max_radius + 1, 0.0f); + std::vector excluded_in_ring(max_radius + 1, 0); + for (int iter = 0; iter < 3; iter++) + ParallelFor(max_radius + 1, nthreads, [&](int r) { + const int lo = rings.offset[r], hi = rings.offset[r + 1]; + const int n = hi - lo; + const int m = excluded_in_ring[r]; + const int avail = n - m; + baseline[r] = (avail <= 0) ? 0.0f : rings.values[lo + m + avail / 2]; + // Counted the same way the per-pixel test is written, so the two agree bit for bit; the + // ring is sorted, so this is the length of a prefix. + const float d = std::max(baseline[r], 1e-6f); + int excl = 0; + while (excl < n && rings.values[lo + excl] / d < SHADOW_RATIO) + excl++; + excluded_in_ring[r] = excl; + }); + return baseline; +} + +// The median of `values` in each sector of each radial band (band * HARMONIC_SECTORS + sector), over the +// pixels `use` marks; -1 for a sector with fewer than MIN_SECTOR_PIXELS of them. The sector is the +// pixel's azimuth in the lab frame (ring_position). +std::vector sector_medians(const Plane &values, const Plane &use, const Plane &radius, + int W, float centre_x, float centre_y, const float *det_matrix, + float distance_px, int n_bands, size_t nthreads) { + const size_t n_sectors = static_cast(n_bands) * HARMONIC_SECTORS; + const int n_pixels = static_cast(values.size()); + // Gathered by blocks of rows in parallel and joined in block order. Only the median of each + // sector is read, and that is the same whatever order its values were gathered in. + constexpr int SECTOR_BLOCKS = 64; + std::vector>> block_values(SECTOR_BLOCKS); + ParallelFor(SECTOR_BLOCKS, nthreads, [&](int b) { + auto &v = block_values[b]; + v.resize(n_sectors); + const int lo = static_cast(static_cast(n_pixels) * b / SECTOR_BLOCKS); + const int hi = static_cast(static_cast(n_pixels) * (b + 1) / SECTOR_BLOCKS); + for (int i = lo; i < hi; i++) { + if (!use[i]) + continue; + const RingPosition ring = ring_position(det_matrix, distance_px, + static_cast(i % W) - centre_x, + static_cast(i / W) - centre_y); + const double phi = std::atan2(ring.ly, ring.lx) + std::numbers::pi; + const int sector = std::min(HARMONIC_SECTORS - 1, static_cast(phi / (2.0 * std::numbers::pi) * HARMONIC_SECTORS)); + v[static_cast(radius[i] / HARMONIC_BAND_PX) * HARMONIC_SECTORS + sector].push_back(values[i]); + } + }); + std::vector> sector_values(n_sectors); + ParallelFor(static_cast(n_sectors), nthreads, [&](int k) { + for (const auto &v : block_values) + sector_values[k].insert(sector_values[k].end(), v[k].begin(), v[k].end()); + }); + std::vector median(n_sectors, -1.0); + ParallelFor(static_cast(n_sectors), nthreads, [&](int k) { + auto &v = sector_values[k]; + if (v.size() >= MIN_SECTOR_PIXELS) { + std::nth_element(v.begin(), v.begin() + v.size() / 2, v.end()); + median[k] = v[v.size() / 2]; + } + }); + return median; +} + } // namespace ShadowFinder::ShadowFinder(const DiffractionExperiment &experiment, const PixelMask &mask) @@ -407,7 +469,6 @@ ShadowFinder::ShadowFinder(const DiffractionExperiment &experiment, const PixelM beam_x(experiment.GetBeamX_pxl()), beam_y(experiment.GetBeamY_pxl()), geometry(experiment.GetDiffractionGeometry()), - polarization(experiment.GetPolarizationFactor()), pixel_mask(mask.GetMask(experiment)), modules(experiment.GetModules()) { if (pixel_mask.size() != static_cast(width) * height) @@ -669,8 +730,6 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { setup.det_matrix[k] = rot[k]; setup.pixel_size_mm = geometry.GetPixelSize_mm(); setup.distance_mm = geometry.GetDetectorDistance_mm(); - setup.has_polarization = polarization.has_value(); - setup.polarization = polarization.value_or(0.0f); ScopedGpuPin pin(gpu->GetDevice()); return gpu->Mask(setup, pixel_mask); } @@ -718,19 +777,16 @@ std::vector ShadowFinder::MaskOn(const Projection &p, const std::vecto return mask; // mean projection, usable pixels and the ring each pixel is on (ring_position). The mean is - // divided by the polarization factor and by the obliquity of a tilted detector, so that what is - // left varies around a ring only where something is in the way; `pol` keeps the product because - // the counts the Poisson test is made of are the ones that were recorded, not these. Both are - // read off the experiment's own geometry - which follows the detector's tilt, quarter turns and - // in-plane rotation - about the centre the caller named. - // Following Kahn, Fourme, Gadet, Janin, Dumas & Andre (1982) J. Appl. Cryst. 15, 330-337 - DiffractionGeometry pol_geometry = geometry; - pol_geometry.BeamX_pxl(centre_x).BeamY_pxl(centre_y); + // divided by the obliquity of a tilted detector here and by the background's azimuthal modulation + // below, so that what is left varies around a ring only where something is in the way; `shape` + // keeps the product because the counts the Poisson test is made of are the ones that were + // recorded, not these. The ring follows the experiment's own geometry - the detector's tilt, + // quarter turns and in-plane rotation - about the centre the caller named. const std::vector det_matrix = geometry.GetDetectorMatrix().arr(); const float distance_px = geometry.GetDetectorDistance_mm() / geometry.GetPixelSize_mm(); Plane mean(n_pixels); - Plane pol(n_pixels); + Plane shape(n_pixels); Plane valid(n_pixels); Plane radius(n_pixels); std::atomic max_radius_atomic{0}; @@ -740,21 +796,14 @@ std::vector ShadowFinder::MaskOn(const Projection &p, const std::vecto for (int x = 0; x < W; x++) { const int i = y * W + x; const RingPosition ring = ring_position(det_matrix.data(), distance_px, x - centre_x, y - centre_y); - pol[i] = ring.obliquity - * (polarization.has_value() - ? pol_geometry.CalcAzIntPolarizationCorr(static_cast(x), - static_cast(y), - polarization.value()) - : 1.0f); - // pol > 0 also rejects a NaN factor. The correction is exactly zero only where - // the polarization fully suppresses the scattering, so such a pixel carries no - // usable ring information anyway - and divided by zero it would put an inf into - // the pooled means, which the running box sums then turn into NaN for a whole row. + shape[i] = ring.obliquity; + // shape > 0 also rejects a NaN obliquity, which divided into the pooled means would + // turn a whole row of the running box sums into NaN. mean[i] = 0.0f; valid[i] = 0; - if (valid_count[i] > 0 && known_mask[i] == 0 && pol[i] > 0.0f) { + if (valid_count[i] > 0 && known_mask[i] == 0 && shape[i] > 0.0f) { mean[i] = static_cast(static_cast(sum_value[i]) - / valid_count[i] / pol[i]); + / valid_count[i] / shape[i]); valid[i] = 1; } radius[i] = static_cast(std::lround(ring.radius)); @@ -766,6 +815,7 @@ std::vector ShadowFinder::MaskOn(const Projection &p, const std::vecto while (prev < local_max && !max_radius_atomic.compare_exchange_weak(prev, local_max)) {} }); const int max_radius = max_radius_atomic.load(); + const int n_bands = max_radius / HARMONIC_BAND_PX + 1; // Pool the background over a small box before testing it. A background of a fraction of @@ -793,37 +843,57 @@ std::vector ShadowFinder::MaskOn(const Projection &p, const std::vecto }); + // The background is not flat around a ring even with nothing in the beam: a polarized source + // suppresses it in its own plane, by a factor of three at 2 theta = 55 deg and four at 70 deg - + // several times the dip being looked for, so on a short-distance geometry the two lobes of every + // outer ring read as shadow. Whatever the degree of polarization, the plane it lies in and the + // 2 theta, that modulation is a second harmonic in azimuth, cos 2(phi - phi0), so it is measured + // here rather than predicted: each radial band's harmonic is fitted to the medians of its sectors, + // taken relative to their rings, and the background is divided by it. A prediction needs numbers + // the data does not confirm - the degree of polarization, the plane, and on a stills run a + // detector distance that is not known yet when the pre-scan runs - and where they were wrong the + // division left a ring less flat than it found it, with large dim lobes sitting on the thresholds + // below, so that a change of p in the third decimal moved thousands of pixels in and out of the + // mask. Hardware is not a second harmonic of the ring: it stands in a few sectors, which the fit + // leaves out. + { + const auto level = ring_baseline(bin_by_ring(pooled, valid, radius, max_radius, nthreads), max_radius, + nthreads); + Plane relative(n_pixels); + ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) { + for (int i = lo; i < hi; i++) + relative[i] = valid[i] ? pooled[i] / std::max(level[radius[i]], 1e-6f) : 0.0f; + }); + std::vector shape_c, shape_s; + AzimuthalShapeFit(sector_medians(relative, valid, radius, W, centre_x, centre_y, det_matrix.data(), + distance_px, n_bands, nthreads), + n_bands, shape_c, shape_s); + ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) { + for (int i = lo; i < hi; i++) { + if (!valid[i]) + continue; + const RingPosition ring = ring_position(det_matrix.data(), distance_px, + static_cast(i % W) - centre_x, + static_cast(i / W) - centre_y); + const float a = azimuthal_shape(shape_c.data(), shape_s.data(), n_bands, radius[i], ring.lx, ring.ly); + // Zero only where the source suppresses the scattering completely, so such a pixel + // carries no usable ring information. + if (a > 0.0f) { + pooled[i] /= a; + shape[i] *= a; + } else + valid[i] = 0; + } + }); + } + + // Azimuthal comparison: the median of the ring, iterated so the shadow stays out of the // baseline it is measured against. Comparing a pixel's background against the background at its // own resolution is how DEFPIX recognises a shaded detector region. // Following Kabsch (2010) Acta Cryst. D66, 125-132 - // - // The iteration only ever excludes pixels whose pooled background is below a cut, and dividing by - // a positive baseline is monotone - so the excluded pixels of a ring are exactly the lowest ones, - // and the pixels the next median is taken over are exactly the rest. Each iteration's median is - // therefore an order statistic of the ring's values, which do not change: bin and sort the rings - // once, then each iteration picks a rank and counts a prefix. That replaces nine full-image - // passes (use / ratio / excluded, three times) with one, and three re-binnings with none. const RingValues rings = bin_by_ring(pooled, valid, radius, max_radius, nthreads); - std::vector baseline(max_radius + 1, 0.0f); - { - std::vector excluded_in_ring(max_radius + 1, 0); - for (int iter = 0; iter < 3; iter++) - ParallelFor(max_radius + 1, nthreads, [&](int r) { - const int lo = rings.offset[r], hi = rings.offset[r + 1]; - const int n = hi - lo; - const int m = excluded_in_ring[r]; - const int avail = n - m; - baseline[r] = (avail <= 0) ? 0.0f : rings.values[lo + m + avail / 2]; - // Counted the same way the per-pixel test below is written, so the two agree bit for - // bit; the ring is sorted, so this is the length of a prefix. - const float d = std::max(baseline[r], 1e-6f); - int excl = 0; - while (excl < n && rings.values[lo + excl] / d < SHADOW_RATIO) - excl++; - excluded_in_ring[r] = excl; - }); - } + const std::vector baseline = ring_baseline(rings, max_radius, nthreads); Plane ratio(n_pixels); ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) { for (int i = lo; i < hi; i++) @@ -867,7 +937,7 @@ std::vector ShadowFinder::MaskOn(const Projection &p, const std::vecto // background and what the ring expects of it are both put back on the detector's own // scale before they are compared - the same factor on each, so it is the units of the // comparison that change and not its answer. - const double counted = static_cast(frames) * pooled_count[i] * pol[i]; + const double counted = static_cast(frames) * pooled_count[i] * shape[i]; deficit[i] = static_cast(poisson_deficit_sigma(pooled[i] * counted, baseline[radius[i]] * counted)); low[i] = ratio[i] < SHADOW_RATIO && deficit[i] > MIN_DEFICIT_SIGMA; @@ -970,76 +1040,13 @@ std::vector ShadowFinder::MaskOn(const Projection &p, const std::vecto // 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; - const size_t n_sectors = static_cast(n_bands) * HARMONIC_SECTORS; - // Gathered by blocks of rows in parallel and joined in block order. Only the median of each - // sector is read, and that is the same whatever order its values were gathered in. - constexpr int SECTOR_BLOCKS = 64; - std::vector>> block_values(SECTOR_BLOCKS); - ParallelFor(SECTOR_BLOCKS, nthreads, [&](int b) { - auto &values = block_values[b]; - values.resize(n_sectors); - const int lo = static_cast(static_cast(n_pixels) * b / SECTOR_BLOCKS); - const int hi = static_cast(static_cast(n_pixels) * (b + 1) / SECTOR_BLOCKS); - for (int i = lo; i < hi; i++) { - if (!valid[i] || region[i]) - continue; - const RingPosition ring = ring_position(det_matrix.data(), distance_px, - static_cast(i % W) - centre_x, - static_cast(i / W) - centre_y); - const double phi = std::atan2(ring.ly, ring.lx) + std::numbers::pi; - const int sector = std::min(HARMONIC_SECTORS - 1, static_cast(phi / (2.0 * std::numbers::pi) * HARMONIC_SECTORS)); - values[static_cast(radius[i] / HARMONIC_BAND_PX) * HARMONIC_SECTORS + sector].push_back(ratio[i]); - } - }); - std::vector> sector_values(n_sectors); - ParallelFor(static_cast(n_sectors), nthreads, [&](int k) { - for (const auto &values : block_values) - sector_values[k].insert(sector_values[k].end(), values[k].begin(), values[k].end()); - }); - block_values.clear(); - // The median of each sector with enough pixels to have one. - std::vector sector_median(n_sectors, -1.0); - ParallelFor(static_cast(n_sectors), nthreads, [&](int k) { - auto &v = sector_values[k]; - if (v.size() >= MIN_SECTOR_PIXELS) { - std::nth_element(v.begin(), v.begin() + v.size() / 2, v.end()); - sector_median[k] = v[v.size() / 2]; - } - }); - std::vector harm_c, harm_s; - HarmonicFit(sector_median, n_bands, harm_c, harm_s); - - Plane dim = filled_plane(n_pixels, 0, nthreads); + // Dim is judged against the ring with its azimuthal modulation divided out (above): left in, the + // dim lobes of a polarized background sit below PENUMBRA_RATIO over a quarter of the detector at + // high angle, deep in places from counting noise, and joined they make one piece that covers it. + Plane dim(n_pixels); 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 RingPosition ring = ring_position(det_matrix.data(), distance_px, - static_cast(i % W) - centre_x, - static_cast(i / W) - centre_y); - const float dx = ring.lx, dy = ring.ly; - 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; - } + for (int i = lo; i < hi; i++) + dim[i] = valid[i] && !region[i] && ratio[i] < PENUMBRA_RATIO && deficit[i] > MIN_DEFICIT_SIGMA; }); const Plane joined = bridge_gaps(dilate(dim, W, H, BRIDGE_PX, nthreads), valid, W, H, nthreads); const auto pieces = label_components(joined, W, H, nthreads); @@ -1100,9 +1107,10 @@ int BlockedOutTo(const std::vector &baseline, const std::vector &rin } void HarmonicFit(const std::vector §or_median, int n_bands, - std::vector &harm_c, std::vector &harm_s) { + std::vector &harm_c, std::vector &harm_s, std::vector &fitted) { harm_c.assign(n_bands, 0.0f); harm_s.assign(n_bands, 0.0f); + fitted.assign(n_bands, 0); for (int band = 0; band < n_bands; band++) { std::vector med(HARMONIC_SECTORS), c(HARMONIC_SECTORS), s(HARMONIC_SECTORS); for (int k = 0; k < HARMONIC_SECTORS; k++) { @@ -1128,6 +1136,7 @@ void HarmonicFit(const std::vector §or_median, int n_bands, 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; + fitted[band] = 0; break; } const double m = (y * (scc * sss - scs * scs) - sc * (yc * sss - scs * ys) + ss * (yc * scs - scc * ys)) / det; @@ -1135,14 +1144,37 @@ void HarmonicFit(const std::vector §or_median, int n_bands, 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; + fitted[band] = 0; break; } harm_c[band] = static_cast(p / m); harm_s[band] = static_cast(q / m); + fitted[band] = 1; 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]); } } } +void AzimuthalShapeFit(const std::vector §or_median, int n_bands, + std::vector &harm_c, std::vector &harm_s) { + std::vector c, s; + std::vector fitted; + HarmonicFit(sector_median, n_bands, c, s, fitted); + harm_c.assign(n_bands, 0.0f); + harm_s.assign(n_bands, 0.0f); + for (int band = 0; band < n_bands; band++) { + // The nearest fitted band, the inner one of two equally near. + for (int d = 0; d < n_bands; d++) { + const int k = band - d >= 0 && fitted[band - d] ? band - d + : band + d < n_bands && fitted[band + d] ? band + d : -1; + if (k >= 0) { + harm_c[band] = c[k]; + harm_s[band] = s[k]; + break; + } + } + } +} + } // namespace shadow_finder diff --git a/image_analysis/beam_stop/ShadowFinder.h b/image_analysis/beam_stop/ShadowFinder.h index 1c844f5cd..35515844f 100644 --- a/image_analysis/beam_stop/ShadowFinder.h +++ b/image_analysis/beam_stop/ShadowFinder.h @@ -48,19 +48,11 @@ class ShadowFinder { float beam_x; float beam_y; - // What the background owes to the source rather than to the hardware. The scattered background - // is not flat around a ring: a polarized source suppresses it in its own plane by a factor that - // reaches three at the 2 theta a short detector distance puts in a corner - several times the - // dip this class is looking for - so the comparison divides it out before it compares. The only - // other correction that varies along a ring is the solid angle of a tilted detector, divided out - // with it (ring_position); detector and air absorption are functions of 2 theta alone, nearly - // enough, and the ring's own median absorbs them. - // The geometry is kept whole rather than reduced to a distance and a pixel size because the - // azimuth is the whole point: detector tilt, a quarter-turned image and an in-plane rotation - // all move the polarization plane across the stored image, and the geometry already knows where - // it lies. Only the centre is replaced, by the one the caller measured. + // The geometry the rings are drawn in. A ring is a cone of constant 2 theta, so it is kept whole + // rather than reduced to a distance and a pixel size: a detector tilt makes a ring a circle about + // no pixel at all, and the azimuth the background's modulation is measured in is the lab's. Only + // the centre is replaced, by the one the caller measured. const DiffractionGeometry geometry; - const std::optional polarization; // unset leaves the background as it was measured std::vector pixel_mask; // pixels already masked carry no background to test // The panels the image is made of. With any, the mask is made on the detector plane (GetMask). diff --git a/image_analysis/beam_stop/ShadowFinderInternal.h b/image_analysis/beam_stop/ShadowFinderInternal.h index 6b1f0d55d..c3c2900d4 100644 --- a/image_analysis/beam_stop/ShadowFinderInternal.h +++ b/image_analysis/beam_stop/ShadowFinderInternal.h @@ -28,7 +28,7 @@ namespace shadow_finder { // its offset from the beam in that frame, which gives its azimuth. A tilt also makes a pixel's solid // angle vary along the ring, by tens of per cent at that same tilt; `obliquity` is the pixel's solid // angle over the one the untilted detector has at the same 2theta, (distance / lz)^3, and the -// background is divided by it as by the polarization factor. +// background is divided by it as by its azimuthal modulation. // (dx, dy) is the pixel's offset from the geometry's beam centre, `det_matrix` the detector matrix // (row major) and `distance_px` the detector distance in pixels. Untilted, (lx, ly) is (dx, dy) and // the obliquity is 1, exactly. @@ -87,10 +87,9 @@ inline constexpr int MIN_SHADOW_PIXELS = 2000; // and deep nowhere. inline 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 +// A ring of background is not flat around the beam: a polarized source suppresses it in its own plane, +// and the modulation 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. inline constexpr int HARMONIC_BAND_PX = 64; inline constexpr int HARMONIC_SECTORS = 24; @@ -116,8 +115,29 @@ int BlockedOutTo(const std::vector &baseline, const std::vector &rin // The second harmonic in azimuth of each radial band, relative to its level, fitted to the medians of // its HARMONIC_SECTORS sectors (sector_median[band * HARMONIC_SECTORS + k], negative where the sector -// has too few pixels). See GetMask. +// has too few pixels). A band whose sectors cannot be fitted gets 0 and fitted[band] = 0. See GetMask. void HarmonicFit(const std::vector §or_median, int n_bands, - std::vector &harm_c, std::vector &harm_s); + std::vector &harm_c, std::vector &harm_s, std::vector &fitted); + +// The background's own azimuthal modulation, band by band: HarmonicFit, with a band that could not be +// fitted - inside the stop, or at a detector's corners where a band is only an arc - taking the nearest +// fitted band's harmonic. 0 everywhere when no band could be fitted. See GetMask. +void AzimuthalShapeFit(const std::vector §or_median, int n_bands, + std::vector &harm_c, std::vector &harm_s); + +// That modulation at a pixel: 1 + c cos 2phi + s sin 2phi, the band harmonics interpolated linearly in +// radius between band centres. (lx, ly) and the radius are the pixel's ring_position. +SHADOW_FINDER_HD float azimuthal_shape(const float *harm_c, const float *harm_s, int n_bands, + int radius, float lx, float ly) { + const float t = fminf(fmaxf(static_cast(radius) / HARMONIC_BAND_PX - 0.5f, 0.0f), + static_cast(n_bands - 1)); + const int b0 = static_cast(t); + const int b1 = b0 + 1 < n_bands ? b0 + 1 : b0; + const float w = t - static_cast(b0); + const float c = (1.0f - w) * harm_c[b0] + w * harm_c[b1]; + const float s = (1.0f - w) * harm_s[b0] + w * harm_s[b1]; + const float r2 = fmaxf(lx * lx + ly * ly, 1e-6f); + return 1.0f + c * (lx * lx - ly * ly) / r2 + s * 2.0f * lx * ly / r2; +} } // namespace shadow_finder diff --git a/image_analysis/beam_stop/ShadowMaskGPU.cu b/image_analysis/beam_stop/ShadowMaskGPU.cu index b57efc671..a3dd6d2c8 100644 --- a/image_analysis/beam_stop/ShadowMaskGPU.cu +++ b/image_analysis/beam_stop/ShadowMaskGPU.cu @@ -53,24 +53,6 @@ __device__ double poisson_deficit_sigma(double observed, double expected) { return ll > 0.0 ? sqrt(ll) : 0.0; } -// DiffractionGeometry::CalcAzIntPolarizationCorr about the centre the rings are drawn about. -__device__ float polarization_factor(const ShadowMaskSetup &s, float x, float y) { - const float u = (x - s.beam_x) * s.pixel_size_mm; - const float v = (y - s.beam_y) * s.pixel_size_mm; - const float *m = s.det_matrix; - const float lx = m[0] * u + m[1] * v + m[2] * s.distance_mm; - const float ly = m[3] * u + m[4] * v + m[5] * s.distance_mm; - const float lz = m[6] * u + m[7] * v + m[8] * s.distance_mm; - const float two_theta = atan2f(sqrtf(lx * lx + ly * ly), lz); - float phi = atan2f(ly, lx); - if (phi < 0) - phi += 2.0f * PI; - const float cos_2theta = cosf(two_theta); - const float cos_2theta_2 = cos_2theta * cos_2theta; - const float cos_2phi = cosf(2.0f * phi); - return 0.5f * (1.0f + cos_2theta_2 - s.polarization * cos_2phi * (1.0f - cos_2theta_2)); -} - __global__ void setup_kernel(ShadowMaskSetup s, const uint32_t *__restrict__ pixel_mask, const int64_t *__restrict__ sum_value, const uint32_t *__restrict__ valid_count, float *__restrict__ pol, char *__restrict__ valid, int *__restrict__ radius, @@ -81,8 +63,7 @@ __global__ void setup_kernel(ShadowMaskSetup s, const uint32_t *__restrict__ pix return; const int x = static_cast(i % s.width), y = static_cast(i / s.width); const RingPosition ring = ring_position(s.det_matrix, s.distance_mm / s.pixel_size_mm, x - s.beam_x, y - s.beam_y); - const float p = ring.obliquity - * (s.has_polarization ? polarization_factor(s, static_cast(x), static_cast(y)) : 1.0f); + const float p = ring.obliquity; pol[i] = p; float mean = 0.0f; char v = 0; @@ -197,6 +178,37 @@ __global__ void baseline_kernel(int rings, const uint64_t *__restrict__ keys, co baseline[r] = b; } +// The pooled background relative to its ring's baseline, 0 where unusable. +__global__ void relative_kernel(size_t n, const char *__restrict__ valid, const float *__restrict__ pooled, + const int *__restrict__ radius, const float *__restrict__ baseline, + float *__restrict__ relative) { + const size_t i = blockIdx.x * static_cast(blockDim.x) + threadIdx.x; + if (i >= n) return; + relative[i] = valid[i] ? pooled[i] / fmaxf(baseline[radius[i]], 1e-6f) : 0.0f; +} + +// The pooled background divided by its azimuthal modulation (GetMask), the factor kept with the +// obliquity in `pol`, and the ring keys of the result. +__global__ void shape_kernel(ShadowMaskSetup s, int n_bands, const float *__restrict__ shape_c, + const float *__restrict__ shape_s, const int *__restrict__ radius, char *__restrict__ valid, + float *__restrict__ pooled, float *__restrict__ pol, uint64_t *__restrict__ ring_key) { + const size_t n = static_cast(s.width) * s.height; + const size_t i = blockIdx.x * static_cast(blockDim.x) + threadIdx.x; + if (i >= n) return; + if (valid[i]) { + const RingPosition ring = ring_position(s.det_matrix, s.distance_mm / s.pixel_size_mm, + static_cast(i % s.width) - s.beam_x, + static_cast(i / s.width) - s.beam_y); + const float a = azimuthal_shape(shape_c, shape_s, n_bands, radius[i], ring.lx, ring.ly); + if (a > 0.0f) { + pooled[i] /= a; + pol[i] *= a; + } else + valid[i] = 0; + } + ring_key[i] = valid[i] ? (static_cast(radius[i]) << 32) | float_key(pooled[i]) : UINT64_MAX; +} + __global__ void low_kernel(size_t n, uint32_t frames, const char *__restrict__ valid, const float *__restrict__ pooled, const int32_t *__restrict__ pooled_count, const float *__restrict__ pol, const int *__restrict__ radius, const float *__restrict__ baseline, @@ -378,36 +390,11 @@ __global__ void sector_median_kernel(int sectors, const uint64_t *__restrict__ k median[k] = n >= MIN_SECTOR_PIXELS ? key_float(static_cast(keys[offset[k] + n / 2])) : -1.0; } -__global__ void dim_kernel(ShadowMaskSetup s, uint32_t frames, const char *__restrict__ valid, - const char *__restrict__ region, const float *__restrict__ ratio, - const float *__restrict__ deficit, const int *__restrict__ radius, - const float *__restrict__ harm_c, const float *__restrict__ harm_s, - const int32_t *__restrict__ pooled_count, const float *__restrict__ pol, - const float *__restrict__ pooled, const float *__restrict__ baseline, - char *__restrict__ dim) { - const size_t n = static_cast(s.width) * s.height; +__global__ void dim_kernel(size_t n, const char *__restrict__ valid, const char *__restrict__ region, + const float *__restrict__ ratio, const float *__restrict__ deficit, char *__restrict__ dim) { const size_t i = blockIdx.x * static_cast(blockDim.x) + threadIdx.x; if (i >= n) return; - char d = 0; - if (valid[i] && !region[i] && ratio[i] < PENUMBRA_RATIO) { - // cos 2phi and sin 2phi from the offset to the beam. - const RingPosition ring = ring_position(s.det_matrix, s.distance_mm / s.pixel_size_mm, - static_cast(i % s.width) - s.beam_x, - static_cast(i / s.width) - s.beam_y); - const float dx = ring.lx, dy = ring.ly; - const float r2 = fmaxf(dx * dx + dy * dy, 1e-6f); - const int band = radius[i] / HARMONIC_BAND_PX; - const float model = fminf(1.0f, 1.0f + harm_c[band] * (dx * dx - dy * dy) / r2 - + harm_s[band] * 2.0f * dx * dy / r2); - if (model >= 1.0f) { - d = deficit[i] > MIN_DEFICIT_SIGMA; - } else { - const double counted = static_cast(frames) * pooled_count[i] * pol[i]; - d = ratio[i] < PENUMBRA_RATIO * model - && poisson_deficit_sigma(pooled[i] * counted, baseline[radius[i]] * model * counted) > MIN_DEFICIT_SIGMA; - } - } - dim[i] = d; + dim[i] = valid[i] && !region[i] && ratio[i] < PENUMBRA_RATIO && deficit[i] > MIN_DEFICIT_SIGMA; } // Join a region across the module gaps it crosses, one line per thread (bridge_gaps in @@ -524,7 +511,7 @@ std::vector ShadowMaskOnDevice(const ShadowMaskSetup &setup, const std CudaDevicePtr d_pixel_mask(n), mask(n); e.Upload(d_pixel_mask.get(), pixel_mask); - // Mean projection over the polarization factor, usable pixels and radius from the beam centre. + // Mean projection over the obliquity, usable pixels and radius from the beam centre. CudaDevicePtr pol(n), pooled(n), ratio(n), deficit(n); CudaDevicePtr valid(n), low(n), region(n), a(n), b(n), c(n); CudaDevicePtr radius(n), root(n), count_a(n), count_b(n), max_radius(1); @@ -546,16 +533,43 @@ std::vector ShadowMaskOnDevice(const ShadowMaskSetup &setup, const std const double *pooled_sum = num; const int32_t *pooled_count = den; - // The rings, each sorted once; the baseline is an order statistic of them. + // The rings, each sorted once; the baseline is an order statistic of them. Made twice: once to measure + // the background's azimuthal modulation against, and again once that is divided out. CudaDevicePtr keys(n), sorted(n); + CudaDevicePtr ring_offset(rings + 1); + CudaDevicePtr baseline(rings); + const auto ring_baseline = [&] { + e.SortKeys(keys, sorted); + group_offsets<<>>(sorted, n, rings, ring_offset); + baseline_kernel<<>>(rings, sorted, ring_offset, baseline); + e.Check("baseline"); + }; pooled_kernel<<>>(n, pooled_sum, pooled_count, valid, radius, pooled, keys); e.Check("pooled"); + ring_baseline(); + + // The background's azimuthal modulation, measured and divided out (GetMask): the medians of each + // band's sectors of the background relative to its ring, over every valid pixel. + const int n_bands = max_r / HARMONIC_BAND_PX + 1; + const int n_sectors = n_bands * HARMONIC_SECTORS; + relative_kernel<<>>(n, valid, pooled, radius, baseline, ratio); + check(cudaMemsetAsync(region.get(), 0, n, stream), "memset"); // no pixel excluded + sector_key_kernel<<>>(setup, valid, region, radius, ratio, keys); + e.Check("sectors"); e.SortKeys(keys, sorted); - CudaDevicePtr ring_offset(rings + 1); - group_offsets<<>>(sorted, n, rings, ring_offset); - CudaDevicePtr baseline(rings); - baseline_kernel<<>>(rings, sorted, ring_offset, baseline); - e.Check("baseline"); + CudaDevicePtr sector_offset(n_sectors + 1); + group_offsets<<>>(sorted, n, n_sectors, sector_offset); + CudaDevicePtr sector_median(n_sectors); + sector_median_kernel<<>>(n_sectors, sorted, sector_offset, sector_median); + e.Check("sector medians"); + std::vector shape_c, shape_s; + AzimuthalShapeFit(e.Download(sector_median.get(), n_sectors), n_bands, shape_c, shape_s); + CudaDevicePtr d_shape_c(n_bands), d_shape_s(n_bands); + e.Upload(d_shape_c.get(), shape_c); + e.Upload(d_shape_s.get(), shape_s); + shape_kernel<<>>(setup, n_bands, d_shape_c, d_shape_s, radius, valid, pooled, pol, keys); + e.Check("shape"); + ring_baseline(); const auto host_baseline = e.Download(baseline.get(), rings); const auto offsets = e.Download(ring_offset.get(), rings + 1); std::vector ring_pixels(rings); @@ -594,24 +608,8 @@ std::vector ShadowMaskOnDevice(const ShadowMaskSetup &setup, const std final_kernel<<>>(n, root, count_a, c, region, mask); e.Check("fill"); - // The arm search: sector medians of the ratio, the harmonic of each band, the dim pixels. - const int n_bands = max_r / HARMONIC_BAND_PX + 1; - const int n_sectors = n_bands * HARMONIC_SECTORS; - sector_key_kernel<<>>(setup, valid, region, radius, ratio, keys); - e.Check("sectors"); - e.SortKeys(keys, sorted); - CudaDevicePtr sector_offset(n_sectors + 1); - group_offsets<<>>(sorted, n, n_sectors, sector_offset); - CudaDevicePtr sector_median(n_sectors); - sector_median_kernel<<>>(n_sectors, sorted, sector_offset, sector_median); - e.Check("sector medians"); - std::vector harm_c, harm_s; - HarmonicFit(e.Download(sector_median.get(), n_sectors), n_bands, harm_c, harm_s); - CudaDevicePtr d_harm_c(n_bands), d_harm_s(n_bands); - e.Upload(d_harm_c.get(), harm_c); - e.Upload(d_harm_s.get(), harm_s); - dim_kernel<<>>(setup, frames, valid, region, ratio, deficit, radius, d_harm_c, d_harm_s, - pooled_count, pol, pooled, baseline, a); + // The arm search: the dim pixels. + dim_kernel<<>>(n, valid, region, ratio, deficit, a); e.Check("dim"); e.Dilate(a, b, c, BRIDGE_PX); check(cudaMemcpyAsync(c.get(), b.get(), n, cudaMemcpyDeviceToDevice, stream), "copy"); diff --git a/image_analysis/beam_stop/ShadowMaskGPU.h b/image_analysis/beam_stop/ShadowMaskGPU.h index 16b641cbf..741118aaf 100644 --- a/image_analysis/beam_stop/ShadowMaskGPU.h +++ b/image_analysis/beam_stop/ShadowMaskGPU.h @@ -10,24 +10,21 @@ #include -// What the mask is drawn about: the detector, the beam centre, and the geometry the rings and the -// polarization factor are read off (ShadowFinder keeps it as a DiffractionGeometry; the device takes -// it as numbers). +// What the mask is drawn about: the detector, the beam centre, and the geometry the rings are read off +// (ShadowFinder keeps it as a DiffractionGeometry; the device takes it as numbers). struct ShadowMaskSetup { int width = 0, height = 0; float beam_x = 0.0f, beam_y = 0.0f; float det_matrix[9] = {}; // row major float pixel_size_mm = 0.0f; float distance_mm = 0.0f; - bool has_polarization = false; - float polarization = 0.0f; }; // ShadowFinder::GetMask on the device, from the projection ShadowAccumulatorGPU holds there. Step for // step the host's algorithm - the same pooling, ring medians, components, morphology and arm search - // and the same answer wherever the arithmetic is exact: every integer, comparison, sort and component -// is. What is not is the floating point the two compilers evaluate differently - the polarization -// factor's trigonometry, the Poisson test's logarithm and the azimuth of the arm search - so a pixel +// is. What is not is the floating point the two compilers evaluate differently - the ring geometry, the +// Poisson test's logarithm and the azimuth of the sectors - so a pixel // within a rounding of one of those thresholds can come out the other way. std::vector ShadowMaskOnDevice(const ShadowMaskSetup &setup, const std::vector &pixel_mask, const int64_t *max_value, const int64_t *sum_value, diff --git a/tests/ShadowFinderTest.cpp b/tests/ShadowFinderTest.cpp index 626fa902b..f253d036e 100644 --- a/tests/ShadowFinderTest.cpp +++ b/tests/ShadowFinderTest.cpp @@ -329,11 +329,13 @@ TEST_CASE("ShadowFinder_TheRingsFollowTheCentreTheCallerNames", "[ShadowFinder]" CHECK(at_measured[I((PATCH_X0 + PATCH_X1) / 2, (PATCH_Y0 + PATCH_Y1) / 2)] == 1); // About the true centre the rings are flat and only the patch and its penumbra come back; - // about the file's they cut across the fall-off and a large part of the detector does. + // about the file's they cut across the fall-off and more of the detector does. (Less than it + // used to: the part of the fall-off that is a second harmonic about the wrong centre is taken out + // with the background's azimuthal modulation.) const auto masked_at_file = std::count(at_file.begin(), at_file.end(), 1u); const auto masked_at_measured = std::count(at_measured.begin(), at_measured.end(), 1u); CHECK(masked_at_measured < 2 * PATCH_PX); - CHECK(masked_at_file > 3 * PATCH_PX); + CHECK(masked_at_file > masked_at_measured); } // A detector on a swung 2theta arm. The beam no longer lands on the geometry's centre (the point of @@ -377,8 +379,7 @@ TEST_CASE("ShadowFinder_ATiltedDetectorIsNotAShadow", "[ShadowFinder]") { // several times the dip this class looks for, so with the modulation left in, half of every outer // ring reads as shadow with nothing in the beam. The scene here is exactly that and nothing else - // a background carrying the Kahn factor of a horizontally polarized source, no hardware anywhere - -// and it must come back empty when the experiment states its polarization, and does not when it -// says nothing, which is what the correction is for. +// and it must come back empty. TEST_CASE("ShadowFinder_APolarizedBackgroundIsNotAShadow", "[ShadowFinder]") { // Short enough that this small detector reaches 2 theta = 70 degrees in its corner, where // the source suppresses its own plane to a fifth. A real geometry gets there with a short @@ -410,7 +411,11 @@ TEST_CASE("ShadowFinder_APolarizedBackgroundIsNotAShadow", "[ShadowFinder]") { frame[I(x, y)] = static_cast(std::lround(LEVEL * factor / at_peak)); } - auto mask_for = [&](const std::optional &polarization) { + // The modulation is measured from the background, so what the experiment states about the source + // - the beam's own value, another degree, none at all, or the plane turned by 90 degrees - is not + // read and cannot turn the lobes into a shadow. + for (const std::optional polarization : {std::optional(POLARIZATION), std::optional(0.9f), + std::optional(), std::optional(-POLARIZATION)}) { DiffractionExperiment x = TestExperiment(); x.DetectorDistance_mm(DISTANCE_MM).PolarizationFactor(polarization); const PixelMask pixel_mask(x); @@ -422,20 +427,65 @@ TEST_CASE("ShadowFinder_APolarizedBackgroundIsNotAShadow", "[ShadowFinder]") { finder.AddImage(msg, buffer); } const auto mask = finder.GetMask(); - return std::count(mask.begin(), mask.end(), 1u); - }; - - CHECK(mask_for(POLARIZATION) == 0); - CHECK(mask_for(std::nullopt) > static_cast(W) * H / 20); + CHECK(std::count_if(mask.begin(), mask.end(), [](uint32_t v) { return v != 0; }) == 0); + } } -// 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. +// A beam stop with four arms on a polarized background at a short distance: the stop has to come back, +// two of its arms inside the background's dim lobes, the lobes themselves must not, and the mask must not move when the experiment states another +// degree of polarization or none. The degree cannot be told from the data, so a mask that followed +// it would follow a guess. +TEST_CASE("ShadowFinder_TheMaskDoesNotDependOnTheStatedPolarization", "[ShadowFinder]") { + constexpr float DISTANCE_MM = 8.0f; + constexpr double TRUE_POLARIZATION = -0.998; // the plane along the image's y + constexpr int LEVEL = 24; + + const float pixel_mm = TestExperiment().GetPixelSize_mm(); + const auto stop = Scene(/*cross=*/true, false); + std::vector frame(static_cast(W) * H); + for (int y = 0; y < H; y++) + for (int xi = 0; xi < W; xi++) { + const double dx = xi - C, dy = y - C; + const double rho = std::hypot(dx, dy) * pixel_mm; + const double cos_2theta_sq = DISTANCE_MM * DISTANCE_MM / (DISTANCE_MM * DISTANCE_MM + rho * rho); + const double rr = dx * dx + dy * dy; + const double cos_2phi = rr > 0.0 ? (dx * dx - dy * dy) / rr : 0.0; + const double factor = 0.5 * (1.0 + cos_2theta_sq - TRUE_POLARIZATION * cos_2phi * (1.0 - cos_2theta_sq)); + frame[I(xi, y)] = stop[I(xi, y)] == 0 ? 0 : static_cast(std::lround(LEVEL * factor)); + } + + auto mask_for = [&](const std::optional &polarization) { + DiffractionExperiment x = TestExperiment(); + x.DetectorDistance_mm(DISTANCE_MM).PolarizationFactor(polarization); + 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); + } + return finder.GetMask(); + }; + + const auto mask = mask_for(0.99f); + CHECK(mask[I(C, C)] == ShadowFinder::SHADOW); // the stop + CHECK(mask[I(W - 3, C)] == ShadowFinder::SHADOW); // and its arms, along x + CHECK(mask[I(C, 3)] == ShadowFinder::SHADOW); // and along y, inside the suppressed lobe + CHECK(mask[I(C + 20, 3)] == 0); // but not the lobe itself + CHECK(mask[I(C - 20, H - 4)] == 0); + CHECK(mask_for(0.998f) == mask); + CHECK(mask_for(std::nullopt) == mask); + CHECK(mask_for(-0.998f) == mask); +} + +// An unpolarized background in an experiment that states 0.99. When the finder divided by the stated +// factor, every outer ring came out 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 modulation is +// measured from the background now, so there is no lobe, 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; @@ -554,9 +604,9 @@ namespace { } // The mask is made on the GPU wherever the projection is there. On scenes where no pixel sits within a -// rounding of a threshold the two must agree to the pixel: every step but the polarization factor, -// the Poisson test's logarithm and the arm search's azimuth is exact on both, and these scenes are -// built so that none of the three decides anything at an edge. +// rounding of a threshold the two must agree to the pixel: every step but the ring geometry, the +// Poisson test's logarithm and the sectors' azimuth is exact on both, and these scenes are built so +// that none of the three decides anything at an edge. TEST_CASE("ShadowFinder_DeviceMaskMatchesHost", "[ShadowFinder]") { if (get_gpu_count() == 0) SKIP("no GPU");