cdt.hpp

Go to the documentation of this file.

 
#pragma once
#include "meshkit/detail/storage.hpp"
#include "meshkit/terrainpolicy.hpp"
#include "meshkit/progress.hpp"
#include <functional>
#include <limits>
 
namespace meshkit
{
 
template <class Storage> class BasicConstrainedDelaunay
{
public:
    using Point = typename Storage::Point;
    using Box = typename Storage::Box;
    using String = typename Storage::String;
    using PointInput = typename Storage::PointInput;
    template <class T> using Vector = typename Storage::template Vector<T>;
    template <class K, class V> using Hash = typename Storage::template Hash<K, V>;
    template <class K, class V> using MultiHash = typename Storage::template MultiHash<K, V>;
    template <class K> using Set = typename Storage::template Set<K>;
    struct Triangle {
        int v[3];   
        int adj[3]; 
        bool constrained[3]; 
        bool alive;
    };
 
    [[nodiscard]] static double orientExact(const Point &a, const Point &b, const Point &c);
 
    bool build(const PointInput &points, Vector<int> *vertexOfPoint = nullptr,
               const std::function<bool()> &cancelled = {},
               const std::function<void(std::uint64_t)> &progress = {});
 
    bool insertConstraint(int a, int b);
 
    int insertPoint(const Point &p);
 
    bool removeVertex(int v);
    [[nodiscard]] bool isVertexActive(int v) const
    {
        return v >= 0 && v < m_vertexTri.size() && m_vertexTri[v] >= 0;
    }
 
    void removeExterior();
    void removeSuperTriangles();
    int removeRegionLeftOf(int a, int b);
    [[nodiscard]] int inputVertexCount() const { return m_superBase; }
 
    int removeRegionAt(const Point &p); 
 
    struct QualityOptions {
        std::function<double(double, double)> hAt;
        std::function<double(const Point *)> maximumAreaForTriangle;
        double maximumAdjacentRatio = 0;
        std::function<double(int, int)> fixedNeighbourSize;
        int sizingRepairPasses = 0, maximumSizingRepairInsertions = 200000;
        std::function<double(const Point *, const double *, Point *)> terrainError;
        std::function<double(double, double)> terrainElevationAt;
        double terrainTolerance = 0.0;
        std::function<double(const Box &)>
            terrainToleranceForFootprint; 
 
        double terrainMinSpacing = 0.0;
        bool terrainWorstFirst = false;
        bool terrainFinalCheck = true;
        int smoothingPasses = 0;
        int coarsenPasses = 0; 
        CoarseningOrder coarseningOrder = CoarseningOrder::LegacyFlatness;
        std::function<double(const Box &)>
            coarseningTargetSize; 
        double coarseningReferenceSize = 10, normalVerticalToXY = 1;
 
        std::function<bool(const Point &)> protectedLocation;
        Set<int> protectedVertices; 
        std::function<double(const Point *, const double *)> terrainErrorBound;
        bool prioritizeQuality = true;
        int maxTriangles = std::numeric_limits<int>::max();
        double minAngleDeg = 30.0;
        int maxInsertions = std::numeric_limits<int>::max();
        double minEdge = 0.0;
        std::function<bool()> cancelled;
        std::function<void(const GenerationProgress &)> onProgress;
    };
    struct QualityReport {
        int sizingRepairInserted = 0, sizingRepairPasses = 0;
        bool sizingRepairLimited = false;
        int inserted = 0; 
        int segmentSplits = 0;
        int blockedByFixed =
            0; 
        int sizeInsertions = 0, qualityInsertions = 0, terrainInsertions = 0;
        int terrainUnresolved = 0, terrainUnknown = 0;
        double maxTerrainError = 0.0;
        bool capped = false; 
        bool cancelled = false;
        int coarsenRemoved = 0, coarsenRejectedAngle = 0, coarsenRejectedSize = 0;
        int coarsenRejectedTerrain = 0, coarsenRejectedEdge = 0, coarsenAttempts = 0;
        int coarsenHeightRejects = 0, coarsenBoundAccepts = 0, coarsenExactQueries = 0;
        bool coarsenSkippedTerrain = false;
        double coarsenSeconds = 0;
        std::uint64_t coarsenQueuePeak = 0, coarsenPriorityQueries = 0;
        bool invalidPolicy = false;
    };
    QualityReport refineQuality(const QualityOptions &opt);
    void setFixedConstraint(int a, int b);
    [[nodiscard]] bool smallAngleExempt(int t) const;
    [[nodiscard]] bool touchesFixed(int t) const;
 
    // Optional diagnostic provenance. Bits may overlap at coincident inputs;
    // endpoint associations are not a causal explanation of final cell counts.
    enum Origin : std::uint32_t {
        Domain = 1,
        Hole = 2,
        Line = 4,
        Patch = 8,
        Node = 16,
        Breakline = 32,
        Region = 64,
        Lattice = 128,
        Size = 256,
        Angle = 512,
        Terrain = 1024,
        SegmentSplit = 2048
    };
    void enableProvenance() { m_origins.fill(0, m_pts.size()); }
    void addVertexOrigin(int v, std::uint32_t origin)
    {
        if (!m_origins.isEmpty() && v >= 0 && v < m_origins.size())
            m_origins[v] |= origin;
    }
    [[nodiscard]] const Vector<std::uint32_t> &vertexOrigins() const { return m_origins; }
 
    // ── Read-out ───────────────────────────────────────────────────────
    [[nodiscard]] const Vector<Point> &vertices() const { return m_pts; }
    [[nodiscard]] const Vector<Triangle> &triangles() const { return m_tris; }
    [[nodiscard]] const Vector<double> &terrainElevations() const { return m_terrainElevations; }
    [[nodiscard]] int liveTriangleCount() const;
    [[nodiscard]] int locate(const Point &p) const;
    [[nodiscard]] bool isConstrained(int a, int b) const;
    [[nodiscard]] Vector<int> constrainedChain(int a, int b) const;
    [[nodiscard]] bool isSuperVertex(int v) const
    {
        return v >= m_superBase && v < m_superBase + 3;
    }
    [[nodiscard]] String errorMsg() const { return m_errorMsg; }
 
private:
    struct Removal {
        Vector<int> star, ring;
        Vector<Triangle> replacement;
    };
    bool prepareRemoval(int v, Removal *out) const;
    void commitRemoval(int v, const Removal &removal);
    double orient(int a, int b, int c) const;
    double orient(int a, int b, const Point &p) const;
    bool inCircle(int a, int b, int c, int d) const;
    int findTriangle(const Point &p, int start) const;
    int triangleWithEdge(int a, int b, int *edgeIndexOut) const;
    void legalize(int t, int e);
    void flip(int t, int e);
    void setAdj(int t, int e, int n);
    void splitTriangle(int t, int v);
    void splitEdge(int t, int e, int v);
    int edgeIndex(int t, int v) const;
    void trianglesAround(int v, Vector<int> *out) const;
    void markConstrained(int a, int b);
    bool insertConstraintImpl(int a, int b, std::uint64_t origin);
    int splitSubsegment(int a, int b);
    bool subsegmentEncroached(int a, int b) const;
    bool isFixedSub(int a, int b) const;
    bool exemptShortestEdge(int u, int v) const;
 
    Vector<Point> m_pts;
    Vector<std::uint32_t> m_origins;
    Vector<double> m_terrainElevations;
    Vector<Triangle> m_tris;
    int m_liveTriangles = 0;
    Vector<int> m_vertexTri; 
    int m_superBase = 0;
    int m_lastLocate = 0;
    String m_errorMsg;
    Hash<std::uint64_t, std::uint64_t>
        m_segOrigin; 
    MultiHash<std::uint64_t, std::uint64_t>
        m_segOriginExtra; 
    Hash<std::uint64_t, std::uint64_t>
        m_segPiece; 
    Hash<int, std::uint64_t>
        m_vertexSeg;               
    Set<std::uint64_t> m_fixedSub; 
    // Installed only during refinement. Every geometric edit invalidates the
    // affected work, including flips outside the inserted vertex's final star.
    std::function<void(int)> m_triangleChanged;
};
 
using ConstrainedDelaunay = BasicConstrainedDelaunay<detail::StandardStorage>;
} // namespace meshkit