mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 82 Wiki Activity

Beam stop: rings of constant 2theta and the tilt's solid angle divided out, so a swung detector is not read as shadow - RESULT-CHANGING

Browse Source 
ShadowFinder compared each pixel with the median of a 1-pixel circle about the
geometry's beam centre. On a tilted detector that centre is the point of normal
incidence, not where the beam lands, and the rings of constant 2theta are not
circles about it: on a 2theta arm swung by 5-10 deg at 50 mm the circles cut across
the background's steep fall-off and the far side of the detector read as shadow.
Every tilted home-source SMV set in the open arm masked 14-41 % of its detector.

A pixel's ring is now distance * tan(2theta) from the detector matrix, its azimuth
for the arm search the lab azimuth, and the mean is divided by the obliquity
(distance / l_z)^3 - the part of the solid angle that varies along a 2theta ring -
alongside the polarization factor. Host and device share one helper. Untilted, all
three are exact identities, so untilted sets keep their mask bit for bit.

Named-set before/after on the same commit (16 sets): 5reo, 5uth, 6h2p_1p89A, 4bwl,
lyso_x06da_ref unchanged. Shadow 26.5->5.9 % (3mc4), 40.6->10.4 (3r6o),
39.4->6.4 (5vml), 37.1->9.2 (3meb), 33.1->6.8 (5cc8), 13.8->5.8 (3p85);
completeness 77->91, 43->98, 80->84, 43->78, 45->80, 86->87 %. The 0.2 deg-tilted
in-house PILATUS4 sets move by a few hundred pixels; the 20/30 deg tilted
small-molecule sets mask 0.4->0.8 and 0.1->0.6 %, SHELXL R1 unchanged or better.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4
This commit is contained in:
leonarski_f
2026-10-15 18:13:53 +02:00
co-authored by Claude Opus 5.5
parent 873880bba1
commit c207012a13
7 changed files with 131 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/beam_stop/SHADOW_FINDER.md
+12 -6
View File
@@ -50,13 +50,19 @@ own radius**, not to a global threshold.
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. It is the **only** correction a
ring carries that varies along it — solid angle, detector and air absorption are functions of
`2 theta` alone and the ring's own median absorbs them. The Poisson test of step 5 puts the
factor back, so the significance is still the significance of the counts that were recorded.
Nothing is corrected when the experiment states no polarization.
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 centre. The median ignores the shadow (a minority of each ring) and the spots.
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
arm) a circle about the centre runs across the background's fall-off, and at 5 deg and 50 mm
that put a quarter of the detector in the shadow. The median ignores the shadow (a minority of each ring) and the spots.
Iterate 3×, excluding the shadow it finds, so the baseline is not pulled down by the shadow.
The rings are drawn about the centre the caller names with `BeamCenter()`, which defaults to the
file's. Tens of pixels of error cost nothing, but a couple of hundred draws the rings across the
image_analysis/beam_stop/ShadowFinder.cpp
+24 -15
View File
@@ -694,14 +694,17 @@ std::vector<uint32_t> ShadowFinder::GetMask(size_t nthreads) const {
|| std::fabs(beam_x - W * 0.5f) > 4.0f * diag || std::fabs(beam_y - H * 0.5f) > 4.0f * diag)
return mask;
// mean projection, usable pixels and radius from the beam centre. The mean is divided by the
// polarization factor, so that what is left varies around a ring only where something is in the
// way; `pol` is kept because the counts the Poisson test is made of are the ones that were
// recorded, not these. The factor is read off the experiment's own geometry - which follows the
// detector's tilt, quarter turns and in-plane rotation - about the centre the rings use.
// mean projection, usable pixels and the ring each pixel is on (ring_position). The mean is
// divided by the polarization factor and by the obliquity of a tilted detector, so that what is
// left varies around a ring only where something is in the way; `pol` keeps the product because
// the counts the Poisson test is made of are the ones that were recorded, not these. Both are
// read off the experiment's own geometry - which follows the detector's tilt, quarter turns and
// in-plane rotation - about the centre the caller named.
// Following Kahn, Fourme, Gadet, Janin, Dumas & Andre (1982) J. Appl. Cryst. 15, 330-337
DiffractionGeometry pol_geometry = geometry;
pol_geometry.BeamX_pxl(beam_x).BeamY_pxl(beam_y);
const std::vector<float> det_matrix = geometry.GetDetectorMatrix().arr();
const float distance_px = geometry.GetDetectorDistance_mm() / geometry.GetPixelSize_mm();
Plane<float> mean(n_pixels);
Plane<float> pol(n_pixels);
@@ -713,12 +716,13 @@ std::vector<uint32_t> ShadowFinder::GetMask(size_t nthreads) const {
for (int y = ylo; y < yhi; y++)
for (int x = 0; x < W; x++) {
const int i = y * W + x;
const float dx = x - beam_x, dy = y - beam_y;
pol[i] = polarization.has_value()
? pol_geometry.CalcAzIntPolarizationCorr(static_cast<float>(x),
static_cast<float>(y),
polarization.value())
: 1.0f;
const RingPosition ring = ring_position(det_matrix.data(), distance_px, x - beam_x, y - beam_y);
pol[i] = ring.obliquity
* (polarization.has_value()
? pol_geometry.CalcAzIntPolarizationCorr(static_cast<float>(x),
static_cast<float>(y),
polarization.value())
: 1.0f);
// pol > 0 also rejects a NaN factor. The correction is exactly zero only where
// the polarization fully suppresses the scattering, so such a pixel carries no
// usable ring information anyway - and divided by zero it would put an inf into
@@ -730,7 +734,7 @@ std::vector<uint32_t> ShadowFinder::GetMask(size_t nthreads) const {
/ valid_count[i] / pol[i]);
valid[i] = 1;
}
radius[i] = static_cast<int>(std::lround(std::sqrt(dx * dx + dy * dy)));
radius[i] = static_cast<int>(std::lround(ring.radius));
local_max = std::max(local_max, radius[i]);
}
// max is associative, so folding the per-worker maxima gives the same answer whatever
@@ -965,8 +969,10 @@ std::vector<uint32_t> ShadowFinder::GetMask(size_t nthreads) const {
for (int i = lo; i < hi; i++) {
if (!valid[i] || region[i])
continue;
const float dx = static_cast<float>(i % W) - beam_x, dy = static_cast<float>(i / W) - beam_y;
const double phi = std::atan2(dy, dx) + std::numbers::pi;
const RingPosition ring = ring_position(det_matrix.data(), distance_px,
static_cast<float>(i % W) - beam_x,
static_cast<float>(i / W) - beam_y);
const double phi = std::atan2(ring.ly, ring.lx) + std::numbers::pi;
const int sector = std::min(HARMONIC_SECTORS - 1, static_cast<int>(phi / (2.0 * std::numbers::pi) * HARMONIC_SECTORS));
values[static_cast<size_t>(radius[i] / HARMONIC_BAND_PX) * HARMONIC_SECTORS + sector].push_back(ratio[i]);
}
@@ -995,7 +1001,10 @@ std::vector<uint32_t> ShadowFinder::GetMask(size_t nthreads) const {
if (!valid[i] || region[i] || ratio[i] >= PENUMBRA_RATIO)
continue;
// cos 2phi and sin 2phi from the offset to the beam.
const float dx = static_cast<float>(i % W) - beam_x, dy = static_cast<float>(i / W) - beam_y;
const RingPosition ring = ring_position(det_matrix.data(), distance_px,
static_cast<float>(i % W) - beam_x,
static_cast<float>(i / W) - beam_y);
const float dx = ring.lx, dy = ring.ly;
const float r2 = std::max(dx * dx + dy * dy, 1e-6f);
const int band = radius[i] / HARMONIC_BAND_PX;
const float model = std::min(1.0f, 1.0f + harm_c[band] * (dx * dx - dy * dy) / r2
image_analysis/beam_stop/ShadowFinder.h
+9 -7
View File
@@ -26,7 +26,7 @@
// and is 1 where the beam stop shadows the detector.
//
// The shadow is a place where the background is missing, so it is found by comparing each
// pixel's mean against the typical background at the same radius - the median over its ring,
// pixel's mean against the typical background at the same 2 theta - the median over its ring,
// taken over the pixels not already known to be shadowed. That comparison holds wherever the
// ring still has unshadowed pixels to measure. Where it does not - a ring lying wholly inside
// the stop - there is nothing to compare against, and such a ring is shadow in its entirety.
@@ -51,9 +51,10 @@ class ShadowFinder {
// 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. Of the
// corrections a ring carries this is the only one that varies along it; solid angle, detector
// and air absorption are all functions of 2 theta alone and the ring's own median absorbs them.
// 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
@@ -126,9 +127,10 @@ public:
ShadowFinder(const DiffractionExperiment &experiment, const PixelMask &mask);
~ShadowFinder();
// The centre the rings are drawn about. It starts as the file's, which is the only one there
// is when the finder is built; a caller that has measured one replaces it before reading the
// mask. A centre far from the truth draws the rings across the background's own radial
// The beam centre of the geometry the rings are drawn in - the point of normal incidence, which
// is where the beam lands only on an untilted detector. It starts as the file's, which is the
// only one there is when the finder is built; a caller that has measured one replaces it before
// reading the mask. A centre far from the truth draws the rings across the background's own radial
// fall-off instead of along it, and the comparison then describes the fall-off rather than the
// hardware. The projection is not centred on anything, so this may be set after the frames.
void BeamCenter(float x, float y);
image_analysis/beam_stop/ShadowFinderInternal.h
+35 -2
View File
@@ -4,14 +4,47 @@
#pragma once
// What ShadowFinder::GetMask shares with its device twin (ShadowMaskGPU): the constants it tests
// against, and the two small fits made over whole rings and sectors, which stay on the host on both
// paths.
// against, where a pixel sits on its ring, and the two small fits made over whole rings and sectors,
// which stay on the host on both paths.
#include <cmath>
#include <cstdint>
#include <vector>
#ifdef __CUDACC__
#define SHADOW_FINDER_HD __host__ __device__ inline
#else
#define SHADOW_FINDER_HD inline
#endif
namespace shadow_finder {
// Where a pixel sits on the ring it is compared along. A ring is a cone of constant 2theta about the
// beam. On an untilted detector that is a circle about the beam centre, but on a tilted one - a
// detector on a swung 2theta arm - it is a circle about no pixel at all, and a circle drawn about the
// centre runs across the background's own fall-off: 5 degrees of tilt at 50 mm put a quarter of a
// detector in the shadow. So a pixel's ring is its 2theta, in pixels: `radius` is distance *
// tan(2theta), the radius it would have on an untilted detector at the same distance, and (lx, ly) is
// its offset from the beam in that frame, which gives its azimuth. A tilt also makes a pixel's solid
// angle vary along the ring, by tens of per cent at that same tilt; `obliquity` is the pixel's solid
// angle over the one the untilted detector has at the same 2theta, (distance / lz)^3, and the
// background is divided by it as by the polarization factor.
// (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.
struct RingPosition {
float lx, ly, radius, obliquity;
};
SHADOW_FINDER_HD RingPosition ring_position(const float *det_matrix, float distance_px, float dx, float dy) {
const float *m = det_matrix;
const float lx = m[0] * dx + m[1] * dy + m[2] * distance_px;
const float ly = m[3] * dx + m[4] * dy + m[5] * distance_px;
const float lz = m[6] * dx + m[7] * dy + m[8] * distance_px;
const float stretch = distance_px / lz;
return {lx, ly, sqrtf(lx * lx + ly * ly) * stretch, stretch * stretch * stretch};
}
// The values of ShadowFinder::SHADOW and ShadowFinder::TRANSMITTING, for the device code, which does
// not include ShadowFinder.h.
inline constexpr uint32_t MASK_SHADOW = 1;
image_analysis/beam_stop/ShadowMaskGPU.cu
+12 -6
View File
@@ -80,8 +80,9 @@ __global__ void setup_kernel(ShadowMaskSetup s, const uint32_t *__restrict__ pix
if (i >= n)
return;
const int x = static_cast<int>(i % s.width), y = static_cast<int>(i / s.width);
const float dx = x - s.beam_x, dy = y - s.beam_y;
const float p = s.has_polarization ? polarization_factor(s, static_cast<float>(x), static_cast<float>(y)) : 1.0f;
const RingPosition ring = ring_position(s.det_matrix, s.distance_mm / s.pixel_size_mm, x - s.beam_x, y - s.beam_y);
const float p = ring.obliquity
* (s.has_polarization ? polarization_factor(s, static_cast<float>(x), static_cast<float>(y)) : 1.0f);
pol[i] = p;
float mean = 0.0f;
char v = 0;
@@ -92,7 +93,7 @@ __global__ void setup_kernel(ShadowMaskSetup s, const uint32_t *__restrict__ pix
valid[i] = v;
num[i] = v ? mean : 0.0;
den[i] = v ? 1 : 0;
const int r = static_cast<int>(lroundf(sqrtf(dx * dx + dy * dy)));
const int r = static_cast<int>(lroundf(ring.radius));
radius[i] = r;
atomicMax(max_radius, r);
}
@@ -360,8 +361,10 @@ __global__ void sector_key_kernel(ShadowMaskSetup s, const char *__restrict__ va
key[i] = UINT64_MAX;
return;
}
const float dx = static_cast<float>(i % s.width) - s.beam_x, dy = static_cast<float>(i / s.width) - s.beam_y;
const double phi = atan2f(dy, dx) + PI;
const RingPosition ring = ring_position(s.det_matrix, s.distance_mm / s.pixel_size_mm,
static_cast<float>(i % s.width) - s.beam_x,
static_cast<float>(i / s.width) - s.beam_y);
const double phi = atan2f(ring.ly, ring.lx) + PI;
const int sector = min(HARMONIC_SECTORS - 1, static_cast<int>(phi / (2.0 * PI) * HARMONIC_SECTORS));
const uint64_t k = static_cast<uint64_t>(radius[i] / HARMONIC_BAND_PX) * HARMONIC_SECTORS + sector;
key[i] = (k << 32) | float_key(ratio[i]);
@@ -388,7 +391,10 @@ __global__ void dim_kernel(ShadowMaskSetup s, uint32_t frames, const char *__res
char d = 0;
if (valid[i] && !region[i] && ratio[i] < PENUMBRA_RATIO) {
// cos 2phi and sin 2phi from the offset to the beam.
const float dx = static_cast<float>(i % s.width) - s.beam_x, dy = static_cast<float>(i / s.width) - s.beam_y;
const RingPosition ring = ring_position(s.det_matrix, s.distance_mm / s.pixel_size_mm,
static_cast<float>(i % s.width) - s.beam_x,
static_cast<float>(i / s.width) - s.beam_y);
const float dx = ring.lx, dy = ring.ly;
const float r2 = fmaxf(dx * dx + dy * dy, 1e-6f);
const int band = radius[i] / HARMONIC_BAND_PX;
const float model = fminf(1.0f, 1.0f + harm_c[band] * (dx * dx - dy * dy) / r2
image_analysis/beam_stop/ShadowMaskGPU.h
+3 -3
View File
@@ -10,9 +10,9 @@
#include <cuda_runtime.h>
// What the mask is drawn about: the detector, the centre the rings are drawn about, and the geometry
// the polarization factor is 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;
tests/ShadowFinderTest.cpp
+36
View File
@@ -336,6 +336,42 @@ TEST_CASE("ShadowFinder_TheRingsFollowTheCentreTheCallerNames", "[ShadowFinder]"
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
// normal incidence) and the rings of constant 2 theta are no circle about it, nor about any pixel;
// the solid angle of a pixel, which falls away from the normal, also varies along each of them. A
// background that depends on 2 theta alone, seen through that solid angle, has nothing in the
// beam - and circles drawn about the centre cut across its fall-off and read a large part of the
// detector as shadow.
TEST_CASE("ShadowFinder_ATiltedDetectorIsNotAShadow", "[ShadowFinder]") {
// Short enough, and tilted enough, that the beam lands 47 px from the centre.
constexpr float DISTANCE_MM = 20.0f;
constexpr float TILT_RAD = 0.174533f; // 10 degrees
DiffractionExperiment x = TestExperiment();
x.DetectorDistance_mm(DISTANCE_MM).PoniRot1_rad(TILT_RAD);
const DiffractionGeometry geom = x.GetDiffractionGeometry();
// Below MIN_REFLECTION, so no pixel is exempt as a recorded reflection.
std::vector<int32_t> frame(static_cast<size_t>(W) * H);
for (int y = 0; y < H; y++)
for (int xi = 0; xi < W; xi++) {
const float two_theta = geom.TwoTheta_rad(static_cast<float>(xi), static_cast<float>(y));
const float solid_angle = geom.CalcAzIntSolidAngleCorr(static_cast<float>(xi), static_cast<float>(y));
frame[I(xi, y)] = static_cast<int32_t>(std::lround(24.0f * solid_angle * std::exp(-two_theta / 0.15f)));
}
const PixelMask pixel_mask(x);
ShadowFinder finder(x, pixel_mask);
std::vector<uint8_t> buffer;
for (int f = 0; f < NFRAMES; f++) {
DataMessage msg{};
msg.image = CompressedImage(frame, W, H);
finder.AddImage(msg, buffer);
}
const auto mask = finder.GetMask();
CHECK(std::count_if(mask.begin(), mask.end(), [](uint32_t v) { return v != 0; }) == 0);
}
// The background is not flat around a ring: a polarized source suppresses it in its own plane, by
// a factor that reaches four at the 2 theta a short detector distance puts in a corner. That is
// several times the dip this class looks for, so with the modulation left in, half of every outer