mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 78 Wiki Activity
Files
b3c9987ed7396d2adf65abdbdbfbf1f52ed22fef
Jungfraujoch/image_analysis/beam_stop/ShadowFinder.cpp
T
leonarski_fandClaude Opus 5 a7ed4e2d34 beam stop: a bright ring is not a beam stop 
Inside the stop a whole ring is blocked, so its own median is blocked too and the per-pixel
comparison has nothing to work with. The walk that covers that case declared a ring to be inside the
stop when its background fell below a third of the LARGEST background of any ring further out - and
on a sample whose background peaks in a strong ring away from the beam, the ordinary background
inside that ring is legitimately below a third of the peak. The walk then runs out to the ring and
returns a filled disk of good detector: on one corpus dataset 16 % of the area, with diffraction
rings plainly visible inside the disk it masked, and on another 2.5 %.

The defect predates the branch; what changed is that the flood from the beam centre used to discard
the disk whenever its seeds landed on invalid pixels, which is how the first of those two datasets
came back with an empty mask instead of a wrong one. Dropping that anchor was right, and it made
this visible.

A ring is now compared against what this detector's background typically is - the median over the
rings the walk is willing to judge - which is robust to a bright ring and to a corner ring of a
handful of pixels alike, and is less code. Evaluated over 151 corpus datasets: the number returning
a disk larger than 0.2 % of the detector falls from 12 to 3, the two pathological cases collapse
(radius 974 -> 68 px and 188 -> 12 px), and genuine stops move by a few pixels at most (152 -> 132,
122 -> 120, 70 -> 64). The new test fails on the old walk and passes on this one.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-09-13 07:36:57 +02:00

