diff --git a/image_analysis/indexing/MultiLatticeSearch.cpp b/image_analysis/indexing/MultiLatticeSearch.cpp index 78681e55..07da2edb 100644 --- a/image_analysis/indexing/MultiLatticeSearch.cpp +++ b/image_analysis/indexing/MultiLatticeSearch.cpp @@ -3,6 +3,7 @@ #include "MultiLatticeSearch.h" +#include #include namespace { @@ -69,13 +70,9 @@ std::vector MultiLatticeSearch(const std::vector MultiLatticeSearch(const std::vector(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 - }); + std::cout << rotation_vector << std::endl; + std::cout << rotation_vector.Length() * 180.0 / M_PI << " deg\n" << std::endl; + + if (found) + continue; + ret.push_back({rotation_vector}); found = true; } } diff --git a/image_analysis/indexing/MultiLatticeSearch.h b/image_analysis/indexing/MultiLatticeSearch.h index e5f341ae..af369bb5 100644 --- a/image_analysis/indexing/MultiLatticeSearch.h +++ b/image_analysis/indexing/MultiLatticeSearch.h @@ -7,11 +7,8 @@ #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] }; @@ -19,7 +16,7 @@ struct MultiLatticeSearchResult { // 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. +// proper rotation R mapping reference -> lattice std::vector MultiLatticeSearch(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/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index 7da8584e..8c865ad4 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -1211,10 +1211,10 @@ TEST_CASE("JFJochReader_Spots_OldMasterFormat", "[HDF5][Full]") { std::vector spots; spots.push_back(SpotToSave{ .x = 1, .y = 2, .intensity = 376, - .ice_ring = true, - .indexed = true, .h = 11, .k = -3, .l = -5, - .dist_ewald_sphere = 0.1234f + .dist_ewald_sphere = 0.1234f, + .ice_ring = true, + .indexed = true }); spots.push_back(SpotToSave{ .x = 7, .y = -3, .intensity = 0.156f, @@ -1308,10 +1308,10 @@ TEST_CASE("JFJochReader_Spots_VDS", "[HDF5][Full]") { std::vector spots; spots.push_back(SpotToSave{ .x = 1, .y = 2, .intensity = 376, - .ice_ring = true, - .indexed = true, .h = 11, .k = -3, .l = -5, - .dist_ewald_sphere = 0.1234f + .dist_ewald_sphere = 0.1234f, + .ice_ring = true, + .indexed = true }); spots.push_back(SpotToSave{ .x = 7, .y = -3, .intensity = 0.156f, @@ -2127,10 +2127,10 @@ static std::vector MakeTestSpots(int i) { return { SpotToSave{ .x = 1, .y = 2, .intensity = 376, - .ice_ring = true, - .indexed = true, .h = 11, .k = -3, .l = -5, - .dist_ewald_sphere = 0.1234f + .dist_ewald_sphere = 0.1234f, + .ice_ring = true, + .indexed = true }, SpotToSave{ .x = 7, .y = static_cast(-3 - i), .intensity = 0.156f, diff --git a/tests/MultiLatticeSearchTest.cpp b/tests/MultiLatticeSearchTest.cpp index 6e5df735..4f94b868 100644 --- a/tests/MultiLatticeSearchTest.cpp +++ b/tests/MultiLatticeSearchTest.cpp @@ -5,17 +5,6 @@ #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); @@ -30,36 +19,24 @@ TEST_CASE("MultiLatticeSearch_RecoversRotation") { 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); + REQUIRE(result.size() == 1); // 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(result[0].rotation_vector.Length() == Catch::Approx(angle).margin(1e-4)); + const Coord recovered_axis = result[0].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)); + REQUIRE(result.empty()); } TEST_CASE("MultiLatticeSearch_Empty") { @@ -67,6 +44,8 @@ TEST_CASE("MultiLatticeSearch_Empty") { CHECK(result.empty()); } +#include + TEST_CASE("MultiLatticeSearch_EP") { // Real EP case CrystalLattice cell1(Coord(-13.2, -30.0, -29.6), @@ -77,7 +56,10 @@ TEST_CASE("MultiLatticeSearch_EP") { Coord(1.8, 45.4, -23.7) ); - auto result = MultiLatticeSearch({cell1, cell2}, - 0.1, 5); - REQUIRE(result.size() == 2); + auto result = MultiLatticeSearch({cell1, cell2}, 0.1, 3); + + + REQUIRE(result.size() == 1); + std::cout << result[0].rotation_vector << std::endl; + std::cout << result[0].rotation_vector.Length() * 180.0 / M_PI << std::endl; } \ No newline at end of file