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 <noreply@anthropic.com>

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

image_analysis/scale_merge/CellPopulation.cpp

@@ -0,0 +1,77 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "CellPopulation.h"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+
+std::pair<double, double> GivenCellPopulation(const std::vector<double> &relative_volume) {
+    constexpr double everything = std::numeric_limits<double>::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<double> 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<double> 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<double>(n));
+    if (!(bin > 0))
+        return {-everything, everything};
+
+    const double lo = v.front();
+    const int nbins = static_cast<int>((v.back() - lo) / bin) + 1;
+    std::vector<double> count(nbins, 0.0);
+    for (double x : v)
+        count[std::min(nbins - 1, static_cast<int>((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<double> 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<int>(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};
+}

image_analysis/scale_merge/CellPopulation.h

@@ -0,0 +1,19 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <utility>
+#include <vector>
+
+// 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<double, double> GivenCellPopulation(const std::vector<double> &relative_volume);

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

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;

rugnux/RugnuxScaleMerge.cpp

@@ -83,6 +83,7 @@
 #include <array>
 #include <map>
 #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<double> 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();

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<std::pair<float, float>> 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);

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

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

tests/CellPopulationTest.cpp

@@ -0,0 +1,50 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <random>
+
+#include <catch2/catch_all.hpp>
+
+#include "../image_analysis/scale_merge/CellPopulation.h"
+
+TEST_CASE("CellPopulation_OnePopulationKeepsEverything", "[CellPopulation][portable]") {
+    std::mt19937 rng(7);
+    std::normal_distribution<double> d(0.004, 0.01);
+    std::vector<double> 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<double> small(0.0, 0.012), large(0.09, 0.008);
+    std::vector<double> 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<double> &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<double> 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));
+}