UT_TriangleMesh.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * NAME:        UT_TriangleMesh.h (UT Library, C++)
 *
 * COMMENTS: A class to represents a constant (unmodifiable) abstract
 * (topological) triangulation of a 2-manifold (surface) possibly with
 * boundary, providing the usual traversal methods. Note that no geometric
 * information is kept in this class.
 *
 * Be particularly careful about what a manifold "boundary point" is! The
 * neighbourhood of a boundary point must be homeomorphic to a half-disc
 * (in particular, the neighbourhood minus the point itself must be connected),
 * meaning pinched boundaries (such as what you get from two triangles sharing
 * a single point) cannot be represented by this class.
 *
 */

#ifndef __UT_TriangleMesh__
#define __UT_TriangleMesh__

#include "UT_API.h"
#include "UT_ValArray.h"
#include "UT_ParallelUtil.h"
#include "UT_Quaternion.h"
#include <complex>

// We use the usual Houdini terminology of "points", "vertices", and half-edges
// (hedges for short): a vertex belongs to a single triangle and is mapped to a
// point which exists independently. A half-edge is the share of a single
// triangle of an edge and is oriented consistently with triangle's winding
// with a source and a destination point (and vertex). A neighbouring triangle
// sharing the same edge has its own corresponding half-edge with the same
// endpoints. The two half-edges are said to be *equivalent* and point to each
// other through the sym() method. Boundary half-edges have negative sym().
// Manifold requirement dictates that each non-boundary half-edge has exactly
// one sym() which in turn would have that half-edge as its sym(). Furthermore,
// the source of a half-edge is always the destination of its sym() and
// vice versa.
//
// Triangles are numbered between 0 and numTriangles() - 1, points are numbered
// 0 to numPoints() - 1, and vertices (or equivalently half-edges) are numbered
// 0..numHedges() - 1. Note that numHedges() is the same as numVertices() and
// is 3 * numTriangles().
//
// NB. Any correspondence between elements of the mesh and the world outside
// should be tracked externally and should never be added to this class!

class UT_API UT_TriangleMesh
{
public:
                         UT_TriangleMesh() = default;

    // Build a triangulation from a list of triangle points. The length
    // of `tri_pts` must be a multiple of 3 and should represent a manifold
    // surface. Edges appearing to be shared by more than two triangles are
    // automatically treated as cut-edges. Triangle edges with the same
    // endpoints are identified as paired half-edges, unless forced apart
    // by passing `boundary_edge_pts`, which is an array of even length broken
    // into pairs of consecutive point numbers each representing an edge of a
    // triangle that should not be identified with any other triangle edges
    // and therefore forced to be a boundary edge.
    //
    // The created triangle mesh will have as many triangles as indicated by
    // tri_pts (length of tri_pts divided by 3). The "adopted" number of mesh
    // points will be one plus the largest integer in the tri_pts array
    // although any number in the range 0 to numPoints() - 1 that does not
    // appear in tri_pts, has no manifestation in the actual mesh. In other
    // words there can be unused points if tri_pts array does not cover a
    // contiguous list starting from 0.
    //
    // If `reindex_points` is set to true, the points are re-indexed to to the
    // range from 0 to number of distinct entries in `tri_pts` minus one. The
    // index of a point will agree with its rank `tri_pts` viewed as a set:
    // the smallest number in `tri_pts` will have rank 0, second smallest 1,
    // etc.
    //
    // If `hedge_pair_class` is passed, it should have the same number of
    // entries as `tri_pts`, each corresponding to the half-edge implicitly
    // defined by that array in each position. Two half-edges are then allowed
    // to be paired as sym() to each other only if they have matching
    // hedge_pair_classes. For example, if the triangles are coming from
    // triangulation of larger polygons, putting original hedges and those
    // created by triangulation in different classes prevents them from being
    // identified even when endpoints match.

    void                 build(const UT_IntArray &tri_pts,
                                const UT_IntArray *boundary_edge_pts = nullptr,
                                bool reindex_points = false,
                                const UT_IntArray *hedge_pair_class = nullptr);

    SYS_FORCE_INLINE
    int                  numTriangles() const
                            { return int(myVertexPoint.size()) / 3; }

    SYS_FORCE_INLINE
    int                  numPoints() const
                            { return int(myPointVertex.size()); }

