ConformalLabpp/code/include/stereographic_layout.hpp

// Copyright (c) 2024-2026 Tarik Moussa.
// SPDX-License-Identifier: MIT

// stereographic_layout.hpp
//
// Phase 9d.3 — Stereographic projection for spherical DCE output.
//
// Converts a spherical layout (points on S²) to a 2-D conformal map via:
// 1. Stereographic projection: S² → ℂ ∪ {∞}, mapping the sphere to the complex plane.
// 2. Möbius centring: centres the resulting point cloud for canonical position.
//
// Mathematical reference:
//   Stereographic projection from the north pole (0,0,1):
//     (x,y,z) ↦ (x/(1-z), y/(1-z)) in ℂ (complex coordinate u+iv).
//   North pole (0,0,1) maps to ∞ (removed from the layout).
//   South pole (0,0,-1) maps to (0,0) in ℂ.
//   The projection is conformal (angle-preserving).
//
// Möbius centring: apply a Möbius transformation to centre the layout
//   (e.g. shift the centroid to the origin, possibly scale/rotate).
//
// Java source (ported from):
//   unwrapper/StereographicUnwrapper.java (266 lines)
//   The supporting math/CP1 + ComplexUtility.stereographic operations
//   (deliberately NOT ported — redundant with std::complex).

#pragma once

#include "conformal_mesh.hpp"
#include "layout.hpp"
#include <complex>
#include <vector>
#include <cmath>
#include <array>

namespace conformallab {

// ────────────────────────────────────────────────────────────────────────────
// Stereographic Projection: S² → ℂ
// ────────────────────────────────────────────────────────────────────────────

/// Stereographic projection from the north pole (0, 0, 1).
/// Maps a point on the unit sphere S² to the complex plane ℂ.
/// North pole (0,0,1) projects to ∞ (not representable; returns NaN).
/// South pole (0,0,-1) projects to 0+0i.
///
/// Formula: (x,y,z) ↦ x/(1-z) + i·y/(1-z)
inline std::complex<double> stereographic_project(double x, double y, double z)
{
    const double denom = 1.0 - z;
    if (std::abs(denom) < 1e-15) {
        // North pole (z ≈ 1) — maps to ∞.
        // Return NaN to signal infinity.
        return std::complex<double>(std::nan(""), std::nan(""));
    }
    return std::complex<double>(x / denom, y / denom);
}

/// Stereographic projection of a 3-D point (as Point3).
inline std::complex<double> stereographic_project(const Point3& p)
{
    return stereographic_project(p.x(), p.y(), p.z());
}

// ────────────────────────────────────────────────────────────────────────────
// Möbius Centring
// ────────────────────────────────────────────────────────────────────────────

/// Simple centring: translate the point cloud so that its centroid
/// is at the origin (u+iv = 0).
inline void centre_at_origin(std::vector<std::complex<double>>& points)
{
    if (points.empty()) return;

    // Compute centroid.
    std::complex<double> centroid(0.0, 0.0);
    int n_valid = 0;
    for (const auto& z : points) {
        if (std::isfinite(z.real()) && std::isfinite(z.imag())) {
            centroid += z;
            n_valid++;
        }
    }
    if (n_valid <= 0) return;

    centroid /= static_cast<double>(n_valid);

    // Translate: z' = z - centroid.
    for (auto& z : points) {
        if (std::isfinite(z.real()) && std::isfinite(z.imag())) {
            z -= centroid;
        }
    }
}

// ────────────────────────────────────────────────────────────────────────────
// Stereographic Layout: S² → ℂ (2-D)
// ────────────────────────────────────────────────────────────────────────────

/// Convert a spherical layout (3-D points on S²) to a 2-D conformal map
/// via stereographic projection.
///
/// Output: a Layout2D where:
///   - uv[v.idx()] = (Re, Im) of the stereographic projection of the 3-D point.
///   - The north pole is excluded (uv[v] = NaN for projections at ∞).
///
/// Möbius centring: the resulting layout is centred at the origin.
///
/// \param mesh    Input surface mesh.
/// \param layout  Input spherical layout (3-D points on S²).
/// \return        Output Layout2D in the complex plane (ℂ).
inline Layout2D stereographic_layout(
    const ConformalMesh&          mesh,
    const Layout3D&               layout)
{
    Layout2D result;
    result.uv.resize(mesh.number_of_vertices());
    result.halfedge_uv.resize(mesh.number_of_halfedges());

    // Step 1: Stereographic projection for each vertex.
    std::vector<std::complex<double>> complex_points;
    complex_points.reserve(mesh.number_of_vertices());

    for (auto v : mesh.vertices()) {
        const auto& p3d = layout.pos[v.idx()];
        // Convert Eigen::Vector3d to Point3-like coordinates.
        double x = p3d[0], y = p3d[1], z = p3d[2];
        auto z_complex = stereographic_project(x, y, z);
        complex_points.push_back(z_complex);

        // Store as Eigen::Vector2d (Re, Im).
        result.uv[v.idx()] = Eigen::Vector2d(z_complex.real(), z_complex.imag());
    }

    // Step 2: Möbius centring.
    centre_at_origin(complex_points);

    // Update uv after centring.
    for (auto v : mesh.vertices()) {
        const auto& z = complex_points[v.idx()];
        result.uv[v.idx()] = Eigen::Vector2d(z.real(), z.imag());
    }

    // Step 3: Halfedge UV (for texture atlasing).
    // Copy the primary vertex UV to each halfedge's source.
    for (auto h : mesh.halfedges()) {
        Vertex_index src = mesh.source(h);
        result.halfedge_uv[h.idx()] = result.uv[src.idx()];
    }

    return result;
}

// ────────────────────────────────────────────────────────────────────────────
// Inverse Stereographic Projection: ℂ → S²
// ────────────────────────────────────────────────────────────────────────────

/// Inverse stereographic projection: ℂ → S².
/// Given a complex number z = u + iv, recover the 3-D point on the unit sphere.
///
/// Formula: (u,v) ↦ (2u/(1+u²+v²), 2v/(1+u²+v²), (u²+v²-1)/(u²+v²+1))
/// Inverse of: (x,y,z) ↦ (x/(1-z), y/(1-z)).
///
/// The origin (u,v) = (0,0) maps back to (0,0,-1) (south pole).
inline Point3 inverse_stereographic_project(std::complex<double> z)
{
    double u = z.real();
    double v = z.imag();

    double u2_plus_v2 = u * u + v * v;
    double denom = 1.0 + u2_plus_v2;

    double x = 2.0 * u / denom;
    double y = 2.0 * v / denom;
    double zz = (u2_plus_v2 - 1.0) / denom;

    return Point3(x, y, zz);
}

/// Inverse stereographic projection from a 2-D layout point.
inline Point3 inverse_stereographic_project(const Eigen::Vector2d& uv)
{
    return inverse_stereographic_project(std::complex<double>(uv.x(), uv.y()));
}

/// Round-trip validation: project a 3-D point to 2-D and back.
/// Returns the error (distance on S²) between the original and recovered point.
inline double stereographic_roundtrip_error(const Point3& original)
{
    auto z = stereographic_project(original);
    if (!std::isfinite(z.real()) || !std::isfinite(z.imag())) {
        return std::numeric_limits<double>::infinity();  // north pole
    }
    auto recovered = inverse_stereographic_project(z);

    // Distance on the unit sphere: ‖p - q‖.
    double dx = original.x() - recovered.x();
    double dy = original.y() - recovered.y();
    double dz = original.z() - recovered.z();
    return std::sqrt(dx * dx + dy * dy + dz * dz);
}

}  // namespace conformallab