feat(geometry): pn_geometry.hpp — Pn projective-metric substrate (jReality port)
Ports the small but pervasive de.jreality.math.Pn surface that the Java
algorithmic core uses across every not-yet-ported geometric phase:
HyperbolicLayout, SphericalLayout, FundamentalPolygon (9c), KoebePolyhedron
(10c'), quasi-isothermic (10e), hyperelliptic theta, CircleDomain (11c).
Porting it once unblocks all downstream consumers instead of re-deriving
the metric ad hoc per phase.
pn_geometry.hpp — header-only, Eigen, ~120 LOC:
PnMetric { EUCLIDEAN=0, ELLIPTIC=+1, HYPERBOLIC=-1 } matching jReality.
pn_inner_product — bilinear form per signature.
pn_norm / pn_dehomogenize / pn_set_to_length / pn_normalize.
pn_distance_between — Euclidean (spatial), elliptic (acos), hyperbolic (acosh).
pn_linear_interpolation — affine (E) and slerp (S/H constant-speed geodesic).
HYPERBOLIC uses the timelike-positive ("upper-sheet") convention, identical to
the already-verified projective_math.hpp::hyperbolicDistance; pn_distance_between
is regression-anchored against it in the tests.
test_pn_geometry.cpp — 6 tests, all GREEN:
InnerProductSignatures, EuclideanDistance, EllipticDistanceIsAngle,
HyperbolicDistanceClosedFormAndAnchor (+ projective_math.hpp anchor),
NormAndScaling, LinearInterpolationGeodesic.
doc/roadmap/java-parity.md — new "Infrastructure / support layers" section:
de.jreality.math.Pn → pn_geometry.hpp (partial, 6 fns covering core usage)
de.jreality.math.Rn → Eigen (no separate port needed)
MatrixBuilder → not yet (only needed for hyperelliptic theta)
conformallab XML types → not planned (GUI persistence, not algorithmic)
246/246 cgal tests pass.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Tarik Moussa
156
code/include/pn_geometry.hpp
Normal file
156
code/include/pn_geometry.hpp
Normal file
@@ -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
|
||||
Reference in New Issue
Block a user