diff --git a/image_analysis/geom_refinement/GeometryRefiner.cpp b/image_analysis/geom_refinement/GeometryRefiner.cpp index 642db86d3..d18ddba7f 100644 --- a/image_analysis/geom_refinement/GeometryRefiner.cpp +++ b/image_analysis/geom_refinement/GeometryRefiner.cpp @@ -9,6 +9,7 @@ #include #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 ScanModuleShifts(const DiffractionGeometry &geom, const std::vector &frames, const std::vector &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 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 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(offsets.size()); + std::vector 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]; } } } diff --git a/image_analysis/geom_refinement/GeometryRefiner.h b/image_analysis/geom_refinement/GeometryRefiner.h index cb1d79aa1..f624b4a2b 100644 --- a/image_analysis/geom_refinement/GeometryRefiner.h +++ b/image_analysis/geom_refinement/GeometryRefiner.h @@ -110,5 +110,5 @@ struct ModuleScan { std::vector ScanModuleShifts(const DiffractionGeometry &geom, const std::vector &frames, const std::vector &modules, - double range_px); + double range_px, size_t nthreads = 1); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index eaa72cb87..8efde9269 100644 --- a/rugnux/Rugnux.cpp +++ b/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; + 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_. diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index b0d485b71..d7d61a719 100644 --- a/rugnux/Rugnux.h +++ b/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 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; diff --git a/rugnux/RugnuxStills.cpp b/rugnux/RugnuxStills.cpp index 361bcd1ed..92bb1b67c 100644 --- a/rugnux/RugnuxStills.cpp +++ b/rugnux/RugnuxStills.cpp @@ -27,6 +27,7 @@ #include #include #include +#include #include #include @@ -233,12 +234,14 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ } }; + TimingMark sampling_mark(logger, "geometry refinement sample"); std::vector> futures; futures.reserve(config_.nthreads); for (int i = 0; i < config_.nthreads; ++i) futures.push_back(std::async(std::launch::async, worker)); for (auto &fut : futures) fut.get(); + sampling_mark.Stop(); if (cancelled_) return; @@ -310,8 +313,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 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 +325,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 +337,53 @@ 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 on a fifth of the threads. The placement is computed whether or not the held-out + // residuals then call for it. + TimingMark placement_mark(logger, "module placement fits"); + const int fit_threads = std::max(1, config_.nthreads / 5); + GeometryRefinerSettings st = gr_settings; + st.modules = start; + st.refine_modules = true; + st.num_threads = fit_threads; + 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 &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::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::infinity(); - sum += e.capped_rms_px; - } - return sum / 2.0; + auto held_out = [&](const std::vector &modules, bool refine, int fit) { + GeometryRefinerSettings sf = gr_settings; + sf.modules = modules; + sf.refine_modules = refine; + sf.num_threads = fit_threads; + const auto f = RefineGlobalGeometry(nominal, *cell, half[fit], sf); + if (!f.ok) + return std::numeric_limits::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; + ev.num_threads = fit_threads; + const auto e = RefineGlobalGeometry(g, f.cell, half[1 - fit], ev); + if (!e.ok) + return std::numeric_limits::infinity(); + return e.capped_rms_px; }; - const double held = held_out(centres, false); - const double moved = held_out(start, true); + std::future 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]; @@ -374,6 +394,7 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ // so every spot found and every prediction integrated from here on uses it. The image is // not touched. if (config_.module_diagnostics) { + st.num_threads = gr_settings.num_threads; // its fits run one at a time ModuleDiagnosticsInput diag{nominal, *cell, placing, centres, start, st, rm, experiment_.GetDetectorSetup().GetSensorMaterial(), experiment_.GetDetectorSetup().GetSensorThickness_um()}; @@ -524,7 +545,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 *>, int64_t> tried; const auto trial = [&](float x, float y, float distance, const std::vector &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 +560,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 +689,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 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 &mask = pixel_mask_.GetMask(); std::mutex m; @@ -712,17 +741,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 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) diff --git a/tests/ModuleRefinementTest.cpp b/tests/ModuleRefinementTest.cpp index 5b4886055..ede657447 100644 --- a/tests/ModuleRefinementTest.cpp +++ b/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 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 FramesOnMovedModules(const DiffractionExperiment &exp, const CrystalLattice &latt, + const std::vector &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 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 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 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 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 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); + } +}