791 lines
36 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "ShadowFinder.h"
#include <algorithm>
#include <cmath>
#include <limits>
#include <atomic>
#include <future>
#include <thread>
#include <queue>
#include <type_traits>
#include <spdlog/spdlog.h>
#include "../../common/CUDAWrapper.h"
#include "../../common/ParallelFor.h"
#include "../../common/JFJochException.h"
// A pixel is shadow when its background is below this fraction of the background it is
// compared against.
constexpr float SHADOW_RATIO = 0.50f;
// The boundary grows outward into partially shadowed pixels down to this fraction, but no
// further than PENUMBRA_MAX_PX from the core. A pin or a loop casts a wide half-shadow, so the
// reach is a good deal more than the beam stop's own edge needs.
constexpr float PENUMBRA_RATIO = 0.75f;
constexpr int PENUMBRA_MAX_PX = 30;
// Bridge module gaps and small breaks that the holder arm crosses.
constexpr int BRIDGE_PX = 6;
// A pixel whose maximum reaches this recorded a real reflection and is never masked - a
// beam stop cannot block a reflection that was measured.
constexpr int64_t MIN_REFLECTION = 25;
// How far below the background it is compared against a pixel must sit before the dip is
// believed, in standard deviations of the counts that back it. The counts are photons, so their
// scatter is Poisson and the deficit is measured against it rather than against a fixed number:
// on a well-exposed sweep a third of the background missing is overwhelming, and on a handful of
// low-background frames the same third is noise. Without this a six-frame pre-scan of a
// low-background sweep masks three quarters of the detector.
constexpr double MIN_DEFICIT_SIGMA = 6.0;
// Smallest region the per-pixel test may return. A shadow is cast by something physical and is
// correspondingly large; an isolated patch this small is the background wandering, not hardware.
// This is what keeps the test specific now that a shadow no longer has to touch the direct beam.
constexpr int MIN_SHADOW_PIXELS = 2000;
// Side of the box the background is pooled over before testing. Its area is how many pixels back
// a ring's countability test, which decides where an azimuthal comparison is possible at all.
constexpr int POOL_PX = 5;
constexpr double MEAN_POOLED_PIXELS = POOL_PX * POOL_PX;
// A ring with fewer valid pixels than this says nothing about whether it was counted.
constexpr int MIN_RING_PIXELS = 32;
// A ring lies wholly inside the stop when its background is below this fraction of the background
// further out. This asks about a whole ring rather than about a pixel, so it keeps a threshold of
// its own and does not follow SHADOW_RATIO.
constexpr float BLOCKED_RING_RATIO = 0.35f;
// Binary-image helpers on a width*height frame stored row-major as char (0/1). All run once,
// at GetMask() time; the BFS forms keep them O(pixels) rather than O(pixels * radius).
namespace {
// 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) {
if (expected <= 0.0 || observed >= expected)
return 0.0;
const double ll = 2.0 * (expected - observed
+ (observed > 0.0 ? observed * std::log(observed / expected) : 0.0));
return ll > 0.0 ? std::sqrt(ll) : 0.0;
}
// 8-connected dilation by `r` pixels, i.e. every pixel within Chebyshev distance r of a set one.
//
// This was a multi-source BFS, which is what the distance is defined by - but on a full rectangle
// with no obstacles the 8-connected graph distance IS the Chebyshev distance (a path stepping
// towards the target never has to leave the frame), so the result is a dilation by the (2r+1) square
// clipped to the frame. A square is the product of a horizontal and a vertical segment and max is
// associative, so it separates into one pass along x and one along y - O(1) per pixel whatever r is,
// no queue, and no 4-bytes-per-pixel distance array. The values are 0/1, so the running "is any set"
// is just a count of set pixels in the window.
std::vector<char> dilate(const std::vector<char> &in, int W, int H, int r, size_t nthreads) {
if (r <= 0)
return in;
std::vector<char> tmp(in.size()), out(in.size());
ParallelChunks(H, nthreads, [&](int ylo, int yhi) {
for (int y = ylo; y < yhi; y++) {
const char *src = in.data() + static_cast<size_t>(y) * W;
char *dst = tmp.data() + static_cast<size_t>(y) * W;
int count = 0;
for (int x = 0; x <= std::min(r, W - 1); x++)
count += src[x];
for (int x = 0; x < W; x++) {
dst[x] = count > 0;
if (x + r + 1 < W) count += src[x + r + 1];
if (x - r >= 0) count -= src[x - r];
}
}
});
// The vertical pass walks a column, which strides by a whole row. Take a strip of columns at a
// time so each row it touches is read contiguously instead of one byte per cache line.
constexpr int STRIP = 64;
ParallelChunks((W + STRIP - 1) / STRIP, nthreads, [&](int slo, int shi) {
std::vector<int> count(STRIP);
for (int st = slo; st < shi; st++) {
const int x0 = st * STRIP, xn = std::min(STRIP, W - x0);
std::fill(count.begin(), count.begin() + xn, 0);
for (int y = 0; y <= std::min(r, H - 1); y++)
for (int i = 0; i < xn; i++)
count[i] += tmp[static_cast<size_t>(y) * W + x0 + i];
for (int y = 0; y < H; y++) {
for (int i = 0; i < xn; i++)
out[static_cast<size_t>(y) * W + x0 + i] = count[i] > 0;
if (y + r + 1 < H)
for (int i = 0; i < xn; i++)
count[i] += tmp[static_cast<size_t>(y + r + 1) * W + x0 + i];
if (y - r >= 0)
for (int i = 0; i < xn; i++)
count[i] -= tmp[static_cast<size_t>(y - r) * W + x0 + i];
}
}
});
return out;
}
// Erosion by `r` = dilation of the complement. The dilation is clipped to the frame and cannot seed
// outside it, so outside the frame contributes nothing - a pixel within r of the edge is eroded only
// by what the frame actually holds.
std::vector<char> erode(const std::vector<char> &in, int W, int H, int r, size_t nthreads) {
std::vector<char> comp(in.size());
ParallelChunks(static_cast<int>(in.size()), nthreads, [&](int lo, int hi) {
for (int i = lo; i < hi; i++) comp[i] = !in[i];
});
const auto grown = dilate(comp, W, H, r, nthreads);
std::vector<char> out(in.size());
ParallelChunks(static_cast<int>(out.size()), nthreads, [&](int lo, int hi) {
for (int i = lo; i < hi; i++) out[i] = !grown[i];
});
return out;
}
// Fill holes: background not reachable from the image border becomes region.
//
// The flood is run over the bounding box of `region` grown by one, not the whole detector. Outside
// that box every pixel is background and the box's surrounding ring is background too, so the whole
// outside is one border-connected component: a background pixel inside the box is border-connected
// exactly when it reaches the ring. The beam stop occupies a small part of a detector, so this is
// the same answer over a fraction of the pixels.
std::vector<char> fill_holes(const std::vector<char> &region, int W, int H) {
int x0 = W, x1 = -1, y0 = H, y1 = -1;
for (int y = 0; y < H; y++)
for (int x = 0; x < W; x++)
if (region[static_cast<size_t>(y) * W + x]) {
x0 = std::min(x0, x); x1 = std::max(x1, x);
y0 = std::min(y0, y); y1 = std::max(y1, y);
}
if (x1 < 0)
return region; // nothing to enclose
x0 = std::max(0, x0 - 1); x1 = std::min(W - 1, x1 + 1);
y0 = std::max(0, y0 - 1); y1 = std::min(H - 1, y1 + 1);
const int BW = x1 - x0 + 1, BH = y1 - y0 + 1;
std::vector<char> bg_visited(static_cast<size_t>(BW) * BH, 0);
std::queue<int> q; // indices into the box
auto push = [&](int bx, int by) {
const int j = by * BW + bx;
if (!region[static_cast<size_t>(by + y0) * W + bx + x0] && !bg_visited[j]) {
bg_visited[j] = 1; q.push(j);
}
};
for (int bx = 0; bx < BW; bx++) { push(bx, 0); push(bx, BH - 1); }
for (int by = 0; by < BH; by++) { push(0, by); push(BW - 1, by); }
while (!q.empty()) {
const int i = q.front(); q.pop();
const int by = i / BW, bx = i % BW;
for (int dy = -1; dy <= 1; dy++)
for (int dx = -1; dx <= 1; dx++) {
const int yy = by + dy, xx = bx + dx;
if (yy < 0 || yy >= BH || xx < 0 || xx >= BW)
continue;
const int j = yy * BW + xx;
if (!region[static_cast<size_t>(yy + y0) * W + xx + x0] && !bg_visited[j]) {
bg_visited[j] = 1; q.push(j);
}
}
}
std::vector<char> out = region;
for (int by = 0; by < BH; by++)
for (int bx = 0; bx < BW; bx++) {
const size_t i = static_cast<size_t>(by + y0) * W + bx + x0;
if (!region[i] && !bg_visited[by * BW + bx])
out[i] = 1;
}
return out;
}
// Sum of `in` over the k x k box centred on each pixel, zero outside the frame.
//
// Each row's and each column's running sum keeps exactly the terms it had, in the order it had them,
// so the floating-point form rounds identically - only the traversal changes. Walking one column at
// a time strides a whole row per step and misses on every access, so the vertical pass takes a strip
// of columns together and reads each row it touches contiguously. Rows, and strips, are independent.
template <typename T>
std::vector<T> box_sum(const std::vector<T> &in, int W, int H, int k, size_t nthreads) {
const int half = k / 2;
std::vector<T> row(in.size()), out(in.size());
ParallelChunks(H, nthreads, [&](int ylo, int yhi) {
for (int y = ylo; y < yhi; y++) {
const T *src = in.data() + static_cast<size_t>(y) * W;
T *dst = row.data() + static_cast<size_t>(y) * W;
T s = 0;
for (int x = 0; x <= std::min(half, W - 1); x++)
s += src[x];
for (int x = 0; x < W; x++) {
dst[x] = s;
if (x + half + 1 < W) s += src[x + half + 1];
if (x - half >= 0) s -= src[x - half];
}
}
});
constexpr int STRIP = 64;
ParallelChunks((W + STRIP - 1) / STRIP, nthreads, [&](int slo, int shi) {
std::vector<T> s(STRIP);
for (int st = slo; st < shi; st++) {
const int x0 = st * STRIP, xn = std::min(STRIP, W - x0);
std::fill(s.begin(), s.begin() + xn, T{0});
for (int y = 0; y <= std::min(half, H - 1); y++)
for (int i = 0; i < xn; i++)
s[i] += row[static_cast<size_t>(y) * W + x0 + i];
for (int y = 0; y < H; y++) {
for (int i = 0; i < xn; i++)
out[static_cast<size_t>(y) * W + x0 + i] = s[i];
if (y + half + 1 < H)
for (int i = 0; i < xn; i++)
s[i] += row[static_cast<size_t>(y + half + 1) * W + x0 + i];
if (y - half >= 0)
for (int i = 0; i < xn; i++)
s[i] -= row[static_cast<size_t>(y - half) * W + x0 + i];
}
}
});
return out;
}
// The values of each ring, laid end to end, with the ring's slice given by offset[r]..offset[r+1].
// radius and pooled do not change over the three baseline iterations, so this is built once and
// every iteration is an order statistic of the same, already-sorted, ring.
struct RingValues {
std::vector<float> values;
std::vector<int> offset;
};
RingValues bin_by_ring(const std::vector<float> &values, const std::vector<char> &valid,
const std::vector<int> &radius, int max_radius, size_t nthreads) {
RingValues rv;
rv.offset.assign(max_radius + 2, 0);
for (size_t i = 0; i < values.size(); i++)
if (valid[i])
rv.offset[radius[i] + 1]++;
for (int r = 0; r <= max_radius; r++)
rv.offset[r + 1] += rv.offset[r];
rv.values.resize(rv.offset[max_radius + 1]);
std::vector<int> cursor(rv.offset.begin(), rv.offset.end() - 1);
for (size_t i = 0; i < values.size(); i++)
if (valid[i])
rv.values[cursor[radius[i]]++] = values[i];
// Sorted once; the three iterations then only pick a rank and count a prefix.
ParallelFor(max_radius + 1, nthreads, [&](int r) {
std::sort(rv.values.begin() + rv.offset[r], rv.values.begin() + rv.offset[r + 1]);
});
return rv;
}
} // namespace
ShadowFinder::ShadowFinder(const DiffractionExperiment &experiment, const PixelMask &mask)
: width(static_cast<int>(experiment.GetXPixelsNumConv())),
height(static_cast<int>(experiment.GetYPixelsNumConv())),
beam_x(experiment.GetBeamX_pxl()),
beam_y(experiment.GetBeamY_pxl()),
geometry(experiment.GetDiffractionGeometry()),
polarization(experiment.GetPolarizationFactor()),
pixel_mask(mask.GetMask(experiment)) {
if (pixel_mask.size() != static_cast<size_t>(width) * height)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"ShadowFinder: pixel mask does not match the detector");
SetShardCount(1);
#ifdef JFJOCH_USE_CUDA
if (get_gpu_count() > 0) {
const size_t npixels = static_cast<size_t>(width) * height;
gpu_pending = std::async(std::launch::async, [npixels] {
return std::make_unique<ShadowAccumulatorGPU>(npixels);
});
}
#endif
}
void ShadowFinder::BeamCenter(float x, float y) {
std::unique_lock ul(m);
beam_x = x;
beam_y = y;
}
// A shard's accumulators are allocated when a frame is first added to it, not here: with a GPU they
// are never used at all, and on a 16 Mpx detector eight of them are 2.9 GB to allocate and clear -
// which measured 0.8 s of the pre-scan, all of it wasted.
void ShadowFinder::SetShardCount(size_t n) {
shards.clear();
shards.resize(std::max<size_t>(1, n));
}
ShadowFinder::Projection ShadowFinder::Reduce() const {
#ifdef JFJOCH_USE_CUDA
// The device holds its own projection. Bring it back and let it take part in the fold below as
// one more shard; when every frame went to the GPU it is the whole answer.
Projection device;
if (Gpu() && gpu->GetFrameCount() > 0) {
gpu->Download(device.max_value, device.sum_value, device.valid_count);
device.frames = gpu->GetFrameCount();
bool host_empty = true;
for (const auto &p : shards)
host_empty = host_empty && (p.frames == 0);
if (host_empty)
return device;
}
#endif
// Only when that shard actually holds something: its accumulators are allocated on first use, so
// an unused shard is empty rather than zeroed, and returning it would hand the callers below a
// projection they index by pixel.
if (shards.size() == 1 && shards[0].frames > 0
#ifdef JFJOCH_USE_CUDA
&& !(gpu && gpu->GetFrameCount() > 0)
#endif
)
return shards[0];
Projection out;
const size_t npixels = static_cast<size_t>(width) * height;
out.max_value.assign(npixels, 0);
out.sum_value.assign(npixels, 0);
out.valid_count.assign(npixels, 0);
for (const auto &p : shards)
out.frames += p.frames;
#ifdef JFJOCH_USE_CUDA
out.frames += device.frames;
#endif
// Each worker owns a slice of the pixels and folds every shard into it. The sums and counts are
// integers and a pixel is touched by one worker only, so the result is the same as folding them
// one shard at a time on one thread - this is several hundred megabytes per shard and is limited
// by memory rather than by arithmetic.
const size_t nthreads = std::max<size_t>(1, std::min<size_t>(std::thread::hardware_concurrency(),
shards.size() * 2));
const size_t chunk = (npixels + nthreads - 1) / nthreads;
std::vector<std::future<void>> futures;
futures.reserve(nthreads);
for (size_t t = 0; t < nthreads; t++) {
const size_t lo = t * chunk, hi = std::min(npixels, lo + chunk);
if (lo >= hi) break;
futures.emplace_back(std::async(std::launch::async, [&, lo, hi] {
#ifdef JFJOCH_USE_CUDA
const Projection *extra[1] = {&device};
for (const auto *pp : extra) {
const auto &p = *pp;
if (p.frames == 0) continue;
for (size_t i = lo; i < hi; i++) {
if (p.valid_count[i] == 0)
continue;
if (out.valid_count[i] == 0 || p.max_value[i] > out.max_value[i])
out.max_value[i] = p.max_value[i];
out.sum_value[i] += p.sum_value[i];
out.valid_count[i] += p.valid_count[i];
}
}
#endif
for (const auto &p : shards) {
if (p.frames == 0) continue; // never used, and its accumulators were never allocated
for (size_t i = lo; i < hi; i++) {
if (p.valid_count[i] == 0)
continue;
if (out.valid_count[i] == 0 || p.max_value[i] > out.max_value[i])
out.max_value[i] = p.max_value[i];
out.sum_value[i] += p.sum_value[i];
out.valid_count[i] += p.valid_count[i];
}
}
}));
}
for (auto &f : futures) f.get();
return out;
}
template<class T>
void ShadowFinder::Add(const T *ptr, Projection &p) {
// The pixel type's sentinel extreme marks "no data" (module gap / masked): the
// preprocessor/writer stores INT*_MIN for signed and UINT*_MAX for unsigned. For signed
// types the opposite extreme is a genuine saturated value and is kept, so a saturated
// reflection still registers as bright.
T masked;
if constexpr (std::is_signed_v<T>)
masked = std::numeric_limits<T>::min();
else
masked = std::numeric_limits<T>::max();
for (size_t i = 0; i < p.max_value.size(); i++) {
const T v = ptr[i];
if (v == masked)
continue;
const int64_t vi = static_cast<int64_t>(v);
if (p.valid_count[i] == 0 || vi > p.max_value[i])
p.max_value[i] = vi;
p.sum_value[i] += vi;
p.valid_count[i]++;
}
p.frames++;
}
void ShadowFinder::AddImage(const DataMessage &data, std::vector<uint8_t> &buffer, size_t shard) {
if (static_cast<size_t>(data.image.GetWidth()) * data.image.GetHeight()
!= static_cast<size_t>(width) * height)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"ShadowFinder: image size does not match the detector");
if (shard >= shards.size())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"ShadowFinder: shard out of range");
#ifdef JFJOCH_USE_CUDA
// One device, so the frames queue here - but each is only a chunk upload plus two kernels, and
// nothing was decompressed on the host to get this far.
if (ShadowAccumulatorGPU *acc = ShadowAccumulatorGPU::Supports(data.image) ? Gpu() : nullptr) {
std::unique_lock ul(gpu_mutex);
try {
acc->Add(data.image);
return;
} catch (const std::exception &e) {
spdlog::warn("Beam stop: GPU accumulate failed ({}), falling back to the host", e.what());
}
}
#endif
Projection &p = shards[shard];
if (p.max_value.empty()) {
const size_t npixels = static_cast<size_t>(width) * height;
p.max_value.assign(npixels, 0);
p.sum_value.assign(npixels, 0);
p.valid_count.assign(npixels, 0);
}
const auto ptr = data.image.GetUncompressedPtr(buffer);
switch (data.image.GetMode()) {
case CompressedImageMode::Int8: Add(reinterpret_cast<const int8_t *>(ptr), p); break;
case CompressedImageMode::Uint8: Add(reinterpret_cast<const uint8_t *>(ptr), p); break;
case CompressedImageMode::Int16: Add(reinterpret_cast<const int16_t *>(ptr), p); break;
case CompressedImageMode::Uint16: Add(reinterpret_cast<const uint16_t *>(ptr), p); break;
case CompressedImageMode::Int32: Add(reinterpret_cast<const int32_t *>(ptr), p); break;
case CompressedImageMode::Uint32: Add(reinterpret_cast<const uint32_t *>(ptr), p); break;
default:
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"ShadowFinder: unsupported image mode");
}
}
#ifdef JFJOCH_USE_CUDA
ShadowAccumulatorGPU *ShadowFinder::Gpu() const {
std::unique_lock ul(gpu_mutex);
if (gpu_pending.valid()) {
try {
gpu = gpu_pending.get();
} catch (const std::exception &e) {
// A GPU that cannot hold the projection is not a reason to fail: the host path gives
// the same answer, only slower.
spdlog::warn("Beam stop: GPU projection unavailable ({}), accumulating on the host",
e.what());
gpu.reset();
}
}
return gpu.get();
}
#endif
uint32_t ShadowFinder::GetFrameCount() const {
std::unique_lock ul(m);
uint32_t frames = 0;
for (const auto &p : shards) frames += p.frames;
#ifdef JFJOCH_USE_CUDA
if (gpu) frames += gpu->GetFrameCount();
#endif
return frames;
}
std::vector<float> ShadowFinder::GetMeanProjection() const {
std::unique_lock ul(m);
const Projection p = Reduce();
const auto &sum_value = p.sum_value;
const auto &valid_count = p.valid_count;
std::vector<float> mean(static_cast<size_t>(width) * height, NAN);
for (size_t i = 0; i < mean.size(); i++)
if (valid_count[i] > 0 && pixel_mask[i] == 0)
mean[i] = static_cast<float>(static_cast<double>(sum_value[i]) / valid_count[i]);
return mean;
}
std::vector<uint32_t> ShadowFinder::GetMask(size_t nthreads) const {
std::unique_lock ul(m);
if (nthreads == 0)
nthreads = std::max(1u, std::thread::hardware_concurrency());
const Projection p = Reduce();
const auto &max_value = p.max_value;
const auto &sum_value = p.sum_value;
const auto &valid_count = p.valid_count;
const uint32_t frames = p.frames;
const int W = width, H = height;
const int n_pixels = W * H;
std::vector<uint32_t> mask(n_pixels, 0);
if (frames == 0)
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.
// 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);
std::vector<float> mean(n_pixels, 0.0f);
std::vector<float> pol(n_pixels, 1.0f);
std::vector<char> valid(n_pixels, 0);
std::vector<int> radius(n_pixels, 0);
std::atomic<int> 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 float dx = x - beam_x, dy = y - beam_y;
if (polarization.has_value())
pol[i] = pol_geometry.CalcAzIntPolarizationCorr(static_cast<float>(x),
static_cast<float>(y),
polarization.value());
if (valid_count[i] > 0 && pixel_mask[i] == 0) {
mean[i] = static_cast<float>(static_cast<double>(sum_value[i])
/ valid_count[i] / pol[i]);
valid[i] = 1;
}
radius[i] = static_cast<int>(std::lround(std::sqrt(dx * dx + dy * dy)));
local_max = std::max(local_max, radius[i]);
}
// max is associative, so folding the per-worker maxima gives the same answer whatever
// order they finish in.
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();
// Pool the background over a small box before testing it. A background of a fraction of
// a count per pixel per frame gives no single pixel enough counts to tell a shadow from
// a Poisson hole; the stop and its arm are wider than the box, so pooling costs no
// resolution that matters and multiplies the statistics by the pixels in the box.
// The count is a count: at most 25 pixels, so an integer box sum is exact and costs half the
// memory of the floating-point one it replaces. The background itself stays in double - its
// running sum adds and subtracts across a whole row, and in float the rounding of the two would
// not cancel, which moves pixels across the shadow threshold below.
std::vector<double> num(n_pixels);
std::vector<int32_t> den(n_pixels);
ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
for (int i = lo; i < hi; i++) {
num[i] = valid[i] ? mean[i] : 0.0;
den[i] = valid[i] ? 1 : 0;
}
});
const auto pooled_sum = box_sum(num, W, H, POOL_PX, nthreads);
const auto pooled_count = box_sum(den, W, H, POOL_PX, nthreads);
std::vector<float> pooled(n_pixels, 0.0f);
ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
for (int i = lo; i < hi; i++)
if (pooled_count[i] > 0)
pooled[i] = static_cast<float>(pooled_sum[i] / pooled_count[i]);
});
// Azimuthal comparison: the median of the ring, iterated so the shadow stays out of the
// baseline it is measured against. Comparing a pixel's background against the background at its
// own resolution is how DEFPIX recognises a shaded detector region.
// Following Kabsch (2010) Acta Cryst. D66, 125-132
//
// The iteration only ever excludes pixels whose pooled background is below a cut, and dividing by
// a positive baseline is monotone - so the excluded pixels of a ring are exactly the lowest ones,
// and the pixels the next median is taken over are exactly the rest. Each iteration's median is
// therefore an order statistic of the ring's values, which do not change: bin and sort the rings
// once, then each iteration picks a rank and counts a prefix. That replaces nine full-image
// passes (use / ratio / excluded, three times) with one, and three re-binnings with none.
const RingValues rings = bin_by_ring(pooled, valid, radius, max_radius, nthreads);
std::vector<float> baseline(max_radius + 1, 0.0f);
{
std::vector<int> 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;
});
}
std::vector<float> ratio(n_pixels, 1.0f);
ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
for (int i = lo; i < hi; i++)
if (valid[i])
ratio[i] = pooled[i] / std::max(baseline[radius[i]], 1e-6f);
});
// A ring whose background was never counted carries no information to test a pixel against.
// Walking outward, every ring before the first countable one lies wholly inside the stop - a
// ring fully within the disk has no unshadowed pixel for the median to find, which is exactly
// where an azimuthal comparison must fail. Those rings are shadow in their entirety.
// Innermost rings hold only a handful of pixels, too few to judge, so they are stepped over
// rather than allowed to end the walk.
std::vector<int> ring_pixels(max_radius + 1, 0);
for (int i = 0; i < n_pixels; i++)
if (valid[i])
ring_pixels[radius[i]]++;
// A ring lies inside the stop when its background is a fraction of what this detector's
// background typically is. Counting statistics cannot decide this: on a bright dataset the
// shadow is still well counted. The comparison used to be against the LARGEST background of any
// ring further out, and that reads a sample whose background peaks in a strong ring away from
// the beam - a powder standard, a strong solvent ring - as a beam stop the size of that ring:
// the ordinary background inside it is legitimately below a third of the peak. On one corpus
// dataset it declared 16 % of the detector to be stop, with diffraction rings visible inside the
// disk it masked. The median over the rings this walk is willing to judge is what "typically"
// means, and it is robust from both sides - to a bright ring, and to a corner ring of a handful
// of pixels whose median is one pixel's mean.
std::vector<float> judgeable;
for (int rad = 0; rad <= max_radius; rad++)
if (ring_pixels[rad] >= MIN_RING_PIXELS)
judgeable.push_back(baseline[rad]);
float typical_background = 0.0f;
if (!judgeable.empty()) {
const auto middle = judgeable.begin() + judgeable.size() / 2;
std::nth_element(judgeable.begin(), middle, judgeable.end());
typical_background = *middle;
}
int blocked_out_to = -1;
for (int rad = 0; rad <= max_radius; rad++) {
if (ring_pixels[rad] < MIN_RING_PIXELS)
continue;
if (baseline[rad] >= BLOCKED_RING_RATIO * typical_background)
break;
blocked_out_to = rad;
}
// The counts a pixel's pooled background is made of, and the counts the ring says it should
// have had. The test is on the deficit between them, in units of its own Poisson scatter.
std::vector<float> deficit(n_pixels, 0.0f);
std::vector<char> low(n_pixels, 0);
ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
for (int i = lo; i < hi; i++) {
if (!valid[i])
continue;
// The deficit is significant or not in the counts that were recorded, so the pooled
// 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<double>(frames) * pooled_count[i] * pol[i];
deficit[i] = static_cast<float>(poisson_deficit_sigma(pooled[i] * counted,
baseline[radius[i]] * counted));
low[i] = ratio[i] < SHADOW_RATIO && deficit[i] > MIN_DEFICIT_SIGMA;
}
});
// The shadow is a low region large enough to have been cast by something. Nothing anchors it to
// the beam centre: hardware that shadows the detector need not touch the direct beam, and a pin
// or a loop typically does not - it begins some way out in radius, with lit detector between it
// and the stop. What keeps the test specific instead is size, since the background wanders by a
// pixel or two at a time and hardware does not.
const std::vector<char> bridged = dilate(low, W, H, BRIDGE_PX, nthreads);
std::vector<char> region(n_pixels, 0);
{
std::vector<char> seen(n_pixels, 0);
std::vector<int> component;
std::queue<int> q;
for (int start = 0; start < n_pixels; start++) {
if (!bridged[start] || seen[start])
continue;
component.clear();
int n_low = 0;
seen[start] = 1;
q.push(start);
while (!q.empty()) {
const int i = q.front(); q.pop();
component.push_back(i);
n_low += low[i];
const int y = i / W, x = i % W;
for (int dy = -1; dy <= 1; dy++)
for (int dx = -1; dx <= 1; dx++) {
const int yy = y + dy, xx = x + dx;
if (yy < 0 || yy >= H || xx < 0 || xx >= W)
continue;
const int j = yy * W + xx;
if (bridged[j] && !seen[j]) { seen[j] = 1; q.push(j); }
}
}
if (n_low >= MIN_SHADOW_PIXELS)
for (const int i : component)
region[i] = low[i];
}
}
// The rings that lie wholly inside the stop are decided by the ring walk above rather than by
// the per-pixel test, so they join the region after it and are not asked to be large.
for (int i = 0; i < n_pixels; i++)
if (valid[i] && radius[i] <= blocked_out_to)
region[i] = 1;
// Recorded reflections. A small cluster is required so a single-frame zinger does not count.
std::vector<char> lit(n_pixels, 0);
ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
for (int i = lo; i < hi; i++)
lit[i] = (valid_count[i] > 0) && (max_value[i] >= MIN_REFLECTION);
});
std::vector<char> reflection(n_pixels, 0);
ParallelChunks(H, nthreads, [&](int ylo, int yhi) {
for (int y = ylo; y < yhi; y++)
for (int x = 0; x < W; x++) {
const int i = y * W + x;
if (!lit[i]) continue;
int neighbours = 0;
for (int dy = -1; dy <= 1; dy++)
for (int dx = -1; dx <= 1; dx++) {
const int yy = y + dy, xx = x + dx;
if ((dx || dy) && yy >= 0 && yy < H && xx >= 0 && xx < W && lit[yy * W + xx])
neighbours++;
}
reflection[i] = (neighbours >= 2);
}
});
// Grow the soft boundary, round it and fill the disk interior.
const std::vector<char> penumbra = dilate(region, W, H, PENUMBRA_MAX_PX, nthreads);
for (int i = 0; i < n_pixels; i++)
if (penumbra[i] && valid[i] && ratio[i] < PENUMBRA_RATIO && deficit[i] > MIN_DEFICIT_SIGMA)
region[i] = 1;
region = erode(dilate(region, W, H, 2, nthreads), W, H, 2, nthreads);
region = fill_holes(region, W, H);
// Expose recorded reflections - done last, with no fill afterwards, so a spot the shadow
// still covered is given back rather than re-enclosed.
const std::vector<char> reflection_grown = dilate(reflection, W, H, 1, nthreads);
for (int i = 0; i < n_pixels; i++)
if (reflection_grown[i])
region[i] = 0;
for (int i = 0; i < n_pixels; i++)
mask[i] = region[i] ? 1 : 0;
return mask;
}
Reference in New Issue View Git Blame Copy Permalink