From 380cdf3dc7c6fa4ba095d99eefeb0b6ff6bba74d Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 20 Oct 2026 18:02:10 +0200 Subject: [PATCH] rugnux stills with a given cell: one crystal form merged, distance and centre searched together, a stated distance kept - A given cell keeps only the frames whose refined cell belongs to its population (GivenCellPopulation): the peak of the refined-volume histogram nearest the given cell, out to the valley that separates it from another population, with the bin set by the data's own spread. A sample of two crystal forms indexed against one form's cell otherwise merges both; a run with one population keeps every frame. - The distance grid is run from every centre hypothesis, not only the one decided at the file's distance, and a local grid of centres follows at the distance decided. A wrong header distance and a centre several pixels off (both background estimators misread it on a Rayonix set) otherwise leave each other flat; now the cell alone finds both. - A distance given with --detector-distance is not searched. Co-Authored-By: Claude Opus 5.5 --- image_analysis/scale_merge/CMakeLists.txt | 2 + image_analysis/scale_merge/CellPopulation.cpp | 77 +++++++++++++++++++ image_analysis/scale_merge/CellPopulation.h | 19 +++++ rugnux/Rugnux.cpp | 2 + rugnux/Rugnux.h | 2 + rugnux/RugnuxScaleMerge.cpp | 26 +++++++ rugnux/RugnuxStills.cpp | 68 ++++++++++++---- rugnux/rugnux_cli.cpp | 1 + tests/CMakeLists.txt | 1 + tests/CellPopulationTest.cpp | 50 ++++++++++++ 10 files changed, 233 insertions(+), 15 deletions(-) create mode 100644 image_analysis/scale_merge/CellPopulation.cpp create mode 100644 image_analysis/scale_merge/CellPopulation.h create mode 100644 tests/CellPopulationTest.cpp diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt index bae06d3d2..0fd093890 100644 --- a/image_analysis/scale_merge/CMakeLists.txt +++ b/image_analysis/scale_merge/CMakeLists.txt @@ -11,6 +11,8 @@ ADD_LIBRARY(JFJochScaleMerge Merge.h ScaleOnTheFly.cpp ScaleOnTheFly.h + CellPopulation.cpp + CellPopulation.h StillsPartialityRefine.cpp StillsPartialityRefine.h RotationScaleMerge.cpp diff --git a/image_analysis/scale_merge/CellPopulation.cpp b/image_analysis/scale_merge/CellPopulation.cpp new file mode 100644 index 000000000..ff8f7b145 --- /dev/null +++ b/image_analysis/scale_merge/CellPopulation.cpp @@ -0,0 +1,77 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "CellPopulation.h" + +#include +#include +#include + +std::pair GivenCellPopulation(const std::vector &relative_volume) { + constexpr double everything = std::numeric_limits::infinity(); + // Below this many frames a histogram has no valleys to speak of. + constexpr size_t MIN_FRAMES = 50; + // A valley separates two populations only where it dips below this fraction of the kept peak; a + // shallower dip is noise in one population's tail. + constexpr double VALLEY_FRACTION = 0.5; + + const size_t n = relative_volume.size(); + if (n < MIN_FRAMES) + return {-everything, everything}; + + std::vector v = relative_volume; + std::sort(v.begin(), v.end()); + // The bin: twice a robust standard deviation over n^(1/3) (Freedman-Diaconis with the median + // absolute deviation in place of the interquartile range, which a second population inflates). + const double median = v[n / 2]; + std::vector dev(n); + for (size_t i = 0; i < n; i++) + dev[i] = std::fabs(v[i] - median); + std::nth_element(dev.begin(), dev.begin() + n / 2, dev.end()); + const double bin = 2.0 * 1.4826 * dev[n / 2] / std::cbrt(static_cast(n)); + if (!(bin > 0)) + return {-everything, everything}; + + const double lo = v.front(); + const int nbins = static_cast((v.back() - lo) / bin) + 1; + std::vector count(nbins, 0.0); + for (double x : v) + count[std::min(nbins - 1, static_cast((x - lo) / bin))] += 1.0; + // Three-bin running mean, so a single empty or full bin is not a valley or a peak. + std::vector smooth(nbins, 0.0); + for (int i = 0; i < nbins; i++) { + const int a = std::max(0, i - 1), b = std::min(nbins - 1, i + 1); + double s = 0; + for (int j = a; j <= b; j++) + s += count[j]; + smooth[i] = s / (b - a + 1); + } + + // The peak nearest the given cell: from its bin, uphill. + int peak = std::clamp(static_cast(std::floor(-lo / bin)), 0, nbins - 1); + while (true) { + if (peak + 1 < nbins && smooth[peak + 1] > smooth[peak]) peak++; + else if (peak > 0 && smooth[peak - 1] > smooth[peak]) peak--; + else break; + } + + // Out from the peak on each side: the valley is the lowest point, below VALLEY_FRACTION of the peak, + // before the histogram rises again to twice its height by more than counting noise - the next + // population, however small. None: the population runs out to the end of the data. + const auto valley = [&](int step) { + int lowest = peak; + for (int i = peak + step; i >= 0 && i < nbins; i += step) { + if (smooth[i] < smooth[lowest]) + lowest = i; + const double low = smooth[lowest], here = smooth[i]; + // A three-bin mean of Poisson counts has a variance of a third of its value. + const bool rises = here >= 2.0 * low && here - low > 3.0 * std::sqrt((here + low) / 3.0); + if (low < VALLEY_FRACTION * smooth[peak] && rises) + return lowest; + } + return -1; + }; + const int left = valley(-1), right = valley(+1); + return {left < 0 ? -everything : lo + (left + 0.5) * bin, + right < 0 ? everything : lo + (right + 0.5) * bin}; +} diff --git a/image_analysis/scale_merge/CellPopulation.h b/image_analysis/scale_merge/CellPopulation.h new file mode 100644 index 000000000..09f326375 --- /dev/null +++ b/image_analysis/scale_merge/CellPopulation.h @@ -0,0 +1,19 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +// Serial stills given a cell: which frames belong to that crystal. A sample can hold two crystal forms +// (a polymorph, a different hydration state), and a known-cell indexer fits both to the given cell's +// orientation search - the refined cells then fall into two volume populations, and merging both mixes +// two structures. The frames kept are the population the given cell sits in: the peak of the histogram +// of v / v_given - 1 nearest to 0, out to the valley on either side that separates it from another +// population. The histogram bin follows the spread of the data, so neither the +// width of a population nor its distance to the next is a fixed number. A run with one population has +// no such valley and keeps every frame. +// +// relative_volume: v / v_given - 1 per indexed frame. Returns the [low, high] range of it to keep. +std::pair GivenCellPopulation(const std::vector &relative_volume); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index ed5a51d6d..6f1db8cd2 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -849,6 +849,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // sample of strongly-indexed frames and apply it to experiment_, so everything below (azimuthal // mapping, the main index+integrate pass, scaling) uses the refined geometry. Rotation is skipped // inside (it has its own two-pass). + if (!experiment_.IsRotationIndexing()) + stills_user_cell_ = experiment_.GetUnitCell(); if (full && config_.refine_geometry.has_value()) RefineStillsGeometry(start_image, end_image, images_to_process, observer); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 23da36c49..14c37f309 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -164,6 +164,8 @@ struct ProcessConfig { // two answers - the one comparison an on-failure trigger cannot make, because a centre wrong // along the spindle does not fail. bool beam_center_check = true; + // The detector distance was given on the command line: stills do not search it (DecideStillsGeometry). + bool detector_distance_given = false; // Whether a saturation level the frames show below the file's limit replaces it (DetectSaturation): // for every input but a Jungfraujoch-written file, whose stated limit is trusted. bool detect_saturation = false; diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 661941763..49dd2d3b4 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -83,6 +83,7 @@ #include #include #include "RugnuxPipeline.h" +#include "../image_analysis/scale_merge/CellPopulation.h" using namespace rugnux_internal; @@ -512,6 +513,31 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { logger.Info("Sweep ends: left out {} frame(s) at the start and {} at the end that recorded no " "diffraction - a shutter or beam that was not on them", at_start, at_end); } + // Stills given a cell: only the frames whose refined cell belongs to its population are merged + // (see GivenCellPopulation) - a second crystal form, fitted to the same cell, is left out. + if (!experiment_.IsRotationIndexing() && stills_user_cell_) { + auto &outcomes = indexer->GetIntegrationOutcome(); + const double given = CrystalLattice(*stills_user_cell_).CalcVolume(); + std::vector rel; + for (const auto &o : outcomes) + if (!o.reflections.empty()) + rel.push_back(o.latt.CalcVolume() / given - 1.0); + const auto [low, high] = GivenCellPopulation(rel); + size_t left_out = 0; + for (auto &o : outcomes) { + if (o.reflections.empty()) + continue; + const double r = o.latt.CalcVolume() / given - 1.0; + if (r < low || r > high) { + o.reflections.clear(); + left_out++; + } + } + if (left_out > 0) + logger.Info("Given cell: {} of {} indexed frames refined to a cell outside its population " + "({:+.1f}% to {:+.1f}% in volume) and are left out of the merge", + left_out, rel.size(), 100 * low, 100 * high); + } size_t ice_flagged = 0, ice_total = 0; { auto &outcomes = indexer->GetIntegrationOutcome(); diff --git a/rugnux/RugnuxStills.cpp b/rugnux/RugnuxStills.cpp index 02d516729..a7a5c3ced 100644 --- a/rugnux/RugnuxStills.cpp +++ b/rugnux/RugnuxStills.cpp @@ -95,8 +95,6 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ return; } // The bundle anchors a known cell; without one there is nothing to refine against. - if (!stills_user_cell_) - stills_user_cell_ = experiment_.GetUnitCell(); const auto cell = stills_user_cell_; if (!cell.has_value()) { logger.Warning("--refine-geometry needs a reference unit cell (-C / -S / reference MTZ). Skipping."); @@ -446,30 +444,41 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { // and only then is the distance a hypothesis of its own: without a cell a wrong distance indexes as a // cell of the wrong size, which only the curvature far from the beam tells apart. The stills geometry // refinement takes it from there (RunPipeline). + // A distance the user gave on the command line is taken as given: a sample of mixed crystal forms + // can put more spots on one given cell at a distance that is not the detector's. float distance = file_distance; - if (experiment_.GetUnitCell().has_value()) { - float grid_best = file_distance; - int64_t grid_score = trial(x, y, file_distance, strongest); - for (float f = MIN_FACTOR; f <= MAX_FACTOR; f *= STEP) { - if (std::fabs(f - 1.0f) < 0.5f * (STEP - 1.0f)) - continue; // the file's own, scored above - const int64_t score = trial(x, y, file_distance * f, strongest); - if (score > grid_score) { - grid_score = score; - grid_best = file_distance * f; + if (experiment_.GetUnitCell().has_value() && !config_.detector_distance_given) { + // Every centre hypothesis gets its own grid: where the file's distance is wrong, the centre + // decided at it had little indexing to go on, and a wrong centre in turn leaves the grid flat. + std::vector> centres{{x, y}, {file_x, file_y}}; + centres.insert(centres.end(), candidates.begin(), candidates.end()); + float grid_best = file_distance, grid_x = x, grid_y = y; + int64_t grid_score = -1; + for (const auto &[cx, cy] : centres) + for (float f = MIN_FACTOR; f <= MAX_FACTOR; f *= STEP) { + const int64_t score = trial(cx, cy, file_distance * f, strongest); + if (score > grid_score) { + grid_score = score; + grid_best = file_distance * f; + grid_x = cx; + grid_y = cy; + } } - } - float best_distance = file_distance; + float best_distance = file_distance, best_x = x, best_y = y; int64_t best = score_here; for (const float f : {0.98f, 0.99f, 1.0f, 1.01f, 1.02f}) { - const int64_t score = trial(x, y, grid_best * f, frames); + const int64_t score = trial(grid_x, grid_y, grid_best * f, frames); if (score > best) { best = score; best_distance = grid_best * f; + best_x = grid_x; + best_y = grid_y; } } if (score_here < FILE_KEPT_ABOVE * best) { distance = best_distance; + x = best_x; + y = best_y; stills_distance_from_cell_ = true; if (stills_user_cell_) experiment_.SetUnitCell(*stills_user_cell_); @@ -482,6 +491,35 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { decide_centre(distance); } + // With a cell, also a local grid of centres about the one decided, at the distance decided: the two + // background estimators can both be several pixels off (a shadowed or uneven background), and the + // geometry refinement only moves the centre a few pixels. Coarse, then fine about the coarse best. + if (experiment_.GetUnitCell().has_value()) { + float gx = x, gy = y; + int64_t g_best = trial(x, y, distance, strongest); + for (const auto &[step, reach] : {std::pair{3.0f, 4}, std::pair{1.0f, 2}}) { + const float cx0 = gx, cy0 = gy; + for (int i = -reach; i <= reach; i++) + for (int j = -reach; j <= reach; j++) { + if (i == 0 && j == 0) + continue; + const int64_t score = trial(cx0 + step * i, cy0 + step * j, distance, strongest); + if (score > g_best) { + g_best = score; + gx = cx0 + step * i; + gy = cy0 + step * j; + } + } + } + const int64_t here = trial(x, y, distance, frames), there = trial(gx, gy, distance, frames); + logger.Info("Stills beam centre grid at {:.2f} mm: ({:.2f},{:.2f}) puts {} sampled spots on the lattice, " + "({:.2f},{:.2f}) {}", distance, gx, gy, there, x, y, here); + if (here < FILE_KEPT_ABOVE * there) { + x = gx; + y = gy; + } + } + logger.Info("Stills beam centre: the run uses ({:.2f},{:.2f}){}", x, y, x == file_x && y == file_y ? ", the file's" : ""); experiment_.BeamX_pxl(x).BeamY_pxl(y).DetectorDistance_mm(distance); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 1dcf98c1c..ff8363290 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -3207,6 +3207,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { config.two_pass_rotation = two_pass_rotation; config.detect_beam_stop = detect_beam_stop; config.beam_center_check = beam_center_check; + config.detector_distance_given = detector_distance_mm.has_value(); config.detect_saturation = !saturation_value && (input_is_cbf || input_is_marccd || input_is_smv || input_is_crystfel || input_is_cxi || input_is_xtc); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b440124ba..9d1805532 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -12,6 +12,7 @@ ADD_EXECUTABLE(jfjoch_portable_test EXCLUDE_FROM_ALL JFJochReaderTest.cpp CrystFELGeometryTest.cpp SaturationCeilingTest.cpp + CellPopulationTest.cpp SweepLayoutTest.cpp CBORTest.cpp ImageSpotFinderCPUTest.cpp diff --git a/tests/CellPopulationTest.cpp b/tests/CellPopulationTest.cpp new file mode 100644 index 000000000..e7fc0f248 --- /dev/null +++ b/tests/CellPopulationTest.cpp @@ -0,0 +1,50 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include + +#include "../image_analysis/scale_merge/CellPopulation.h" + +TEST_CASE("CellPopulation_OnePopulationKeepsEverything", "[CellPopulation][portable]") { + std::mt19937 rng(7); + std::normal_distribution d(0.004, 0.01); + std::vector r(2000); + for (auto &x : r) x = d(rng); + const auto [low, high] = GivenCellPopulation(r); + for (double x : r) + CHECK((x >= low && x <= high)); +} + +TEST_CASE("CellPopulation_SecondFormLeftOut", "[CellPopulation][portable]") { + // Two crystal forms 9% apart in volume, the given cell one of them; whichever is given, its own + // population is kept and the other left out. + std::mt19937 rng(11); + std::normal_distribution small(0.0, 0.012), large(0.09, 0.008); + std::vector r; + for (int i = 0; i < 450; i++) r.push_back(small(rng)); + for (int i = 0; i < 150; i++) r.push_back(large(rng)); + + // Kept frames of each form, against either cell. + const auto kept = [](const std::vector &x, size_t from, size_t to) { + const auto [low, high] = GivenCellPopulation(x); + size_t n = 0; + for (size_t i = from; i < to; i++) + n += (x[i] >= low && x[i] <= high); + return n; + }; + CHECK(kept(r, 0, 450) > 440); + CHECK(kept(r, 450, 600) < 3); + + std::vector r2; + for (double x : r) r2.push_back((1.0 + x) / 1.09 - 1.0); + CHECK(kept(r2, 450, 600) > 145); + CHECK(kept(r2, 0, 450) < 5); +} + +TEST_CASE("CellPopulation_TooFewFrames", "[CellPopulation][portable]") { + const auto [low, high] = GivenCellPopulation({0.0, 0.1, 0.2}); + CHECK(std::isinf(low)); + CHECK(std::isinf(high)); +}