From 8999287ba459f2a47f4a6c3e488e05dbef23c7a9 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 10:01:49 +0200 Subject: [PATCH 001/179] Test: a failed calibration now leaves the machine in Error, not Inactive JFJochStateMachine_CalibrationFailure was written in rc.164 against the behaviour of the time. rc.165 moved every calibration FAILURE path - CalibrateDetector's catch, TakeDarkMaskInternal's mask-size check and the three ImportPedestal checks - from Inactive to Error, so that /wait_till_done and /wait_until_running answer with the reason instead of the bodiless 502 that a deliberate Deactivate() also leaves; only a CANCELLED calibration still ends in Inactive. The test kept asserting the old state and started failing. Assert Error. The rest of the case is unaffected: the message and severity are what the new SetState writes, and Start() still throws from Error because it requires Idle, so "a data collection is refused on it" keeps its meaning. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H --- tests/JFJochStateMachineTest.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/JFJochStateMachineTest.cpp b/tests/JFJochStateMachineTest.cpp index dd5f9ea5a..3b140ceec 100644 --- a/tests/JFJochStateMachineTest.cpp +++ b/tests/JFJochStateMachineTest.cpp @@ -336,10 +336,11 @@ TEST_CASE("JFJochStateMachine_CalibrationFailure") { REQUIRE_NOTHROW(state_machine.Initialize()); REQUIRE_NOTHROW(state_machine.WaitTillMeasurementDone()); - // Inactive, not Idle and not Error: the calibration is undefined, so the detector has to be - // initialized again rather than looking ready to measure. + // Error, not Idle: the calibration is undefined, so the detector has to be initialized again + // rather than looking ready to measure. Error and not Inactive, so that the wait calls report + // the reason instead of the bodiless answer a deliberate power-off leaves behind. auto status = state_machine.GetStatus(); - REQUIRE(status.state == JFJochState::Inactive); + REQUIRE(status.state == JFJochState::Error); REQUIRE(status.message_severity == BrokerStatus::MessageSeverity::Error); REQUIRE(status.message == "Pedestal not collected properly"); -- 2.54.0 From a6fbd1d195aeeda2206f653cf6eaf55017a88d24 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 10:18:19 +0200 Subject: [PATCH 002/179] Viewer: Alt and the wheel step through the images The wheel already zooms, and with Ctrl or Shift it moves the foreground; Alt now moves through the dataset, one image per notch, wheel up forward - the direction QAbstractSlider's own wheel handling uses, which is what the toolbar's scrub slider follows. The view does not know which image it is showing, so it emits the step and the navigation toolbar applies it: the same loadImage() every other control ends in, with the same clamp and the same Sum setting, so nothing about loading is duplicated. An empty dataset is left alone, because loadImage(-1) does not mean "before the first image" but "the latest one" when the viewer is following a running collection over HTTP. The signal is on JFJochImage, so the other views emit it too; only the diffraction view is connected, which leaves them as they were. Note for a desktop where Alt+wheel does nothing: many window managers grab Alt-modified mouse events before the application sees them, and that is a window-manager setting, not something the viewer can take back. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H --- docs/CHANGELOG.md | 3 +++ viewer/JFJochViewerWindow.cpp | 2 ++ viewer/image_viewer/JFJochImage.cpp | 10 ++++++++++ viewer/image_viewer/JFJochImage.h | 2 ++ viewer/toolbar/JFJochViewerToolbarImage.cpp | 10 ++++++++++ viewer/toolbar/JFJochViewerToolbarImage.h | 1 + 6 files changed, 28 insertions(+) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index a568fc030..14505b792 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,5 +1,8 @@ # Changelog ## 1.0.0 +### 1.0.0-rc.166 +* In `jfjoch_viewer`, Alt and the mouse wheel step through the dataset one image at a time. + ### 1.0.0-rc.165 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. diff --git a/viewer/JFJochViewerWindow.cpp b/viewer/JFJochViewerWindow.cpp index 3c8cfabb5..215557691 100644 --- a/viewer/JFJochViewerWindow.cpp +++ b/viewer/JFJochViewerWindow.cpp @@ -224,6 +224,8 @@ JFJochViewerWindow::JFJochViewerWindow(QWidget *parent, bool dbus, const QString connect(toolBarImage, &JFJochViewerToolbarImage::loadImage, reading_worker, &JFJochImageReadingWorker::LoadImage); + connect(viewer, &JFJochDiffractionImage::stepImage, toolBarImage, &JFJochViewerToolbarImage::stepImage); + connect(toolBarDisplay, &JFJochViewerToolbarDisplay::setForeground, viewer, &JFJochDiffractionImage::changeForeground); diff --git a/viewer/image_viewer/JFJochImage.cpp b/viewer/image_viewer/JFJochImage.cpp index e06eb3133..39e4c21ad 100644 --- a/viewer/image_viewer/JFJochImage.cpp +++ b/viewer/image_viewer/JFJochImage.cpp @@ -143,6 +143,16 @@ void JFJochImage::setFeatureColor(QColor input) { void JFJochImage::wheelEvent(QWheelEvent *event) { if (!scene()) return; + // Alt+wheel steps through the dataset instead of zooming. The view does not know which image + // it shows, so the step is emitted and the navigation toolbar applies it. + if (event->modifiers() & Qt::AltModifier) { + const int delta = event->angleDelta().y(); + if (delta != 0) + emit stepImage(delta > 0 ? 1 : -1); + event->accept(); + return; + } + const double zoomFactor = 1.15; // Zoom factor // Get the position of the mouse in scene coordinates diff --git a/viewer/image_viewer/JFJochImage.h b/viewer/image_viewer/JFJochImage.h index 8b43034ca..b8cb6ea67 100644 --- a/viewer/image_viewer/JFJochImage.h +++ b/viewer/image_viewer/JFJochImage.h @@ -230,6 +230,8 @@ signals: void hoverScenePos(QPointF scenePos); // A new frame has been rendered into Frame(). Follower views repaint on this. void frameRendered(); + // Alt+wheel asks for a move through the dataset: +1 one image forward, -1 one back. + void stepImage(int steps); private slots: void onScroll(int value); public slots: diff --git a/viewer/toolbar/JFJochViewerToolbarImage.cpp b/viewer/toolbar/JFJochViewerToolbarImage.cpp index 1fbddc77e..8ab28efd2 100644 --- a/viewer/toolbar/JFJochViewerToolbarImage.cpp +++ b/viewer/toolbar/JFJochViewerToolbarImage.cpp @@ -2,6 +2,7 @@ // SPDX-License-Identifier: GPL-3.0-only #include "JFJochViewerToolbarImage.h" +#include #include #include "../widgets/ToolbarIcons.h" @@ -176,6 +177,15 @@ void JFJochViewerToolbarImage::leftmostButtonPressed() { emit loadImage(0, sum); } +void JFJochViewerToolbarImage::stepImage(int steps) { + if (image_count_in_dataset == 0) + return; + const auto target = std::clamp(curr_image + steps, 0, + static_cast(image_count_in_dataset) - 1); + if (target != curr_image) + emit loadImage(target, sum); +} + void JFJochViewerToolbarImage::imageNumberSliderPressed() { image_number_slider_manual = true; } diff --git a/viewer/toolbar/JFJochViewerToolbarImage.h b/viewer/toolbar/JFJochViewerToolbarImage.h index 93792ef2b..88b9b9eb8 100644 --- a/viewer/toolbar/JFJochViewerToolbarImage.h +++ b/viewer/toolbar/JFJochViewerToolbarImage.h @@ -57,6 +57,7 @@ public: explicit JFJochViewerToolbarImage(QWidget *parent = nullptr); public slots: void setImageNumber(int64_t total_images, int64_t current_image); + void stepImage(int steps); void setAutoloadMode(JFJochImageReadingWorker::AutoloadMode input); void setHttpConnection(bool connected, QString addr); private slots: -- 2.54.0 From ac202a55a161b57db39aecb6c15744e4a32bfdda Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 10:38:10 +0200 Subject: [PATCH 003/179] Spot finding: the strong-pixel limit follows the detector An image with 65535 or more strong pixels was given up on and reported ZERO spots - silently, no log line, indistinguishable from a frame that did not diffract. 65535 is one pixel in 64 of the JUNGFRAU 4M the number was written for; left fixed while the detectors grew it became one in 276 of an 18-megapixel EIGER, which a strongly diffracting crystal passes on its best frames. On the strong rotation set just added to the battery it cost 767 of 1800 images: peakCountUnfiltered 0 and resolutionEstimate NaN across two blocks of the sweep, the two where the crystal diffracts hardest. Make the bar one pixel in 64 everywhere, and never below the value that stood here, so no smaller detector loses ground. It lived in three places - the host extractor, StrongPixelSet, and SpotExtractorGPU's buffer capacity - now one function. The bar was there for a reason and raising it alone would not have been safe. sparseccl walks a sliding window of the last two lines and tests every pixel in it, which is quadratic in how many strong pixels a line pair holds: a handful for the silicon-tracker hits upstream wrote it for, four thousand for a lit detector line, and 76 seconds for a fully lit frame. But the pixels arrive in raster order, so the window need not be walked at all - a pixel's earlier 8-neighbours are the one to its left and the at most three above it, which is what the GPU extractor already finds by binary search. Keeping the previous line's range and a forward-only cursor gives the same edge set and the same unions in the same order, so the labels are identical, and the fully lit frame now takes 0.16 s. Verified bit-identical on real frames, on fully dense frames, across occupancy 1e-5 to 5e-2, and on 4000 randomised images including ones with blank lines; SpotExtractorGPU's host-vs-device parity test passes untouched. ImagePreprocessorBufferGPU's gather staging was sized to the old constant, with a comment tying it to the caller's give-up. Raising that give-up without it would have run the gather off the end of the device buffer, so it follows the same limit now. Byte-identical .hkl on three battery crystals that never reach the bar. On the strong set, with symmetry, cell and geometry pinned so only the spot list moves: better in every resolution shell, CC1/2 97.8 -> 98.5%, R_meas 30.5 -> 28.4%, ISa 3.36 -> 3.58, indexing rate 0.772 -> 0.824. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H --- docs/CHANGELOG.md | 2 + .../ImagePreprocessorBufferGPU.cu | 5 +- .../ImagePreprocessorBufferGPU.h | 9 ++- .../spot_finding/ImageSpotFinder.cpp | 6 +- .../spot_finding/SpotExtractorGPU.cu | 29 ++++---- .../spot_finding/SpotExtractorGPU.h | 7 +- .../spot_finding/SpotFindingSettings.h | 16 +++++ .../spot_finding/StrongPixelSet.cpp | 66 +++++++++++-------- image_analysis/spot_finding/StrongPixelSet.h | 1 + 9 files changed, 89 insertions(+), 52 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 14505b792..a1dac33e0 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,8 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* Spot finding no longer throws away a whole image when it holds many strong pixels: the limit follows the detector (one pixel in 64) instead of standing at the 65535 that suited a 4-megapixel detector, which a strongly diffracting crystal on an 18-megapixel one passes on its best frames. +* The connected-component search is linear in the strong pixels rather than quadratic in how many of them a detector line holds; a fully lit image is labelled in 0.16 s instead of 76 s, and the spots it finds are unchanged. * In `jfjoch_viewer`, Alt and the mouse wheel step through the dataset one image at a time. ### 1.0.0-rc.165 diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu index 54265f5c3..e45f97dd3 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu +++ b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu @@ -16,8 +16,9 @@ ImagePreprocessorBufferGPU::ImagePreprocessorBufferGPU(size_t npixel, bool host_ gpu_image(npixel), // A no-op when the mirror was not allocated: CudaRegisteredVector skips an empty vector. buffer_reg(buffer), - gpu_gather_index(MAX_GATHER), - gpu_gather_value(MAX_GATHER) { + max_gather(StrongPixelLimit(npixel)), + gpu_gather_index(max_gather), + gpu_gather_value(max_gather) { } int32_t *ImagePreprocessorBufferGPU::getGPUBuffer() { diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h index dd223be15..45fbb7bc6 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h +++ b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h @@ -7,15 +7,18 @@ #include "ImagePreprocessorBuffer.h" #include "../indexing/CUDAMemHelpers.h" +#include "../spot_finding/SpotFindingSettings.h" class ImagePreprocessorBufferGPU : public ImagePreprocessorBuffer { CudaDevicePtr gpu_image; CudaRegisteredVector buffer_reg; // Staging for Gather(). Its only caller is ImageSpotFinder::ExtractSpots, which gives up on a frame - // with UINT16_MAX or more strong pixels (the connected-component search rejects it anyway), so that - // is the largest gather that can be asked for. - static constexpr size_t MAX_GATHER = UINT16_MAX; + // with StrongPixelLimit or more strong pixels (the connected-component search rejects it anyway), + // so that is the largest gather that can be asked for. It follows the detector, so this has to as + // well - sized to a constant while the caller's bar was raised, the gather would run off the end + // of the buffer rather than merely lose the frame. + const size_t max_gather; CudaDevicePtr gpu_gather_index; CudaDevicePtr gpu_gather_value; // Own stream: every analysis engine synchronises its own stream before it returns, so the device diff --git a/image_analysis/spot_finding/ImageSpotFinder.cpp b/image_analysis/spot_finding/ImageSpotFinder.cpp index 47a44d3f7..8fe18f98a 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.cpp +++ b/image_analysis/spot_finding/ImageSpotFinder.cpp @@ -71,9 +71,9 @@ void ImageSpotFinder::ExtractComponentsHost(const ImagePreprocessorBuffer &image } components.clear(); - // The connected-component search rejects a frame with this many strong pixels, so their values are - // of no use. - if (strong_pixel.size() >= UINT16_MAX) + // The connected-component search gives up on a frame with this many strong pixels, so their values + // are of no use - not even worth gathering off the device. + if (strong_pixel.size() >= StrongPixelLimit(static_cast(width) * height)) return; image.Gather(strong_pixel, strong_pixel_value); diff --git a/image_analysis/spot_finding/SpotExtractorGPU.cu b/image_analysis/spot_finding/SpotExtractorGPU.cu index 2a5e9c38c..2795de1b9 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.cu +++ b/image_analysis/spot_finding/SpotExtractorGPU.cu @@ -185,7 +185,7 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 const int n = static_cast(min(*nstrong, capacity)); if (threadIdx.x == 0) *nout = 0; __syncthreads(); - // Same give-up as StrongPixelSet::FindComponentsImage - except that here the count is known + // The same give-up the host makes at StrongPixelLimit - except that here the count is known // before a single pixel has been written anywhere, so the frame costs nothing to reject. if (n == 0 || static_cast(n) >= capacity) return; @@ -275,16 +275,17 @@ SpotExtractorGPU::SpotExtractorGPU(int32_t in_width, int32_t in_height, std::sha : stream(std::move(in_stream)), width(in_width), nwords((static_cast(in_width) * in_height + 31) / 32), + max_strong(StrongPixelLimit(static_cast(in_width) * in_height)), gpu_res_mask(nwords), gpu_nstrong(1), - gpu_index(MAX_STRONG), - gpu_value(MAX_STRONG), - gpu_parent(MAX_STRONG), - gpu_root(MAX_STRONG), - gpu_label(MAX_STRONG), - gpu_count(MAX_STRONG), - gpu_spot(MAX_STRONG), - gpu_spot_out(MAX_STRONG), + gpu_index(max_strong), + gpu_value(max_strong), + gpu_parent(max_strong), + gpu_root(max_strong), + gpu_label(max_strong), + gpu_count(max_strong), + gpu_spot(max_strong), + gpu_spot_out(max_strong), gpu_nspot(1), host_nspot(1), host_spot(SPOT_PREFIX) { @@ -326,15 +327,15 @@ void SpotExtractorGPU::Extract(const uint32_t *gpu_strong, const int32_t *gpu_im scan_block_counts<<<1, FINISH_THREADS, 0, *stream>>>(gpu_block_count, gpu_block_offset, gpu_nstrong, compact_blocks); scatter_bits<<>>(gpu_strong, gpu_res_mask, gpu_block_offset, gpu_image, - gpu_index, gpu_value, nwords, MAX_STRONG); + gpu_index, gpu_value, nwords, max_strong); // Fixed grids reading the strong-pixel count from device memory: the host never learns it, so it // never has to synchronise in the middle of the frame. - init_parent<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_nstrong, MAX_STRONG); - union_neighbours<<<512, THREADS, 0, *stream>>>(gpu_index, gpu_parent, gpu_nstrong, MAX_STRONG, width); - resolve_roots<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_root, gpu_nstrong, MAX_STRONG); + init_parent<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_nstrong, max_strong); + union_neighbours<<<512, THREADS, 0, *stream>>>(gpu_index, gpu_parent, gpu_nstrong, max_strong, width); + resolve_roots<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_root, gpu_nstrong, max_strong); finish_components<<<1, FINISH_THREADS, 0, *stream>>>(gpu_index, gpu_value, gpu_root, gpu_label, gpu_count, gpu_spot, gpu_spot_out, gpu_nspot, - gpu_nstrong, MAX_STRONG, width, max_pix); + gpu_nstrong, max_strong, width, max_pix); cuda_err(cudaMemcpyAsync(host_nspot, gpu_nspot, sizeof(uint32_t), cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaMemcpyAsync(host_spot, gpu_spot_out, SPOT_PREFIX * sizeof(SpotExtractorGPUSpot), cudaMemcpyDeviceToHost, *stream)); diff --git a/image_analysis/spot_finding/SpotExtractorGPU.h b/image_analysis/spot_finding/SpotExtractorGPU.h index e00a86cc6..4cbc93735 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.h +++ b/image_analysis/spot_finding/SpotExtractorGPU.h @@ -51,9 +51,10 @@ class SpotExtractorGPU { const int32_t width; const size_t nwords; - // The connected-component search gives up above this many strong pixels (see - // StrongPixelSet::FindComponentsImage), so nothing larger is ever built. - static constexpr uint32_t MAX_STRONG = UINT16_MAX; + // Strong pixels this engine's buffers hold, and above which the extraction gives up on the frame - + // StrongPixelLimit, so it follows the detector rather than standing at a constant. 104 bytes of + // device memory apiece, 29 MB on an 18-megapixel detector. + const uint32_t max_strong; // Spots copied back together with their count in one transfer. A frame with more than this many // surviving spots - far past anything indexable - simply takes a second copy. static constexpr uint32_t SPOT_PREFIX = 4096; diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index 5ecc300d7..a2bdde060 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/image_analysis/spot_finding/SpotFindingSettings.h @@ -3,9 +3,25 @@ #pragma once +#include #include #include +// Strong pixels above which the connected-component search gives up on a frame, unlabelled: an image +// with this much of the detector over threshold is not a diffraction pattern, and it is what the +// spot extractor's buffers are sized to. It is no longer a time limit - the search is linear in the +// strong pixels either way, and a fully lit 18-megapixel frame labels in 0.16 s. +// +// The bar has to be a FRACTION of the detector. It stood at a fixed 65535, which is one pixel in 64 +// of the JUNGFRAU 4M it was written for; left fixed while the detectors grew it became one pixel in +// 276 of an 18-megapixel EIGER - a bar a strongly diffracting crystal clears on its best frames, +// which were then dropped whole, and in silence. One in 64 everywhere, and never below the value +// that used to stand here, so no smaller detector loses ground. +constexpr uint32_t StrongPixelLimit(size_t pixel_count) { + const auto limit = static_cast(pixel_count / 64); + return limit > UINT16_MAX ? limit : UINT16_MAX; +} + struct SpotFindingSettings { bool enable = true; float signal_to_noise_threshold = 4.0; // STRONG_PIXEL in XDS diff --git a/image_analysis/spot_finding/StrongPixelSet.cpp b/image_analysis/spot_finding/StrongPixelSet.cpp index 7af8acfba..6c631236a 100644 --- a/image_analysis/spot_finding/StrongPixelSet.cpp +++ b/image_analysis/spot_finding/StrongPixelSet.cpp @@ -4,6 +4,9 @@ // SparseCCL code taken from https://github.com/acts-project/traccc/blob/main/core/include/traccc/clusterization/detail/sparse_ccl.hpp // (c) 2021-2022 CERN for the benefit of the ACTS project // Mozilla Public License Version 2.0 +// +// The union-find and the two-scan structure are theirs. How a pixel's earlier neighbours are FOUND +// is not: see sparseccl below. #include @@ -18,16 +21,6 @@ void StrongPixelSet::AddStrongPixel(uint16_t col, uint16_t line, int32_t photons ++strong_pixel_count; } -bool is_far_enough(strong_pixel pixel0, strong_pixel pixel1) { - return (pixel1.line - pixel0.line) > 1; -} - -bool is_adjacent(strong_pixel pixel0, strong_pixel pixel1) { - auto line_diff = pixel0.line - pixel1.line; - auto col_diff = pixel0.col - pixel1.col; - return line_diff <= 1 && line_diff >= -1 && col_diff <= 1 && col_diff >= -1; -} - uint32_t StrongPixelSet::find_root(uint32_t e) { uint32_t r = e; while (L[r] != r) @@ -52,18 +45,38 @@ std::vector StrongPixelSet::sparseccl() { unsigned int labels = 0; - // first scan: pixel association - uint32_t start_j = 0; + // First scan: pixel association. The pixels arrive in raster order - line ascending, column + // ascending within a line - which upstream uses to walk a sliding window of the last two lines, + // testing every pixel in it for adjacency. That is quadratic in how many strong pixels a line + // pair holds: fine for the silicon-tracker hits it was written for, but a flooded detector line + // holds four thousand of them, and labelling a fully lit frame took 76 seconds. + // + // Since the columns ascend, the window need not be walked. A pixel's earlier 8-neighbours are + // exactly the one to its left and the at most three above it, so keep the previous line's range + // and a cursor into it that only ever moves forward - the same four neighbours the GPU extractor + // finds by binary search. Same edge set, same unions in the same order, therefore the same + // labels; the flooded frame now takes 0.16 s. + uint32_t line_begin = 0; // first pixel of the line being scanned + uint32_t prev_begin = 0, prev_end = 0; // the pixels of the line above it + uint32_t up = 0; // cursor into [prev_begin, prev_end) for (uint32_t i = 0; i < pixels.size(); ++i) { L[i] = i; - uint32_t ai = i; - for (uint32_t j = start_j; j < i; ++j) { - if (is_adjacent(pixels[i], pixels[j])) { - ai = make_union(ai, find_root(j)); - } else if (is_far_enough(pixels[j], pixels[i])) { - ++start_j; - } + if (i > 0 && pixels[i].line != pixels[i - 1].line) { + // The line above is the previous one only if it really is the line above: a line with no + // strong pixel at all leaves nothing to join to. + prev_begin = (pixels[i].line == pixels[i - 1].line + 1) ? line_begin : i; + prev_end = i; + line_begin = i; + up = prev_begin; } + uint32_t ai = i; + while (up < prev_end && pixels[up].col + 1 < pixels[i].col) + ++up; + for (uint32_t j = up; j < prev_end && pixels[j].col <= pixels[i].col + 1; ++j) + ai = make_union(ai, find_root(j)); + // The pixel to the left comes last, as it did when the window was walked in order. + if (i > line_begin && pixels[i - 1].col + 1 == pixels[i].col) + ai = make_union(ai, find_root(i - 1)); } // second scan: transitive closure @@ -85,19 +98,18 @@ std::vector StrongPixelSet::sparseccl() { void StrongPixelSet::FindComponentsImage(const SpotFindingSettings &settings, std::vector &spots) { - // Avoid spot finding, when more than 65536 strong pixel count (will be super slow) - if (!pixels.empty() && (strong_pixel_count < UINT16_MAX)) { - for (const auto &spot: sparseccl()) { - if (spot.PixelCount() <= settings.max_pix_per_spot) - spots.push_back(spot); - } + // No StrongPixelLimit test here: the caller knows how big the image is and has already applied it. + for (const auto &spot: sparseccl()) { + if (spot.PixelCount() <= settings.max_pix_per_spot) + spots.push_back(spot); } } void StrongPixelSet::FindSpots(const DiffractionExperiment &experiment, const SpotFindingSettings &settings, std::vector &spots, uint16_t module_number) { - // Avoid spot finding, when more than 65536 strong pixel count (will be super slow) - if (!pixels.empty() && (strong_pixel_count < UINT16_MAX)) { + // Per module, so the bar is the module's own - and ReadFPGAOutput has already refused anything + // past max_strong_pixel_per_module, far below it. + if (!pixels.empty() && (strong_pixel_count < StrongPixelLimit(RAW_MODULE_SIZE))) { for (const auto &spot: sparseccl()) { if ((spot.PixelCount() <= settings.max_pix_per_spot) && (spot.PixelCount() >= settings.min_pix_per_spot.value_or(2))) { diff --git a/image_analysis/spot_finding/StrongPixelSet.h b/image_analysis/spot_finding/StrongPixelSet.h index 2e53c28c2..5026122d3 100644 --- a/image_analysis/spot_finding/StrongPixelSet.h +++ b/image_analysis/spot_finding/StrongPixelSet.h @@ -38,6 +38,7 @@ public: // Every connected component of at most max-pix pixels. min-pix is deliberately NOT applied: it is // the only spot setting that changes between the passes of the per-image min-pix search, so one // connected-component search serves all of them (see ImageSpotFinder::ExtractComponents). + // The caller is responsible for not handing over more than StrongPixelLimit pixels. void FindComponentsImage(const SpotFindingSettings &settings, std::vector &spots); uint32_t GetStrongPixelCount() const; }; -- 2.54.0 From 1388b16d7a2a4d6985fcabcbc938679b703f5c05 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 10:41:34 +0200 Subject: [PATCH 004/179] Resolution estimate: predict past the edge of the detector The spot-finding resolution estimate was clamped so it could never beat the detector corner. On a crystal that diffracts past the corner that reports where the DETECTOR stops, which is the one thing this number is not for - it is meant to say how far a merge of data like these would reach, a property of the crystal and the exposure. The clamp also hid the interesting case: an estimate finer than what the run actually merged is the statement "this run was detector-limited", and there was no way to make it. The statistic already extrapolates. Its quantile sits in the middle of the fall-off, well inside what the detector records, so it goes on measuring the crystal's own decay when the detector cuts that decay short. Measured by truncating the spot lists of 31 battery crystals at an artificial detector edge and scoring the unclamped answer against each crystal's own measured CC1/2 = 0.30 crossing, it holds its 8-9% floor out to about 1.7x past the cut and only then drifts pessimistic, which is the safe direction. Every genuinely detector-limited crystal in the battery needs between 1.10x and 1.63x. Against a truth corrected for censoring - the six crystals whose merge is cut off by their own detector cannot have a measured crossing, so theirs is extrapolated from multiplicity-corrected and anchored on the 25 where both exist: symmetric-log RMS 13.5 -> 9.4% over 37 crystals, 25 -> 28 within 0.2 A. On the six detector-limited ones 26.1 -> 9.8% and the bias goes +19 -> -3%; on the 31 that are not, 9.23 -> 9.32%, i.e. it costs them nothing. The 0.30 tail fraction and the 2.25 reach were refit by leave-one-out against that truth and did not move. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H --- docs/CHANGELOG.md | 1 + docs/CPU_DATA_ANALYSIS.md | 2 +- image_analysis/spot_finding/SpotUtils.cpp | 12 ++++++------ image_analysis/spot_finding/SpotUtils.h | 8 +++++--- rugnux/ResultReport.cpp | 4 +++- tests/SpotUtilsTest.cpp | 7 +++++-- 6 files changed, 21 insertions(+), 13 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index a1dac33e0..ecbec9c6e 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,7 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* The spot-finding resolution estimate is no longer capped at the corner of the detector, so a crystal that diffracts past the edge is reported as reaching past it, and a run where the estimate is finer than what was merged is a detector-limited run. * Spot finding no longer throws away a whole image when it holds many strong pixels: the limit follows the detector (one pixel in 64) instead of standing at the 65535 that suited a 4-megapixel detector, which a strongly diffracting crystal on an 18-megapixel one passes on its best frames. * The connected-component search is linear in the strong pixels rather than quadratic in how many of them a detector line holds; a fully lit image is labelled in 0.16 s instead of 76 s, and the spots it finds are unchanged. * In `jfjoch_viewer`, Alt and the mouse wheel step through the dataset one image at a time. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index b53529b9e..29650cc68 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -389,7 +389,7 @@ is kept; the frame is then integrated once at that min-pix. The fraction factor ### 3.6 Predicting the resolution the merged data will reach -A per-image **resolution estimate** is read off the finished spot list. It predicts how far the *merged* data will reach, not how far the furthest spot on this image lies. Each non-ice spot is weighted by $\sqrt{I}$ — the intensity is a summed photon count, so $\sqrt{I}$ is its Poisson significance — the $1/d^2$ is found beyond which a fraction $f=0.30$ of that weight lies, and the estimate is that resolution taken $2.25\times$ further in $1/d$, clamped so it can never beat the corner of the detector. The dataset value is the median over images. +A per-image **resolution estimate** is read off the finished spot list. It predicts how far the *merged* data will reach, not how far the furthest spot on this image lies. Each non-ice spot is weighted by $\sqrt{I}$ — the intensity is a summed photon count, so $\sqrt{I}$ is its Poisson significance — the $1/d^2$ is found beyond which a fraction $f=0.30$ of that weight lies, and the estimate is that resolution taken $2.25\times$ further in $1/d$. It is deliberately **not** limited to what the detector records: the quantile sits in the middle of the fall-off, well inside the recorded range, so it goes on measuring the crystal where the detector stops before the diffraction does, and on such a run it reads finer than the detector corner. The dataset value is the median over images. Both constants carry a mechanism. A quantile from the middle of the distribution measures the *shape* of the fall-off, which is the crystal's own $\exp(-B/2d^{2})$, where the extreme end of it measures where detection stops — a threshold that moves with the exposure and with how many reflections the unit cell puts on a frame. And merging averages many observations of each reflection, so intensities go on being measurable a fixed factor in $1/d$ past the point at which one image's spot finder still detects them; that factor is the $2.25$. Both are calibrated on rotation data against the resolution at which per-shell CC1/2 falls through 0.30, and the estimate is good to about 0.2 Å there. It is a prediction and not a measurement of what a run achieved: nothing downstream is cut on it, and it is reported alone (rugnux `SPOT_RESOLUTION_ESTIMATE`, and per image in the stream, the plots and HDF5). diff --git a/image_analysis/spot_finding/SpotUtils.cpp b/image_analysis/spot_finding/SpotUtils.cpp index 04fe89661..433eb6610 100644 --- a/image_analysis/spot_finding/SpotUtils.cpp +++ b/image_analysis/spot_finding/SpotUtils.cpp @@ -129,7 +129,7 @@ namespace { constexpr size_t SPOT_RESOLUTION_MIN_SPOTS = 4; } -std::optional GetResolution(const std::vector &spots, float detector_d_min_A) { +std::optional GetResolution(const std::vector &spots) { // Each spot enters weighted by its own signal-to-noise. The intensity is a summed photon count, so // it is Poisson and its significance is sqrt(I): that keeps a marginal high-resolution detection // from counting for as much as a real reflection, without letting the handful of very strong @@ -161,10 +161,10 @@ std::optional GetResolution(const std::vector &spots, float d break; } - const float d_A = 1.0f / (SPOT_RESOLUTION_MERGE_REACH * std::sqrt(one_over_d2)); - - // However far the crystal diffracts, no merge reaches past the corner of the detector. - return detector_d_min_A > 0.0f ? std::max(d_A, detector_d_min_A) : d_A; + // Not clamped at the corner of the detector. The quantile is read from the middle of the + // fall-off, so it still measures the crystal where the detector cuts that fall-off short; + // clamping reported where the detector stops instead, which is the one thing this is not for. + return 1.0f / (SPOT_RESOLUTION_MERGE_REACH * std::sqrt(one_over_d2)); } void GenerateSpotPlot(DataMessage &msg, const std::vector &spots, float d_min_A) { @@ -232,7 +232,7 @@ void SpotAnalyze(const DiffractionExperiment &experiment, GenerateSpotPlot(output, spots_out, spot_d_min.value_or(0.0f) > 0 ? *spot_d_min : experiment.GetDetectorMaxResolution_A()); - output.resolution_estimate = GetResolution(spots_out, experiment.GetDetectorMaxResolution_A()); + output.resolution_estimate = GetResolution(spots_out); // One decision drives both: if indexing is to use the ice-band spots, the spot budget must not // throw them away before it gets the chance. diff --git a/image_analysis/spot_finding/SpotUtils.h b/image_analysis/spot_finding/SpotUtils.h index f0b1aa3d3..1eaa5f03e 100644 --- a/image_analysis/spot_finding/SpotUtils.h +++ b/image_analysis/spot_finding/SpotUtils.h @@ -37,9 +37,11 @@ void FilterSpuriousHighResolutionSpots(std::vector &spots, float thr // factor further in 1/d than the quantile, because averaging many observations goes on measuring // intensities that one image cannot detect. // -// detector_d_min_A is the corner of the detector, which the answer is never allowed to beat; pass 0 to -// leave it unclamped. Returns nothing when the image has too few spots to have a fall-off at all. -std::optional GetResolution(const std::vector &spots, float detector_d_min_A = 0.0f); +// The answer is deliberately NOT limited to what this detector records. The quantile sits in the +// middle of the fall-off, well inside the recorded range, so it goes on measuring the crystal when +// the detector stops before the diffraction does - which is the case the number is most wanted for. +// Returns nothing when the image has too few spots to have a fall-off at all. +std::optional GetResolution(const std::vector &spots); void SpotAnalyze(const DiffractionExperiment &experiment, const SpotFindingSettings &settings, diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 8b6ef13fb..876121f8d 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -102,7 +102,9 @@ std::string RenderResultReport(const std::string &output_prefix, << " How far the merged data are expected to reach, read off the found spots alone - no\n" << " lattice, no integration, no merge. It is a prediction, good to about 0.2 A on the\n" << " rotation data it was calibrated on, and it is not what the run achieved: compare it\n" - << " with INCLUDE_RESOLUTION_RANGE in section 5.\n"; + << " with INCLUDE_RESOLUTION_RANGE in section 5. It is not limited to what this detector\n" + << " records: where it reads finer than the high-resolution end of that range, the crystal\n" + << " diffracts past the corner and the run is detector-limited.\n"; } if (result.pass_count > 1) { diff --git a/tests/SpotUtilsTest.cpp b/tests/SpotUtilsTest.cpp index 17f70b701..0a804c9f0 100644 --- a/tests/SpotUtilsTest.cpp +++ b/tests/SpotUtilsTest.cpp @@ -39,8 +39,11 @@ TEST_CASE("GetResolution") { REQUIRE(d.has_value()); CHECK(*d == Catch::Approx(1.0 / (2.25 * std::sqrt(0.8))).epsilon(1e-4)); - // The merged data cannot beat the corner of the detector. - CHECK(*GetResolution(spots, 2.0f) == Catch::Approx(2.0)); + // The answer is not limited to what a detector records. Keeping only the five spots a detector + // reaching 1/d^2 = 0.5 would have recorded leaves the quantile at 0.4, and the estimate still + // extrapolates 2.25x past it instead of stopping at the cut. + const std::vector cut(spots.begin(), spots.begin() + 5); + CHECK(*GetResolution(cut) == Catch::Approx(1.0 / (2.25 * std::sqrt(0.4))).epsilon(1e-4)); // Ice-flagged spots take no part, however strong they are. std::vector with_ice = spots; -- 2.54.0 From 42add7c0e215e96c904252709ad9039ccb8afcc3 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 10:46:47 +0200 Subject: [PATCH 005/179] Report the strong-pixel count off the FPGA path too DataMessage::strong_pixel_count is written all the way to CBOR, to /entry/MX/strongPixels, to the receiver plots and to the frontend - and only the FPGA receiver ever set it. Every EIGER dataset, and every rugnux run, stored zeros. It is the one number that distinguishes an image the spot finder gave up on from an image that did not diffract, and its absence is why a defect that cost 767 of 1800 images their spots looked like a crystal that stopped diffracting for two blocks of the sweep. Fill it on the CPU/GPU path as well: the finders know the count, the GPU extractor already had it on the device, and it rides back with the spot count on the frame's one synchronisation, so it costs nothing. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H --- docs/CHANGELOG.md | 1 + image_analysis/MXAnalysisWithoutFPGA.cpp | 5 +++++ image_analysis/spot_finding/AdaptiveSpotFinderGPU.h | 1 + image_analysis/spot_finding/ImageSpotFinder.cpp | 1 + image_analysis/spot_finding/ImageSpotFinder.h | 6 ++++++ image_analysis/spot_finding/ImageSpotFinderGPU.h | 1 + image_analysis/spot_finding/SpotExtractorGPU.cu | 4 ++++ image_analysis/spot_finding/SpotExtractorGPU.h | 6 ++++++ 8 files changed, 25 insertions(+) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index ecbec9c6e..01fb728d7 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,7 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* `/entry/MX/strongPixels` and the strong-pixel plot are filled on the CPU/GPU analysis path, not only behind the FPGA, so an image the spot finder gave up on can be told from one that did not diffract. * The spot-finding resolution estimate is no longer capped at the corner of the detector, so a crystal that diffracts past the edge is reported as reaching past it, and a run where the estimate is finer than what was merged is a detector-limited run. * Spot finding no longer throws away a whole image when it holds many strong pixels: the limit follows the detector (one pixel in 64) instead of standing at the 65535 that suited a 4-megapixel detector, which a strongly diffracting crystal on an 18-megapixel one passes on its best frames. * The connected-component search is linear in the strong pixels rather than quadratic in how many of them a detector line holds; a fully lit image is labelled in 0.16 s instead of 76 s, and the spots it finds are unchanged. diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 0c964db70..57849fb15 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -251,6 +251,11 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, indexer.ProcessImage(output, spot_finding_settings, *prediction, integrate_fn); } + // Recorded whichever way the frame went. A frame holding StrongPixelLimit of them is given up + // on and reports no spots at all, and this is the only thing that tells such a frame from one + // that did not diffract. + output.strong_pixel_count = finder.StrongPixelCount(); + #ifdef JFJOCH_USE_CUDA if (fused) { // Lift the azimuthal profile the fused engine computed in the same detection pass; its azint diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h index 24e8cb74e..c15d48c0e 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h @@ -115,6 +115,7 @@ public: void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; void SetResolutionMaskBits(const std::vector &packed_mask) override; + [[nodiscard]] uint32_t StrongPixelCount() const override { return extractor.StrongPixelCount(); } const std::vector &ExtractComponents(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; diff --git a/image_analysis/spot_finding/ImageSpotFinder.cpp b/image_analysis/spot_finding/ImageSpotFinder.cpp index 8fe18f98a..7eab95ccf 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.cpp +++ b/image_analysis/spot_finding/ImageSpotFinder.cpp @@ -70,6 +70,7 @@ void ImageSpotFinder::ExtractComponentsHost(const ImagePreprocessorBuffer &image } } + strong_pixel_count = static_cast(strong_pixel.size()); components.clear(); // The connected-component search gives up on a frame with this many strong pixels, so their values // are of no use - not even worth gathering off the device. diff --git a/image_analysis/spot_finding/ImageSpotFinder.h b/image_analysis/spot_finding/ImageSpotFinder.h index 6566a84d4..14ec681bf 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.h +++ b/image_analysis/spot_finding/ImageSpotFinder.h @@ -15,6 +15,7 @@ class ImageSpotFinder { // value. Kept as members only to reuse the allocation from image to image. std::vector strong_pixel; std::vector strong_pixel_value; + uint32_t strong_pixel_count = 0; protected: const int32_t width, height; std::vector output_buffer; @@ -45,6 +46,11 @@ public: // components from those pixels. virtual void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) = 0; + // Strong pixels the last extraction saw, after the resolution mask. Reported whether or not the + // frame was given up on for holding StrongPixelLimit of them, which is the point of it: such a + // frame yields no spots at all, and without this nothing distinguishes it from a blank one. + [[nodiscard]] virtual uint32_t StrongPixelCount() const { return strong_pixel_count; } + // Peak-excluded per-ring background of the last Detect(), in the bins of the azimuthal-integration // mapping and in raw photon counts. Only the adaptive finders build one (it is what sets their // threshold); empty for everyone else, and for a frame with nothing valid to reduce. diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.h b/image_analysis/spot_finding/ImageSpotFinderGPU.h index df9e76621..d63f4acd1 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.h +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.h @@ -26,6 +26,7 @@ public: void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; void SetResolutionMaskBits(const std::vector &packed_mask) override; + [[nodiscard]] uint32_t StrongPixelCount() const override { return extractor.StrongPixelCount(); } const std::vector &ExtractComponents(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; }; diff --git a/image_analysis/spot_finding/SpotExtractorGPU.cu b/image_analysis/spot_finding/SpotExtractorGPU.cu index 2795de1b9..0cc6624b5 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.cu +++ b/image_analysis/spot_finding/SpotExtractorGPU.cu @@ -287,6 +287,7 @@ SpotExtractorGPU::SpotExtractorGPU(int32_t in_width, int32_t in_height, std::sha gpu_spot(max_strong), gpu_spot_out(max_strong), gpu_nspot(1), + host_nstrong(1), host_nspot(1), host_spot(SPOT_PREFIX) { // One block per few hundred words: enough blocks to fill the device, few enough that the serial @@ -336,6 +337,9 @@ void SpotExtractorGPU::Extract(const uint32_t *gpu_strong, const int32_t *gpu_im finish_components<<<1, FINISH_THREADS, 0, *stream>>>(gpu_index, gpu_value, gpu_root, gpu_label, gpu_count, gpu_spot, gpu_spot_out, gpu_nspot, gpu_nstrong, max_strong, width, max_pix); + // Rides along with the spot count on the frame's one synchronisation, so knowing how many strong + // pixels there were costs nothing. + cuda_err(cudaMemcpyAsync(host_nstrong, gpu_nstrong, sizeof(uint32_t), cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaMemcpyAsync(host_nspot, gpu_nspot, sizeof(uint32_t), cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaMemcpyAsync(host_spot, gpu_spot_out, SPOT_PREFIX * sizeof(SpotExtractorGPUSpot), cudaMemcpyDeviceToHost, *stream)); diff --git a/image_analysis/spot_finding/SpotExtractorGPU.h b/image_analysis/spot_finding/SpotExtractorGPU.h index 4cbc93735..415c69e94 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.h +++ b/image_analysis/spot_finding/SpotExtractorGPU.h @@ -75,6 +75,7 @@ class SpotExtractorGPU { CudaDevicePtr gpu_spot_out; CudaDevicePtr gpu_nspot; + CudaHostPtr host_nstrong; CudaHostPtr host_nspot; CudaHostPtr host_spot; // SPOT_PREFIX entries, pinned std::vector overflow_spot; // only for a frame with more spots than that @@ -89,4 +90,9 @@ public: // extractor would. void Extract(const uint32_t *gpu_strong, const int32_t *gpu_image, const SpotFindingSettings &settings, std::vector &spots); + + // Strong pixels the last Extract() saw, after the resolution mask. Reported whether or not the + // frame was given up on, which is the point of it: a frame at or above max_strong yields no spots + // at all, and this is what says so. + [[nodiscard]] uint32_t StrongPixelCount() const { return *host_nstrong.get(); } }; -- 2.54.0 From fe1f3ba96e75ae664ac2a843285c142a10106a8d Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 10:49:18 +0200 Subject: [PATCH 006/179] rugnux: report the geometry in the form the broker takes it in A rotation run post-refines the detector distance and the beam centre, and that is usually the best measurement of them anyone has - but the only way back into the instrument was to read two numbers off the report by eye and retype them into a collection. Write them once more, as the object jfjoch_broker takes: the four required properties of dataset_settings in broker/jfjoch_api.yaml, spelled the way the API spells them, on one line of valid JSON that a script can lift with a grep and POST. The existing DETECTOR_DISTANCE and BEAM_CENTRE keys are unchanged and stay the ones a person reads. REPORT_VERSION is not bumped: a new key breaks no consumer of the old ones. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H --- docs/CHANGELOG.md | 1 + docs/RUGNUX.md | 17 ++++++++++++++++- rugnux/ResultReport.cpp | 15 +++++++++++++++ 3 files changed, 32 insertions(+), 1 deletion(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 01fb728d7..b0d50b5f0 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,7 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* The rugnux results report carries `JFJOCH_DATASET_SETTINGS=`, the geometry the run integrated at written as one line of JSON in the form `jfjoch_broker` takes it, so a refined beam centre and distance can be carried back to the instrument. * `/entry/MX/strongPixels` and the strong-pixel plot are filled on the CPU/GPU analysis path, not only behind the FPGA, so an image the spot finder gave up on can be told from one that did not diffract. * The spot-finding resolution estimate is no longer capped at the corner of the detector, so a crystal that diffracts past the edge is reported as reaching past it, and a run where the estimate is finer than what was merged is a detector-limited run. * Spot finding no longer throws away a whole image when it holds many strong pixels: the limit follows the detector (one pixel in 64) instead of standing at the 65535 that suited a 4-megapixel detector, which a strongly diffracting crystal on an 18-megapixel one passes on its best frames. diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index c8ed84f1d..77508371e 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -564,7 +564,22 @@ keep their numbers. **`SPOT_RESOLUTION_ESTIMATE=`** in section 1 is how far the merged data are expected to reach, read off the found spots alone — no lattice, no integration, no merge — so it is there on a run that never merges, and on a run that does it can be read against `INCLUDE_RESOLUTION_RANGE` in section 5. It is a -prediction, good to about 0.2 Å on rotation data; nothing is cut on it. +prediction, good to about 0.2 Å on rotation data; nothing is cut on it. It is **not** limited to what +the detector records: where it reads finer than the high-resolution end of `INCLUDE_RESOLUTION_RANGE`, +the crystal diffracts past the corner and the run was detector-limited. + +**`JFJOCH_DATASET_SETTINGS=`** in section 1 is the geometry the run integrated at — on a rotation run +the post-refined one — written as the object `jfjoch_broker` takes it in: the four required properties +of `dataset_settings` in `broker/jfjoch_api.yaml`, on one line of valid JSON, so a refined beam centre +and distance can go back to the instrument for the next collection without anyone retyping them. + +``` +JFJOCH_DATASET_SETTINGS= {"beam_x_pxl": 2078.24, "beam_y_pxl": 2233.92, "detector_distance_mm": 190.311, "incident_energy_keV": 12.4000} +``` + +```bash +grep '^JFJOCH_DATASET_SETTINGS=' out_report.txt | cut -d' ' -f2- > geometry.json +``` **Which pass.** A rotation run integrates twice — once at the geometry in the input file, then again at the post-refined geometry — and can integrate a third time if a guard rejects the second pass. diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 876121f8d..07822fe5a 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -95,6 +95,21 @@ std::string RenderResultReport(const std::string &output_prefix, << " The distance and beam centre above are the ones this result was integrated at, which on\n" << " a rotation run is the post-refined geometry rather than the values in the input file.\n"; + // The same geometry once more, as the object jfjoch_broker takes it in: the four required + // properties of dataset_settings in broker/jfjoch_api.yaml, spelled the way the API spells them. + // A run that refined the geometry is usually the best measurement of it anyone has, and without + // this the only way back into the instrument is to read two numbers off this report by eye and + // retype them. One line, valid JSON, so a script can lift it with a grep and POST it. + os << "\n"; + Key(os, "JFJOCH_DATASET_SETTINGS", + fmt::format(R"({{"beam_x_pxl": {:.2f}, "beam_y_pxl": {:.2f}, "detector_distance_mm": {:.3f}, )" + R"("incident_energy_keV": {:.4f}}})", + result.used_beam_x_pxl, result.used_beam_y_pxl, result.used_distance_mm, + experiment.GetDatasetSettings().GetPhotonEnergy_keV())); + os << "\n" + << " The geometry above as jfjoch_broker's dataset_settings, to carry a refined beam centre and\n" + << " distance back to the instrument for the next collection.\n"; + if (result.spot_resolution_estimate_A.has_value()) { os << "\n"; Key(os, "SPOT_RESOLUTION_ESTIMATE", fmt::format("{:.2f}", *result.spot_resolution_estimate_A)); -- 2.54.0 From 6516bc96af2c88a286736da092b3018d8b8e0060 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 11:27:28 +0200 Subject: [PATCH 007/179] Ice rings: carry the list past 1.5 A, where ice does not stop The eleven measured bands end at 1.522 A because their source says so in its own words - "pure hexagonal ice has 11 diffraction rings between 4 and 1.5 A resolution" - and its subject was detecting ice in deposited data, not masking it. On a detector that reaches further, the rings it does not list are the ones left in the data: on the strong rotation set just added to the battery, 44% of every image's spots sit in ice bands, and beyond 1.5 A the spot list is ice and nothing else, which is why the resolution estimate read the ice rather than the crystal. There is nothing measured to copy below 1.522 A, so the eight added bands are calculated. Enumerating hkl is not enough and the code already said so: ice Ih is P6_3/mmc with O on 4f, and most of what enumeration emits is extinguished by the OXYGEN SUBLATTICE rather than by the space group - which is why (004) at 1.830 A and (104) at 1.657 A are missing from the measured list although they sit inside its range and its reflection conditions allow them (for (00l) the structure factor goes as cos(2*pi*l*z), and z ~ 1/16 kills l = 4). So compute structure factors - oxygen only, the hydrogens being half-occupancy disordered and weak to X-rays - and keep the lines reaching 3% of the strongest. That rule REPRODUCES THE MEASURED ELEVEN EXACTLY and every line it drops inside their range computes to zero, which is what makes it trustworthy below 1.522 A. It stops at 1.170 A: below that the real lines fall to 2-3% while the extinct ones rise to about 1%, and an oxygen-only calculation cannot separate them honestly. Every added band was independently confirmed in the data - the spot-count histogram of the strong set peaks at each of them and is empty between - and every line the rule calls extinct is absent there too. Costs, measured. The bands are inert above 1.6 A: on 38 of 39 battery sets the profile ice score does not move at all, and the one that appeared to (a jet set, 1.25 -> 2.60) does not on the peak-excluded profile the score actually uses - that was Bragg peaks in the plain profile, which is what the peak exclusion is for. Where a detector does reach past 1.5 A the bands cover more of reciprocal space: unchanged at 1.6 A, +7.4 points at 1.4 A, +16.5 at 1.18 A. On the strong set that is 17% -> 27% of reflections held out of the scale fit, and it shows: the spot resolution estimate improves from 1.33 to 1.46 A against a truth near 1.42, while CC1/2 falls 98.5 -> 97.6% and ISa 3.58 -> 3.37. Ice handling only runs at all on a run that trips the ice gate, so a clean crystal pays nothing. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H --- common/Definitions.h | 26 ++++++++++++++++++--- docs/ACKNOWLEDGEMENT.md | 15 ++++++++---- docs/CHANGELOG.md | 1 + docs/CPU_DATA_ANALYSIS.md | 12 +++++++--- image_analysis/geom_refinement/Calibrants.h | 5 ++-- tests/CalibrationTest.cpp | 2 +- 6 files changed, 48 insertions(+), 13 deletions(-) diff --git a/common/Definitions.h b/common/Definitions.h index 08de6e7f8..140d48836 100644 --- a/common/Definitions.h +++ b/common/Definitions.h @@ -57,11 +57,31 @@ constexpr size_t PEDESTAL_G0_WRONG_GAIN_ALLOWED_COUNT = 2; constexpr size_t MESSAGE_SIZE_FOR_START_END = (256*1024*1024); // pessimistic highest value constexpr float LAB6_CELL_A = 4.156468f; -// Ice ring resolution taken from: +// Hexagonal-ice ring positions. The first eleven are the measured ones from: // Moreau, Atakisi, Thorne, Acta Cryst D77, 2021, 540,554 // https://journals.iucr.org/d/issues/2021/04/00/tz5104/index.html - -constexpr std::array ICE_RING_RES_A = {3.895, 3.661, 3.438, 2.667, 2.249, 2.068, 1.947, 1.916, 1.882, 1.719, 1.522}; +// That paper's Table 1 stops at 1.522 A because it says so in its own words - "pure hexagonal ice has +// 11 diffraction rings between 4 and 1.5 A resolution" - and its subject was detecting ice in the PDB, +// not masking it. Ice does not stop there, and on a detector that reaches past 1.5 A the rings it does +// not list are the ones left in the data: one strongly diffracting set had 44% of every frame's spots +// in ice bands, and beyond 1.5 A its spots were ice and nothing else. +// +// The rest are calculated, because past that paper there is nothing measured to copy. Ice Ih is +// P6_3/mmc with O on 4f, so enumerating hkl is not enough - most of what it emits is extinguished by +// the oxygen sublattice rather than by the space group, which is why (004) at 1.830 A and (104) at +// 1.657 A are absent from the measured list although they lie inside its range and its reflection +// conditions allow them. Structure factors were computed instead (oxygen only - the hydrogens are +// half-occupancy disordered and scatter X-rays weakly) and the lines kept are those reaching 3% of the +// strongest. That reproduces the measured eleven exactly, and every line it drops in their range +// computes to zero, which is what makes the same rule trustworthy below 1.522 A. +// Following Roettger, Endriss, Ihringer, Doyle & Kuhs (1994) Acta Cryst. B50, 644-648 for the cell. +// +// It stops at 1.170 A: below that the calculated real lines fall to 2-3% while the extinct ones rise +// to about 1%, and an oxygen-only calculation cannot separate them honestly any further. +constexpr std::array ICE_RING_RES_A = {3.895, 3.661, 3.438, 2.667, 2.249, 2.068, 1.947, + 1.916, 1.882, 1.719, 1.522, + 1.472, 1.443, 1.371, 1.366, 1.298, 1.261, 1.224, + 1.170}; // True when resolution d (Angstrom) sits within half_width of a hexagonal-ice powder ring, in the // q = 2*pi/d units the spot-finder uses (ice_ring_width_Q_recipA). Used to drop ice-contaminated diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index a09ebc6dc..3ee56f9df 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -98,12 +98,19 @@ structure-factor machinery, and MTZ / XDS_ASCII I/O. Vendored in `gemmi_gph/`, s licence obligation. M. Wojdyr, "GEMMI: A library for structural biology" (2022), J. Open Source Softw. 7, 4200 [doi:10.21105/joss.04200](https://doi.org/10.21105/joss.04200). -**Hexagonal-ice ring positions** — the eleven measured ring $d$ spacings the ice-ring score, the -ice-ring flagging and the ice calibrant are all built on are taken from the measurements of, not -enumerated from a cell. D. W. Moreau, H. Atakisi and R. E. Thorne, "Ice in biomolecular -cryocrystallography" (2021), Acta Cryst. D77, 540-554 +**Hexagonal-ice ring positions** — the eleven ring $d$ spacings from 3.895 to 1.522 Å that the +ice-ring score, the ice-ring flagging and the ice calibrant are all built on are taken from the +measurements of, not enumerated from a cell. D. W. Moreau, H. Atakisi and R. E. Thorne, "Ice in +biomolecular cryocrystallography" (2021), Acta Cryst. D77, 540-554 [doi:10.1107/S2059798321001170](https://doi.org/10.1107/S2059798321001170). +That list ends at 1.522 Å by its own scope, so the eight bands below it are calculated here rather +than taken from anyone: ice Ih structure factors on the oxygen sublattice, kept where they reach 3% of +the strongest line, which reproduces the eleven measured positions exactly. The lattice constants are +Röttger and co-workers'. A. Röttger, A. Endriss, J. Ihringer, S. Doyle and W. F. Kuhs, "Lattice +constants and thermal expansion of H2O and D2O ice Ih between 10 and 265 K" (1994), Acta Cryst. B50, +644-648 [doi:10.1107/S0108768194004933](https://doi.org/10.1107/S0108768194004933). + **Diffraction anisotropy** — the description of the overall fall-off by a single anisotropic displacement tensor, its symmetry constraints, and the fact that only its deviatoric part is determined (the isotropic part being degenerate with the overall scale) are Sheriff and Hendrickson's. diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index b0d50b5f0..e6074e4cf 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,7 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* Eight further hexagonal-ice bands between 1.472 and 1.170 A are recognised. The measured list they extend stops at 1.522 A because its source went no further, not because ice does, and on a detector reaching past 1.5 A the unlisted rings were left in the data - on a strongly diffracting set they were 44% of every image's spots. Ice handling still only applies where a run trips the ice gate. * The rugnux results report carries `JFJOCH_DATASET_SETTINGS=`, the geometry the run integrated at written as one line of JSON in the form `jfjoch_broker` takes it, so a refined beam centre and distance can be carried back to the instrument. * `/entry/MX/strongPixels` and the strong-pixel plot are filled on the CPU/GPU analysis path, not only behind the FPGA, so an image the spot finder gave up on can be told from one that did not diffract. * The spot-finding resolution estimate is no longer capped at the corner of the detector, so a crystal that diffracts past the edge is reported as reaching past it, and a run where the estimate is finer than what was merged is a detector-limited run. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 29650cc68..ece6feb53 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -344,13 +344,19 @@ Because detection reads the pixel's ring, a pixel that falls outside the azimuth Spot finding can be restricted to a resolution range $[d_\mathrm{high}, d_\mathrm{low}]$ by masking pixels outside the range. Optionally, spots in identified ice-ring regions can be tagged so that subsequent indexing/refinement may include or exclude them (see §4 and §6). -A single per-image **ice-ring score** is derived from a radial profile: for each hexagonal-ice powder ring (positions $d$ from Moreau *et al.*, Acta Cryst D77, 2021), the profile intensity at the ring is divided by a smooth background estimated from the *whole* profile — a running median of the non-ice bins, interpolated under each ring — and the strongest ring's ratio is reported (1 = no ice, $>1$ = ice above background). A whole-profile background is used rather than a couple of adjacent shoulder bins so the estimate is robust to the radial binning: at a coarse Q-spacing a local shoulder can be only ~1 bin and would double-count the ring's own edge (offline processing defaults to a fine 0.01 1/Å spacing, `--azim-q-spacing`, so the rings are well resolved). The reported quantity is the ice *magnitude* rather than a significance: with many photons any real ice ring is statistically significant, so significance does not discriminate. +A single per-image **ice-ring score** is derived from a radial profile: for each hexagonal-ice powder ring (see *Where the ring positions come from* below), the profile intensity at the ring is divided by a smooth background estimated from the *whole* profile — a running median of the non-ice bins, interpolated under each ring — and the strongest ring's ratio is reported (1 = no ice, $>1$ = ice above background). A whole-profile background is used rather than a couple of adjacent shoulder bins so the estimate is robust to the radial binning: at a coarse Q-spacing a local shoulder can be only ~1 bin and would double-count the ring's own edge (offline processing defaults to a fine 0.01 1/Å spacing, `--azim-q-spacing`, so the rings are well resolved). The reported quantity is the ice *magnitude* rather than a significance: with many photons any real ice ring is statistically significant, so significance does not discriminate. The profile the score is read off is the **peak-excluded** one, not the plain azimuthal integration: where adaptive spot finding runs (§3.2 — the offline and viewer default), the score uses the sigma-clipped per-resolution-ring background that finder already computes for its threshold. This matters more than it sounds. A plain azimuthal profile is a per-ring *mean*, so a few strong low-resolution reflections landing in a ring's bin raise it exactly as ice would; measured over a rotation battery, that alone ranked ice-free crystals above crystals that really are iced. An ice ring is azimuthally smooth and survives the sigma clip, while Bragg peaks do not, so on the clipped profile ice-free crystals sit near 1 and crystals with confirmed ice above 2. Only where no adaptive finder ran (the FPGA workflow) does the score fall back to the plain profile. The radial profile sees ice only as a **smooth powder ring**. Ice in large crystallites diffracts as discrete spots, leaves the profile flat, and is invisible to the score above, so a second channel is read from the spot list itself: the spots found on the ice bands are counted against the spots found in the ice-free flanks $[w,2w)$ either side of each band, rescaled to the bands' own $q$ width (a flank landing on another ring is dropped with its width). The indicator is the ratio pooled over the run — per image the control is a handful of spots and the ratio means nothing — and it is taken before the spot-count filter, which orders ice spots last and would discard them first. The two channels barely overlap: smooth ice reads high on the profile and ~1 on the spots, textured ice the reverse, and a clean crystal ~1 on both. Both counts are stored per image (`spot_count_ice_rings`, `spot_count_ice_control`; HDF5 `/entry/MX/peakCountIceRingRes`, `/entry/MX/peakCountIceRingControl`). -Both channels are used offline as a **gate** on ice handling: unless the run reaches `--ice-min-score` (default 1.5) on the profile or `--ice-min-spot-ratio` (default 2.0) on the spots — 0 disables a channel — ice-ring flagging and the exclusion from the scale fit (§10.10) are skipped. The eleven fixed bands cover 16–26 % of the unique reflections at typical resolutions whether or not the crystal has ice, so handling ice on a clean crystal is a pure loss. +Both channels are used offline as a **gate** on ice handling: unless the run reaches `--ice-min-score` (default 1.5) on the profile or `--ice-min-spot-ratio` (default 2.0) on the spots — 0 disables a channel — ice-ring flagging and the exclusion from the scale fit (§10.10) are skipped. The fixed bands cover 16–26 % of the unique reflections at typical resolutions whether or not the crystal has ice, and more than that on a detector reaching past 1.5 Å, so handling ice on a clean crystal is a pure loss. + +**Where the ring positions come from.** The eleven bands from 3.895 to 1.522 Å are the *measured* positions of Moreau, Atakisi & Thorne (Acta Cryst D77, 2021, 540–554). That list ends at 1.522 Å by its own scope — the paper states that hexagonal ice "has 11 diffraction rings between 4 and 1.5 Å resolution", and its subject was detecting ice in deposited data rather than masking it — not because ice stops there. On a detector that reaches further, the rings it does not list are the ones left in the data. + +The eight bands below it are **calculated**, since past that paper there is nothing measured to copy. Enumerating $hkl$ from the ice Ih cell is not enough: ice Ih is $P6_3/mmc$ with oxygen on $4f$, and most of what enumeration emits is extinguished by the *oxygen sublattice* rather than by the space group — which is why (004) at 1.830 Å and (104) at 1.657 Å are missing from the measured list even though they lie inside its range and its reflection conditions allow them (for $(00l)$ the structure factor goes as $\cos 2\pi l z$, and $z \approx 1/16$ kills $l=4$). Structure factors are computed instead — oxygen only, the hydrogens being half-occupancy disordered and weak to X-rays — and the lines kept are those reaching 3 % of the strongest. That rule **reproduces the measured eleven exactly**, and every line it drops inside their range computes to zero, which is what makes it trustworthy below 1.522 Å. The cell is that of Röttger *et al.* (Acta Cryst B50, 1994, 644–648). The list stops at 1.170 Å because below it the calculated real lines fall to 2–3 % while the extinct ones rise to about 1 %, and an oxygen-only calculation cannot separate them any further. + +One consequence is worth stating: the profile score is the **strongest** ring's ratio, a maximum over the bands, so a longer list can only raise it. The gate at 1.5 is therefore read against a list of this length, and lengthening it again would need the gate re-checked. A further optional safeguard removes isolated high-resolution “spur” spots by detecting large gaps in $1/d$ (or $q$) space and discarding spots beyond the gap. This is intended for macromolecular diffraction where edge-of-detector backgrounds can be extremely low. @@ -568,7 +574,7 @@ Everything above fits the geometry to *Bragg* data, where the beam centre is the The ring positions are matched to the observed rings and the geometry is refined (Ceres, five parameters: beam centre, distance, and the two detector tilts) so that the $|s|$ predicted at each observed ring point matches the ring it belongs to. -**Calibrants.** LaB₆, silver behenate, CeO₂ and silicon are held as unit cells and their rings enumerated from them. Ice is held as the eleven **measured** hexagonal-ice ring positions of §3.3 instead, because hexagonal ice is $P6_3/mmc$ and enumerating $hkl$ from its cell would emit rings that are systematically absent. A calibrant is therefore a list of ring $q$ values throughout, not a cell. +**Calibrants.** LaB₆, silver behenate, CeO₂ and silicon are held as unit cells and their rings enumerated from them. Ice is held as the hexagonal-ice ring positions of §3.3 instead — measured to 1.522 Å, calculated below it — because hexagonal ice is $P6_3/mmc$ with oxygen on $4f$ and enumerating $hkl$ from its cell would emit rings the oxygen sublattice extinguishes. A calibrant is therefore a list of ring $q$ values throughout, not a cell. **What a ring can and cannot determine.** A ring is a conic centred on the beam, so a wrong centre makes its apparent radius oscillate once per turn, $r(\phi)=R+\delta_x\cos\phi+\delta_y\sin\phi$, with the **same amplitude on every ring**. A detector tilt $\beta$ produces a $\cos\phi$ term too — not the $\cos2\phi$ one might expect — but one that grows as the ring's radius *squared*, $r(\phi)=R+(R^2/F)(\beta_x\cos\phi+\beta_y\sin\phi)$; the true $\cos2\phi$ term is $O(R^3\beta^2/F^2)$, hundredths of a pixel. The two are therefore separated by how the amplitude scales with radius, which needs **at least two rings** — on a single ring they are exactly degenerate. None of this uses the calibrant's $d$-spacings, so the centre is determined without assuming anything about the standard. diff --git a/image_analysis/geom_refinement/Calibrants.h b/image_analysis/geom_refinement/Calibrants.h index 899cdae3d..7f1f0b9f7 100644 --- a/image_analysis/geom_refinement/Calibrants.h +++ b/image_analysis/geom_refinement/Calibrants.h @@ -26,8 +26,9 @@ const std::vector &Calibrants(); // RingOptimizerInput::q_expected expects - sorted ascending. Empty if the name is not a calibrant. // // A ring LIST rather than a UnitCell because ice cannot be given as one: hexagonal ice is P6_3/mmc, so -// enumerating hkl from its cell would emit rings that are systematically absent, while the eleven -// entries of ICE_RING_RES_A are measured ring positions. +// enumerating hkl from its cell would emit rings the oxygen sublattice extinguishes, while the entries +// of ICE_RING_RES_A are ring positions that carry intensity - measured to 1.522 A, and calculated with +// that extinction taken out below it. std::vector CalibrantRings(const std::string &name); // "lab6, agbh, ceo2, si, ice" - for the error message when a name is not one of them. diff --git a/tests/CalibrationTest.cpp b/tests/CalibrationTest.cpp index 00cb8f7dc..135f3f62d 100644 --- a/tests/CalibrationTest.cpp +++ b/tests/CalibrationTest.cpp @@ -65,7 +65,7 @@ TEST_CASE("Calibrants_CentredStandardsOmitTheExtinctRings", "[DetGeomCalib]") { // Ice is the reason the calibrant abstraction is a ring list and not a UnitCell: its entries are // measured ring positions, and enumerating hkl from the hexagonal cell would add rings that are // systematically absent in P6_3/mmc. -TEST_CASE("Calibrants_IceIsTheMeasuredRingList", "[DetGeomCalib]") { +TEST_CASE("Calibrants_IceIsTheRingList", "[DetGeomCalib]") { const auto q = CalibrantRings("ice"); REQUIRE(q.size() == ICE_RING_RES_A.size()); CHECK(std::is_sorted(q.begin(), q.end())); -- 2.54.0 From 0af109b838a49f2a60ac439c1a14c8ff67bb640c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 11:49:49 +0200 Subject: [PATCH 008/179] Spot finding: bound how large a spot may be, not how bright A component of more than 50 connected pixels was discarded. Under the self-calibrating threshold that is an INTENSITY CEILING, not a size bound: the threshold is an absolute per-ring contour, so a component's area above it grows as sigma^2*ln(A/T), without bound in the peak amplitude. Measured on the strong rotation set, footprints run 3 px at 30-100 counts to 50 px above 10000 - a slope of 4.4 px per ln(peak) against the fixed local-box test's 0.8, which saturates at 13 px and never reaches the bound at all. So the brighter a reflection, the more certainly it was thrown away: 59% of the box finder's d>3 A spots were missing from the adaptive finder's list, including every one of its ten strongest, at an intensity ratio of 1.085 for those that did match. Indexed spots per image collapsed from 220 to 7. Raise the bound to 200 - CrystFEL peakfinder8's --max-pix-count, the only directly comparable number in the field; XDS has no such parameter and guards on shape instead - and ask a component above 50 pixels to be COMPACT: it must fill a fifth of the square its bounding box fits inside. A Bragg reflection is round and fills about half of that square however bright it is; an ice arc, a cosmic-ray track or a lit detector row fills a fifth or less, and those are what an upper bound was ever protecting against. Below 50 nothing is asked of the shape, so every component accepted before still is. Integer arithmetic on both sides, and on the GPU the bounding side fits in what was padding, so the device struct does not grow. The shape test is what makes the raise safe. With a flat 200 alone, two battery crystals moved: one lost a little I/sigma, and the other's de-novo lattice was NOT MONOTONE in the bound - correct at 50, 100 and 200, wrong at 150 and at 250 and above - so 200 was partly luck. With the shape test both are unchanged to three decimals. De novo at defaults on the strong set: indexing rate 0.574 -> 0.804, completeness 97.3 -> 99.6%, 3.42 -> 6.07, R_meas 0.299 -> 0.249, CC1/2 0.947 -> 0.961, ISa 3.29 -> 4.13. The full 37-crystal battery, scored per shell, is still owed. Also documents what StrongPixelSet::AddStrongPixel has required since the component search became linear - pixels in raster order - and puts the existing test's insertion order into it. Both callers scan a bitmap in ascending flat index and always satisfied it; the test did not, and was the only thing that did not. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H --- docs/CHANGELOG.md | 2 + docs/CPU_DATA_ANALYSIS.md | 2 + etc/broker_local.json | 2 +- .../spot_finding/SpotExtractorGPU.cu | 25 ++++++++-- .../spot_finding/SpotExtractorGPU.h | 4 +- .../spot_finding/SpotFindingSettings.h | 24 +++++++++- .../spot_finding/StrongPixelSet.cpp | 41 +++++++++++----- image_analysis/spot_finding/StrongPixelSet.h | 9 +++- tests/StrongPixelSetTest.cpp | 47 +++++++++++++++---- 9 files changed, 127 insertions(+), 29 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index e6074e4cf..b2940fb9c 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,8 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* A spot may hold up to 200 connected pixels rather than 50. Under the self-calibrating threshold a reflection's footprint grows with its brightness, so the old bound acted as an intensity ceiling and discarded the strongest reflections of a strongly diffracting crystal - on one such set, every one of the ten brightest on an image. +* A spot larger than 50 pixels must also be compact, filling a fifth of the square its bounding box fits inside. An ice arc, a cosmic-ray track or a lit detector row is refused however many pixels it has, which is what the size bound used to do and now does without capping brightness. * Eight further hexagonal-ice bands between 1.472 and 1.170 A are recognised. The measured list they extend stops at 1.522 A because its source went no further, not because ice does, and on a detector reaching past 1.5 A the unlisted rings were left in the data - on a strongly diffracting set they were 44% of every image's spots. Ice handling still only applies where a run trips the ice gate. * The rugnux results report carries `JFJOCH_DATASET_SETTINGS=`, the geometry the run integrated at written as one line of JSON in the form `jfjoch_broker` takes it, so a refined beam centre and distance can be carried back to the instrument. * `/entry/MX/strongPixels` and the strong-pixel plot are filled on the CPU/GPU analysis path, not only behind the FPGA, so an image the spot finder gave up on can be told from one that did not diffract. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index ece6feb53..256cb58fe 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -372,6 +372,8 @@ Strong pixels are grouped into connected components (adjacent strong pixels) usi Spot-level filters include minimum/maximum pixel count and resolution limits. +**The upper bound is on how large a *round* spot may be, not on how bright one may be.** Under the self-calibrating threshold (§3.2) a component's area above the contour grows as $\sigma^2\ln(A/T)$ — without bound in the peak amplitude $A$ — so an upper bound on pixel count alone becomes an *intensity ceiling*: measured on a strongly diffracting crystal, footprints run 3 px at 30–100 counts to 50 px above 10 000, four times the slope the fixed local-box test gives, and the old bound of 50 discarded the brightest reflections on every image. The bound is therefore 200 (the same as CrystFEL `peakfinder8`'s `--max-pix-count`; XDS has no such parameter at all and guards on shape instead, with `SPOT_MAXIMUM-CENTROID`), and a component above 50 pixels must in addition **fill at least a fifth of the square its bounding box fits inside**. A Bragg reflection is round and fills about half of that square however bright it is; an ice arc, a cosmic-ray track or a lit detector row fills a fifth or less, which is what an upper bound was ever protecting against. Below 50 pixels no shape is asked for, so everything accepted before still is. The test is integer arithmetic, so the host and the GPU extractor agree by construction. + The host implementation (`StrongPixelSet::sparseccl`) is the SparseCCL of the ACTS/traccc project: it runs over the strong pixels sorted row-major, uses a sliding window over the previous line and a union-find whose root is each component's lowest index. On the GPU the same labelling runs **on the diff --git a/etc/broker_local.json b/etc/broker_local.json index 417b14e4a..ce9845eee 100644 --- a/etc/broker_local.json +++ b/etc/broker_local.json @@ -20,7 +20,7 @@ "signal_to_noise_threshold":4.0, "photon_count_threshold":10, "low_resolution_limit":50.0, - "max_pix_per_spot":50, + "max_pix_per_spot":200, "min_pix_per_spot":2, "indexing":true, "quick_integration":true diff --git a/image_analysis/spot_finding/SpotExtractorGPU.cu b/image_analysis/spot_finding/SpotExtractorGPU.cu index 0cc6624b5..d286c5586 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.cu +++ b/image_analysis/spot_finding/SpotExtractorGPU.cu @@ -181,7 +181,7 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 int32_t *__restrict__ count, SpotExtractorGPUSpot *__restrict__ scratch, SpotExtractorGPUSpot *__restrict__ out, uint32_t *__restrict__ nout, const uint32_t *__restrict__ nstrong, uint32_t capacity, - int width, int max_pix) { + int width, int max_pix, int shape_free_pix, int min_fill_percent) { const int n = static_cast(min(*nstrong, capacity)); if (threadIdx.x == 0) *nout = 0; __syncthreads(); @@ -231,10 +231,14 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 if (want > max_pix) continue; long long x = 0, y = 0; long long photons = 0, max_photons = LLONG_MIN; + int min_col = INT_MAX, max_col = INT_MIN, min_line = INT_MAX, max_line = INT_MIN; int found = 0; for (int j = i; j < n && found < want; j++) { if (root[j] != static_cast(i)) continue; const long long counts = value[j]; + const int col = static_cast(index[j] % width), line = static_cast(index[j] / width); + min_col = min(min_col, col); max_col = max(max_col, col); + min_line = min(min_line, line); max_line = max(max_line, line); // Integers, exactly as DiffractionSpot::AddPixel does them, so host and device agree by // construction - no rounding mode to match and nothing for either compiler to contract. x += static_cast(index[j] % width) * counts; @@ -247,14 +251,23 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 scratch[l].y = y; scratch[l].photons = photons; scratch[l].max_photons = max_photons; + scratch[l].bbox_side = max(max_col - min_col, max_line - min_line) + 1; } __syncthreads(); - // 4) max-pix filter, compacted by another prefix sum so the surviving spots keep their order + // 4) size and shape filter, compacted by another prefix sum so the surviving spots keep their + // order. The test is SpotShapeAccepted written out - the constants come in as arguments rather + // than being included here, so there is exactly one definition of them. + auto keep = [&](const SpotExtractorGPUSpot &s) { + if (s.pixel_count > max_pix) return false; + if (s.pixel_count <= shape_free_pix) return true; + return static_cast(s.pixel_count) * 100 + >= static_cast(min_fill_percent) * s.bbox_side * s.bbox_side; + }; const int label_chunk = (nlabel + nthreads - 1) / nthreads; const int label_lo = min(t * label_chunk, nlabel), label_hi = min(label_lo + label_chunk, nlabel); uint32_t nkeep = 0; - for (int i = label_lo; i < label_hi; i++) nkeep += (scratch[i].pixel_count <= max_pix) ? 1u : 0u; + for (int i = label_lo; i < label_hi; i++) nkeep += keep(scratch[i]) ? 1u : 0u; shared[t] = nkeep; __syncthreads(); for (int d = 1; d < nthreads; d <<= 1) { @@ -265,7 +278,7 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 } uint32_t pos = shared[t] - nkeep; for (int i = label_lo; i < label_hi; i++) - if (scratch[i].pixel_count <= max_pix) out[pos++] = scratch[i]; + if (keep(scratch[i])) out[pos++] = scratch[i]; if (t == nthreads - 1) *nout = shared[t]; } @@ -336,7 +349,9 @@ void SpotExtractorGPU::Extract(const uint32_t *gpu_strong, const int32_t *gpu_im resolve_roots<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_root, gpu_nstrong, max_strong); finish_components<<<1, FINISH_THREADS, 0, *stream>>>(gpu_index, gpu_value, gpu_root, gpu_label, gpu_count, gpu_spot, gpu_spot_out, gpu_nspot, - gpu_nstrong, max_strong, width, max_pix); + gpu_nstrong, max_strong, width, max_pix, + static_cast(SPOT_SHAPE_FREE_PIXELS), + static_cast(SPOT_MIN_FILL_PERCENT)); // Rides along with the spot count on the frame's one synchronisation, so knowing how many strong // pixels there were costs nothing. cuda_err(cudaMemcpyAsync(host_nstrong, gpu_nstrong, sizeof(uint32_t), cudaMemcpyDeviceToHost, *stream)); diff --git a/image_analysis/spot_finding/SpotExtractorGPU.h b/image_analysis/spot_finding/SpotExtractorGPU.h index 415c69e94..4bef4f31d 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.h +++ b/image_analysis/spot_finding/SpotExtractorGPU.h @@ -43,7 +43,9 @@ struct SpotExtractorGPUSpot { int64_t photons; int64_t max_photons; int32_t pixel_count; - int32_t padding; + // Longer side of the component's bounding box, for SpotShapeAccepted. It fits in what used to be + // padding, so carrying it costs nothing. + int32_t bbox_side; }; class SpotExtractorGPU { diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index a2bdde060..d610277bc 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/image_analysis/spot_finding/SpotFindingSettings.h @@ -22,6 +22,26 @@ constexpr uint32_t StrongPixelLimit(size_t pixel_count) { return limit > UINT16_MAX ? limit : UINT16_MAX; } +// A connected component this size or smaller is judged on size alone. That is where the size bound +// stood before it was raised, so every component that used to be accepted still is. +constexpr int64_t SPOT_SHAPE_FREE_PIXELS = 50; +// A larger one must also be COMPACT: it has to fill at least this percent of the square its bounding +// box fits inside. A Bragg reflection is round and fills about half of that square however bright it +// is; an ice arc, a cosmic-ray track or a lit detector row fills a fifth or less, and those are what +// an upper bound on spot size was ever for. Measured on the strong rotation set, this keeps 97% of the +// genuine wide reflections a raised bound admits and none of the ice arcs. +// A shape test rather than a size test is what XDS and DISTL bound spots with; neither uses this +// statistic, but the choice of shape over size is theirs. +// Following Kabsch (2010) Acta Cryst. D66, 125-132 and Zhang, Sauter et al. (2006) J. Appl. Cryst. 39, 112-119 +constexpr int64_t SPOT_MIN_FILL_PERCENT = 20; + +// Integer arithmetic throughout, so the host and the GPU extractor agree by construction. +constexpr bool SpotShapeAccepted(int64_t pixel_count, int64_t bbox_side) { + if (pixel_count <= SPOT_SHAPE_FREE_PIXELS) + return true; + return pixel_count * 100 >= SPOT_MIN_FILL_PERCENT * bbox_side * bbox_side; +} + struct SpotFindingSettings { bool enable = true; float signal_to_noise_threshold = 4.0; // STRONG_PIXEL in XDS @@ -31,7 +51,9 @@ struct SpotFindingSettings { // is kept (see MXAnalysisWithoutFPGA::Analyze); a value fixes it. Defaults to a concrete value, so // the online receiver and the FPGA path keep the single-pass fixed behaviour unless set otherwise. std::optional min_pix_per_spot = 2; - int64_t max_pix_per_spot = 50; // Maximum pixels per spot + // Maximum pixels per spot. A component above SPOT_SHAPE_FREE_PIXELS must also pass + // SpotShapeAccepted, so this bounds how large a ROUND spot may be, not how bright. + int64_t max_pix_per_spot = 200; // High-resolution limit for spot finding [A]. std::nullopt = as far as the detector reaches, i.e. no // resolution clipping of the detection at all (DiffractionExperiment::GetDetectorMaxResolution_A // supplies the number where one is needed, e.g. for the spot plot's shells). diff --git a/image_analysis/spot_finding/StrongPixelSet.cpp b/image_analysis/spot_finding/StrongPixelSet.cpp index 6c631236a..0c12cc241 100644 --- a/image_analysis/spot_finding/StrongPixelSet.cpp +++ b/image_analysis/spot_finding/StrongPixelSet.cpp @@ -8,6 +8,7 @@ // The union-find and the two-scan structure are theirs. How a pixel's earlier neighbours are FOUND // is not: see sparseccl below. +#include #include #include "StrongPixelSet.h" @@ -40,7 +41,7 @@ uint32_t StrongPixelSet::make_union(uint32_t e1, uint32_t e2) { return e; } -std::vector StrongPixelSet::sparseccl() { +std::vector StrongPixelSet::sparseccl(const SpotFindingSettings &settings) { L.resize(pixels.size()); unsigned int labels = 0; @@ -89,20 +90,39 @@ std::vector StrongPixelSet::sparseccl() { } std::vector spots(labels); + // The bounding box travels with the accumulation rather than on DiffractionSpot: it is wanted + // only to judge the shape here, and carrying it further would leave a member that + // ConvertToImageCoordinates silently invalidates. + std::vector min_col(labels, UINT16_MAX), max_col(labels, 0); + std::vector min_line(labels, UINT16_MAX), max_line(labels, 0); - for (uint32_t i = 0; i < L.size(); i++) - spots[L[i]].AddPixel(pixels[i].col, pixels[i].line, pixels[i].counts); + for (uint32_t i = 0; i < L.size(); i++) { + const uint32_t l = L[i]; + spots[l].AddPixel(pixels[i].col, pixels[i].line, pixels[i].counts); + min_col[l] = std::min(min_col[l], pixels[i].col); + max_col[l] = std::max(max_col[l], pixels[i].col); + min_line[l] = std::min(min_line[l], pixels[i].line); + max_line[l] = std::max(max_line[l], pixels[i].line); + } - return spots; + std::vector out; + for (uint32_t l = 0; l < labels; l++) { + if (spots[l].PixelCount() > settings.max_pix_per_spot) + continue; + const int64_t w = static_cast(max_col[l]) - min_col[l] + 1; + const int64_t h = static_cast(max_line[l]) - min_line[l] + 1; + if (SpotShapeAccepted(spots[l].PixelCount(), std::max(w, h))) + out.push_back(spots[l]); + } + return out; } void StrongPixelSet::FindComponentsImage(const SpotFindingSettings &settings, std::vector &spots) { // No StrongPixelLimit test here: the caller knows how big the image is and has already applied it. - for (const auto &spot: sparseccl()) { - if (spot.PixelCount() <= settings.max_pix_per_spot) - spots.push_back(spot); - } + // Size and shape are applied by sparseccl. + for (const auto &spot: sparseccl(settings)) + spots.push_back(spot); } void StrongPixelSet::FindSpots(const DiffractionExperiment &experiment, const SpotFindingSettings &settings, @@ -110,9 +130,8 @@ void StrongPixelSet::FindSpots(const DiffractionExperiment &experiment, const Sp // Per module, so the bar is the module's own - and ReadFPGAOutput has already refused anything // past max_strong_pixel_per_module, far below it. if (!pixels.empty() && (strong_pixel_count < StrongPixelLimit(RAW_MODULE_SIZE))) { - for (const auto &spot: sparseccl()) { - if ((spot.PixelCount() <= settings.max_pix_per_spot) - && (spot.PixelCount() >= settings.min_pix_per_spot.value_or(2))) { + for (const auto &spot: sparseccl(settings)) { + if (spot.PixelCount() >= settings.min_pix_per_spot.value_or(2)) { auto s = spot; s.ConvertToImageCoordinates(experiment, module_number); spots.push_back(s); diff --git a/image_analysis/spot_finding/StrongPixelSet.h b/image_analysis/spot_finding/StrongPixelSet.h index 5026122d3..7004dd510 100644 --- a/image_analysis/spot_finding/StrongPixelSet.h +++ b/image_analysis/spot_finding/StrongPixelSet.h @@ -26,12 +26,19 @@ class StrongPixelSet { uint32_t find_root(uint32_t e); uint32_t make_union(uint32_t e1, uint32_t e2); - std::vector sparseccl(); + // Labels the strong pixels and returns the components that pass BOTH the size bound and + // SpotShapeAccepted. min-pix is deliberately not applied here (see FindComponentsImage). + std::vector sparseccl(const SpotFindingSettings &settings); public: void ReadFPGAOutput(const DiffractionExperiment& experiment, const DeviceOutput& output); StrongPixelSet(); + // Pixels MUST arrive in raster order - line ascending, and column ascending within a line. The + // component search walks the previous line with a forward-only cursor and would miss neighbours + // handed to it out of order. Both callers satisfy this by construction, each scanning a bitmap in + // ascending flat index: ImageSpotFinder::ExtractComponentsHost over the image, ReadFPGAOutput over + // a module (whose row length is a multiple of 32, so no word straddles two lines). void AddStrongPixel(uint16_t col, uint16_t line, int32_t photons = 1); void FindSpots(const DiffractionExperiment &experiment, const SpotFindingSettings &settings, std::vector &spots, uint16_t module_number); diff --git a/tests/StrongPixelSetTest.cpp b/tests/StrongPixelSetTest.cpp index c2243e100..97649d25d 100644 --- a/tests/StrongPixelSetTest.cpp +++ b/tests/StrongPixelSetTest.cpp @@ -29,17 +29,12 @@ TEST_CASE("StrongPixelSet_BuildSpots","[StrongPixelSet]") { settings.min_pix_per_spot = 3; settings.max_pix_per_spot = 200; + // In raster order, which is what AddStrongPixel requires and what both callers produce. std::vector spots; StrongPixelSet strong_pixel_set; - strong_pixel_set.AddStrongPixel(7,105); - strong_pixel_set.AddStrongPixel(7,106); - strong_pixel_set.AddStrongPixel(7,104); - strong_pixel_set.AddStrongPixel(6,106); - strong_pixel_set.AddStrongPixel(6,105); - strong_pixel_set.AddStrongPixel(6,104); - strong_pixel_set.AddStrongPixel(8,106); - strong_pixel_set.AddStrongPixel(8,105); - strong_pixel_set.AddStrongPixel(8,104); + for (uint16_t line = 104; line <= 106; line++) + for (uint16_t col = 6; col <= 8; col++) + strong_pixel_set.AddStrongPixel(col, line); strong_pixel_set.FindSpots(experiment, settings, spots, 0); REQUIRE(spots.size() == 1); @@ -49,6 +44,40 @@ TEST_CASE("StrongPixelSet_BuildSpots","[StrongPixelSet]") { REQUIRE(spots[0].RawCoord().y == Catch::Approx(105.0)); } +TEST_CASE("StrongPixelSet_LargeSpotMustBeCompact","[StrongPixelSet]") { + SpotFindingSettings settings; // max_pix_per_spot 200, so all three pass on size alone + + // A filled 9x9 square: 81 pixels in a 9-wide box, so it fills all of it. A Bragg reflection is + // this shape however bright it is, and this is what a raised size bound is meant to admit. + StrongPixelSet square; + for (uint16_t line = 100; line < 109; line++) + for (uint16_t col = 200; col < 209; col++) + square.AddStrongPixel(col, line); + std::vector spots; + square.FindComponentsImage(settings, spots); + REQUIRE(spots.size() == 1); + CHECK(spots[0].PixelCount() == 81); + + // The same 81 pixels drawn as a line: an ice arc, a cosmic-ray track or a lit detector row. It + // fills 1/81 of its bounding square, so it is refused however many pixels it has. + StrongPixelSet streak; + for (uint16_t col = 200; col < 281; col++) + streak.AddStrongPixel(col, 100); + spots.clear(); + streak.FindComponentsImage(settings, spots); + CHECK(spots.empty()); + + // Shape is asked of large components only. A short line of the same one-pixel width is below + // SPOT_SHAPE_FREE_PIXELS and is kept, exactly as it was before the bound was raised. + StrongPixelSet short_streak; + for (uint16_t col = 200; col < 210; col++) + short_streak.AddStrongPixel(col, 100); + spots.clear(); + short_streak.FindComponentsImage(settings, spots); + REQUIRE(spots.size() == 1); + CHECK(spots[0].PixelCount() == 10); +} + /* TEST_CASE("StrongPixelSet_ReadFPGAOutput_1","[StrongPixelSet]") { DiffractionExperiment experiment(DetectorGeometry(8, 2, 8, 36, true)); -- 2.54.0 From 844b461232936b5e029be74aee4d1746f9fcb876 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 15:04:05 +0200 Subject: [PATCH 009/179] Spot finding: say honestly that the compactness gate is inert The commit that added it claimed the shape test was what made raising the size bound safe - "with a flat 200 alone, two battery crystals moved; with the shape test both are unchanged to three decimals". That is wrong. Those two runs came from different builds, and the comparison never isolated the gate at all. Isolated properly, gate on against gate off in ONE build: byte-identical .hkl on three rotation crystals, the strongly diffracting one included. Not merely close - identical, every digit of every merged reflection. On that crystal the gate rejects 5 to 22 of the 110 to 440 oversize components a frame, so essentially all of what that commit bought, it bought by raising the bound from 50 to 200. Keep it, but for the reason that survives rather than the one that did not. The raised bound admits components up to 200 pixels, and on data carrying ice arcs, cosmic-ray tracks or a lit detector row those are what arrives; none of the crystals measured here carry that population in quantity, which is exactly why the gate reads inert on all of them. So it bounds the SHAPE of what the larger size bound now lets through, and no experiment on these crystals can validate it - an inert gate scores identically whatever its threshold. The comment also carried a "keeps 97% of the genuine wide reflections and none of the ice arcs" that came from a simulation over component lists, not from the pipeline. Removed rather than restated: the pipeline measurement says no reflection anywhere changed. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H --- docs/CPU_DATA_ANALYSIS.md | 2 +- image_analysis/spot_finding/SpotFindingSettings.h | 10 ++++++++-- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 256cb58fe..5165d9545 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -372,7 +372,7 @@ Strong pixels are grouped into connected components (adjacent strong pixels) usi Spot-level filters include minimum/maximum pixel count and resolution limits. -**The upper bound is on how large a *round* spot may be, not on how bright one may be.** Under the self-calibrating threshold (§3.2) a component's area above the contour grows as $\sigma^2\ln(A/T)$ — without bound in the peak amplitude $A$ — so an upper bound on pixel count alone becomes an *intensity ceiling*: measured on a strongly diffracting crystal, footprints run 3 px at 30–100 counts to 50 px above 10 000, four times the slope the fixed local-box test gives, and the old bound of 50 discarded the brightest reflections on every image. The bound is therefore 200 (the same as CrystFEL `peakfinder8`'s `--max-pix-count`; XDS has no such parameter at all and guards on shape instead, with `SPOT_MAXIMUM-CENTROID`), and a component above 50 pixels must in addition **fill at least a fifth of the square its bounding box fits inside**. A Bragg reflection is round and fills about half of that square however bright it is; an ice arc, a cosmic-ray track or a lit detector row fills a fifth or less, which is what an upper bound was ever protecting against. Below 50 pixels no shape is asked for, so everything accepted before still is. The test is integer arithmetic, so the host and the GPU extractor agree by construction. +**The upper bound is on how large a *round* spot may be, not on how bright one may be.** Under the self-calibrating threshold (§3.2) a component's area above the contour grows as $\sigma^2\ln(A/T)$ — without bound in the peak amplitude $A$ — so an upper bound on pixel count alone becomes an *intensity ceiling*: measured on a strongly diffracting crystal, footprints run 3 px at 30–100 counts to 50 px above 10 000, four times the slope the fixed local-box test gives, and the old bound of 50 discarded the brightest reflections on every image. The bound is therefore 200 (the same as CrystFEL `peakfinder8`'s `--max-pix-count`; XDS has no such parameter at all and guards on shape instead, with `SPOT_MAXIMUM-CENTROID`), and a component above 50 pixels must in addition **fill at least a fifth of the square its bounding box fits inside**. A Bragg reflection is round and fills about half of that square however bright it is; an ice arc, a cosmic-ray track or a lit detector row fills a fifth or less, which is what an upper bound was ever protecting against. Below 50 pixels no shape is asked for, so everything accepted before still is. The shape test is inert on every dataset it has been measured on — it exists to bound the *shape* of what the larger size bound now admits, on data carrying arcs or tracks, not because the crystals measured here needed it. The test is integer arithmetic, so the host and the GPU extractor agree by construction. The host implementation (`StrongPixelSet::sparseccl`) is the SparseCCL of the ACTS/traccc project: it runs over the strong pixels sorted row-major, uses a sliding window over the previous line and a diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index d610277bc..8b79aab6c 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/image_analysis/spot_finding/SpotFindingSettings.h @@ -28,8 +28,14 @@ constexpr int64_t SPOT_SHAPE_FREE_PIXELS = 50; // A larger one must also be COMPACT: it has to fill at least this percent of the square its bounding // box fits inside. A Bragg reflection is round and fills about half of that square however bright it // is; an ice arc, a cosmic-ray track or a lit detector row fills a fifth or less, and those are what -// an upper bound on spot size was ever for. Measured on the strong rotation set, this keeps 97% of the -// genuine wide reflections a raised bound admits and none of the ice arcs. +// an upper bound on spot size was ever for. +// +// It is INERT on every dataset it has been measured on: gate on and gate off give a byte-identical +// merge on three rotation crystals including the strongly diffracting one the raised bound was written +// for. What it guards against is a population none of them carry in quantity - the raised bound admits +// components up to 200 px, and on ice-arc or cosmic-track-heavy data those are what arrives. So it is +// kept as a bound on the SHAPE of what the larger size bound now lets through, not because any +// measurement here needed it. // A shape test rather than a size test is what XDS and DISTL bound spots with; neither uses this // statistic, but the choice of shape over size is theirs. // Following Kabsch (2010) Acta Cryst. D66, 125-132 and Zhang, Sauter et al. (2006) J. Appl. Cryst. 39, 112-119 -- 2.54.0 From 9f2696c1acc4f7c859c2f018fe0322083cb269c4 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 15:07:56 +0200 Subject: [PATCH 010/179] Adaptive spot finding: the ring threshold replaces the floor, not the local test The self-calibrating finder was meant to replace the classic finder's FIXED PHOTON FLOOR with a per-resolution-ring threshold read off the image's own noise. As written it replaced the local-box SNR test as well, and that is the defect: a whole-ring threshold is an ABSOLUTE contour with no feedback from a pixel's own surroundings, so the area a spot puts above it grows as sigma^2*ln(peak/threshold) and never saturates. Measured on a strongly diffracting rotation set, the detected footprint grows by +8.05 pixels per e-fold of peak, so the brightest reflections came out as 100-500 pixel blobs and were then discarded for exceeding the size bound - every one of the ten strongest on an image. Intersecting with the local box gives -0.24 pixels per e-fold, the classic finder's own number to two decimals. WHY the local box is the right partner, rather than merely the incumbent: it is a prominence rule whose reference level is a 961-pixel mean. A spot inflates the box's own variance and the peak divides out of the acceptance test, so it cuts at a fixed FRACTION of the spot's own height. Referring that level to fewer pixels makes it inherit their shot noise - at FIXED footprint, estimating the level from 961 pixels, from 25, and from the single maximum gives centroid residuals of 0.524, 0.539 and 0.656 - so flat growth and a stable centroid turn out to be two ends of one dial. A contour on the bare maximum has the flattest growth of anything tried (+0.1) and merges worst. The two arms bind in different regimes, which is why intersecting beats choosing: on serial stills the ring threshold is 0.6x the classic floor, on this rotation sweep 2.3-6.0x. Stills are a strict no-op - 175 components against 175, identical per frame - so the +40% in stills indexing that the adaptive threshold was introduced for is untouched. What it buys, stated as one fact rather than two. Across five geometry pins spanning 1.1 mm it indexes the most frames of any arm tried, 0.831 against 0.803, and integrates 3.04 to 5.76% more observations - but those are the SAME number: regressing observation count on indexing rate over four arms leaves residuals of +/-0.7 percentage points against swings of -7 to +4.5%, so the extra observations ARE the extra indexed frames, not better data per frame. CC1/2, the only statistic here carrying per-observation quality, is +0.66 at one pin and -0.06 at the other: not harmed, not improved. , ISa and R_meas cannot arbitrate on this data - across those pins each crosses zero as a monotone function of the pin. WHY an absolute contour indexes fewer frames, when its spot list is equal or better on every axis measured - recall, top-1000 recall, centroid, ice fraction, component count - is the interesting part, and it is not a detection effect at all: ITS OWN SIZE BOUND DELETES THE BRIGHTEST REFLECTIONS ON THE FRAME. A component is discarded because it grew past 200 px, and it grew past 200 px because it was bright, so the deletions are drawn from the head of the indexing budget rather than uniformly from it: they are 11x enriched in the top 250 of the thousand spots handed to the indexer, and the bound's own real deletions sit at MEDIAN RANK 12. Turning the bound off recovers 66% and 50% of the deficit at the two pins, against a bar registered at 33% before the run. Three of us dismissed this for most of a day on the grounds that the gates delete only ~4% of what is detected. That arithmetic was right and the denominator was wrong - a rate is not an impact when the thing being lost is selected for the property that makes it matter. Reworking the bound instead was measured and rejected: it recovers half the deficit, and it cannot be done without re-admitting what the bound is for - 68 components past 200 px, of which 8 are real and 60 are junk, where the intersect gets the 8 without the 60. The residual once the bound is off, +1.08%/+1.70%, is the contour itself. Component merging is ruled out separately: geometrically impossible here, 33.9 px minimum reflection separation against components spanning 10 px. So is a ranking effect - the intersect's lead runs +0.06% at --max-spots 250, +3.46% at 1000 and +14.26% at 2000, which is backwards for a selection artefact. Costs 0.48 ms per image in the finder, and 0.044 px of bright-spot centroid precision - measured convention-free, by fitting a line to a reflection's own centroid across five frames, after an XDS-referenced figure proved to be four fifths aperture convention. It also makes the compactness gate above it safe. On the absolute contour that gate is net damage, deleting 37 genuine reflections per ten frames; once the footprint stops growing nothing reaches its threshold at all. Also fixes a real but unexercised defect in PoissonThreshold, where the exact tail handed over to a normal approximation with a step. It changes nothing here: the clipped ring sigma is over-dispersed 1.2-4.9x against sqrt(mu) because it still contains diffraction, so the Gaussian arm wins every ring above mu=50 and none of the 522 thresholds move. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H --- docs/CHANGELOG.md | 4 +- .../spot_finding/AdaptiveSpotFinderCPU.cpp | 34 +++++++--- .../spot_finding/AdaptiveSpotFinderCPU.h | 33 ++++++++-- .../spot_finding/AdaptiveSpotFinderGPU.cu | 51 ++++++++------- .../spot_finding/AdaptiveSpotFinderGPU.h | 29 ++++----- .../spot_finding/AdaptiveThreshold.h | 23 ++++++- .../spot_finding/ImageSpotFinderGPU.h | 6 ++ .../spot_finding/SpotFindingSettings.h | 17 +++-- tests/AdaptiveSpotFinderCPUTest.cpp | 62 +++++++++++++++++++ tests/AdaptiveThresholdTest.cpp | 17 ++++- 10 files changed, 215 insertions(+), 61 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index b2940fb9c..37030a2aa 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,8 +1,10 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* Self-calibrating spot detection intersects its per-resolution-ring threshold with the local signal-to-noise test again, instead of replacing it. The ring threshold takes the place of the fixed photon floor and nothing else; standing alone it followed a bright reflection's skirt outwards, so on a strongly diffracting crystal the brightest reflections were detected as 100-500 pixel blobs and then discarded for being too large. +* The self-calibrating threshold no longer steps where it switches from the exact Poisson tail to a normal approximation. * A spot may hold up to 200 connected pixels rather than 50. Under the self-calibrating threshold a reflection's footprint grows with its brightness, so the old bound acted as an intensity ceiling and discarded the strongest reflections of a strongly diffracting crystal - on one such set, every one of the ten brightest on an image. -* A spot larger than 50 pixels must also be compact, filling a fifth of the square its bounding box fits inside. An ice arc, a cosmic-ray track or a lit detector row is refused however many pixels it has, which is what the size bound used to do and now does without capping brightness. +* A spot larger than 50 pixels must also be compact, filling a fifth of the square its bounding box fits inside, so that an ice arc, a cosmic-ray track or a lit detector row is refused however many pixels it has. It bounds the shape of what the raised size limit admits; on data without that population it changes nothing. * Eight further hexagonal-ice bands between 1.472 and 1.170 A are recognised. The measured list they extend stops at 1.522 A because its source went no further, not because ice does, and on a detector reaching past 1.5 A the unlisted rings were left in the data - on a strongly diffracting set they were 44% of every image's spots. Ice handling still only applies where a run trips the ice gate. * The rugnux results report carries `JFJOCH_DATASET_SETTINGS=`, the geometry the run integrated at written as one line of JSON in the form `jfjoch_broker` takes it, so a refined beam centre and distance can be carried back to the instrument. * `/entry/MX/strongPixels` and the strong-pixel plot are filled on the CPU/GPU analysis path, not only behind the FPGA, so an image the spot finder gave up on can be told from one that did not diffract. diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp index a60e78327..9f7cd16ad 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp @@ -9,8 +9,8 @@ #include "AdaptiveThreshold.h" AdaptiveSpotFinderCPU::AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &in_mapping) - : ImageSpotFinder(static_cast(in_mapping.GetWidth()), - static_cast(in_mapping.GetHeight())), + : ImageSpotFinderCPU(static_cast(in_mapping.GetWidth()), + static_cast(in_mapping.GetHeight())), mapping(in_mapping) { const size_t nbins = mapping.GetBinNumber(); ring_sum.assign(nbins, 0); @@ -20,6 +20,7 @@ AdaptiveSpotFinderCPU::AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping & ring_sigma.assign(nbins, 0.0f); ring_thr.assign(nbins, 0.0f); ring_bkg.assign(nbins, NAN); + ring_bits.assign(OutputSize(), 0); } // Accumulate per-ring mean/variance from the raw (photon) image. clip_k <= 0 -> use every valid @@ -61,9 +62,7 @@ void AdaptiveSpotFinderCPU::AccumulateRings(const ImagePreprocessorBuffer &image void AdaptiveSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) { - const auto &pixel_to_bin = mapping.GetPixelToBin(); const size_t nbins = ring_sum.size(); - const size_t npix = static_cast(width) * height; // --- Stage A: robust per-ring background (one plain pass + two sigma-clip passes) --- AccumulateRings(image, 0.0f); @@ -104,9 +103,28 @@ void AdaptiveSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image, ? g_thr : adaptive_threshold::RingThreshold(ring_mean[b], ring_sigma[b], p, z); - // --- Stage C: flag strong pixels into the bit buffer (value >= ring threshold) --- + // --- Stage C: the ring threshold, intersected with the classic local-box SNR test --- + FlagRings(image); + + if (settings.signal_to_noise_threshold <= 0.0f) { + // No local test asked for: the ring threshold alone decides, as the fixed photon floor + // alone would in the classic finder. + output_buffer = ring_bits; + return; + } + + // The ring threshold IS the photon floor here, so the local pass must not apply another one. + SpotFindingSettings local = settings; + local.photon_count_threshold = 0; + ImageSpotFinderCPU::Detect(image, local); for (size_t i = 0; i < OutputSize(); ++i) - output_buffer[i] = 0; + output_buffer[i] &= ring_bits[i]; +} + +void AdaptiveSpotFinderCPU::FlagRings(const ImagePreprocessorBuffer &image) { + const auto &pixel_to_bin = mapping.GetPixelToBin(); + const size_t nbins = ring_thr.size(); + const size_t npix = static_cast(width) * height; std::bitset<32> out = 0; for (size_t pxl = 0; pxl < npix; ++pxl) { @@ -122,10 +140,10 @@ void AdaptiveSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image, if (strong) out.set(bit); if (bit == 31) { - output_buffer[pxl / 32] = out.to_ulong(); + ring_bits[pxl / 32] = out.to_ulong(); out.reset(); } } if (npix % 32 != 0) - output_buffer[OutputSize() - 1] = out.to_ulong(); + ring_bits[OutputSize() - 1] = out.to_ulong(); } diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h index c78d732b0..3604d1b9c 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h @@ -5,7 +5,7 @@ #include -#include "ImageSpotFinder.h" +#include "ImageSpotFinderCPU.h" #include "SpotFindingSettings.h" #include "../../common/AzimuthalIntegrationMapping.h" @@ -22,9 +22,29 @@ // whose (peak-excluded) background mean is mu, the threshold is the smallest count whose Poisson // upper tail is <= p = E / N_pixels, max'd with a Gaussian arm mu + z*sigma to absorb read/flat-field // excess. Because it is set from the image's own noise, the SAME E lands at ~12 photons on the first -// set and ~5 on the weaker one with no user input. Detection then is simply value > ring_threshold, -// fed to the same connected-component builder as the classic finder. -class AdaptiveSpotFinderCPU : public ImageSpotFinder { +// set and ~5 on the weaker one with no user input. +// +// The ring threshold replaces the floor and ONLY the floor: the classic finder's local-box SNR test +// still has to pass, which is why this engine runs it (ImageSpotFinderCPU) and intersects the two +// masks. +// +// A whole-ring threshold is an ABSOLUTE contour with no feedback from the pixel's own surroundings, +// so the area a spot puts above it grows as sigma^2 * ln(peak/threshold) and never saturates: on a +// strongly diffracting rotation set the detected footprint grows by 8 pixels per e-fold of peak, so +// the brightest reflections came out as 100-500 pixel blobs. The local box has no such contour. The +// spot inflates the box's own variance, and the peak divides out of the acceptance test, so the box +// cuts every spot at roughly a fixed FRACTION of its own height - a peak-relative contour. Measured +// on the same set, the footprint then grows by -0.2 pixels per e-fold, i.e. not at all, and lands on +// the classic finder's own number to two decimals. +// +// The two arms bind in different regimes, which is the point of intersecting rather than choosing. +// On serial stills the ring background is a fraction of a count and the ring threshold lands BELOW +// the fixed floor the classic finder would use, so the ring arm decides and the local box passes +// everything - which is the whole reason this engine exists. On a bright rotation set the ring +// background is tens of counts, the ring threshold lands several times ABOVE that floor, and the +// local box decides. The engine is therefore never worse than the classic finder on footprint, and +// never worse than a fixed floor on a weak background. +class AdaptiveSpotFinderCPU : public ImageSpotFinderCPU { const AzimuthalIntegrationMapping &mapping; // per-ring scratch, sized to the mapping's bin count @@ -40,8 +60,13 @@ class AdaptiveSpotFinderCPU : public ImageSpotFinder { // ring_mean of the last Detect(), NaN where the ring holds too few pixels to be its own background. // Kept separately because ring_mean carries the previous frame's value for an empty ring. std::vector ring_bkg; + // Pixels at or above their ring's threshold, packed like output_buffer. Intersected with the + // local-box mask that ImageSpotFinderCPU::Detect leaves in output_buffer. + std::vector ring_bits; void AccumulateRings(const ImagePreprocessorBuffer &image, float clip_k); + // Fill ring_bits from the thresholds of the current frame. + void FlagRings(const ImagePreprocessorBuffer &image); public: explicit AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &mapping); diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu index 65ed5123b..2cafc5c06 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu @@ -276,14 +276,20 @@ __global__ void flag_strong(const int32_t *__restrict__ image, } } +// out &= mask, over the packed bit buffer. One word per thread; the buffer is one bit per pixel, so +// this is 1/32 of an image pass. +__global__ void and_bits(uint32_t *__restrict__ out, const uint32_t *__restrict__ mask, size_t nwords) { + for (size_t i = blockIdx.x * blockDim.x + threadIdx.x; i < nwords; i += blockDim.x * gridDim.x) + out[i] &= mask[i]; +} + } // namespace AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &in_mapping, std::shared_ptr in_stream) - : ImageSpotFinder(static_cast(in_mapping.GetWidth()), - static_cast(in_mapping.GetHeight()), false), + : ImageSpotFinderGPU(static_cast(in_mapping.GetWidth()), + static_cast(in_mapping.GetHeight()), std::move(in_stream)), mapping(in_mapping), - stream(in_stream), nbins(in_mapping.GetBinNumber()), npix(in_mapping.GetPixelToBin().size()), gpu_sum(nbins), @@ -294,7 +300,7 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping & gpu_sum_corr(nbins), gpu_sum2_corr(nbins), gpu_thr(nbins), - gpu_strong(OutputSize()), + gpu_ring(OutputSize()), host_sum(nbins), host_sum2(nbins), host_count(nbins), @@ -308,8 +314,6 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping & prof_sum_reg(prof_sum), prof_sum2_reg(prof_sum2), prof_count_reg(prof_count), - extractor(static_cast(in_mapping.GetWidth()), - static_cast(in_mapping.GetHeight()), std::move(in_stream)), last_profile(in_mapping) { // The current device, not device 0: callers round-robin engines across GPUs, so device 0's shared @@ -464,27 +468,30 @@ void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image, if (host_thr.empty()) { // Nothing valid to threshold against: leave no strong pixels for the extractor to build on. - cuda_err(cudaMemsetAsync(gpu_strong, 0, OutputByteSize(), *stream)); + cuda_err(cudaMemsetAsync(gpu_out_1, 0, OutputByteSize(), *stream)); cuda_err(cudaStreamSynchronize(*stream)); return; } - // --- Stage C: flag strong pixels into the bit buffer (value >= ring threshold) --- + // --- Stage C: the ring threshold, intersected with the classic local-box SNR test --- cuda_err(cudaMemcpyAsync(gpu_thr, host_thr.data(), sizeof(float) * nbins, cudaMemcpyHostToDevice, *stream)); - cuda_err(cudaMemsetAsync(gpu_strong, 0, OutputByteSize(), *stream)); + cuda_err(cudaMemsetAsync(gpu_ring, 0, OutputByteSize(), *stream)); flag_strong<<>>( - image.getGPUBuffer(), gpu_pixel_to_bin->get(), gpu_thr, gpu_strong, npix, nbins); + image.getGPUBuffer(), gpu_pixel_to_bin->get(), gpu_thr, gpu_ring, npix, nbins); + + if (settings.signal_to_noise_threshold <= 0.0f) { + // No local test asked for: the ring threshold alone decides, as the fixed photon floor + // alone would in the classic finder. + cuda_err(cudaMemcpyAsync(gpu_out_1, gpu_ring, OutputByteSize(), cudaMemcpyDeviceToDevice, *stream)); + cuda_err(cudaStreamSynchronize(*stream)); + return; + } + + // The ring threshold IS the photon floor here, so the local pass must not apply another one. + SpotFindingSettings local = settings; + local.photon_count_threshold = 0; + ImageSpotFinderGPU::Detect(image, local); + and_bits<<>>(gpu_out_1, gpu_ring, OutputSize()); // The bit buffer stays on the device and ExtractComponents reads it there, on this same stream, - // so the ordering already guarantees flag_strong has finished. Waiting here only idled the host. -} - -void AdaptiveSpotFinderGPU::SetResolutionMaskBits(const std::vector &packed_mask) { - ImageSpotFinder::SetResolutionMaskBits(packed_mask); - extractor.SetResolutionMask(res_mask_bits); -} - -const std::vector &AdaptiveSpotFinderGPU::ExtractComponents(const ImagePreprocessorBuffer &image, - const SpotFindingSettings &settings) { - extractor.Extract(gpu_strong, image.getGPUBuffer(), settings, components); - return components; + // so the ordering already guarantees and_bits has finished. Waiting here only idled the host. } diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h index c15d48c0e..0a3bcb545 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h @@ -12,20 +12,23 @@ // - the azimuthal-integration profile (mean intensity per ring, in flat-field-corrected space), and // - the per-ring background (mean, sigma, peak-excluded via two sigma-clip passes) that sets the // self-calibrating spot-detection threshold (in raw photon counts). -// It then flags strong pixels (value >= ring threshold) into a packed bit buffer and hands that -// buffer - still on the device - to SpotExtractorGPU, which builds the spots there. +// The ring threshold is the photon-count FLOOR and nothing more: this engine derives from the classic +// ImageSpotFinderGPU and intersects the ring mask with that engine's local-box SNR mask, exactly as +// AdaptiveSpotFinderCPU does on the host (the reason the local test is not optional is written out +// there). The result is left in the inherited bit buffer, still on the device, and the inherited +// SpotExtractorGPU builds the spots from it. // // Numerically it reproduces AdaptiveSpotFinderCPU: the same three-pass robust background, the same // per-ring threshold formula (shared via AdaptiveThreshold.h, computed on the host once per frame), -// and the same raw-count detection test. The only differences from the CPU are those inherent to a -// GPU reduction (float per-ring accumulation in atomic order vs the CPU's serial double sums), which -// shift a handful of borderline pixels at most. The corrected sums for the azint profile are +// and the same intersection with the local test. The only differences from the CPU are those inherent +// to a GPU reduction (float per-ring accumulation in atomic order vs the CPU's serial double sums), +// which shift a handful of borderline pixels at most. The corrected sums for the azint profile are // accumulated in the SAME plain first pass, so one reduction feeds both products. #include #include -#include "ImageSpotFinder.h" +#include "ImageSpotFinderGPU.h" #include "SpotExtractorGPU.h" #include "SpotFindingSettings.h" #include "../../common/AzimuthalIntegrationProfile.h" @@ -33,9 +36,8 @@ #include "../indexing/CUDAMemHelpers.h" #include "../indexing/CudaSharedTables.h" -class AdaptiveSpotFinderGPU : public ImageSpotFinder { +class AdaptiveSpotFinderGPU : public ImageSpotFinderGPU { const AzimuthalIntegrationMapping &mapping; - std::shared_ptr stream; const int nbins; const size_t npix; @@ -68,9 +70,10 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder { CudaDevicePtr gpu_sum_corr; CudaDevicePtr gpu_sum2_corr; - // Per-ring detection threshold (host-computed, uploaded) and the strong-pixel bit buffer. + // Per-ring detection threshold (host-computed, uploaded) and the ring-threshold mask that is + // intersected into the inherited bit buffer. CudaDevicePtr gpu_thr; - CudaDevicePtr gpu_strong; + CudaDevicePtr gpu_ring; // Host mirrors of the small per-ring transfers. std::vector host_sum; // clipped raw sum } input to the host threshold computation @@ -94,8 +97,6 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder { CudaRegisteredVector prof_sum2_reg; CudaRegisteredVector prof_count_reg; - SpotExtractorGPU extractor; // builds the spots from gpu_strong without it leaving the device - AzimuthalIntegrationProfile last_profile; // filled every Run(), retrievable via GetProfile() // One reduction pass over the image into the raw accumulators. clip_k <= 0 -> plain pass (all @@ -114,10 +115,6 @@ public: AdaptiveSpotFinderGPU &operator=(const AdaptiveSpotFinderGPU &) = delete; void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; - void SetResolutionMaskBits(const std::vector &packed_mask) override; - [[nodiscard]] uint32_t StrongPixelCount() const override { return extractor.StrongPixelCount(); } - const std::vector &ExtractComponents(const ImagePreprocessorBuffer &image, - const SpotFindingSettings &settings) override; // The azimuthal profile computed as a byproduct of the last Detect() - lets this engine replace the // separate azint pass in the analysis pipeline. diff --git a/image_analysis/spot_finding/AdaptiveThreshold.h b/image_analysis/spot_finding/AdaptiveThreshold.h index 047ee4698..745c64d27 100644 --- a/image_analysis/spot_finding/AdaptiveThreshold.h +++ b/image_analysis/spot_finding/AdaptiveThreshold.h @@ -61,10 +61,27 @@ inline double NormalQuantile(double p) { // significance floor while the background is countable (it carries the sqrt(mu) shot-noise // implicitly, so a bright low-resolution ring gets a high threshold). It DEGENERATES at mu -> 0 // (a single photon on a zero background is "significant"), which is why it is max'd with a -// read-noise-floored Gaussian arm by the caller. Short-circuits to Gaussian for large mu. +// read-noise-floored Gaussian arm by the caller. +// +// Past the summation limit the quantile is taken from the Cornish-Fisher expansion +// (Cornish and Fisher (1938) Rev. Int. Stat. Inst. 5, 307-320), whose skewness +// term (z^2-1)/6 is what a plain mu + z*sqrt(mu) leaves out. At the 4-6 sigma this operating point +// works at, that term is 3-6 counts, so the Gaussian form alone stood BELOW the true Poisson +// quantile - and it did so with a step at the switch, since below it the exact quantile was used. +// Cornish-Fisher is within one count of the exact value at every mu, so the two arms now join +// smoothly. +// +// How much this is worth depends on which arm of RingThreshold wins, and on measured data it is +// often neither: where the ring background is over-dispersed (a clipped ring sigma that still +// carries the ring's own azimuthal structure, 1.2-4.9x sqrt(mu) on a strongly diffracting rotation +// set) the Gaussian arm is the larger of the two on every ring above mu = 50 and this correction +// changes no threshold at all. It is the right value to return regardless: a caller that ever sees +// the Poisson arm win up there would otherwise get a bar that jumps at mu = 50. inline float PoissonThreshold(double mu, double p, double z) { - if (mu > 50.0) - return static_cast(mu + z * std::sqrt(mu)); + // The summation below needs k up to about mu + z*sqrt(mu), and exp(-mu) has to stay normal. + constexpr double SUM_LIMIT = 200.0; + if (mu > SUM_LIMIT) + return static_cast(mu + z * std::sqrt(mu) + (z * z - 1.0) / 6.0); if (mu < 1e-6) mu = 1e-6; const double target = 1.0 - p; double pmf = std::exp(-mu); diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.h b/image_analysis/spot_finding/ImageSpotFinderGPU.h index d63f4acd1..1d9382771 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.h +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.h @@ -11,12 +11,18 @@ #include "../indexing/CUDAMemHelpers.h" class ImageSpotFinderGPU : public ImageSpotFinder { +protected: + // Protected rather than private because AdaptiveSpotFinderGPU derives from this engine: it is + // this same local-box detection with the fixed photon floor replaced by a per-resolution-ring + // one, so it reuses the stream, the bit buffers and the extractor rather than owning a second + // set of them. std::shared_ptr stream; CudaDevicePtr gpu_out_0; CudaDevicePtr gpu_out_1; // holds the strong-pixel bits after Detect() SpotExtractorGPU extractor; +private: const int numberOfCudaThreads = 128; // #threads per block that should work well for Nvidia L4 const int numberOfWaves = 32; // #waves that should work well for Nvidia L4 const int windowSizeLimit = 32; // limit on the window size (2nby+1, 2nbx+1) to prevent shared memory problems diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index 8b79aab6c..897408920 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/image_analysis/spot_finding/SpotFindingSettings.h @@ -32,10 +32,19 @@ constexpr int64_t SPOT_SHAPE_FREE_PIXELS = 50; // // It is INERT on every dataset it has been measured on: gate on and gate off give a byte-identical // merge on three rotation crystals including the strongly diffracting one the raised bound was written -// for. What it guards against is a population none of them carry in quantity - the raised bound admits -// components up to 200 px, and on ice-arc or cosmic-track-heavy data those are what arrives. So it is -// kept as a bound on the SHAPE of what the larger size bound now lets through, not because any -// measurement here needed it. +// for. The raised bound is what did nearly all of the work there - of the components an absolute ring +// contour pushes past 50 pixels, this shape test rejects a few percent and the raise re-admits the +// rest. +// +// Both only matter while a spot's footprint can grow with its brightness, and since the adaptive +// finder intersects its ring threshold with the local-box SNR test the footprint is peak-relative and +// no component on that set reaches 50 pixels at all. +// +// The shape test is kept anyway because it is the one of the two that does not go stale: a size bound +// is a bet on how large spots are, and the detectors and the detection rule both move underneath it - +// under every peak-relative detector tried, the 200-pixel bound became unreachable while this test +// still fired. The bound is kept alongside it because the cap the local box imposes scales with the +// spot WIDTH, so a set with wider spots than the ones measured here will reach past 50 pixels again. // A shape test rather than a size test is what XDS and DISTL bound spots with; neither uses this // statistic, but the choice of shape over size is theirs. // Following Kabsch (2010) Acta Cryst. D66, 125-132 and Zhang, Sauter et al. (2006) J. Appl. Cryst. 39, 112-119 diff --git a/tests/AdaptiveSpotFinderCPUTest.cpp b/tests/AdaptiveSpotFinderCPUTest.cpp index 819ccbf1b..0db0270a6 100644 --- a/tests/AdaptiveSpotFinderCPUTest.cpp +++ b/tests/AdaptiveSpotFinderCPUTest.cpp @@ -1,6 +1,7 @@ // SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include #include #include @@ -190,3 +191,64 @@ TEST_CASE("AdaptiveSpotFinderCPU_RingBackgroundExcludesPeaks", "[AdaptiveSpotFin for (size_t b = 0; b < clean.size(); b++) CHECK(std::isfinite(clean[b]) == (pixels_per_bin[b] >= 40)); } + +// The footprint of a bright reflection must not grow with its brightness. A whole-ring threshold is +// an absolute contour, so the area a Gaussian puts above it grows as sigma^2 ln(peak/threshold) - the +// same spot detected a hundred times brighter comes back tens of pixels larger, and an upper bound on +// spot size becomes an upper bound on spot INTENSITY. Intersecting with the local-box test removes +// that: the spot inflates the box's own variance, the peak divides out of the acceptance test, and the +// contour lands at a fixed fraction of the spot's own height whatever that height is. +TEST_CASE("AdaptiveSpotFinderCPU_FootprintDoesNotGrowWithBrightness", "[AdaptiveSpotFinder]") { + DiffractionExperiment x(DetJF4M()); + x.DetectorDistance_mm(80).BeamX_pxl(1030).BeamY_pxl(1080); + x.QSpacingForAzimInt_recipA(0.05).QRangeForAzimInt_recipA(0.05, 5.0); + x.GeometryTransformation(false); + + PixelMask pixel_mask(x); + AzimuthalIntegrationMapping mapping(x, pixel_mask); + const auto &pixel_to_bin = mapping.GetPixelToBin(); + + const size_t w = x.GetXPixelsNum(); + const size_t h = x.GetYPixelsNum(); + + size_t spot_row = 0, spot_col = 0; + for (size_t row = 400; row < h - 400 && spot_row == 0; row++) + for (size_t col = 400; col < w - 400; col++) + if (pixel_to_bin[row * w + col] != UINT16_MAX) { + spot_row = row; + spot_col = col; + break; + } + REQUIRE(spot_row > 0); + + std::vector res_mask(x.GetPixelsNum(), false); + auto settings = AdaptiveSettings(); + settings.min_pix_per_spot = 2; + settings.max_pix_per_spot = 100000; // no bound, so the footprint itself is what is measured + + // One Gaussian of width 1.5 px on a flat background of 10, at three amplitudes a hundred apart. + auto footprint = [&](double amplitude) { + ImagePreprocessorBuffer buffer(x.GetPixelsNum()); + for (size_t i = 0; i < w * h; i++) + buffer[i] = 10; + constexpr double sigma = 1.5; + for (int dr = -12; dr <= 12; dr++) + for (int dc = -12; dc <= 12; dc++) { + const double r2 = dr * dr + dc * dc; + buffer[(spot_row + dr) * w + spot_col + dc] = + 10 + static_cast(amplitude * std::exp(-r2 / (2 * sigma * sigma))); + } + AdaptiveSpotFinderCPU finder(mapping); + finder.SetResolutionMask(res_mask); + const auto spots = finder.Run(buffer, settings); + REQUIRE(spots.size() == 1); + return spots[0].PixelCount(); + }; + + const int64_t small = footprint(300.0); + const int64_t large = footprint(30000.0); + CHECK(small > 0); + // A hundredfold in peak is 4.6 e-folds. An absolute contour would add sigma^2 ln(100) ~ 10 pixels + // per e-fold of AREA here, tens of pixels in all; a peak-relative one adds nothing. + CHECK(large - small <= 4); +} diff --git a/tests/AdaptiveThresholdTest.cpp b/tests/AdaptiveThresholdTest.cpp index ecd13da64..dc5c88257 100644 --- a/tests/AdaptiveThresholdTest.cpp +++ b/tests/AdaptiveThresholdTest.cpp @@ -51,7 +51,7 @@ TEST_CASE("AdaptiveThreshold_PoissonThreshold", "[SpotFinding]") { const float z = static_cast(NormalQuantile(1.0 - p)); // The defining property: the returned count is the SMALLEST whose upper tail is within p. - for (const double mu: {1e-6, 0.1, 1.0, 3.0, 10.0, 40.0}) { + for (const double mu: {1e-6, 0.1, 1.0, 3.0, 10.0, 40.0, 60.0, 120.0}) { const int thr = static_cast(PoissonThreshold(mu, p, z)); CHECK(PoissonUpperTail(mu, thr) <= p); CHECK(PoissonUpperTail(mu, thr - 1) > p); @@ -65,8 +65,19 @@ TEST_CASE("AdaptiveThreshold_PoissonThreshold", "[SpotFinding]") { prev = thr; } - // Above mu = 50 it short-circuits to the Gaussian form mu + z sqrt(mu). - CHECK(PoissonThreshold(100.0, p, z) == Catch::Approx(100.0 + z * 10.0).epsilon(1e-5)); + // Above the summation limit the Cornish-Fisher form takes over, and it has to stay a POISSON + // quantile: the skewness term (z^2-1)/6 is what a plain mu + z*sqrt(mu) leaves out, and at this + // z the Gaussian form alone lets through several times the tail asked for. + for (const double mu: {250.0, 400.0}) { + const float thr = PoissonThreshold(mu, p, z); + CHECK(thr > mu + z * std::sqrt(mu)); + CHECK(PoissonUpperTail(mu, static_cast(thr)) <= 2 * p); + CHECK(PoissonUpperTail(mu, static_cast(mu + z * std::sqrt(mu))) > 2 * p); + } + + // And it joins the exact quantile smoothly at the switch - no step for a ring whose background + // drifts across it from frame to frame. + CHECK(std::fabs(PoissonThreshold(200.5, p, z) - PoissonThreshold(199.5, p, z)) < 2.0f); // A tighter operating point (smaller p) can only raise the threshold. CHECK(PoissonThreshold(5.0, 1e-8, static_cast(NormalQuantile(1.0 - 1e-8))) -- 2.54.0 From df98bb878951257c4baaf6e2008e6d750d379ab2 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 17:20:36 +0200 Subject: [PATCH 011/179] Space-group search: ignore the frames the crystal barely diffracted on A crystal that repeatedly leaves the beam over a full turn was read as P2(1) where it is P4(3)2(1)2, on data that merge at CC1/2 96-98% once that symmetry is assumed. The cell was right, the indexing was right - 83% of images - and every 422 operator correlation still came out between 0.08 and 0.38 against a gate of 0.30. The per-frame scale enters the intensities as 1/G. On this crystal the fitted scale spans 356x between its 5th and 95th percentiles and 15% of frames sit below a tenth of the run median, so those frames arrive amplified ten- to a hundredfold - and their sigmas are amplified by exactly the same factor, so nothing weighted by sigma can see it. What saves the PRODUCTION merge is multiplicity: in the determined symmetry a reflection is measured ten or twenty times and --reject-outliers throws the amplified observation out. The search merges in P1, where a reflection has two or three observations and no majority exists to call any of them an outlier, so the bad value lands in the merged intensity at full weight and the operator correlations pay for it. So drop those frames from the search merge only, beside the |zeta| filter and through the same snapshot that undoes it before the production merge, which keeps every frame. A tenth of the median because that is where the two populations sit: over the 39-crystal battery 33 crystals have NOT ONE frame below it, so the filter is inert on them by construction; of the six that do, five spend 0.1-3.5% of their frames there against this crystal's 15%. A twentieth leaves it under the gate; a fifth starts costing frames the battery says are real, the smallest legitimate min(G)/median(G) measured being 0.070. Battery 35/39 -> 36/39 space groups matching XDS, with 38 of the 39 rows character-identical: the only crystal that moves is this one. On it the operators go to 0.37-0.62, and every one improves in every resolution shell, most at low resolution - the signature of a scale error removed rather than noise removed. The run then merges to 1.41 A instead of 1.67, at CC1/2 97.7% and completeness 99.7%. Four other explanations were measured and refuted first: the search resolution cut (the correlations are flat in resolution, and a battery crystal with a finer detector and no cut at all gets 422 without trouble), the space group being inherited from the pass-1 geometry (handing the refined geometry in from the start changes nothing), the ordering of the corr snapshot around the collapsed-scale guard (0.02 in correlation, nothing in space group), and the ice-ring flagging (removing it entirely leaves the answer at P2(1)). Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H --- docs/CHANGELOG.md | 1 + docs/CPU_DATA_ANALYSIS.md | 2 + .../scale_merge/RotationScaleMerge.cpp | 63 ++++++++++++++++++- 3 files changed, 65 insertions(+), 1 deletion(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 37030a2aa..50e64e775 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,7 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* The de-novo space-group search ignores the frames whose fitted per-frame scale came out below a tenth of the run median. Its merge is in P1, where a reflection has too few observations for anything to catch an intensity that scale amplified; the production merge still keeps every frame. * Self-calibrating spot detection intersects its per-resolution-ring threshold with the local signal-to-noise test again, instead of replacing it. The ring threshold takes the place of the fixed photon floor and nothing else; standing alone it followed a bright reflection's skirt outwards, so on a strongly diffracting crystal the brightest reflections were detected as 100-500 pixel blobs and then discarded for being too large. * The self-calibrating threshold no longer steps where it switches from the exact Poisson tail to a normal approximation. * A spot may hold up to 200 connected pixels rather than 50. Under the self-calibrating threshold a reflection's footprint grows with its brightness, so the old bound acted as an intensity ceiling and discarded the strongest reflections of a strongly diffracting crystal - on one such set, every one of the ten brightest on an image. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 5165d9545..68ed963a0 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -1004,6 +1004,8 @@ Several space groups may share an absence pattern exactly. Where they do, the se The Lorentz factor $\zeta$ (§8.3) governs how well a reflection can be measured, so when the spindle lies in a plane of the lattice, an operator permuting the two in-plane axes samples a different mixture of measurement qualities than one that only flips signs. The search is therefore run a second time on a merge of only the well-measured observations (`--search-min-zeta`, rotation default 0.85), both answers are reported, and **where they disagree the merge of all the observations decides**. The filter discards 40–80 % of the observations, which can starve an operator correlation the full merge confirms and can equally leave an operator confirmed that the full merge refuses, so the decision — the point group as well as the absences, which live in the weak reflections the filter removes — rests on the arm with every observation behind it. A tie (same order, different symmetry) is reported with both candidates named, for trying in molecular replacement. +Both search merges also drop the **frames whose fitted per-frame scale came out below a tenth of the run median** — the stretches where the crystal was barely in the beam. The scale enters as $1/G$, so such a frame's intensities arrive amplified tenfold or more, with their $\sigma$ amplified by the identical factor and the scale's own error nowhere in it; the production merge survives that because a reflection in the determined symmetry is measured ten or twenty times and `--reject-outliers` removes the amplified observation, but the search merges in $P1$, where a reflection has two or three observations and no majority exists to call any of them an outlier. On a crystal that repeatedly left the beam over a full turn this cost the $422$ point group outright, its operators reading $0.08$–$0.38$ on a merge that gives $\mathrm{CC}_{1/2} = 96\,\%$ in that same point group once it is assumed. Like the $\zeta$ filter, this is a filter of the search passes alone — the production merge keeps every frame. + **Centering** is accepted when the systematically-absent class is weak relative to the present one by *either* of two floor-independent tests: its mean signed $I/\sigma$ well below the present mean, *or* its rate of individually-significant reflections well below the present class's own significant rate. The second test covers weak and low-energy data, where a positive intensity floor (background and profile leakage) lifts the absent class's mean $I/\sigma$ well above zero and, when the present class is itself weak, carries the plain mean ratio past its bound; a false centering fails both tests, its absent class being as strong as the present one. When several centerings pass, they are ranked by their **net** systematic absences (absent minus violating), not the gross absent count, so a super-centering (e.g. $F$ over a true $C$) whose extra, only-half-populated absent class dilutes the strength ratio does not out-rank the correct lower centering. ### 13.2 Twinning check diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 6679fd6c0..32e744a9a 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -48,6 +48,27 @@ namespace { // failure. constexpr double MIN_CREDIBLE_SCALE_RATIO = 0.02; + // The same quantity, for the P1 merge that feeds the de-novo space-group search only. That merge + // needs a much tighter floor than the production one, because what protects the production merge + // from a badly-scaled frame is multiplicity: in the determined symmetry a reflection is measured + // ten or twenty times and --reject-outliers throws the amplified observation out, while the search + // merges in P1, where a reflection has two or three observations and no majority exists to call any + // of them an outlier. The frame's own scale error is not in its sigma either (see above), so the + // bad value lands in the merged intensity at full weight. It is the operator correlations that pay: + // on a crystal that repeatedly left the beam over a full turn, 15% of frames scaled below a tenth + // of the run median and every 422 operator read 0.08-0.38 on a merge whose own data give CC1/2 96% + // and R_meas 21-30% in that very point group once it is assumed. + // + // A tenth, because that is where the two populations sit. Over the 39-crystal rotation battery, 33 + // crystals have NOT ONE frame below a tenth of their median scale, so the filter is inert on them + // by construction; of the six that do, five spend 0.1-3.5% of their frames there and the crystal + // above spends 15%. Its next stop down, a twentieth, leaves that crystal's operators at 0.13-0.41 - + // still under the gate - and a fifth buys another 0.09 of correlation for frames the battery says + // are real (the smallest legitimate min(G)/median(G) measured is 0.070). This is a filter of the + // search pass alone: like the |zeta| filter beside it, it is undone before the production merge, + // which keeps every frame. + constexpr double SEARCH_MIN_SCALE_RATIO = 0.1; + // --- Sweep-quality diagnostic (MeasureSweepQuality) --- // A stretch is reported only when BOTH per-frame channels are down: the scale (how much the crystal // diffracted) and the CC to merge (whether what it diffracted is still usable). The CC channel is what @@ -3578,7 +3599,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st // this replaces was five passes over the whole partial array - one to scatter the download across the // 80-byte Obs, one to gather the save back out, one per filter, and one to gather the result for the // re-upload - for a value the device never stopped holding. - const bool pass_filters = (for_search && search_min_zeta > 0.0) || min_cc_for_image > 0.0; + const bool pass_filters = for_search || min_cc_for_image > 0.0; if (pass_filters) { corr_before_pass_filters.resize(partials.size()); bool saved_from_gpu = false; @@ -3629,6 +3650,46 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st n_dropped, search_min_zeta); } + // --- 2a'. Search pass only: drop the frames the crystal barely diffracted on, whose 1/G would + // amplify their intensities into a P1 merge with no multiplicity to catch them. --- + if (for_search) { + std::vector fitted; + for (int f = 0; f < n_frames; ++f) + if (frame_scaled_scratch[f] && std::isfinite(g_partial[f]) && g_partial[f] > 0.0) + fitted.push_back(g_partial[f]); + if (!fitted.empty()) { + const size_t mid = fitted.size() / 2; + std::nth_element(fitted.begin(), fitted.begin() + mid, fitted.end()); + const double floor_g = fitted[mid] * SEARCH_MIN_SCALE_RATIO; + std::vector reject(n_frames, 0); + int n_rejected = 0; + for (int f = 0; f < n_frames; ++f) + if (frame_scaled_scratch[f] && std::isfinite(g_partial[f]) && g_partial[f] > 0.0 + && g_partial[f] < floor_g) { + reject[f] = 1; + ++n_rejected; + } + if (n_rejected > 0) { + bool reject_on_gpu = false; +#ifdef JFJOCH_USE_CUDA + if (gpu_active_) { + gpu_->FilterCorrByFrame(reject.data()); + reject_on_gpu = true; + } +#endif + if (!reject_on_gpu) + ParallelChunks(static_cast(partials.size()), + ThreadsForWork(partials.size(), nthreads), [&](int lo, int hi) { + for (int i = lo; i < hi; ++i) + if (reject[partials[i].frame]) partials[i].corr = 0.0f; + }); + logger.Info("Space-group search: ignoring {} of {} frames whose scale came out below " + "1/{:.0f} of the run median (the crystal barely diffracted there)", + n_rejected, n_frames, 1.0 / SEARCH_MIN_SCALE_RATIO); + } + } + } + // --- 2b. Drop frames that do not agree with the merged reference (--min-image-cc). --- if (min_cc_for_image > 0.0) { std::vector reject(n_frames, 0); -- 2.54.0 From 49515383680bd0ef686d3e526b4a7b6022df8593 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 19:02:02 +0200 Subject: [PATCH 012/179] rugnux: measure the spot budget instead of taking it as given --max-spots was a fixed 1000, and on a rotation sweep it sets the DENOMINATOR of the per-frame acceptance gate, which admits a frame when at least a fifth of its spots index. Detections are not all reflections: on a strongly diffracting crystal with heavy solvent background there are 2372 a frame, 1416 of them on ice bands, and only 178 index - so a budget that takes essentially the whole list puts the entire frame population on the gate (median indexed fraction 0.271, tenth percentile 0.218) and the run integrates 83% of its images. The same run with a smaller budget gains 211 frames and loses none, and the frames it gains are the ones with the MOST detections. That is how a larger budget integrates fewer images. So measure it: over the first pass's validation frames, with the sweep's lattice known, tally each spot by rank as +1 if it lies on the lattice and take the budget at argmax over N of n_indexed(N) - 0.20 * n_counted(N) which rises exactly while spots at that depth index better than the gate's own floor and falls after. The 0.20 is that floor, not a new constant - it is lifted out of the function-local it already lived in. It means: as deep into the intensity-ordered list as the image is still showing reflections of THIS crystal. The argmax alone would not do. Under the null that spots index at the same rate at every depth the tally is a driftless random walk, whose maximum is positive whatever the data, so a bare argmax shortens every dataset. The budget therefore has to clear the walk's own noise: the quantity it acts on is the fall from the peak to the end of the list, which is that walk read backwards, and the reflection principle gives its null law in closed form - P(fall > z*sqrt(g(1-g)T)) = 2(1-Phi(z)). That already pays for the search over ranks, so nothing further is owed to multiple comparisons. One false cut in a thousand measurements - a twelfth of one over a 39-crystal two-pass corpus - fixes z at 3.29. The level was chosen before the rule was written and was not revisited afterwards. Battery, same build, one changed default: 36/39 space groups in both arms, none lost, and 36 of the 39 crystals BIT-IDENTICAL. Two move materially - the strong crystal by +17.4% observations and +35.9 points of CC1/2 at 1.58 A, its indexing rate 83.4 -> 95.1%, and another by +93.5% observations with completeness 76.6 -> 98.4%. The significance requirement is what makes that list clean: it removed the one crystal the unguarded rule regressed, and with it two of the five gains, whose peaks do not clear 3.29 sigma on 60 frames. The lever for those is the sample and not the threshold - power grows as the root of the frame count while the bar stays where it is. Rotation only, and the first-pass lattice search always sees the full list: a fixed --max-spots 250 would have been much cheaper to write and fails a battery crystal outright, starving the de-novo FFT into a wrong cell that indexes 7 of 60 frames. Stills never reach the code path - verified bit-identical merged reflections on a serial set - and --max-spots N still pins it, as does the library default the broker and the FPGA use. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H --- docs/CHANGELOG.md | 1 + docs/CPU_DATA_ANALYSIS.md | 14 +++++ docs/RUGNUX.md | 2 +- image_analysis/indexing/AnalyzeIndexing.cpp | 50 ++++++++++++++--- image_analysis/indexing/AnalyzeIndexing.h | 39 ++++++++++++++ rugnux/Rugnux.cpp | 45 ++++++++++++++++ rugnux/Rugnux.h | 7 +++ rugnux/rugnux_cli.cpp | 4 +- tests/SpotUtilsTest.cpp | 60 +++++++++++++++++++++ 9 files changed, 212 insertions(+), 10 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 50e64e775..9b198fa66 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -2,6 +2,7 @@ ## 1.0.0 ### 1.0.0-rc.166 * The de-novo space-group search ignores the frames whose fitted per-frame scale came out below a tenth of the run median. Its merge is in P1, where a reflection has too few observations for anything to catch an intensity that scale amplified; the production merge still keeps every frame. +* rugnux measures its own spot budget on rotation data instead of keeping the strongest 1000 spots of every image: the first pass reads how deep into an image's spot list the spots still lie on the lattice it found, and the run keeps that many. `--max-spots` still pins the budget, and stills are unchanged. * Self-calibrating spot detection intersects its per-resolution-ring threshold with the local signal-to-noise test again, instead of replacing it. The ring threshold takes the place of the fixed photon floor and nothing else; standing alone it followed a bright reflection's skirt outwards, so on a strongly diffracting crystal the brightest reflections were detected as 100-500 pixel blobs and then discarded for being too large. * The self-calibrating threshold no longer steps where it switches from the exact Poisson tail to a normal approximation. * A spot may hold up to 200 connected pixels rather than 50. Under the self-calibrating threshold a reflection's footprint grows with its brightness, so the old bound acted as an intensity ceiling and discarded the strongest reflections of a strongly diffracting crystal - on one such set, every one of the ten brightest on an image. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 68ed963a0..481c85a71 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -434,6 +434,20 @@ A spot is indexed if $\delta^2 < \tau^2$, where $\tau$ is the configured toleran For indexed spots, the reciprocal lattice point $\mathbf{p} = h\mathbf{a}^*+k\mathbf{b}^*+l\mathbf{c}^*$ is used to compute $\Delta_\mathrm{Ewald}(\mathbf{p})$ (stored as a diagnostic and later used in profile-radius estimation). +A frame is taken to be this crystal's when at least a fraction $g = 0.20$ of its in-resolution, non-ice spots index. On rotation data that decision is what admits the frame to integration, so its denominator matters: every spot handed to it that is not a reflection of this crystal argues against the frame. + +### 4.2 The spot budget + +Only the strongest `--max-spots` spots of an image are kept (`FilterSpotsByCount`), and that budget therefore sets the denominator above. Detections are not all reflections — background structure, unlisted ice and detector artefacts are found too — so a budget deeper than an image's reflections makes the test above a measurement of the background rather than of the crystal, and a *larger* budget can integrate *fewer* images. + +rugnux measures the budget instead of fixing it. With the sweep's lattice in hand, the first pass tallies the spots of a sample of frames by their rank in the intensity-ordered list: how many images carried a spot at that rank, and on how many of them it indexed. Weighting each indexed spot by $1-g$ and each unindexed one by $-g$ — the same weighing the frame test applies to the list as a whole — the running sum over ranks + +$$ E(N) = n_\mathrm{indexed}(N) - g\, n_\mathrm{counted}(N) $$ + +rises exactly while the spots at that depth lie on the lattice more often than $g$, and falls after. The budget is $\arg\max_N E(N)$. Its meaning is "as deep into the list as the image is still showing reflections of this crystal": deeper spots cannot help the frame test and can only push a frame towards rejection. On crystals whose spot lists are reflections all the way down the maximum is at the end of the list and the budget is unchanged. + +The peak has to be one. Under the null — the spots lie on the lattice at the same rate at every depth — $E$ is a driftless random walk in the counted spots, with per-spot variance $g(1-g)$, and the maximum of such a walk is positive whatever the data; an $\arg\max$ taken on its own would shorten every dataset, including one with nothing to shorten. What the budget acts on is the fall from the peak to the end of the list, $E(N^*) - E(L)$, which is the maximum of the same walk read backwards; by the reflection principle its null law is $P(\mathrm{fall} > z\sqrt{g(1-g)T}) = 2(1-\Phi(z))$ over $T$ counted spots in all, so the search over ranks is already accounted for and no further multiple-comparison correction applies. The budget is taken only where the fall clears that bar at $z = 3.29$, one false shortening in a thousand measurements; otherwise the whole list is kept. + --- ## 5. FFT indexing (unknown unit cell) diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 77508371e..e91a44de2 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -957,7 +957,7 @@ Spot finding: | `--spot-high-resolution ` | High-resolution limit for spot finding, Å. Omitted (or 0): no resolution clipping — spot finding extends as far as the detector reaches, for rotation data as well as stills | | `--spot-low-resolution ` | Low-resolution limit for spot finding, Å (default: 50; lower it, e.g. 24, to exclude the direct-beam halo on weak serial data; 0 removes the limit) | | `--min-pix-per-spot ` | Minimum connected strong pixels per spot. **If omitted, min-pix is chosen per image** (stills indexing): the frame is indexed at min-pix 3/2/1 and the one maximising indexed-spot count × indexed fraction is kept. Give an explicit value to force a fixed min-pix instead. | -| `--max-spots ` | Maximum spots kept per image (the strongest ones) and handed to indexing (default: 1000) | +| `--max-spots ` | Maximum spots kept per image (the strongest ones) and handed to indexing. **If omitted, the budget is measured** on rotation data: the first pass reads how deep into an image's spot list its spots still lie on the lattice it found, and the run keeps that many (never more than 1000). Give a value to pin it. Stills always use the fixed 1000. | | `--detect-ice-rings[=on\|off]` | Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling. Default: the master file's `detect_ice_rings`, or — where the file carries no such key — **on for rotation and off for stills** | Azimuthal integration (the radial profile behind the per-image ice-ring score): diff --git a/image_analysis/indexing/AnalyzeIndexing.cpp b/image_analysis/indexing/AnalyzeIndexing.cpp index 056b8fe95..a7adc108b 100644 --- a/image_analysis/indexing/AnalyzeIndexing.cpp +++ b/image_analysis/indexing/AnalyzeIndexing.cpp @@ -420,14 +420,13 @@ bool AnalyzeIndexing(DataMessage &message, int64_t indexing_lattice_count = 0; bool outcome = false; - // Minimum fraction of the in-resolution spots a candidate lattice must index to be accepted. - // Lowering it admits weaker/sparser crystals (more real ones on flooded XFEL frames, but also more - // spurious lattices that a downstream merge-consistency gate must remove). The gate is a stills - // notion (CrystFEL, White et al., J. Appl. Cryst. 45, 335-341 (2012)): XDS, MOSFLM and DIALS index - // once over the sweep and then integrate every frame from that lattice - none of them re-decides - // per frame whether a frame may be integrated. - constexpr float min_frac = 0.20f; - const bool lattice_fits = nspots_indexed >= std::lround(min_frac * nspots_ref); + // LATTICE_MIN_INDEXED_FRACTION is the minimum fraction of the in-resolution spots a candidate + // lattice must index to be accepted. Lowering it admits weaker/sparser crystals (more real ones on + // flooded XFEL frames, but also more spurious lattices that a downstream merge-consistency gate must + // remove). The gate is a stills notion (CrystFEL, White et al., J. Appl. Cryst. 45, 335-341 (2012)): + // XDS, MOSFLM and DIALS index once over the sweep and then integrate every frame from that lattice - + // none of them re-decides per frame whether a frame may be integrated. + const bool lattice_fits = nspots_indexed >= std::lround(LATTICE_MIN_INDEXED_FRACTION * nspots_ref); // Two different questions. "Does this frame index?" - reported as the indexing rate, and what the // rotation first pass scores candidate lattices on - needs the absolute floor too, because a // handful of spots sit on almost any lattice by chance. "Is this frame worth integrating?" needs @@ -522,3 +521,38 @@ bool AnalyzeIndexing(DataMessage &message, message.indexing_result = outcome && frame_indexes; return outcome; } + +void AddSpotBudgetEvidence(const std::vector &spots, bool index_ice_rings, + std::vector &indexed, std::vector &counted) { + const size_t n = std::min({spots.size(), indexed.size(), counted.size()}); + for (size_t i = 0; i < n; i++) { + if (!index_ice_rings && spots[i].ice_ring) + continue; + counted[i]++; + if (spots[i].indexed) + indexed[i]++; + } +} + +int64_t SpotBudgetFromEvidence(const std::vector &indexed, const std::vector &counted) { + int64_t indexed_sum = 0; + int64_t counted_sum = 0; + int64_t budget = 0; + float best = 0.0f; + float at_end = 0.0f; + for (size_t i = 0; i < std::min(indexed.size(), counted.size()); i++) { + indexed_sum += indexed[i]; + counted_sum += counted[i]; + at_end = static_cast(indexed_sum) + - LATTICE_MIN_INDEXED_FRACTION * static_cast(counted_sum); + if (at_end > best) { + best = at_end; + budget = static_cast(i) + 1; + } + } + + const float g = LATTICE_MIN_INDEXED_FRACTION; + const float noise = SPOT_BUDGET_SIGNIFICANCE_Z + * std::sqrt(g * (1.0f - g) * static_cast(counted_sum)); + return (best - at_end > noise) ? budget : 0; +} diff --git a/image_analysis/indexing/AnalyzeIndexing.h b/image_analysis/indexing/AnalyzeIndexing.h index 1e612343d..3fb2d1f1b 100644 --- a/image_analysis/indexing/AnalyzeIndexing.h +++ b/image_analysis/indexing/AnalyzeIndexing.h @@ -7,6 +7,45 @@ #include "../../common/DiffractionExperiment.h" #include "../../common/JFJochMessages.h" +// Minimum fraction of a frame's in-resolution spots that must lie on a candidate lattice for the +// frame to be that crystal's. See the frame gate in AnalyzeIndexing, which is where it is applied. +constexpr float LATTICE_MIN_INDEXED_FRACTION = 0.20f; + +// Tally one image's spots by their rank in its intensity-ordered spot list: how many images had a spot +// at that rank at all (`counted`) and on how many of them it lay on the lattice (`indexed`). Ice spots +// are skipped, as they are in the frame gate. Both are added to, and their length bounds the ranks +// considered. Counts rather than weights so that the tally is exact whatever order the images are +// summed in, which is what makes the budget below independent of the thread schedule. +void AddSpotBudgetEvidence(const std::vector &spots, bool index_ice_rings, + std::vector &indexed, std::vector &counted); + +// How far the fall from the peak must exceed the counting noise of the spots for the peak to be one. +// +// Under the null - the spots lie on the lattice at the same rate at every depth - the running sum +// below is a driftless random walk in the counted spots: each is worth 1 - g with probability g and +// -g otherwise, so its step has mean zero and variance g(1-g). The MAXIMUM of such a walk is positive +// whatever the data, so an argmax taken on its own cuts every dataset, including one with nothing to +// cut. What is acted on is the FALL from the peak to the end of the list, which is the maximum of the +// same walk read backwards from the end, and the reflection principle gives that maximum's null law +// exactly: P(fall > z sqrt(g(1-g)T)) = 2(1 - Phi(z)) over T counted spots in all. The search over the +// ranks is therefore already paid for and no further multiple-comparison correction is due. z is set +// for one false cut in a thousand measurements, which over a corpus the size of a rotation test set +// (tens of crystals, a measurement per pass) expects none at all. +constexpr float SPOT_BUDGET_SIGNIFICANCE_Z = 3.29f; // 2(1 - Phi(z)) = 0.001 + +// The spot budget those tallies support: the rank at which indexed - LATTICE_MIN_INDEXED_FRACTION * +// counted, summed over the ranks down to it, peaks. Each spot that lies on the lattice is worth +// 1 - LATTICE_MIN_INDEXED_FRACTION and each one that does not costs LATTICE_MIN_INDEXED_FRACTION - the +// same weighing the frame gate applies to a spot list as a whole - so the sum rises exactly while the +// spots at that depth are on the lattice more often than the gate's floor. Deeper than the peak they +// are not: they are no longer this crystal's reflections, and they can only push a frame towards +// rejection while adding nothing the lattice recognises. +// +// Zero - keep the whole list - when the fall from that peak to the end of the list is no larger than +// the counting noise above, which is the case whenever the spots go on lying on the lattice at the +// same rate all the way down, and the case a bare argmax gets wrong. +int64_t SpotBudgetFromEvidence(const std::vector &indexed, const std::vector &counted); + bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experiment, const CrystalLattice &latt, diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index ed1e51b63..db0f96f86 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -31,6 +31,7 @@ #include "../image_analysis/MXAnalysisWithoutFPGA.h" #include "../image_analysis/beam_stop/ShadowFinder.h" #include "../image_analysis/IndexAndRefine.h" +#include "../image_analysis/indexing/AnalyzeIndexing.h" #include "../image_analysis/geom_refinement/BeamCenterFromBackground.h" #include "../image_analysis/geom_refinement/BeamCenterFromSpots.h" #include "../image_analysis/geom_refinement/GeometryRefiner.h" @@ -996,6 +997,7 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { Logger logger("Rugnux"); const std::string base_prefix = config_.output_prefix; const auto gonio_snapshot = experiment_.GetGoniometer(); + const int64_t max_spot_count_snapshot = experiment_.GetDatasetSettings().GetMaxSpotCount(); prepass_detector_geometry_.reset(); prepass_rotation_scale_.reset(); prepass_result_.reset(); @@ -1014,6 +1016,10 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { pass1.pass_count = 2; if (cancelled_) { config_.output_prefix = base_prefix; return pass1; } if (gonio_snapshot) experiment_.Goniometer(*gonio_snapshot); // undo the pre-pass goniometer shift + // A measured spot budget is the pass's own, not the run's: give the second pass the same list to + // measure from, so the two passes cannot ratchet each other down. The canonical pass then reads + // the budget off the refined geometry, where the spots that do lie on the lattice actually do. + experiment_.MaxSpotCount(max_spot_count_snapshot); // Apply the post-refined detector geometry for the second pass, keeping the header geometry so // the run can go back to it if the refined pass turns out worse (see the quality guard below). @@ -1662,6 +1668,32 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b return count.load(); }; + // How deep into an image's intensity-ordered spot list this lattice is still being seen - the + // measured spot budget. Same frames, same per-image path as count_indexed above, but scoring the + // spots rather than the frames: each frame's spots are tallied by rank and the run keeps the + // ranks the tally supports (AddSpotBudgetEvidence / SpotBudgetFromEvidence). A budget deeper than + // that admits detections that are not this crystal's reflections, and the frame gate then counts + // them against the frame - which is how a larger budget can integrate FEWER images. + auto measure_budget = [&](IndexAndRefine &idx, const RotationIndexerResult &r) -> int64_t { + idx.ForceRotationIndexerResult(r); + prefetch_spots(validation); + std::mutex evidence_mutex; + std::vector indexed(experiment_.GetMaxSpotCount(), 0); + std::vector counted(indexed.size(), 0); + const bool index_ice_rings = experiment_.GetIndexingSettings().GetIndexIceRings(); + ParallelFor(static_cast(validation.size()), + std::min(validation.size(), config_.nthreads), [&](int i) { + DataMessage m{}; + m.number = validation[i]; + m.spots = spot_cache.at(validation[i]); + if (!idx.IndexFrameOnly(m, validation_settings)) + return; + std::unique_lock ul(evidence_mutex); + AddSpotBudgetEvidence(m.spots, index_ice_rings, indexed, counted); + }); + return SpotBudgetFromEvidence(indexed, counted); + }; + // Feed one first-pass scheme (a set of image ordinals) into its own rotation indexer, ready to // be indexed. Spots are pulled from the cache here; the FFT + refinement runs separately // (RunIndexing) so the two schemes' indexing can overlap. @@ -1922,6 +1954,19 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b logger.Info("Two-pass rotation indexing found lattice (scheme '{}': {}/{} validation frames)", best.name, best.score, static_cast(validation.size())); + // The spots in hand were already cut to the budget in force, so the measurement can only + // shorten a budget, never lengthen one - which is why RunAllPasses hands the second pass the + // list the first one started from rather than the list the first one settled on. + if (config_.measure_spot_budget) { + const int64_t budget = measure_budget(*indexer, *best.result); + if (budget > 0 && budget < experiment_.GetMaxSpotCount()) { + logger.Info("Spot budget: keeping the strongest {} spots per image " + "(measured on {} frames; was {})", budget, + static_cast(validation.size()), experiment_.GetMaxSpotCount()); + experiment_.MaxSpotCount(budget); + } + } + // Second pass: compare the de-novo lattice with pass 1's HERE, before integrating every image // with it. A markedly larger cell is the bistable supercell collapse, and a centring the group // carried over from pass 1 cannot describe is the same disagreement seen from the symmetry side. diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 4b49e0a12..1e863dc5c 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -68,6 +68,13 @@ struct ProcessConfig { CalibrationMethod calibration_method = CalibrationMethod::Rings; std::vector calibrant_ring_q; + // Measure the spot budget instead of taking experiment.GetMaxSpotCount() as given (FullAnalysis, + // rotation only; the rugnux CLI sets it unless --max-spots pinned the budget). The first pass reads + // how deep into each image's spot list the spots still lie on the lattice it found, and the run + // keeps that many - see AddSpotBudgetEvidence. Off by default, so the online receiver and the viewer + // keep the fixed budget they are configured with. + bool measure_spot_budget = false; + // Rotation indexing (FullAnalysis) bool rotation_indexing = false; bool two_pass_rotation = true; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 790e7d03a..5019e0a05 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -116,7 +116,7 @@ void print_usage() { std::cout << " --spot-false-pixels Adaptive detection operating point: expected noise pixels tolerated per frame (default: 100; implies --adaptive-spots)" << std::endl; std::cout << " --spot-high-resolution High resolution limit for spot finding. If omitted (or 0), spot finding is not clipped in resolution and extends as far as the detector reaches" << std::endl; std::cout << " --spot-low-resolution Low resolution limit for spot finding, in A (default: 50; lower it, e.g. 24, to exclude the direct-beam halo on weakly-diffracting serial data; 0 removes the limit)" << std::endl; - std::cout << " --max-spots Max spot count per image, the strongest ones, handed to indexing (default: 1000)" << std::endl; + std::cout << " --max-spots Max spot count per image, the strongest ones, handed to indexing. Default: measured from the data on rotation (as deep into each image's spot list as its spots still lie on the lattice), capped at 1000" << std::endl; std::cout << " --detect-ice-rings[=on|off] Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling. Default: the master file's setting, or - where the file says nothing - on for rotation and off for stills" << std::endl; std::cout << std::endl; @@ -2312,6 +2312,8 @@ static int RunRugnux(int argc, char **argv) { config.fit_spindle = fit_spindle; config.adaptive_integration_radius = adaptive_integration_radius; config.rotation_postrefine_geometry = rotation_postrefine_geometry; + // Measure the spot budget from the data unless the user pinned it. + config.measure_spot_budget = !max_spot_count_override.has_value(); config.rotation_scale = rotation_scale; config.rotation_indexing_image_count = rotation_indexing_image_count; config.forced_rotation_lattice = forced_rotation_lattice; diff --git a/tests/SpotUtilsTest.cpp b/tests/SpotUtilsTest.cpp index 0a804c9f0..28abb93c5 100644 --- a/tests/SpotUtilsTest.cpp +++ b/tests/SpotUtilsTest.cpp @@ -4,6 +4,7 @@ #include #include "../image_analysis/spot_finding/SpotUtils.h" +#include "../image_analysis/indexing/AnalyzeIndexing.h" TEST_CASE("FilterSpuriousHighResolutionSpots") { std::vector spots; @@ -59,3 +60,62 @@ TEST_CASE("GetResolution") { // Too few spots to have a fall-off at all. CHECK_FALSE(GetResolution(std::vector(3)).has_value()); } + +TEST_CASE("SpotBudgetFromEvidence") { + // One image's worth of spots, repeated over 60 frames as the first pass does: the first 100 index + // and the next 100 do not. Every indexed spot adds 1 - 0.2 and every unindexed one takes 0.2 away, + // so the running tally rises to rank 100 and falls after it. + std::vector spots(200); + for (size_t i = 0; i < spots.size(); i++) + spots[i].indexed = i < 100; + + constexpr int frames = 60; + std::vector indexed(spots.size(), 0), counted(spots.size(), 0); + for (int f = 0; f < frames; f++) + AddSpotBudgetEvidence(spots, false, indexed, counted); + CHECK(SpotBudgetFromEvidence(indexed, counted) == 100); + + // Spots that go on indexing all the way down: the tally never falls, so there is nothing to cut. + for (auto &s: spots) + s.indexed = true; + std::vector all_hit(spots.size(), 0), all_seen(spots.size(), 0); + for (int f = 0; f < frames; f++) + AddSpotBudgetEvidence(spots, false, all_hit, all_seen); + CHECK(SpotBudgetFromEvidence(all_hit, all_seen) == 0); + + // A budget already cut to its peak has no fall left in it, so a second measurement takes nothing + // further off: the rule does not ratchet down on repetition. + CHECK(SpotBudgetFromEvidence({indexed.begin(), indexed.begin() + 100}, + {counted.begin(), counted.begin() + 100}) == 0); + + // Ice-flagged spots take no part, so a run of them neither ends the budget nor moves it: the peak + // stays at the last indexed non-ice rank before them. + std::vector with_ice(300); + for (size_t i = 0; i < with_ice.size(); i++) { + with_ice[i].ice_ring = (i >= 100 && i < 160); + with_ice[i].indexed = i < 100; + } + std::vector ice_indexed(with_ice.size(), 0), ice_counted(with_ice.size(), 0); + for (int f = 0; f < frames; f++) + AddSpotBudgetEvidence(with_ice, false, ice_indexed, ice_counted); + CHECK(SpotBudgetFromEvidence(ice_indexed, ice_counted) == 100); + + // Nothing indexes: no rank carries evidence and there is no budget to report. + for (auto &s: spots) + s.indexed = false; + std::vector none_indexed(spots.size(), 0), none_counted(spots.size(), 0); + for (int f = 0; f < frames; f++) + AddSpotBudgetEvidence(spots, false, none_indexed, none_counted); + CHECK(SpotBudgetFromEvidence(none_indexed, none_counted) == 0); + + // The case a bare argmax gets wrong: the spots index at exactly the gate's own fraction at every + // depth, so there is no depth at which the list stops being reflections. The tally still has a + // maximum - it always does - but the fall from it is inside the counting noise, and nothing is cut. + std::vector flat(1000); + for (size_t i = 0; i < flat.size(); i++) + flat[i].indexed = (i % 5 == 0); + std::vector flat_indexed(flat.size(), 0), flat_counted(flat.size(), 0); + for (int f = 0; f < frames; f++) + AddSpotBudgetEvidence(flat, false, flat_indexed, flat_counted); + CHECK(SpotBudgetFromEvidence(flat_indexed, flat_counted) == 0); +} -- 2.54.0 From e78fbd55b7222530ed1a0b34515e2e5abd5ef334 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 20:06:54 +0200 Subject: [PATCH 013/179] Indexing: let the FFT search reach past 500 A, and stop clipping the peak background fft_max_unit_cell_A was both the default and an enforced check_max, so 500 A was the longest basis vector the FFT could ever return: FFTIndexer sizes its projected histogram from that value and the transform's last usable bin IS that length. Of the PDB's 206950 X-ray entries, 1091 (0.53%) have an axis longer than that and were unindexable by construction. The accepted range now goes to 1200 A, which leaves 8. The DEFAULT is unchanged at 500 - the histogram is sized from the value in use, so nothing pays for the wider range unless a caller asks for it. The peak picker's running-mean background was truncated at the ends of the spectrum rather than slid inward, so a peak within bg_half (~15 A) of either end - which is exactly where the longest cells sit - was judged on a one-sided background, biasing its prominence by however much the spectrum sloped there. Keep the window a constant width and slide it. Both bounds stay monotonically non-decreasing in j, so the GPU kernel's running sum is still valid. Co-Authored-By: Claude Opus 5 (1M context) --- common/IndexingSettings.cpp | 8 +++++++- common/IndexingSettings.h | 8 ++++++++ image_analysis/indexing/FFTIndexerCPU.cpp | 9 +++++++-- image_analysis/indexing/FFTIndexerGPU.cu | 13 +++++++++++-- 4 files changed, 33 insertions(+), 5 deletions(-) diff --git a/common/IndexingSettings.cpp b/common/IndexingSettings.cpp index 401e9b050..ce5ad34d6 100644 --- a/common/IndexingSettings.cpp +++ b/common/IndexingSettings.cpp @@ -44,10 +44,16 @@ IndexingSettings &IndexingSettings::Algorithm(IndexingAlgorithmEnum input) { return *this; } +// The accepted range, not the default (which stays 500 A - see IndexingSettings.h). The FFT can +// only recover a basis vector up to this length, so the ceiling is exactly the longest cell the +// indexer can ever find; at 500 it excluded 1091 of the PDB's 206950 X-ray entries (0.53%) outright. +// 1200 A leaves 8. Nothing pays for the wider range: the histogram is sized from the VALUE in use, +// and only a caller that asks for more - a given cell that needs it, or the long-axis rescue - gets +// a longer transform. IndexingSettings &IndexingSettings::FFT_MaxUnitCell_A(float input) { check_finite("FFT indexing max unit cell (A)", input); check_min("FFT indexing max unit cell (A)", input, 50); - check_max("FFT indexing max unit cell (A)", input, 500); + check_max("FFT indexing max unit cell (A)", input, fft_max_unit_cell_limit_A); fft_max_unit_cell_A = input; return *this; } diff --git a/common/IndexingSettings.h b/common/IndexingSettings.h index b387c3106..34e2d864b 100644 --- a/common/IndexingSettings.h +++ b/common/IndexingSettings.h @@ -26,6 +26,14 @@ class IndexingSettings { float max_angle_from_ewald_deg = 2.0; float unit_cell_dist_tolerance_vs_reference = 0.05; // relative static constexpr float unit_cell_angle_tolerance_deg = 5.0; // degree +public: + // The longest cell the FFT search can be asked to reach, and so the longest it can ever find: + // FFTIndexer sizes its histogram from fft_max_unit_cell_A and the transform's last usable bin IS + // that length. Callers that widen the bound (a given cell, the long-axis rescue) must clamp to + // this rather than let the setter throw - a rescue that recovers an implausible axis must not + // take the whole run down with it. + static constexpr float fft_max_unit_cell_limit_A = 1200.0; +private: int64_t indexing_threads = 4; // Threads splitting the candidate-cell refinement WITHIN one indexer call. 1 (the default) is the // right answer whenever indexers already run one per image across all workers; it is raised only diff --git a/image_analysis/indexing/FFTIndexerCPU.cpp b/image_analysis/indexing/FFTIndexerCPU.cpp index 92b0bd55f..8061b62f5 100644 --- a/image_analysis/indexing/FFTIndexerCPU.cpp +++ b/image_analysis/indexing/FFTIndexerCPU.cpp @@ -110,8 +110,13 @@ void FFTIndexerCPU::ExecuteFFT(const std::vector &coord, size_t nspots) { double len = len_coeff * static_cast(j); if (len <= static_cast(min_length_A)) continue; - const int lo = std::max(0, j - bg_half); - const int hi = std::min(out_len, j + bg_half + 1); + // Slide the background window inward at the ends rather than truncating it. A peak within + // bg_half bins of either end - which is where the LONGEST cells sit, the last usable bin + // being max_length_A itself - otherwise gets its background from a one-sided window, and + // the prominence it is judged on is biased by however much the spectrum slopes there. + int lo = j - bg_half, hi = j + bg_half + 1; + if (lo < 0) { hi = std::min(out_len, hi - lo); lo = 0; } + if (hi > out_len) { lo = std::max(0, lo - (hi - out_len)); hi = out_len; } const double background = (pref[hi] - pref[lo]) / static_cast(hi - lo); const double prominence = mag[j] - background; diff --git a/image_analysis/indexing/FFTIndexerGPU.cu b/image_analysis/indexing/FFTIndexerGPU.cu index e6b011f2a..22b0b1400 100644 --- a/image_analysis/indexing/FFTIndexerGPU.cu +++ b/image_analysis/indexing/FFTIndexerGPU.cu @@ -40,9 +40,18 @@ __global__ void calculate_fft_result( FFTResult result{.magnitude = 0.0f, .direction = i, .length = -1}; for (int j = 0; j < out_len; ++j) { - const int want_hi = min(j + bg_half, out_len - 1); + // Constant-width window, slid inward at the ends instead of truncated (see + // FFTIndexerCPU): a peak within bg_half of either end - where the LONGEST cells sit - + // otherwise gets a one-sided background and a biased prominence. Both bounds stay + // monotonically non-decreasing in j, so the running sum below is still valid. + int want_lo = j - bg_half; + int want_hi = j + bg_half; + if (want_lo < 0) { want_hi = min(out_len - 1, want_hi - want_lo); want_lo = 0; } + if (want_hi > out_len - 1) { + want_lo = max(0, want_lo - (want_hi - (out_len - 1))); + want_hi = out_len - 1; + } while (whi < want_hi) { ++whi; winsum += complex_abs(d_output[offset + whi]); } - const int want_lo = max(0, j - bg_half); while (wlo < want_lo) { winsum -= complex_abs(d_output[offset + wlo]); ++wlo; } const float len = len_coeff * static_cast(j); -- 2.54.0 From dc16a00271b97d787196cbab8195e140df87a76f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 20:07:10 +0200 Subject: [PATCH 014/179] reader: read a PILATUS miniCBF sweep natively, without libcbf Most facilities still archive rotation data as a directory of miniCBF frames, which until now had to be converted to HDF5 before rugnux could see it. Nothing in that format needs a CIF parser or a library: it is an ASCII header, four separator bytes, then one byte-offset compressed image, and every value the reader wants sits on a "# " comment line or a MIME line. MiniCBF holds the format itself - header parse and the byte-offset decoder, which is a running value with deltas stored smallest-container-first. Verified byte-exact against dxtbx on PILATUS 6M, 6M-F, 300K, silicon and CdTe sensors, and three sensor thicknesses. JFJochCBFReader is a sibling of JFJochHDF5Reader under the JFJochReader base. NAMING ANY FRAME READS ITS WHOLE SWEEP: the sweep is identified by the template (prefix + digit count) the named frame belongs to, not by "every .cbf in the directory", so a directory holding two sweeps does not splice two crystals together. Naming a directory takes the sweep with the most frames in it. Images decode on demand, one per call, so any number of workers can read at once - there is no global lock as there is on the HDF5 path, HDF5 not being thread-safe. A raw CBF carries no analysis results, so the dataset it builds is the geometry, the mask and nothing else, exactly as a plain DECTRIS file with no /entry/MX gives. Two header quirks are handled because real files have them: the sensor material is written "Silicon" where the rest of the code compares against "CdTe", and the thickness unit is sometimes omitted. Headers are not a fixed size either - one set carries 6335 bytes - so the parse runs to the binary separator rather than over a fixed prefix. Co-Authored-By: Claude Opus 5 (1M context) --- reader/CMakeLists.txt | 4 + reader/JFJochCBFReader.cpp | 235 +++++++++++++++++++++++++++++++++++++ reader/JFJochCBFReader.h | 54 +++++++++ reader/MiniCBF.cpp | 215 +++++++++++++++++++++++++++++++++ reader/MiniCBF.h | 64 ++++++++++ 5 files changed, 572 insertions(+) create mode 100644 reader/JFJochCBFReader.cpp create mode 100644 reader/JFJochCBFReader.h create mode 100644 reader/MiniCBF.cpp create mode 100644 reader/MiniCBF.h diff --git a/reader/CMakeLists.txt b/reader/CMakeLists.txt index a5168f185..5fb20d52a 100644 --- a/reader/CMakeLists.txt +++ b/reader/CMakeLists.txt @@ -2,6 +2,10 @@ ADD_LIBRARY(JFJochReader STATIC JFJochReader.cpp JFJochReader.h JFJochHDF5Reader.cpp JFJochHDF5Reader.h + MiniCBF.cpp + MiniCBF.h + JFJochCBFReader.cpp + JFJochCBFReader.h HDF5ImageLocator.cpp HDF5ImageLocator.h HDF5ImageSource.cpp diff --git a/reader/JFJochCBFReader.cpp b/reader/JFJochCBFReader.cpp new file mode 100644 index 000000000..cd404f14d --- /dev/null +++ b/reader/JFJochCBFReader.cpp @@ -0,0 +1,235 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "JFJochCBFReader.h" + +#include +#include +#include +#include +#include +#include + +#include "../common/JFJochException.h" +#include "../common/JFJochMath.h" + +namespace { + +bool HasCBFExtension(const std::filesystem::path &p) { + std::string ext = p.extension().string(); + std::transform(ext.begin(), ext.end(), ext.begin(), [](unsigned char c) { return std::tolower(c); }); + return ext == ".cbf"; +} + +// The sweep a file belongs to, as a template: everything before the trailing run of digits, the +// number of digits, and the extension. "o8_1_0042.cbf" -> {"o8_1_", 4}. A directory can hold several +// sweeps ("o8_1_*" beside "o8_2_*"), so collecting every .cbf in it would silently splice two +// crystals together; matching the template is what makes "point at any frame" safe. +struct Template { + std::string prefix; + size_t digits = 0; + + bool Matches(const std::string &name) const { + if (name.size() != prefix.size() + digits + 4) // + ".cbf" + return false; + if (name.compare(0, prefix.size(), prefix) != 0) + return false; + for (size_t i = 0; i < digits; i++) + if (!std::isdigit(static_cast(name[prefix.size() + i]))) + return false; + return true; + } +}; + +std::optional