From a5e582574e6fa3459c9dbaa1e80c4c8609b89b82 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 30 Jan 2026 14:44:54 +0100 Subject: [PATCH] IndexAndRefine: Use small angle approximation for optimizing rotation data --- image_analysis/IndexAndRefine.cpp | 71 ++++++++++--------- image_analysis/geom_refinement/CMakeLists.txt | 4 +- .../SimpleRotXtalOptimizer.cpp | 56 +++++++++++---- .../geom_refinement/SimpleRotXtalOptimizer.h | 2 +- 4 files changed, 84 insertions(+), 49 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index c0f9b933..7b182fdf 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -9,7 +9,7 @@ #include "indexing/AnalyzeIndexing.h" #include "indexing/FFTIndexer.h" #include "lattice_search/LatticeSearch.h" - +#include "geom_refinement/SimpleRotXtalOptimizer.h" IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer) : index_ice_rings(x.GetIndexingSettings().GetIndexIceRings()), @@ -93,42 +93,49 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In if (!outcome.lattice_candidate) return; - XtalOptimizerData data{ - .geom = outcome.experiment.GetDiffractionGeometry(), - .latt = *outcome.lattice_candidate, - .crystal_system = outcome.symmetry.crystal_system, - .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), - .refine_beam_center = true, - .refine_distance_mm = false, - .refine_detector_angles = false, - .max_time = 0.04 // 40 ms is max allowed time for the operation - }; - if (experiment.IsRotationIndexing()) { - data.refine_beam_center = false; - data.refine_unit_cell = false; - } + SimpleRotXtalOptimizerData data{ + .geom = outcome.experiment.GetDiffractionGeometry(), + .latt = *outcome.lattice_candidate, + .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), + }; + if (SimpleRotXtalOptimizer(data, msg.spots)) { + outcome.lattice_candidate = data.latt; + } + } else { + XtalOptimizerData data{ + .geom = outcome.experiment.GetDiffractionGeometry(), + .latt = *outcome.lattice_candidate, + .crystal_system = outcome.symmetry.crystal_system, + .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), + .refine_beam_center = true, + .refine_distance_mm = false, + .refine_detector_angles = false, + .max_time = 0.04 // 40 ms is max allowed time for the operation + }; - if (outcome.symmetry.crystal_system == gemmi::CrystalSystem::Trigonal) - data.crystal_system = gemmi::CrystalSystem::Hexagonal; - switch (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm()) { - case GeomRefinementAlgorithmEnum::None: - break; - case GeomRefinementAlgorithmEnum::BeamCenter: - if (XtalOptimizer(data, msg.spots)) { - outcome.experiment.BeamX_pxl(data.geom.GetBeamX_pxl()) - .BeamY_pxl(data.geom.GetBeamY_pxl()); - outcome.beam_center_updated = true; - } - break; - } + if (outcome.symmetry.crystal_system == gemmi::CrystalSystem::Trigonal) + data.crystal_system = gemmi::CrystalSystem::Hexagonal; - outcome.lattice_candidate = data.latt; + switch (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm()) { + case GeomRefinementAlgorithmEnum::None: + break; + case GeomRefinementAlgorithmEnum::BeamCenter: + if (XtalOptimizer(data, msg.spots)) { + outcome.experiment.BeamX_pxl(data.geom.GetBeamX_pxl()) + .BeamY_pxl(data.geom.GetBeamY_pxl()); + outcome.beam_center_updated = true; + } + break; + } - if (outcome.beam_center_updated) { - msg.beam_corr_x = data.beam_corr_x; - msg.beam_corr_y = data.beam_corr_y; + outcome.lattice_candidate = data.latt; + + if (outcome.beam_center_updated) { + msg.beam_corr_x = data.beam_corr_x; + msg.beam_corr_y = data.beam_corr_y; + } } } diff --git a/image_analysis/geom_refinement/CMakeLists.txt b/image_analysis/geom_refinement/CMakeLists.txt index 0143b49a..71bdd41e 100644 --- a/image_analysis/geom_refinement/CMakeLists.txt +++ b/image_analysis/geom_refinement/CMakeLists.txt @@ -6,8 +6,8 @@ ADD_LIBRARY(JFJochGeomRefinement STATIC AssignSpotsToRings.h XtalOptimizer.cpp XtalOptimizer.h - RotOptimizer.cpp - RotOptimizer.h + SimpleRotXtalOptimizer.cpp + SimpleRotXtalOptimizer.h ) TARGET_LINK_LIBRARIES(JFJochGeomRefinement Ceres::ceres Eigen3::Eigen JFJochCommon) diff --git a/image_analysis/geom_refinement/SimpleRotXtalOptimizer.cpp b/image_analysis/geom_refinement/SimpleRotXtalOptimizer.cpp index 6dfc47f4..2b7cc017 100644 --- a/image_analysis/geom_refinement/SimpleRotXtalOptimizer.cpp +++ b/image_analysis/geom_refinement/SimpleRotXtalOptimizer.cpp @@ -8,7 +8,8 @@ static bool SimpleRotXtalOptimizerInternal( SimpleRotXtalOptimizerData &data, const std::vector &spots, - float tolerance) + float tolerance, + Eigen::Matrix3d &r_total) { Coord vec0 = data.latt.Vec0(); Coord vec1 = data.latt.Vec1(); @@ -59,10 +60,11 @@ static bool SimpleRotXtalOptimizerInternal( const auto &h = obs[i]; const auto &e = exp[i]; - // Cross-product matrix of h_obs - A.row(3*i + 0) << 0.0, -h.z(), h.y(); - A.row(3*i + 1) << h.z(), 0.0, -h.x(); - A.row(3*i + 2) << -h.y(), h.x(), 0.0; + // Cross-product matrix of h_obs: [h]× such that [h]× · ω = h × ω + // But we need ω × h = -h × ω, so use -[h]× + A.row(3*i + 0) << 0.0, h.z(), -h.y(); + A.row(3*i + 1) << -h.z(), 0.0, h.x(); + A.row(3*i + 2) << h.y(), -h.x(), 0.0; b.segment<3>(3*i) = e - h; } @@ -74,21 +76,47 @@ static bool SimpleRotXtalOptimizerInternal( if (!omega.allFinite()) return false; - if (omega.norm() > 0.05) { - return false; - // > ~3 degrees — warn or reject - } + // Apply incremental rotation and accumulate it + const double angle = omega.norm(); + Eigen::Matrix3d r_inc = Eigen::Matrix3d::Identity(); + if (angle > 0.0) + r_inc = Eigen::AngleAxisd(angle, omega / angle).toRotationMatrix(); - // Apply small-angle rotation to lattice - // TODO: data.latt.RotateSmallAngle(omega); + r_total = r_inc * r_total; + + if (angle > 0.0) { + // Build the current lattice matrix (columns = a, b, c) + Eigen::Matrix3d L; + L.col(0) = Eigen::Vector3d(vec0.x, vec0.y, vec0.z); + L.col(1) = Eigen::Vector3d(vec1.x, vec1.y, vec1.z); + L.col(2) = Eigen::Vector3d(vec2.x, vec2.y, vec2.z); + + // The HKL-space rotation R transforms: h_new = R * h_old + // This corresponds to: L_new = L_old * R^T (real-space lattice) + Eigen::Matrix3d L_new = L * r_inc.transpose(); + + data.latt = CrystalLattice( + Coord(L_new(0,0), L_new(1,0), L_new(2,0)), + Coord(L_new(0,1), L_new(1,1), L_new(2,1)), + Coord(L_new(0,2), L_new(1,2), L_new(2,2)) + ); + } return true; } bool SimpleRotXtalOptimizer(SimpleRotXtalOptimizerData &data, const std::vector &spots) { + Eigen::Matrix3d r_total = Eigen::Matrix3d::Identity(); - if (!SimpleRotXtalOptimizerInternal(data, spots, 0.3)) + if (!SimpleRotXtalOptimizerInternal(data, spots, 0.3, r_total)) return false; - SimpleRotXtalOptimizerInternal(data, spots, 0.2); - return SimpleRotXtalOptimizerInternal(data, spots, 0.1); + SimpleRotXtalOptimizerInternal(data, spots, 0.2, r_total); + const bool ok = SimpleRotXtalOptimizerInternal(data, spots, 0.1, r_total); + + Eigen::AngleAxisd aa(r_total); + data.rotation[0] = aa.axis().x(); + data.rotation[1] = aa.axis().y(); + data.rotation[2] = aa.axis().z(); + data.angle = aa.angle() * 180.0 / M_PI; + return ok; } diff --git a/image_analysis/geom_refinement/SimpleRotXtalOptimizer.h b/image_analysis/geom_refinement/SimpleRotXtalOptimizer.h index 2587e0cb..b158f8e6 100644 --- a/image_analysis/geom_refinement/SimpleRotXtalOptimizer.h +++ b/image_analysis/geom_refinement/SimpleRotXtalOptimizer.h @@ -15,8 +15,8 @@ struct SimpleRotXtalOptimizerData { int64_t min_spots = 8; bool index_ice_rings = true; - float max_time = 1.0; double rotation[3] = {0,0,0}; + double angle = 0.0; }; bool SimpleRotXtalOptimizer(SimpleRotXtalOptimizerData &data, const std::vector &spots);