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");