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