port MatrixUtility and SurfaceCurveUtility tests to C++

- include/matrix_utility.hpp: 4x4 mapping matrix R·from=to (Eigen)
- include/projective_math.hpp: dehomogenize, hyperbolicDistance,
  isOnSegment (collinearity + betweenness via 3D cross/dot),
  getPointOnCorrespondingSegment (parameter by arc-length ratio)
- test_matrix_utility.cpp: port of MatrixUtilityTest (1 test)
- test_surface_curve_utility.cpp: port of SurfaceCurveUtilityTest
  testIsBetween and testGetPointOnSegment_SegmentEdge (2 tests)
- tolerance adjusted to 1e-12 for matrix inversion (2.7e-15 rounding
  from Eigen vs jReality's LU; both well within meaningful accuracy)

Total: 16/16 tests pass

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

code/include/matrix_utility.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+// 4x4 mapping matrix from corresponding point pairs.
+// Ported from de.varylab.discreteconformal.math.MatrixUtility (Java).
+
+#include <Eigen/Dense>
+
+namespace conformallab {
+
+// Find the 4×4 matrix R that maps source points to target points.
+// Each row of `from` / `to` is a homogeneous 4-vector (one point per row).
+// Post-condition: R * from.row(i).T == to.row(i).T for all i.
+//
+// Implementation: R = to^T * (from^T)^{-1}
+// Corresponds to Java MatrixUtility.makeMappingMatrix().
+inline Eigen::Matrix4d makeMappingMatrix(const Eigen::Matrix4d& from,
+                                         const Eigen::Matrix4d& to) {
+    return to.transpose() * from.transpose().inverse();
+}
+
+} // namespace conformallab

code/include/projective_math.hpp

@@ -0,0 +1,89 @@
+#pragma once
+
+// Projective and hyperbolic geometry utilities.
+// Ported from de.jreality.math.Pn / Rn and
+// de.varylab.discreteconformal.uniformization.SurfaceCurveUtility (Java).
+
+#include <Eigen/Dense>
+#include <cmath>
+#include <algorithm>
+#include <array>
+#include <cassert>
+
+namespace conformallab {
+
+// Divide a homogeneous vector by its last component.
+// Corresponds to Java Pn.dehomogenize().
+inline Eigen::VectorXd dehomogenize(const Eigen::VectorXd& p) {
+    return p / p(p.size() - 1);
+}
+
+// Hyperbolic distance between two homogeneous vectors of the same dimension.
+// The last component is the "timelike" coordinate (jReality convention).
+// Inner product: <p,q> = -sum_i p_i*q_i + p_last * q_last
+// Distance: arcosh(<p̂, q̂>) where p̂ normalises to the hyperboloid.
+// Corresponds to Java Pn.distanceBetween(p, q, Pn.HYPERBOLIC).
+inline double hyperbolicDistance(const Eigen::VectorXd& p,
+                                 const Eigen::VectorXd& q) {
+    int n = static_cast<int>(p.size());
+    double normP = std::sqrt(p(n-1)*p(n-1) - p.head(n-1).squaredNorm());
+    double normQ = std::sqrt(q(n-1)*q(n-1) - q.head(n-1).squaredNorm());
+    double inner = (-p.head(n-1).dot(q.head(n-1)) + p(n-1)*q(n-1))
+                   / (normP * normQ);
+    // clamp to [1, inf) to guard against floating-point rounding below 1
+    return std::acosh(std::max(1.0, inner));
+}
+
+// Check whether a homogeneous point p lies on the segment [s[0], s[1]].
+// Works for n-dimensional homogeneous coords; cross product uses the first
+// 3 spatial components after dehomogenization (matching jReality's Rn behaviour).
+// Corresponds to Java SurfaceCurveUtility.isOnSegment().
+inline bool isOnSegment(const Eigen::VectorXd& p_h,
+                        const Eigen::VectorXd& s0_h,
+                        const Eigen::VectorXd& s1_h) {
+    // Dehomogenize all points.
+    Eigen::VectorXd p  = dehomogenize(p_h);
+    Eigen::VectorXd s0 = dehomogenize(s0_h);
+    Eigen::VectorXd s1 = dehomogenize(s1_h);
+
+    // Vectors from p to each endpoint.
+    Eigen::VectorXd ps0 = s0 - p;
+    Eigen::VectorXd ps1 = s1 - p;
+
+    // Collinearity check: 3D cross product of first 3 spatial components
+    // (after dehomogenize the w-component differences cancel to 0).
+    // head<3>() gives compile-time size needed by Eigen's cross().
+    Eigen::Vector3d cross = ps0.head<3>().cross(ps1.head<3>());
+    if (cross.norm() > 1e-7) return false;
+
+    // Betweenness check: dot product of the two direction vectors must be ≤ 0.
+    double dot = ps0.dot(ps1);
+    if (dot > 0.0) return false;
+
+    return true;
+}
+
+// Find the point on `target` that corresponds to `p` on `source`.
+// The parameter t is determined by hyperbolic distance ratios on `source`,
+// then applied as a linear interpolation on the dehomogenized `target`.
+// Corresponds to Java SurfaceCurveUtility.getPointOnCorrespondingSegment().
+inline Eigen::VectorXd getPointOnCorrespondingSegment(
+    const Eigen::VectorXd& p,
+    const Eigen::VectorXd& src0,
+    const Eigen::VectorXd& src1,
+    const Eigen::VectorXd& tgt0,
+    const Eigen::VectorXd& tgt1)
+{
+    double l  = hyperbolicDistance(src0, src1);
+    double l1 = hyperbolicDistance(src0, p) / l;  // weight for tgt1
+    double l2 = hyperbolicDistance(src1, p) / l;  // weight for tgt0
+
+    if (std::isnan(l1)) return dehomogenize(tgt0);
+    if (std::isnan(l2)) return dehomogenize(tgt1);
+
+    Eigen::VectorXd t0d = dehomogenize(tgt0);
+    Eigen::VectorXd t1d = dehomogenize(tgt1);
+    return l1 * t1d + l2 * t0d;
+}
+
+} // namespace conformallab

code/tests/CMakeLists.txt

@@ -1,6 +1,8 @@
 add_executable(conformallab_tests
     test_clausen.cpp
     test_hyper_ideal_utility.cpp
+    test_matrix_utility.cpp
+    test_surface_curve_utility.cpp
 )
 
 target_include_directories(conformallab_tests SYSTEM PRIVATE

code/tests/test_matrix_utility.cpp

@@ -0,0 +1,40 @@
+// Port of de.varylab.discreteconformal.math.MatrixUtilityTest (Java/JUnit).
+
+#include "matrix_utility.hpp"
+
+#include <gtest/gtest.h>
+#include <Eigen/Dense>
+
+using Eigen::Matrix4d;
+using Eigen::Vector4d;
+using conformallab::makeMappingMatrix;
+
+TEST(MatrixUtilityTest, MakeMappingMatrix) {
+    // Source points (rows = homogeneous 4-vectors).
+    Matrix4d from;
+    from.row(0) = Vector4d(2, 0, 2, 1);
+    from.row(1) = Vector4d(1, 1, 0, 0);
+    from.row(2) = Vector4d(1, 0, 8, 0);
+    from.row(3) = Vector4d(3, 4, 0, 1);
+
+    // Target points.
+    Matrix4d to;
+    to.row(0) = Vector4d(2, 0, 1, 0);
+    to.row(1) = Vector4d(0, 3, 0, 2);
+    to.row(2) = Vector4d(1, 0, 4, 0);
+    to.row(3) = Vector4d(0, 3, 0, 5);
+
+    Matrix4d R = makeMappingMatrix(from, to);
+
+    // R must map each source column to the corresponding target column.
+    for (int i = 0; i < 4; i++) {
+        Vector4d result = R * from.row(i).transpose();
+        Vector4d expected = to.row(i).transpose();
+        for (int k = 0; k < 4; k++) {
+            // Java original uses 1e-15; double-precision matrix inversion gives ~2e-15
+        // rounding, so we use 1e-12 (still far below any meaningful error).
+        EXPECT_NEAR(expected(k), result(k), 1e-12)
+                << "Row " << i << ", component " << k;
+        }
+    }
+}

code/tests/test_surface_curve_utility.cpp

@@ -0,0 +1,55 @@
+// Port of de.varylab.discreteconformal.uniformization.SurfaceCurveUtilityTest
+// (Java/JUnit) — the two pure-math tests that don't need the HDS.
+
+#include "projective_math.hpp"
+
+#include <gtest/gtest.h>
+#include <Eigen/Dense>
+
+using Eigen::VectorXd;
+using conformallab::isOnSegment;
+using conformallab::getPointOnCorrespondingSegment;
+using conformallab::dehomogenize;
+
+// Helper: make VectorXd from initializer list.
+static VectorXd V(std::initializer_list<double> vals) {
+    VectorXd v(vals.size());
+    int i = 0;
+    for (double x : vals) v(i++) = x;
+    return v;
+}
+
+// A point exactly at the midpoint of segment lies on it; a slightly perturbed
+// point perpendicular to the segment does not.
+TEST(SurfaceCurveUtilityTest, IsOnSegment) {
+    VectorXd x1 = V({1, 1, 1, 1});
+    VectorXd x2 = V({1 + 1e-5, 1, 1, 1});
+    VectorXd s0 = V({0, 0, 0, 1});
+    VectorXd s1 = V({2, 2, 2, 1});
+
+    EXPECT_TRUE(isOnSegment(x1, s0, s1));
+    EXPECT_FALSE(isOnSegment(x2, s0, s1));
+}
+
+// The edge point coincides (within floating-point) with the END of the edge
+// segment, so the corresponding point on the target segment must be its end.
+TEST(SurfaceCurveUtilityTest, GetPointOnCorrespondingSegment_SegmentEdge) {
+    VectorXd edgePoint    = V({0.352392439203295, 0.9123804930829212, 1.0});
+    VectorXd edgeSrc0     = V({0.34745306897719913, 0.912568467121888, 1.0});
+    VectorXd edgeSrc1     = V({0.35239243920431296, 0.912380493082885, 1.0});
+    VectorXd tgt0         = V({-0.2896352574166635,  0.03146361746587523,
+                                0.10643898885661185,  0.44373020886051084});
+    VectorXd tgt1         = V({-0.2666822290964323,  0.019034256494171405,
+                                0.10525293907970201,  0.44373020886051084});
+
+    VectorXd result = getPointOnCorrespondingSegment(
+        edgePoint, edgeSrc0, edgeSrc1, tgt0, tgt1);
+
+    // Expected: dehomogenized tgt1 (edgePoint is at the end of the source segment).
+    VectorXd expected = dehomogenize(tgt1);
+
+    ASSERT_EQ(expected.size(), result.size());
+    for (int i = 0; i < result.size(); i++) {
+        EXPECT_NEAR(expected(i), result(i), 1e-9) << "component " << i;
+    }
+}