Jungfraujoch/tests/MultiLatticeSearchTest.cpp

// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include <catch2/catch_all.hpp>

#include "../image_analysis/indexing/MultiLatticeSearch.h"

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() == 1);

    // Second entry: angle and axis recovered
    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));
}

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.empty());
}

TEST_CASE("MultiLatticeSearch_Empty") {
    auto result = MultiLatticeSearch({});
    CHECK(result.empty());
}

TEST_CASE("MultiLatticeSearch_EP") {
    // Real EP case
    CrystalLattice cell1(Coord(-13.2, -30.0, -29.6),
        Coord(70.1, -12.16, -18.6),
        Coord(7.6, -23.9, 44.8));
    CrystalLattice cell2(Coord(-13.2, -29.9, -29.5),
        Coord(-70.1, 12.15, 18.8),
        Coord(1.8, 45.4, -23.7)
    );

    auto result = MultiLatticeSearch({cell1, cell2}, 0.1, 3);

    REQUIRE(result.size() == 1);

    RotMatrix matrix(result[0].rotation_vector.Length(), result[0].rotation_vector.Normalize());
    auto cell2_rot = cell1.Multiply(matrix);

    CHECK((cell2_rot.Vec0() - cell2.Vec0()).Length() < cell2.Vec0().Length() / 100);
    CHECK((cell2_rot.Vec1() - cell2.Vec1()).Length() < cell2.Vec1().Length() / 100);
    CHECK((cell2_rot.Vec2() - cell2.Vec2()).Length() < cell2.Vec2().Length() / 100);
}