    SYS_FORCE_INLINE
    int                  numHedges() const { return int(myVertexPoint.size()); }

    SYS_FORCE_INLINE
    int                  numVertices() const { return numHedges(); }

    // The jth hedge of triangle t.
    SYS_FORCE_INLINE
    int                  triangleHedge(int t, int j = 0) const
                            { return 3 * t + j % 3; }

    // The jth vertex of triangle t (same as jth hedge).
    SYS_FORCE_INLINE
    int                  triangleVertex(int t, int j = 0) const
                            { return triangleHedge(t, j); }

    // Triangle of hedge h.
    SYS_FORCE_INLINE
    int                  hedgeTriangle(int h) const { return h / 3; }

    // Triangle of vertex v.
    SYS_FORCE_INLINE
    int                  vertexTriangle(int v) const
                            { return hedgeTriangle(v); }

    // The jth point of triangle t.
    SYS_FORCE_INLINE
    int                  trianglePoint(int t, int j = 0) const
                            { return myVertexPoint(3 * t + j % 3); }

    SYS_FORCE_INLINE
    int                  srcVertex(int h) const { return h; }

    SYS_FORCE_INLINE
    int                  dstVertex(int h) const
                            { return (h % 3 == 2) ? h - 2 : h + 1; }

    SYS_FORCE_INLINE
    int                  apxVertex(int h) const
                            { return (h % 3 == 0) ? h + 2 : h - 1; }

    SYS_FORCE_INLINE
    int                  onext(int h) const { return sym(lprev(h)); }

    SYS_FORCE_INLINE
    int                  vertexPoint(int v) const
                            { return myVertexPoint(v); }

    SYS_FORCE_INLINE
    int                  srcPoint(int h) const
                            { return myVertexPoint(srcVertex(h)); }

    SYS_FORCE_INLINE
    int                  dstPoint(int h) const
                            { return myVertexPoint(dstVertex(h)); }

    SYS_FORCE_INLINE
    int                  apxPoint(int h) const
                            { return myVertexPoint(apxVertex(h)); }

    SYS_FORCE_INLINE
    int                  lnext(int h) const { return dstVertex(h); }

    SYS_FORCE_INLINE
    int                  lprev(int h) const { return apxVertex(h); }

    SYS_FORCE_INLINE
    int                  sym(int h) const { return myHedgeSym(h); }

    SYS_FORCE_INLINE
    int                  pointVertex(int p) const { return myPointVertex(p); }

    SYS_FORCE_INLINE
    int                  pointHedge(int p) const { return myPointVertex(p); }

    SYS_FORCE_INLINE
    bool                 isBoundaryHedge(int h) const { return sym(h) < 0; }

    SYS_FORCE_INLINE
    bool                 isBoundaryPoint(int p) const
                            { return isBoundaryHedge(pointHedge(p)); }

    // Number of boundary components
    SYS_FORCE_INLINE
    int                  numBoundaries() const
                            { return int(myBoundaryRepHedges.size()); }

    // Return a first (representative) hedge for component comp.
    SYS_FORCE_INLINE
    int                  boundaryHedge(int comp) const
                            { return myBoundaryRepHedges(comp); }

    // Returns the boundary component number to which h belongs, or -1 if h
    // is an interior hedge.
    SYS_FORCE_INLINE
    int                  hedgeBoundary(int h) const
                            { return -SYSmin(0, sym(h)) - 1; }

    // Get all the boundary component hedges in order.
    void                 traceBoundary(int h0, UT_IntArray &bd) const;

    int                  numInteriorHedges() const
                            { return myNumInteriorHedges; }

    int                  numBoundaryHedges() const
                            { return numHedges() - numInteriorHedges(); }

    int                  numBoundaryPoints() const
                            { return numBoundaryHedges(); }

    int                  numInteriorPoints() const
                            { return numPoints() - numBoundaryPoints(); }

    int                  numEdges() const
                            { return numBoundaryHedges()
                                        + numInteriorHedges() / 2; }

    int                  eulerCharacteristic() const
                            { return numPoints() + numTriangles()
                                        - numEdges(); }

private:
    UT_IntArray          myVertexPoint; // Points of triangles (3 per tri)
    UT_IntArray          myHedgeSym;    // Syms of triangle hedges (3 per tri)
    UT_IntArray          myPointVertex; // First hedge of each point

