diff --git a/code/include/matrix_utility.hpp b/code/include/matrix_utility.hpp new file mode 100644 index 0000000..f18aee6 --- /dev/null +++ b/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 + +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 diff --git a/code/include/projective_math.hpp b/code/include/projective_math.hpp new file mode 100644 index 0000000..2c8f297 --- /dev/null +++ b/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 +#include +#include +#include +#include + +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: = -sum_i p_i*q_i + p_last * q_last +// Distance: arcosh() 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(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 diff --git a/code/tests/CMakeLists.txt b/code/tests/CMakeLists.txt index f8e25d7..d379268 100644 --- a/code/tests/CMakeLists.txt +++ b/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 diff --git a/code/tests/test_matrix_utility.cpp b/code/tests/test_matrix_utility.cpp new file mode 100644 index 0000000..c300c8f --- /dev/null +++ b/code/tests/test_matrix_utility.cpp @@ -0,0 +1,40 @@ +// Port of de.varylab.discreteconformal.math.MatrixUtilityTest (Java/JUnit). + +#include "matrix_utility.hpp" + +#include +#include + +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; + } + } +} diff --git a/code/tests/test_surface_curve_utility.cpp b/code/tests/test_surface_curve_utility.cpp new file mode 100644 index 0000000..85db577 --- /dev/null +++ b/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 +#include + +using Eigen::VectorXd; +using conformallab::isOnSegment; +using conformallab::getPointOnCorrespondingSegment; +using conformallab::dehomogenize; + +// Helper: make VectorXd from initializer list. +static VectorXd V(std::initializer_list 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; + } +}