leonarski_fandClaude Opus 5.5
b7e4ea23a3
A pixel stuck at the saturation limit was never masked: the defective-pixel scan counted only the error value as "invalid on most frames", and a saturated pixel merely as lit, which the hot rule does not judge on stills. A Rayonix stills set has 3,627 pixels at 65535 on every frame, at high resolution; each flagged every reflection whose disk it fell in as overloaded, and the spot finder found them on every frame. The saturation limit now counts with the error value (CPU and CUDA), so such a pixel is masked like one holding the error value. On that set the indexed fraction rises from 10% to 33% (and 7% to 30% on the second lattice of the same run), R_split 0.120 -> 0.106 and the unrefined R-free against the deposition -0.001 / -0.024. Rotation merges checked bit-identical (three in-house sets with overloads). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
241 lines
12 KiB
C++
241 lines
12 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <catch2/catch_all.hpp>
|
|
|
|
#include <cmath>
|
|
#include <algorithm>
|
|
#include <random>
|
|
#include <tuple>
|
|
#include <vector>
|
|
|
|
#include "../common/CUDAWrapper.h"
|
|
#include "../common/DetectorSetup.h"
|
|
#include "../common/DiffractionExperiment.h"
|
|
#include "../common/PixelMask.h"
|
|
#include "../image_analysis/hot_pixels/HotPixels.h"
|
|
|
|
namespace {
|
|
constexpr int W = 257, H = 257, C = 128;
|
|
constexpr int NFRAMES = 60;
|
|
constexpr double OSC_DEG = 0.1, SPACING_DEG = 6.0; // 60 frames spread over a full turn
|
|
|
|
constexpr size_t I(int x, int y) { return static_cast<size_t>(y) * W + x; }
|
|
}
|
|
|
|
// A Poisson background with a powder ring, and on it four kinds of pixel: one that reads 60 counts
|
|
// high on every frame, one only 20 high - persistent, but too weak to make an outlier - one that holds
|
|
// the error value on every frame, one stuck at the saturation limit, and one crossed by a genuine
|
|
// reflection, saturated on a few frames. The last sits close to the rotation axis, where one reflection stays on a pixel for a
|
|
// long stretch of rotation - here nine consecutive sampled frames, more than the chance bound allows
|
|
// but fewer than the one-reflection bound at that zeta. Only the first, third and fourth are masked - and
|
|
// a 3x3 patch reading high on every frame, which is stationary in the lab and so no reflection either.
|
|
// Three dark pixels besides: one dead on the powder ring, which is masked, one on the ring reading 70%
|
|
// of it, and one dead on the 5-count background, whose deficit is too small to make an outlier.
|
|
TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") {
|
|
DiffractionExperiment x(DetDECTRIS(W, H, "Test detector", ""));
|
|
x.IncidentEnergy_keV(WVL_1A_IN_KEV).DetectorDistance_mm(10.0f);
|
|
x.BeamX_pxl(static_cast<float>(C)).BeamY_pxl(static_cast<float>(C));
|
|
x.Goniometer(GoniometerAxis("omega", 0.0f, static_cast<float>(OSC_DEG), Coord(1, 0, 0), std::nullopt));
|
|
const PixelMask pixel_mask(x);
|
|
HotPixelFinder finder(x, pixel_mask, 4);
|
|
|
|
constexpr int HOT_X = 200, HOT_Y = 200, WARM_X = 190, WARM_Y = 60, ERR_X = 60, ERR_Y = 190;
|
|
constexpr int BRAGG_X = 40, BRAGG_Y = 115, STUCK_X = 100, STUCK_Y = 30;
|
|
constexpr int PATCH_X = 150, PATCH_Y = 225;
|
|
constexpr int DEAD_X = C + 82, DEAD_Y = C, DIM_X = C, DIM_Y = C + 82, DEAD_BKG_X = 20, DEAD_BKG_Y = 20;
|
|
std::mt19937 rng(1);
|
|
// The powder ring is a smooth radial profile, as a real one is: 45 counts over the background at
|
|
// 82 px, 4 px sigma.
|
|
std::vector<int32_t> frame(static_cast<size_t>(W) * H), scratch;
|
|
for (int f = 0; f < NFRAMES; f++) {
|
|
for (int y = 0; y < H; y++)
|
|
for (int x_ = 0; x_ < W; x_++) {
|
|
const double r = std::hypot(x_ - C, y - C);
|
|
const double mean = 5.0 + 45.0 * std::exp(-0.5 * (r - 82.0) * (r - 82.0) / 16.0);
|
|
frame[I(x_, y)] = std::poisson_distribution<int>(mean)(rng);
|
|
}
|
|
frame[I(HOT_X, HOT_Y)] += 60;
|
|
frame[I(WARM_X, WARM_Y)] += 20;
|
|
frame[I(ERR_X, ERR_Y)] = INT32_MIN;
|
|
frame[I(STUCK_X, STUCK_Y)] = INT32_MAX;
|
|
for (int dy = -1; dy <= 1; dy++)
|
|
for (int dx = -1; dx <= 1; dx++)
|
|
frame[I(PATCH_X + dx, PATCH_Y + dy)] += 60;
|
|
if (f >= 20 && f < 29)
|
|
frame[I(BRAGG_X, BRAGG_Y)] = f < 23 ? INT32_MAX : frame[I(BRAGG_X, BRAGG_Y)] + 500;
|
|
frame[I(DEAD_X, DEAD_Y)] = 0;
|
|
frame[I(DIM_X, DIM_Y)] = static_cast<int32_t>(0.7 * frame[I(DIM_X, DIM_Y)]);
|
|
frame[I(DEAD_BKG_X, DEAD_BKG_Y)] = 0;
|
|
finder.AddImage(frame.data(), scratch);
|
|
}
|
|
|
|
const auto result = finder.GetMask(OSC_DEG, SPACING_DEG, 4);
|
|
CHECK(result.frames == NFRAMES);
|
|
CHECK(result.mask[I(HOT_X, HOT_Y)] == 1); // 1 = hot, 2 = error value, 3 = dead
|
|
CHECK(result.mask[I(ERR_X, ERR_Y)] == 2);
|
|
CHECK(result.mask[I(STUCK_X, STUCK_Y)] == 2);
|
|
CHECK(result.mask[I(WARM_X, WARM_Y)] == 0);
|
|
CHECK(result.mask[I(BRAGG_X, BRAGG_Y)] == 0);
|
|
for (int dy = -1; dy <= 1; dy++)
|
|
for (int dx = -1; dx <= 1; dx++)
|
|
CHECK(result.mask[I(PATCH_X + dx, PATCH_Y + dy)] == 1);
|
|
CHECK(result.mask[I(DEAD_X, DEAD_Y)] == 3);
|
|
CHECK(result.mask[I(DIM_X, DIM_Y)] == 0);
|
|
CHECK(result.mask[I(DEAD_BKG_X, DEAD_BKG_Y)] == 0);
|
|
CHECK(result.hot == 10);
|
|
CHECK(result.error == 2);
|
|
CHECK(result.dead == 1);
|
|
}
|
|
|
|
#ifdef JFJOCH_USE_CUDA
|
|
// The device path against the host one, on frames that take every branch of both: a background bright
|
|
// enough on one annulus that the host reads its rings' statistics off the exact selection rather than
|
|
// its histogram, negative counts, saturated and error pixels scattered at random, masked pixels, and a
|
|
// few hundred planted pixels whose excess and persistence straddle every threshold, and as many read at
|
|
// fractions of their counts that straddle the dead rule's, so that a level or a threshold one count off
|
|
// would move some of them across. The masks must be identical.
|
|
TEST_CASE("HotPixelFinder_DeviceMatchesHost", "[HotPixelFinder]") {
|
|
if (get_gpu_count() == 0) {
|
|
WARN("No CUDA GPU present. Skipping HotPixelFinder_DeviceMatchesHost");
|
|
return;
|
|
}
|
|
DiffractionExperiment x(DetDECTRIS(W, H, "Test detector", ""));
|
|
x.IncidentEnergy_keV(WVL_1A_IN_KEV).DetectorDistance_mm(10.0f);
|
|
x.BeamX_pxl(static_cast<float>(C) + 0.3f).BeamY_pxl(static_cast<float>(C) - 0.6f);
|
|
x.Goniometer(GoniometerAxis("omega", 0.0f, static_cast<float>(OSC_DEG), Coord(1, 0, 0), std::nullopt));
|
|
PixelMask pixel_mask(x);
|
|
std::vector<uint32_t> masked(static_cast<size_t>(W) * H, 0);
|
|
for (int y = 100; y < 110; y++)
|
|
masked[I(30, y)] = 1;
|
|
pixel_mask.LoadUserMask(x, masked);
|
|
|
|
HotPixelFinder host(x, pixel_mask, 4), device(x, pixel_mask, 4);
|
|
HotPixelFinderGPU::Frame frame(std::make_shared<CudaStream>());
|
|
CudaDevicePtr<int32_t> device_image(static_cast<size_t>(W) * H);
|
|
|
|
std::mt19937 rng(7);
|
|
std::uniform_int_distribution<int> coord(0, W - 1);
|
|
struct Planted { size_t i; int excess; double rate; };
|
|
std::vector<Planted> planted;
|
|
for (int p = 0; p < 300; p++)
|
|
planted.push_back({I(coord(rng), coord(rng)), 5 + p / 3, 0.2 + 0.8 * (p % 7) / 6.0});
|
|
std::vector<std::pair<size_t, double>> dark;
|
|
for (int p = 0; p < 300; p++)
|
|
dark.push_back({I(coord(rng), coord(rng)), 0.3 + 0.3 * (p % 13) / 12.0});
|
|
|
|
std::vector<int32_t> image(static_cast<size_t>(W) * H), scratch;
|
|
std::uniform_real_distribution<double> u(0.0, 1.0);
|
|
for (int f = 0; f < NFRAMES; f++) {
|
|
for (int y = 0; y < H; y++)
|
|
for (int x_ = 0; x_ < W; x_++) {
|
|
const double r = std::hypot(x_ - C, y - C);
|
|
const double mean = 5.0 + 45.0 * std::exp(-0.5 * (r - 82.0) * (r - 82.0) / 16.0)
|
|
+ (r > 40.0 && r < 50.0 ? 3000.0 : 0.0);
|
|
int32_t v = std::poisson_distribution<int>(mean)(rng) - (r > 110.0 ? 3 : 0);
|
|
const double d = u(rng);
|
|
if (d < 0.002) v = INT32_MIN;
|
|
else if (d < 0.003) v = INT32_MAX;
|
|
image[I(x_, y)] = v;
|
|
}
|
|
for (const auto &p : planted)
|
|
if (u(rng) < p.rate && image[p.i] != INT32_MIN && image[p.i] != INT32_MAX)
|
|
image[p.i] += p.excess;
|
|
for (const auto &[i, fraction] : dark)
|
|
if (image[i] != INT32_MIN && image[i] != INT32_MAX)
|
|
image[i] = static_cast<int32_t>(fraction * image[i]);
|
|
for (size_t i = 0; i < image.size(); i++)
|
|
if (pixel_mask.GetMask()[i] != 0)
|
|
image[i] = INT32_MIN; // as the preprocessor leaves a masked pixel
|
|
|
|
host.AddImage(image.data(), scratch);
|
|
REQUIRE(cudaMemcpy(device_image, image.data(), image.size() * sizeof(int32_t), cudaMemcpyHostToDevice)
|
|
== cudaSuccess);
|
|
device.AddDeviceImage(device_image, frame);
|
|
// The frame is queued, not processed: the next copy must not overwrite it before its kernels ran.
|
|
REQUIRE(cudaStreamSynchronize(*frame.stream) == cudaSuccess);
|
|
}
|
|
|
|
const auto expected = host.GetMask(OSC_DEG, SPACING_DEG, 4);
|
|
const auto result = device.GetMask(OSC_DEG, SPACING_DEG, 4);
|
|
CHECK(expected.hot > 10);
|
|
CHECK(expected.dead > 3);
|
|
CHECK(result.dead == expected.dead);
|
|
CHECK(result.frames == expected.frames);
|
|
CHECK(result.hot == expected.hot);
|
|
CHECK(result.error == expected.error);
|
|
CHECK(result.mask == expected.mask);
|
|
}
|
|
|
|
// The pre-scan keeps one device half per card and folds the cards together when the mask is read. Two
|
|
// halves, each given part of the frames, folded into one must give exactly the chance counts and the
|
|
// candidates of one half given every frame.
|
|
TEST_CASE("HotPixelFinderGPU_FoldedCardsMatchOneCard", "[HotPixelFinder]") {
|
|
if (get_gpu_count() == 0) {
|
|
WARN("No CUDA GPU present. Skipping HotPixelFinderGPU_FoldedCardsMatchOneCard");
|
|
return;
|
|
}
|
|
// Four rings of two sectors: a key by band of rows and half of the detector, one column masked.
|
|
constexpr int NRINGS = 4, SECTORS = 2;
|
|
const size_t npixels = static_cast<size_t>(W) * H;
|
|
std::vector<int32_t> key(npixels);
|
|
std::vector<uint32_t> key_begin(NRINGS * SECTORS + 1, 0);
|
|
for (int y = 0; y < H; y++)
|
|
for (int x_ = 0; x_ < W; x_++) {
|
|
const int k = x_ == 7 ? -1 : (y * NRINGS / H) * SECTORS + x_ * SECTORS / W;
|
|
key[I(x_, y)] = k;
|
|
if (k >= 0) key_begin[k + 1]++;
|
|
}
|
|
for (int k = 0; k < NRINGS * SECTORS; k++)
|
|
key_begin[k + 1] += key_begin[k];
|
|
const HotPixelLevelRules rules{HotPixelFinder::MIN_SECTOR_PIXELS, HotPixelFinder::MIN_RING_PIXELS,
|
|
HotPixelFinder::LIT_NSIGMA, HotPixelFinder::LIT_OFFSET};
|
|
HotPixelFinderGPU one(key.data(), npixels, key_begin, NRINGS, SECTORS, rules);
|
|
HotPixelFinderGPU first(key.data(), npixels, key_begin, NRINGS, SECTORS, rules);
|
|
HotPixelFinderGPU second(key.data(), npixels, key_begin, NRINGS, SECTORS, rules);
|
|
HotPixelFinderGPU::Frame frame(std::make_shared<CudaStream>());
|
|
CudaDevicePtr<int32_t> device_image(npixels);
|
|
|
|
std::mt19937 rng(11);
|
|
std::uniform_real_distribution<double> u(0.0, 1.0);
|
|
std::vector<int32_t> image(npixels);
|
|
for (int f = 0; f < NFRAMES; f++) {
|
|
for (size_t i = 0; i < npixels; i++) {
|
|
const double d = u(rng);
|
|
image[i] = d < 0.01 ? INT32_MIN : d < 0.015 ? INT32_MAX : std::poisson_distribution<int>(8.0)(rng);
|
|
}
|
|
for (int p = 0; p < 40; p++)
|
|
image[I(20 + p, 30 + p)] += f % 3 == 0 ? 0 : 40; // lit on most frames
|
|
REQUIRE(cudaMemcpy(device_image, image.data(), npixels * sizeof(int32_t), cudaMemcpyHostToDevice)
|
|
== cudaSuccess);
|
|
one.Add(device_image, frame);
|
|
(f % 3 == 0 ? second : first).Add(device_image, frame);
|
|
REQUIRE(cudaStreamSynchronize(*frame.stream) == cudaSuccess);
|
|
}
|
|
first.Fold(second);
|
|
|
|
std::vector<int64_t> lit_one, seen_one, lit_folded, seen_folded;
|
|
one.ChanceCounts(lit_one, seen_one);
|
|
first.ChanceCounts(lit_folded, seen_folded);
|
|
CHECK(lit_folded == lit_one);
|
|
CHECK(seen_folded == seen_one);
|
|
const std::vector<int> k_chance(NRINGS, 5);
|
|
const auto expected = one.GetCandidates(NFRAMES, 10, true, k_chance);
|
|
const auto result = first.GetCandidates(NFRAMES, 10, true, k_chance);
|
|
CHECK(expected.index.size() >= 40);
|
|
// The candidates are written by atomic slot, so their order is the device's; compare them by pixel.
|
|
const auto by_pixel = [](const HotPixelFinderGPU::Candidates &c) {
|
|
std::vector<std::tuple<uint32_t, uint16_t, uint16_t, uint16_t, int64_t, int64_t>> v;
|
|
for (size_t j = 0; j < c.index.size(); j++)
|
|
v.emplace_back(c.index[j], c.n_lit[j], c.n_error[j], c.n_error_ring_ok[j], c.sum_value[j],
|
|
c.error_level_sum[j]);
|
|
std::sort(v.begin(), v.end());
|
|
return v;
|
|
};
|
|
CHECK(by_pixel(result) == by_pixel(expected));
|
|
CHECK(result.key_frames == expected.key_frames);
|
|
CHECK(result.key_level_sum == expected.key_level_sum);
|
|
CHECK(result.key_local_level_sum == expected.key_local_level_sum);
|
|
}
|
|
#endif
|