From e97b902603b85643da2c0f46e08dee68110219ee Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 25 Oct 2026 02:45:20 +0200 Subject: [PATCH] Beam stop: a reflection whose signal disk reaches the shadow is not integrated; no-data pixels are dead, not shadow The profile fit keeps a signal disk that lost pixels to the mask and fits the pixels that remain. That holds for a gap or a dead pixel, but not for the beam-stop shadow: the mask is drawn where the background is dimmed, and the Bragg half-shadow does not end there, so the rest of such a disk is attenuated too. Measured on summed rocking events against their symmetry mates (stratified by resolution and zeta, the detector edge excluded), observations whose disk touches the shadow read low by 12-60 % on 5uth, 7qij, 5nw5, 3r6o, 3meb (CCD and PILATUS alike) and level on the in-house photon-counting sweeps; observations with a clean disk are unbiased except on 5nw5, whose half-shadow reaches ~20 px past the mask. The test sits in MXAnalysisWithoutFPGA::IntegrateOnImage, so rotation, stills, the _process.h5 and --mode scale see the same reflections; with no shadow in the mask (the broker) it does nothing. The shadow finder also reports pixels holding a detector's no-data mark (PILATUS CBF module gaps at -1): they are the darkest pixels there are. The pre-scan now masks those as dead (bit 10) instead of beam-stop shadow (bit 9), so that the shadow test does not drop every reflection clipped by a module gap and the reported shadow share is the shadow's. A mean of zero is not taken as dead: a deep shadow on a still counts nothing in the projected frames. CPU subset battery vs rc180 (REFMAC 0-cycle R-free): 3meb -0.0004, 5uth -0.0004, 6h2p -0.0003, 6oel -0.0002, 3r6o +0.0001, 7qij/7rji 0, 5nw5 +0.002 (its lowest-shell CC1/2 0.78 -> 0.89, ISa 8.4 -> 9.3); in-house refres ISa up 0.1-0.3 on lyso/thau/insu/cytc/myob, half_image refres R_meas 75 -> 64 %, no set lost quality. Stills (lcls242_s/t, lcls222, euxfel83, sacla236): R-free within +-0.002, the noise of a perturbed run. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu (cherry picked from commit a7610d847694b5540dacc9ef609344ea7c7ae8fd) --- docs/RUGNUX_ADVANCED.md | 2 +- image_analysis/MXAnalysisWithoutFPGA.cpp | 35 ++++++++++++++++++++++-- image_analysis/MXAnalysisWithoutFPGA.h | 6 ++++ rugnux/RugnuxPreScan.cpp | 23 ++++++++++++++++ rugnux/rugnux_cli.cpp | 2 +- tests/BraggIntegrationEngineCPUTest.cpp | 32 ++++++++++++++++++++++ 6 files changed, 96 insertions(+), 4 deletions(-) diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index 70c0c85a3..637637fb2 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -353,7 +353,7 @@ Detector mask: | Option | Description | | --- | --- | -| `--detect-beam-stop[=N\|off]` | Find the beam stop and its holder in a projection of N images and add them to the pixel mask as bit 9, so nothing shadowed by them is integrated. **On by default** (60 images); `=off` disables. Reflections behind the stop are attenuated but not flagged, so they integrate low with a plausible sigma and no existing rejection catches them | +| `--detect-beam-stop[=N\|off]` | Find the beam stop and its holder in a projection of N images and add them to the pixel mask as bit 9, so nothing shadowed by them is integrated. **On by default** (60 images); `=off` disables. A reflection whose signal disk reaches into the shadow is not integrated: the rest of such a disk lies in the stop's half-shadow and reads low. Pixels the finder reports that hold the detector's no-data mark (a PILATUS gap) are masked as dead instead | Geometry: diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index e56230cd4..66976fc29 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -535,8 +535,15 @@ void MXAnalysisWithoutFPGA::UpdateMaskResolution(const SpotFindingSettings &sett std::vector MXAnalysisWithoutFPGA::IntegrateOnImage(const std::vector &predicted, size_t npredicted, int64_t image_number) { const auto &modules = experiment.GetModules(); - if (modules.empty()) - return bragg_engine->Run(*preprocessor_buffer, predicted, npredicted, image_number); + const auto drop_shadow_clipped = [this](std::vector &reflections) { + DropShadowClipped(reflections, mask.GetMask(), static_cast(experiment.GetXPixelsNumConv()), + static_cast(experiment.GetYPixelsNumConv()), experiment.GetBraggIntegrationSettings().GetR1()); + }; + if (modules.empty()) { + auto ret = bragg_engine->Run(*preprocessor_buffer, predicted, npredicted, image_number); + drop_shadow_clipped(ret); + return ret; + } // Predictions are where the reflections fall on the detector plane; the engine reads the image, so // each goes to where its panel put those pixels, and comes back to the plane once integrated. One in @@ -557,6 +564,7 @@ std::vector MXAnalysisWithoutFPGA::IntegrateOnImage(const std::vecto on_image.push_back(r); } auto ret = bragg_engine->Run(*preprocessor_buffer, on_image, on_image.size(), image_number); + drop_shadow_clipped(ret); for (auto &r: ret) { const int64_t panel = modules.ModuleOfImagePixel(r.predicted_x, r.predicted_y); const auto [px, py] = modules.ImageToPlane(panel, r.predicted_x, r.predicted_y); @@ -570,3 +578,26 @@ std::vector MXAnalysisWithoutFPGA::IntegrateOnImage(const std::vecto } return ret; } + +// A signal disk that lost pixels to the mask is otherwise kept, its profile fitted over the pixels that +// remain. That holds for a gap or a dead pixel, which take nothing from their neighbours, but not for the +// beam-stop shadow: the mask is drawn where the BACKGROUND is dimmed, and the Bragg signal's half-shadow +// does not end where the background's does, so the rest of such a disk is attenuated too. Measured against +// their symmetry mates, these observations came out 12-60 % low on the sweeps whose shadow edge crosses +// busy low-resolution reflections (CCD and pixel detectors alike), and level with them elsewhere. +void DropShadowClipped(std::vector &reflections, const std::vector &mask, + int W, int H, float r1) { + std::erase_if(reflections, [&](const Reflection &r) { + const int x0 = std::max(0, static_cast(std::floor(r.predicted_x - r1))); + const int x1 = std::min(W - 1, static_cast(std::ceil(r.predicted_x + r1))); + const int y0 = std::max(0, static_cast(std::floor(r.predicted_y - r1))); + const int y1 = std::min(H - 1, static_cast(std::ceil(r.predicted_y + r1))); + for (int y = y0; y <= y1; y++) + for (int x = x0; x <= x1; x++) { + const float dx = static_cast(x) - r.predicted_x, dy = static_cast(y) - r.predicted_y; + if (dx * dx + dy * dy < r1 * r1 && (mask[static_cast(y) * W + x] & (1u << PixelMask::BeamStopPixelBit))) + return true; + } + return false; + }); +} diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index 3d3b6e802..27cdff2e5 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/image_analysis/MXAnalysisWithoutFPGA.h @@ -26,6 +26,12 @@ class CudaStream; class AdaptiveSpotFinderGPU; +// Leave out the reflections whose signal disk - radius r1 about the predicted position, in image +// coordinates - reaches a pixel of the beam-stop shadow (PixelMask::BeamStopPixelBit of `mask`, a +// width x height image). +void DropShadowClipped(std::vector &reflections, const std::vector &mask, + int width, int height, float r1); + // MXAnalysisWithoutFPGA is not thread safe - it has to owned by a single thread class MXAnalysisWithoutFPGA { const DiffractionExperiment &experiment; diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index 9534da056..9cafe9646 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -978,6 +978,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // Where the scattered background puts the beam, when the projection exists to say so. std::optional> ring_centre; std::vector opaque_shadow; + std::vector dead_pixels; // reported by the shadow finder, counted nothing (below) // The mean projection, taken once and shared read-only by everything below that reads it - the // ring-centre fit here, and the diagnostic picture and the beam-centre capture on threads of // their own. It is 4 bytes a pixel, and a copy for each of them was a copy on this thread. @@ -1052,6 +1053,22 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru Note("polarization plane taken from the background, turned by 90 degrees"); experiment_.PolarizationFactor(polarization); } + // A pixel whose mean is negative holds the detector's own mark for no data - a module gap a PILATUS + // CBF leaves at -1 - and is the darkest pixel there is, so the shadow finder reports it. It is dead + // rather than shadowed, and goes to the defective pixels (after MaskDefectivePixels below), so that + // the shadow bit holds only what something in the beam shadowed: what the share reported here + // counts, and what the integration's shadow test reads (MXAnalysisWithoutFPGA::DropShadowClipped). + // A mean of zero is not enough: a deep shadow on a low background counts nothing in the frames + // projected, and on a still it often does. + dead_pixels.assign(shadow.size(), 0); + for (size_t i = 0; i < shadow.size(); i++) + if (shadow[i] != 0 && !((*projection)[i] > -0.5f)) { + shadow[i] = 0; + dead_pixels[i] = 1; + } + if (const auto dead = std::count(dead_pixels.begin(), dead_pixels.end(), 1u); dead > 0) + logger.Info("Beam stop: {} pixels reading the detector's no-data mark in {} images are masked as dead, " + "not as shadow", dead, shadow_frames); const auto shadowed = std::count_if(shadow.begin(), shadow.end(), [](uint32_t v) { return v != 0; }); pixel_mask_.LoadBeamStopMask(experiment_, shadow); opaque_shadow = shadow; @@ -1107,6 +1124,12 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru measured.emplace(); MaskDefectivePixels(start_image, shadow_sample, ring_centre, frame, measured->defects); } + if (!dead_pixels.empty()) { + const auto &m = pixel_mask_.GetMask(); + for (size_t i = 0; i < m.size(); i++) + dead_pixels[i] |= (m[i] >> PixelMask::HotPixelBit) & 1u; + pixel_mask_.LoadHotPixelMask(experiment_, dead_pixels); + } // The centre is fitted with the shadow masked but NOT the transmitting pieces the shadow finder // adds after it. Their dimming is multiplicative in azimuth, which the fit's per-sector amplitude diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 58f89b955..be30055f0 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -167,7 +167,7 @@ void print_usage() { std::cout << " Detector mask" << std::endl; std::cout << " --mask-good [:][@] HDF5/CXI stills: a mask array in the frames' own layout (as stored, before the panels are placed) - in , or without it in each data file (a per-frame mask stored beside the frames); a pixel is usable only where it is non-zero, or with only where it has all of them. May be given more than once" << std::endl; std::cout << " --mask-bad [:][@] As --mask-good, but a pixel is unusable where the array is non-zero, or with where it has any of them (e.g. @0xfeff)" << std::endl; - std::cout << " --detect-beam-stop[=N|off] Find the beam stop and its holder in a projection of N images and add them to the pixel mask (bit 9), so nothing shadowed by them is integrated. ON by default (60 images); =off disables. Reflections behind the stop are attenuated but not flagged, so they are integrated low with a plausible sigma and no existing rejection catches them" << std::endl; + std::cout << " --detect-beam-stop[=N|off] Find the beam stop and its holder in a projection of N images and add them to the pixel mask (bit 9), so nothing shadowed by them is integrated. ON by default (60 images); =off disables. A reflection whose signal disk reaches into the shadow is not integrated: the rest of such a disk lies in the stop's half-shadow and reads low. Pixels the finder reports that hold the detector's no-data mark (a PILATUS gap) are masked as dead instead" << std::endl; std::cout << " --beam-center-check[=off] Measure the beam centre from the isotropy of the scattered background on EVERY run, report how far the file's value is from it, and index a second first pass with it to see whether the two centres give the same lattice. The fit reads the projection --detect-beam-stop already builds, so it costs no extra frames. On a run that indexes, the measured centre is adopted where the file's centre indexes nothing and the measured one indexes a majority, and where the two centres return cells related by an integer volume factor, in either direction, and the measured centre's lattice carries materially more of the pooled spots against its own chance level - the axis harmonic a centre error along the spindle produces, which nothing downstream repairs. Where the indexing cannot tell the two centres apart - the same cell with a different metric symmetry, or the same lattice found at both only by reading less of each frame - the first pass is run at both and the run keeps the one that merges better (on the same lattice: the one that measures more reflections at I/sigma >= 2). ON by default; =off disables On stills, where there is no first pass, the file's centre, the measured one and the point-symmetry capture's are each indexed on a sample of frames and the run takes the one that puts materially more spots on the lattice of their frame; with a reference cell (-C) the detector distance is searched too, over 30-3000 mm at every centre, the file's value being one hypothesis, before the stills geometry refinement" << std::endl; std::cout << std::endl; diff --git a/tests/BraggIntegrationEngineCPUTest.cpp b/tests/BraggIntegrationEngineCPUTest.cpp index c89320728..46d2f30b5 100644 --- a/tests/BraggIntegrationEngineCPUTest.cpp +++ b/tests/BraggIntegrationEngineCPUTest.cpp @@ -13,6 +13,8 @@ #include "../common/DiffractionExperiment.h" #include "../common/JFJochMath.h" #include "../common/Reflection.h" +#include "../common/PixelMask.h" +#include "../image_analysis/MXAnalysisWithoutFPGA.h" #include "../image_analysis/bragg_integration/BraggIntegrationEngineCPU.h" #include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" @@ -540,3 +542,33 @@ TEST_CASE("BraggIntegrationEngineCPU_AduIntegratesAsPhotons", "[Integration][por } } } + +// A reflection is left out only when its signal disk reaches the beam-stop shadow: a shadow pixel just +// outside the disk, or a pixel masked for any other reason inside it, keeps it. +TEST_CASE("DropShadowClipped_SignalDiskTouchingTheShadow", "[Integration][portable]") { + const int W = 64, H = 32; + const float r1 = 4.0f; + std::vector mask(W * H, 0); + for (int y = 0; y < H; y++) + mask[y * W + 30] = 1u << PixelMask::BeamStopPixelBit; // a shadow column at x = 30 + mask[16 * W + 10] = 1u << PixelMask::ModuleGapPixelBit; // a gap pixel at (10, 16) + + auto at = [](float x, float y) { + Reflection r{}; + r.predicted_x = x; + r.predicted_y = y; + return r; + }; + std::vector refl = {at(10.0f, 16.0f), // a gap pixel at its centre: kept + at(26.5f, 16.0f), // 3.5 px from the shadow: dropped + at(25.0f, 16.0f), // 5 px from it: kept + at(33.9f, 5.0f), // 3.9 px on the other side: dropped + at(34.0f, 5.0f), // 4 px, the disk's own edge: kept + at(1.0f, 1.0f)}; // near the image corner: kept + DropShadowClipped(refl, mask, W, H, r1); + REQUIRE(refl.size() == 4); + CHECK(refl[0].predicted_x == 10.0f); + CHECK(refl[1].predicted_x == 25.0f); + CHECK(refl[2].predicted_x == 34.0f); + CHECK(refl[3].predicted_x == 1.0f); +}