Merge fixed-cost into stills-easy: stills geometry and pre-scan run in parallel, no output change

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu

image_analysis/geom_refinement/GeometryRefiner.cpp

@@ -9,6 +9,7 @@
 #include <vector>
 
 #include "../../common/JFJochMath.h"   // PI
+#include "../../common/ParallelFor.h"
 #include "LatticeReduction.h"
 #include "XtalResidual.h"
 #include "ceres/ceres.h"
@@ -560,7 +561,7 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge
 std::vector<ModuleScan> ScanModuleShifts(const DiffractionGeometry &geom,
                                          const std::vector<GeomRefineFrame> &frames,
                                          const std::vector<ModuleShift> &modules,
-                                         double range_px) {
+                                         double range_px, size_t nthreads) {
     // A spot is on its lattice when every fractional index is within this of a whole number.
     constexpr double ON_LATTICE = 0.15;
 
@@ -602,14 +603,21 @@ std::vector<ModuleScan> ScanModuleShifts(const DiffractionGeometry &geom,
             const double step = pass == 0 ? 1.0 : 0.5;
             const double half = pass == 0 ? range_px : 1.0;
             const double cx = r.dx, cy = r.dy;
-            for (double dy = -half; dy <= half; dy += step)
-                for (double dx = -half; dx <= half; dx += step) {
-                    const int n = score(on_module[m], modules[m], cx + dx, cy + dy);
-                    if (n > r.on_lattice_at_best) {
-                        r.on_lattice_at_best = n;
-                        r.dx = cx + dx;
-                        r.dy = cy + dy;
-                    }
+            std::vector<double> offsets;
+            for (double d = -half; d <= half; d += step)
+                offsets.push_back(d);
+            // Every shift of the grid is scored on the workers, then the grid is walked in order: the first
+            // shift with the most spots wins, as it would on one thread.
+            const int n = static_cast<int>(offsets.size());
+            std::vector<int> on_lattice(n * n);
+            ParallelFor(n * n, nthreads, [&](int k) {
+                on_lattice[k] = score(on_module[m], modules[m], cx + offsets[k % n], cy + offsets[k / n]);
+            });
+            for (int k = 0; k < n * n; k++)
+                if (on_lattice[k] > r.on_lattice_at_best) {
+                    r.on_lattice_at_best = on_lattice[k];
+                    r.dx = cx + offsets[k % n];
+                    r.dy = cy + offsets[k / n];
                 }
         }
     }

image_analysis/geom_refinement/GeometryRefiner.h

@@ -110,5 +110,5 @@ struct ModuleScan {
 std::vector<ModuleScan> ScanModuleShifts(const DiffractionGeometry &geom,
                                          const std::vector<GeomRefineFrame> &frames,
                                          const std::vector<ModuleShift> &modules,
-                                         double range_px);
+                                         double range_px, size_t nthreads = 1);
 

rugnux/Rugnux.cpp

@@ -864,6 +864,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
     // inside (it has its own two-pass).
     if (!experiment_.IsRotationIndexing())
         stills_user_cell_ = experiment_.GetUnitCell();
+    // The saturation histogram reads only frames and the mask: it is taken beside the geometry refinement and
+    // decided where it always was, after it (DetectSaturation).
+    std::future<saturation_ceiling::Histogram> saturation;
+    if (full && !experiment_.IsRotationIndexing() && !config_.live)
+        saturation = std::async(std::launch::async, [&] { return SampleSaturation(start_image, images_to_process); });
     if (full && config_.refine_geometry.has_value())
         RefineStillsGeometry(start_image, end_image, images_to_process, observer);
 
@@ -873,8 +878,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
     const int prepass_images = (prepass_end - start_image + config_.stride - 1) / config_.stride;
     const int early_images = experiment_.IsRotationIndexing() ? prepass_images : images_to_process;
 
-    if (full && !experiment_.IsRotationIndexing() && !config_.live)
-        DetectSaturation(start_image, images_to_process);
+    if (saturation.valid())
+        DetectSaturation(saturation.get(), images_to_process);
 
     // After any geometry refinement, so the shadow is found about the beam centre actually used,
     // and before the azimuthal mapping and the output mask below, which both read pixel_mask_.

rugnux/Rugnux.h

@@ -44,6 +44,7 @@
 #include "../image_analysis/geom_refinement/PowderCalibration.h"   // CalibrationMethod, CalibrationResult
 #include "../image_analysis/WriteReflections.h"   // OutputProvenance
 #include "../image_analysis/bragg_integration/SpotWidth.h"       // FluxCurve
+#include "../image_analysis/image_preprocessing/SaturationCeiling.h"   // saturation_ceiling::Histogram
 #include "../image_analysis/bragg_integration/SpotFootprint.h"   // FootprintSpot
 #include "DiagnosticOutput.h"                                     // PerImagePlot
 
@@ -1415,7 +1416,9 @@ class Rugnux {
     std::optional<UnitCell> stills_user_cell_;
     bool stills_distance_from_cell_ = false;
     // Stills: where the frames' pixel values pile up, against the file's saturation limit (see there).
-    void DetectSaturation(int start_image, int images_to_process);
+    static constexpr int SATURATION_SAMPLE_FRAMES = 40;
+    saturation_ceiling::Histogram SampleSaturation(int start_image, int images_to_process);
+    void DetectSaturation(const saturation_ceiling::Histogram &histogram, int images_to_process);
 
     // One full run of the pipeline (index + integrate + optional scale/merge/space-group + output).
     // write_output == false makes no outputs (ProcessResult::outputs) and streams nothing to the image sink;
@@ -1436,6 +1439,9 @@ class Rugnux {
     // How many workers with GPU engines of their own - at most `most` - one card takes: decided from
     // the card's total memory and the detector, never from what is free (RugnuxImagePass.cpp).
     [[nodiscard]] int GPUWorkersPerCard(int most) const;
+    // How many workers a loop analysing whole frames on the cards runs - the image pass, and the stills
+    // samples ahead of it (RugnuxImagePass.cpp).
+    [[nodiscard]] int ImageWorkers() const;
     // The FFT indexers of a first pass's pool: two, for its two schemes at once; one for a probe pass on a
     // card too small for it beside the main line, which then indexes the schemes one after the other.
     // An FFT indexer holds 0.3-0.6 GB of the card, and which one indexes a frame does not change what it

rugnux/RugnuxImagePass.cpp

@@ -95,6 +95,21 @@ int Rugnux::GPUWorkersPerCard(int most) const {
     return static_cast<int>(std::clamp<int64_t>(fit, 1, most));
 }
 
+// Every worker submits its own kernels to a card, and a card runs out of room to accept them long before it
+// runs out of work: measured on two GPUs, the image loop's own time falls to sixteen workers and then rises
+// again, so forty-eight is slower than sixteen on a 16 Mpx set. Sixteen is where that measurement put the
+// minimum, and a later one on a single card agrees - eight leaves the loop waiting on the queue rather than
+// the card. The cap is per card because that is what the queue depth belongs to, and it only applies when -N
+// was left at its default - an explicit -N is a deliberate instruction and is obeyed, which is what a
+// previous attempt at this got wrong. A worker also holds GPU memory, so a small card gets fewer of them
+// (GPUWorkersPerCard).
+int Rugnux::ImageWorkers() const {
+    const int gpus = get_gpu_count();
+    if (!config_.nthreads_auto || gpus == 0)
+        return config_.nthreads;
+    return std::min(config_.nthreads, GPUWorkersPerCard(16) * gpus);
+}
+
 void Rugnux::ImagePass(PipelineLocals &p) {
     auto &logger = p.logger;
     auto &result = p.result;
@@ -374,22 +389,11 @@ void Rugnux::ImagePass(PipelineLocals &p) {
 
     std::function<void()> worker = per_image_analysis ? std::function<void()>(full_worker)
                                                       : std::function<void()>(azint_worker);
-    // How many workers the loop actually wants. Every one of them submits its own kernels to a card,
-    // and a card runs out of room to accept them long before it runs out of work: measured on two
-    // GPUs, the loop's own time falls to sixteen workers and then rises again, so forty-eight is
-    // slower than sixteen on a 16 Mpx set. Sixteen is where that measurement put the minimum, and a
-    // later one on a single card agrees - eight leaves the loop waiting on the queue rather than the
-    // card. The cap is per card because that is what the queue depth belongs to, and it only applies
-    // when -N was left at its default - an explicit -N is a deliberate instruction and is obeyed,
-    // which is what a previous attempt at this got wrong.
-    // Only the full-analysis worker drives a card; the azimuthal one preprocesses and integrates on
-    // the CPU and wants every thread it can have, so the cap must not reach it.
-    // A worker also holds GPU memory, so a small card gets fewer of them (GPUWorkersPerCard).
+    // How many workers the loop actually wants (ImageWorkers). Only the full-analysis worker drives a
+    // card; the azimuthal one preprocesses and integrates on the CPU and wants every thread it can have,
+    // so the cap must not reach it.
     const int gpus = get_gpu_count();
-    const int workers_per_card = gpus > 0 ? GPUWorkersPerCard(16) : 16;
-    const int image_workers = (per_image_analysis && config_.nthreads_auto && gpus > 0)
-            ? std::min(config_.nthreads, workers_per_card * gpus)
-            : config_.nthreads;
+    const int image_workers = per_image_analysis ? ImageWorkers() : config_.nthreads;
     if (image_workers < config_.nthreads)
         logger.Info("Per-image loop: {} of {} threads ({} GPUs) - past what a card can take queue "
                     "work from, more workers make it slower, not faster; pass -N to override",

rugnux/RugnuxStills.cpp

@@ -27,6 +27,7 @@
 #include <numeric>
 #include <map>
 #include <set>
+#include <tuple>
 #include <sstream>
 #include <thread>
 
@@ -233,12 +234,15 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
         }
     };
 
+    // As many workers as the image pass: each drives the card the same way, and how many there are does not
+    // change which frames are examined (the stripes above).
+    TimingMark sampling_mark(logger, "geometry refinement sample");
     std::vector<std::future<void>> futures;
-    futures.reserve(config_.nthreads);
-    for (int i = 0; i < config_.nthreads; ++i)
+    for (int i = 0; i < ImageWorkers(); ++i)
         futures.push_back(std::async(std::launch::async, worker));
     for (auto &fut : futures)
         fut.get();
+    sampling_mark.Stop();
 
     if (cancelled_)
         return;
@@ -310,8 +314,9 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
         // went unindexed: each module is first searched for over +-SCAN_RANGE_PX on the bundle's frames,
         // and the placement starts from the shift that puts the most of its spots on the lattices.
         constexpr double SCAN_RANGE_PX = 25.0;
+        TimingMark scan_mark(logger, "module scan");
         std::vector<ModuleShift> start = centres;
-        const auto scan = ScanModuleShifts(nominal, strong, centres, SCAN_RANGE_PX);
+        const auto scan = ScanModuleShifts(nominal, strong, centres, SCAN_RANGE_PX, config_.nthreads);
         for (size_t m = 0; m < start.size(); m++)
             if (scan[m].on_lattice_at_best > scan[m].on_lattice_at_zero) {
                 start[m].dx = scan[m].dx;
@@ -321,6 +326,7 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
                                 "({:+.1f}, {:+.1f}) px", m, start[m].x0, start[m].y0, scan[m].on_lattice_at_zero,
                                 scan[m].on_lattice_at_best, scan[m].dx, scan[m].dy);
             }
+        scan_mark.Stop();
         // Once the modules are placed, spots that went unindexed may sit on the lattice: the placement
         // takes every spot and lets the fit sort them.
         for (auto &f: placing) {
@@ -332,38 +338,50 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
         for (size_t i = 0; i < placing.size(); i++)
             half[i % 2].push_back(placing[i]);
 
+        // The four held-out fits and the placement on every frame do not depend on each other, and one fit
+        // keeps only five or six cores busy however many threads it is given (measured on 32): all five run
+        // side by side, each with every thread - on a busy machine that measured faster than a fifth each
+        // (JF16M: 10.9 s against 14.6 s, 18.9 s one after the other). The placement is computed whether or
+        // not the held-out residuals then call for it.
+        TimingMark placement_mark(logger, "module placement fits");
+        GeometryRefinerSettings st = gr_settings;
+        st.modules = start;
+        st.refine_modules = true;
+        auto placement = std::async(std::launch::async, [&] {
+            return RefineGlobalGeometry(nominal, *cell, placing, st);
+        });
+
         // Held-out residual of the geometry fitted on half `fit`, scored on every spot of the other half.
-        auto held_out = [&](const std::vector<ModuleShift> &start, bool refine) {
-            double sum = 0;
-            for (int fit = 0; fit < 2; fit++) {
-                GeometryRefinerSettings st = gr_settings;
-                st.modules = start;
-                st.refine_modules = refine;
-                const auto f = RefineGlobalGeometry(nominal, *cell, half[fit], st);
-                if (!f.ok)
-                    return std::numeric_limits<double>::infinity();
-                DiffractionGeometry g = nominal;
-                g.BeamX_pxl(f.beam_x_px).BeamY_pxl(f.beam_y_px).DetectorDistance_mm(f.distance_mm);
-                GeometryRefinerSettings ev = gr_settings;
-                ev.refine_beam_distance = false;
-                ev.modules = f.modules;
-                const auto e = RefineGlobalGeometry(g, f.cell, half[1 - fit], ev);
-                if (!e.ok)
-                    return std::numeric_limits<double>::infinity();
-                sum += e.capped_rms_px;
-            }
-            return sum / 2.0;
+        auto held_out = [&](const std::vector<ModuleShift> &modules, bool refine, int fit) {
+            GeometryRefinerSettings sf = gr_settings;
+            sf.modules = modules;
+            sf.refine_modules = refine;
+            const auto f = RefineGlobalGeometry(nominal, *cell, half[fit], sf);
+            if (!f.ok)
+                return std::numeric_limits<double>::infinity();
+            DiffractionGeometry g = nominal;
+            g.BeamX_pxl(f.beam_x_px).BeamY_pxl(f.beam_y_px).DetectorDistance_mm(f.distance_mm);
+            GeometryRefinerSettings ev = gr_settings;
+            ev.refine_beam_distance = false;
+            ev.modules = f.modules;
+            const auto e = RefineGlobalGeometry(g, f.cell, half[1 - fit], ev);
+            if (!e.ok)
+                return std::numeric_limits<double>::infinity();
+            return e.capped_rms_px;
         };
-        const double held = held_out(centres, false);
-        const double moved = held_out(start, true);
+        std::future<double> held_fits[2], moved_fits[2];
+        for (int fit = 0; fit < 2; fit++) {
+            held_fits[fit] = std::async(std::launch::async, held_out, std::cref(centres), false, fit);
+            moved_fits[fit] = std::async(std::launch::async, held_out, std::cref(start), true, fit);
+        }
+        const double held = (held_fits[0].get() + held_fits[1].get()) / 2.0;
+        const double moved = (moved_fits[0].get() + moved_fits[1].get()) / 2.0;
+        const GeometryRefinerResult rm = placement.get();
+        placement_mark.Stop();
         logger.Info("Module placement: {} modules; held-out spot residual (capped rms) {:.3f} px as given, "
                     "{:.3f} px with each placed", centres.size(), held, moved);
 
         if (moved < held) {
-            GeometryRefinerSettings st = gr_settings;
-            st.modules = start;
-            st.refine_modules = true;
-            const auto rm = RefineGlobalGeometry(nominal, *cell, placing, st);
             if (rm.ok) {
                 for (size_t m = 0; m < rm.modules.size(); m++) {
                     const auto &ms = rm.modules[m];
@@ -438,7 +456,7 @@ std::vector<DataMessage> Rugnux::FindStillsSampleSpots(const std::vector<int> &o
         }
     };
     std::vector<std::future<void>> futures;
-    for (int i = 0; i < config_.nthreads; ++i)
+    for (int i = 0; i < ImageWorkers(); ++i)
         futures.push_back(std::async(std::launch::async, worker));
     for (auto &f : futures)
         f.get();
@@ -524,7 +542,13 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) {
 
     const float file_x = experiment_.GetBeamX_pxl(), file_y = experiment_.GetBeamY_pxl();
     const float file_distance = experiment_.GetDetectorDistance_mm();
+    // A hypothesis is indexed once: the searches below meet the same one more than once - the decided centre
+    // is the file's whenever the file's is kept - and it gives the same count every time.
+    std::map<std::tuple<float, float, float, const std::vector<DataMessage> *>, int64_t> tried;
     const auto trial = [&](float x, float y, float distance, const std::vector<DataMessage> &on) {
+        const auto key = std::make_tuple(x, y, distance, &on);
+        if (const auto it = tried.find(key); it != tried.end())
+            return it->second;
         DiffractionExperiment e = experiment_;
         e.BeamX_pxl(x).BeamY_pxl(y).DetectorDistance_mm(distance);
         if (stills_user_cell_)
@@ -533,6 +557,7 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) {
         logger.Debug("Stills geometry trial: centre ({:.2f},{:.2f}) distance {:.2f} mm - {}/{} frames indexed, "
                      "{} of their {} spots on the lattice", x, y, distance, t.indexed, t.frames,
                      t.on_lattice, t.spots);
+        tried[key] = t.on_lattice;
         return t.on_lattice;
     };
 
@@ -661,12 +686,13 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) {
 // Jungfraujoch wrote (config_.detect_saturation) - a file that states no limit, or the ADC's full
 // range where each pixel's dark has already been taken off, otherwise integrates clipped peaks as if
 // they were measured, and the strongest reflections at low resolution are the ones it spoils.
-void Rugnux::DetectSaturation(int start_image, int images_to_process) {
+// The histogram reads only the frames and the mask, so it is taken beside the geometry refinement; the
+// limit is decided once both are done (DetectSaturation), where the refinement leaves the experiment.
+saturation_ceiling::Histogram Rugnux::SampleSaturation(int start_image, int images_to_process) {
     Logger logger("Rugnux");
     TimingMark mark(logger, "saturation ceiling");
-    constexpr int SAMPLE_FRAMES = 40;
     const std::vector<int> sample = select_equally_spaced_image_ordinals(
-            images_to_process, std::min(images_to_process, SAMPLE_FRAMES));
+            images_to_process, std::min(images_to_process, SATURATION_SAMPLE_FRAMES));
     const std::vector<uint32_t> &mask = pixel_mask_.GetMask();
 
     std::mutex m;
@@ -712,17 +738,22 @@ void Rugnux::DetectSaturation(int start_image, int images_to_process) {
         futures.push_back(std::async(std::launch::async, worker));
     for (auto &f : futures)
         f.get();
+    return total;
+}
 
+void Rugnux::DetectSaturation(const saturation_ceiling::Histogram &total, int images_to_process) {
+    Logger logger("Rugnux");
+    const int frames = std::min(images_to_process, SATURATION_SAMPLE_FRAMES);
     const int64_t stated = experiment_.GetSaturationLimit();
     const std::optional<int64_t> ceiling = saturation_ceiling::Ceiling(total);
     if (!ceiling) {
         logger.Info("Saturation: the values of {} sampled frames fall off without piling up below any level - "
-                    "nothing clips them under the file's limit of {}", sample.size(), stated);
+                    "nothing clips them under the file's limit of {}", frames, stated);
         return;
     }
     const bool adopt = config_.detect_saturation && *ceiling < stated;
     logger.Info("Saturation: the values of {} sampled frames pile up from {} - the file's limit is {}{}",
-                sample.size(), *ceiling, stated,
+                frames, *ceiling, stated,
                 adopt ? ", so pixels from there on are treated as overloaded"
                       : *ceiling < stated ? "; the stated limit is kept (written by Jungfraujoch, or given with --saturation)" : "");
     if (adopt)

tests/ModuleRefinementTest.cpp

@@ -9,27 +9,14 @@
 #include "../image_analysis/geom_refinement/GeometryRefiner.h"
 #include "../image_analysis/bragg_prediction/BraggPrediction.h"
 
-// A detector of four modules (the quadrants around the beam), each displaced by a known rigid move: the
-// bundle adjustment with module placement must find the moves, and the fit without it must be worse.
-TEST_CASE("RefineGlobalGeometry_ModulePlacement") {
-    DiffractionExperiment exp;
-    exp.IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(150);
-    const DiffractionGeometry geom = exp.GetDiffractionGeometry();
+namespace {
 
-    // Moves that hold none of what the gauge takes away: zero in sum, in turn and outwards.
-    std::vector<ModuleShift> truth(4);
-    const double centre[4][2] = {{500, 500}, {1500, 500}, {500, 1500}, {1500, 1500}};
-    const double move[4][3] = {{2.0, -2.0, 0.2}, {2.0, 2.0, -0.2}, {-2.0, -2.0, 0.1}, {-2.0, 2.0, -0.1}};
-    for (int m = 0; m < 4; m++) {
-        truth[m].x0 = centre[m][0];
-        truth[m].y0 = centre[m][1];
-        truth[m].dx = move[m][0];
-        truth[m].dy = move[m][1];
-        truth[m].rot_rad = move[m][2] * PI / 180.0;
-    }
+// Frames of a 50/60/70 A orthorhombic crystal in random orientations on a detector of four modules (the
+// quadrants around the beam at 1000, 1000), each displaced by the rigid move `truth` gives it: the spots are
+// where the image shows them, each indexed as predicted.
+std::vector<GeomRefineFrame> FramesOnMovedModules(const DiffractionExperiment &exp, const CrystalLattice &latt,
+                                                  const std::vector<ModuleShift> &truth) {
     auto module_of = [](double x, double y) { return (x < 1000 ? 0 : 1) + (y < 1000 ? 0 : 2); };
-
-    const CrystalLattice latt(50, 60, 70, 90, 90, 90);
     BraggPredictionSettings settings{.high_res_A = 2.0, .ewald_dist_cutoff = 0.001};
     std::mt19937 rng(7);
     std::uniform_real_distribution<float> u(-1.0f, 1.0f);
@@ -62,6 +49,31 @@ TEST_CASE("RefineGlobalGeometry_ModulePlacement") {
         if (frame.spots.size() >= 20)
             frames.push_back(frame);
     }
+    return frames;
+}
+
+}  // namespace
+
+// A detector of four modules (the quadrants around the beam), each displaced by a known rigid move: the
+// bundle adjustment with module placement must find the moves, and the fit without it must be worse.
+TEST_CASE("RefineGlobalGeometry_ModulePlacement") {
+    DiffractionExperiment exp;
+    exp.IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(150);
+    const DiffractionGeometry geom = exp.GetDiffractionGeometry();
+
+    // Moves that hold none of what the gauge takes away: zero in sum, in turn and outwards.
+    std::vector<ModuleShift> truth(4);
+    const double centre[4][2] = {{500, 500}, {1500, 500}, {500, 1500}, {1500, 1500}};
+    const double move[4][3] = {{2.0, -2.0, 0.2}, {2.0, 2.0, -0.2}, {-2.0, -2.0, 0.1}, {-2.0, 2.0, -0.1}};
+    for (int m = 0; m < 4; m++) {
+        truth[m].x0 = centre[m][0];
+        truth[m].y0 = centre[m][1];
+        truth[m].dx = move[m][0];
+        truth[m].dy = move[m][1];
+        truth[m].rot_rad = move[m][2] * PI / 180.0;
+    }
+    const CrystalLattice latt(50, 60, 70, 90, 90, 90);
+    const std::vector<GeomRefineFrame> frames = FramesOnMovedModules(exp, latt, truth);
     REQUIRE(frames.size() > 30);
 
     GeometryRefinerSettings st;
@@ -92,3 +104,48 @@ TEST_CASE("RefineGlobalGeometry_ModulePlacement") {
     CHECK(placed.beam_y_px == Catch::Approx(1000).margin(0.1));
     CHECK(placed.distance_mm == Catch::Approx(150).margin(0.05));
 }
+
+// The coarse search ahead of the placement: a module off by several pixels is found near where it is, and the
+// search gives the same answer on any number of threads.
+TEST_CASE("ScanModuleShifts_FindsShiftedModule") {
+    DiffractionExperiment exp;
+    exp.IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(150);
+    const DiffractionGeometry geom = exp.GetDiffractionGeometry();
+
+    std::vector<ModuleShift> truth(4), given(4);
+    const double centre[4][2] = {{500, 500}, {1500, 500}, {500, 1500}, {1500, 1500}};
+    for (int m = 0; m < 4; m++) {
+        truth[m].x0 = given[m].x0 = centre[m][0];
+        truth[m].y0 = given[m].y0 = centre[m][1];
+    }
+    truth[1].dx = 6.0;
+    truth[1].dy = -4.0;
+
+    const CrystalLattice latt(50, 60, 70, 90, 90, 90);
+    std::vector<GeomRefineFrame> frames = FramesOnMovedModules(exp, latt, truth);
+    REQUIRE(frames.size() > 30);
+    // The module's spots did not index where it was given: the search sees them as unindexed.
+    for (auto &f: frames)
+        for (auto &s: f.spots)
+            if (s.module == 1) {
+                s.h = s.k = s.l = 0;
+                f.unindexed.push_back(s);
+            }
+    for (auto &f: frames)
+        std::erase_if(f.spots, [](const GeomRefineSpot &s) { return s.module == 1; });
+
+    const auto one = ScanModuleShifts(geom, frames, given, 25.0, 1);
+    const auto many = ScanModuleShifts(geom, frames, given, 25.0, 4);
+    REQUIRE(one.size() == 4);
+    // Coarse: a spot is on its lattice within 0.15 of an index, here ~4 px of the spacing, so the search lands
+    // within that of the true shift and the placement takes it from there.
+    CHECK(one[1].dx == Catch::Approx(6.0).margin(4.0));
+    CHECK(one[1].dy == Catch::Approx(-4.0).margin(4.0));
+    CHECK(one[1].on_lattice_at_best > 2 * one[1].on_lattice_at_zero);
+    for (int m = 0; m < 4; m++) {
+        CHECK(many[m].dx == one[m].dx);
+        CHECK(many[m].dy == one[m].dy);
+        CHECK(many[m].on_lattice_at_best == one[m].on_lattice_at_best);
+        CHECK(many[m].on_lattice_at_zero == one[m].on_lattice_at_zero);
+    }
+}