    // One half-edge representing each boundary component.
    UT_IntArray          myBoundaryRepHedges;
    int                  myNumInteriorHedges                    = 0;
};

// An embedded triangle mesh will support a point position for each point
// in the mesh. The entire implementation is left to the template parameter.
// This includes the return type of the point position.

template <typename E>
class UT_EmbeddedTriangleMesh : public UT_TriangleMesh
{
public:

    using PointCoords = typename E::value_type;
    using Real = typename PointCoords::value_type;
    using Complex = typename std::complex<Real>;
    using Vector3 = UT_Vector3T<Real>;

    explicit             UT_EmbeddedTriangleMesh() = default;

    SYS_FORCE_INLINE
    auto                 pointPosition(int pt) const
                            { return myEmbedding->pointPosition(pt); }

    void                 setMeshPointPositions(const E *ptpos)
                            { myEmbedding = ptpos; }

    bool                 hasEmbedding() const
                            { return myEmbedding != nullptr; }

    // Returns the barycenter of t.
    SYS_FORCE_INLINE
    auto                 triangleBarycenter(int t) const;

    // Returns the normal of the triangle t. If normalize is set to false,
    // the normal vector will have magnitude equal to twice the area of t.
    //
    // NB. The normals used here are based on COUNTERCLOCKWISE orientation
    // of triangles UNLIKE Houdini polygon normals which are backwards. Based
    // on what these meshes are intended for, it's sheer asininity to adopt
    // Houdini's convention here and people in disagreement should take/retake
    // a bunch of math classes.

    SYS_FORCE_INLINE
    auto                 triangleNormal(int t, bool normalize = true) const;

    // Returns the area of triangle t.
    SYS_FORCE_INLINE
    auto                 triangleArea(int t) const
                            { return 0.5 * triangleNormal(t, false).length(); }

    // Returns the (integral of) gradient of the linear function defined on a
    // triangle t by values f0, f1, and f2 on its points 0, 1, and 2,
    // respectively, over t.
    template <typename T>
    auto                 triangleGradient(int t, T f0, T f1, T f2) const;

    // Returns the (integral of) gradient of the function defined on dual
    // vertices (triangle centroids) over the dual cell of a point p.
    template <typename F, typename A>
    Vector3              pointGradient(int p,
                                const F &tri_fn,
                                const A &pt_atlas) const;

    // Calculates the (integral of) Laplacian of a function defined on dual
    // vertices over triangle t. The function value at triangles is supplied by
    // tri_fn.
    template <typename F>
    auto                 triangleLaplacian(int t, const F &tri_fn) const;

    // Calculates the (integral of) Laplacian of a function defined on dual
    // vertices over the dual cell of point p. The function value at a point
    // is supplied by pt_fn.
    template <typename F>
    auto                 pointLaplacian(int p, const F &pt_fn) const;

    // Returns the vector from the src to dst of half-edge h.
    SYS_FORCE_INLINE
    auto                 hedgeVector(int h) const;

    // Returns the length of a half-edge h.
    SYS_FORCE_INLINE
    auto                 hedgeLength(int h) const
                            { return hedgeVector(h).length(); }

    // The cotan-based estimated length of the dual to h.
    SYS_FORCE_INLINE
    auto                 dualHedgeLength(int h) const;

    // Returns the dihedral angle between the two triangles incident to h.
    // This is angles between the normals with a positive sign if the edge
    // convex and negative otherwise.
    //
    // NB. These dihedral angles are calculating based on counterclockwise
    // triangle normals. So, the sign must be flipped if the mesh is built
    // consistently with a Houdini polygon geometry to give the expected
    // dihedral angle. (See notes above `triangleNormal()`).
    SYS_FORCE_INLINE
    auto                 hedgeDihedralAngle(int h) const;

    // Returns the cotan of the opposite angle to h.
    SYS_FORCE_INLINE
    auto                 hedgeCotan(int h) const;

    // Returns the angle (in radians) of the angle opposite to h.
    SYS_FORCE_INLINE
    auto                 hedgeAngle(int h) const
                            { return SYSatan2(Real(1.0), hedgeCotan(h)); }

