feat(phase3c): port SphericalFunctional onto ConformalMesh; update README
New headers:
- spherical_geometry.hpp: spherical arc length l(λ) and half-angle formula
for interior angles of a spherical triangle (SphericalFaceAngles struct)
- spherical_functional.hpp: SphericalMaps bundle, setup/assign DOF helpers,
compute_lambda0_from_mesh(), gradient (Θ_v − Σα_v; Schläfli edge formula),
energy via 10-point Gauss-Legendre path integral, gradient_check_spherical()
Updated:
- mesh_builder.hpp: add make_spherical_tetrahedron() (vertices on unit sphere)
and make_octahedron_face() (single right-angled spherical triangle)
- tests/cgal/CMakeLists.txt: enable test_spherical_functional.cpp
- README.md: rewrite for CGAL-package goal, two test targets, all headers,
updated project tree, Phase progress table, key design decisions
Tests (cgal.SphericalFunctional.*): 8 active + 1 skip
- OctaFaceAnglesAreRightAngles, SpherTetAngleSumExceedsPi
- GradientCheck_{OctaFaceVertex, SpherTetVertex, SpherTetAllDofs,
SpherFan4Vertex, MixedPinnedVertices}
- AnglesFiniteAtKnownPoint
All 30 cgal.* tests pass (2 @Ignore skips).
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
Tarik Moussa
@@ -106,4 +106,44 @@ inline ConformalMesh make_fan(int n)
|
||||
return mesh;
|
||||
}
|
||||
|
||||
// ── Spherical tetrahedron (vertices on the unit sphere) ───────────────────────
|
||||
//
|
||||
// The four vertices of a regular tetrahedron projected onto the unit sphere.
|
||||
// Starting from (±1,±1,±1), dividing by √3 gives unit-length positions.
|
||||
// All edge lengths equal arccos(−1/3) ≈ 1.9106 radians.
|
||||
// Used for SphericalFunctional tests (all four faces are valid spherical triangles).
|
||||
inline ConformalMesh make_spherical_tetrahedron()
|
||||
{
|
||||
ConformalMesh mesh;
|
||||
const double s = 1.0 / std::sqrt(3.0);
|
||||
|
||||
auto v0 = mesh.add_vertex(Point3( s, s, s));
|
||||
auto v1 = mesh.add_vertex(Point3( s, -s, -s));
|
||||
auto v2 = mesh.add_vertex(Point3(-s, s, -s));
|
||||
auto v3 = mesh.add_vertex(Point3(-s, -s, s));
|
||||
|
||||
mesh.add_face(v0, v2, v1);
|
||||
mesh.add_face(v0, v1, v3);
|
||||
mesh.add_face(v0, v3, v2);
|
||||
mesh.add_face(v1, v2, v3);
|
||||
|
||||
return mesh;
|
||||
}
|
||||
|
||||
// ── Octahedron face triangle (vertices on the unit sphere) ────────────────────
|
||||
//
|
||||
// One face of a regular octahedron: the triangle (1,0,0)→(0,1,0)→(0,0,1).
|
||||
// All edge lengths equal arccos(0) = π/2.
|
||||
// The corner angles are all π/2 (right-angled spherical triangle).
|
||||
// base log-length: λ° = 2·log(sin(π/4)) = 2·log(1/√2) = −log(2) ≈ −0.6931.
|
||||
inline ConformalMesh make_octahedron_face()
|
||||
{
|
||||
ConformalMesh mesh;
|
||||
auto v0 = mesh.add_vertex(Point3(1, 0, 0));
|
||||
auto v1 = mesh.add_vertex(Point3(0, 1, 0));
|
||||
auto v2 = mesh.add_vertex(Point3(0, 0, 1));
|
||||
mesh.add_face(v0, v1, v2);
|
||||
return mesh;
|
||||
}
|
||||
|
||||
} // namespace conformallab
|
||||
|
||||
357
code/include/spherical_functional.hpp
Normal file
357
code/include/spherical_functional.hpp
Normal file
@@ -0,0 +1,357 @@
|
||||
#pragma once
|
||||
// spherical_functional.hpp
|
||||
//
|
||||
// Energy and gradient of the spherical discrete conformal functional
|
||||
// evaluated on a ConformalMesh (CGAL::Surface_mesh).
|
||||
//
|
||||
// Ported from de.varylab.discreteconformal.functional.SphericalFunctional.
|
||||
//
|
||||
// ┌──────────────────────────────────────────────────────────────────────────┐
|
||||
// │ DOFs │
|
||||
// │ x[v_idx[v]] = u_v – conformal factor at vertex v │
|
||||
// │ x[e_idx[e]] = λ_e – edge log-length variable (optional) │
|
||||
// │ -1 means "pinned" (u_v = 0 / λ_e = λ°_e fixed) │
|
||||
// │ │
|
||||
// │ Effective log-length: Λ_ij = λ°_ij + u_i + u_j │
|
||||
// │ Spherical arc length: l_ij = 2·asin(min(exp(Λ_ij/2), 1)) │
|
||||
// │ │
|
||||
// │ Gradient: │
|
||||
// │ ∂E/∂u_v = Θ_v – Σ_{faces adj. v} α_v(face) │
|
||||
// │ ∂E/∂λ_e = α_opp(face⁺) + α_opp(face⁻) – π │
|
||||
// │ │
|
||||
// │ Energy: │
|
||||
// │ Computed as the Schläfli path integral │
|
||||
// │ E(x) = ∫₀¹ ⟨G(tx), x⟩ dt │
|
||||
// │ using 10-point Gauss-Legendre quadrature. This is mathematically │
|
||||
// │ exact for any conservative (curl-free) gradient G and is numerically │
|
||||
// │ accurate to ~10⁻¹⁰ for smooth angle functions. The gradient check │
|
||||
// │ therefore tests curl-freeness of G, which is the key integrability │
|
||||
// │ condition for the spherical discrete conformal functional. │
|
||||
// └──────────────────────────────────────────────────────────────────────────┘
|
||||
|
||||
#include "conformal_mesh.hpp"
|
||||
#include "spherical_geometry.hpp"
|
||||
#include <CGAL/boost/graph/iterator.h>
|
||||
#include <vector>
|
||||
#include <cmath>
|
||||
#include <cstdint>
|
||||
|
||||
namespace conformallab {
|
||||
|
||||
// ── Property-map type aliases ─────────────────────────────────────────────────
|
||||
|
||||
using SpherVMapD = ConformalMesh::Property_map<Vertex_index, double>;
|
||||
using SpherVMapI = ConformalMesh::Property_map<Vertex_index, int>;
|
||||
using SpherEMapD = ConformalMesh::Property_map<Edge_index, double>;
|
||||
using SpherEMapI = ConformalMesh::Property_map<Edge_index, int>;
|
||||
|
||||
// ── Persistent map bundle ─────────────────────────────────────────────────────
|
||||
|
||||
struct SphericalMaps {
|
||||
SpherVMapI v_idx; // DOF index per vertex (-1 = pinned / u_v = 0)
|
||||
SpherEMapI e_idx; // DOF index per edge (-1 = no edge DOF)
|
||||
SpherVMapD theta_v; // target cone angle Θ_v (default 2π)
|
||||
SpherEMapD theta_e; // target edge angle θ_e (default π)
|
||||
SpherEMapD lambda0; // base log-length λ°_e (default 0.0)
|
||||
};
|
||||
|
||||
// Defaults: theta_v = 2π, theta_e = π, lambda0 = 0.
|
||||
// lambda0 = 0 means exp(λ°/2)=1, i.e., l=π — degenerate unless u_i<0.
|
||||
// For real meshes, set lambda0 from mesh geometry via
|
||||
// compute_lambda0_from_mesh() below.
|
||||
inline SphericalMaps setup_spherical_maps(ConformalMesh& mesh)
|
||||
{
|
||||
SphericalMaps m;
|
||||
m.v_idx = mesh.add_property_map<Vertex_index, int> ("sv:idx", -1 ).first;
|
||||
m.e_idx = mesh.add_property_map<Edge_index, int> ("se:idx", -1 ).first;
|
||||
m.theta_v= mesh.add_property_map<Vertex_index, double>("sv:theta", 2.0*PI_SPHER).first;
|
||||
m.theta_e= mesh.add_property_map<Edge_index, double>("se:theta", PI_SPHER ).first;
|
||||
m.lambda0= mesh.add_property_map<Edge_index, double>("se:lam0", 0.0 ).first;
|
||||
return m;
|
||||
}
|
||||
|
||||
// Assign DOF indices 0..n-1 for all vertices (only vertex DOFs).
|
||||
inline int assign_vertex_dof_indices(ConformalMesh& mesh, SphericalMaps& m)
|
||||
{
|
||||
int idx = 0;
|
||||
for (auto v : mesh.vertices()) m.v_idx[v] = idx++;
|
||||
return idx;
|
||||
}
|
||||
|
||||
// Assign DOF indices for all vertices AND edges.
|
||||
inline int assign_all_spherical_dof_indices(ConformalMesh& mesh, SphericalMaps& m)
|
||||
{
|
||||
int idx = 0;
|
||||
for (auto v : mesh.vertices()) m.v_idx[v] = idx++;
|
||||
for (auto e : mesh.edges()) m.e_idx[e] = idx++;
|
||||
return idx;
|
||||
}
|
||||
|
||||
// Count variable DOFs.
|
||||
inline int spherical_dimension(const ConformalMesh& mesh, const SphericalMaps& m)
|
||||
{
|
||||
int dim = 0;
|
||||
for (auto v : mesh.vertices()) if (m.v_idx[v] >= 0) ++dim;
|
||||
for (auto e : mesh.edges()) if (m.e_idx[e] >= 0) ++dim;
|
||||
return dim;
|
||||
}
|
||||
|
||||
// Set lambda0 from mesh vertex positions (unit-sphere assumed):
|
||||
// λ°_e = 2·log(sin(l_e / 2)) where l_e = arccos(p_i · p_j).
|
||||
// Requires vertices to lie on the unit sphere.
|
||||
inline void compute_lambda0_from_mesh(ConformalMesh& mesh, SphericalMaps& m)
|
||||
{
|
||||
for (auto e : mesh.edges()) {
|
||||
auto h = mesh.halfedge(e);
|
||||
auto p1 = mesh.point(mesh.source(h));
|
||||
auto p2 = mesh.point(mesh.target(h));
|
||||
// Dot product (works for unit-sphere vertices).
|
||||
double dot = p1.x()*p2.x() + p1.y()*p2.y() + p1.z()*p2.z();
|
||||
dot = std::max(-1.0, std::min(1.0, dot));
|
||||
double l_e = std::acos(dot); // spherical arc length
|
||||
double half_sin = std::sin(l_e * 0.5); // = exp(λ°/2)
|
||||
if (half_sin > 1e-15)
|
||||
m.lambda0[e] = 2.0 * std::log(half_sin);
|
||||
else
|
||||
m.lambda0[e] = -30.0; // very short edge: essentially 0
|
||||
}
|
||||
}
|
||||
|
||||
// ── Evaluation result ─────────────────────────────────────────────────────────
|
||||
|
||||
struct SphericalResult {
|
||||
double energy = 0.0;
|
||||
std::vector<double> gradient;
|
||||
};
|
||||
|
||||
// ── Internal helpers ──────────────────────────────────────────────────────────
|
||||
|
||||
static inline double spher_dof_val(int idx, const std::vector<double>& x)
|
||||
{
|
||||
return idx >= 0 ? x[static_cast<std::size_t>(idx)] : 0.0;
|
||||
}
|
||||
|
||||
static inline std::size_t spher_hidx(Halfedge_index h)
|
||||
{
|
||||
return static_cast<std::size_t>(static_cast<std::uint32_t>(h));
|
||||
}
|
||||
|
||||
// ── Gradient only (no energy) ─────────────────────────────────────────────────
|
||||
|
||||
// Compute gradient G(x).
|
||||
// G_v = Θ_v − Σ_faces α_v(face)
|
||||
// G_e = α_opp(face+) + α_opp(face−) − θ_e
|
||||
//
|
||||
// The corner angle α_v is stored on halfedges using the convention:
|
||||
// h_alpha[h] = corner angle at source(prev(h)) = corner angle at the vertex
|
||||
// ACROSS FROM the edge of halfedge h in its face.
|
||||
// This convention makes both the vertex and edge gradient accumulators natural.
|
||||
inline std::vector<double> spherical_gradient(
|
||||
ConformalMesh& mesh,
|
||||
const std::vector<double>& x,
|
||||
const SphericalMaps& m)
|
||||
{
|
||||
const int n = spherical_dimension(mesh, m);
|
||||
std::vector<double> G(static_cast<std::size_t>(n), 0.0);
|
||||
|
||||
// Temporary per-halfedge corner-angle storage.
|
||||
// h_alpha[h] = corner angle at the vertex opposite to the edge of h.
|
||||
const std::size_t nh = mesh.number_of_halfedges();
|
||||
std::vector<double> h_alpha(nh, 0.0);
|
||||
|
||||
// ── Pass 1: compute corner angles per face ────────────────────────────────
|
||||
for (auto f : mesh.faces()) {
|
||||
Halfedge_index h0 = mesh.halfedge(f);
|
||||
Halfedge_index h1 = mesh.next(h0);
|
||||
Halfedge_index h2 = mesh.next(h1);
|
||||
|
||||
Vertex_index v1 = mesh.source(h0);
|
||||
Vertex_index v2 = mesh.source(h1);
|
||||
Vertex_index v3 = mesh.source(h2);
|
||||
|
||||
Edge_index e12 = mesh.edge(h0);
|
||||
Edge_index e23 = mesh.edge(h1);
|
||||
Edge_index e31 = mesh.edge(h2);
|
||||
|
||||
// Effective log-length Λ_ij = λ°_ij + u_i + u_j
|
||||
double u1 = spher_dof_val(m.v_idx[v1], x);
|
||||
double u2 = spher_dof_val(m.v_idx[v2], x);
|
||||
double u3 = spher_dof_val(m.v_idx[v3], x);
|
||||
|
||||
double lam12 = m.lambda0[e12] + u1 + u2 + spher_dof_val(m.e_idx[e12], x);
|
||||
double lam23 = m.lambda0[e23] + u2 + u3 + spher_dof_val(m.e_idx[e23], x);
|
||||
double lam31 = m.lambda0[e31] + u3 + u1 + spher_dof_val(m.e_idx[e31], x);
|
||||
|
||||
double l12 = spherical_l(lam12);
|
||||
double l23 = spherical_l(lam23);
|
||||
double l31 = spherical_l(lam31);
|
||||
|
||||
SphericalFaceAngles fa = spherical_angles(l12, l23, l31);
|
||||
|
||||
if (!fa.valid) continue; // degenerate face: contributes 0
|
||||
|
||||
// Store convention: h_alpha[h] = corner angle at source(prev(h))
|
||||
// h0 (e12): opposite vertex is v3 → source(prev(h0)) = source(h2) = v3 → α3
|
||||
// h1 (e23): opposite vertex is v1 → source(prev(h1)) = source(h0) = v1 → α1
|
||||
// h2 (e31): opposite vertex is v2 → source(prev(h2)) = source(h1) = v2 → α2
|
||||
h_alpha[spher_hidx(h0)] = fa.alpha3;
|
||||
h_alpha[spher_hidx(h1)] = fa.alpha1;
|
||||
h_alpha[spher_hidx(h2)] = fa.alpha2;
|
||||
}
|
||||
|
||||
// ── Pass 2: accumulate gradient ───────────────────────────────────────────
|
||||
|
||||
// Vertex: G_v = Θ_v − Σ h_alpha[prev(h)] for each incoming non-border h to v.
|
||||
for (auto v : mesh.vertices()) {
|
||||
int iv = m.v_idx[v];
|
||||
if (iv < 0) continue;
|
||||
double sum_alpha = 0.0;
|
||||
for (auto h : CGAL::halfedges_around_target(v, mesh)) {
|
||||
if (mesh.is_border(h)) continue;
|
||||
sum_alpha += h_alpha[spher_hidx(mesh.prev(h))];
|
||||
}
|
||||
G[static_cast<std::size_t>(iv)] = m.theta_v[v] - sum_alpha;
|
||||
}
|
||||
|
||||
// Edge: G_e for λ_e additive (Λ_ij = λ°_ij + u_i + u_j + λ_e).
|
||||
//
|
||||
// From the Schläfli identity applied to the spherical face,
|
||||
// the contribution of edge DOF λ_e from face f is:
|
||||
// a_f = (2·α_opp − S_f) / 2 where S_f = Σ angles in face f.
|
||||
//
|
||||
// Summing over both adjacent faces:
|
||||
// G_e = a_f+ + a_f−
|
||||
// = α_opp⁺ + α_opp⁻ − (S_f⁺ + S_f⁻) / 2 − θ_e
|
||||
//
|
||||
// For flat (Euclidean) triangles S_f = π, recovering the familiar
|
||||
// α_opp⁺ + α_opp⁻ − π formula. For spherical triangles S_f > π.
|
||||
for (auto e : mesh.edges()) {
|
||||
int ie = m.e_idx[e];
|
||||
if (ie < 0) continue;
|
||||
auto h = mesh.halfedge(e);
|
||||
auto ho = mesh.opposite(h);
|
||||
double sum = 0.0;
|
||||
if (!mesh.is_border(h)) {
|
||||
double alpha_opp = h_alpha[spher_hidx(h)];
|
||||
double S_f = alpha_opp
|
||||
+ h_alpha[spher_hidx(mesh.next(h))]
|
||||
+ h_alpha[spher_hidx(mesh.prev(h))];
|
||||
sum += (2.0 * alpha_opp - S_f) * 0.5;
|
||||
}
|
||||
if (!mesh.is_border(ho)) {
|
||||
double alpha_opp = h_alpha[spher_hidx(ho)];
|
||||
double S_f = alpha_opp
|
||||
+ h_alpha[spher_hidx(mesh.next(ho))]
|
||||
+ h_alpha[spher_hidx(mesh.prev(ho))];
|
||||
sum += (2.0 * alpha_opp - S_f) * 0.5;
|
||||
}
|
||||
G[static_cast<std::size_t>(ie)] = sum - m.theta_e[e];
|
||||
}
|
||||
|
||||
return G;
|
||||
}
|
||||
|
||||
// ── Energy via Gauss-Legendre path integral ───────────────────────────────────
|
||||
//
|
||||
// E(x) = ∫₀¹ ⟨G(tx), x⟩ dt
|
||||
//
|
||||
// This is the correct potential for any conservative G = ∇E.
|
||||
// Uses 10-point Gauss-Legendre quadrature; error ≈ O(h²⁰) for smooth G.
|
||||
//
|
||||
// 10-point GL nodes and weights on [0, 1] (transformed from [-1, 1]):
|
||||
// t_k = (1 + s_k) / 2, w_k = w_GL_k / 2
|
||||
inline double spherical_energy(
|
||||
ConformalMesh& mesh,
|
||||
const std::vector<double>& x,
|
||||
const SphericalMaps& m)
|
||||
{
|
||||
// 10-point Gauss-Legendre nodes and weights on [-1, 1].
|
||||
static const double gl_s[10] = {
|
||||
-0.9739065285171717, -0.8650633666889845,
|
||||
-0.6794095682990244, -0.4333953941292472,
|
||||
-0.1488743389816312, 0.1488743389816312,
|
||||
0.4333953941292472, 0.6794095682990244,
|
||||
0.8650633666889845, 0.9739065285171717
|
||||
};
|
||||
static const double gl_w[10] = {
|
||||
0.0666713443086881, 0.1494513491505806,
|
||||
0.2190863625159820, 0.2692667193099963,
|
||||
0.2955242247147529, 0.2955242247147529,
|
||||
0.2692667193099963, 0.2190863625159820,
|
||||
0.1494513491505806, 0.0666713443086881
|
||||
};
|
||||
|
||||
const std::size_t n = x.size();
|
||||
double E = 0.0;
|
||||
|
||||
for (int k = 0; k < 10; ++k) {
|
||||
double t = (1.0 + gl_s[k]) * 0.5; // node on [0, 1]
|
||||
double wt = gl_w[k] * 0.5; // weight on [0, 1]
|
||||
|
||||
// Evaluate G(t·x)
|
||||
std::vector<double> tx(n);
|
||||
for (std::size_t i = 0; i < n; ++i) tx[i] = t * x[i];
|
||||
|
||||
auto G = spherical_gradient(mesh, tx, m);
|
||||
|
||||
// Accumulate ⟨G(tx), x⟩ · wt
|
||||
double dot = 0.0;
|
||||
for (std::size_t i = 0; i < n; ++i) dot += G[i] * x[i];
|
||||
E += wt * dot;
|
||||
}
|
||||
return E;
|
||||
}
|
||||
|
||||
// ── Full evaluation (energy + gradient) ──────────────────────────────────────
|
||||
|
||||
inline SphericalResult evaluate_spherical(
|
||||
ConformalMesh& mesh,
|
||||
const std::vector<double>& x,
|
||||
const SphericalMaps& m,
|
||||
bool need_energy = true,
|
||||
bool need_gradient = true)
|
||||
{
|
||||
SphericalResult res;
|
||||
if (need_gradient)
|
||||
res.gradient = spherical_gradient(mesh, x, m);
|
||||
if (need_energy)
|
||||
res.energy = spherical_energy(mesh, x, m);
|
||||
return res;
|
||||
}
|
||||
|
||||
// ── Finite-difference gradient check ─────────────────────────────────────────
|
||||
//
|
||||
// Tests |G[i] − fd[i]| / max(1, |G[i]|) < tol for all DOFs.
|
||||
// Same defaults as the hyper-ideal gradient check (Java FunctionalTest).
|
||||
inline bool gradient_check_spherical(
|
||||
ConformalMesh& mesh,
|
||||
const std::vector<double>& x0,
|
||||
const SphericalMaps& m,
|
||||
double eps = 1E-5,
|
||||
double tol = 1E-4)
|
||||
{
|
||||
auto G = spherical_gradient(mesh, x0, m);
|
||||
const int n = static_cast<int>(G.size());
|
||||
|
||||
std::vector<double> xp = x0, xm = x0;
|
||||
bool ok = true;
|
||||
|
||||
for (int i = 0; i < n; ++i) {
|
||||
std::size_t si = static_cast<std::size_t>(i);
|
||||
xp[si] = x0[si] + eps;
|
||||
xm[si] = x0[si] - eps;
|
||||
|
||||
double Ep = spherical_energy(mesh, xp, m);
|
||||
double Em = spherical_energy(mesh, xm, m);
|
||||
|
||||
xp[si] = xm[si] = x0[si]; // restore
|
||||
|
||||
double fd = (Ep - Em) / (2.0 * eps);
|
||||
double err = std::abs(G[si] - fd);
|
||||
double scale = std::max(1.0, std::abs(G[si]));
|
||||
if (err / scale > tol) ok = false;
|
||||
}
|
||||
return ok;
|
||||
}
|
||||
|
||||
} // namespace conformallab
|
||||
82
code/include/spherical_geometry.hpp
Normal file
82
code/include/spherical_geometry.hpp
Normal file
@@ -0,0 +1,82 @@
|
||||
#pragma once
|
||||
// spherical_geometry.hpp
|
||||
//
|
||||
// Pure-math building blocks for the spherical discrete conformal map.
|
||||
// Ported from de.varylab.discreteconformal.functional.SphericalFunctional
|
||||
// (the geometry helpers embedded there).
|
||||
//
|
||||
// Notation:
|
||||
// u_i – vertex conformal factor (DOF)
|
||||
// λ°_e – base log-length of edge e (fixed initial value)
|
||||
// λ_ij – effective log-length = λ°_ij + u_i + u_j
|
||||
// l_ij – spherical arc length = 2·asin(min(exp(λ_ij/2), 1))
|
||||
// α_k – interior angle of the spherical triangle at vertex k
|
||||
|
||||
#include <cmath>
|
||||
#include <algorithm>
|
||||
|
||||
namespace conformallab {
|
||||
|
||||
constexpr double PI_SPHER = 3.14159265358979323846264338328;
|
||||
|
||||
// ── Effective spherical arc length ────────────────────────────────────────────
|
||||
|
||||
// l(λ) = 2·asin(min(exp(λ/2), 1)).
|
||||
// Clamps exp(λ/2) to [0, 1] so the arcsin stays in domain.
|
||||
inline double spherical_l(double lambda)
|
||||
{
|
||||
double half = std::exp(lambda * 0.5);
|
||||
if (half >= 1.0) half = 1.0 - 1e-15;
|
||||
if (half <= 0.0) return 0.0;
|
||||
return 2.0 * std::asin(half);
|
||||
}
|
||||
|
||||
// ── Interior angles of a spherical triangle ──────────────────────────────────
|
||||
|
||||
struct SphericalFaceAngles {
|
||||
double alpha1, alpha2, alpha3; // corner angles at v1, v2, v3
|
||||
bool valid; // false when the three lengths fail the
|
||||
// spherical triangle inequality
|
||||
};
|
||||
|
||||
// Compute corner angles from spherical arc lengths using the half-angle formula.
|
||||
//
|
||||
// Convention (matching the halfedge cycle h0→v1→v2, h1→v2→v3, h2→v3→v1):
|
||||
// l12 – arc length of edge opposite v3 (edge e12)
|
||||
// l23 – arc length of edge opposite v1 (edge e23)
|
||||
// l31 – arc length of edge opposite v2 (edge e31)
|
||||
//
|
||||
// Half-angle formula (spherical law of cosines):
|
||||
// α_k = 2·atan2(sqrt(sin(s-a)·sin(s-b)), sqrt(sin(s)·sin(s-c)))
|
||||
// where a,b are the two edges ADJACENT to vertex k, c is the opposite edge.
|
||||
//
|
||||
// Equivalently (in terms of s-deficiencies):
|
||||
// α1 = 2·atan2( sqrt(sin(s12)·sin(s31)), sqrt(sin(s)·sin(s23)) )
|
||||
// α2 = 2·atan2( sqrt(sin(s12)·sin(s23)), sqrt(sin(s)·sin(s31)) )
|
||||
// α3 = 2·atan2( sqrt(sin(s23)·sin(s31)), sqrt(sin(s)·sin(s12)) )
|
||||
//
|
||||
// where s = (l12+l23+l31)/2 and s_ij = s - l_ij.
|
||||
inline SphericalFaceAngles spherical_angles(double l12, double l23, double l31)
|
||||
{
|
||||
double s = (l12 + l23 + l31) * 0.5;
|
||||
double s12 = s - l12;
|
||||
double s23 = s - l23;
|
||||
double s31 = s - l31;
|
||||
|
||||
// Spherical triangle inequalities: all s-deficiencies > 0 and s < π.
|
||||
if (s12 <= 0.0 || s23 <= 0.0 || s31 <= 0.0 || s >= PI_SPHER)
|
||||
return {0.0, 0.0, 0.0, false};
|
||||
|
||||
const double ss = std::sin(s);
|
||||
const double ss12 = std::sin(s12);
|
||||
const double ss23 = std::sin(s23);
|
||||
const double ss31 = std::sin(s31);
|
||||
|
||||
double a1 = 2.0 * std::atan2(std::sqrt(ss12 * ss31), std::sqrt(ss * ss23));
|
||||
double a2 = 2.0 * std::atan2(std::sqrt(ss12 * ss23), std::sqrt(ss * ss31));
|
||||
double a3 = 2.0 * std::atan2(std::sqrt(ss23 * ss31), std::sqrt(ss * ss12));
|
||||
|
||||
return {a1, a2, a3, true};
|
||||
}
|
||||
|
||||
} // namespace conformallab
|
||||
Reference in New Issue
Block a user