diff --git a/image_analysis/image_preprocessing/CMakeLists.txt b/image_analysis/image_preprocessing/CMakeLists.txt index fe6dd3ea9..7fa931935 100644 --- a/image_analysis/image_preprocessing/CMakeLists.txt +++ b/image_analysis/image_preprocessing/CMakeLists.txt @@ -2,7 +2,9 @@ ADD_LIBRARY(JFJochImagePreprocessing STATIC ImagePreprocessor.cpp ImagePreproces ImagePreprocessorCPU.cpp ImagePreprocessorCPU.h ImagePreprocessorBuffer.cpp - ImagePreprocessorBuffer.h) + ImagePreprocessorBuffer.h + SaturationCeiling.cpp + SaturationCeiling.h) TARGET_LINK_LIBRARIES(JFJochImagePreprocessing JFJochCommon) diff --git a/image_analysis/image_preprocessing/SaturationCeiling.cpp b/image_analysis/image_preprocessing/SaturationCeiling.cpp new file mode 100644 index 000000000..7d07dce25 --- /dev/null +++ b/image_analysis/image_preprocessing/SaturationCeiling.cpp @@ -0,0 +1,46 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "SaturationCeiling.h" + +#include +#include + +namespace saturation_ceiling { + +namespace { + +int Bin(int64_t value) { + return static_cast(std::floor(std::log2(static_cast(value) / FLOOR) * BINS_PER_OCTAVE)); +} + +int64_t LowerEdge(int bin) { + return static_cast(std::ceil(FLOOR * std::exp2(static_cast(bin) / BINS_PER_OCTAVE))); +} + +} // namespace + +void Add(Histogram &h, int64_t value) { + if (value < FLOOR) + return; + const int b = Bin(value); + if (b < N_BINS) + h[b]++; +} + +void Add(Histogram &h, const std::vector &values, const std::vector &mask) { + for (size_t i = 0; i < values.size(); i++) + if (i >= mask.size() || mask[i] == 0) + Add(h, values[i]); +} + +std::optional Ceiling(const Histogram &h) { + for (int b = WINDOW; b < N_BINS; b++) { + const uint64_t below = *std::max_element(h.begin() + b - WINDOW, h.begin() + b); + if (h[b] >= MIN_PIXELS && static_cast(h[b]) > RISE * static_cast(below)) + return LowerEdge(b); + } + return std::nullopt; +} + +} // namespace saturation_ceiling diff --git a/image_analysis/image_preprocessing/SaturationCeiling.h b/image_analysis/image_preprocessing/SaturationCeiling.h new file mode 100644 index 000000000..f61ecbb25 --- /dev/null +++ b/image_analysis/image_preprocessing/SaturationCeiling.h @@ -0,0 +1,39 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include + +// Where a detector's pixel values stop, read off the values themselves. The bright tail of a frame - +// the cores of Bragg peaks - falls off steadily with value; a detector that saturates clips everything +// beyond its ceiling down to it, so the values pile up just below where they stop. Where each pixel has +// its own ceiling (an ADC limit less that pixel's dark), the pile-up is a band rather than one value, +// and its lower edge is where clipping starts. A file's own statement of the limit is often missing +// or wrong for detectors outside the facility that wrote it. +namespace saturation_ceiling { + +// Values are binned on a log scale, BINS_PER_OCTAVE to a factor of two, from FLOOR up: a bin of the +// tail is then about 1% wide, narrow enough to place the edge, while below FLOOR a bin is narrower +// than one integer and the counts jump between neighbouring bins for no reason in the data. +constexpr int BINS_PER_OCTAVE = 64; +constexpr int64_t FLOOR = 256; +constexpr int N_BINS = BINS_PER_OCTAVE * 23; // up to 2^31 +// A pile-up bin holds more than RISE times the most any bin in the WINDOW below it held - a tail that +// only falls never does that - and at least MIN_PIXELS pixel values. +constexpr double RISE = 2.0; +constexpr int WINDOW = 8; +constexpr uint64_t MIN_PIXELS = 50; + +using Histogram = std::array; + +void Add(Histogram &h, const std::vector &values, const std::vector &mask); +void Add(Histogram &h, int64_t value); + +// The lower edge of the lowest bin that rises above the tail below it, if any does. +std::optional Ceiling(const Histogram &h); + +} // namespace saturation_ceiling diff --git a/reader/JFJochCrystFELReader.cpp b/reader/JFJochCrystFELReader.cpp index cb8d5e0f3..62fcf5d7c 100644 --- a/reader/JFJochCrystFELReader.cpp +++ b/reader/JFJochCrystFELReader.cpp @@ -76,6 +76,11 @@ void JFJochCrystFELReader::ReadFiles(const std::string &path, const std::string // A Cheetah CXI file carries what one flat panel needs - the frames at /entry_1/data_1/data, the pixel // size, the photon energy and a camera length, all in SI units - but not where the beam is. The panel is // put with its middle on the beam, which is as good as no centre at all: the run measures it. +// +// Only a square frame is taken this way. A Cheetah file from a tiled detector (CSPAD, AGIPD, ePix, +// JUNGFRAU, MPCCD) holds its panels stacked in a slab, which read as one flat panel would place every +// pixel wrongly and say nothing; those slabs are not square, the single-panel CCDs are, and a tiled +// detector needs its geometry file (--crystfel-geometry) anyway. bool JFJochCrystFELReader::CanReadCXI(const std::string &path) { try { if (!std::filesystem::exists(path)) @@ -83,7 +88,10 @@ bool JFJochCrystFELReader::CanReadCXI(const std::string &path) { const auto files = ListFiles(path); std::unique_lock ul(hdf5_mutex); HDF5ReadOnlyFile file(files[0]); - return file.Exists("/entry_1/data_1/data") && file.GetDimension("/entry_1/data_1/data").size() == 3 + if (!file.Exists("/entry_1/data_1/data")) + return false; + const auto dim = file.GetDimension("/entry_1/data_1/data"); + return dim.size() == 3 && dim[1] == dim[2] && file.Exists("/entry_1/instrument_1/detector_1/x_pixel_size") && file.Exists("/entry_1/instrument_1/source_1/energy"); } catch (...) { @@ -260,6 +268,12 @@ void JFJochCrystFELReader::Open(const std::string &path, const std::string &geom detector.PixelSize_um(static_cast(1e6 / p0.res)); // Nothing a parallax or sensor-efficiency correction needs is stated in a geometry file. detector.SensorThickness_um(0); + // Where the geometry states no limit, the stored integer type still has one - a Cheetah file keeps + // its frames as int16, and a value clipped at 32767 is no measurement either. + if (!saturation.has_value() && !stored_float_ && stored_size_ > 0) { + const double type_max = std::exp2(8.0 * static_cast(stored_size_) - (stored_signed_ ? 1 : 0)) - 1; + saturation = type_max / *std::max_element(adu_per_photon_.begin(), adu_per_photon_.end()); + } if (saturation.has_value()) detector.SaturationLimit(static_cast(std::floor(saturation.value()))); detector.BitDepthImage(32); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 23e053bc3..c81c21744 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -858,6 +858,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const int prepass_images = (prepass_end - start_image + config_.stride - 1) / config_.stride; const int early_images = experiment_.IsRotationIndexing() ? prepass_images : images_to_process; + if (full && !experiment_.IsRotationIndexing() && !config_.live) + DetectSaturation(start_image, images_to_process); + // After any geometry refinement, so the shadow is found about the beam centre actually used, // and before the azimuthal mapping and the output mask below, which both read pixel_mask_. if (config_.detect_beam_stop.has_value() || config_.adaptive_integration_radius) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 49a9d1e26..4cac16c57 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -164,6 +164,9 @@ struct ProcessConfig { // two answers - the one comparison an on-failure trigger cannot make, because a centre wrong // along the spindle does not fail. bool beam_center_check = true; + // Whether a saturation level the frames show below the file's limit replaces it (DetectSaturation): + // for every input but a Jungfraujoch-written file, whose stated limit is trusted. + bool detect_saturation = false; // Adaptive integration radius (--adaptive-integration-radius). When set, the pre-scan measures how // wide the recorded spots are - r80, the radius holding 80 % of a spot's flux, over a fixed @@ -1378,6 +1381,8 @@ class Rugnux { const std::vector &frames); // Returns whether the run's geometry was moved off the file's. bool DecideStillsGeometry(int start_image, int images_to_process); + // Stills: where the frames' pixel values pile up, against the file's saturation limit (see there). + void DetectSaturation(int start_image, int images_to_process); // One full run of the pipeline (index + integrate + optional scale/merge/space-group + output). // write_output == false makes no outputs (ProcessResult::outputs) and streams nothing to the image sink; diff --git a/rugnux/RugnuxStills.cpp b/rugnux/RugnuxStills.cpp index d353fb020..440ddbcec 100644 --- a/rugnux/RugnuxStills.cpp +++ b/rugnux/RugnuxStills.cpp @@ -52,6 +52,7 @@ #include "../image_analysis/azint/AzIntEngineCPU.h" #include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" #include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" +#include "../image_analysis/image_preprocessing/SaturationCeiling.h" #ifdef JFJOCH_USE_CUDA #include "../image_analysis/image_preprocessing/ImagePreprocessorGPU.h" #include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" @@ -468,3 +469,77 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { experiment_.BeamX_pxl(x).BeamY_pxl(y).DetectorDistance_mm(distance); return x != file_x || y != file_y || distance != file_distance; } + +// Where this run's pixel values stop, read off the frames (saturation_ceiling): a sample spread over +// the run, every unmasked pixel of it. Measured on every stills run and reported beside what the file +// states. It replaces the file's limit only where the frames pile up below it and the file is not one +// Jungfraujoch wrote (config_.detect_saturation) - a file that states no limit, or the ADC's full +// range where each pixel's dark has already been taken off, otherwise integrates clipped peaks as if +// they were measured, and the strongest reflections at low resolution are the ones it spoils. +void Rugnux::DetectSaturation(int start_image, int images_to_process) { + Logger logger("Rugnux"); + TimingMark mark(logger, "saturation ceiling"); + constexpr int SAMPLE_FRAMES = 40; + const std::vector sample = select_equally_spaced_image_ordinals( + images_to_process, std::min(images_to_process, SAMPLE_FRAMES)); + const std::vector &mask = pixel_mask_.GetMask(); + + std::mutex m; + saturation_ceiling::Histogram total{}; + std::atomic next = 0; + auto worker = [&]() { + saturation_ceiling::Histogram h{}; + JFJochReaderRawImage img; + std::vector buffer; + std::vector values; + for (size_t i = next.fetch_add(1); i < sample.size() && !cancelled_; i = next.fetch_add(1)) { + if (!reader_.ReadRawImage(start_image + sample[i] * config_.stride, img)) + continue; + const size_t n = img.image.GetWidth() * img.image.GetHeight(); + const uint8_t *ptr = img.image.GetUncompressedPtr(buffer); + values.resize(n); + switch (img.image.GetMode()) { + case CompressedImageMode::Int32: + std::memcpy(values.data(), ptr, n * sizeof(int32_t)); + break; + case CompressedImageMode::Uint32: + for (size_t j = 0; j < n; j++) + values[j] = static_cast(std::min( + reinterpret_cast(ptr)[j], INT32_MAX)); + break; + case CompressedImageMode::Uint16: + std::copy_n(reinterpret_cast(ptr), n, values.begin()); + break; + case CompressedImageMode::Int16: + std::copy_n(reinterpret_cast(ptr), n, values.begin()); + break; + default: + continue; // a float or colour image has no integer ceiling to find + } + saturation_ceiling::Add(h, values, mask); + } + std::lock_guard lock(m); + for (size_t b = 0; b < h.size(); b++) + total[b] += h[b]; + }; + std::vector> futures; + for (int i = 0; i < config_.nthreads; ++i) + futures.push_back(std::async(std::launch::async, worker)); + for (auto &f : futures) + f.get(); + + const int64_t stated = experiment_.GetSaturationLimit(); + const std::optional ceiling = saturation_ceiling::Ceiling(total); + if (!ceiling) { + logger.Info("Saturation: the values of {} sampled frames fall off without piling up below any level - " + "nothing clips them under the file's limit of {}", sample.size(), stated); + return; + } + const bool adopt = config_.detect_saturation && *ceiling < stated; + logger.Info("Saturation: the values of {} sampled frames pile up from {} - the file's limit is {}{}", + sample.size(), *ceiling, stated, + adopt ? ", so pixels from there on are treated as overloaded" + : *ceiling < stated ? "; the file's is kept (written by Jungfraujoch)" : ""); + if (adopt) + experiment_.Detector().SaturationLimit(*ceiling); +} diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index ade463023..98bf3feac 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -151,7 +151,7 @@ void print_usage() { std::cout << " Detector geometry" << std::endl; std::cout << " --crystfel-geometry Read - an HDF5/CXI file, or a directory of them - as CrystFEL does with this geometry file: frames from its data path, panels placed at its corners, camera length, pixel size, gain, photon energy and bad regions from it. One wavelength and one camera length for the run; panel positions are taken as given" << std::endl; - std::cout << " (A Cheetah CXI file, or a directory of them, is read without one: one flat panel, with the pixel size, photon energy and camera length the file states and the beam centre measured; --detector-distance overrides a camera length the file has wrong)" << std::endl; + std::cout << " (A Cheetah CXI file with a square frame, or a directory of them, is read without one: one flat panel, with the pixel size, photon energy and camera length the file states and the beam centre measured; --detector-distance overrides a camera length the file has wrong)" << std::endl; std::cout << std::endl; std::cout << " Detector mask" << std::endl; @@ -3182,6 +3182,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { config.two_pass_rotation = two_pass_rotation; config.detect_beam_stop = detect_beam_stop; config.beam_center_check = beam_center_check; + config.detect_saturation = input_is_cbf || input_is_marccd || input_is_smv || input_is_crystfel || input_is_cxi; config.adaptive_integration_radius = adaptive_integration_radius; config.counts_in_adu = input_is_smv || input_is_marccd || input_is_xtc || ((input_is_crystfel || input_is_cxi) && !crystfel_reader.CountsInPhotons()); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b7eb6dae8..b440124ba 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -11,6 +11,7 @@ ADD_EXECUTABLE(jfjoch_portable_test EXCLUDE_FROM_ALL ZSTDCompressorTest.cpp JFJochReaderTest.cpp CrystFELGeometryTest.cpp + SaturationCeilingTest.cpp SweepLayoutTest.cpp CBORTest.cpp ImageSpotFinderCPUTest.cpp diff --git a/tests/SaturationCeilingTest.cpp b/tests/SaturationCeilingTest.cpp new file mode 100644 index 000000000..1a5f16cab --- /dev/null +++ b/tests/SaturationCeilingTest.cpp @@ -0,0 +1,52 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include +#include + +#include "../image_analysis/image_preprocessing/SaturationCeiling.h" + +// A bright tail, P(> v) ~ 1/v, as the cores of Bragg peaks give. +static saturation_ceiling::Histogram Tail(int64_t clip_from, int64_t clip_to) { + saturation_ceiling::Histogram h{}; + std::mt19937 rng(1); + std::uniform_real_distribution u(0.0, 1.0); + for (int i = 0; i < 2000000; i++) { + auto v = static_cast(std::min(1e9, 100.0 / (1.0 - u(rng)))); // P(> v) ~ 1/v above 100 + if (v >= clip_from) { + // Clipped at a ceiling of its own, somewhere in [clip_from, clip_to]. + v = clip_from + static_cast(u(rng) * static_cast(clip_to - clip_from)); + } + saturation_ceiling::Add(h, v); + } + return h; +} + +TEST_CASE("SaturationCeiling_NoClipping", "[portable]") { + CHECK(!saturation_ceiling::Ceiling(Tail(INT32_MAX, INT32_MAX)).has_value()); +} + +TEST_CASE("SaturationCeiling_OneLevel", "[portable]") { + const auto c = saturation_ceiling::Ceiling(Tail(12000, 12000)); + REQUIRE(c.has_value()); + CHECK(*c <= 12000); + CHECK(*c > 12000 * 0.98); +} + +TEST_CASE("SaturationCeiling_BandPerPixel", "[portable]") { + // Each pixel stops at its own level - the ADC's limit less that pixel's dark. + const auto c = saturation_ceiling::Ceiling(Tail(12400, 15600)); + REQUIRE(c.has_value()); + CHECK(*c <= 12400); + CHECK(*c > 12400 * 0.98); +} + +TEST_CASE("SaturationCeiling_MaskedPixelsLeftOut", "[portable]") { + saturation_ceiling::Histogram h{}; + const std::vector values(1000, 5000); + saturation_ceiling::Add(h, values, std::vector(1000, 1)); + CHECK(!saturation_ceiling::Ceiling(h).has_value()); + saturation_ceiling::Add(h, values, {}); + CHECK(saturation_ceiling::Ceiling(h).has_value()); +}