Jungfraujoch/tests/LatticeReductionTest.cpp

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

#include <catch2/catch_all.hpp>

#include <Eigen/Dense>

#include "../image_analysis/geom_refinement/LatticeReduction.h"


namespace {
    Eigen::Vector3d to_eigen(const Coord& v) {
        return {v[0], v[1], v[2]};
    }

    double angle_rad(const Eigen::Vector3d& a, const Eigen::Vector3d& b) {
        const double na = a.norm();
        const double nb = b.norm();
        if (na == 0.0 || nb == 0.0)
            return 0.0;
        double c = a.dot(b) / (na * nb);
        c = std::max(-1.0, std::min(1.0, c));
        return std::acos(c);
    }

    // Compare two lattices up to a global rotation: compare Gram matrices G = L^T L (rotation-invariant).
    Eigen::Matrix3d gram(const CrystalLattice& latt) {
        const Eigen::Vector3d A = to_eigen(latt.Vec0());
        const Eigen::Vector3d B = to_eigen(latt.Vec1());
        const Eigen::Vector3d C = to_eigen(latt.Vec2());
        Eigen::Matrix3d G;
        G(0,0) = A.dot(A); G(0,1) = A.dot(B); G(0,2) = A.dot(C);
        G(1,0) = B.dot(A); G(1,1) = B.dot(B); G(1,2) = B.dot(C);
        G(2,0) = C.dot(A); G(2,1) = C.dot(B); G(2,2) = C.dot(C);
        return G;
    }

    void check_gram_close(const CrystalLattice& a,
                          const CrystalLattice& b,
                          double abs_eps,
                          double rel_eps) {
        const Eigen::Matrix3d Ga = gram(a);
        const Eigen::Matrix3d Gb = gram(b);

        for (int r = 0; r < 3; ++r) {
            for (int c = 0; c < 3; ++c) {
                const double va = Ga(r, c);
                const double vb = Gb(r, c);

                // Scale for relative error; avoid blowing up around zero.
                const double scale = std::max({1.0, std::abs(va), std::abs(vb)});
                const double tol = std::max(abs_eps, rel_eps * scale);

                INFO("G(" << r << "," << c << ") va=" << va << " vb=" << vb
                          << " scale=" << scale << " tol=" << tol);

                CHECK(va == Catch::Approx(vb).margin(tol));
            }
        }
    }
} // namespace


TEST_CASE("LatticeToRodrigues") {
    double rod[3];
    double lengths[3];

    CrystalLattice latt_i(40,50,80,90,90,90);

    LatticeToRodriguesAndLengths_GS(latt_i, rod, lengths);
    CHECK(lengths[0] == Catch::Approx(40.0));
    CHECK(lengths[1] == Catch::Approx(50.0));
    CHECK(lengths[2] == Catch::Approx(80.0));
    CHECK(fabs(rod[0]) < 0.001);
    CHECK(fabs(rod[1]) < 0.001);
    CHECK(fabs(rod[2]) < 0.001);

    auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, M_PI/2, M_PI/2);

    CHECK(latt_o.Vec0().Length() == Catch::Approx(40.0));
    CHECK(latt_o.Vec1().Length() == Catch::Approx(50.0));
    CHECK(latt_o.Vec2().Length() == Catch::Approx(80.0));
}


TEST_CASE("LatticeToRodrigues_irregular") {
    double rod[3];
    double lengths[3];

    CrystalLattice latt_i(Coord(40,0,0),
                          Coord(0, 50 / sqrt(2), -50 / sqrt(2)),
                          Coord(0, 80 / sqrt(2), 80 / sqrt(2)));

    LatticeToRodriguesAndLengths_GS(latt_i, rod, lengths);
    CHECK(lengths[0] == Catch::Approx(40.0));
    CHECK(lengths[1] == Catch::Approx(50.0));
    CHECK(lengths[2] == Catch::Approx(80.0));

    auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, M_PI/2, M_PI/2);

    CHECK(latt_o.Vec0().Length() == Catch::Approx(40.0));
    CHECK(latt_o.Vec1().Length() == Catch::Approx(50.0));
    CHECK(latt_o.Vec2().Length() == Catch::Approx(80.0));
}

TEST_CASE("LatticeToRodrigues_Hex") {
    double rod[3];
    double lengths[3];

    Coord a = Coord(40,0,0);
    Coord b = Coord(40 / 2, 40 * sqrt(3)/ 2.0, 0);
    Coord c = Coord(0, 0, 70);

    RotMatrix R(1.0, Coord(0,1,1));
    CrystalLattice latt_i(R*a,R*b,R*c);

    LatticeToRodriguesAndLengths_Hex(latt_i, rod, lengths);
    CHECK(lengths[0] == Catch::Approx(40.0));
    CHECK(lengths[1] == Catch::Approx(40.0));
    CHECK(lengths[2] == Catch::Approx(70.0));

    auto latt_o = AngleAxisAndCellToLattice(rod, lengths,M_PI / 2.0, M_PI / 2.0, 2.0 * M_PI / 3.0);
    auto uc_o = latt_o.GetUnitCell();
    CHECK(uc_o.a == Catch::Approx(40.0));
    CHECK(uc_o.b == Catch::Approx(40.0));
    CHECK(uc_o.c == Catch::Approx(70.0));
    CHECK(uc_o.alpha == Catch::Approx(90.0));
    CHECK(uc_o.beta == Catch::Approx(90.0));
    CHECK(uc_o.gamma == Catch::Approx(120.0));
}


