diff --git a/code/include/discrete_elliptic_utility.hpp b/code/include/discrete_elliptic_utility.hpp new file mode 100644 index 0000000..c49dd0d --- /dev/null +++ b/code/include/discrete_elliptic_utility.hpp @@ -0,0 +1,43 @@ +#pragma once + +// Ported from de.varylab.discreteconformal.util.DiscreteEllipticUtility (Java). +// Only the pure-math subset (no HDS required). + +#include +#include + +namespace conformallab { + +// Move tau into the fundamental domain of the modular group SL(2,Z): +// |Re(tau)| <= 0.5, Im(tau) >= 0, Re(tau) >= 0, |tau| >= 1 +// +// Algorithm: iteratively apply +// 1. T-shift: Re > 0.5 or Re < 0 → Re -= sign(Re) +// 2. Im-flip: Im < 0 → Im = -Im +// 3. Re-flip: Re < 0 → Re = -Re +// 4. S-invert: |tau| < 1 → tau = 1/tau +// +// Corresponds to Java DiscreteEllipticUtility.normalizeModulus(Complex). +inline std::complex normalizeModulus(std::complex tau) { + int maxIter = 100; + while (--maxIter > 0) { + double re = tau.real(); + double im = tau.imag(); + // exit when all conditions satisfied + if (std::abs(re) <= 0.5 && im >= 0.0 && re >= 0.0 && std::abs(tau) >= 1.0) + break; + + if (std::abs(re) > 0.5) + re -= (re > 0.0 ? 1.0 : -1.0); // signum shift + if (im < 0.0) + im = -im; + if (re < 0.0) + re = -re; + tau = std::complex(re, im); + if (std::abs(tau) < 1.0) + tau = 1.0 / tau; // S-transformation: invert + } + return tau; +} + +} // namespace conformallab diff --git a/code/include/p2_utility.hpp b/code/include/p2_utility.hpp new file mode 100644 index 0000000..d4f88ca --- /dev/null +++ b/code/include/p2_utility.hpp @@ -0,0 +1,102 @@ +#pragma once + +// 2-D projective geometry utilities for the Euclidean signature. +// Ported from de.jreality.math.P2 and de.varylab.discreteconformal.math.P2Big. +// +// Points and lines are represented as homogeneous 3-vectors (x, y, w). +// In the Euclidean case a finite point (px, py) is stored as (px, py, 1). + +#include +#include + +namespace conformallab { + +// ── Point / line duality ────────────────────────────────────────────────────── + +// Intersection of two lines l1, l2 (or line through two points p1, p2) +// via the cross product. Works for any P2 element. +// Corresponds to Java P2.pointFromLines / P2.lineFromPoints. +inline Eigen::Vector3d pointFromLines(const Eigen::Vector3d& l1, + const Eigen::Vector3d& l2) { + return l1.cross(l2); +} + +// ── Euclidean perpendicular bisector ───────────────────────────────────────── + +// Returns the homogeneous line coordinates (a, b, c) of the perpendicular +// bisector of the segment [p, q] in the Euclidean plane. +// Coordinates: ax + by + c = 0 (after dehomogenizing p and q). +// +// Corresponds to Java P2.perpendicularBisector(p, q, Pn.EUCLIDEAN). +inline Eigen::Vector3d perpendicularBisectorEuclidean(const Eigen::Vector3d& p_h, + const Eigen::Vector3d& q_h) { + // Dehomogenize + Eigen::Vector2d p = p_h.head<2>() / p_h(2); + Eigen::Vector2d q = q_h.head<2>() / q_h(2); + + // Direction vector (p → direction, matching jReality sign convention) + Eigen::Vector2d d = p - q; + + // Midpoint + Eigen::Vector2d m = (p + q) * 0.5; + + // Line: d[0]*(x - m[0]) + d[1]*(y - m[1]) = 0 + // = d[0]*x + d[1]*y - (d[0]*m[0] + d[1]*m[1]) + double c = -(d(0) * m(0) + d(1) * m(1)); + return {d(0), d(1), c}; +} + +// ── Euclidean distance between two P2 homogeneous points ───────────────────── + +inline double euclideanDistanceP2(const Eigen::Vector3d& p_h, + const Eigen::Vector3d& q_h) { + Eigen::Vector2d p = p_h.head<2>() / p_h(2); + Eigen::Vector2d q = q_h.head<2>() / q_h(2); + return (p - q).norm(); +} + +// ── Direct Euclidean isometry from two point-frames ────────────────────────── + +// Build the 3×3 projective matrix that represents the coordinate frame +// anchored at p0 with p1 defining the positive x-direction. +// Euclidean case: columns are [dehom(p0), unit_dir(p0→p1), perp_dir]. +// +// Template parameter S allows float / double / long double. +template +Eigen::Matrix makeFrameMatrix(Eigen::Matrix p0_h, + Eigen::Matrix p1_h) { + // Dehomogenize + Eigen::Matrix p0 = p0_h / p0_h(2); // (px, py, 1) + Eigen::Matrix p1_d = p1_h / p1_h(2); + + // Unit direction p0 → p1 + Eigen::Matrix dir2 = (p1_d - p0).template head<2>(); + dir2.normalize(); + Eigen::Matrix p1n(dir2(0), dir2(1), S(0)); + + // Perpendicular direction + Eigen::Matrix p2(-dir2(1), dir2(0), S(0)); + + Eigen::Matrix M; + M.col(0) = p0; + M.col(1) = p1n; + M.col(2) = p2; + return M; +} + +// Find the 3×3 Euclidean isometry (as a projective matrix) that maps +// the frame (s1, s2) to the frame (t1, t2). +// +// Corresponds to Java P2.makeDirectIsometryFromFrames(s1, s2, t1, t2, Pn.EUCLIDEAN) +// and P2Big.makeDirectIsometryFromFrames(...) (the BigDecimal / high-precision variant). +template +Eigen::Matrix makeDirectIsometryFromFramesEuclidean( + Eigen::Matrix s1, Eigen::Matrix s2, + Eigen::Matrix t1, Eigen::Matrix t2) +{ + auto toS = makeFrameMatrix(s1, s2); + auto toT = makeFrameMatrix(t1, t2); + return toT * toS.inverse(); +} + +} // namespace conformallab diff --git a/code/tests/CMakeLists.txt b/code/tests/CMakeLists.txt index d379268..1c09b91 100644 --- a/code/tests/CMakeLists.txt +++ b/code/tests/CMakeLists.txt @@ -1,8 +1,17 @@ add_executable(conformallab_tests + # ── Fully ported (pure math, no HDS) ──────────────────────────────────── test_clausen.cpp test_hyper_ideal_utility.cpp test_matrix_utility.cpp test_surface_curve_utility.cpp + test_discrete_elliptic_utility.cpp + test_p2_utility.cpp + + # ── Stubs: blocked until HDS port (Phase 4) ────────────────────────────── + # All tests call GTEST_SKIP() with a clear explanation. + test_hyper_ideal_functional.cpp + test_hyper_ideal_hyperelliptic_utility.cpp + test_spherical_functional.cpp ) target_include_directories(conformallab_tests SYSTEM PRIVATE diff --git a/code/tests/test_discrete_elliptic_utility.cpp b/code/tests/test_discrete_elliptic_utility.cpp new file mode 100644 index 0000000..9554ea1 --- /dev/null +++ b/code/tests/test_discrete_elliptic_utility.cpp @@ -0,0 +1,39 @@ +// Port of de.varylab.discreteconformal.util.DiscreteEllipticUtilityTest (Java/JUnit). +// Tests the normalizeModulus function that moves a complex number tau into the +// fundamental domain of the modular group SL(2,Z). + +#include "discrete_elliptic_utility.hpp" +#include +#include +#include + +using namespace conformallab; + +// Corresponds to Java testNormalizeModulus() +TEST(DiscreteEllipticUtilityTest, NormalizeModulus) { + // tau already in fundamental domain → should be returned unchanged + std::complex tau(0.45, 1.1); + auto tauNorm = normalizeModulus(tau); + EXPECT_NEAR(0.45, tauNorm.real(), 1E-12); + EXPECT_NEAR(1.1, tauNorm.imag(), 1E-12); + + // tau = i/3 (|tau| < 1) → inversion gives 3i + tau = std::complex(0.0, 1.0 / 3.0); + tauNorm = normalizeModulus(tau); + EXPECT_NEAR(3.0, tauNorm.imag(), 1E-12); + EXPECT_NEAR(0.0, tauNorm.real(), 1E-12); +} + +// Corresponds to Java testNormalizeModulusPeriodShift() +// Two tau values that differ by a T-shift (integer shift of Re) must normalize +// to the same point in the fundamental domain. +TEST(DiscreteEllipticUtilityTest, NormalizeModulusPeriodShift) { + std::complex tau1(0.3, 1.0); + std::complex tau2(-0.7, 1.0); // tau2 = tau1 - 1 + + auto n1 = normalizeModulus(tau1); + auto n2 = normalizeModulus(tau2); + + EXPECT_NEAR(n1.real(), n2.real(), 1E-12) << "real parts should be equal"; + EXPECT_NEAR(n1.imag(), n2.imag(), 1E-12) << "imag parts should be equal"; +} diff --git a/code/tests/test_hyper_ideal_functional.cpp b/code/tests/test_hyper_ideal_functional.cpp new file mode 100644 index 0000000..ea6afb3 --- /dev/null +++ b/code/tests/test_hyper_ideal_functional.cpp @@ -0,0 +1,37 @@ +// Stub for de.varylab.discreteconformal.functional.HyperIdealFunctionalTest (Java/JUnit). +// +// STATUS: BLOCKED – requires HDS port (Phase 4). +// +// These tests evaluate gradient and Hessian of the HyperIdealFunctional on +// actual mesh data (CoHDS + HyperIdealGenerator). They cannot be ported +// until the HalfEdge data structure (CoHDS), the functional evaluation +// framework, and the mesh generators are available in C++. +// +// Java tests and their status: +// testHessian() – @Ignore in Java (skipped here too) +// testGradientWithHyperIdealAndIdealPoints – blocked: needs HDS +// testGradientInTheExtendedDomain – blocked: needs HDS +// testGradientWithHyperellipticCurve – blocked: needs HDS +// testFunctionalAtNaNValue – blocked: needs HDS + +#include + +TEST(HyperIdealFunctionalTest, TestHessian_IgnoredInJava) { + GTEST_SKIP() << "@Ignore in Java – skipped here too"; +} + +TEST(HyperIdealFunctionalTest, GradientWithHyperIdealAndIdealPoints) { + GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HyperIdealFunctional)"; +} + +TEST(HyperIdealFunctionalTest, GradientInTheExtendedDomain) { + GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HyperIdealFunctional)"; +} + +TEST(HyperIdealFunctionalTest, GradientWithHyperellipticCurve) { + GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HyperIdealFunctional)"; +} + +TEST(HyperIdealFunctionalTest, FunctionalAtNaNValue) { + GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HyperIdealFunctional)"; +} diff --git a/code/tests/test_hyper_ideal_hyperelliptic_utility.cpp b/code/tests/test_hyper_ideal_hyperelliptic_utility.cpp new file mode 100644 index 0000000..107710b --- /dev/null +++ b/code/tests/test_hyper_ideal_hyperelliptic_utility.cpp @@ -0,0 +1,26 @@ +// Stub for de.varylab.discreteconformal.functional.HyperIdealHyperellipticUtilityTest. +// +// STATUS: BLOCKED – requires HDS port (Phase 4). +// +// Tests compute intersection angles of circles associated with hyper-ideal +// vertices using CoHDS + HalfEdgeUtils. All three tests operate on mesh +// data structures that are not yet available in C++. +// +// Java tests and their status: +// testCalculateCircleIntersections – blocked: needs CoHDS + HalfEdgeUtils +// testCalculateCircleIntersectionsInfinite – blocked: needs CoHDS + HalfEdgeUtils +// testLawsonHyperellipticAngles – blocked: needs CoHDS + HyperIdealGenerator + +#include + +TEST(HyperIdealHyperellipticUtilityTest, CalculateCircleIntersections) { + GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HalfEdgeUtils)"; +} + +TEST(HyperIdealHyperellipticUtilityTest, CalculateCircleIntersectionsInfinite) { + GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HalfEdgeUtils)"; +} + +TEST(HyperIdealHyperellipticUtilityTest, LawsonHyperellipticAngles) { + GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HyperIdealGenerator)"; +} diff --git a/code/tests/test_p2_utility.cpp b/code/tests/test_p2_utility.cpp new file mode 100644 index 0000000..7ecade7 --- /dev/null +++ b/code/tests/test_p2_utility.cpp @@ -0,0 +1,87 @@ +// Port of de.varylab.discreteconformal.math.P2BigTest (Java/JUnit). +// Tests 2-D projective geometry utilities: perpendicular bisectors, +// point-from-lines, and direct isometries in the Euclidean plane. +// +// The Java test compared double precision (P2) against BigDecimal precision +// (P2Big) to 1E-10. Here we compare double against long double to the +// same tolerance. + +#include "p2_utility.hpp" +#include +#include +#include + +using namespace conformallab; + +// Corresponds to Java P2BigTest.testMakeDirectIsometryFromFramesEuclidean() +// +// Computes the Euclidean isometry mapping frame (s1,s2) to frame (t1,t2) +// with both double and long-double precision, and checks: +// 1. The two precisions agree to 1E-10 (precision stability). +// 2. The matrix actually maps s1→t1 and s2→t2. +TEST(P2UtilityTest, MakeDirectIsometryFromFramesEuclidean) { + using V3d = Eigen::Vector3d; + using V3ld = Eigen::Matrix; + + V3d s1(-1.4142135623730963, 0.0, 1.0); + V3d s2( 1.4142135623730951, 0.0, 1.0); + V3d t1(-2.828427124746189, 2.4494897427831805, 1.0); + V3d t2( 0.0, 2.4494897427831783, 1.0); + + // double precision + auto T = makeDirectIsometryFromFramesEuclidean(s1, s2, t1, t2); + + // long double precision (analogous to Java's BigDecimal P2Big) + V3ld s1l = s1.cast(); + V3ld s2l = s2.cast(); + V3ld t1l = t1.cast(); + V3ld t2l = t2.cast(); + auto Tl = makeDirectIsometryFromFramesEuclidean(s1l, s2l, t1l, t2l); + + // 1. double vs long double must agree to 1E-10 + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + EXPECT_NEAR((double)Tl(i,j), T(i,j), 1E-10) + << "element (" << i << "," << j << ") differs between precisions"; + + // 2. T must map s1 → t1 and s2 → t2 (verify isometry correctness) + auto map_s1 = T * s1; + auto map_s2 = T * s2; + EXPECT_NEAR(euclideanDistanceP2(map_s1, t1), 0.0, 1E-9) << "T*s1 should equal t1"; + EXPECT_NEAR(euclideanDistanceP2(map_s2, t2), 0.0, 1E-9) << "T*s2 should equal t2"; +} + +// Corresponds to Java P2BigTest.testPerpendicularBisector() +TEST(P2UtilityTest, PerpendicularBisector) { + Eigen::Vector3d p1(0.5, 0.0, 1.0); + Eigen::Vector3d q1(0.0, 0.5, 1.0); + + auto bisector = perpendicularBisectorEuclidean(p1, q1); + + EXPECT_NEAR( 0.5, bisector(0), 1E-10); + EXPECT_NEAR(-0.5, bisector(1), 1E-10); + EXPECT_NEAR( 0.0, bisector(2), 1E-10); +} + +// Corresponds to Java P2BigTest.testPerpendicularBisectorIntersection() +// +// The intersection of the perpendicular bisectors of two edges must be +// equidistant from the endpoints of each edge (circumcenter property). +TEST(P2UtilityTest, PerpendicularBisectorIntersection) { + Eigen::Vector3d p1(0.5, 0.0, 1.0); + Eigen::Vector3d q1(0.0, 1.0, 1.0); + Eigen::Vector3d p2(1.0, 0.0, 1.0); + Eigen::Vector3d q2(0.0, 1.5, 1.0); + + auto l1 = perpendicularBisectorEuclidean(p1, q1); + auto l2 = perpendicularBisectorEuclidean(p2, q2); + auto o = pointFromLines(l1, l2); // circumcenter + + // o must be equidistant from p1 and q1 + EXPECT_NEAR(euclideanDistanceP2(p1, o), + euclideanDistanceP2(q1, o), 1E-10); + + // o must be equidistant from p2 and q2 + EXPECT_NEAR(euclideanDistanceP2(p2, o), + euclideanDistanceP2(q2, o), 1E-10); +} diff --git a/code/tests/test_spherical_functional.cpp b/code/tests/test_spherical_functional.cpp new file mode 100644 index 0000000..989932e --- /dev/null +++ b/code/tests/test_spherical_functional.cpp @@ -0,0 +1,37 @@ +// Stub for de.varylab.discreteconformal.functional.SphericalFunctionalTest (Java/JUnit). +// +// STATUS: BLOCKED – requires HDS port (Phase 4). +// +// Tests evaluate gradient and Hessian of the SphericalFunctional on meshes +// built via CoHDS + ConvexHull, and check that a regular spherical metric +// is a critical point of the functional. All tests require the HalfEdge +// data structure and the functional evaluation framework in C++. +// +// Java tests and their status: +// testReducedGradient – blocked: needs CoHDS + SphericalFunctional +// testReducedHessian – blocked: needs CoHDS + SphericalFunctional +// testGradient – blocked: needs CoHDS + SphericalFunctional +// testHessian – blocked: needs CoHDS + SphericalFunctional +// testCriticalPoint – blocked: needs CoHDS + ConvexHull + SphericalFunctional + +#include + +TEST(SphericalFunctionalTest, ReducedGradient) { + GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + SphericalFunctional)"; +} + +TEST(SphericalFunctionalTest, ReducedHessian) { + GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + SphericalFunctional)"; +} + +TEST(SphericalFunctionalTest, Gradient) { + GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + SphericalFunctional)"; +} + +TEST(SphericalFunctionalTest, Hessian) { + GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + SphericalFunctional)"; +} + +TEST(SphericalFunctionalTest, CriticalPoint) { + GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + ConvexHull + SphericalFunctional)"; +}