    // Returns the dual area of a point which is the a third of the sum of
    // areas of triangles incident to it.
    SYS_FORCE_INLINE
    auto                 pointDualArea(int p) const;

    // Returns the sum of angles of vertices incident to p.
    SYS_FORCE_INLINE
    auto                 pointAngleSum(int p) const;

    // Returns a unit normal for a point p. By default the point normal is the
    // normalized angle-weighted sum of the normals of the incident triangles.
    // If `area_weighted` is set to true, incident triangle areas are used
    // instead of angles.
    //
    // NB. Like triangle normals, point normals are based on counterclockwise
    // winding of triangles. See the notes above `triangleNormal()`.
    SYS_FORCE_INLINE
    auto                 pointNormal(int p, bool area_weighted = false) const;

    // Returns the integral of the gaussian curvature in the mesh dual cell
    // of a point p.
    SYS_FORCE_INLINE
    auto                 pointGaussCurvature(int p) const;

    // Returns the integral of the mean curvature in the mesh dual cell
    // of a point p.
    SYS_FORCE_INLINE
    auto                 pointMeanCurvature(int p) const;

    // Returns the integral of the gaussian curvature over a mesh triangle.
    // Note that the curvature is calculated based on the Gaussian curvature
    // of the incident vertices and therefore depends on all the angles
    // incident to three corners of the tringle.
    SYS_FORCE_INLINE
    auto                 triangleGaussCurvature(int t) const;

    // Returns the integral of the mean curvature in the mesh dual cell
    // of a point p.
    SYS_FORCE_INLINE
    auto                 triangleMeanCurvature(int t) const;

    // Returns maximum principal direction as a complex number. The complex
    // plane based on which the number should be interpreted is the plane
    // of the triangle with the triangles first half-edge as the real axis.
    SYS_FORCE_INLINE
    Complex              triangleMaxPrincipalDirection(int t) const;


    // Given the polar coordinates of a tangent vector in the exponential
    // mapping of the star of point p, generate an extrinsic (vector 3)
    // tangent vector to the embedding of the mesh.
    Vector3              extrinsicTriangleTangent(int p, Real radial,
                                Real angular) const;

    // Same as the abvoe method, using a complex number.
    SYS_FORCE_INLINE
    Vector3              extrinsicTriangleTangent(int t, Complex z) const;

private:
    const E             *myEmbedding                            = nullptr;
};

template <typename E>
auto
UT_EmbeddedTriangleMesh<E>::triangleBarycenter(int t) const
{
    return (pointPosition(trianglePoint(t, 0))
            + pointPosition(trianglePoint(t, 1))
            + pointPosition(trianglePoint(t, 2))) / 3.0;
}

template <typename E>
auto
UT_EmbeddedTriangleMesh<E>::triangleNormal(int t, bool normalize) const
{
    auto a = pointPosition(trianglePoint(t, 0));
    auto b = pointPosition(trianglePoint(t, 1));
    auto c = pointPosition(trianglePoint(t, 2));

    auto nml = cross(b - a, c - a);
    if (normalize)
        nml.normalize();

    return nml;
}

template <typename E>
template <typename T>
auto
UT_EmbeddedTriangleMesh<E>::triangleGradient(int t, T f0, T f1, T f2) const
{
    auto a = pointPosition(trianglePoint(t, 0));
    auto b = pointPosition(trianglePoint(t, 1));
    auto c = pointPosition(trianglePoint(t, 2));

    auto nml = cross(b - a, c - a);
    auto areax2 = nml.normalize();

    return (f2 * cross(nml, b - a) + f0 * cross(nml, c - b)
            + f1 * cross(nml, a - c)) / areax2;
}

template <typename E>
template <typename F>
auto
UT_EmbeddedTriangleMesh<E>::triangleLaplacian(int t, const F &tri_fn) const
{
    auto f = tri_fn(t);
    auto res = 0.0 * f;

    for (int i = 0; i < 3; i++)
    {
        auto h = triangleHedge(t, i);
        auto h_sym = sym(h);
        if (h_sym < 0)
            continue;

        auto t_sym = hedgeTriangle(h_sym);
        auto f_sym = tri_fn(t_sym);

#if 0
        auto wt = 1.0 / (hedgeCotan(h) + hedgeCotan(h_sym));
        res += (f_sym - f) * wt;
#else
        res += (f_sym - f);
#endif

    }
    return 0.5 * res;
}

