diff --git a/image_analysis/indexing/CMakeLists.txt b/image_analysis/indexing/CMakeLists.txt index bea0431b..f40c1492 100644 --- a/image_analysis/indexing/CMakeLists.txt +++ b/image_analysis/indexing/CMakeLists.txt @@ -13,8 +13,10 @@ ADD_LIBRARY(JFJochIndexing STATIC FFTResult.h FFTIndexer.cpp FFTIndexer.h - PostIndexingRefinement.cpp) -TARGET_LINK_LIBRARIES(JFJochIndexing JFJochCommon) + PostIndexingRefinement.cpp + MultiLatticeSearch.cpp + MultiLatticeSearch.h) +TARGET_LINK_LIBRARIES(JFJochIndexing JFJochCommon JFJochLatticeSearch) IF (JFJOCH_CUDA_AVAILABLE) FetchContent_Declare( diff --git a/image_analysis/indexing/MultiLatticeSearch.cpp b/image_analysis/indexing/MultiLatticeSearch.cpp new file mode 100644 index 00000000..c799d9ab --- /dev/null +++ b/image_analysis/indexing/MultiLatticeSearch.cpp @@ -0,0 +1,99 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "MultiLatticeSearch.h" + +#include + +namespace { + // Direct lattice vectors as matrix columns: M = [a | b | c] + Eigen::Matrix3d LatticeMatrix(const CrystalLattice &latt) { + const Coord a = latt.Vec0(); + const Coord b = latt.Vec1(); + const Coord c = latt.Vec2(); + + Eigen::Matrix3d M; + M.col(0) = Eigen::Vector3d(a.x, a.y, a.z); + M.col(1) = Eigen::Vector3d(b.x, b.y, b.z); + M.col(2) = Eigen::Vector3d(c.x, c.y, c.z); + return M; + } + + // Proper rotation mapping reference -> target (target = R * reference). + // R = Mtarget * Mref^-1, projected to the nearest rotation via SVD. + Eigen::Matrix3d RotationRefToTarget(const CrystalLattice &reference, const CrystalLattice &target) { + Eigen::Matrix3d R = LatticeMatrix(target) * LatticeMatrix(reference).inverse(); + + Eigen::JacobiSVD svd(R, Eigen::ComputeFullU | Eigen::ComputeFullV); + R = svd.matrixU() * svd.matrixV().transpose(); + if (R.determinant() < 0.0) { + Eigen::Matrix3d U = svd.matrixU(); + U.col(2) *= -1.0; + R = U * svd.matrixV().transpose(); + } + return R; + } +} + +std::vector MultiLatticeSearch(const std::vector &lattices, + float dist_tolerance, + float angle_tolerance_deg) { + std::vector ret; + if (lattices.empty()) + return ret; + + const CrystalLattice &reference = lattices[0]; + const UnitCell ref_cell = reference.GetUnitCell(); + + ret.push_back({reference, reference, Coord(0, 0, 0)}); + + for (size_t i = 1; i < lattices.size(); i++) { + const CrystalLattice &latt = lattices[i]; + + if (!latt.GetUnitCell().is_close(ref_cell, dist_tolerance, angle_tolerance_deg)) + continue; + + const Eigen::Matrix3d R = RotationRefToTarget(reference, latt); + + const Eigen::AngleAxisd aa(R); + const Eigen::Vector3d rod = aa.angle() * aa.axis(); + const Coord rotation_vector(static_cast(rod.x()), + static_cast(rod.y()), + static_cast(rod.z())); + + // output_lattice = R * reference, built straight from the Eigen matrix + // so it is exactly a proper rotation of the reference (full double precision). + const Eigen::Matrix3d out = R * LatticeMatrix(reference); + ret.push_back({ + latt, + CrystalLattice(Coord(out(0, 0), out(1, 0), out(2, 0)), + Coord(out(0, 1), out(1, 1), out(2, 1)), + Coord(out(0, 2), out(1, 2), out(2, 2))), + rotation_vector + }); + } + + return ret; +} + +std::vector MultiLatticeSearchReduced(const std::vector &lattices, + float dist_tolerance, + float angle_tolerance_deg) { + if (lattices.empty()) + return {}; + + // Reduce lattice 0 to its primitive Niggli setting; that becomes the reference basis. + std::vector reduced; + reduced.reserve(lattices.size()); + + // Reduce every candidate the same way so metric comparison and rotation + // are all done in the reduced setting. + for (const auto & lattice : lattices) + reduced.push_back(lattice.NiggliReduce()); + + auto ret = MultiLatticeSearch(reduced, dist_tolerance, angle_tolerance_deg); + + // Keep input_lattice pointing at the original (un-reduced) lattices for debugging. + + return ret; +} \ No newline at end of file diff --git a/image_analysis/indexing/MultiLatticeSearch.h b/image_analysis/indexing/MultiLatticeSearch.h new file mode 100644 index 00000000..09eddabf --- /dev/null +++ b/image_analysis/indexing/MultiLatticeSearch.h @@ -0,0 +1,33 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "../common/CrystalLattice.h" +#include "../common/Coord.h" +#include "../lattice_search/LatticeSearch.h" + +struct MultiLatticeSearchResult { + CrystalLattice input_lattice; // lattice as it came in (kept for debugging only) + CrystalLattice output_lattice; // R * reference - a proper rotation of the first lattice + Coord rotation_vector; // Rodrigues vector: axis * angle, magnitude == angle [rad] +}; + +// Input: lattices assumed to describe the same crystal (same unit cell), differing +// only by orientation, and describing non-overlapping sets of spots. +// The first lattice is the reference. A subsequent lattice is kept only if its unit +// cell matches the reference (within tolerances); for each kept lattice we find the +// proper rotation R mapping reference -> lattice and store output_lattice = R * reference. +std::vector MultiLatticeSearch(const std::vector &lattices, + float dist_tolerance = 0.03f, + float angle_tolerance_deg = 3.0f); + +// Same as above, but first runs LatticeSearch on lattice 0 to obtain a +// Niggli-reduced / conventional reference, and reduces every candidate the +// same way before comparing/rotating. Rotations are returned in the reduced +// reference setting. +std::vector MultiLatticeSearchReduced(const std::vector &lattices, + float dist_tolerance = 0.03f, + float angle_tolerance_deg = 3.0f); \ No newline at end of file diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 951350cd..f4868b12 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -70,6 +70,7 @@ ADD_EXECUTABLE(jfjoch_test MergeScaleTest.cpp UnitCellTest.cpp CCTest.cpp + MultiLatticeSearchTest.cpp ) target_link_libraries(jfjoch_test Catch2WithMain JFJochBroker JFJochReceiver JFJochReader JFJochWriter diff --git a/tests/MultiLatticeSearchTest.cpp b/tests/MultiLatticeSearchTest.cpp new file mode 100644 index 00000000..b40cc537 --- /dev/null +++ b/tests/MultiLatticeSearchTest.cpp @@ -0,0 +1,68 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include "../image_analysis/indexing/MultiLatticeSearch.h" + +namespace { + // Rotation-invariant comparison via Gram matrix G = L^T L + void check_same_cell(const CrystalLattice &a, const CrystalLattice &b, double margin) { + const Coord av[3] = {a.Vec0(), a.Vec1(), a.Vec2()}; + const Coord bv[3] = {b.Vec0(), b.Vec1(), b.Vec2()}; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + CHECK((av[i] * av[j]) == Catch::Approx(bv[i] * bv[j]).margin(margin)); + } +} + +TEST_CASE("MultiLatticeSearch_RecoversRotation") { + CrystalLattice reference(40, 50, 80, 90, 95, 90); + + // Known rotation: 0.4 rad about a tilted axis + const Coord axis = Coord(0.3f, 0.9f, 0.1f).Normalize(); + const float angle = 0.4f; + const RotMatrix R(angle, axis); + + const CrystalLattice rotated(R * reference.Vec0(), + R * reference.Vec1(), + R * reference.Vec2()); + + auto result = MultiLatticeSearch({reference, rotated}); + + REQUIRE(result.size() == 2); + + // First entry: identity + CHECK(result[0].rotation_vector.Length() == Catch::Approx(0.0).margin(1e-6)); + check_same_cell(result[0].output_lattice, reference, 1e-3); + + // Second entry: angle and axis recovered + CHECK(result[1].rotation_vector.Length() == Catch::Approx(angle).margin(1e-4)); + const Coord recovered_axis = result[1].rotation_vector.Normalize(); + CHECK(recovered_axis.x == Catch::Approx(axis.x).margin(1e-3)); + CHECK(recovered_axis.y == Catch::Approx(axis.y).margin(1e-3)); + CHECK(recovered_axis.z == Catch::Approx(axis.z).margin(1e-3)); + + // output_lattice is a proper rotation of the reference => same metric + check_same_cell(result[1].output_lattice, reference, 1e-2); + + // and output equals the rotated input (same orientation) + check_same_cell(result[1].output_lattice, rotated, 1e-2); +} + +TEST_CASE("MultiLatticeSearch_SkipsDifferentCell") { + CrystalLattice reference(40, 50, 80, 90, 90, 90); + CrystalLattice other_cell(45, 50, 80, 90, 90, 90); // a differs by 5 A + + + auto result = MultiLatticeSearch({reference, other_cell}); + + // Only the reference survives + REQUIRE(result.size() == 1); + CHECK(result[0].rotation_vector.Length() == Catch::Approx(0.0).margin(1e-6)); +} + +TEST_CASE("MultiLatticeSearch_Empty") { + auto result = MultiLatticeSearch({}); + CHECK(result.empty()); +} \ No newline at end of file