From 89d827a24f70f9cf077ebf77d824bbc71036f305 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 19 Feb 2026 10:29:38 +0100 Subject: [PATCH 01/52] jfjoch_process: Add resolution limit --- image_analysis/scale_merge/ScaleAndMerge.cpp | 7 +++++++ image_analysis/scale_merge/ScaleAndMerge.h | 1 + tools/jfjoch_process.cpp | 14 ++++++++++---- 3 files changed, 18 insertions(+), 4 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 8c6c63ca..30733018 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -268,6 +268,9 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &ob if (!std::isfinite(r.I)) continue; + if (opt.d_min_limit_A > 0.0 && d < opt.d_min_limit_A) + continue; + if (!std::isfinite(r.zeta) || r.zeta <= 0.0f) continue; if (!std::isfinite(r.rlp) || r.rlp == 0.0f) @@ -582,6 +585,8 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &ob for (int h = 0; h < nhkl; ++h) { const auto d = static_cast(out.merged[h].d); if (std::isfinite(d) && d > 0.0f) { + if (opt.d_min_limit_A > 0.0 && d < static_cast(opt.d_min_limit_A)) + continue; stat_d_min = std::min(stat_d_min, d); stat_d_max = std::max(stat_d_max, d); } @@ -600,6 +605,8 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &ob for (int h = 0; h < nhkl; ++h) { const auto d = static_cast(out.merged[h].d); if (std::isfinite(d) && d > 0.0f) { + if (opt.d_min_limit_A > 0.0 && d < static_cast(opt.d_min_limit_A)) + continue; auto s = stat_shells.GetShell(d); if (s.has_value()) hkl_shell[h] = s.value(); diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h index 7cf9a3c4..439ec25d 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.h +++ b/image_analysis/scale_merge/ScaleAndMerge.h @@ -48,6 +48,7 @@ struct ScaleMergeOptions { bool smoothen_g = true; bool smoothen_mos = true; + double d_min_limit_A = 0.0; }; struct MergedReflection { diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 646392f0..208a5723 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -39,6 +39,7 @@ void print_usage(Logger &logger) { logger.Info(" -F Use FFT indexing algorithm (default: Auto)"); logger.Info(" -x No least-square beam center refinement"); logger.Info(" -d High resolution limit for spot finding (default: 1.5)"); + logger.Info(" -D High resolution limit for scaling/merging (default: 0.0; no limit)"); logger.Info(" -S Space group number"); logger.Info(" -M Scale and merge (refine mosaicity) and write scaled.hkl + image.dat"); logger.Info(" -L Use log-scaling residual"); @@ -67,6 +68,8 @@ int main(int argc, char **argv) { bool anomalous_mode = false; std::optional space_group_number; + double resolution_limit = 0.0; + bool log_residual = false; enum class MosaicityRefinementMode { None, Fixed, Image }; MosaicityRefinementMode mosaicity_refinement_mode = MosaicityRefinementMode::Image; @@ -79,7 +82,7 @@ int main(int argc, char **argv) { } int opt; - while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:MLm:A")) != -1) { + while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:MLm:AD:")) != -1) { switch (opt) { case 'o': output_prefix = optarg; @@ -96,6 +99,9 @@ int main(int argc, char **argv) { case 'v': verbose = true; break; + case 'd': + resolution_limit = atof(optarg); + break; case 'R': rotation_indexing = true; if (optarg) rotation_indexing_range = atof(optarg); @@ -106,8 +112,9 @@ int main(int argc, char **argv) { case 'x': refine_beam_center = false; break; - case 'd': + case 'D': d_high = atof(optarg); + logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_high); break; case 'S': space_group_number = atoi(optarg); @@ -425,6 +432,7 @@ int main(int argc, char **argv) { scale_opts.max_solver_time_s = 240.0; // generous cutoff for now scale_opts.log_scaling_residual = log_residual; scale_opts.merge_friedel = !anomalous_mode; + scale_opts.d_min_limit_A = d_high; switch (mosaicity_refinement_mode) { case MosaicityRefinementMode::None: @@ -450,8 +458,6 @@ int main(int argc, char **argv) { double scale_time = std::chrono::duration(scale_end - scale_start).count(); if (scale_result) { - -// ... existing code ... logger.Info("Scaling completed in {:.2f} s ({} unique reflections, {} images)", scale_time, scale_result->merged.size(), -- 2.49.1 From eefea7f36bad578f7b4407f99fc4cb30bb9712c7 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 20 Feb 2026 09:23:34 +0100 Subject: [PATCH 02/52] XtalOptimizer: Don't regularize lattices at input, as this can give problems to downstream procedures. We assume that SearchLattice gave reasonable response. Need to better handle this in the future. --- image_analysis/geom_refinement/XtalOptimizer.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp index 3809d314..0a4685df 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.cpp +++ b/image_analysis/geom_refinement/XtalOptimizer.cpp @@ -427,8 +427,6 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, const std::vector &spots, const float tolerance) { try { - data.latt.Regularize(data.crystal_system); - Coord vec0 = data.latt.Vec0(); Coord vec1 = data.latt.Vec1(); Coord vec2 = data.latt.Vec2(); @@ -635,9 +633,6 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, // Triclinic via the same generic builder data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, latt_vec2[0], latt_vec2[1], latt_vec2[2]); } - - data.latt.Regularize(data.crystal_system); - return true; } catch (...) { // Convergence problems, likely not updated -- 2.49.1 From 43c340219892901ec33c4bb00d196d30605050b8 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 20 Feb 2026 09:28:50 +0100 Subject: [PATCH 03/52] ScaleMerge: Simplify by removing log residual --- image_analysis/scale_merge/ScaleAndMerge.cpp | 95 ++------------------ image_analysis/scale_merge/ScaleAndMerge.h | 1 - tools/jfjoch_process.cpp | 25 +++--- 3 files changed, 18 insertions(+), 103 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 30733018..9aeab336 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -115,59 +115,6 @@ namespace { return key; } - /// CrystFEL-like log-scaling residual - /// - /// residual = w * [ ln(I_obs) - ln(G) - ln(partiality) - ln(lp) - ln(I_true) ] - /// - /// Only observations with I_obs > 0 should be fed in (the caller skips the rest). - /// G and I_true are constrained to be positive via Ceres lower bounds. - struct LogIntensityResidual { - LogIntensityResidual(const Reflection &r, double sigma_obs, double wedge_deg, bool refine_partiality) - : log_Iobs_(std::log(std::max(static_cast(r.I), 1e-30))), - weight_(SafeInv(sigma_obs / std::max(static_cast(r.I), 1e-30), 1.0)), - delta_phi_(r.delta_phi_deg), - log_lp_(std::log(std::max(SafeInv(static_cast(r.rlp), 1.0), 1e-30))), - half_wedge_(wedge_deg / 2.0), - c1_(r.zeta / std::sqrt(2.0)), - partiality_(r.partiality), - refine_partiality_(refine_partiality) { - } - - template - bool operator()(const T *const G, - const T *const mosaicity, - const T *const Itrue, - T *residual) const { - T partiality; - if (refine_partiality_ && mosaicity[0] != 0.0) { - const T arg_plus = T(delta_phi_ + half_wedge_) * T(c1_) / mosaicity[0]; - const T arg_minus = T(delta_phi_ - half_wedge_) * T(c1_) / mosaicity[0]; - partiality = (ceres::erf(arg_plus) - ceres::erf(arg_minus)) / T(2.0); - } else { - partiality = T(partiality_); - } - - // Clamp partiality away from zero so log is safe - const T min_p = T(1e-30); - if (partiality < min_p) - partiality = min_p; - - // ln(I_pred) = ln(G) + ln(partiality) + ln(lp) + ln(Itrue) - const T log_Ipred = ceres::log(G[0]) + ceres::log(partiality) + T(log_lp_) + ceres::log(Itrue[0]); - residual[0] = (log_Ipred - T(log_Iobs_)) * T(weight_); - return true; - } - - double log_Iobs_; - double weight_; // w_i ≈ I_obs / sigma_obs (relative weight in log-space) - double delta_phi_; - double log_lp_; - double half_wedge_; - double c1_; - double partiality_; - bool refine_partiality_; - }; - struct IntensityResidual { IntensityResidual(const Reflection &r, double sigma_obs, double wedge_deg, bool refine_partiality) : Iobs_(static_cast(r.I)), @@ -362,32 +309,15 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &ob for (const auto &o: obs) { size_t mos_slot = opt.per_image_mosaicity ? o.img_slot : 0; + auto *cost = new ceres::AutoDiffCostFunction( + new IntensityResidual(*o.r, o.sigma, opt.wedge_deg, refine_partiality)); - if (opt.log_scaling_residual) { - // Log residual requires positive I_obs - if (o.r->I <= 0.0f) - continue; - - auto *cost = new ceres::AutoDiffCostFunction( - new LogIntensityResidual(*o.r, o.sigma, opt.wedge_deg, refine_partiality)); - - problem.AddResidualBlock(cost, - nullptr, - &g[o.img_slot], - &mosaicity[mos_slot], - &Itrue[o.hkl_slot]); - is_valid_hkl_slot[o.hkl_slot] = true; - } else { - auto *cost = new ceres::AutoDiffCostFunction( - new IntensityResidual(*o.r, o.sigma, opt.wedge_deg, refine_partiality)); - - problem.AddResidualBlock(cost, - nullptr, - &g[o.img_slot], - &mosaicity[mos_slot], - &Itrue[o.hkl_slot]); - is_valid_hkl_slot[o.hkl_slot] = true; - } + problem.AddResidualBlock(cost, + nullptr, + &g[o.img_slot], + &mosaicity[mos_slot], + &Itrue[o.hkl_slot]); + is_valid_hkl_slot[o.hkl_slot] = true; } if (opt.smoothen_g) { @@ -418,15 +348,6 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &ob for (int i = 0; i < nimg; ++i) problem.SetParameterLowerBound(&g[i], 0, 1e-12); - - // For log residual, Itrue must stay positive - if (opt.log_scaling_residual) { - for (int h = 0; h < nhkl; ++h) { - if (is_valid_hkl_slot[h]) - problem.SetParameterLowerBound(&Itrue[h], 0, 1e-12); - } - } - // Mosaicity refinement + bounds if (!opt.refine_mosaicity) { for (int i = 0; i < (opt.per_image_mosaicity ? nimg : 1); ++i) diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h index 439ec25d..8d766785 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.h +++ b/image_analysis/scale_merge/ScaleAndMerge.h @@ -42,7 +42,6 @@ struct ScaleMergeOptions { // --- Optional: regularize per-image scale k towards 1 (Kabsch-like) --- bool regularize_scale_to_one = true; double scale_regularization_sigma = 0.05; - bool log_scaling_residual = false; double min_partiality_for_merge = 0.2; bool smoothen_g = true; diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 208a5723..b24b6af2 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -42,7 +42,6 @@ void print_usage(Logger &logger) { logger.Info(" -D High resolution limit for scaling/merging (default: 0.0; no limit)"); logger.Info(" -S Space group number"); logger.Info(" -M Scale and merge (refine mosaicity) and write scaled.hkl + image.dat"); - logger.Info(" -L Use log-scaling residual"); logger.Info(" -m Mosaicity refinement none|fixed|image (default: image)"); logger.Info(" -A Anomalous mode (don't merge Friedel pairs)"); } @@ -68,13 +67,11 @@ int main(int argc, char **argv) { bool anomalous_mode = false; std::optional space_group_number; - double resolution_limit = 0.0; - - bool log_residual = false; enum class MosaicityRefinementMode { None, Fixed, Image }; MosaicityRefinementMode mosaicity_refinement_mode = MosaicityRefinementMode::Image; - float d_high = 1.5; + float d_min_spot_finding = 1.5; + std::optional d_min_scale_merge; if (argc == 1) { print_usage(logger); @@ -82,7 +79,7 @@ int main(int argc, char **argv) { } int opt; - while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:MLm:AD:")) != -1) { + while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:Mm:AD:")) != -1) { switch (opt) { case 'o': output_prefix = optarg; @@ -100,7 +97,7 @@ int main(int argc, char **argv) { verbose = true; break; case 'd': - resolution_limit = atof(optarg); + d_min_spot_finding = atof(optarg); break; case 'R': rotation_indexing = true; @@ -113,8 +110,8 @@ int main(int argc, char **argv) { refine_beam_center = false; break; case 'D': - d_high = atof(optarg); - logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_high); + d_min_scale_merge = atof(optarg); + logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_min_spot_finding); break; case 'S': space_group_number = atoi(optarg); @@ -122,9 +119,6 @@ int main(int argc, char **argv) { case 'M': run_scaling = true; break; - case 'L': - log_residual = true; - break; case 'A': anomalous_mode = true; break; @@ -213,7 +207,9 @@ int main(int argc, char **argv) { SpotFindingSettings spot_settings; spot_settings.enable = true; spot_settings.indexing = true; - spot_settings.high_resolution_limit = d_high; + spot_settings.high_resolution_limit = d_min_spot_finding; + if (d_min_scale_merge > 0) + spot_settings.high_resolution_limit = d_min_spot_finding; // Initialize Analysis Components PixelMask pixel_mask = dataset->pixel_mask; @@ -430,9 +426,8 @@ int main(int argc, char **argv) { scale_opts.refine_mosaicity = true; scale_opts.max_num_iterations = 500; scale_opts.max_solver_time_s = 240.0; // generous cutoff for now - scale_opts.log_scaling_residual = log_residual; scale_opts.merge_friedel = !anomalous_mode; - scale_opts.d_min_limit_A = d_high; + scale_opts.d_min_limit_A = d_min_scale_merge.value_or(0.0); switch (mosaicity_refinement_mode) { case MosaicityRefinementMode::None: -- 2.49.1 From e6e0f5b2f0d420cfaddbc01e6f07da9d0ce7dd99 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 20 Feb 2026 11:37:59 +0100 Subject: [PATCH 04/52] ScaleAndMerge: clean-up (no log residual, no fixed mosaicity, return per-image values) --- image_analysis/IndexAndRefine.cpp | 26 +-- image_analysis/IndexAndRefine.h | 2 +- image_analysis/scale_merge/ScaleAndMerge.cpp | 220 +++++++++---------- image_analysis/scale_merge/ScaleAndMerge.h | 6 +- tools/jfjoch_process.cpp | 24 +- 5 files changed, 128 insertions(+), 150 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index fa45a846..e5b20981 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -17,6 +17,7 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool indexer_(indexer) { if (indexer && x.IsRotationIndexing()) rotation_indexer = std::make_unique(x, *indexer); + reflections.resize(x.GetImageNum()); } IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(DataMessage &msg) { @@ -198,16 +199,15 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, [](const Reflection& a, const Reflection& b) { return a.d < b.d; }); } + { + std::unique_lock ul(reflections_mutex); + reflections[msg.number] = refl_ret; // Image is not processed twice, so thread-safe in principle, but better safe than sorry :) + } + msg.reflections = std::move(refl_ret); CalcISigma(msg); CalcWilsonBFactor(msg); - - // Append reflections to the class-wide reflections buffer (thread-safe) - { - std::unique_lock ul(reflections_mutex); - reflections.insert(reflections.end(), msg.reflections.begin(), msg.reflections.end()); - } } void IndexAndRefine::ProcessImage(DataMessage &msg, @@ -242,15 +242,13 @@ std::optional IndexAndRefine::Finalize() { } std::optional IndexAndRefine::ScaleRotationData(const ScaleMergeOptions &opts) const { - std::vector snapshot; - { - std::unique_lock ul(reflections_mutex); - snapshot = reflections; // cheap copy under lock - } + size_t nrefl = 0; + for (const auto &i: reflections) + nrefl += i.size(); // Need a reasonable number of reflections to make refinement meaningful constexpr size_t kMinReflections = 20; - if (snapshot.size() < kMinReflections) + if (nrefl < kMinReflections) return std::nullopt; // Build options focused on mosaicity refinement but allow caller override @@ -267,5 +265,5 @@ std::optional IndexAndRefine::ScaleRotationData(const ScaleMer options.space_group = *sg; } - return ScaleAndMergeReflectionsCeres(snapshot, options); -} \ No newline at end of file + return ScaleAndMergeReflectionsCeres(reflections, options); +} diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 15fb2b13..a04592b8 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -45,7 +45,7 @@ class IndexAndRefine { }; mutable std::mutex reflections_mutex; - std::vector reflections; + std::vector> reflections; IndexingOutcome DetermineLatticeAndSymmetry(DataMessage &msg); void RefineGeometryIfNeeded(DataMessage &msg, IndexingOutcome &outcome); diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 9aeab336..e4f1e27b 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -187,98 +187,86 @@ namespace { }; } // namespace -ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &observations, +ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector> &observations, const ScaleMergeOptions &opt) { - ScaleMergeResult out; + const bool refine_partiality = opt.wedge_deg > 0.0; + ceres::Problem problem; struct ObsRef { const Reflection *r = nullptr; int img_id = 0; - int img_slot = -1; int hkl_slot = -1; double sigma = 0.0; }; - std::vector obs; - obs.reserve(observations.size()); + size_t nrefl = 0; + for (const auto &i: observations) + nrefl += i.size(); - std::unordered_map imgIdToSlot; - imgIdToSlot.reserve(256); + std::vector obs; + obs.reserve(nrefl); std::unordered_map hklToSlot; - hklToSlot.reserve(observations.size()); + hklToSlot.reserve(nrefl); - for (const auto &r: observations) { - const double d = SafeD(r.d); - if (!std::isfinite(d)) - continue; - if (!std::isfinite(r.I)) - continue; + std::vector g(observations.size(), 1.0); + std::vector mosaicity(observations.size(), opt.mosaicity_init_deg); + std::vector image_slot_used(observations.size(), 0); - if (opt.d_min_limit_A > 0.0 && d < opt.d_min_limit_A) - continue; + for (int i = 0; i < observations.size(); i++) { + for (const auto &r: observations[i]) { + const double d = SafeD(r.d); + if (!std::isfinite(d)) + continue; + if (!std::isfinite(r.I)) + continue; - if (!std::isfinite(r.zeta) || r.zeta <= 0.0f) - continue; - if (!std::isfinite(r.rlp) || r.rlp == 0.0f) - continue; + if (opt.d_min_limit_A > 0.0 && d < opt.d_min_limit_A) + continue; - const double sigma = SafeSigma(static_cast(r.sigma), opt.min_sigma); + if (!std::isfinite(r.zeta) || r.zeta <= 0.0f) + continue; + if (!std::isfinite(r.rlp) || r.rlp == 0.0f) + continue; - const int img_id = RoundImageId(r.image_number, opt.image_number_rounding); + const double sigma = SafeSigma(r.sigma, opt.min_sigma); - int img_slot; - { - auto it = imgIdToSlot.find(img_id); - if (it == imgIdToSlot.end()) { - img_slot = static_cast(imgIdToSlot.size()); - imgIdToSlot.emplace(img_id, img_slot); - } else { - img_slot = it->second; + const int img_id = i; + + image_slot_used[img_id] = 1; + + int hkl_slot; + try { + const HKLKey key = CanonicalizeHKLKey(r, opt); + auto it = hklToSlot.find(key); + if (it == hklToSlot.end()) { + hkl_slot = static_cast(hklToSlot.size()); + hklToSlot.emplace(key, hkl_slot); + } else { + hkl_slot = it->second; + } + } catch (...) { + continue; } - } - int hkl_slot; - try { - const HKLKey key = CanonicalizeHKLKey(r, opt); - auto it = hklToSlot.find(key); - if (it == hklToSlot.end()) { - hkl_slot = static_cast(hklToSlot.size()); - hklToSlot.emplace(key, hkl_slot); - } else { - hkl_slot = it->second; - } - } catch (...) { - continue; + ObsRef o; + o.r = &r; + o.img_id = img_id; + o.hkl_slot = hkl_slot; + o.sigma = sigma; + obs.push_back(o); } - - ObsRef o; - o.r = &r; - o.img_id = img_id; - o.img_slot = img_slot; - o.hkl_slot = hkl_slot; - o.sigma = sigma; - obs.push_back(o); } - const int nimg = static_cast(imgIdToSlot.size()); const int nhkl = static_cast(hklToSlot.size()); - out.image_scale_g.assign(nimg, 1.0); - out.image_ids.assign(nimg, 0); - - for (const auto &kv: imgIdToSlot) { - out.image_ids[kv.second] = kv.first; - } - - std::vector g(nimg, 1.0); std::vector Itrue(nhkl, 0.0); - // Mosaicity: always per-image - std::vector mosaicity(nimg, opt.mosaicity_init_deg); - for (int i = 0; i < nimg && i < opt.mosaicity_init_deg_vec.size(); ++i) { - if (opt.mosaicity_init_deg_vec[i] > 0.0) - mosaicity[i] = opt.mosaicity_init_deg_vec[i]; + for (int i = 0; i < observations.size(); i++) { + if (!image_slot_used[i]) { + mosaicity[i] = NAN; + g[i] = NAN; + } } // Initialize Itrue from per-HKL median of observed intensities @@ -301,61 +289,73 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &ob } } - ceres::Problem problem; - - const bool refine_partiality = opt.wedge_deg > 0.0; - std::vector is_valid_hkl_slot(nhkl, false); for (const auto &o: obs) { - size_t mos_slot = opt.per_image_mosaicity ? o.img_slot : 0; auto *cost = new ceres::AutoDiffCostFunction( new IntensityResidual(*o.r, o.sigma, opt.wedge_deg, refine_partiality)); - problem.AddResidualBlock(cost, nullptr, - &g[o.img_slot], - &mosaicity[mos_slot], + &g[o.img_id], + &mosaicity[o.img_id], &Itrue[o.hkl_slot]); is_valid_hkl_slot[o.hkl_slot] = true; } - - if (opt.smoothen_g) { - for (int i = 0; i < nimg - 2; ++i) { - auto *cost = new ceres::AutoDiffCostFunction( - new SmoothnessRegularizationResidual(0.05)); - - problem.AddResidualBlock(cost, nullptr, - &g[i], - &g[i + 1], - &g[i + 2]); + for (int i = 0; i < g.size(); ++i) { + if (image_slot_used[i]) { + auto *cost = new ceres::AutoDiffCostFunction( + new ScaleRegularizationResidual(0.05)); + problem.AddResidualBlock(cost, nullptr, &g[i]); } } - if (opt.smoothen_mos && opt.per_image_mosaicity) { - for (int i = 0; i < nimg - 2; ++i) { - auto *cost = new ceres::AutoDiffCostFunction( - new SmoothnessRegularizationResidual(0.05)); - problem.AddResidualBlock(cost, nullptr, - &mosaicity[i], - &mosaicity[i + 1], - &mosaicity[i + 2]); + if (opt.smoothen_g) { + for (int i = 0; i < g.size() - 2; ++i) { + if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) { + auto *cost = new ceres::AutoDiffCostFunction( + new SmoothnessRegularizationResidual(0.05)); + + problem.AddResidualBlock(cost, nullptr, + &g[i], + &g[i + 1], + &g[i + 2]); + } + } + } + + if (opt.smoothen_mos && opt.refine_mosaicity) { + for (int i = 0; i < mosaicity.size() - 2; ++i) { + if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) { + auto *cost = new ceres::AutoDiffCostFunction( + new SmoothnessRegularizationResidual(0.05)); + + problem.AddResidualBlock(cost, nullptr, + &mosaicity[i], + &mosaicity[i + 1], + &mosaicity[i + 2]); + } } } // Scaling factors must be always positive - for (int i = 0; i < nimg; ++i) - problem.SetParameterLowerBound(&g[i], 0, 1e-12); + for (int i = 0; i < g.size(); i++) { + if (image_slot_used[i]) + problem.SetParameterLowerBound(&g[i], 0, 1e-12); + } // Mosaicity refinement + bounds if (!opt.refine_mosaicity) { - for (int i = 0; i < (opt.per_image_mosaicity ? nimg : 1); ++i) - problem.SetParameterBlockConstant(&mosaicity[i]); + for (int i = 0; i < mosaicity.size(); ++i) { + if (image_slot_used[i]) + problem.SetParameterBlockConstant(&mosaicity[i]); + } } else { - for (int i = 0; i < (opt.per_image_mosaicity ? nimg : 1); ++i) { - problem.SetParameterLowerBound(&mosaicity[i], 0, opt.mosaicity_min_deg); - problem.SetParameterUpperBound(&mosaicity[i], 0, opt.mosaicity_max_deg); + for (int i = 0; i < mosaicity.size(); ++i) { + if (image_slot_used[i]) { + problem.SetParameterLowerBound(&mosaicity[i], 0, opt.mosaicity_min_deg); + problem.SetParameterUpperBound(&mosaicity[i], 0, opt.mosaicity_max_deg); + } } } @@ -375,18 +375,15 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &ob ceres::Solver::Summary summary; ceres::Solve(options, &problem, &summary); - // --- Export per-image results --- - for (int i = 0; i < nimg; ++i) - out.image_scale_g[i] = g[i]; - - out.mosaicity_deg.resize(nimg); - for (int i = 0; i < nimg; ++i) - out.mosaicity_deg[i] = opt.per_image_mosaicity ? mosaicity[i] : mosaicity[0]; + ScaleMergeResult out; + out.image_scale_g = g; + out.mosaicity_deg.resize(observations.size()); + for (int i = 0; i < out.mosaicity_deg.size(); i++) + out.mosaicity_deg[i] = mosaicity[i]; std::vector slotToHKL(nhkl); - for (const auto &kv: hklToSlot) { + for (const auto &kv: hklToSlot) slotToHKL[kv.second] = kv.first; - } out.merged.resize(nhkl); for (int h = 0; h < nhkl; ++h) { @@ -432,15 +429,14 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &ob for (const auto &o: obs) { const Reflection &r = *o.r; const double lp = SafeInv(static_cast(r.rlp), 1.0); - const double G_i = g[o.img_slot]; + const double G_i = g[o.img_id]; // Compute partiality with refined mosaicity double partiality; - const size_t mos_slot = opt.per_image_mosaicity ? o.img_slot : 0; - if (refine_partiality && mosaicity[mos_slot] > 0.0) { + if (refine_partiality && mosaicity[o.img_id] > 0.0) { const double c1 = r.zeta / std::sqrt(2.0); - const double arg_plus = (r.delta_phi_deg + half_wedge) * c1 / mosaicity[mos_slot]; - const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[mos_slot]; + const double arg_plus = (r.delta_phi_deg + half_wedge) * c1 / mosaicity[o.img_id]; + const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[o.img_id]; partiality = (std::erf(arg_plus) - std::erf(arg_minus)) / 2.0; } else { partiality = r.partiality; diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h index 8d766785..917b24fb 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.h +++ b/image_analysis/scale_merge/ScaleAndMerge.h @@ -32,7 +32,6 @@ struct ScaleMergeOptions { // --- Mosaicity (user input in degrees; internally converted to radians) --- bool refine_mosaicity = true; - bool per_image_mosaicity = true; double mosaicity_init_deg = 0.17; // ~0.003 rad double mosaicity_min_deg = 1e-3; @@ -81,14 +80,11 @@ struct MergeStatistics { struct ScaleMergeResult { std::vector merged; std::vector image_scale_g; - std::vector image_ids; - - // One mosaicity value per image (degrees). std::vector mosaicity_deg; /// Per-shell and overall merging statistics (populated after merging) MergeStatistics statistics; }; -ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& observations, +ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector>& observations, const ScaleMergeOptions& opt = {}); \ No newline at end of file diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index b24b6af2..ce4d485e 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -42,7 +42,7 @@ void print_usage(Logger &logger) { logger.Info(" -D High resolution limit for scaling/merging (default: 0.0; no limit)"); logger.Info(" -S Space group number"); logger.Info(" -M Scale and merge (refine mosaicity) and write scaled.hkl + image.dat"); - logger.Info(" -m Mosaicity refinement none|fixed|image (default: image)"); + logger.Info(" -m Mosaicity refinement none|image (default: image)"); logger.Info(" -A Anomalous mode (don't merge Friedel pairs)"); } @@ -67,7 +67,7 @@ int main(int argc, char **argv) { bool anomalous_mode = false; std::optional space_group_number; - enum class MosaicityRefinementMode { None, Fixed, Image }; + enum class MosaicityRefinementMode { None, Image }; MosaicityRefinementMode mosaicity_refinement_mode = MosaicityRefinementMode::Image; float d_min_spot_finding = 1.5; @@ -125,8 +125,6 @@ int main(int argc, char **argv) { case 'm': if (strcmp(optarg, "none") == 0) mosaicity_refinement_mode = MosaicityRefinementMode::None; - else if (strcmp(optarg, "fixed") == 0) - mosaicity_refinement_mode = MosaicityRefinementMode::Fixed; else if (strcmp(optarg, "image") == 0) mosaicity_refinement_mode = MosaicityRefinementMode::Image; else { @@ -433,13 +431,8 @@ int main(int argc, char **argv) { case MosaicityRefinementMode::None: scale_opts.refine_mosaicity = false; break; - case MosaicityRefinementMode::Fixed: - scale_opts.refine_mosaicity = true; - scale_opts.per_image_mosaicity = false; - break; case MosaicityRefinementMode::Image: scale_opts.refine_mosaicity = true; - scale_opts.per_image_mosaicity = true; break; } if (space_group) @@ -453,10 +446,8 @@ int main(int argc, char **argv) { double scale_time = std::chrono::duration(scale_end - scale_start).count(); if (scale_result) { - logger.Info("Scaling completed in {:.2f} s ({} unique reflections, {} images)", - scale_time, - scale_result->merged.size(), - scale_result->image_ids.size()); + logger.Info("Scaling completed in {:.2f} s ({} unique reflections)", + scale_time, scale_result->merged.size()); // Print resolution-shell statistics table { @@ -499,13 +490,10 @@ int main(int argc, char **argv) { logger.Error("Cannot open {} for writing", img_path); } else { img_file << "# image_id mosaicity_deg K\n"; - for (size_t i = 0; i < scale_result->image_ids.size(); ++i) { - img_file << scale_result->image_ids[i] << " " - << scale_result->mosaicity_deg[i] << " " - << scale_result->image_scale_g[i] << "\n"; + for (size_t i = 0; i < scale_result->mosaicity_deg.size(); ++i) { + img_file << i << " " << scale_result->mosaicity_deg[i] << " " << scale_result->image_scale_g[i] << "\n"; } img_file.close(); - logger.Info("Wrote {} image records to {}", scale_result->image_ids.size(), img_path); } } -- 2.49.1 From 3819fe5a86b13a4bfe55115c3a318487cf6a0c42 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 20 Feb 2026 11:47:01 +0100 Subject: [PATCH 05/52] ScaleAndMerge: add option to cluster images for refinement (faster convergence, but pays price in quality) --- image_analysis/scale_merge/ScaleAndMerge.cpp | 30 +++++++++++++------- image_analysis/scale_merge/ScaleAndMerge.h | 2 ++ 2 files changed, 22 insertions(+), 10 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index e4f1e27b..33dee5fb 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -189,6 +189,9 @@ namespace { ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector> &observations, const ScaleMergeOptions &opt) { + if (opt.image_cluster <= 0) + throw std::invalid_argument("image_cluster must be positive"); + const bool refine_partiality = opt.wedge_deg > 0.0; ceres::Problem problem; @@ -209,9 +212,11 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector hklToSlot; hklToSlot.reserve(nrefl); - std::vector g(observations.size(), 1.0); - std::vector mosaicity(observations.size(), opt.mosaicity_init_deg); - std::vector image_slot_used(observations.size(), 0); + size_t n_image_slots = observations.size() / opt.image_cluster + (observations.size() % opt.image_cluster > 0 ? 1 : 0); + + std::vector g(n_image_slots, 1.0); + std::vector mosaicity(n_image_slots, opt.mosaicity_init_deg); + std::vector image_slot_used(n_image_slots, 0); for (int i = 0; i < observations.size(); i++) { for (const auto &r: observations[i]) { @@ -231,8 +236,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector Itrue(nhkl, 0.0); - for (int i = 0; i < observations.size(); i++) { + for (int i = 0; i < n_image_slots; i++) { if (!image_slot_used[i]) { mosaicity[i] = NAN; g[i] = NAN; @@ -376,10 +380,16 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector slotToHKL(nhkl); for (const auto &kv: hklToSlot) diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h index 917b24fb..1943158a 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.h +++ b/image_analysis/scale_merge/ScaleAndMerge.h @@ -47,6 +47,8 @@ struct ScaleMergeOptions { bool smoothen_mos = true; double d_min_limit_A = 0.0; + + int64_t image_cluster = 1; }; struct MergedReflection { -- 2.49.1 From 203eefca66b1c8576819b193a35d357f1b67cdb8 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 20 Feb 2026 11:48:34 +0100 Subject: [PATCH 06/52] jfjoch_process: Save scale factors --- image_analysis/scale_merge/ScaleAndMerge.h | 4 ++-- tools/jfjoch_process.cpp | 2 ++ 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h index 1943158a..2699f898 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.h +++ b/image_analysis/scale_merge/ScaleAndMerge.h @@ -81,8 +81,8 @@ struct MergeStatistics { struct ScaleMergeResult { std::vector merged; - std::vector image_scale_g; - std::vector mosaicity_deg; + std::vector image_scale_g; + std::vector mosaicity_deg; /// Per-shell and overall merging statistics (populated after merging) MergeStatistics statistics; diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index ce4d485e..95c56fbd 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -446,6 +446,8 @@ int main(int argc, char **argv) { double scale_time = std::chrono::duration(scale_end - scale_start).count(); if (scale_result) { + end_msg.scale_factor = scale_result->image_scale_g; + logger.Info("Scaling completed in {:.2f} s ({} unique reflections)", scale_time, scale_result->merged.size()); -- 2.49.1 From c33366f3d6bd697ec1de270c3ef2c97f2c4780ab Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 20 Feb 2026 13:18:16 +0100 Subject: [PATCH 07/52] IndexAndRefine: Use mosaicity from indexing for scaling --- image_analysis/IndexAndRefine.cpp | 9 ++++++-- image_analysis/IndexAndRefine.h | 1 + image_analysis/scale_merge/ScaleAndMerge.cpp | 24 ++++++++------------ 3 files changed, 18 insertions(+), 16 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index e5b20981..6311c596 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -18,6 +18,7 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool if (indexer && x.IsRotationIndexing()) rotation_indexer = std::make_unique(x, *indexer); reflections.resize(x.GetImageNum()); + mosaicity.resize(x.GetImageNum(), NAN); } IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(DataMessage &msg) { @@ -167,8 +168,10 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, if (experiment.GetGoniometer().has_value()) { wedge_deg = experiment.GetGoniometer()->GetWedge_deg() / 2.0; - if (msg.mosaicity_deg) + if (msg.mosaicity_deg) { mos_deg = msg.mosaicity_deg.value(); + mosaicity[msg.number] = mos_deg; + } } const BraggPredictionSettings settings_prediction{ @@ -255,8 +258,10 @@ std::optional IndexAndRefine::ScaleRotationData(const ScaleMer ScaleMergeOptions options = opts; // If the experiment provides a wedge, propagate it - if (experiment.GetGoniometer().has_value()) + if (experiment.GetGoniometer().has_value()) { options.wedge_deg = experiment.GetGoniometer()->GetWedge_deg(); + options.mosaicity_init_deg_vec = mosaicity; + } // If caller left space_group unset, try to pick it from the indexed lattice if (!options.space_group.has_value()) { diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index a04592b8..4b0f893b 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -46,6 +46,7 @@ class IndexAndRefine { mutable std::mutex reflections_mutex; std::vector> reflections; + std::vector mosaicity; IndexingOutcome DetermineLatticeAndSymmetry(DataMessage &msg); void RefineGeometryIfNeeded(DataMessage &msg, IndexingOutcome &outcome); diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 33dee5fb..74f8c667 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -214,8 +214,6 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector 0 ? 1 : 0); - std::vector g(n_image_slots, 1.0); - std::vector mosaicity(n_image_slots, opt.mosaicity_init_deg); std::vector image_slot_used(n_image_slots, 0); for (int i = 0; i < observations.size(); i++) { @@ -262,17 +260,20 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector(hklToSlot.size()); - - std::vector Itrue(nhkl, 0.0); - + std::vector g(n_image_slots, 1.0); + std::vector mosaicity(n_image_slots, opt.mosaicity_init_deg); for (int i = 0; i < n_image_slots; i++) { if (!image_slot_used[i]) { mosaicity[i] = NAN; g[i] = NAN; + } else if (opt.mosaicity_init_deg_vec.size() > i && std::isfinite(opt.mosaicity_init_deg_vec[i])) { + mosaicity[i] = opt.mosaicity_init_deg_vec[i]; } } + const int nhkl = static_cast(hklToSlot.size()); + std::vector Itrue(nhkl, 0.0); + // Initialize Itrue from per-HKL median of observed intensities { std::vector > per_hkl_I(nhkl); @@ -305,6 +306,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector( @@ -320,10 +322,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector( new SmoothnessRegularizationResidual(0.05)); - problem.AddResidualBlock(cost, nullptr, - &g[i], - &g[i + 1], - &g[i + 2]); + problem.AddResidualBlock(cost, nullptr, &g[i], &g[i + 1], &g[i + 2]); } } } @@ -334,10 +333,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector( new SmoothnessRegularizationResidual(0.05)); - problem.AddResidualBlock(cost, nullptr, - &mosaicity[i], - &mosaicity[i + 1], - &mosaicity[i + 2]); + problem.AddResidualBlock(cost, nullptr, &mosaicity[i], &mosaicity[i + 1], &mosaicity[i + 2]); } } } -- 2.49.1 From ff6e02046a93db81d94ab9f8e0c814b375f8b9f7 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 20 Feb 2026 13:40:51 +0100 Subject: [PATCH 08/52] ScaleAndMerge: Fix still case (though it doesn't really work at the moment) --- image_analysis/IndexAndRefine.cpp | 2 +- image_analysis/scale_merge/ScaleAndMerge.cpp | 38 ++++++++++---------- image_analysis/scale_merge/ScaleAndMerge.h | 2 +- 3 files changed, 22 insertions(+), 20 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 6311c596..8605b232 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -258,7 +258,7 @@ std::optional IndexAndRefine::ScaleRotationData(const ScaleMer ScaleMergeOptions options = opts; // If the experiment provides a wedge, propagate it - if (experiment.GetGoniometer().has_value()) { + if (experiment.GetGoniometer().has_value() && rotation_indexer) { options.wedge_deg = experiment.GetGoniometer()->GetWedge_deg(); options.mosaicity_init_deg_vec = mosaicity; } diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 74f8c667..e8bf2f5e 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -192,7 +192,8 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector 0.0; + const bool rotation_crystallography = opt.wedge_deg.has_value(); + ceres::Problem problem; struct ObsRef { @@ -298,7 +299,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector( - new IntensityResidual(*o.r, o.sigma, opt.wedge_deg, refine_partiality)); + new IntensityResidual(*o.r, o.sigma, opt.wedge_deg.value_or(0.0), rotation_crystallography)); problem.AddResidualBlock(cost, nullptr, &g[o.img_id], @@ -315,25 +316,26 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector( + new SmoothnessRegularizationResidual(0.05)); - if (opt.smoothen_g) { - for (int i = 0; i < g.size() - 2; ++i) { - if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) { - auto *cost = new ceres::AutoDiffCostFunction( - new SmoothnessRegularizationResidual(0.05)); - - problem.AddResidualBlock(cost, nullptr, &g[i], &g[i + 1], &g[i + 2]); + problem.AddResidualBlock(cost, nullptr, &g[i], &g[i + 1], &g[i + 2]); + } } } - } - if (opt.smoothen_mos && opt.refine_mosaicity) { - for (int i = 0; i < mosaicity.size() - 2; ++i) { - if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) { - auto *cost = new ceres::AutoDiffCostFunction( - new SmoothnessRegularizationResidual(0.05)); + if (opt.smoothen_mos && opt.refine_mosaicity) { + for (int i = 0; i < mosaicity.size() - 2; ++i) { + if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) { + auto *cost = new ceres::AutoDiffCostFunction( + new SmoothnessRegularizationResidual(0.05)); - problem.AddResidualBlock(cost, nullptr, &mosaicity[i], &mosaicity[i + 1], &mosaicity[i + 2]); + problem.AddResidualBlock(cost, nullptr, &mosaicity[i], &mosaicity[i + 1], &mosaicity[i + 2]); + } } } } @@ -430,7 +432,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector 0.0) { + if (rotation_crystallography && mosaicity[o.img_id] > 0.0) { const double c1 = r.zeta / std::sqrt(2.0); const double arg_plus = (r.delta_phi_deg + half_wedge) * c1 / mosaicity[o.img_id]; const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[o.img_id]; diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h index 2699f898..764ce01c 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.h +++ b/image_analysis/scale_merge/ScaleAndMerge.h @@ -28,7 +28,7 @@ struct ScaleMergeOptions { // --- Kabsch(XDS)-style partiality model --- // Rotation range (wedge) used in partiality calculation. // Set to 0 to disable partiality correction. - double wedge_deg = 1.0; + std::optional wedge_deg; // --- Mosaicity (user input in degrees; internally converted to radians) --- bool refine_mosaicity = true; -- 2.49.1 From b5a7ae14ff335389dfab83fff4885c2d4117fd83 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 20 Feb 2026 14:12:22 +0100 Subject: [PATCH 09/52] jfjoch_tests: Fix ordering of unit cell of lysozyme --- tests/JFJochReceiverProcessingTest.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/JFJochReceiverProcessingTest.cpp b/tests/JFJochReceiverProcessingTest.cpp index 301ba0a0..79b80591 100644 --- a/tests/JFJochReceiverProcessingTest.cpp +++ b/tests/JFJochReceiverProcessingTest.cpp @@ -209,7 +209,7 @@ TEST_CASE("JFJochIntegrationTest_ZMQ_lysozyme_spot_and_index_refinement_tetragon .FilePrefix("lyso_test_refinement_tetragonal").JungfrauConvPhotonCnt(false).SetFileWriterFormat( FileWriterFormat::NXmxVDS).OverwriteExistingFiles(true) .DetectorDistance_mm(75).BeamY_pxl(1136).BeamX_pxl(1090).IncidentEnergy_keV(12.4) - .SetUnitCell(UnitCell{.a = 36.9, .b = 78.95, .c = 78.95, .alpha = 90, .beta = 90, .gamma = 90}) + .SetUnitCell(UnitCell{.a = 78.95, .b = 78.95, .c = 36.9, .alpha = 90, .beta = 90, .gamma = 90}) .GeomRefinementAlgorithm(GeomRefinementAlgorithmEnum::BeamCenter) .SpaceGroupNumber(96); experiment.SampleTemperature_K(123.0).RingCurrent_mA(115); -- 2.49.1 From 6eeb9c0449399d03f3409d8e7642d41b69358d17 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 20 Feb 2026 15:08:40 +0100 Subject: [PATCH 10/52] jfjoch_viewer: Read image scale factor --- reader/JFJochHDF5Reader.cpp | 2 ++ reader/JFJochReaderDataset.h | 2 ++ viewer/JFJochViewerDatasetInfo.cpp | 6 ++++++ 3 files changed, 10 insertions(+) diff --git a/reader/JFJochHDF5Reader.cpp b/reader/JFJochHDF5Reader.cpp index 2761dc06..8c5eef97 100644 --- a/reader/JFJochHDF5Reader.cpp +++ b/reader/JFJochHDF5Reader.cpp @@ -307,6 +307,8 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) { dataset->experiment.IndexingAlgorithm(IndexingAlgorithmEnum::FFT); else if (indexing == "ffbidx") dataset->experiment.IndexingAlgorithm(IndexingAlgorithmEnum::FFBIDX); + + dataset->scale_factor = master_file->ReadOptVector("/entry/MX/imageScaleFactor"); } auto ring_current_A = master_file->GetOptFloat("/entry/source/current"); diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h index 282d91e4..b1268f94 100644 --- a/reader/JFJochReaderDataset.h +++ b/reader/JFJochReaderDataset.h @@ -41,6 +41,8 @@ struct JFJochReaderDataset { std::vector mosaicity_deg; std::vector b_factor; + std::vector scale_factor; + std::vector max_value; std::vector roi; diff --git a/viewer/JFJochViewerDatasetInfo.cpp b/viewer/JFJochViewerDatasetInfo.cpp index 516cf33d..2029b10b 100644 --- a/viewer/JFJochViewerDatasetInfo.cpp +++ b/viewer/JFJochViewerDatasetInfo.cpp @@ -90,6 +90,10 @@ void JFJochViewerDatasetInfo::UpdateLabels() { combo_box->addItem("Mosaicity", 9); } + if (!dataset->scale_factor.empty()) { + combo_box->insertSeparator(1000); + combo_box->addItem("Scale factor", 10); + } for (int i = 0; i < this->dataset->roi.size(); i++) { std::string name = std::string("ROI ") + this->dataset->roi[i]; combo_box->insertSeparator(1000); @@ -167,6 +171,8 @@ void JFJochViewerDatasetInfo::UpdatePlot() { data = dataset->b_factor; } else if (val == 9) { data = dataset->mosaicity_deg; + } else if (val == 10) { + data = dataset->scale_factor; } else if (val >= 100) { int roi_index = (val - 100) / 4; -- 2.49.1 From 745755184ada0ab8fda9ce800f2d0721063e53e7 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 20 Feb 2026 15:27:28 +0100 Subject: [PATCH 11/52] ScaleAndMerge: Split code into functions --- image_analysis/scale_merge/ScaleAndMerge.cpp | 549 ++++++++++--------- 1 file changed, 291 insertions(+), 258 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index e8bf2f5e..90afbd75 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -185,16 +185,6 @@ namespace { double inv_sigma_; }; -} // namespace - -ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector> &observations, - const ScaleMergeOptions &opt) { - if (opt.image_cluster <= 0) - throw std::invalid_argument("image_cluster must be positive"); - - const bool rotation_crystallography = opt.wedge_deg.has_value(); - - ceres::Problem problem; struct ObsRef { const Reflection *r = nullptr; @@ -203,269 +193,130 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector obs; - obs.reserve(nrefl); - - std::unordered_map hklToSlot; - hklToSlot.reserve(nrefl); - - size_t n_image_slots = observations.size() / opt.image_cluster + (observations.size() % opt.image_cluster > 0 ? 1 : 0); - - std::vector image_slot_used(n_image_slots, 0); - - for (int i = 0; i < observations.size(); i++) { - for (const auto &r: observations[i]) { - const double d = SafeD(r.d); - if (!std::isfinite(d)) - continue; - if (!std::isfinite(r.I)) - continue; - - if (opt.d_min_limit_A > 0.0 && d < opt.d_min_limit_A) - continue; - - if (!std::isfinite(r.zeta) || r.zeta <= 0.0f) - continue; - if (!std::isfinite(r.rlp) || r.rlp == 0.0f) - continue; - - const double sigma = SafeSigma(r.sigma, opt.min_sigma); - - const int img_id = i / opt.image_cluster; - image_slot_used[img_id] = 1; - - int hkl_slot; - try { - const HKLKey key = CanonicalizeHKLKey(r, opt); - auto it = hklToSlot.find(key); - if (it == hklToSlot.end()) { - hkl_slot = static_cast(hklToSlot.size()); - hklToSlot.emplace(key, hkl_slot); - } else { - hkl_slot = it->second; - } - } catch (...) { - continue; - } - - ObsRef o; - o.r = &r; - o.img_id = img_id; - o.hkl_slot = hkl_slot; - o.sigma = sigma; - obs.push_back(o); - } - } - - std::vector g(n_image_slots, 1.0); - std::vector mosaicity(n_image_slots, opt.mosaicity_init_deg); - for (int i = 0; i < n_image_slots; i++) { - if (!image_slot_used[i]) { - mosaicity[i] = NAN; - g[i] = NAN; - } else if (opt.mosaicity_init_deg_vec.size() > i && std::isfinite(opt.mosaicity_init_deg_vec[i])) { - mosaicity[i] = opt.mosaicity_init_deg_vec[i]; - } - } - - const int nhkl = static_cast(hklToSlot.size()); - std::vector Itrue(nhkl, 0.0); - - // Initialize Itrue from per-HKL median of observed intensities - { - std::vector > per_hkl_I(nhkl); - for (const auto &o: obs) { - per_hkl_I[o.hkl_slot].push_back(static_cast(o.r->I)); - } - for (int h = 0; h < nhkl; ++h) { - auto &v = per_hkl_I[h]; - if (v.empty()) { - Itrue[h] = std::max(opt.min_sigma, 1e-6); - continue; - } - std::nth_element(v.begin(), v.begin() + static_cast(v.size() / 2), v.end()); - double med = v[v.size() / 2]; - if (!std::isfinite(med) || med <= opt.min_sigma) - med = opt.min_sigma; - Itrue[h] = med; - } - } - - std::vector is_valid_hkl_slot(nhkl, false); - - for (const auto &o: obs) { - auto *cost = new ceres::AutoDiffCostFunction( - new IntensityResidual(*o.r, o.sigma, opt.wedge_deg.value_or(0.0), rotation_crystallography)); - problem.AddResidualBlock(cost, - nullptr, - &g[o.img_id], - &mosaicity[o.img_id], - &Itrue[o.hkl_slot]); - is_valid_hkl_slot[o.hkl_slot] = true; - } - - for (int i = 0; i < g.size(); ++i) { - if (image_slot_used[i]) { - auto *cost = new ceres::AutoDiffCostFunction( - new ScaleRegularizationResidual(0.05)); - problem.AddResidualBlock(cost, nullptr, &g[i]); - } - } - - if (rotation_crystallography) { - if (opt.smoothen_g) { - for (int i = 0; i < g.size() - 2; ++i) { - if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) { - auto *cost = new ceres::AutoDiffCostFunction( - new SmoothnessRegularizationResidual(0.05)); - - problem.AddResidualBlock(cost, nullptr, &g[i], &g[i + 1], &g[i + 2]); - } - } - } - - if (opt.smoothen_mos && opt.refine_mosaicity) { - for (int i = 0; i < mosaicity.size() - 2; ++i) { - if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) { - auto *cost = new ceres::AutoDiffCostFunction( - new SmoothnessRegularizationResidual(0.05)); - - problem.AddResidualBlock(cost, nullptr, &mosaicity[i], &mosaicity[i + 1], &mosaicity[i + 2]); - } - } - } - } - - // Scaling factors must be always positive - for (int i = 0; i < g.size(); i++) { - if (image_slot_used[i]) - problem.SetParameterLowerBound(&g[i], 0, 1e-12); - } - - // Mosaicity refinement + bounds - if (!opt.refine_mosaicity) { - for (int i = 0; i < mosaicity.size(); ++i) { - if (image_slot_used[i]) - problem.SetParameterBlockConstant(&mosaicity[i]); - } - } else { - for (int i = 0; i < mosaicity.size(); ++i) { - if (image_slot_used[i]) { - problem.SetParameterLowerBound(&mosaicity[i], 0, opt.mosaicity_min_deg); - problem.SetParameterUpperBound(&mosaicity[i], 0, opt.mosaicity_max_deg); - } - } - } - - // use all available threads - unsigned int hw = std::thread::hardware_concurrency(); - if (hw == 0) - hw = 1; // fallback - - ceres::Solver::Options options; - - options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY; - options.minimizer_progress_to_stdout = true; - options.max_num_iterations = opt.max_num_iterations; - options.max_solver_time_in_seconds = opt.max_solver_time_s; - options.num_threads = static_cast(hw); - - ceres::Solver::Summary summary; - ceres::Solve(options, &problem, &summary); - - ScaleMergeResult out; - - out.image_scale_g.resize(observations.size(), NAN); - out.mosaicity_deg.resize(observations.size(), NAN); - for (int i = 0; i < observations.size(); i++) { - size_t img_slot = i / opt.image_cluster; - if (image_slot_used[img_slot]) { - out.image_scale_g[i] = g[img_slot]; - out.mosaicity_deg[i] = mosaicity[img_slot]; - } - } - - std::vector slotToHKL(nhkl); - for (const auto &kv: hklToSlot) - slotToHKL[kv.second] = kv.first; - - out.merged.resize(nhkl); - for (int h = 0; h < nhkl; ++h) { - out.merged[h].h = slotToHKL[h].h; - out.merged[h].k = slotToHKL[h].k; - out.merged[h].l = slotToHKL[h].l; - out.merged[h].I = Itrue[h]; - out.merged[h].sigma = 0.0; - out.merged[h].d = 0.0; - } - - // Populate d from median of observations per HKL - { - std::vector > per_hkl_d(nhkl); - for (const auto &o: obs) { - const double d_val = static_cast(o.r->d); - if (std::isfinite(d_val) && d_val > 0.0) - per_hkl_d[o.hkl_slot].push_back(d_val); - } - for (int h = 0; h < nhkl; ++h) { - auto &v = per_hkl_d[h]; - if (!v.empty()) { - std::nth_element(v.begin(), v.begin() + static_cast(v.size() / 2), v.end()); - out.merged[h].d = v[v.size() / 2]; - } - } - } - - std::cout << summary.FullReport() << std::endl; - - // ---- Compute corrected observations once (used for both merging and statistics) ---- struct CorrectedObs { int hkl_slot; double I_corr; double sigma_corr; }; - std::vector corr_obs; - corr_obs.reserve(obs.size()); - { - const double half_wedge = opt.wedge_deg.value_or(0.0) / 2.0; + void scale(const ScaleMergeOptions &opt, + std::vector &g, + std::vector &mosaicity, + const std::vector &image_slot_used, + bool rotation_crystallography, + size_t nhkl, + const std::vector &obs) { + ceres::Problem problem; + + std::vector Itrue(nhkl, 0.0); + + // Initialize Itrue from per-HKL median of observed intensities + { + std::vector > per_hkl_I(nhkl); + for (const auto &o: obs) { + per_hkl_I[o.hkl_slot].push_back(static_cast(o.r->I)); + } + for (int h = 0; h < nhkl; ++h) { + auto &v = per_hkl_I[h]; + if (v.empty()) { + Itrue[h] = std::max(opt.min_sigma, 1e-6); + continue; + } + std::nth_element(v.begin(), v.begin() + static_cast(v.size() / 2), v.end()); + double med = v[v.size() / 2]; + if (!std::isfinite(med) || med <= opt.min_sigma) + med = opt.min_sigma; + Itrue[h] = med; + } + } + + std::vector is_valid_hkl_slot(nhkl, false); for (const auto &o: obs) { - const Reflection &r = *o.r; - const double lp = SafeInv(static_cast(r.rlp), 1.0); - const double G_i = g[o.img_id]; + auto *cost = new ceres::AutoDiffCostFunction( + new IntensityResidual(*o.r, o.sigma, opt.wedge_deg.value_or(0.0), rotation_crystallography)); + problem.AddResidualBlock(cost, + nullptr, + &g[o.img_id], + &mosaicity[o.img_id], + &Itrue[o.hkl_slot]); + is_valid_hkl_slot[o.hkl_slot] = true; + } - // Compute partiality with refined mosaicity - double partiality; - if (rotation_crystallography && mosaicity[o.img_id] > 0.0) { - const double c1 = r.zeta / std::sqrt(2.0); - const double arg_plus = (r.delta_phi_deg + half_wedge) * c1 / mosaicity[o.img_id]; - const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[o.img_id]; - partiality = (std::erf(arg_plus) - std::erf(arg_minus)) / 2.0; - } else { - partiality = r.partiality; + for (int i = 0; i < g.size(); ++i) { + if (image_slot_used[i]) { + auto *cost = new ceres::AutoDiffCostFunction( + new ScaleRegularizationResidual(0.05)); + problem.AddResidualBlock(cost, nullptr, &g[i]); + } + } + + if (rotation_crystallography) { + if (opt.smoothen_g) { + for (int i = 0; i < g.size() - 2; ++i) { + if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) { + auto *cost = new ceres::AutoDiffCostFunction( + new SmoothnessRegularizationResidual(0.05)); + + problem.AddResidualBlock(cost, nullptr, &g[i], &g[i + 1], &g[i + 2]); + } + } } - if (partiality <= opt.min_partiality_for_merge) - continue; - const double correction = G_i * partiality * lp; - if (correction <= 0.0) - continue; + if (opt.smoothen_mos && opt.refine_mosaicity) { + for (int i = 0; i < mosaicity.size() - 2; ++i) { + if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) { + auto *cost = new ceres::AutoDiffCostFunction( + new SmoothnessRegularizationResidual(0.05)); - corr_obs.push_back({ - o.hkl_slot, - static_cast(r.I) / correction, - o.sigma / correction - }); + problem.AddResidualBlock(cost, nullptr, &mosaicity[i], &mosaicity[i + 1], &mosaicity[i + 2]); + } + } + } } + + // Scaling factors must be always positive + for (int i = 0; i < g.size(); i++) { + if (image_slot_used[i]) + problem.SetParameterLowerBound(&g[i], 0, 1e-12); + } + + // Mosaicity refinement + bounds + if (!opt.refine_mosaicity) { + for (int i = 0; i < mosaicity.size(); ++i) { + if (image_slot_used[i]) + problem.SetParameterBlockConstant(&mosaicity[i]); + } + } else { + for (int i = 0; i < mosaicity.size(); ++i) { + if (image_slot_used[i]) { + problem.SetParameterLowerBound(&mosaicity[i], 0, opt.mosaicity_min_deg); + problem.SetParameterUpperBound(&mosaicity[i], 0, opt.mosaicity_max_deg); + } + } + } + + // use all available threads + unsigned int hw = std::thread::hardware_concurrency(); + if (hw == 0) + hw = 1; // fallback + + ceres::Solver::Options options; + + options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY; + options.minimizer_progress_to_stdout = true; + options.max_num_iterations = opt.max_num_iterations; + options.max_solver_time_in_seconds = opt.max_solver_time_s; + options.num_threads = static_cast(hw); + + ceres::Solver::Summary summary; + ceres::Solve(options, &problem, &summary); + std::cout << summary.FullReport() << std::endl; } - // ---- Merge (XDS/XSCALE style: inverse-variance weighted mean) ---- - { + void merge(size_t nhkl, ScaleMergeResult &out, const std::vector &corr_obs) { + // ---- Merge (XDS/XSCALE style: inverse-variance weighted mean) ---- + struct HKLAccum { double sum_wI = 0.0; double sum_w = 0.0; @@ -501,6 +352,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &corr_obs) // ---- Compute per-shell merging statistics ---- { constexpr int kStatShells = 10; @@ -657,6 +509,187 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &g, + std::vector &mosaicity, + bool rotation_crystallography, + const std::vector &obs, + std::vector &corr_obs) { + + // ---- Compute corrected observations once (used for both merging and statistics) ---- + const double half_wedge = opt.wedge_deg.value_or(0.0) / 2.0; + + for (const auto &o: obs) { + const Reflection &r = *o.r; + const double lp = SafeInv(static_cast(r.rlp), 1.0); + const double G_i = g[o.img_id]; + + // Compute partiality with refined mosaicity + double partiality; + if (rotation_crystallography && mosaicity[o.img_id] > 0.0) { + const double c1 = r.zeta / std::sqrt(2.0); + const double arg_plus = (r.delta_phi_deg + half_wedge) * c1 / mosaicity[o.img_id]; + const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[o.img_id]; + partiality = (std::erf(arg_plus) - std::erf(arg_minus)) / 2.0; + } else { + partiality = r.partiality; + } + + if (partiality <= opt.min_partiality_for_merge) + continue; + const double correction = G_i * partiality * lp; + if (correction <= 0.0) + continue; + + corr_obs.push_back({ + o.hkl_slot, + static_cast(r.I) / correction, + o.sigma / correction + }); + } + } + + void proc_obs(const std::vector > &observations, + const ScaleMergeOptions &opt, + std::vector &image_slot_used, + std::vector &obs, + std::unordered_map &hklToSlot + ) { + for (int i = 0; i < observations.size(); i++) { + for (const auto &r: observations[i]) { + const double d = SafeD(r.d); + if (!std::isfinite(d)) + continue; + if (!std::isfinite(r.I)) + continue; + + if (opt.d_min_limit_A > 0.0 && d < opt.d_min_limit_A) + continue; + + if (!std::isfinite(r.zeta) || r.zeta <= 0.0f) + continue; + if (!std::isfinite(r.rlp) || r.rlp == 0.0f) + continue; + + const double sigma = SafeSigma(r.sigma, opt.min_sigma); + + const int img_id = i / opt.image_cluster; + image_slot_used[img_id] = 1; + + int hkl_slot; + try { + const HKLKey key = CanonicalizeHKLKey(r, opt); + auto it = hklToSlot.find(key); + if (it == hklToSlot.end()) { + hkl_slot = static_cast(hklToSlot.size()); + hklToSlot.emplace(key, hkl_slot); + } else { + hkl_slot = it->second; + } + } catch (...) { + continue; + } + + ObsRef o; + o.r = &r; + o.img_id = img_id; + o.hkl_slot = hkl_slot; + o.sigma = sigma; + obs.push_back(o); + } + } + + } +} // namespace + +ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector > &observations, + const ScaleMergeOptions &opt) { + if (opt.image_cluster <= 0) + throw std::invalid_argument("image_cluster must be positive"); + + const bool rotation_crystallography = opt.wedge_deg.has_value(); + + size_t nrefl = 0; + for (const auto &i: observations) + nrefl += i.size(); + + std::vector obs; + std::vector corr_obs; + + obs.reserve(nrefl); + corr_obs.reserve(nrefl); + + std::unordered_map hklToSlot; + hklToSlot.reserve(nrefl); + + size_t n_image_slots = observations.size() / opt.image_cluster + + (observations.size() % opt.image_cluster > 0 ? 1 : 0); + + std::vector image_slot_used(n_image_slots, 0); + + proc_obs(observations, opt, image_slot_used, obs, hklToSlot); + + const int nhkl = static_cast(hklToSlot.size()); + + std::vector g(n_image_slots, 1.0); + std::vector mosaicity(n_image_slots, opt.mosaicity_init_deg); + for (int i = 0; i < n_image_slots; i++) { + if (!image_slot_used[i]) { + mosaicity[i] = NAN; + g[i] = NAN; + } else if (opt.mosaicity_init_deg_vec.size() > i && std::isfinite(opt.mosaicity_init_deg_vec[i])) { + mosaicity[i] = opt.mosaicity_init_deg_vec[i]; + } + } + + scale(opt, g, mosaicity, image_slot_used, rotation_crystallography, nhkl, obs); + + ScaleMergeResult out; + + out.image_scale_g.resize(observations.size(), NAN); + out.mosaicity_deg.resize(observations.size(), NAN); + for (int i = 0; i < observations.size(); i++) { + size_t img_slot = i / opt.image_cluster; + if (image_slot_used[img_slot]) { + out.image_scale_g[i] = g[img_slot]; + out.mosaicity_deg[i] = mosaicity[img_slot]; + } + } + + std::vector slotToHKL(nhkl); + for (const auto &kv: hklToSlot) + slotToHKL[kv.second] = kv.first; + + out.merged.resize(nhkl); + for (int h = 0; h < nhkl; ++h) { + out.merged[h].h = slotToHKL[h].h; + out.merged[h].k = slotToHKL[h].k; + out.merged[h].l = slotToHKL[h].l; + out.merged[h].I = 0.0; + out.merged[h].sigma = 0.0; + out.merged[h].d = 0.0; + } + + // Populate d from median of observations per HKL + { + std::vector > per_hkl_d(nhkl); + for (const auto &o: obs) { + const double d_val = static_cast(o.r->d); + if (std::isfinite(d_val) && d_val > 0.0) + per_hkl_d[o.hkl_slot].push_back(d_val); + } + for (int h = 0; h < nhkl; ++h) { + auto &v = per_hkl_d[h]; + if (!v.empty()) { + std::nth_element(v.begin(), v.begin() + static_cast(v.size() / 2), v.end()); + out.merged[h].d = v[v.size() / 2]; + } + } + } + + calc_obs(opt, g, mosaicity, rotation_crystallography, obs, corr_obs); + merge(nhkl, out, corr_obs); + stats(opt, nhkl, out, corr_obs); return out; } -- 2.49.1 From 5527fd0aa2e901e6d41dd2ed6ce899cf749a3a7a Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Fri, 20 Feb 2026 18:15:17 +0100 Subject: [PATCH 12/52] IndexAndRefine: For the moment - symmetry is enforced only if space group and unit cell are both provided --- image_analysis/IndexAndRefine.cpp | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 8605b232..db41ad43 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -65,17 +65,14 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data auto latt = indexer_result.lattice[0]; auto sg = experiment.GetGemmiSpaceGroup(); - // If space group provided => always enforce symmetry in refinement + // If space group and cell provided => always enforce symmetry in refinement // If space group not provided => guess symmetry - if (sg) { - // If space group provided but no unit cell fixed, it is better to keep triclinic for now - if (experiment.GetUnitCell()) { - outcome.symmetry = LatticeMessage{ - .centering = sg->centring_type(), - .niggli_class = 0, - .crystal_system = sg->crystal_system() - }; - } + if (sg && experiment.GetUnitCell()) { + outcome.symmetry = LatticeMessage{ + .centering = sg->centring_type(), + .niggli_class = 0, + .crystal_system = sg->crystal_system() + }; outcome.lattice_candidate = latt; } else { auto sym_result = LatticeSearch(latt); -- 2.49.1 From 7535ba56e3ee26a46105aeacbef4978ed3b9ad2a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 22 Feb 2026 09:38:00 +0100 Subject: [PATCH 13/52] ScaleAndMerge: Add more partiality models --- image_analysis/scale_merge/ScaleAndMerge.cpp | 79 +++++++++++++++----- image_analysis/scale_merge/ScaleAndMerge.h | 6 +- tools/jfjoch_process.cpp | 31 +++----- 3 files changed, 74 insertions(+), 42 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 90afbd75..e9102f04 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -115,8 +115,8 @@ namespace { return key; } - struct IntensityResidual { - IntensityResidual(const Reflection &r, double sigma_obs, double wedge_deg, bool refine_partiality) + struct IntensityRotResidual { + IntensityRotResidual(const Reflection &r, double sigma_obs, double wedge_deg, bool refine_partiality) : Iobs_(static_cast(r.I)), weight_(SafeInv(sigma_obs, 1.0)), delta_phi_(r.delta_phi_deg), @@ -133,12 +133,12 @@ namespace { const T *const Itrue, T *residual) const { T partiality; - if (refine_partiality_ && mosaicity[0] != 0.0) { + if (mosaicity[0] >= 0.0) { const T arg_plus = T(delta_phi_ + half_wedge_) * T(c1_) / mosaicity[0]; const T arg_minus = T(delta_phi_ - half_wedge_) * T(c1_) / mosaicity[0]; partiality = (ceres::erf(arg_plus) - ceres::erf(arg_minus)) / T(2.0); } else - partiality = T(partiality_); + partiality = T(1.0); const T Ipred = G[0] * partiality * T(lp_) * Itrue[0]; residual[0] = (Ipred - T(Iobs_)) * T(weight_); @@ -155,6 +155,27 @@ namespace { bool refine_partiality_; }; + struct IntensityFixedResidual { + IntensityFixedResidual(const Reflection &r, double sigma_obs, double partiality) + : Iobs_(static_cast(r.I)), + weight_(SafeInv(sigma_obs, 1.0)), + corr_(partiality * SafeInv(r.rlp, 1.0)) + {} + + template + bool operator()(const T *const G, + const T *const Itrue, + T *residual) const { + const T Ipred = T(corr_) * G[0] * Itrue[0]; + residual[0] = (Ipred - T(Iobs_)) * T(weight_); + return true; + } + + double Iobs_; + double weight_; + double corr_; + }; + struct ScaleRegularizationResidual { explicit ScaleRegularizationResidual(double sigma_k) : inv_sigma_(SafeInv(sigma_k, 1.0)) { @@ -231,15 +252,39 @@ namespace { } std::vector is_valid_hkl_slot(nhkl, false); - for (const auto &o: obs) { - auto *cost = new ceres::AutoDiffCostFunction( - new IntensityResidual(*o.r, o.sigma, opt.wedge_deg.value_or(0.0), rotation_crystallography)); - problem.AddResidualBlock(cost, - nullptr, - &g[o.img_id], - &mosaicity[o.img_id], - &Itrue[o.hkl_slot]); + switch (opt.partiality_model) { + case ScaleMergeOptions::PartialityModel::Rotation: { + auto *cost = new ceres::AutoDiffCostFunction( + new IntensityRotResidual(*o.r, o.sigma, opt.wedge_deg.value_or(0.0), rotation_crystallography)); + problem.AddResidualBlock(cost, + nullptr, + &g[o.img_id], + &mosaicity[o.img_id], + &Itrue[o.hkl_slot]); + } + break; + case ScaleMergeOptions::PartialityModel::Unity: { + auto *cost = new ceres::AutoDiffCostFunction( + new IntensityFixedResidual(*o.r, o.sigma, 1.0)); + problem.AddResidualBlock(cost, + nullptr, + &g[o.img_id], + &Itrue[o.hkl_slot]); + is_valid_hkl_slot[o.hkl_slot] = true; + } + break; + case ScaleMergeOptions::PartialityModel::Fixed: { + auto *cost = new ceres::AutoDiffCostFunction( + new IntensityFixedResidual(*o.r, o.sigma, o.r->partiality)); + problem.AddResidualBlock(cost, + nullptr, + &g[o.img_id], + &Itrue[o.hkl_slot]); + is_valid_hkl_slot[o.hkl_slot] = true; + } + break; + } is_valid_hkl_slot[o.hkl_slot] = true; } @@ -263,12 +308,11 @@ namespace { } } - if (opt.smoothen_mos && opt.refine_mosaicity) { + if (opt.smoothen_mos && opt.partiality_model == ScaleMergeOptions::PartialityModel::Rotation) { for (int i = 0; i < mosaicity.size() - 2; ++i) { if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) { auto *cost = new ceres::AutoDiffCostFunction( new SmoothnessRegularizationResidual(0.05)); - problem.AddResidualBlock(cost, nullptr, &mosaicity[i], &mosaicity[i + 1], &mosaicity[i + 2]); } } @@ -282,12 +326,7 @@ namespace { } // Mosaicity refinement + bounds - if (!opt.refine_mosaicity) { - for (int i = 0; i < mosaicity.size(); ++i) { - if (image_slot_used[i]) - problem.SetParameterBlockConstant(&mosaicity[i]); - } - } else { + if (opt.partiality_model == ScaleMergeOptions::PartialityModel::Rotation) { for (int i = 0; i < mosaicity.size(); ++i) { if (image_slot_used[i]) { problem.SetParameterLowerBound(&mosaicity[i], 0, opt.mosaicity_min_deg); diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h index 764ce01c..287188ed 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.h +++ b/image_analysis/scale_merge/ScaleAndMerge.h @@ -31,9 +31,7 @@ struct ScaleMergeOptions { std::optional wedge_deg; // --- Mosaicity (user input in degrees; internally converted to radians) --- - bool refine_mosaicity = true; - - double mosaicity_init_deg = 0.17; // ~0.003 rad + double mosaicity_init_deg = 0.17; double mosaicity_min_deg = 1e-3; double mosaicity_max_deg = 2.0; std::vector mosaicity_init_deg_vec; @@ -49,6 +47,8 @@ struct ScaleMergeOptions { double d_min_limit_A = 0.0; int64_t image_cluster = 1; + + enum class PartialityModel {Fixed, Rotation, Unity} partiality_model = PartialityModel::Fixed; }; struct MergedReflection { diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 95c56fbd..075b847f 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -42,7 +42,7 @@ void print_usage(Logger &logger) { logger.Info(" -D High resolution limit for scaling/merging (default: 0.0; no limit)"); logger.Info(" -S Space group number"); logger.Info(" -M Scale and merge (refine mosaicity) and write scaled.hkl + image.dat"); - logger.Info(" -m Mosaicity refinement none|image (default: image)"); + logger.Info(" -P Partiality refinement fixed|rot|unity (default: fixed)"); logger.Info(" -A Anomalous mode (don't merge Friedel pairs)"); } @@ -67,8 +67,7 @@ int main(int argc, char **argv) { bool anomalous_mode = false; std::optional space_group_number; - enum class MosaicityRefinementMode { None, Image }; - MosaicityRefinementMode mosaicity_refinement_mode = MosaicityRefinementMode::Image; + ScaleMergeOptions::PartialityModel partiality_model = ScaleMergeOptions::PartialityModel::Fixed; float d_min_spot_finding = 1.5; std::optional d_min_scale_merge; @@ -79,7 +78,7 @@ int main(int argc, char **argv) { } int opt; - while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:Mm:AD:")) != -1) { + while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:MP:AD:")) != -1) { switch (opt) { case 'o': output_prefix = optarg; @@ -122,13 +121,15 @@ int main(int argc, char **argv) { case 'A': anomalous_mode = true; break; - case 'm': - if (strcmp(optarg, "none") == 0) - mosaicity_refinement_mode = MosaicityRefinementMode::None; - else if (strcmp(optarg, "image") == 0) - mosaicity_refinement_mode = MosaicityRefinementMode::Image; + case 'P': + if (strcmp(optarg, "unity") == 0) + partiality_model = ScaleMergeOptions::PartialityModel::Unity; + else if (strcmp(optarg, "fixed") == 0) + partiality_model = ScaleMergeOptions::PartialityModel::Fixed; + else if (strcmp(optarg, "rot") == 0) + partiality_model = ScaleMergeOptions::PartialityModel::Rotation; else { - logger.Error("Invalid mosaicity refinement mode: {}", optarg); + logger.Error("Invalid partiality mode: {}", optarg); print_usage(logger); exit(EXIT_FAILURE); } @@ -421,20 +422,12 @@ int main(int argc, char **argv) { logger.Info("Running scaling (mosaicity refinement) ..."); ScaleMergeOptions scale_opts; - scale_opts.refine_mosaicity = true; + scale_opts.partiality_model = partiality_model; scale_opts.max_num_iterations = 500; scale_opts.max_solver_time_s = 240.0; // generous cutoff for now scale_opts.merge_friedel = !anomalous_mode; scale_opts.d_min_limit_A = d_min_scale_merge.value_or(0.0); - switch (mosaicity_refinement_mode) { - case MosaicityRefinementMode::None: - scale_opts.refine_mosaicity = false; - break; - case MosaicityRefinementMode::Image: - scale_opts.refine_mosaicity = true; - break; - } if (space_group) scale_opts.space_group = *space_group; else -- 2.49.1 From 3a33944dc447347eca875bf15e95d8e67d606462 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 22 Feb 2026 10:02:35 +0100 Subject: [PATCH 14/52] ScaleAndMerge: Set iteration limit to 10 --- image_analysis/scale_merge/ScaleAndMerge.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index e9102f04..7a4eba1b 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -347,6 +347,7 @@ namespace { options.max_num_iterations = opt.max_num_iterations; options.max_solver_time_in_seconds = opt.max_solver_time_s; options.num_threads = static_cast(hw); + options.max_num_iterations = 10; // 6 is usually enough to converge ceres::Solver::Summary summary; ceres::Solve(options, &problem, &summary); -- 2.49.1 From f8c326191fc1b40c6f813baca80f375140ee65a6 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Mon, 23 Feb 2026 14:43:06 +0100 Subject: [PATCH 15/52] ScaleAndMerge: Refine partiality model options --- image_analysis/scale_merge/ScaleAndMerge.cpp | 81 ++++++++++++++++---- image_analysis/scale_merge/ScaleAndMerge.h | 2 +- 2 files changed, 66 insertions(+), 17 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 7a4eba1b..e561db16 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -116,15 +116,14 @@ namespace { } struct IntensityRotResidual { - IntensityRotResidual(const Reflection &r, double sigma_obs, double wedge_deg, bool refine_partiality) + IntensityRotResidual(const Reflection &r, double sigma_obs, double wedge_deg) : Iobs_(static_cast(r.I)), weight_(SafeInv(sigma_obs, 1.0)), delta_phi_(r.delta_phi_deg), - lp_(SafeInv(static_cast(r.rlp), 1.0)), + lp_(SafeInv(r.rlp, 1.0)), half_wedge_(wedge_deg / 2.0), c1_(r.zeta / std::sqrt(2.0)), - partiality_(r.partiality), - refine_partiality_(refine_partiality) { + partiality_(r.partiality) { } template @@ -152,7 +151,31 @@ namespace { double half_wedge_; double c1_; double partiality_; - bool refine_partiality_; + }; + + struct IntensityStillResidual { + IntensityStillResidual(const Reflection &r, double sigma_obs) + : Iobs_(static_cast(r.I)), + weight_(SafeInv(sigma_obs, 1.0)), + lp_(SafeInv(r.rlp, 1.0)), + dist_ewald_sq_(r.dist_ewald) { + } + + template + bool operator()(const T *const G, + const T *const R, + const T *const Itrue, + T *residual) const { + const T partiality = ceres::exp(-T(dist_ewald_sq_)/R[0]); + const T Ipred = G[0] * partiality * T(lp_) * Itrue[0]; + residual[0] = (Ipred - T(Iobs_)) * T(weight_); + return true; + } + + double Iobs_; + double weight_; + double lp_; + double dist_ewald_sq_; }; struct IntensityFixedResidual { @@ -256,7 +279,7 @@ namespace { switch (opt.partiality_model) { case ScaleMergeOptions::PartialityModel::Rotation: { auto *cost = new ceres::AutoDiffCostFunction( - new IntensityRotResidual(*o.r, o.sigma, opt.wedge_deg.value_or(0.0), rotation_crystallography)); + new IntensityRotResidual(*o.r, o.sigma, opt.wedge_deg.value_or(0.0))); problem.AddResidualBlock(cost, nullptr, &g[o.img_id], @@ -264,6 +287,15 @@ namespace { &Itrue[o.hkl_slot]); } break; + case ScaleMergeOptions::PartialityModel::Still: { + auto *cost = new ceres::AutoDiffCostFunction( + new IntensityStillResidual(*o.r, o.sigma)); + problem.AddResidualBlock(cost, + nullptr, + &g[o.img_id], + &mosaicity[o.img_id], + &Itrue[o.hkl_slot]); + } case ScaleMergeOptions::PartialityModel::Unity: { auto *cost = new ceres::AutoDiffCostFunction( new IntensityFixedResidual(*o.r, o.sigma, 1.0)); @@ -271,7 +303,6 @@ namespace { nullptr, &g[o.img_id], &Itrue[o.hkl_slot]); - is_valid_hkl_slot[o.hkl_slot] = true; } break; case ScaleMergeOptions::PartialityModel::Fixed: { @@ -281,7 +312,6 @@ namespace { nullptr, &g[o.img_id], &Itrue[o.hkl_slot]); - is_valid_hkl_slot[o.hkl_slot] = true; } break; } @@ -319,6 +349,14 @@ namespace { } } + if (opt.partiality_model == ScaleMergeOptions::PartialityModel::Still) { + for (int i = 0; i < mosaicity.size(); ++i) { + if (image_slot_used[i]) { + problem.SetParameterLowerBound(&mosaicity[i], 0, 0.0); + } + } + } + // Scaling factors must be always positive for (int i = 0; i < g.size(); i++) { if (image_slot_used[i]) @@ -565,18 +603,29 @@ namespace { const double G_i = g[o.img_id]; // Compute partiality with refined mosaicity - double partiality; - if (rotation_crystallography && mosaicity[o.img_id] > 0.0) { - const double c1 = r.zeta / std::sqrt(2.0); - const double arg_plus = (r.delta_phi_deg + half_wedge) * c1 / mosaicity[o.img_id]; - const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[o.img_id]; - partiality = (std::erf(arg_plus) - std::erf(arg_minus)) / 2.0; - } else { - partiality = r.partiality; + double partiality = 1.0; + + switch (opt.partiality_model) { + case ScaleMergeOptions::PartialityModel::Fixed: + partiality = r.partiality; + break; + case ScaleMergeOptions::PartialityModel::Rotation: { + const double c1 = r.zeta / std::sqrt(2.0); + const double arg_plus = (r.delta_phi_deg + half_wedge) * c1 / mosaicity[o.img_id]; + const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[o.img_id]; + partiality = (std::erf(arg_plus) - std::erf(arg_minus)) / 2.0; + } + break; + case ScaleMergeOptions::PartialityModel::Still: + partiality = std::exp(-r.dist_ewald * r.dist_ewald / mosaicity[o.img_id]); + break; + case ScaleMergeOptions::PartialityModel::Unity: + break; } if (partiality <= opt.min_partiality_for_merge) continue; + const double correction = G_i * partiality * lp; if (correction <= 0.0) continue; diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h index 287188ed..9710ce05 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.h +++ b/image_analysis/scale_merge/ScaleAndMerge.h @@ -48,7 +48,7 @@ struct ScaleMergeOptions { int64_t image_cluster = 1; - enum class PartialityModel {Fixed, Rotation, Unity} partiality_model = PartialityModel::Fixed; + enum class PartialityModel {Fixed, Rotation, Unity, Still} partiality_model = PartialityModel::Fixed; }; struct MergedReflection { -- 2.49.1 From 94aea965eb6b6926e5913dd08c0398fdc290601a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 23 Feb 2026 14:53:19 +0100 Subject: [PATCH 16/52] jfjoch_process: Add still partiality model (though why it fails)? --- tools/jfjoch_process.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 075b847f..db751f45 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -128,6 +128,8 @@ int main(int argc, char **argv) { partiality_model = ScaleMergeOptions::PartialityModel::Fixed; else if (strcmp(optarg, "rot") == 0) partiality_model = ScaleMergeOptions::PartialityModel::Rotation; + else if (strcmp(optarg, "still") == 0) + partiality_model = ScaleMergeOptions::PartialityModel::Still; else { logger.Error("Invalid partiality mode: {}", optarg); print_usage(logger); -- 2.49.1 From 64390c71d176947103e6009ab2814f6bd1e4e654 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 23 Feb 2026 15:07:41 +0100 Subject: [PATCH 17/52] ScaleAndMerge: Don't allow zero in still mode (need to figure out proper bounds...) --- image_analysis/scale_merge/ScaleAndMerge.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index e561db16..f543e471 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -352,7 +352,7 @@ namespace { if (opt.partiality_model == ScaleMergeOptions::PartialityModel::Still) { for (int i = 0; i < mosaicity.size(); ++i) { if (image_slot_used[i]) { - problem.SetParameterLowerBound(&mosaicity[i], 0, 0.0); + problem.SetParameterLowerBound(&mosaicity[i], 0, 1e-6); } } } @@ -385,7 +385,6 @@ namespace { options.max_num_iterations = opt.max_num_iterations; options.max_solver_time_in_seconds = opt.max_solver_time_s; options.num_threads = static_cast(hw); - options.max_num_iterations = 10; // 6 is usually enough to converge ceres::Solver::Summary summary; ceres::Solve(options, &problem, &summary); -- 2.49.1 From b00ff5c9c1a6919c759f10c680b28ebfa5e91b89 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Mon, 23 Feb 2026 15:44:12 +0100 Subject: [PATCH 18/52] ScaleAndMerge: Change bounds for R --- image_analysis/scale_merge/ScaleAndMerge.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index f543e471..bc9b931b 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -352,7 +352,8 @@ namespace { if (opt.partiality_model == ScaleMergeOptions::PartialityModel::Still) { for (int i = 0; i < mosaicity.size(); ++i) { if (image_slot_used[i]) { - problem.SetParameterLowerBound(&mosaicity[i], 0, 1e-6); + problem.SetParameterLowerBound(&mosaicity[i], 0, 1e-9); + problem.SetParameterUpperBound(&mosaicity[i], 0, 1.0); } } } @@ -589,7 +590,6 @@ namespace { void calc_obs(const ScaleMergeOptions &opt, std::vector &g, std::vector &mosaicity, - bool rotation_crystallography, const std::vector &obs, std::vector &corr_obs) { @@ -775,7 +775,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector Date: Sat, 28 Feb 2026 13:05:56 +0100 Subject: [PATCH 19/52] StatusVector: Prereserve arrays for better performance --- common/StatusVector.cpp | 12 ++++-- common/StatusVector.h | 5 ++- receiver/JFJochReceiverPlots.cpp | 67 +++++++++++++++++--------------- 3 files changed, 48 insertions(+), 36 deletions(-) diff --git a/common/StatusVector.cpp b/common/StatusVector.cpp index b8765ca9..f5ddb3ab 100644 --- a/common/StatusVector.cpp +++ b/common/StatusVector.cpp @@ -4,9 +4,11 @@ #include "StatusVector.h" -void StatusVector::Clear() { +void StatusVector::Clear(size_t reserve) { std::unique_lock ul(m); content.clear(); + if (reserve > 0) + content.reserve(reserve); mean = NAN; count = 0; @@ -187,15 +189,19 @@ MultiLinePlot StatusVector::GetMaxPlot(int64_t bin_size, float x_start, float x_ return ret; } -void StatusMultiVector::Clear() { +void StatusMultiVector::Clear(size_t reserve) { std::unique_lock ul(m); + if (reserve > 0) + r = reserve; status.clear(); } void StatusMultiVector::AddElement(const std::string &s, uint32_t id, float val) { std::unique_lock ul(m); - if (!status.contains(s)) + if (!status.contains(s)) { status[s] = std::make_unique(); + status[s]->Clear(r); + } status[s]->AddElement(id, val); } diff --git a/common/StatusVector.h b/common/StatusVector.h index 2326223f..65ea4956 100644 --- a/common/StatusVector.h +++ b/common/StatusVector.h @@ -22,7 +22,7 @@ class StatusVector { size_t count = 0; double sum = 0; public: - void Clear(); + void Clear(size_t reserve = 0); void AddElement(uint32_t id, std::optional val); void AddElement(uint32_t id, float val); std::optional GetElement(uint32_t id) const; @@ -44,8 +44,9 @@ public: class StatusMultiVector { std::mutex m; std::map> status; + size_t r = 0; public: - void Clear(); + void Clear(size_t reserve = 0); void AddElement(const std::string& s, uint32_t id, float val); void AddElement(const std::string& s, uint32_t id, std::optional val); diff --git a/receiver/JFJochReceiverPlots.cpp b/receiver/JFJochReceiverPlots.cpp index aada4c8d..9581db1f 100644 --- a/receiver/JFJochReceiverPlots.cpp +++ b/receiver/JFJochReceiverPlots.cpp @@ -13,41 +13,46 @@ void JFJochReceiverPlots::Setup(const DiffractionExperiment &experiment, const A grid_scan = experiment.GetGridScan(); default_binning = experiment.GetDefaultPlotBinning(); + size_t r = experiment.GetImageNum(); // Reset all status vectors xfel_pulse_id.Clear(); + if (experiment.IsPulsedSource()) { + xfel_pulse_id.reserve(r); + xfel_event_code.reserve(r); + } xfel_event_code.Clear(); - bkg_estimate.Clear(); - spot_count.Clear(); - spot_count_low_res.Clear(); - spot_count_indexed.Clear(); - spot_count_ice.Clear(); + bkg_estimate.Clear(r); + spot_count.Clear(r); + spot_count_low_res.Clear(r); + spot_count_indexed.Clear(r); + spot_count_ice.Clear(r); - indexing_solution.Clear(); - indexing_unit_cell_angle.Clear(); - indexing_unit_cell_len.Clear(); - error_pixels.Clear(); - saturated_pixels.Clear(); - strong_pixels.Clear(); - receiver_delay.Clear(); - receiver_free_send_buf.Clear(); - image_collection_efficiency.Clear(); - roi_sum.Clear(); - roi_max_count.Clear(); - roi_pixels.Clear(); - roi_x.Clear(); - roi_y.Clear(); - roi_mean.Clear(); - packets_received.Clear(); - max_value.Clear(); - resolution_estimate.Clear(); - indexing_time.Clear(); - profile_radius.Clear(); - mosaicity_deg.Clear(); - b_factor.Clear(); - processing_time.Clear(); - beam_center_x.Clear(); - beam_center_y.Clear(); - pixel_sum.Clear(); + indexing_solution.Clear(r); + indexing_unit_cell_angle.Clear(r); + indexing_unit_cell_len.Clear(r); + error_pixels.Clear(r); + saturated_pixels.Clear(r); + strong_pixels.Clear(r); + receiver_delay.Clear(r); + receiver_free_send_buf.Clear(r); + image_collection_efficiency.Clear(r); + roi_sum.Clear(r); + roi_max_count.Clear(r); + roi_pixels.Clear(r); + roi_x.Clear(r); + roi_y.Clear(r); + roi_mean.Clear(r); + packets_received.Clear(r); + max_value.Clear(r); + resolution_estimate.Clear(r); + indexing_time.Clear(r); + profile_radius.Clear(r); + mosaicity_deg.Clear(r); + b_factor.Clear(r); + processing_time.Clear(r); + beam_center_x.Clear(r); + beam_center_y.Clear(r); + pixel_sum.Clear(r); } void JFJochReceiverPlots::Add(const DataMessage &msg, const AzimuthalIntegrationProfile &profile) { -- 2.49.1 From 9bb91ee84900ca981dce6ecac7abcf8f1fe0743f Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 13:06:23 +0100 Subject: [PATCH 20/52] ScaleAndMerge: Minor fixes --- image_analysis/scale_merge/ScaleAndMerge.cpp | 21 ++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index bc9b931b..74f339f7 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -158,7 +158,7 @@ namespace { : Iobs_(static_cast(r.I)), weight_(SafeInv(sigma_obs, 1.0)), lp_(SafeInv(r.rlp, 1.0)), - dist_ewald_sq_(r.dist_ewald) { + dist_ewald_sq_(r.dist_ewald * r.dist_ewald) { } template @@ -246,6 +246,7 @@ namespace { void scale(const ScaleMergeOptions &opt, std::vector &g, std::vector &mosaicity, + std::vector &R_sq, const std::vector &image_slot_used, bool rotation_crystallography, size_t nhkl, @@ -293,7 +294,7 @@ namespace { problem.AddResidualBlock(cost, nullptr, &g[o.img_id], - &mosaicity[o.img_id], + &R_sq[o.img_id], &Itrue[o.hkl_slot]); } case ScaleMergeOptions::PartialityModel::Unity: { @@ -350,10 +351,10 @@ namespace { } if (opt.partiality_model == ScaleMergeOptions::PartialityModel::Still) { - for (int i = 0; i < mosaicity.size(); ++i) { + for (int i = 0; i < R_sq.size(); ++i) { if (image_slot_used[i]) { - problem.SetParameterLowerBound(&mosaicity[i], 0, 1e-9); - problem.SetParameterUpperBound(&mosaicity[i], 0, 1.0); + problem.SetParameterLowerBound(&R_sq[i], 0, 1e-9); + problem.SetParameterUpperBound(&R_sq[i], 0, 1.0); } } } @@ -590,6 +591,7 @@ namespace { void calc_obs(const ScaleMergeOptions &opt, std::vector &g, std::vector &mosaicity, + std::vector &R_sq, const std::vector &obs, std::vector &corr_obs) { @@ -616,7 +618,7 @@ namespace { } break; case ScaleMergeOptions::PartialityModel::Still: - partiality = std::exp(-r.dist_ewald * r.dist_ewald / mosaicity[o.img_id]); + partiality = std::exp(-r.dist_ewald * r.dist_ewald / R_sq[o.img_id]); break; case ScaleMergeOptions::PartialityModel::Unity: break; @@ -721,16 +723,19 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector g(n_image_slots, 1.0); std::vector mosaicity(n_image_slots, opt.mosaicity_init_deg); + std::vector R_sq(n_image_slots, 0.001 * 0.001); + for (int i = 0; i < n_image_slots; i++) { if (!image_slot_used[i]) { mosaicity[i] = NAN; g[i] = NAN; + R_sq[i] = NAN; } else if (opt.mosaicity_init_deg_vec.size() > i && std::isfinite(opt.mosaicity_init_deg_vec[i])) { mosaicity[i] = opt.mosaicity_init_deg_vec[i]; } } - scale(opt, g, mosaicity, image_slot_used, rotation_crystallography, nhkl, obs); + scale(opt, g, mosaicity, R_sq, image_slot_used, rotation_crystallography, nhkl, obs); ScaleMergeResult out; @@ -775,7 +780,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector Date: Sat, 28 Feb 2026 15:40:44 +0100 Subject: [PATCH 21/52] JFJochReceiverPlots: Tiny bit safer code --- common/StatusVector.cpp | 1 + receiver/JFJochReceiverPlots.cpp | 11 ++++++++--- receiver/JFJochReceiverPlots.h | 3 ++- 3 files changed, 11 insertions(+), 4 deletions(-) diff --git a/common/StatusVector.cpp b/common/StatusVector.cpp index f5ddb3ab..34ec53b1 100644 --- a/common/StatusVector.cpp +++ b/common/StatusVector.cpp @@ -59,6 +59,7 @@ int32_t StatusVector::GetActualBinning(int32_t bin_size) const { } [[nodiscard]] float StatusVector::Mean() const { + std::unique_lock ul(m); return mean; } diff --git a/receiver/JFJochReceiverPlots.cpp b/receiver/JFJochReceiverPlots.cpp index 9581db1f..3edd7d0c 100644 --- a/receiver/JFJochReceiverPlots.cpp +++ b/receiver/JFJochReceiverPlots.cpp @@ -95,7 +95,8 @@ void JFJochReceiverPlots::Add(const DataMessage &msg, const AzimuthalIntegration { std::unique_lock ul(m); - *az_int_profile += profile; + if (az_int_profile) + *az_int_profile += profile; if (msg.xfel_pulse_id.has_value()) xfel_pulse_id[msg.number] = msg.xfel_pulse_id.value(); if (msg.xfel_event_code.has_value()) @@ -288,10 +289,12 @@ std::optional JFJochReceiverPlots::GetBkgEstimate() const { } void JFJochReceiverPlots::GetXFELPulseID(std::vector &v) const { + std::unique_lock ul(m); v = xfel_pulse_id.vec(); } void JFJochReceiverPlots::GetXFELEventCode(std::vector &v) const { + std::unique_lock ul(m); v = xfel_event_code.vec(); } @@ -383,12 +386,14 @@ void JFJochReceiverPlots::GetPlotRaw(std::vector &v, PlotType type, const break; case PlotType::AzInt: { std::unique_lock ul(m); - v = az_int_profile->GetResult(); + if (az_int_profile) + v = az_int_profile->GetResult(); break; } case PlotType::AzInt1D: { std::unique_lock ul(m); - v = az_int_profile->GetResult1D(); + if (az_int_profile) + v = az_int_profile->GetResult1D(); break; } case PlotType::PacketsReceived: diff --git a/receiver/JFJochReceiverPlots.h b/receiver/JFJochReceiverPlots.h index 70e2bcd9..7c81a87b 100644 --- a/receiver/JFJochReceiverPlots.h +++ b/receiver/JFJochReceiverPlots.h @@ -15,7 +15,8 @@ #include "../common/ScanResult.h" class JFJochReceiverPlots { - mutable std::mutex m; + mutable std::mutex m; // protects xfel_pulse_id, xfel_event_code and az_int_profile + std::optional goniometer; std::optional grid_scan; int64_t default_binning = 1; -- 2.49.1 From 6e97b7c15c65616bbbab29aaaaebafb0975176d2 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 16:09:14 +0100 Subject: [PATCH 22/52] JFJochReceiverPlots: Safer code --- common/StatusVector.cpp | 53 --------- common/StatusVector.h | 16 --- receiver/JFJochReceiverPlots.cpp | 190 ++++++++++++++++++++++--------- receiver/JFJochReceiverPlots.h | 34 ++++-- tests/StatusVectorTest.cpp | 22 ---- 5 files changed, 164 insertions(+), 151 deletions(-) diff --git a/common/StatusVector.cpp b/common/StatusVector.cpp index 34ec53b1..a0af404a 100644 --- a/common/StatusVector.cpp +++ b/common/StatusVector.cpp @@ -189,56 +189,3 @@ MultiLinePlot StatusVector::GetMaxPlot(int64_t bin_size, float x_start, float x_ ret.AddPlot(GetMaxPerBin(bin_size, x_start, x_incr, fill_value)); return ret; } - -void StatusMultiVector::Clear(size_t reserve) { - std::unique_lock ul(m); - if (reserve > 0) - r = reserve; - status.clear(); -} - -void StatusMultiVector::AddElement(const std::string &s, uint32_t id, float val) { - std::unique_lock ul(m); - if (!status.contains(s)) { - status[s] = std::make_unique(); - status[s]->Clear(r); - } - status[s]->AddElement(id, val); -} - -void StatusMultiVector::AddElement(const std::string &s, uint32_t id, std::optional val) { - // no need to lock, as AddElement(string, u32, T) has lock already - if (val.has_value()) - AddElement(s, id, val.value()); -} - -MultiLinePlotStruct StatusMultiVector::GetMeanPerBin(const std::string &in_key, int64_t bin_size, float x_start, - float x_incr, - const std::optional &fill_value) const { - MultiLinePlotStruct ret{}; - for (const auto &[key, value]: status) { - if (key == in_key) { - ret = value->GetMeanPerBin(bin_size, x_start, x_incr, fill_value); - ret.title = key; - } - } - return ret; -} - -MultiLinePlot StatusMultiVector::GetMeanPlot(int64_t bin_size, float x_start, float x_incr, - const std::optional &fill_value) const { - MultiLinePlot ret; - for (const auto &[key, value]: status) { - auto tmp = value->GetMeanPerBin(bin_size, x_start, x_incr, fill_value); - tmp.title = key; - ret.AddPlot(tmp); - } - return ret; -} - -std::vector StatusMultiVector::ExportArray(const std::string &s) const { - auto iter = status.find(s); - if (iter == status.end() || !iter->second) - return {}; - return iter->second->ExportArray(); -} diff --git a/common/StatusVector.h b/common/StatusVector.h index 65ea4956..c8c7e4b7 100644 --- a/common/StatusVector.h +++ b/common/StatusVector.h @@ -41,20 +41,4 @@ public: const std::optional &fill_value = {}) const; }; -class StatusMultiVector { - std::mutex m; - std::map> status; - size_t r = 0; -public: - void Clear(size_t reserve = 0); - void AddElement(const std::string& s, uint32_t id, float val); - void AddElement(const std::string& s, uint32_t id, std::optional val); - - [[nodiscard]] MultiLinePlotStruct GetMeanPerBin(const std::string& in_key, int64_t bin_size, float x_start, float x_incr, - const std::optional &fill_value = {}) const; - [[nodiscard]] MultiLinePlot GetMeanPlot(int64_t bin_size, float x_start, float x_incr, - const std::optional &fill_value = {}) const; - [[nodiscard]] std::vector ExportArray(const std::string& s) const; -}; - #endif //JUNGFRAUJOCH_STATUSVECTOR_H diff --git a/receiver/JFJochReceiverPlots.cpp b/receiver/JFJochReceiverPlots.cpp index 3edd7d0c..5578770d 100644 --- a/receiver/JFJochReceiverPlots.cpp +++ b/receiver/JFJochReceiverPlots.cpp @@ -3,6 +3,44 @@ #include "JFJochReceiverPlots.h" +#include + +MultiLinePlot JFJochReceiverPlots::GetROIPlot(PlotType type, int64_t nbins, float start, float incr, + const std::optional &fill_value) const { + MultiLinePlot ret; + + std::shared_lock sl(roi_m); + for (const auto &[key, roi] : roi_status) { + MultiLinePlotStruct plot; + switch (type) { + case PlotType::ROISum: + plot = roi.sum.GetMeanPerBin(nbins, start, incr, fill_value); + break; + case PlotType::ROIMaxCount: + plot = roi.max_count.GetMeanPerBin(nbins, start, incr, fill_value); + break; + case PlotType::ROIPixels: + plot = roi.pixels.GetMeanPerBin(nbins, start, incr, fill_value); + break; + case PlotType::ROIMean: + plot = roi.mean.GetMeanPerBin(nbins, start, incr, fill_value); + break; + case PlotType::ROIWeightedX: + plot = roi.x.GetMeanPerBin(nbins, start, incr, fill_value); + break; + case PlotType::ROIWeightedY: + plot = roi.y.GetMeanPerBin(nbins, start, incr, fill_value); + break; + default: + continue; + } + plot.title = key; + ret.AddPlot(plot); + } + + return ret; +} + void JFJochReceiverPlots::Setup(const DiffractionExperiment &experiment, const AzimuthalIntegration &mapping) { std::unique_lock ul(m); @@ -14,13 +52,14 @@ void JFJochReceiverPlots::Setup(const DiffractionExperiment &experiment, const A default_binning = experiment.GetDefaultPlotBinning(); size_t r = experiment.GetImageNum(); + // Reset all status vectors xfel_pulse_id.Clear(); + xfel_event_code.Clear(); if (experiment.IsPulsedSource()) { xfel_pulse_id.reserve(r); xfel_event_code.reserve(r); } - xfel_event_code.Clear(); bkg_estimate.Clear(r); spot_count.Clear(r); spot_count_low_res.Clear(r); @@ -28,20 +67,33 @@ void JFJochReceiverPlots::Setup(const DiffractionExperiment &experiment, const A spot_count_ice.Clear(r); indexing_solution.Clear(r); - indexing_unit_cell_angle.Clear(r); - indexing_unit_cell_len.Clear(r); + indexing_uc_a.Clear(r); + indexing_uc_b.Clear(r); + indexing_uc_c.Clear(r); + indexing_uc_alpha.Clear(r); + indexing_uc_beta.Clear(r); + indexing_uc_gamma.Clear(r); error_pixels.Clear(r); saturated_pixels.Clear(r); strong_pixels.Clear(r); receiver_delay.Clear(r); receiver_free_send_buf.Clear(r); image_collection_efficiency.Clear(r); - roi_sum.Clear(r); - roi_max_count.Clear(r); - roi_pixels.Clear(r); - roi_x.Clear(r); - roi_y.Clear(r); - roi_mean.Clear(r); + + { + std::unique_lock roi_lock(roi_m); + roi_status.clear(); + for (const auto &[name, _id] : experiment.ROI().GetROINameMap()) { + auto &entry = roi_status[name]; + entry.sum.Clear(r); + entry.max_count.Clear(r); + entry.pixels.Clear(r); + entry.x.Clear(r); + entry.y.Clear(r); + entry.mean.Clear(r); + } + } + packets_received.Clear(r); max_value.Clear(r); resolution_estimate.Clear(r); @@ -77,12 +129,12 @@ void JFJochReceiverPlots::Add(const DataMessage &msg, const AzimuthalIntegration processing_time.AddElement(msg.number, msg.processing_time_s); if (msg.indexing_unit_cell) { - indexing_unit_cell_len.AddElement("a", msg.number, msg.indexing_unit_cell->a); - indexing_unit_cell_len.AddElement("b", msg.number, msg.indexing_unit_cell->b); - indexing_unit_cell_len.AddElement("c", msg.number, msg.indexing_unit_cell->c); - indexing_unit_cell_angle.AddElement("alpha", msg.number, msg.indexing_unit_cell->alpha); - indexing_unit_cell_angle.AddElement("beta", msg.number, msg.indexing_unit_cell->beta); - indexing_unit_cell_angle.AddElement("gamma", msg.number, msg.indexing_unit_cell->gamma); + indexing_uc_a.AddElement(msg.number, msg.indexing_unit_cell->a); + indexing_uc_b.AddElement(msg.number, msg.indexing_unit_cell->b); + indexing_uc_c.AddElement(msg.number, msg.indexing_unit_cell->c); + indexing_uc_alpha.AddElement(msg.number, msg.indexing_unit_cell->alpha); + indexing_uc_beta.AddElement(msg.number, msg.indexing_unit_cell->beta); + indexing_uc_gamma.AddElement(msg.number, msg.indexing_unit_cell->gamma); } beam_center_x.AddElement(msg.number, msg.beam_corr_x); @@ -104,12 +156,22 @@ void JFJochReceiverPlots::Add(const DataMessage &msg, const AzimuthalIntegration } for (const auto &[key, value] : msg.roi) { - roi_sum.AddElement(key, msg.number, value.sum); - roi_mean.AddElement(key, msg.number, static_cast(value.sum) / static_cast(value.pixels)); - roi_max_count.AddElement(key, msg.number, value.max_count); - roi_pixels.AddElement(key, msg.number, value.pixels); - roi_x.AddElement(key, msg.number, static_cast(value.x_weighted) / static_cast(value.sum)); - roi_y.AddElement(key, msg.number, static_cast(value.y_weighted) / static_cast(value.sum)); + if (value.pixels == 0) + continue; + + std::shared_lock sl(roi_m); + auto it = roi_status.find(key); + if (it == roi_status.end()) + continue; // ROI not configured in setup -> ignore + + it->second.sum.AddElement(msg.number, value.sum); + it->second.mean.AddElement(msg.number, static_cast(value.sum) / static_cast(value.pixels)); + it->second.max_count.AddElement(msg.number, value.max_count); + it->second.pixels.AddElement(msg.number, value.pixels); + if (value.sum > 0) { + it->second.x.AddElement(msg.number, static_cast(value.x_weighted) / static_cast(value.sum)); + it->second.y.AddElement(msg.number, static_cast(value.y_weighted) / static_cast(value.sum)); + } } } @@ -120,7 +182,7 @@ void JFJochReceiverPlots::AddEmptyImage(const DataMessage &msg) { MultiLinePlot JFJochReceiverPlots::GetPlots(const PlotRequest &request) { MultiLinePlot ret; MultiLinePlotUnits units = MultiLinePlotUnits::ImageNumber; - int64_t nbins; + int64_t nbins = 1; std::optional local_grid_scan; float start = 0.0; @@ -137,6 +199,7 @@ MultiLinePlot JFJochReceiverPlots::GetPlots(const PlotRequest &request) { nbins = default_binning; if (request.binning > 0) nbins = request.binning; + nbins = std::max(1, nbins); if (request.experimental_coord && goniometer) { start = goniometer->GetStart_deg(); @@ -211,22 +274,12 @@ MultiLinePlot JFJochReceiverPlots::GetPlots(const PlotRequest &request) { ret = strong_pixels.GetMeanPlot(nbins, start, incr, request.fill_value); break; case PlotType::ROISum: - ret = roi_sum.GetMeanPlot(nbins, start, incr, request.fill_value); - break; case PlotType::ROIMaxCount: - ret = roi_max_count.GetMeanPlot(nbins, start, incr, request.fill_value); - break; case PlotType::ROIPixels: - ret = roi_pixels.GetMeanPlot(nbins, start, incr, request.fill_value); - break; case PlotType::ROIMean: - ret = roi_mean.GetMeanPlot(nbins, start, incr, request.fill_value); - break; case PlotType::ROIWeightedX: - ret = roi_x.GetMeanPlot(nbins, start, incr, request.fill_value); - break; case PlotType::ROIWeightedY: - ret = roi_y.GetMeanPlot(nbins, start, incr, request.fill_value); + ret = GetROIPlot(request.type, nbins, start, incr, request.fill_value); break; case PlotType::AzInt: ret = GetAzIntProfilePlot(false, request.azint_unit); @@ -234,18 +287,36 @@ MultiLinePlot JFJochReceiverPlots::GetPlots(const PlotRequest &request) { case PlotType::AzInt1D: ret = GetAzIntProfilePlot(true, request.azint_unit); break; - case PlotType::IndexingUnitCellLength: - ret = indexing_unit_cell_len.GetMeanPlot(nbins, start, incr, request.fill_value); + case PlotType::IndexingUnitCellLength: { + auto a = indexing_uc_a.GetMeanPerBin(nbins, start, incr, request.fill_value); + auto b = indexing_uc_b.GetMeanPerBin(nbins, start, incr, request.fill_value); + auto c = indexing_uc_c.GetMeanPerBin(nbins, start, incr, request.fill_value); + a.title = "a"; + b.title = "b"; + c.title = "c"; + ret.AddPlot(a); + ret.AddPlot(b); + ret.AddPlot(c); break; - case PlotType::IndexingUnitCellAngle: - ret = indexing_unit_cell_angle.GetMeanPlot(nbins, start, incr, request.fill_value); + } + case PlotType::IndexingUnitCellAngle: { + auto alpha = indexing_uc_alpha.GetMeanPerBin(nbins, start, incr, request.fill_value); + auto beta = indexing_uc_beta.GetMeanPerBin(nbins, start, incr, request.fill_value); + auto gamma = indexing_uc_gamma.GetMeanPerBin(nbins, start, incr, request.fill_value); + alpha.title = "alpha"; + beta.title = "beta"; + gamma.title = "gamma"; + ret.AddPlot(alpha); + ret.AddPlot(beta); + ret.AddPlot(gamma); break; + } case PlotType::PacketsReceived: ret = packets_received.GetMeanPlot(nbins, start, incr, request.fill_value); break; case PlotType::MaxValue: ret = max_value.GetMaxPlot(nbins, start, incr, request.fill_value); - break; // doesn't make sense to give mean here + break; case PlotType::PixelSum: ret = pixel_sum.GetMeanPlot(nbins, start, incr, request.fill_value); break; @@ -259,7 +330,6 @@ MultiLinePlot JFJochReceiverPlots::GetPlots(const PlotRequest &request) { ret = beam_center_y.GetMeanPlot(nbins, start, incr, request.fill_value); break; default: - // Do nothing break; } @@ -367,23 +437,41 @@ void JFJochReceiverPlots::GetPlotRaw(std::vector &v, PlotType type, const v = strong_pixels.ExportArray(); break; case PlotType::ROISum: - v = roi_sum.ExportArray(roi); - break; case PlotType::ROIMaxCount: - v = roi_max_count.ExportArray(roi); - break; case PlotType::ROIPixels: - v = roi_pixels.ExportArray(roi); - break; case PlotType::ROIMean: - v = roi_mean.ExportArray(roi); - break; case PlotType::ROIWeightedX: - v = roi_x.ExportArray(roi); - break; - case PlotType::ROIWeightedY: - v = roi_y.ExportArray(roi); + case PlotType::ROIWeightedY: { + std::shared_lock sl(roi_m); + auto it = roi_status.find(roi); + if (it == roi_status.end()) { + v.clear(); + break; + } + switch (type) { + case PlotType::ROISum: + v = it->second.sum.ExportArray(); + break; + case PlotType::ROIMaxCount: + v = it->second.max_count.ExportArray(); + break; + case PlotType::ROIPixels: + v = it->second.pixels.ExportArray(); + break; + case PlotType::ROIMean: + v = it->second.mean.ExportArray(); + break; + case PlotType::ROIWeightedX: + v = it->second.x.ExportArray(); + break; + case PlotType::ROIWeightedY: + v = it->second.y.ExportArray(); + break; + default: + break; + } break; + } case PlotType::AzInt: { std::unique_lock ul(m); if (az_int_profile) diff --git a/receiver/JFJochReceiverPlots.h b/receiver/JFJochReceiverPlots.h index 7c81a87b..d7989467 100644 --- a/receiver/JFJochReceiverPlots.h +++ b/receiver/JFJochReceiverPlots.h @@ -4,6 +4,10 @@ #ifndef JUNGFRAUJOCH_JFJOCHRECEIVERPLOTS_H #define JUNGFRAUJOCH_JFJOCHRECEIVERPLOTS_H +#include +#include +#include + #include "../common/StatusVector.h" #include "../common/Histogram.h" #include "../common/ADUHistogram.h" @@ -31,8 +35,12 @@ class JFJochReceiverPlots { StatusVector spot_count_ice; StatusVector indexing_solution; - StatusMultiVector indexing_unit_cell_len; - StatusMultiVector indexing_unit_cell_angle; + StatusVector indexing_uc_a; + StatusVector indexing_uc_b; + StatusVector indexing_uc_c; + StatusVector indexing_uc_alpha; + StatusVector indexing_uc_beta; + StatusVector indexing_uc_gamma; StatusVector error_pixels; StatusVector saturated_pixels; StatusVector strong_pixels; @@ -42,12 +50,18 @@ class JFJochReceiverPlots { StatusVector packets_received; StatusVector max_value; StatusVector resolution_estimate; - StatusMultiVector roi_sum; - StatusMultiVector roi_max_count; - StatusMultiVector roi_pixels; - StatusMultiVector roi_x; - StatusMultiVector roi_y; - StatusMultiVector roi_mean; + + struct ROIStatus { + StatusVector sum; + StatusVector max_count; + StatusVector pixels; + StatusVector x; + StatusVector y; + StatusVector mean; + }; + mutable std::shared_mutex roi_m; + std::map roi_status; + StatusVector indexing_time; StatusVector profile_radius; StatusVector mosaicity_deg; @@ -58,6 +72,9 @@ class JFJochReceiverPlots { StatusVector beam_center_y; StatusVector processing_time; + + MultiLinePlot GetROIPlot(PlotType type, int64_t nbins, float start, float incr, + const std::optional &fill_value) const; public: void Setup(const DiffractionExperiment& experiment, const AzimuthalIntegration& mapping); @@ -77,5 +94,4 @@ public: void GetPlotRaw(std::vector &v, PlotType type, const std::string &roi); }; - #endif //JUNGFRAUJOCH_JFJOCHRECEIVERPLOTS_H diff --git a/tests/StatusVectorTest.cpp b/tests/StatusVectorTest.cpp index f919892e..6562ddc6 100644 --- a/tests/StatusVectorTest.cpp +++ b/tests/StatusVectorTest.cpp @@ -272,28 +272,6 @@ TEST_CASE("StatusVector_Plot_NoBinning","[StatusVector]") { REQUIRE(plot_out.GetPlots()[0].y[0] == Catch::Approx(11)); } -TEST_CASE("StatusMultiVector","[StatusMultiVector]") { - StatusMultiVector status_vector; - status_vector.AddElement("plot1", 0, 4); - status_vector.AddElement("plot1", 1, 3); - status_vector.AddElement("plot2", 0, 5); - status_vector.AddElement("plot2", 1, 4); - - auto ret = status_vector.GetMeanPlot(1, 3.0, 5.0); - REQUIRE(ret.GetPlots().size() == 2); - REQUIRE(ret.GetPlots()[0].title == "plot1"); - REQUIRE(ret.GetPlots()[0].x.size() == 2); - REQUIRE(ret.GetPlots()[0].y.size() == 2); - REQUIRE(ret.GetPlots()[0].x[0] == 3.0 + 5.0 * 0.0f); - REQUIRE(ret.GetPlots()[0].y[1] == 3.0f); - - REQUIRE(ret.GetPlots()[1].title == "plot2"); - REQUIRE(ret.GetPlots()[1].x.size() == 2); - REQUIRE(ret.GetPlots()[1].y.size() == 2); - REQUIRE(ret.GetPlots()[1].x[1] == 3.0 + 5.0 * 1.0f); - REQUIRE(ret.GetPlots()[1].y[1] == 4.0f); -} - TEST_CASE("StatusVector_Clear","[StatusVector]") { StatusVector status_vector; status_vector.AddElement(5, 800); -- 2.49.1 From 79cb6d307327d9396798e0c6e18a4de0cdbeb219 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 16:18:05 +0100 Subject: [PATCH 23/52] ImageBuffer: Minor fixes --- common/ImageBuffer.cpp | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/common/ImageBuffer.cpp b/common/ImageBuffer.cpp index c88a9795..06c29a00 100644 --- a/common/ImageBuffer.cpp +++ b/common/ImageBuffer.cpp @@ -15,15 +15,18 @@ ImageBuffer::ImageBuffer(size_t buffer_size_bytes) : buffer_size(buffer_size_bytes) { #ifdef JFJOCH_USE_NUMA - buffer = (uint8_t *) numa_alloc_interleaved(buffer_size); -#else - buffer = (uint8_t *) mmap (nullptr, buffer_size, PROT_READ | PROT_WRITE, - MAP_PRIVATE|MAP_ANONYMOUS, -1, 0) ; -#endif + buffer = static_cast(numa_alloc_interleaved(buffer_size)); if (buffer == nullptr) throw JFJochException(JFJochExceptionCategory::MemAllocFailed, "Failed to allocate image buffer"); + #else + buffer = (uint8_t *) mmap (nullptr, buffer_size, PROT_READ | PROT_WRITE, + MAP_PRIVATE|MAP_ANONYMOUS, -1, 0) ; + if (buffer == MAP_FAILED) + throw JFJochException(JFJochExceptionCategory::MemAllocFailed, + "Failed to allocate image buffer"); +#endif memset(buffer, 0, buffer_size); } @@ -241,6 +244,7 @@ void ImageBuffer::GetStartMessage(std::vector &out_v) { } void ImageBuffer::SaveStartMessage(const std::vector &msg) { + std::unique_lock ul(start_message_mutex); start_message = msg; } -- 2.49.1 From 88020744cf5eaf6db9e4f8c57068cc348d2a2941 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 16:18:25 +0100 Subject: [PATCH 24/52] JFJochReceiverCurrentStatus: Use std::atomic for progress --- receiver/JFJochReceiverCurrentStatus.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/receiver/JFJochReceiverCurrentStatus.h b/receiver/JFJochReceiverCurrentStatus.h index 28dda27e..d7e4789e 100644 --- a/receiver/JFJochReceiverCurrentStatus.h +++ b/receiver/JFJochReceiverCurrentStatus.h @@ -4,6 +4,7 @@ #ifndef JFJOCH_JFJOCHRECEIVERCURRENTSTATUS_H #define JFJOCH_JFJOCHRECEIVERCURRENTSTATUS_H +#include #include #include #include @@ -29,7 +30,7 @@ struct JFJochReceiverStatus { class JFJochReceiverCurrentStatus { mutable std::mutex m; std::optional status; - volatile float progress = NAN; + std::atomic progress = NAN; public: std::optional GetStatus() const; void Clear(); -- 2.49.1 From 6e2c65011e88448346e9c77b4540b890b8db9763 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 16:29:18 +0100 Subject: [PATCH 25/52] Remove volatile, use std::atomic instead --- acquisition_device/FPGAAcquisitionDevice.h | 2 +- broker/JFJochServices.h | 2 +- broker/JFJochStateMachine.h | 4 ++-- image_puller/ZMQImagePuller.h | 2 +- writer/StreamWriter.h | 2 +- writer/jfjoch_writer.cpp | 2 +- 6 files changed, 7 insertions(+), 7 deletions(-) diff --git a/acquisition_device/FPGAAcquisitionDevice.h b/acquisition_device/FPGAAcquisitionDevice.h index 96234e51..f81ed235 100644 --- a/acquisition_device/FPGAAcquisitionDevice.h +++ b/acquisition_device/FPGAAcquisitionDevice.h @@ -26,7 +26,7 @@ class FPGAAcquisitionDevice : public AcquisitionDevice { void ReadWorkCompletionThread(); std::future send_work_request_future; - volatile bool stop_work_requests = false; + std::atomic stop_work_requests = false; void SendWorkRequestThread(); virtual void HW_LoadCalibration(const LoadCalibrationConfig &config) = 0; diff --git a/broker/JFJochServices.h b/broker/JFJochServices.h index 6d18cad2..8e697184 100644 --- a/broker/JFJochServices.h +++ b/broker/JFJochServices.h @@ -18,7 +18,7 @@ class JFJochServices { JFJochReceiverService *receiver = nullptr; std::unique_ptr detector; - volatile bool cannot_stop_detector = false; + std::atomic cannot_stop_detector = false; Logger &logger; public: explicit JFJochServices(Logger &in_logger); diff --git a/broker/JFJochStateMachine.h b/broker/JFJochStateMachine.h index 661bb733..5bd7852c 100644 --- a/broker/JFJochStateMachine.h +++ b/broker/JFJochStateMachine.h @@ -101,8 +101,8 @@ class JFJochStateMachine { // mutex m is protecting: mutable std::mutex m; std::condition_variable c; - volatile JFJochState state = JFJochState::Inactive; // state should not be set directly, but through SetState function - volatile bool cancel_sequence = false; + std::atomic state = JFJochState::Inactive; // state should not be set directly, but through SetState function + std::atomic cancel_sequence = false; std::unique_ptr calibration; PixelMask pixel_mask; int64_t current_detector_setup; // Lock only on change diff --git a/image_puller/ZMQImagePuller.h b/image_puller/ZMQImagePuller.h index 8a70a360..0a10c307 100644 --- a/image_puller/ZMQImagePuller.h +++ b/image_puller/ZMQImagePuller.h @@ -26,7 +26,7 @@ class ZMQImagePuller : public ImagePuller { ZMQSocket socket; std::string addr; - volatile int disconnect = 0; + std::atomic disconnect = 0; std::unique_ptr repub_socket; diff --git a/writer/StreamWriter.h b/writer/StreamWriter.h index 5bebc3d3..045d0d50 100644 --- a/writer/StreamWriter.h +++ b/writer/StreamWriter.h @@ -31,7 +31,7 @@ struct StreamWriterOutput { }; class StreamWriter { - volatile bool abort = false; + std::atomic abort = false; const bool verbose; std::string file_done_address; diff --git a/writer/jfjoch_writer.cpp b/writer/jfjoch_writer.cpp index ca8d9239..5cf54152 100644 --- a/writer/jfjoch_writer.cpp +++ b/writer/jfjoch_writer.cpp @@ -12,7 +12,7 @@ static Logger logger("jfjoch_writer"); static StreamWriter *writer; -volatile static bool quitok = false; +static std::atomic quitok = false; void print_usage() { logger.Info("Usage ./jfjoch_writer {options}
"); -- 2.49.1 From 0ba6f495d1089367ca3d27b42524fec7c5cb52a5 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 16:47:10 +0100 Subject: [PATCH 26/52] ScaleAndMerge: Fix fall through in switch --- image_analysis/scale_merge/ScaleAndMerge.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 74f339f7..14b129cb 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -288,15 +288,16 @@ namespace { &Itrue[o.hkl_slot]); } break; - case ScaleMergeOptions::PartialityModel::Still: { - auto *cost = new ceres::AutoDiffCostFunction( + case ScaleMergeOptions::PartialityModel::Still: { + auto *cost = new ceres::AutoDiffCostFunction( new IntensityStillResidual(*o.r, o.sigma)); - problem.AddResidualBlock(cost, - nullptr, - &g[o.img_id], - &R_sq[o.img_id], - &Itrue[o.hkl_slot]); - } + problem.AddResidualBlock(cost, + nullptr, + &g[o.img_id], + &R_sq[o.img_id], + &Itrue[o.hkl_slot]); + } + break; case ScaleMergeOptions::PartialityModel::Unity: { auto *cost = new ceres::AutoDiffCostFunction( new IntensityFixedResidual(*o.r, o.sigma, 1.0)); -- 2.49.1 From 64067f743067afa8ee21a6d6f6807579e03ccda6 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 16:55:08 +0100 Subject: [PATCH 27/52] CBOR: Handle rotation_lattice having no set value --- frame_serialize/CBORStream2Serializer.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index c34872f7..97a2f8bd 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -675,7 +675,8 @@ void CBORStream2Serializer::SerializeSequenceEnd(const EndMessage& message) { CBOR_ENC(mapEncoder, "bkg_estimate", message.bkg_estimate); CBOR_ENC(mapEncoder, "rotation_lattice_type", message.rotation_lattice_type); - CBOR_ENC(mapEncoder, "rotation_lattice", message.rotation_lattice->GetVector()); + if (message.rotation_lattice.has_value()) + CBOR_ENC(mapEncoder, "rotation_lattice", message.rotation_lattice->GetVector()); CBOR_ENC(mapEncoder, "image_scale_factor", message.scale_factor); cborErr(cbor_encoder_close_container(&encoder, &mapEncoder)); -- 2.49.1 From a9e88a9b6ef93b35e1e792c8b9f4e284dc169464 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 17:00:27 +0100 Subject: [PATCH 28/52] AcquisiotnCounters: Add more synchronization and avoid busy waiting --- acquisition_device/AcquisitionCounters.cpp | 98 +++++++++++----------- 1 file changed, 49 insertions(+), 49 deletions(-) diff --git a/acquisition_device/AcquisitionCounters.cpp b/acquisition_device/AcquisitionCounters.cpp index ae9a84f9..5cb94335 100644 --- a/acquisition_device/AcquisitionCounters.cpp +++ b/acquisition_device/AcquisitionCounters.cpp @@ -7,8 +7,8 @@ #include "../common/JFJochException.h" AcquisitionCounters::AcquisitionCounters() - : curr_frame_number(max_modules, 0), slowest_frame_number(0), fastest_frame_number(0), - total_packets(0), expected_frames(0), acquisition_finished(false) {} + : total_packets(0), slowest_frame_number(0), fastest_frame_number(0), + curr_frame_number(max_modules, 0), acquisition_finished(false), expected_frames(0) {} void AcquisitionCounters::Reset(const DiffractionExperiment &experiment, uint16_t data_stream) { std::unique_lock ul(m); @@ -56,32 +56,31 @@ void AcquisitionCounters::UpdateCounters(const Completion *c) { if (c->frame_number >= expected_frames) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "UpdateCounters frame number is out of bounds"); - else { - if (curr_frame_number.at(c->module_number) < c->frame_number) - curr_frame_number.at(c->module_number) = c->frame_number; - if (c->frame_number > slowest_frame_number) { - slowest_frame_number = curr_frame_number[0]; - for (int i = 1; i < nmodules; i++) - if (curr_frame_number[i] < slowest_frame_number) - slowest_frame_number = curr_frame_number[i]; - } + if (curr_frame_number.at(c->module_number) < c->frame_number) + curr_frame_number.at(c->module_number) = c->frame_number; - if (c->frame_number > fastest_frame_number) { - fastest_frame_number = c->frame_number; - } - - if (fastest_frame_number - slowest_frame_number > ThresholdFramesLost) - slowest_frame_number = fastest_frame_number - ThresholdFramesLost; - - packets_collected.at(c->frame_number * nmodules + c->module_number) = c->packet_count; - handle_for_frame.at(c->frame_number * nmodules + c->module_number) = c->handle; - saved_completions.at(c->frame_number * nmodules + c->module_number) = *c; - - total_packets += c->packet_count; - packets_per_module[c->module_number] += c->packet_count; + if (c->frame_number > slowest_frame_number) { + slowest_frame_number = curr_frame_number[0]; + for (int i = 1; i < nmodules; i++) + if (curr_frame_number[i] < slowest_frame_number) + slowest_frame_number = curr_frame_number[i]; } + if (c->frame_number > fastest_frame_number) { + fastest_frame_number = c->frame_number; + } + + if (fastest_frame_number - slowest_frame_number > ThresholdFramesLost) + slowest_frame_number = fastest_frame_number - ThresholdFramesLost; + + packets_collected.at(c->frame_number * nmodules + c->module_number) = c->packet_count; + handle_for_frame.at(c->frame_number * nmodules + c->module_number) = c->handle; + saved_completions.at(c->frame_number * nmodules + c->module_number) = *c; + + total_packets += c->packet_count; + packets_per_module[c->module_number] += c->packet_count; + data_updated.notify_all(); } @@ -132,44 +131,50 @@ uint64_t AcquisitionCounters::GetCurrFrameNumber(uint16_t module_number) const { if (module_number >= max_modules) throw JFJochException(JFJochExceptionCategory::ArrayOutOfBounds, "GetCurrFrameNumber Wrong module number: " + std::to_string(module_number)); + std::shared_lock sl(m); return curr_frame_number[module_number]; } uint64_t AcquisitionCounters::GetSlowestFrameNumber() const { + std::shared_lock sl(m); return slowest_frame_number; } uint64_t AcquisitionCounters::GetFastestFrameNumber() const { + std::shared_lock sl(m); return fastest_frame_number; } void AcquisitionCounters::WaitForFrame(size_t curr_frame, uint16_t module_number) const { - uint64_t slowest_head_tmp = (module_number == UINT16_MAX) ? GetSlowestFrameNumber() : GetCurrFrameNumber(module_number); if (curr_frame == 0) - curr_frame = 1; // Cannot wait for frame 0, as this is initial value of slowest frame number, so waiting for frame 1 in this case - while (!acquisition_finished && (slowest_head_tmp < curr_frame)) { - std::this_thread::sleep_for(std::chrono::microseconds(100)); - slowest_head_tmp = (module_number == UINT16_MAX) ? GetSlowestFrameNumber() : GetCurrFrameNumber(module_number); - } + curr_frame = 1; // Cannot wait for frame 0, as this is initial value + + std::shared_lock sl(m); + data_updated.wait(sl, [&] { + const uint64_t head = (module_number == UINT16_MAX) + ? slowest_frame_number + : curr_frame_number.at(module_number); + return acquisition_finished || (head >= curr_frame); + }); } int64_t AcquisitionCounters::CalculateDelay(size_t curr_frame, uint16_t module_number) const { - uint64_t slowest_head_tmp; - - if (module_number == UINT16_MAX) - slowest_head_tmp = GetSlowestFrameNumber(); - else - slowest_head_tmp = GetCurrFrameNumber(module_number); - if (slowest_head_tmp < curr_frame) + std::shared_lock sl(m); + const uint64_t head = (module_number == UINT16_MAX) + ? slowest_frame_number + : curr_frame_number.at(module_number); + if (head < curr_frame) return 0; - return slowest_head_tmp - curr_frame; + return static_cast(head - curr_frame); } bool AcquisitionCounters::IsAcquisitionFinished() const { + std::shared_lock sl(m); return acquisition_finished; } bool AcquisitionCounters::IsFullModuleCollected(size_t frame, uint16_t module_number) const { + std::shared_lock sl(m); if (frame >= expected_frames) throw JFJochException(JFJochExceptionCategory::ArrayOutOfBounds, "IsFullModuleCollected Wrong frame number: " + std::to_string(frame)); @@ -183,6 +188,7 @@ bool AcquisitionCounters::IsFullModuleCollected(size_t frame, uint16_t module_nu } bool AcquisitionCounters::IsAnyPacketCollected(size_t frame, uint16_t module_number) const { + std::shared_lock sl(m); if (frame >= expected_frames) throw JFJochException(JFJochExceptionCategory::ArrayOutOfBounds, "IsFullModuleCollected Wrong frame number: " + std::to_string(frame)); @@ -196,19 +202,10 @@ bool AcquisitionCounters::IsAnyPacketCollected(size_t frame, uint16_t module_num } uint64_t AcquisitionCounters::GetTotalPackets() const { + std::shared_lock sl(m); return total_packets; } -uint64_t AcquisitionCounters::GetTotalPackets(uint16_t module_number) const { - std::unique_lock ul(m); - - if (module_number >= nmodules) - throw JFJochException(JFJochExceptionCategory::ArrayOutOfBounds, - "GetTotalPackets Wrong module number: " + std::to_string(module_number)); - - return packets_per_module[module_number]; -} - uint64_t AcquisitionCounters::GetExpectedPacketsPerImage() const { return expected_packets_per_module * nmodules; } @@ -218,13 +215,16 @@ uint64_t AcquisitionCounters::GetTotalExpectedPackets() const { } uint64_t AcquisitionCounters::GetTotalExpectedPacketsPerModule() const { + std::shared_lock sl(m); return expected_frames * expected_packets_per_module; } uint64_t AcquisitionCounters::GetModuleNumber() const { + std::shared_lock sl(m); return nmodules; } uint64_t AcquisitionCounters::GetBytesReceived() const { - return GetTotalPackets() * bytes_per_packet; + std::shared_lock sl(m); + return total_packets * bytes_per_packet; } -- 2.49.1 From ab077625b90be77cb1b27aead1e089fc25c03311 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 17:02:27 +0100 Subject: [PATCH 29/52] StatusVector: Handle (hypothetical) situation when one image is overwritten by itself --- common/StatusVector.cpp | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/common/StatusVector.cpp b/common/StatusVector.cpp index a0af404a..112d60d5 100644 --- a/common/StatusVector.cpp +++ b/common/StatusVector.cpp @@ -25,11 +25,23 @@ void StatusVector::AddElement(uint32_t id, float val) { if (id >= content.size()) { content.resize(id + 1, NAN); } + + const float old_val = content[id]; + const bool had_old = std::isfinite(old_val); + content[id] = val; - sum += val; - count += 1; - mean = sum / count; + if (had_old) { + sum += static_cast(val) - static_cast(old_val); + } else { + sum += val; + count += 1; + } + + if (count > 0) + mean = static_cast(sum / static_cast(count)); + else + mean = NAN; } std::optional StatusVector::GetElement(uint32_t id) const { -- 2.49.1 From fe61cfcfdecfb3424cbdf88e5863d63e7799a689 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 17:32:34 +0100 Subject: [PATCH 30/52] AcquisitionDevice: Fix potential infinite loop --- acquisition_device/AcquisitionDevice.cpp | 57 +++++++++++++----------- 1 file changed, 30 insertions(+), 27 deletions(-) diff --git a/acquisition_device/AcquisitionDevice.cpp b/acquisition_device/AcquisitionDevice.cpp index 14262055..a087b2aa 100644 --- a/acquisition_device/AcquisitionDevice.cpp +++ b/acquisition_device/AcquisitionDevice.cpp @@ -81,44 +81,47 @@ void AcquisitionDevice::WaitForActionComplete() { while (c.type != Completion::Type::End) { DeviceOutput* output; + bool skip_analysis = false; + try { output = GetDeviceOutput(c.handle); } catch (const JFJochException &e) { if (logger) logger->ErrorException(e); - continue; + skip_analysis = true; } - c.module_number = output->module_statistics.module_number; - c.packet_count = output->module_statistics.packet_count; - c.frame_number = output->module_statistics.frame_number; + if (!skip_analysis) { + c.module_number = output->module_statistics.module_number; + c.packet_count = output->module_statistics.packet_count; + c.frame_number = output->module_statistics.frame_number; - if (c.frame_number >= expected_frames) { - Cancel(); - // this frame is not of any interest, therefore its location can be immediately released - SendWorkRequest(c.handle); - } else if (c.module_number >= max_modules) { - // Module number out of bounds, don't process - if (logger != nullptr) - logger->Error("Completion with wrong module number data stream {} completion frame number {} module {} handle {}", - data_stream, c.frame_number, c.module_number, c.handle); - SendWorkRequest(c.handle); - } else if (c.frame_number < counters.GetSlowestFrameNumber()) { - // Module is falling behind, needs to return the handle then - SendWorkRequest(c.handle); - } else { - try { - counters.UpdateCounters(&c); - } catch (const JFJochException &e) { - if (logger) - logger->ErrorException(e); + if (c.frame_number >= expected_frames) { + Cancel(); + // this frame is not of any interest, therefore its location can be immediately released SendWorkRequest(c.handle); + } else if (c.module_number >= max_modules) { + // Module number out of bounds, don't process + if (logger != nullptr) + logger->Error("Completion with wrong module number data stream {} completion frame number {} module {} handle {}", + data_stream, c.frame_number, c.module_number, c.handle); + SendWorkRequest(c.handle); + } else if (c.frame_number < counters.GetSlowestFrameNumber()) { + // Module is falling behind, needs to return the handle then + SendWorkRequest(c.handle); + } else { + try { + counters.UpdateCounters(&c); + } catch (const JFJochException &e) { + if (logger) + logger->ErrorException(e); + SendWorkRequest(c.handle); + } } + if (logger != nullptr) + logger->Debug("Data stream {} completion frame number {} module {} handle {}", + data_stream, c.frame_number, c.module_number, c.handle); } - if (logger != nullptr) - logger->Debug("Data stream {} completion frame number {} module {} handle {}", - data_stream, c.frame_number, c.module_number, c.handle); - c = work_completion_queue.GetBlocking(); } counters.SetAcquisitionFinished(); -- 2.49.1 From 8691d0e83fa2a31e81d68281ad2c57ae05b0fc6c Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 17:32:54 +0100 Subject: [PATCH 31/52] AcquisitionCounters: Fixes --- acquisition_device/AcquisitionCounters.cpp | 10 ++++++++++ acquisition_device/AcquisitionCounters.h | 2 +- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/acquisition_device/AcquisitionCounters.cpp b/acquisition_device/AcquisitionCounters.cpp index 5cb94335..1ef54d1c 100644 --- a/acquisition_device/AcquisitionCounters.cpp +++ b/acquisition_device/AcquisitionCounters.cpp @@ -206,6 +206,16 @@ uint64_t AcquisitionCounters::GetTotalPackets() const { return total_packets; } +uint64_t AcquisitionCounters::GetTotalPackets(uint16_t module_number) const { + std::unique_lock ul(m); + + if (module_number >= nmodules) + throw JFJochException(JFJochExceptionCategory::ArrayOutOfBounds, + "GetTotalPackets Wrong module number: " + std::to_string(module_number)); + + return packets_per_module[module_number]; +} + uint64_t AcquisitionCounters::GetExpectedPacketsPerImage() const { return expected_packets_per_module * nmodules; } diff --git a/acquisition_device/AcquisitionCounters.h b/acquisition_device/AcquisitionCounters.h index 435d9e7a..622f2a06 100644 --- a/acquisition_device/AcquisitionCounters.h +++ b/acquisition_device/AcquisitionCounters.h @@ -16,7 +16,7 @@ // so uses mutex to ensure consistency class AcquisitionCounters { - uint16_t expected_packets_per_module; + uint32_t expected_packets_per_module; constexpr static const uint64_t max_modules = 32; mutable std::shared_mutex m; mutable std::condition_variable_any data_updated; -- 2.49.1 From c2d930a5999e7588c3b35321b3ea73087b28a19e Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 20:07:44 +0100 Subject: [PATCH 32/52] ZMQWrappers: Make SendZeroCopy more robust for error handling --- common/ZMQWrappers.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/common/ZMQWrappers.cpp b/common/ZMQWrappers.cpp index 0378165a..59d83fad 100644 --- a/common/ZMQWrappers.cpp +++ b/common/ZMQWrappers.cpp @@ -91,12 +91,13 @@ void ZMQSocket::SendZeroCopy(void *data, size_t size,void (*callback)(void *, vo zmq_msg_t msg; if (zmq_msg_init_data(&msg, data, size, callback, hint) != 0) { callback(data, hint); - throw JFJochException(JFJochExceptionCategory::ZeroMQ, "zmq_msg_init_data failed"); - } - if (zmq_msg_send(&msg, socket, 0) < 0) { - callback(data, hint); - throw JFJochException(JFJochExceptionCategory::ZeroMQ, "zmq_msg_send failed"); + throw JFJochException(JFJochExceptionCategory::ZeroMQ, "zmq_msg_init_data failed " + std::string(std::strerror(errno))); } + + auto result = zmq_msg_send(&msg, socket, 0); + zmq_msg_close(&msg); + if (result < 0) + throw JFJochException(JFJochExceptionCategory::ZeroMQ, "zmq_msg_send failed " + std::string(std::strerror(errno))); } void ZMQSocket::SetSocketOption(int32_t option_name, int32_t value) { -- 2.49.1 From 7f416b09568ba5e47ad51d050f09e20cd923134d Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 20:31:58 +0100 Subject: [PATCH 33/52] ZMQWrappers: SendZeroCopy can be blocking --- common/ZMQWrappers.cpp | 14 ++++++++++---- common/ZMQWrappers.h | 2 +- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/common/ZMQWrappers.cpp b/common/ZMQWrappers.cpp index 59d83fad..a969838c 100644 --- a/common/ZMQWrappers.cpp +++ b/common/ZMQWrappers.cpp @@ -86,7 +86,7 @@ void ZMQSocket::Send(zmq_msg_t *msg) { throw JFJochException(JFJochExceptionCategory::ZeroMQ, "zmq_msg_send failed"); } -void ZMQSocket::SendZeroCopy(void *data, size_t size,void (*callback)(void *, void *), void *hint) { +bool ZMQSocket::SendZeroCopy(void *data, size_t size,void (*callback)(void *, void *), void *hint, bool blocking) { std::unique_lock ul(m); zmq_msg_t msg; if (zmq_msg_init_data(&msg, data, size, callback, hint) != 0) { @@ -94,10 +94,16 @@ void ZMQSocket::SendZeroCopy(void *data, size_t size,void (*callback)(void *, vo throw JFJochException(JFJochExceptionCategory::ZeroMQ, "zmq_msg_init_data failed " + std::string(std::strerror(errno))); } - auto result = zmq_msg_send(&msg, socket, 0); + auto result = zmq_msg_send(&msg, socket, (blocking ? 0 : ZMQ_DONTWAIT)); zmq_msg_close(&msg); - if (result < 0) - throw JFJochException(JFJochExceptionCategory::ZeroMQ, "zmq_msg_send failed " + std::string(std::strerror(errno))); + + if (result == 0) + return true; + + if ((errno == EAGAIN) || (errno == EINTR)) + return false; + + throw JFJochException(JFJochExceptionCategory::ZeroMQ, "zmq_msg_send failed " + std::string(std::strerror(errno))); } void ZMQSocket::SetSocketOption(int32_t option_name, int32_t value) { diff --git a/common/ZMQWrappers.h b/common/ZMQWrappers.h index bca00748..97730ad8 100644 --- a/common/ZMQWrappers.h +++ b/common/ZMQWrappers.h @@ -67,7 +67,7 @@ public: void Send(); bool Send(const void *buf, size_t buf_size, bool blocking = true, bool multipart = false); - void SendZeroCopy(void *data, size_t size,void (*callback)(void *, void *), void *hint); + bool SendZeroCopy(void *data, size_t size,void (*callback)(void *, void *), void *hint, bool blocking = true); template bool Send(const std::vector &buf) { return Send(buf.data(), buf.size() * sizeof(T)); } -- 2.49.1 From c92afa95fbdc69ea317eef4536c3fd46c0be2eb3 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 20:59:03 +0100 Subject: [PATCH 34/52] ZeroCopyReturnValue: Add GetHandle() --- common/ZeroCopyReturnValue.cpp | 4 ++++ common/ZeroCopyReturnValue.h | 1 + 2 files changed, 5 insertions(+) diff --git a/common/ZeroCopyReturnValue.cpp b/common/ZeroCopyReturnValue.cpp index 9e7f3e85..5f376f7a 100644 --- a/common/ZeroCopyReturnValue.cpp +++ b/common/ZeroCopyReturnValue.cpp @@ -42,3 +42,7 @@ void ZeroCopyReturnValue::SetIndexed(bool input) { bool ZeroCopyReturnValue::IsIndexed() const { return indexed; } + +uint32_t ZeroCopyReturnValue::GetHandle() const { + return handle; +} diff --git a/common/ZeroCopyReturnValue.h b/common/ZeroCopyReturnValue.h index 1f774407..11e5de0d 100644 --- a/common/ZeroCopyReturnValue.h +++ b/common/ZeroCopyReturnValue.h @@ -23,6 +23,7 @@ public: void SetIndexed(bool input); void *GetImage() const; + [[nodiscard]] uint32_t GetHandle() const; [[nodiscard]] size_t GetImageSize() const; [[nodiscard]] int64_t GetImageNumber() const; [[nodiscard]] bool IsIndexed() const; -- 2.49.1 From c897008d6e1f196d7eae6229eb71b1d453c3bc67 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 28 Feb 2026 21:22:58 +0100 Subject: [PATCH 35/52] ImageBuffer: Improved logic for slot management --- broker/OpenAPIConvert.cpp | 6 + broker/gen/model/Image_buffer_status.cpp | 30 +++- broker/gen/model/Image_buffer_status.h | 14 ++ broker/jfjoch_api.yaml | 12 ++ broker/redoc-static.html | 6 +- common/ImageBuffer.cpp | 134 +++++++++++++----- common/ImageBuffer.h | 25 ++-- common/ZeroCopyReturnValue.cpp | 24 +++- common/ZeroCopyReturnValue.h | 5 +- docs/python_client/docs/ImageBufferStatus.md | 2 + .../src/openapi/models/image_buffer_status.ts | 10 ++ receiver/JFJochReceiverFPGA.cpp | 1 + receiver/JFJochReceiverLite.cpp | 1 + tests/ImageBufferTest.cpp | 29 +++- 14 files changed, 242 insertions(+), 57 deletions(-) diff --git a/broker/OpenAPIConvert.cpp b/broker/OpenAPIConvert.cpp index 13576b9d..3b603dcc 100644 --- a/broker/OpenAPIConvert.cpp +++ b/broker/OpenAPIConvert.cpp @@ -753,6 +753,9 @@ org::openapitools::server::model::Image_buffer_status Convert(const ImageBufferS ret.setImageNumbers(input.images_in_the_buffer); ret.setMaxImageNumber(input.max_image_number); ret.setMinImageNumber(input.min_image_number); + ret.setInPreparationSlots(input.preparation_slots); + ret.setInSendingSlots(input.sending_slots); + if (input.current_counter.has_value()) ret.setCurrentCounter(input.current_counter.value()); return ret; @@ -765,6 +768,9 @@ ImageBufferStatus Convert(const org::openapitools::server::model::Image_buffer_s ret.images_in_the_buffer = input.getImageNumbers(); ret.max_image_number = input.getMaxImageNumber(); ret.min_image_number = input.getMinImageNumber(); + ret.sending_slots = input.getInSendingSlots(); + ret.preparation_slots = input.getInPreparationSlots(); + if (input.currentCounterIsSet()) ret.current_counter = input.getCurrentCounter(); return ret; diff --git a/broker/gen/model/Image_buffer_status.cpp b/broker/gen/model/Image_buffer_status.cpp index 198ae9f2..9e415bf2 100644 --- a/broker/gen/model/Image_buffer_status.cpp +++ b/broker/gen/model/Image_buffer_status.cpp @@ -25,6 +25,8 @@ Image_buffer_status::Image_buffer_status() m_Max_image_number = 0L; m_Total_slots = 0L; m_Available_slots = 0L; + m_In_preparation_slots = 0L; + m_In_sending_slots = 0L; m_Current_counter = 0L; m_Current_counterIsSet = false; @@ -98,7 +100,7 @@ bool Image_buffer_status::validate(std::stringstream& msg, const std::string& pa } } - + return success; } @@ -122,6 +124,12 @@ bool Image_buffer_status::operator==(const Image_buffer_status& rhs) const (getAvailableSlots() == rhs.getAvailableSlots()) && + (getInPreparationSlots() == rhs.getInPreparationSlots()) + && + + (getInSendingSlots() == rhs.getInSendingSlots()) + && + ((!currentCounterIsSet() && !rhs.currentCounterIsSet()) || (currentCounterIsSet() && rhs.currentCounterIsSet() && getCurrentCounter() == rhs.getCurrentCounter())) @@ -141,6 +149,8 @@ void to_json(nlohmann::json& j, const Image_buffer_status& o) j["image_numbers"] = o.m_Image_numbers; j["total_slots"] = o.m_Total_slots; j["available_slots"] = o.m_Available_slots; + j["in_preparation_slots"] = o.m_In_preparation_slots; + j["in_sending_slots"] = o.m_In_sending_slots; if(o.currentCounterIsSet()) j["current_counter"] = o.m_Current_counter; @@ -153,6 +163,8 @@ void from_json(const nlohmann::json& j, Image_buffer_status& o) j.at("image_numbers").get_to(o.m_Image_numbers); j.at("total_slots").get_to(o.m_Total_slots); j.at("available_slots").get_to(o.m_Available_slots); + j.at("in_preparation_slots").get_to(o.m_In_preparation_slots); + j.at("in_sending_slots").get_to(o.m_In_sending_slots); if(j.find("current_counter") != j.end()) { j.at("current_counter").get_to(o.m_Current_counter); @@ -201,6 +213,22 @@ void Image_buffer_status::setAvailableSlots(int64_t const value) { m_Available_slots = value; } +int64_t Image_buffer_status::getInPreparationSlots() const +{ + return m_In_preparation_slots; +} +void Image_buffer_status::setInPreparationSlots(int64_t const value) +{ + m_In_preparation_slots = value; +} +int64_t Image_buffer_status::getInSendingSlots() const +{ + return m_In_sending_slots; +} +void Image_buffer_status::setInSendingSlots(int64_t const value) +{ + m_In_sending_slots = value; +} int64_t Image_buffer_status::getCurrentCounter() const { return m_Current_counter; diff --git a/broker/gen/model/Image_buffer_status.h b/broker/gen/model/Image_buffer_status.h index 84337033..e38c6b53 100644 --- a/broker/gen/model/Image_buffer_status.h +++ b/broker/gen/model/Image_buffer_status.h @@ -84,6 +84,16 @@ public: int64_t getAvailableSlots() const; void setAvailableSlots(int64_t const value); /// + /// Number of slots in the image buffer that are currently in preparation for sending. + /// + int64_t getInPreparationSlots() const; + void setInPreparationSlots(int64_t const value); + /// + /// Number of slots in the image buffer that are currently sending/writing data. + /// + int64_t getInSendingSlots() const; + void setInSendingSlots(int64_t const value); + /// /// Counter of changes in the image buffer - either new start message or new image added. For optimization one can only load new images/datasets from the HTTP if this value changes. Counter is optional as it was not implemented in older versions to avoid breaking change /// int64_t getCurrentCounter() const; @@ -104,6 +114,10 @@ protected: int64_t m_Available_slots; + int64_t m_In_preparation_slots; + + int64_t m_In_sending_slots; + int64_t m_Current_counter; bool m_Current_counterIsSet; diff --git a/broker/jfjoch_api.yaml b/broker/jfjoch_api.yaml index 25329822..1d520b2f 100644 --- a/broker/jfjoch_api.yaml +++ b/broker/jfjoch_api.yaml @@ -695,6 +695,8 @@ components: required: - image_numbers - total_slots + - in_preparation_slots + - in_sending_slots - available_slots - max_image_number - min_image_number @@ -725,6 +727,16 @@ components: type: integer format: int64 description: Slots available for the data collection + in_preparation_slots: + type: integer + format: int64 + description: | + Number of slots in the image buffer that are currently in preparation for sending. + in_sending_slots: + type: integer + format: int64 + description: | + Number of slots in the image buffer that are currently sending/writing data. current_counter: type: integer format: int64 diff --git a/broker/redoc-static.html b/broker/redoc-static.html index f66e2520..133d60e3 100644 --- a/broker/redoc-static.html +++ b/broker/redoc-static.html @@ -790,7 +790,7 @@ This can only be done when detector is Idle, Error or
http://localhost:5232/config/roi

Request samples

Content type
application/json
{
  • "box": {
    },
  • "circle": {
    },
  • "azim": {
    }
}

Response samples

Content type
application/json
{
  • "msg": "Detector in wrong state",
  • "reason": "WrongDAQState"
}

Get general statistics

query Parameters
compression
boolean
Default: false

Enable DEFLATE compression of output data.

Responses

Response samples

Content type
application/json
{
  • "detector": {
    },
  • "detector_list": {
    },
  • "detector_settings": {
    },
  • "image_format_settings": {
    },
  • "instrument_metadata": {
    },
  • "file_writer_settings": {
    },
  • "data_processing_settings": {
    },
  • "measurement": {
    },
  • "broker": {
    },
  • "fpga": [
    ],
  • "calibration": [
    ],
  • "zeromq_preview": {
    },
  • "zeromq_metadata": {
    },
  • "dark_mask": {
    },
  • "pixel_mask": {
    },
  • "roi": {
    },
  • "az_int": {
    },
  • "buffer": {
    },
  • "indexing": {
    }
}

Get data collection statistics

Results of the last data collection

+
http://localhost:5232/statistics

Response samples

Content type
application/json
{
  • "detector": {
    },
  • "detector_list": {
    },
  • "detector_settings": {
    },
  • "image_format_settings": {
    },
  • "instrument_metadata": {
    },
  • "file_writer_settings": {
    },
  • "data_processing_settings": {
    },
  • "measurement": {
    },
  • "broker": {
    },
  • "fpga": [
    ],
  • "calibration": [
    ],
  • "zeromq_preview": {
    },
  • "zeromq_metadata": {
    },
  • "dark_mask": {
    },
  • "pixel_mask": {
    },
  • "roi": {
    },
  • "az_int": {
    },
  • "buffer": {
    },
  • "indexing": {
    }
}

Get data collection statistics

Results of the last data collection

Responses

Response samples

Content type
application/json
{
  • "min_image_number": 0,
  • "max_image_number": 0,
  • "image_numbers": [
    ],
  • "total_slots": 0,
  • "available_slots": 0,
  • "current_counter": 0
}

Get Jungfraujoch version of jfjoch_broker

Responses

Response samples

Content type
application/json
{
  • "min_image_number": 0,
  • "max_image_number": 0,
  • "image_numbers": [
    ],
  • "total_slots": 0,
  • "available_slots": 0,
  • "in_preparation_slots": 0,
  • "in_sending_slots": 0,
  • "current_counter": 0
}

Get Jungfraujoch version of jfjoch_broker

Responses