template <typename E>
template <typename F>
auto
UT_EmbeddedTriangleMesh<E>::pointLaplacian(int p, const F &pt_fn) const
{
    auto f = pt_fn(p);

    auto res = 0.0 * f;
    auto ring_area = 0.0;

    auto h0 = pointHedge(p);
    auto h = h0;
    auto h_sym = -1;
    do
    {
        // auto cot = SYSmax(1e-7, hedgeCotan(h));
        auto cot = hedgeCotan(h);
        h_sym = sym(h);
        if (h_sym >= 0)
            // cot += SYSmax(1e-7, hedgeCotan(h_sym));
            cot += hedgeCotan(h_sym);


        auto f_dst = pt_fn(dstPoint(h));
        res += cot * (f_dst - f);
        ring_area += triangleArea(hedgeTriangle(h));
    } while (h = onext(h), h >= 0 && h != h0);

    if (h < 0)
    {
        h_sym = h_sym >= 0 ? lprev(sym(h_sym)) : lprev(h0);
        auto cot = hedgeCotan(h_sym);
        auto f_dst = pt_fn(srcPoint(h_sym));
        res += cot * (f_dst - f);
    }

    return 0.5 * res;
}

template <typename E>
template <typename F, typename A>
typename UT_EmbeddedTriangleMesh<E>::Vector3
UT_EmbeddedTriangleMesh<E>::pointGradient(int p, const F &tri_fn,
                                          const A &pt_atlas) const
{
    auto hedge_contribution = [&](int h)
    {
        auto i = Complex(0.0, 1.0);
        auto t = hedgeTriangle(h);
        auto s = tri_fn(t);
        auto c = hedgeCotan(h);
        auto h_z = pt_atlas.hedgeComplexCoord(h);
        auto z = s * c * h_z * i;

        auto h_sym = sym(h);
        if (h_sym >= 0)
        {
            auto t_sym = hedgeTriangle(h_sym);
            auto s_sym = tri_fn(t_sym);
            z -= s_sym * c * h_z * i;
        }

        return z;
    };

    Complex z(0, 0);
    auto h0 = pointHedge(p);
    auto h = h0;
    do
    {
        z += hedge_contribution(h);
    } while (h = onext(h), h >= 0 && h != h0);

    return pt_atlas.extrinsicPointTangent(p, z) / pointDualArea(p);
}

template <typename E>
auto
UT_EmbeddedTriangleMesh<E>::hedgeVector(int h) const
{
    return (pointPosition(dstPoint(h)) - pointPosition(srcPoint(h)));
}

template <typename E>
auto
UT_EmbeddedTriangleMesh<E>::dualHedgeLength(int h) const
{
    auto cotsum = hedgeCotan(h);
    auto h_sym = sym(h);
    if (h_sym >= 0)
        cotsum += hedgeCotan(h_sym);

    return 0.5 * hedgeLength(h) * cotsum;
}

template <typename E>
auto
UT_EmbeddedTriangleMesh<E>::hedgeDihedralAngle(int h) const
{
    auto h_sym = sym(h);
    if (h_sym < 0)
        return Real(0.0);

    auto n0 = triangleNormal(hedgeTriangle(h));
    auto n1 = triangleNormal(hedgeTriangle(h_sym));
    auto e = hedgeVector(h);
    e.normalize();
    return SYSatan(dot(e, cross(n0, n1)), dot(n0, n1));
}

template <typename E>
auto
UT_EmbeddedTriangleMesh<E>::hedgeCotan(int h) const
{
    auto a = pointPosition(dstPoint(lnext(h)));
    auto b = pointPosition(srcPoint(h));
    auto c = pointPosition(dstPoint(h));

    auto u = b - a;
    auto v = c - a;

    return dot(u, v) / cross(u, v).length();
}

template <typename E>
auto
UT_EmbeddedTriangleMesh<E>::pointDualArea(int p) const
{
    auto h = pointHedge(p);
    auto h0 = h;
    auto sum = triangleArea(hedgeTriangle(h));

    while (h = onext(h), h >= 0 && h != h0)
        sum += triangleArea(hedgeTriangle(h));

    return sum / Real(3.0);
}

