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); +}