Merge pull request 'feat(geometry): pn_geometry.hpp — Pn projective-metric substrate (jReality port)' (#31) from feat/pn-geometry-substrate into main

code/include/pn_geometry.hpp

@@ -0,0 +1,156 @@
+#pragma once
+// Copyright (c) 2024-2026 Tarik Moussa.
+// SPDX-License-Identifier: MIT
+
+// pn_geometry.hpp
+//
+// Projective-metric geometry substrate — faithful port of the
+// `de.jreality.math.Pn` surface that the Java algorithmic core
+// uses pervasively:
+//   HyperbolicLayout, SphericalLayout, FundamentalPolygon,
+//   KoebePolyhedron, quasi-isothermic utilities, hyperelliptic θ.
+//
+// Porting this ~30-function surface once unblocks every downstream
+// geometric phase (9c / 10c' / 10e / 11c) instead of re-deriving the
+// metric ad hoc per phase.
+//
+// ┌──────────────────────────────────────────────────────────────────┐
+// │  Metric signature convention (matches jReality exactly)          │
+// │                                                                  │
+// │  EUCLIDEAN  = 0   ⟨u,v⟩ = Σ_{i<last} uᵢvᵢ         (spatial)  │
+// │  ELLIPTIC   = +1  ⟨u,v⟩ = Σ_{all} uᵢvᵢ             (S^n)     │
+// │  HYPERBOLIC = -1  ⟨u,v⟩ = u_last·v_last − Σ_{i<last} uᵢvᵢ    │
+// │                           (Minkowski; last coord = timelike)    │
+// │                                                                  │
+// │  HYPERBOLIC uses the timelike-positive ("upper sheet") convention│
+// │  so a normalised hyperbolic point satisfies ⟨p,p⟩ = +1 and     │
+// │  ⟨p,q⟩ ≥ 1 → d(p,q) = acosh(⟨p,q⟩) directly.                 │
+// │  This is identical to projective_math.hpp::hyperbolicDistance,  │
+// │  against which pn_distance_between(..., HYPERBOLIC) is anchored. │
+// └──────────────────────────────────────────────────────────────────┘
+//
+// Vectors are Eigen column vectors (homogeneous, length = dim+1).
+
+#include <Eigen/Core>
+#include <cmath>
+#include <algorithm>
+
+namespace conformallab {
+
+/// Metric signatures matching `de.jreality.math.Pn`.
+enum PnMetric : int {
+    PN_EUCLIDEAN  =  0,
+    PN_ELLIPTIC   = +1,
+    PN_HYPERBOLIC = -1
+};
+
+// ── Inner product ─────────────────────────────────────────────────────────────
+
+/// Bilinear form ⟨u,v⟩ for the given metric.
+/// Spatial part = all but the last coordinate; the last coord is the
+/// homogeneous/timelike one.
+inline double pn_inner_product(const Eigen::VectorXd& u,
+                               const Eigen::VectorXd& v,
+                               int metric)
+{
+    const int n       = static_cast<int>(u.size());
+    const double spat = u.head(n - 1).dot(v.head(n - 1));
+    const double last = u(n - 1) * v(n - 1);
+    switch (metric) {
+        case PN_HYPERBOLIC: return  last - spat;
+        case PN_ELLIPTIC:   return  last + spat;
+        case PN_EUCLIDEAN:
+        default:            return  spat;
+    }
+}
+
+// ── Norm / scale ──────────────────────────────────────────────────────────────
+
+/// Metric norm √|⟨p,p⟩|.  Absolute value guards against tiny negative
+/// round-off in the Euclidean / hyperbolic degenerate cases.
+inline double pn_norm(const Eigen::VectorXd& p, int metric)
+{
+    return std::sqrt(std::abs(pn_inner_product(p, p, metric)));
+}
+
+/// Dehomogenise: divide by the last component (jReality `Pn.dehomogenize`).
+inline Eigen::VectorXd pn_dehomogenize(const Eigen::VectorXd& p)
+{
+    return p / p(p.size() - 1);
+}
+
+/// Scale `p` to metric norm `length` (jReality `Pn.setToLength`).
+/// Returns `p` unchanged when its norm is numerically zero.
+inline Eigen::VectorXd pn_set_to_length(const Eigen::VectorXd& p,
+                                        double length, int metric)
+{
+    const double nrm = pn_norm(p, metric);
+    if (nrm < 1e-300) return p;
+    return p * (length / nrm);
+}
+
+/// Normalise to unit metric norm (jReality `Pn.normalize`).
+inline Eigen::VectorXd pn_normalize(const Eigen::VectorXd& p, int metric)
+{
+    return pn_set_to_length(p, 1.0, metric);
+}
+
+// ── Distance ──────────────────────────────────────────────────────────────────
+
+/// Geodesic distance between two points (jReality `Pn.distanceBetween`).
+///   EUCLIDEAN:  ‖p̂_spatial − q̂_spatial‖ (dehomogenised)
+///   ELLIPTIC:   acos(⟨p̂,q̂⟩)  (spherical angle, clamped to [−1,1])
+///   HYPERBOLIC: acosh(⟨p̂,q̂⟩)  (clamped ≥ 1)
+inline double pn_distance_between(const Eigen::VectorXd& p,
+                                  const Eigen::VectorXd& q,
+                                  int metric)
+{
+    if (metric == PN_EUCLIDEAN) {
+        const auto   pd = pn_dehomogenize(p);
+        const auto   qd = pn_dehomogenize(q);
+        const int    n  = static_cast<int>(pd.size());
+        return (pd.head(n - 1) - qd.head(n - 1)).norm();
+    }
+    const double np = pn_norm(p, metric);
+    const double nq = pn_norm(q, metric);
+    double c = pn_inner_product(p, q, metric) / (np * nq);
+    if (metric == PN_HYPERBOLIC)
+        return std::acosh(std::max(1.0, c));
+    // ELLIPTIC
+    c = std::clamp(c, -1.0, 1.0);
+    return std::acos(c);
+}
+
+// ── Geodesic interpolation ────────────────────────────────────────────────────
+
+/// Constant-speed geodesic interpolation at parameter t ∈ [0,1]
+/// (jReality `Pn.linearInterpolation`).
+///   EUCLIDEAN:  affine blend of the dehomogenised points.
+///   ELLIPTIC:   spherical slerp:
+///                 r(t) = (sin((1−t)d)·p̂ + sin(td)·q̂) / sin(d)
+///   HYPERBOLIC: hyperbolic slerp:
+///                 r(t) = (sinh((1−t)d)·p̂ + sinh(td)·q̂) / sinh(d)
+///   where d = pn_distance_between(p,q,metric) and p̂,q̂ are unit vectors.
+inline Eigen::VectorXd pn_linear_interpolation(const Eigen::VectorXd& p,
+                                               const Eigen::VectorXd& q,
+                                               double t, int metric)
+{
+    if (metric == PN_EUCLIDEAN) {
+        const auto pd = pn_dehomogenize(p);
+        const auto qd = pn_dehomogenize(q);
+        return (1.0 - t) * pd + t * qd;
+    }
+    const auto   ph = pn_normalize(p, metric);
+    const auto   qh = pn_normalize(q, metric);
+    const double d  = pn_distance_between(ph, qh, metric);
+    if (d < 1e-12) return ph;          // coincident — return either endpoint
+    if (metric == PN_ELLIPTIC) {
+        const double s = std::sin(d);
+        return (std::sin((1.0 - t) * d) * ph + std::sin(t * d) * qh) / s;
+    }
+    // HYPERBOLIC
+    const double s = std::sinh(d);
+    return (std::sinh((1.0 - t) * d) * ph + std::sinh(t * d) * qh) / s;
+}
+
+} // namespace conformallab