template <typename E>
auto
UT_EmbeddedTriangleMesh<E>::pointAngleSum(int p) const
{
    auto h = pointHedge(p);
    auto h0 = h;
    auto sum = hedgeAngle(lnext(h));

    while (h = onext(h), h >= 0 && h != h0)
        sum += hedgeAngle(lnext(h));

    return sum;
}

template <typename E>
auto
UT_EmbeddedTriangleMesh<E>::pointNormal(int p, bool area_weighted) const
{
    auto h = pointHedge(p);
    auto h0 = h;

    if (area_weighted)
    {
        auto nml = triangleNormal(hedgeTriangle(h), false);
        while (h = onext(h), h >= 0 && h != h0)
            nml += triangleNormal(hedgeTriangle(h), false);

        nml.normalize();
        return nml;
    }

    auto nml = triangleNormal(hedgeTriangle(h)) * hedgeAngle(lnext(h));
    while (h = onext(h), h >= 0 && h != h0)
        nml += triangleNormal(hedgeTriangle(h)) * hedgeAngle(lnext(h));

    nml.normalize();
    return nml;
}

template <typename E>
auto
UT_EmbeddedTriangleMesh<E>::pointGaussCurvature(int p) const
{
    return (isBoundaryPoint(p) ? Real(M_PI) : Real(2.0 * M_PI))
           - pointAngleSum(p);
}

template <typename E>
auto
UT_EmbeddedTriangleMesh<E>::pointMeanCurvature(int p) const
{
    auto h0 = pointHedge(p);
    auto h = h0;
    auto b = Real(0.0);
    do
    {
        b += hedgeLength(h) * hedgeDihedralAngle(h);
    } while (h = onext(h), h >= 0 && h != h0);
    return Real(0.25) * b;
}

template <typename E>
auto
UT_EmbeddedTriangleMesh<E>::triangleGaussCurvature(int t) const
{
    auto K = Real(-M_PI);
    for (int i = 0; i < 3; i++)
    {
        int h = triangleHedge(t, i);
        int p = apxPoint(h);
        auto a = hedgeAngle(h);
        K += (isBoundaryPoint(p) ? 1 : 2) * M_PI * a / pointAngleSum(p);
    }

    return K;
}

template <typename E>
auto
UT_EmbeddedTriangleMesh<E>::triangleMeanCurvature(int t) const
{
    auto H = Real(0);
    for (int i = 0; i < 3; i++)
    {
        auto h = triangleHedge(t, i);
        H += 0.25 * hedgeLength(h) * hedgeDihedralAngle(h);
    }
    return H;
}

template<typename E>
typename UT_EmbeddedTriangleMesh<E>::Complex
UT_EmbeddedTriangleMesh<E>::triangleMaxPrincipalDirection(int t) const
{
    Complex z(0, 0);
    auto partial_angle_sum = Real(0.0);

    for (int i = 0; i < 3; i++)
    {
        auto h = triangleHedge(t, i);
        auto l = hedgeLength(h);
        auto alpha = hedgeDihedralAngle(h);

        auto theta = partial_angle_sum;

        // Multiply by theta by 2.0 to go from vectors to directions.
        Complex r(SYScos(2.0 * theta), SYSsin(2.0 * theta));

        z += l * alpha * r;
        partial_angle_sum += (M_PI - hedgeAngle(lprev(h)));
    }

    return std::sqrt(z);
}

template<typename M>
typename UT_EmbeddedTriangleMesh<M>::Vector3
UT_EmbeddedTriangleMesh<M>::extrinsicTriangleTangent(int p, Complex z) const
{
    return extrinsicTriangleTangent(p, std::abs(z), std::arg(z));
}

template<typename M>
typename UT_EmbeddedTriangleMesh<M>::Vector3
UT_EmbeddedTriangleMesh<M>::extrinsicTriangleTangent(int t, Real radial,
                                                     Real angular) const
{
    auto n = triangleNormal(t);
    auto x = hedgeVector(triangleHedge(t, 0));
    x.normalize();
    x *= radial;
    auto Jx = cross(n, x);
    return x * SYScos(angular) + Jx * SYSsin(angular);
}


