From fc86fd4a10e5d3c6e084d49681439943ec868161 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 23 Oct 2026 06:41:10 +0200 Subject: [PATCH] ShadowFinder: back to the predicted polarization factor, in the plane the background shows The measured-modulation finder of the previous three commits is withdrawn. On a CCD sweep with a holder arm along x it lost the arm: a small rise in the deep pixels' count promoted part of the arm to shadow, which left the rest of the arm too few deep pixels to pass as transmitting hardware, and the mask fell from 2.5 % to 0.4 % of the detector. On a 16-Mpixel stills detector the per-band fit went wrong in the outer bands and masked 3.6 % of the detector instead of 0.03 %. The shadow finder is rc178's again. The stills defect it was meant to fix came from dividing by the factor of the wrong plane: the background of that run is suppressed along the image's y axis, not lab x, and the division left dim lobes over the corners on the arm search's thresholds, so that the third decimal of p moved thousands of pixels. The degree cannot be told from the data, but the plane can. Before the correction the mean projection is compared with the profile a polarized source gives each radial band (rings relative to their upper decile, each sector the brighter of itself and the one opposite, bands seen over three quarters of the azimuth); the stated degree is applied in the plane the background follows when it follows one clearly (correlation > 0.5, at least 0.3 of a fully polarized source's modulation), and as stated otherwise. Over 46 rotation sweeps and 11 stills runs only that one run is turned, so every other mask is unchanged; there the mask is the 5248-pixel stop at p = 0.99, 0.995 and 0.998. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu --- image_analysis/beam_stop/SHADOW_FINDER.md | 77 ++-- image_analysis/beam_stop/ShadowFinder.cpp | 402 ++++++++++-------- image_analysis/beam_stop/ShadowFinder.h | 16 +- .../beam_stop/ShadowFinderInternal.h | 49 +-- image_analysis/beam_stop/ShadowMaskGPU.cu | 104 ++--- image_analysis/beam_stop/ShadowMaskGPU.h | 11 +- tests/ShadowFinderTest.cpp | 102 +++-- 7 files changed, 385 insertions(+), 376 deletions(-) diff --git a/image_analysis/beam_stop/SHADOW_FINDER.md b/image_analysis/beam_stop/SHADOW_FINDER.md index e3ccca671..998947ea8 100644 --- a/image_analysis/beam_stop/SHADOW_FINDER.md +++ b/image_analysis/beam_stop/SHADOW_FINDER.md @@ -44,30 +44,20 @@ 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. **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's upper decile (the median of a ring half under a shadow sits on the shadow's - edge), each 64 px radial band takes the median of each of its 24 sectors, and each sector is - replaced by the brighter of itself and the sector opposite - the modulation is the same at `phi` and - `phi + 180 deg`, a shadow almost never is, and shadows only darken. A band seen over at least three - quarters of the azimuth that way is fitted with `1 + c cos 2phi + s sin 2phi` (sectors dim against - the fit 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 or beside a wide shadow - taking 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. +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. 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 @@ -119,14 +109,15 @@ 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. - "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 - and then against the ring as it still varies in azimuth: each band's - sector medians of the ratio, the region already found left out, are fitted with a second harmonic - and a pixel is dim below `penumbra_ratio` times that model where the model is below 1. A wide, - smooth shadow over a large arc on one side is explained by it rather than added as a piece. The - model is capped at 1, so it only ever removes dim pixels. + 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. 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 @@ -206,11 +197,23 @@ 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 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. +- **The polarization degree is taken from the experiment; the plane, where the background shows one, + from the background.** 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 degree cannot be told from the data, but the plane can, and it is not always the + one assumed: a detector mounted a quarter turn from another puts the suppressed direction along the + image's other axis. Divided by the factor of the wrong plane, the background comes out less flat + than it went in - on one stills run the outer rings carried dim lobes over the detector's corners, + sitting on the arm search's thresholds, so that the third decimal of the stated degree moved + thousands of pixels in and out of the mask. So before the correction the mean projection is + compared ring by ring with the profile a polarized source gives each 64 px band (each ring relative + to its upper decile, each sector the brighter of itself and the one opposite, since a shadow is + rarely on both sides; a band counted where three quarters of the azimuth are seen). The stated + degree is applied in the plane the background follows when it follows one clearly - a correlation + above 0.5 and at least 0.3 of a fully polarized source's modulation - and as stated otherwise. Over + the rotation and stills corpora only the run whose plane lies along the image's y axis is turned; on + every other run the correction, and so the mask, is what it was. An experiment that states no + polarization is left uncorrected, as before. - **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 6911195e3..a042a40c6 100644 --- a/image_analysis/beam_stop/ShadowFinder.cpp +++ b/image_analysis/beam_stop/ShadowFinder.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include @@ -36,6 +37,16 @@ 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) { @@ -389,76 +400,84 @@ 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 plane the background shows (PolarizationPlane): +1, -1, or 0 where it does not show one. +// The degree of polarization cannot be told from the data, but the plane can, and it is not always the +// one the experiment assumes - a detector mounted a quarter turn from another on the same kind of +// beamline puts the suppressed direction along the image's other axis. Divided by the factor of the +// wrong plane, the background is made less flat rather than more: on a stills run the outer rings +// then carried dim lobes over the detector's corners that sat on the arm search's thresholds, so that +// the third decimal of p moved thousands of pixels in and out of the mask. +// `mean` is the mean projection (W x H, NaN where unusable), compared ring by ring about the centre. +float measured_plane(const std::vector &mean, int W, int H, float centre_x, float centre_y, + const float *det_matrix, float distance_px, size_t nthreads) { + const int n_pixels = W * H; + Plane value(n_pixels); + Plane valid(n_pixels); + Plane radius(n_pixels); + std::atomic max_radius_atomic{0}; + ParallelChunks(H, nthreads, [&](int ylo, int yhi) { + int local_max = 0; + for (int y = ylo; y < yhi; y++) + for (int x = 0; x < W; x++) { + const int i = y * W + x; + const RingPosition ring = ring_position(det_matrix, distance_px, x - centre_x, y - centre_y); + valid[i] = std::isfinite(mean[i]) && ring.obliquity > 0.0f; + value[i] = valid[i] ? mean[i] / ring.obliquity : 0.0f; + radius[i] = static_cast(std::lround(ring.radius)); + local_max = std::max(local_max, radius[i]); + } + int prev = max_radius_atomic.load(); + 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; -// 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) { + // Each ring is taken relative to its upper decile rather than its median: where a shadow covers + // half of a ring the median sits on the shadow's edge and jumps from ring to ring. + const RingValues rings = bin_by_ring(value, valid, radius, max_radius, nthreads); + std::vector level(max_radius + 1, 0.0f); + for (int r = 0; r <= max_radius; r++) { + const int n = rings.offset[r + 1] - rings.offset[r]; + if (n > 0) + level[r] = rings.values[rings.offset[r] + n * 9 / 10]; + } + + // The median of each sector of each band, relative to its rings. 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); + for (int i = 0; i < n_pixels; i++) { + if (!valid[i] || level[radius[i]] <= 0.0f) + 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)); + sector_values[static_cast(radius[i] / HARMONIC_BAND_PX) * HARMONIC_SECTORS + sector] + .push_back(value[i] / level[radius[i]]); + } + 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()); - median[k] = v[v.size() / 2]; + sector_median[k] = v[v.size() / 2]; } }); - return median; + + // The modulation a fully polarized source gives each band, at the 2 theta of its centre. + std::vector amplitude(n_bands); + for (int band = 0; band < n_bands; band++) { + const double t = (band + 0.5) * HARMONIC_BAND_PX / distance_px; // tan 2theta + const double cos2 = 1.0 / (1.0 + t * t); // cos^2 2theta + amplitude[band] = static_cast((1.0 - cos2) / (1.0 + cos2)); + } + return PolarizationPlane(sector_median, amplitude); +} + +// The stated polarization in the plane the background shows; as stated where it shows none. A negative +// coefficient is a vertically polarized source (CalcAzIntPolarizationCorr). +float PlaneApplied(float stated, float plane) { + return plane == 0.0f ? stated : std::fabs(stated) * plane; } } // namespace @@ -469,6 +488,7 @@ 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) @@ -731,6 +751,12 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { setup.pixel_size_mm = geometry.GetPixelSize_mm(); setup.distance_mm = geometry.GetDetectorDistance_mm(); ScopedGpuPin pin(gpu->GetDevice()); + setup.has_polarization = polarization.has_value(); + if (polarization) + setup.polarization = PlaneApplied(*polarization, + measured_plane(gpu->MeanProjection(pixel_mask), width, height, beam_x, + beam_y, setup.det_matrix, + setup.distance_mm / setup.pixel_size_mm, nthreads)); return gpu->Mask(setup, pixel_mask); } #endif @@ -777,16 +803,32 @@ 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 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. + // 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); const std::vector det_matrix = geometry.GetDetectorMatrix().arr(); const float distance_px = geometry.GetDetectorDistance_mm() / geometry.GetPixelSize_mm(); + // The stated degree, in the plane the background shows (measured_plane). + std::optional polarization_applied; + if (polarization) { + std::vector projection(n_pixels); + ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) { + for (int i = lo; i < hi; i++) + projection[i] = valid_count[i] > 0 && known_mask[i] == 0 + ? static_cast(static_cast(sum_value[i]) / valid_count[i]) : NAN; + }); + polarization_applied = PlaneApplied(*polarization, measured_plane(projection, W, H, centre_x, centre_y, + det_matrix.data(), distance_px, nthreads)); + } + Plane mean(n_pixels); - Plane shape(n_pixels); + Plane pol(n_pixels); Plane valid(n_pixels); Plane radius(n_pixels); std::atomic max_radius_atomic{0}; @@ -796,14 +838,21 @@ 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); - 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. + pol[i] = ring.obliquity + * (polarization_applied.has_value() + ? pol_geometry.CalcAzIntPolarizationCorr(static_cast(x), + static_cast(y), + polarization_applied.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. mean[i] = 0.0f; valid[i] = 0; - if (valid_count[i] > 0 && known_mask[i] == 0 && shape[i] > 0.0f) { + if (valid_count[i] > 0 && known_mask[i] == 0 && pol[i] > 0.0f) { mean[i] = static_cast(static_cast(sum_value[i]) - / valid_count[i] / shape[i]); + / valid_count[i] / pol[i]); valid[i] = 1; } radius[i] = static_cast(std::lround(ring.radius)); @@ -815,7 +864,6 @@ 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 @@ -843,65 +891,37 @@ 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. - { - // Each ring is taken relative to its upper decile rather than its median: where a shadow covers - // half of a ring the median sits on the shadow's edge and jumps from ring to ring, and the - // rings of one band then disagree on what 1 means. - const RingValues level_rings = bin_by_ring(pooled, valid, radius, max_radius, nthreads); - std::vector level(max_radius + 1, 0.0f); - for (int r = 0; r <= max_radius; r++) { - const int n = level_rings.offset[r + 1] - level_rings.offset[r]; - if (n > 0) - level[r] = level_rings.values[level_rings.offset[r] + n * 9 / 10]; - } - 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); - const std::vector baseline = ring_baseline(rings, 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; + }); + } Plane ratio(n_pixels); ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) { for (int i = lo; i < hi; i++) @@ -945,7 +965,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] * shape[i]; + const double counted = static_cast(frames) * pooled_count[i] * pol[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; @@ -1048,30 +1068,58 @@ 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. The - // background's modulation is divided out above, but a wide shadow that dims a ring smoothly over a - // large arc - the shadow of a goniometer turning into the beam on one side of the detector - leaves - // dim sectors that are not hardware standing in a few of them; joined, their edge made one piece - // over a large part of the detector. A second harmonic fitted to the sector medians of the ratio, - // the sectors already in the region left out, explains such an arc. The harmonic is only ever + // 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. - Plane outside_region(n_pixels); - ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) { - for (int i = lo; i < hi; i++) - outside_region[i] = valid[i] && !region[i]; + 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; - std::vector harm_fitted; - HarmonicFit(sector_medians(ratio, outside_region, radius, W, centre_x, centre_y, det_matrix.data(), distance_px, - n_bands, nthreads), - n_bands, harm_c, harm_s, harm_fitted); + HarmonicFit(sector_median, n_bands, harm_c, harm_s); - Plane dim(n_pixels); + Plane dim = filled_plane(n_pixels, 0, nthreads); ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) { for (int i = lo; i < hi; i++) { - dim[i] = 0; - if (!outside_region[i] || ratio[i] >= PENUMBRA_RATIO) + 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, @@ -1086,7 +1134,7 @@ std::vector ShadowFinder::MaskOn(const Projection &p, const std::vecto dim[i] = deficit[i] > MIN_DEFICIT_SIGMA; continue; } - const double counted = static_cast(frames) * pooled_count[i] * shape[i]; + 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; } @@ -1150,10 +1198,9 @@ 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 &fitted) { + std::vector &harm_c, std::vector &harm_s) { 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++) { @@ -1179,7 +1226,6 @@ 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; @@ -1187,29 +1233,25 @@ 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) { - // A polarized background is the same at phi and at phi + 180 deg; what shadows the detector rarely - // is - an arm, a pin or a goniometer stands on one side of the beam. Shadows only darken, so each - // sector is taken as the brighter of itself and the sector opposite, and a shadow that covers a - // wide arc on one side - a third of the ring at a detector's edge on some sweeps - is not fitted - // as a harmonic of the background. - // A band is fitted only where three quarters of the azimuth are seen this way: an arc - the corner - // of a detector, or what a wide shadow leaves - says too little about a harmonic to be believed. - std::vector folded(sector_median); +float PolarizationPlane(const std::vector §or_median, const std::vector &litude) { + const int n_bands = static_cast(amplitude.size()); + // The modulation is the same at phi and at phi + 180 deg; what shadows the detector rarely is - an + // arm, a pin or a goniometer stands on one side of the beam - and shadows only darken. So each + // sector is taken as the brighter of itself and the sector opposite, and a band is fitted only where + // three quarters of the azimuth are seen that way. constexpr int HALF = HARMONIC_SECTORS / 2; + std::vector folded(sector_median); + std::vector seen_enough(n_bands); for (int band = 0; band < n_bands; band++) { int seen = 0; for (int k = 0; k < HALF; k++) { @@ -1217,27 +1259,31 @@ void AzimuthalShapeFit(const std::vector §or_median, int n_bands, folded[a] = folded[a + HALF] = std::max(sector_median[a], sector_median[a + HALF]); seen += folded[a] >= 0; } - if (seen < HALF * 3 / 4) - std::fill(folded.begin() + static_cast(band) * HARMONIC_SECTORS, - folded.begin() + static_cast(band + 1) * HARMONIC_SECTORS, -1.0); + seen_enough[band] = seen >= HALF * 3 / 4; } std::vector c, s; - std::vector fitted; - HarmonicFit(folded, n_bands, c, s, fitted); - harm_c.assign(n_bands, 0.0f); - harm_s.assign(n_bands, 0.0f); + HarmonicFit(folded, n_bands, c, s); + // A horizontally polarized source makes the background dim along lab x, c < 0, by more where the + // modulation it can give is larger; a vertically polarized one makes it dim along lab y, c > 0. A + // plane is taken only where the bands' harmonics follow its profile (PLANE_CORRELATION) and amount + // to a good part of what a fully polarized source gives (PLANE_MODULATION): a background with + // little modulation of its own decides nothing. With no modulation at all both are NaN, which + // decides nothing too. + double ca = 0.0, cc = 0.0, aa = 0.0; 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; - } - } + if (!seen_enough[band]) + continue; + ca += c[band] * amplitude[band]; + cc += c[band] * c[band]; + aa += amplitude[band] * amplitude[band]; } + const double correlation = ca / std::sqrt(cc * aa); + const double modulation = ca / aa; + if (correlation > PLANE_CORRELATION && modulation > PLANE_MODULATION) + return -1.0f; + if (correlation < -PLANE_CORRELATION && modulation < -PLANE_MODULATION) + return 1.0f; + return 0.0f; } } // namespace shadow_finder diff --git a/image_analysis/beam_stop/ShadowFinder.h b/image_analysis/beam_stop/ShadowFinder.h index 35515844f..1c844f5cd 100644 --- a/image_analysis/beam_stop/ShadowFinder.h +++ b/image_analysis/beam_stop/ShadowFinder.h @@ -48,11 +48,19 @@ class ShadowFinder { float beam_x; float beam_y; - // 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. + // 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. 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 ac2813f98..af89c80bd 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 its azimuthal modulation. +// background is divided by it as by the polarization factor. // (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,9 +87,10 @@ inline constexpr int MIN_SHADOW_PIXELS = 2000; // and deep nowhere. inline constexpr int MIN_CORE_PIXELS = 200; -// 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 +// The arm search compares a pixel with its ring as the ring actually varies around the beam. A ring +// of background is not flat once divided by the polarization factor when that factor is not the +// beam's: the remainder is a second harmonic in azimuth, cos 2phi, which reaches tens of percent at +// high angle. It is measured over radial bands of HARMONIC_BAND_PX, from the median of each of // HARMONIC_SECTORS sectors that holds at least MIN_SECTOR_PIXELS pixels. inline constexpr int HARMONIC_BAND_PX = 64; inline constexpr int HARMONIC_SECTORS = 24; @@ -115,30 +116,26 @@ 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). A band whose sectors cannot be fitted gets 0 and fitted[band] = 0. See GetMask. +// has too few pixels). See GetMask. void HarmonicFit(const std::vector §or_median, int n_bands, - std::vector &harm_c, std::vector &harm_s, std::vector &fitted); + std::vector &harm_c, std::vector &harm_s); -// The background's own azimuthal modulation, band by band: HarmonicFit to the brighter of each sector -// and the sector opposite, with a band that could not be fitted - inside the stop, or where a band is -// only an arc, at a detector's corners or beside a wide shadow - 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); +// How closely, over the radial bands, the background's cos 2phi harmonic must follow the profile a +// polarized source gives it (the correlation), and how large it must be against the profile of a fully +// polarized source (the least-squares scale), before the plane is taken from the background rather +// than as stated. On 46 rotation sweeps and 11 stills runs whose plane is horizontal, the largest values +// towards the vertical plane are 0.48 and 0.18, the latter on a background with almost no modulation; +// on a stills run whose plane lies along the image's y axis they are 0.998 and 0.47 - the scale reduced +// there by the placeholder distance the pre-scan reads it at. +inline constexpr double PLANE_CORRELATION = 0.5; +inline constexpr double PLANE_MODULATION = 0.3; -// 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; -} +// The plane the source is polarized in, as the background shows it: +1 for the plane containing lab x +// (a horizontally polarized source), -1 for the plane turned by 90 degrees, 0 where the background does +// not show one clearly. `sector_median` is laid +// out as for HarmonicFit, each band's background relative to its rings; `amplitude[band]` is the +// relative modulation a fully polarized source gives the band, sin^2 2theta / (1 + cos^2 2theta). See +// GetMask. +float PolarizationPlane(const std::vector §or_median, const std::vector &litude); } // namespace shadow_finder diff --git a/image_analysis/beam_stop/ShadowMaskGPU.cu b/image_analysis/beam_stop/ShadowMaskGPU.cu index d0ced685b..b57efc671 100644 --- a/image_analysis/beam_stop/ShadowMaskGPU.cu +++ b/image_analysis/beam_stop/ShadowMaskGPU.cu @@ -53,6 +53,24 @@ __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, @@ -63,7 +81,8 @@ __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; + const float p = ring.obliquity + * (s.has_polarization ? polarization_factor(s, static_cast(x), static_cast(y)) : 1.0f); pol[i] = p; float mean = 0.0f; char v = 0; @@ -178,46 +197,6 @@ __global__ void baseline_kernel(int rings, const uint64_t *__restrict__ keys, co baseline[r] = b; } -// The ring's upper decile, what the background's modulation is measured relative to (GetMask). -__global__ void level_kernel(int rings, const uint64_t *__restrict__ keys, const int *__restrict__ offset, - float *__restrict__ level) { - const int r = blockIdx.x * blockDim.x + threadIdx.x; - if (r >= rings) return; - const int n = offset[r + 1] - offset[r]; - level[r] = n > 0 ? key_float(static_cast(keys[offset[r] + n * 9 / 10])) : 0.0f; -} - -// The pooled background relative to its ring's level, 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__ level, - 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(level[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, @@ -545,7 +524,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 obliquity, usable pixels and radius from the beam centre. + // Mean projection over the polarization factor, 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); @@ -567,42 +546,14 @@ 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 level the background's azimuthal modulation is measured against, - // and the baseline once that is divided out, are order statistics of them. + // The rings, each sorted once; the baseline is an order statistic of them. CudaDevicePtr keys(n), sorted(n); - CudaDevicePtr ring_offset(rings + 1); - CudaDevicePtr baseline(rings); pooled_kernel<<>>(n, pooled_sum, pooled_count, valid, radius, pooled, keys); e.Check("pooled"); e.SortKeys(keys, sorted); + CudaDevicePtr ring_offset(rings + 1); group_offsets<<>>(sorted, n, rings, ring_offset); - CudaDevicePtr level(rings); - level_kernel<<>>(rings, sorted, ring_offset, level); - e.Check("level"); - - // The background's azimuthal modulation, measured and divided out (GetMask): the medians of each - // band's sectors of the background relative to its ring's level, 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, level, 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 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"); - e.SortKeys(keys, sorted); - group_offsets<<>>(sorted, n, rings, ring_offset); + CudaDevicePtr baseline(rings); baseline_kernel<<>>(rings, sorted, ring_offset, baseline); e.Check("baseline"); const auto host_baseline = e.Download(baseline.get(), rings); @@ -644,15 +595,18 @@ std::vector ShadowMaskOnDevice(const ShadowMaskSetup &setup, const std 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; - std::vector harm_fitted; - HarmonicFit(e.Download(sector_median.get(), n_sectors), n_bands, harm_c, harm_s, harm_fitted); + 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); diff --git a/image_analysis/beam_stop/ShadowMaskGPU.h b/image_analysis/beam_stop/ShadowMaskGPU.h index 741118aaf..16b641cbf 100644 --- a/image_analysis/beam_stop/ShadowMaskGPU.h +++ b/image_analysis/beam_stop/ShadowMaskGPU.h @@ -10,21 +10,24 @@ #include -// 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). +// 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). 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 ring geometry, the -// Poisson test's logarithm and the azimuth of the sectors - so a pixel +// 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 // 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 acc9bc8ed..8966e4e43 100644 --- a/tests/ShadowFinderTest.cpp +++ b/tests/ShadowFinderTest.cpp @@ -331,13 +331,11 @@ 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 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.) + // about the file's they cut across the fall-off and a large part of the detector does. 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 > masked_at_measured); + CHECK(masked_at_file > 3 * PATCH_PX); } // A detector on a swung 2theta arm. The beam no longer lands on the geometry's centre (the point of @@ -381,7 +379,8 @@ 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. +// 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. 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 @@ -413,11 +412,7 @@ TEST_CASE("ShadowFinder_APolarizedBackgroundIsNotAShadow", "[ShadowFinder]") { frame[I(x, y)] = static_cast(std::lround(LEVEL * factor / at_peak)); } - // 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)}) { + auto mask_for = [&](const std::optional &polarization) { DiffractionExperiment x = TestExperiment(); x.DetectorDistance_mm(DISTANCE_MM).PolarizationFactor(polarization); const PixelMask pixel_mask(x); @@ -429,15 +424,19 @@ TEST_CASE("ShadowFinder_APolarizedBackgroundIsNotAShadow", "[ShadowFinder]") { finder.AddImage(msg, buffer); } const auto mask = finder.GetMask(); - CHECK(std::count_if(mask.begin(), mask.end(), [](uint32_t v) { return v != 0; }) == 0); - } + return std::count(mask.begin(), mask.end(), 1u); + }; + + CHECK(mask_for(POLARIZATION) == 0); + CHECK(mask_for(std::nullopt) > static_cast(W) * H / 20); } -// 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]") { +// A beam stop with four arms on a background polarized in the plane along the image's y, at a short +// distance, in an experiment that states the usual horizontal plane. Divided by the factor of the +// stated plane the background's lobes read as shadow, and how much of them did followed the third +// decimal of the stated degree. The plane is taken from the background, so the stop and its arms come +// back, the lobes do not, and the mask is the same for any stated degree and either stated plane. +TEST_CASE("ShadowFinder_ThePolarizationPlaneIsTheOneTheBackgroundShows", "[ShadowFinder]") { constexpr float DISTANCE_MM = 8.0f; constexpr double TRUE_POLARIZATION = -0.998; // the plane along the image's y constexpr int LEVEL = 24; @@ -477,44 +476,43 @@ TEST_CASE("ShadowFinder_TheMaskDoesNotDependOnTheStatedPolarization", "[ShadowFi 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); } -// The modulation is measured from the sector medians of each band. A wide shadow on one side of the -// beam - a third of the ring dimmed to 70 %, too little to be left out as hardware - must not be -// fitted as part of it: the modulation is the same on opposite sides of the beam and the shadow is -// not. A band seen only as an arc is not fitted and takes its inner neighbour's harmonic. -TEST_CASE("ShadowFinder_AOneSidedShadowIsNotPartOfTheModulation", "[ShadowFinder]") { +// The plane is decided from the sector medians of each band: the cos 2phi harmonic has to follow the +// profile a polarized source gives it, with a one-sided shadow over a third of the ring left out by +// comparing opposite sectors. A background without modulation decides nothing. +TEST_CASE("ShadowFinder_PolarizationPlane", "[ShadowFinder]") { using namespace shadow_finder; - constexpr double C2 = -0.3; // the background's own modulation: dim along x - const int n_bands = 3; - std::vector medians(static_cast(n_bands) * HARMONIC_SECTORS); - for (int band = 0; band < n_bands; band++) - for (int k = 0; k < HARMONIC_SECTORS; k++) { - const double phi = (k + 0.5) * 2.0 * std::numbers::pi / HARMONIC_SECTORS - std::numbers::pi; - double m = 1.0 + C2 * std::cos(2.0 * phi); - if (k >= 8 && k < 16) - m *= 0.7; // the shadow, a third of the ring around phi = 0 - if (band == 2 && (k < 6 || k >= 12)) - m = -1.0; // the outer band is only an arc - medians[static_cast(band) * HARMONIC_SECTORS + k] = m; - } - std::vector c, s; - AzimuthalShapeFit(medians, n_bands, c, s); - for (int band = 0; band < n_bands; band++) { - CHECK(c[band] == Catch::Approx(C2).margin(0.01)); - CHECK(s[band] == Catch::Approx(0.0).margin(0.01)); - } + constexpr int N_BANDS = 8; + std::vector amplitude(N_BANDS); + for (int band = 0; band < N_BANDS; band++) + amplitude[band] = 0.05f * static_cast(band + 1); + auto medians = [&](double sign, bool shadow) { + std::vector m(static_cast(N_BANDS) * HARMONIC_SECTORS); + for (int band = 0; band < N_BANDS; band++) + for (int k = 0; k < HARMONIC_SECTORS; k++) { + const double phi = (k + 0.5) * 2.0 * std::numbers::pi / HARMONIC_SECTORS - std::numbers::pi; + double v = 1.0 + sign * amplitude[band] * std::cos(2.0 * phi); + if (shadow && k >= 8 && k < 16) + v *= 0.7; // a third of the ring around phi = 0 + m[static_cast(band) * HARMONIC_SECTORS + k] = v; + } + return m; + }; + CHECK(PolarizationPlane(medians(-1.0, false), amplitude) == 1.0f); // dim along x: horizontal + CHECK(PolarizationPlane(medians(+1.0, false), amplitude) == -1.0f); // dim along y: vertical + CHECK(PolarizationPlane(medians(+1.0, true), amplitude) == -1.0f); // the shadow along x does not hide it + CHECK(PolarizationPlane(medians(0.0, false), amplitude) == 0.0f); // no modulation: no decision } -// 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. +// The other side of that correction: a polarization factor that is not the beam's. Here the +// background is unpolarized and the experiment states 0.99, so after the division every outer ring is +// bright along x and dim along y - a smooth lobe below PENUMBRA_RATIO across the top and bottom of the +// detector with nothing in the beam. On its own the lobe is never deep, but a small absorbing speck +// inside it is, and the transmitting-hardware step joined the two into one piece "dim everywhere and +// deep in places" and masked the whole lobe. The step compares a pixel with its ring's own second +// harmonic, so the lobe is explained away and the speck, far below a shadow's size, stays unmasked. TEST_CASE("ShadowFinder_AWrongPolarizationFactorIsNotAnArm", "[ShadowFinder]") { constexpr float DISTANCE_MM = 6.0f; constexpr int LEVEL = 24; @@ -633,9 +631,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 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. +// 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. TEST_CASE("ShadowFinder_DeviceMaskMatchesHost", "[ShadowFinder]") { if (get_gpu_count() == 0) SKIP("no GPU");