From e6e0f5b2f0d420cfaddbc01e6f07da9d0ce7dd99 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 20 Feb 2026 11:37:59 +0100 Subject: [PATCH] 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); } }