// UT_PointAtlas is constructs a local polar coordinate system around each
// point of an embedded triangle mesh to support computations on a complex
// plane corresponding to (the exponential map of) the ring of each point
// in which the real axis correspond to the half-edge returned by
// `pointHedge()` method of the mesh.

template <typename M>
class UT_PointAtlas
{
public:
    using Real = typename M::Real;
    using Complex = typename M::Complex;
    using Vector3 = typename M::Vector3;

    explicit             UT_PointAtlas(const M &mesh);


    // Returns the exponential map scale factor at p which is 2π divided by the
    // sum of the mesh angles incident to p.
    SYS_FORCE_INLINE
    Real                 pointAngleScale(int p) const
                            { return myPointScale(p); }

    // Returns the *raw* angle sum between h0 = mesh.pointHedge(p), where p
    // is the source point of h, and h in a counterclockwise scale (using
    // mesh.onext()).
    SYS_FORCE_INLINE
    Real                 hedgeAngularOffset(int h) const;

    // Returns the angular component of the polar coordinates of the vector
    // represented by h in the exponential mapping of star of the source point
    // of h onto a planar patch. This is the angularOffset of h scaled by point
    // angle scale at the source of h. The output will be in range [0, 2π).
    SYS_FORCE_INLINE
    Real                 hedgeAngularCoord(int h) const
                            { return pointAngleScale(myMesh.srcPoint(h))
                                        * hedgeAngularOffset(h); }

    // Returns the radial component of the polar coordinates of the vector
    // represented by h (or equivalently the destination point of h) in the
    // exponential mapping of star of the source point of h onto a planar
    // patch. This is simply the length of h!
    SYS_FORCE_INLINE
    Real                 radialCoord(int h) const
                            { return myMesh.hedgeLength(h); }

    // Returns the polar coordinates of the vector represented by h in the
    // exponential mapping of the star of the source point of h as a complex
    // number with the same polar coordinates.
    SYS_FORCE_INLINE
    Complex              hedgeComplexCoord(int h) const;

    SYS_FORCE_INLINE
    Complex              hedgeComplexDirection(int h, int n = 1) const;

    // Returns the the angular polar component of the complex number z. This
    // is the phase of z, implemented using std::arg(z), in range [0, 2π).
    SYS_FORCE_INLINE
    Real                 complexToAngular(Complex z) const;

    // Returns the radial component of the complex number z.
    SYS_FORCE_INLINE
    Real                 complexToRadial(Complex z) const
                            { return std::abs(z); }

    // Returns the outgoing half-edge at p with the largest angular coordinate
    // less than or equal to `angle`.
    //
    // Note: This is done using a linear scan. We could be brought down to
    // O(log(valence(p))) if we maintained a random-access list of incident
    // hedge angular values. Since we like to be able to lookup hedge angular
    // coordinates in constant time, this would require de facto repeated
    // storage of the values with a rearrangement. Given that average point
    // valence in a triangle mesh is 6, this would be excessive.
    int                  angularLowerBoundHedge(int p, Real angle) const;

    // Given the polar coordinates of a tangent vector in the exponential
    // mapping of the star of point p, generate an extrinsic (vector 3)
    // tangent vector to the embedding of the mesh.
    Vector3              extrinsicPointTangent(int p, Real radial,
                                Real angular) const;

    // Same as the abvoe method, using a complex number.
    SYS_FORCE_INLINE
    Vector3              extrinsicPointTangent(int p, Complex z) const;

    SYS_FORCE_INLINE
    Complex              pointMaxPrincipalDirection(int p) const;

private:
    const M             &myMesh;
    UT_Array<Real>       myPointScale;
    UT_Array<Real>       myHedgeAngle;
};

template<typename M>
UT_PointAtlas<M>::UT_PointAtlas(const M &mesh) :
        myMesh(mesh)
{
    myPointScale.setSizeNoInit(mesh.numPoints());
    myHedgeAngle.setSizeNoInit(mesh.numHedges());

    UTparallelFor(UT_BlockedRange<int>(0, mesh.numPoints()),
    [&](const UT_BlockedRange<int> &range)
    {
        for (int p = range.begin(), pe = range.end(); p != pe; p++)
        {
            auto h0 = mesh.pointHedge(p);
            auto h = h0;
            Real angle_sum = Real(0);
            do
            {
                // Note that myHedgeAngle(h) carries the sum of source angles
                // from h0 to h (inclusive). The last of these should be
                // excluded when calculating the angle coordinate of h.
                // We include it here in the sum because this is the only
                // place to keep track that angle semantically (think of
                // boundary points).
                angle_sum += mesh.hedgeAngle(mesh.lnext(h));
                myHedgeAngle(h) = angle_sum;
            } while (h = mesh.onext(h), h >= 0 && h != h0);
            myPointScale(p) = 2 * M_PI / angle_sum;
        }
    });
}