TEST_CASE("LatticeReduction Lattice param roundtrip (GS) preserves Gram matrix") {
    // Non-orthogonal, irregular basis, but still a valid lattice
    CrystalLattice latt_i(Coord(40,  1,  2),
                          Coord( 3, 50, -4),
                          Coord(-5,  6, 80));

    double rod[3]{}, lengths[3]{};
    LatticeToRodriguesAndLengths_GS(latt_i, rod, lengths);
    auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, M_PI/2, M_PI/2);

    // This parametrization only keeps "lengths + rotation", i.e. it cannot reproduce shear.
    // So we *do not* compare Gram matrices here.
    // Instead, sanity-check: reconstructed vectors have the requested lengths.
    CHECK(latt_o.Vec0().Length() == Catch::Approx(lengths[0]).margin(1e-9));
    CHECK(latt_o.Vec1().Length() == Catch::Approx(lengths[1]).margin(1e-9));
    CHECK(latt_o.Vec2().Length() == Catch::Approx(lengths[2]).margin(1e-9));
}

TEST_CASE("LatticeReduction Lattice param roundtrip (Hex) preserves unit cell") {
    Coord a = Coord(40, 0, 0);
    Coord b = Coord(40 / 2.0, 40 * std::sqrt(3) / 2.0, 0);
    Coord c = Coord(0, 0, 70);

    // Apply an arbitrary rotation to ensure the rod extraction is meaningful
    RotMatrix R(1.0, Coord(0, 1, 1));
    CrystalLattice latt_i(R * a, R * b, R * c);

    double rod[3]{}, ac[3]{};
    LatticeToRodriguesAndLengths_Hex(latt_i, rod, ac);
    auto latt_o = AngleAxisAndCellToLattice(rod, ac, M_PI/2, M_PI/2, 2*M_PI/3);

    auto uc_o = latt_o.GetUnitCell();
    CHECK(uc_o.a == Catch::Approx(40.0).margin(1e-6));
    CHECK(uc_o.b == Catch::Approx(40.0).margin(1e-6));
    CHECK(uc_o.c == Catch::Approx(70.0).margin(1e-6));
    CHECK(uc_o.alpha == Catch::Approx(90.0).margin(1e-6));
    CHECK(uc_o.beta  == Catch::Approx(90.0).margin(1e-6));
    CHECK(uc_o.gamma == Catch::Approx(120.0).margin(1e-6));
}

TEST_CASE("LatticeReduction Monoclinic param roundtrip") {
    struct Case { double beta_deg; };
    const std::vector<Case> cases = {
        {60.0}, {75.0}, {115.0}, {130.0}
    };

    for (const auto& cs : cases) {
        INFO("beta_deg=" << cs.beta_deg);

        // Start from a clean monoclinic cell in its conventional setting (unique axis b).
        CrystalLattice latt0(50, 60, 70, 90, cs.beta_deg, 90);

        // Now apply a TRUE global rotation: rotate each basis vector (left-multiply).
        RotMatrix R(0.7, Coord(0.3, 0.9, 0.1));
        CrystalLattice latt_i(R * latt0.Vec0(),
                              R * latt0.Vec1(),
                              R * latt0.Vec2());

        double rod[3]{}, lengths[3]{}, beta_rad = 0.0;
        LatticeToRodriguesLengthsBeta_Mono(latt_i, rod, lengths, beta_rad);

        // Basic sanity
        CHECK(lengths[0] == Catch::Approx(50.0).margin(1e-6));
        CHECK(lengths[1] == Catch::Approx(60.0).margin(1e-6));
        CHECK(lengths[2] == Catch::Approx(70.0).margin(1e-6));
        CHECK(beta_rad * 180.0 / M_PI == Catch::Approx(cs.beta_deg).margin(1e-6));

        auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, beta_rad, M_PI/2);

        // Rotation-invariant check: Gram matrices match.
        check_gram_close(latt_i, latt_o, /*abs_eps=*/5e-4, /*rel_eps=*/1e-10);

        // Also check the unit-cell angles we expect for monoclinic(unique b).
        auto uc_o = latt_o.GetUnitCell();
        CHECK(uc_o.alpha == Catch::Approx(90.0).margin(1e-4));
        CHECK(uc_o.gamma == Catch::Approx(90.0).margin(1e-4));
        CHECK(uc_o.beta  == Catch::Approx(cs.beta_deg).margin(1e-4));
    }
}

TEST_CASE("LatticeReduction monoclinic") {
    // This isolates only the beta definition, independent of other choices.
    CrystalLattice latt_i(50, 60, 70, 90, 130, 90);

    const Eigen::Vector3d A = to_eigen(latt_i.Vec0());
    const Eigen::Vector3d C = to_eigen(latt_i.Vec2());
    const double beta_geom = angle_rad(A, C);

    double rod[3]{}, lengths[3]{}, beta_rad = 0.0;
    LatticeToRodriguesLengthsBeta_Mono(latt_i, rod, lengths, beta_rad);

    CHECK(beta_rad == Catch::Approx(beta_geom).margin(1e-12));
}