ConformalLabpp/code/include/pn_geometry.hpp

#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