template<typename M>
typename UT_PointAtlas<M>::Real
UT_PointAtlas<M>::hedgeAngularOffset(int h) const
{
    if (h == myMesh.pointHedge(myMesh.srcPoint(h)))
        return 0.0;

    return myHedgeAngle(myMesh.lnext(myMesh.sym(h)));
}

template<typename M>
typename UT_PointAtlas<M>::Complex
UT_PointAtlas<M>::hedgeComplexCoord(int h) const
{
    auto a = hedgeAngularCoord(h);
    return radialCoord(h) * Complex(SYScos(a), SYSsin(a));
}

template<typename M>
typename UT_PointAtlas<M>::Complex
UT_PointAtlas<M>::hedgeComplexDirection(int h, int n) const
{
    auto a = hedgeAngularCoord(h);
    return Complex(SYScos(Real(n) * a), SYSsin(Real(n) * a));
}

template<typename M>
typename UT_PointAtlas<M>::Real
UT_PointAtlas<M>::complexToAngular(Complex z) const
{
    auto t = std::arg(z);
    return t - 2.0 * M_PI * SYSfloor(t / (2.0 * M_PI));
}

template<typename M>
int
UT_PointAtlas<M>::angularLowerBoundHedge(int p, Real angle) const
{
    // This can in theory be made to run in time logarithmic in valence
    // of p; but this requires random access into the orderd list incident
    // outgoing hedges at p. Unless we want to do lots of these queries on
    // very irregular meshes this would be excessive.

    angle = angle - 2.0 * M_PI * SYSfloor(angle / (2.0 * M_PI));
    angle /= pointAngleScale(p);

    auto h0 = myMesh.pointHedge(p);
    auto h = h0;
    auto last_h = h0;
    do
    {
        if (myHedgeAngle(h) > angle)
            return h;
        last_h = h;
    } while (h = myMesh.onext(h), h >= 0 && h != h0);

    // We have zoned the angle above to [0..2π) range, therefore getting
    // here means that when dividing by pointAngleScale(p) our corrected
    // angle ended up going over the angle sum at p which could only be
    // due to numerical error. So, we just drop down to the last outgoing
    // half-edge before getting back to h0.
    return last_h;
}

template<typename M>
typename UT_PointAtlas<M>::Vector3
UT_PointAtlas<M>::extrinsicPointTangent(int p, Real radial, Real angular) const
{
    auto h = angularLowerBoundHedge(p, angular);
    auto n = myMesh.triangleNormal(myMesh.hedgeTriangle(h));
    auto a = (angular - hedgeAngularCoord(h)) / pointAngleScale(p);

    auto v = myMesh.hedgeVector(h);
    v.normalize();
    v *= radial;

    if (a > 0.0)
    {
        UT_QuaternionR q;
        q.updateFromAngleAxis(a, n);
        v = q.rotate(v);
    }
    return v;
}

template<typename M>
typename UT_PointAtlas<M>::Vector3
UT_PointAtlas<M>::extrinsicPointTangent(int p, Complex z) const
{
    return extrinsicPointTangent(p, complexToRadial(z), complexToAngular(z));
}

template<typename M>
typename UT_PointAtlas<M>::Complex
UT_PointAtlas<M>::pointMaxPrincipalDirection(int p) const
{
    Complex z(0, 0);
    auto h0 = myMesh.pointHedge(p);
    auto h = h0;
    do
    {
        auto l = myMesh.hedgeLength(h);
        auto alpha = myMesh.hedgeDihedralAngle(h);
        z += l * alpha * hedgeComplexDirection(h, 2);
    } while (h = myMesh.onext(h), h >= 0 && h != h0);

    return std::sqrt(z);
}

#endif