You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1771 lines
67 KiB
1771 lines
67 KiB
// David Eberly, Geometric Tools, Redmond WA 98052
|
|
// Copyright (c) 1998-2021
|
|
// Distributed under the Boost Software License, Version 1.0.
|
|
// https://www.boost.org/LICENSE_1_0.txt
|
|
// https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
|
|
// Version: 4.8.2021.04.22
|
|
|
|
#pragma once
|
|
|
|
// Incremental insertion and removal of vertices in a Delaunay triangulation.
|
|
// The triangles are counterclockwise ordered. The insertion code is that of
|
|
// template <typename T> class Delaunay2<T> in Delaunay2.h. For now, the code
|
|
// was copy-pasted here. Any changes made to this code must be made in both
|
|
// places. The duplication is
|
|
// using Triangle = VETManifoldMesh::Triangle;
|
|
// using DirectedEdgeKeySet = std::set<EdgeKey<true>>;
|
|
// using TrianglePtrSet = std::set<std::shared_ptr<Triangle>>;
|
|
// VETManifoldMesh mGraph;
|
|
// std::array<std::array<size_t, 2>, 3> const mIndex;
|
|
// bool GetContainingTriangle(size_t, std::shared_ptr<Triangle>&) const;
|
|
// void GetAndRemoveInsertionPolygon(size_t,
|
|
// TrianglePtrSet&, DirectedEdgeKeySet&);
|
|
// void Update(size_t);
|
|
// static ComputeRational const& Copy(InputRational const&,
|
|
// ComputeRational&);
|
|
// int32_t ToLine(size_t, size_t, size_t) const;
|
|
// int32_t ToTriangle(size_t, size_t, size_t, size_t) const;
|
|
// int32_t ToCircumcircle(size_t, size_t, size_t, size_t) const;
|
|
//
|
|
// The removal code is an implementation of the algorithm in
|
|
// Olivier Devillers,
|
|
// "On Deletion in Delaunay Triangulations",
|
|
// International Journal of Computational Geometry and Applications,
|
|
// World Scientific Publishing, 2002, 12, pp. 193-205.
|
|
// https://hal.inria.fr/inria-00167201/document
|
|
// The weight function for the priority queue, implemented as a min-heap, is
|
|
// the negative of the function power(p,circle(q0,q1,q2)) function described
|
|
// in the paper.
|
|
//
|
|
// The paper appears to assume that the removal point is an interior point of
|
|
// the trianglation. Just as the insertion algorithms are different for
|
|
// interior points and for boundary points, the removal algorithms are
|
|
// different for interior points and for boundary points.
|
|
//
|
|
// The paper mentions that degeneracies (colinear points, cocircular points)
|
|
// are handled by jittering. Although one hopes that jittering prevents
|
|
// degeneracies--and perhaps probabilistically this is acceptable, the only
|
|
// guarantee for a correct result is to use exact arithmetic on the input
|
|
// points. The implementation here uses a blend of interval and rational
|
|
// arithmetic for exactness; the input points are not jittered.
|
|
//
|
|
// The details of the algorithms and implementation are provided in
|
|
// https://www.geometrictools.com/Documentation/IncrementalDelaunayTriangulation.pdf
|
|
|
|
#include <Mathematics/ArbitraryPrecision.h>
|
|
#include <Mathematics/MinHeap.h>
|
|
#include <Mathematics/SWInterval.h>
|
|
#include <Mathematics/Vector2.h>
|
|
#include <Mathematics/VETManifoldMesh.h>
|
|
#include <functional>
|
|
#include <map>
|
|
#include <set>
|
|
|
|
namespace gte
|
|
// The input type must be 'float' or 'double'. The compute type is defined
|
|
// internally and has enough bits of precision to handle any
|
|
// floating-point inputs.
|
|
{
|
|
template <typename T>
|
|
class IncrementalDelaunay2
|
|
{
|
|
public:
|
|
// Construction and destruction. A bounding rectangle for the input
|
|
// points must be specified.
|
|
IncrementalDelaunay2(T const& xMin, T const& yMin, T const& xMax, T const& yMax)
|
|
:
|
|
mXMin(xMin),
|
|
mYMin(yMin),
|
|
mXMax(xMax),
|
|
mYMax(yMax),
|
|
mCRPool(maxNumCRPool),
|
|
mGraph{},
|
|
mIndex{ { { 0, 1 }, { 1, 2 }, { 2, 0 } } },
|
|
mTriangles{},
|
|
mAdjacencies{},
|
|
mTrianglesAndAdjacenciesNeedUpdate(true),
|
|
mQueryPoint{},
|
|
mIRQueryPoint{}
|
|
{
|
|
static_assert(
|
|
std::is_floating_point<T>::value,
|
|
"Invalid floating-point type.");
|
|
|
|
LogAssert(
|
|
mXMin < mXMax && mYMin < mYMax,
|
|
"Invalid bounding rectangle.");
|
|
|
|
mToLineWrapper = [this](size_t vPrev, size_t vCurr, size_t vNext)
|
|
{
|
|
return ToLine(vPrev, vCurr, vNext);
|
|
};
|
|
|
|
// Create a supertriangle that contains the input rectangle.
|
|
T xDelta = mXMax - mXMin;
|
|
T yDelta = mYMax - mYMin;
|
|
T x0 = mXMin - xDelta;
|
|
T y0 = mYMin - yDelta;
|
|
T x1 = mXMin + static_cast<T>(5) * xDelta;
|
|
T y1 = y0;
|
|
T x2 = x0;
|
|
T y2 = mYMin + static_cast<T>(5) * yDelta;
|
|
Vector2<T> supervertex0{ x0, y0 };
|
|
Vector2<T> supervertex1{ x1, y1 };
|
|
Vector2<T> supervertex2{ x2, y2 };
|
|
|
|
// Insert the supertriangle vertices into the vertex storage.
|
|
mVertexIndexMap.emplace(supervertex0, 0);
|
|
mVertexIndexMap.emplace(supervertex1, 1);
|
|
mVertexIndexMap.emplace(supervertex2, 2);
|
|
mVertices.emplace_back(supervertex0);
|
|
mVertices.emplace_back(supervertex1);
|
|
mVertices.emplace_back(supervertex2);
|
|
mIRVertices.emplace_back(IRVector{ x0, y0 });
|
|
mIRVertices.emplace_back(IRVector{ x1, y1 });
|
|
mIRVertices.emplace_back(IRVector{ x2, y2 });
|
|
|
|
// Insert the supertriangle into the triangulation.
|
|
auto inserted = mGraph.Insert(0, 1, 2);
|
|
LogAssert(
|
|
inserted != nullptr,
|
|
"Failed to insert supertriangle.");
|
|
}
|
|
|
|
~IncrementalDelaunay2() = default;
|
|
|
|
// Insert a point into the triangulation. The return value is the
|
|
// index associated with the vertex in the vertex map. The
|
|
// supertriangle vertices are at indices 0, 1, and 2. If the input
|
|
// point already exists, its vertex-map index is simply returned. If
|
|
// the position is outside the domain specified in the constructor,
|
|
// an exception is thrown.
|
|
size_t Insert(Vector2<T> const& position)
|
|
{
|
|
mTrianglesAndAdjacenciesNeedUpdate = true;
|
|
|
|
LogAssert(
|
|
mXMin <= position[0] && position[0] <= mXMax &&
|
|
mYMin <= position[1] && position[1] <= mYMax,
|
|
"The position is outside the domain specified in the constructor.");
|
|
|
|
auto iter = mVertexIndexMap.find(position);
|
|
if (iter != mVertexIndexMap.end())
|
|
{
|
|
// The vertex already exists.
|
|
return iter->second;
|
|
}
|
|
|
|
// Store the position in the various pools.
|
|
size_t posIndex = mVertices.size();
|
|
mVertexIndexMap.emplace(position, posIndex);
|
|
mVertices.emplace_back(position);
|
|
mIRVertices.emplace_back(IRVector{ position[0], position[1] });
|
|
|
|
Update(posIndex);
|
|
return posIndex;
|
|
}
|
|
|
|
// Remove a point from the triangulation. The return value is the index
|
|
// associated with the vertex in the vertex map when that vertex exists.
|
|
// If the vertex does not exist, the return value is
|
|
// std::numeric_limit<size_t>::max().
|
|
size_t Remove(Vector2<T> const& position)
|
|
{
|
|
mTrianglesAndAdjacenciesNeedUpdate = true;
|
|
|
|
auto iter = mVertexIndexMap.find(position);
|
|
if (iter == mVertexIndexMap.end())
|
|
{
|
|
// The position is not a vertex of the triangulation.
|
|
return invalid;
|
|
}
|
|
int32_t vRemovalIndex = static_cast<int32_t>(iter->second);
|
|
|
|
if (mVertexIndexMap.size() == 4)
|
|
{
|
|
// Only a single point has been inserted previously into the
|
|
// triangulation.
|
|
for (int32_t i0 = 2, i1 = 0; i1 < 3; i0 = i1++)
|
|
{
|
|
auto removed = mGraph.Remove(vRemovalIndex, i0, i1);
|
|
LogAssert(
|
|
removed,
|
|
"Unexpected removal failure.");
|
|
}
|
|
|
|
auto inserted = mGraph.Insert(0, 1, 2);
|
|
LogAssert(
|
|
inserted != nullptr,
|
|
"Failed to insert supertriangle.");
|
|
|
|
mVertexIndexMap.erase(iter);
|
|
return static_cast<size_t>(vRemovalIndex);
|
|
}
|
|
|
|
// Locate the position in the vertices of the graph.
|
|
auto const& vMap = mGraph.GetVertices();
|
|
auto vIter = vMap.find(vRemovalIndex);
|
|
LogAssert(
|
|
vIter != vMap.end(),
|
|
"Expecting to find the to-be-removed vertex in the triangulation.");
|
|
|
|
bool removalPointOnBoundary = false;
|
|
for (auto vIndex : vIter->second->VAdjacent)
|
|
{
|
|
if (IsSupervertex(vIndex))
|
|
{
|
|
// The triangle has a supervertex, so the removal point
|
|
// is on the boundary of the Delaunay triangulation.
|
|
removalPointOnBoundary = true;
|
|
break;
|
|
}
|
|
}
|
|
|
|
auto const& adjacents = vIter->second->TAdjacent;
|
|
std::vector<int32_t> polygon;
|
|
DeleteRemovalPolygon(vRemovalIndex, adjacents, polygon);
|
|
|
|
if (removalPointOnBoundary)
|
|
{
|
|
RetriangulateBoundaryRemovalPolygon(vRemovalIndex, polygon);
|
|
}
|
|
else
|
|
{
|
|
RetriangulateInteriorRemovalPolygon(vRemovalIndex, polygon);
|
|
}
|
|
|
|
mVertexIndexMap.erase(iter);
|
|
return static_cast<size_t>(vRemovalIndex);
|
|
}
|
|
|
|
// Get the current triangulation including the supervertices and
|
|
// triangles containing supervertices.
|
|
void GetTriangulation(std::vector<Vector2<T>>& vertices,
|
|
std::vector<std::array<size_t, 3>>& triangles)
|
|
{
|
|
vertices.resize(mVertices.size());
|
|
std::copy(mVertices.begin(), mVertices.end(), vertices.begin());
|
|
|
|
auto const& tMap = mGraph.GetTriangles();
|
|
triangles.reserve(tMap.size());
|
|
triangles.clear();
|
|
for (auto const& tri : tMap)
|
|
{
|
|
auto const& tKey = tri.first;
|
|
triangles.push_back({
|
|
static_cast<size_t>(tKey.V[0]),
|
|
static_cast<size_t>(tKey.V[1]),
|
|
static_cast<size_t>(tKey.V[2]) });
|
|
}
|
|
}
|
|
|
|
// Get the current graph, which includes all triangles whether
|
|
// Delaunay or those containing a supervertex.
|
|
inline VETManifoldMesh const& GetGraph() const
|
|
{
|
|
return mGraph;
|
|
}
|
|
|
|
// Queries associated with the mesh of Delaunay triangles. The
|
|
// triangles containing a supervertex are not included in these
|
|
// queries.
|
|
inline size_t GetNumVertices() const
|
|
{
|
|
return mVertices.size();
|
|
}
|
|
|
|
inline std::vector<Vector2<T>> const& GetVertices() const
|
|
{
|
|
return mVertices;
|
|
}
|
|
|
|
size_t GetNumTriangles() const
|
|
{
|
|
if (mTrianglesAndAdjacenciesNeedUpdate)
|
|
{
|
|
UpdateTrianglesAndAdjacencies();
|
|
mTrianglesAndAdjacenciesNeedUpdate = false;
|
|
}
|
|
|
|
return mTriangles.size();
|
|
}
|
|
|
|
std::vector<std::array<size_t, 3>> const& GetTriangles() const
|
|
{
|
|
if (mTrianglesAndAdjacenciesNeedUpdate)
|
|
{
|
|
UpdateTrianglesAndAdjacencies();
|
|
mTrianglesAndAdjacenciesNeedUpdate = false;
|
|
}
|
|
|
|
return mTriangles;
|
|
}
|
|
|
|
std::vector<std::array<size_t, 3>> const& GetAdjacencies() const
|
|
{
|
|
if (mTrianglesAndAdjacenciesNeedUpdate)
|
|
{
|
|
UpdateTrianglesAndAdjacencies();
|
|
mTrianglesAndAdjacenciesNeedUpdate = false;
|
|
}
|
|
|
|
return mAdjacencies;
|
|
}
|
|
|
|
// Get the vertex indices for triangle t. The function returns 'true'
|
|
// when t is a valid triangle index, in which case 'triangle' is valid;
|
|
// otherwise, the function returns 'false' and 'triangle' is invalid.
|
|
bool GetTriangle(size_t t, std::array<size_t, 3>& triangle) const
|
|
{
|
|
if (mTrianglesAndAdjacenciesNeedUpdate)
|
|
{
|
|
UpdateTrianglesAndAdjacencies();
|
|
mTrianglesAndAdjacenciesNeedUpdate = false;
|
|
}
|
|
|
|
if (t < mTriangles.size())
|
|
{
|
|
triangle = mTriangles[t];
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
// Get the indices for triangles adjacent to triangle t. The function
|
|
// returns 'true' when t is a valid triangle index, in which case
|
|
// 'adjacent' is valid; otherwise, the function returns 'false' and
|
|
// 'adjacent' is invalid. When valid, triangle t has ordered vertices
|
|
// <V[0], V[1], V[2]>. The value adjacent[0] is the index for the
|
|
// triangle adjacent to edge <V[0], V[1]>, adjacent[1] is the index
|
|
// for the triangle adjacent to edge <V[1], V[2]>, and adjacent[2] is
|
|
// the index for the triangle adjacent to edge <V[2], V[0]>.
|
|
bool GetAdjacent(size_t t, std::array<size_t, 3>& adjacent) const
|
|
{
|
|
if (mTrianglesAndAdjacenciesNeedUpdate)
|
|
{
|
|
UpdateTrianglesAndAdjacencies();
|
|
mTrianglesAndAdjacenciesNeedUpdate = false;
|
|
}
|
|
|
|
if (t < mAdjacencies.size())
|
|
{
|
|
adjacent = mAdjacencies[t];
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
// Get the convex polygon that is the hull of the Delaunay triangles.
|
|
// The polygon is counterclockwise ordered with vertices V[hull[0]],
|
|
// V[hull[1]], ..., V[hull.size()-1].
|
|
void GetHull(std::vector<size_t>& hull) const
|
|
{
|
|
if (mTrianglesAndAdjacenciesNeedUpdate)
|
|
{
|
|
UpdateTrianglesAndAdjacencies();
|
|
mTrianglesAndAdjacenciesNeedUpdate = false;
|
|
}
|
|
|
|
// The hull edges are shared by the triangles with exactly one
|
|
// supervertex.
|
|
std::map<size_t, size_t> edges;
|
|
auto const& vmap = mGraph.GetVertices();
|
|
for (int32_t v = 0; v < 3; ++v)
|
|
{
|
|
auto vIter = vmap.find(v);
|
|
LogAssert(
|
|
vIter != vmap.end(),
|
|
"Expecting the supervertices to exist in the graph.");
|
|
|
|
for (auto const& adj : vIter->second->TAdjacent)
|
|
{
|
|
for (size_t i0 = 1, i1 = 2, i2 = 0; i2 < 3; i0 = i1, i1 = i2, ++i2)
|
|
{
|
|
if (adj->V[i0] == v)
|
|
{
|
|
if (IsDelaunayVertex(adj->V[i1]) && IsDelaunayVertex(adj->V[i2]))
|
|
{
|
|
edges.insert(std::make_pair(adj->V[i2], adj->V[i1]));
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Repackage the edges into a convex polygon with vertices ordered
|
|
// counterclockwise.
|
|
size_t numEdges = edges.size();
|
|
hull.resize(numEdges);
|
|
auto eIter = edges.begin();
|
|
size_t vStart = eIter->first;
|
|
size_t vNext = eIter->second;
|
|
size_t i = 0;
|
|
hull[0] = vStart;
|
|
while (vNext != vStart)
|
|
{
|
|
hull[++i] = vNext;
|
|
eIter = edges.find(vNext);
|
|
LogAssert(
|
|
eIter != edges.end(),
|
|
"Expecting to find a hull edge.");
|
|
vNext = eIter->second;
|
|
}
|
|
}
|
|
|
|
// Support for searching the Delaunay triangles that contains the
|
|
// point p. If there is a containing triangle, the returned value is
|
|
// a triangle index i with 0 <= i < GetNumTriangles(). If there is
|
|
// not a containing triangle, 'invalid' is returned. The computations
|
|
// are performed using exact rational arithmetic.
|
|
//
|
|
// The SearchInfo input stores information about the triangle search
|
|
// when looking for the triangle (if any) that contains p. The first
|
|
// triangle searched is 'initialTriangle'. On return 'path' stores the
|
|
// ordered triangle indices visited during the search. The last
|
|
// visited triangle has index 'finalTriangle' and vertex indices
|
|
// 'finalV[0,1,2]', stored in counterclockwise order. The last edge of
|
|
// the search is <finalV[0], finalV[1]>. For spatially coherent inputs
|
|
// p for numerous calls to this function, you will want to specify
|
|
// 'finalTriangle' from the previous call as 'initialTriangle' for the
|
|
// next call, which should reduce search times.
|
|
|
|
static size_t constexpr invalid = std::numeric_limits<size_t>::max();
|
|
|
|
struct SearchInfo
|
|
{
|
|
SearchInfo()
|
|
:
|
|
initialTriangle(invalid),
|
|
finalTriangle(invalid),
|
|
finalV{ invalid, invalid, invalid },
|
|
numPath(0),
|
|
path{}
|
|
{
|
|
}
|
|
|
|
size_t initialTriangle;
|
|
size_t finalTriangle;
|
|
std::array<size_t, 3> finalV;
|
|
size_t numPath;
|
|
std::vector<size_t> path;
|
|
};
|
|
|
|
size_t GetContainingTriangle(Vector2<T> const& p, SearchInfo& info) const
|
|
{
|
|
if (mTrianglesAndAdjacenciesNeedUpdate)
|
|
{
|
|
UpdateTrianglesAndAdjacencies();
|
|
mTrianglesAndAdjacenciesNeedUpdate = false;
|
|
}
|
|
|
|
mQueryPoint = p;
|
|
mIRQueryPoint = IRVector{ p[0], p[1] };
|
|
|
|
size_t const numTriangles = mTriangles.size();
|
|
info.path.resize(numTriangles);
|
|
info.numPath = 0;
|
|
size_t tIndex;
|
|
if (info.initialTriangle < numTriangles)
|
|
{
|
|
tIndex = info.initialTriangle;
|
|
}
|
|
else
|
|
{
|
|
info.initialTriangle = 0;
|
|
tIndex = 0;
|
|
}
|
|
|
|
for (size_t t = 0; t < numTriangles; ++t)
|
|
{
|
|
auto const& v = mTriangles[tIndex];
|
|
auto const& adj = mAdjacencies[tIndex];
|
|
|
|
info.finalTriangle = tIndex;
|
|
info.finalV = v;
|
|
info.path[info.numPath++] = tIndex;
|
|
|
|
size_t i0, i1, i2;
|
|
for (i0 = 1, i1 = 2, i2 = 0; i2 < 3; i0 = i1, i1 = i2++)
|
|
{
|
|
// ToLine(pIndex, v0Index, v1Index) uses mQueryPoint when
|
|
// pIndex is set to 'invalid'.
|
|
if (ToLine(invalid, v[i0], v[i1]) > 0)
|
|
{
|
|
tIndex = adj[i0];
|
|
if (tIndex == invalid)
|
|
{
|
|
info.finalV[0] = v[i0];
|
|
info.finalV[1] = v[i1];
|
|
info.finalV[2] = v[i2];
|
|
return invalid;
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
if (i2 == 3)
|
|
{
|
|
return tIndex;
|
|
}
|
|
}
|
|
return invalid;
|
|
}
|
|
|
|
private:
|
|
// The minimum-size rational type of the input points.
|
|
static int32_t constexpr InputNumWords = std::is_same<T, float>::value ? 2 : 4;
|
|
using InputRational = BSNumber<UIntegerFP32<InputNumWords>>;
|
|
using IRVector = Vector2<InputRational>;
|
|
|
|
// The compute type used for exact sign classification.
|
|
static int32_t constexpr ComputeNumWords = std::is_same<T, float>::value ? 36 : 264;
|
|
using ComputeRational = BSNumber<UIntegerFP32<ComputeNumWords>>;
|
|
using CRVector = Vector2<ComputeRational>;
|
|
|
|
// The rectangular domain in which all input points live.
|
|
T mXMin, mYMin, mXMax, mYMax;
|
|
|
|
// The current vertices.
|
|
std::map<Vector2<T>, size_t> mVertexIndexMap;
|
|
std::vector<Vector2<T>> mVertices;
|
|
std::vector<IRVector> mIRVertices;
|
|
|
|
// Sufficient storage for the expression trees related to computing
|
|
// the exact signs in ToLine(...) and ToCircumcircle(...).
|
|
static size_t constexpr maxNumCRPool = 43;
|
|
mutable std::vector<ComputeRational> mCRPool;
|
|
|
|
// Support for inserting a point into the triangulation. The graph
|
|
// is the current triangulation. The mIndex array provides indexing
|
|
using Triangle = VETManifoldMesh::Triangle;
|
|
using DirectedEdgeKeySet = std::set<EdgeKey<true>>;
|
|
using TrianglePtrSet = std::set<Triangle*>;
|
|
VETManifoldMesh mGraph;
|
|
|
|
// Indexing for the vertices of the triangle adjacent to a vertex.
|
|
// The edge adjacent to vertex j is <mIndex[j][0], mIndex[j][1]> and
|
|
// is listed so that the triangle interior is to your left as you walk
|
|
// around the edges.
|
|
std::array<std::array<size_t, 2>, 3> const mIndex;
|
|
|
|
// Wrap the ToLine function for use in retriangulating the removal
|
|
// polygon.
|
|
std::function<int32_t(size_t, size_t, size_t)> mToLineWrapper;
|
|
|
|
|
|
template <typename IntegerType>
|
|
inline bool IsDelaunayVertex(IntegerType vIndex) const
|
|
{
|
|
return vIndex >= 3;
|
|
}
|
|
|
|
template <typename IntegerType>
|
|
inline bool IsSupervertex(IntegerType vIndex) const
|
|
{
|
|
return vIndex < 3;
|
|
}
|
|
|
|
bool GetContainingTriangle(size_t pIndex, Triangle*& tri) const
|
|
{
|
|
size_t const numTriangles = mGraph.GetTriangles().size();
|
|
for (size_t t = 0; t < numTriangles; ++t)
|
|
{
|
|
size_t j;
|
|
for (j = 0; j < 3; ++j)
|
|
{
|
|
size_t v0Index = static_cast<size_t>(tri->V[mIndex[j][0]]);
|
|
size_t v1Index = static_cast<size_t>(tri->V[mIndex[j][1]]);
|
|
if (ToLine(pIndex, v0Index, v1Index) > 0)
|
|
{
|
|
// Point i sees edge <v0,v1> from outside the triangle.
|
|
auto adjTri = tri->T[j];
|
|
if (adjTri)
|
|
{
|
|
// Traverse to the triangle sharing the face.
|
|
tri = adjTri;
|
|
break;
|
|
}
|
|
else
|
|
{
|
|
// We reached a hull edge, so the point is outside
|
|
// the hull.
|
|
return false;
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
if (j == 3)
|
|
{
|
|
// The point is inside all three edges, so the point is
|
|
// inside a triangle.
|
|
return true;
|
|
}
|
|
}
|
|
|
|
LogError(
|
|
"Unexpected termination of loop while searching for a triangle.");
|
|
}
|
|
|
|
void GetAndRemoveInsertionPolygon(size_t pIndex,
|
|
TrianglePtrSet& candidates, DirectedEdgeKeySet& boundary)
|
|
{
|
|
// Locate the triangles that make up the insertion polygon.
|
|
ETManifoldMesh polygon;
|
|
while (candidates.size() > 0)
|
|
{
|
|
Triangle* tri = *candidates.begin();
|
|
candidates.erase(candidates.begin());
|
|
|
|
for (size_t j = 0; j < 3; ++j)
|
|
{
|
|
auto adj = tri->T[j];
|
|
if (adj && candidates.find(adj) == candidates.end())
|
|
{
|
|
size_t v0Index = adj->V[0];
|
|
size_t v1Index = adj->V[1];
|
|
size_t v2Index = adj->V[2];
|
|
if (ToCircumcircle(pIndex, v0Index, v1Index, v2Index) <= 0)
|
|
{
|
|
// Point P is in the circumcircle.
|
|
candidates.insert(adj);
|
|
}
|
|
}
|
|
}
|
|
|
|
auto inserted = polygon.Insert(tri->V[0], tri->V[1], tri->V[2]);
|
|
LogAssert(
|
|
inserted != nullptr,
|
|
"Unexpected insertion failure.");
|
|
|
|
auto removed = mGraph.Remove(tri->V[0], tri->V[1], tri->V[2]);
|
|
LogAssert(
|
|
removed,
|
|
"Unexpected removal failure.");
|
|
}
|
|
|
|
// Get the boundary edges of the insertion polygon.
|
|
for (auto const& element : polygon.GetTriangles())
|
|
{
|
|
Triangle* tri = element.second.get();
|
|
for (size_t j = 0; j < 3; ++j)
|
|
{
|
|
if (!tri->T[j])
|
|
{
|
|
EdgeKey<true> ekey(tri->V[mIndex[j][0]], tri->V[mIndex[j][1]]);
|
|
boundary.insert(ekey);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void Update(size_t pIndex)
|
|
{
|
|
auto const& tmap = mGraph.GetTriangles();
|
|
Triangle* tri = tmap.begin()->second.get();
|
|
if (GetContainingTriangle(pIndex, tri))
|
|
{
|
|
// The point is inside the convex hull. The insertion polygon
|
|
// contains only triangles in the current triangulation; the
|
|
// hull does not change.
|
|
|
|
// Use a depth-first search for those triangles whose
|
|
// circumcircles contain point P.
|
|
TrianglePtrSet candidates;
|
|
candidates.insert(tri);
|
|
|
|
// Get the boundary of the insertion polygon C that contains
|
|
// the triangles whose circumcircles contain point P. Polygon
|
|
// Polygon C contains this point.
|
|
DirectedEdgeKeySet boundary;
|
|
GetAndRemoveInsertionPolygon(pIndex, candidates, boundary);
|
|
|
|
// The insertion polygon consists of the triangles formed by
|
|
// point P and the faces of C.
|
|
for (auto const& key : boundary)
|
|
{
|
|
size_t v0Index = static_cast<size_t>(key.V[0]);
|
|
size_t v1Index = static_cast<size_t>(key.V[1]);
|
|
if (ToLine(pIndex, v0Index, v1Index) < 0)
|
|
{
|
|
auto inserted = mGraph.Insert(static_cast<int32_t>(pIndex),
|
|
key.V[0], key.V[1]);
|
|
LogAssert(
|
|
inserted != nullptr,
|
|
"Unexpected insertion failure.");
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// The point is outside the convex hull. The insertion
|
|
// polygon is formed by point P and any triangles in the
|
|
// current triangulation whose circumcircles contain point P.
|
|
|
|
// Locate the convex hull of the triangles.
|
|
DirectedEdgeKeySet hull;
|
|
for (auto const& element : tmap)
|
|
{
|
|
Triangle* t = element.second.get();
|
|
for (size_t j = 0; j < 3; ++j)
|
|
{
|
|
if (!t->T[j])
|
|
{
|
|
hull.insert(EdgeKey<true>(t->V[mIndex[j][0]], t->V[mIndex[j][1]]));
|
|
}
|
|
}
|
|
}
|
|
|
|
// Iterate over all the hull edges and use the ones visible to
|
|
// point P to locate the insertion polygon.
|
|
auto const& emap = mGraph.GetEdges();
|
|
TrianglePtrSet candidates;
|
|
DirectedEdgeKeySet visible;
|
|
for (auto const& key : hull)
|
|
{
|
|
size_t v0Index = static_cast<size_t>(key.V[0]);
|
|
size_t v1Index = static_cast<size_t>(key.V[1]);
|
|
if (ToLine(pIndex, v0Index, v1Index) > 0)
|
|
{
|
|
auto iter = emap.find(EdgeKey<false>(key.V[0], key.V[1]));
|
|
if (iter != emap.end() && iter->second->T[1] == nullptr)
|
|
{
|
|
auto adj = iter->second->T[0];
|
|
if (adj && candidates.find(adj) == candidates.end())
|
|
{
|
|
size_t a0Index = static_cast<size_t>(adj->V[0]);
|
|
size_t a1Index = static_cast<size_t>(adj->V[1]);
|
|
size_t a2Index = static_cast<size_t>(adj->V[2]);
|
|
if (ToCircumcircle(pIndex, a0Index, a1Index, a2Index) <= 0)
|
|
{
|
|
// Point P is in the circumcircle.
|
|
candidates.insert(adj);
|
|
}
|
|
else
|
|
{
|
|
// Point P is not in the circumcircle but
|
|
// the hull edge is visible.
|
|
visible.insert(key);
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
LogError(
|
|
"This condition should not occur for rational arithmetic.");
|
|
}
|
|
}
|
|
}
|
|
|
|
// Get the boundary of the insertion subpolygon C that
|
|
// contains the triangles whose circumcircles contain point P.
|
|
DirectedEdgeKeySet boundary;
|
|
GetAndRemoveInsertionPolygon(pIndex, candidates, boundary);
|
|
|
|
// The insertion polygon P consists of the triangles formed by
|
|
// point i and the back edges of C and by the visible edges of
|
|
// mGraph-C.
|
|
for (auto const& key : boundary)
|
|
{
|
|
size_t v0Index = static_cast<size_t>(key.V[0]);
|
|
size_t v1Index = static_cast<size_t>(key.V[1]);
|
|
if (ToLine(pIndex, v0Index, v1Index) < 0)
|
|
{
|
|
// This is a back edge of the boundary.
|
|
auto inserted = mGraph.Insert(static_cast<int32_t>(pIndex),
|
|
key.V[0], key.V[1]);
|
|
LogAssert(
|
|
inserted != nullptr,
|
|
"Unexpected insertion failure.");
|
|
}
|
|
}
|
|
for (auto const& key : visible)
|
|
{
|
|
auto inserted = mGraph.Insert(static_cast<int32_t>(pIndex),
|
|
key.V[1], key.V[0]);
|
|
LogAssert(
|
|
inserted != nullptr,
|
|
"Unexpected insertion failure.");
|
|
}
|
|
}
|
|
}
|
|
|
|
static ComputeRational const& Copy(InputRational const& source,
|
|
ComputeRational& target)
|
|
{
|
|
target.SetSign(source.GetSign());
|
|
target.SetBiasedExponent(source.GetBiasedExponent());
|
|
target.GetUInteger().CopyFrom(source.GetUInteger());
|
|
return target;
|
|
}
|
|
|
|
// Given a line with origin V0 and direction <V0,V1> and a query
|
|
// point P, ToLine returns
|
|
// +1, P on right of line
|
|
// -1, P on left of line
|
|
// 0, P on the line
|
|
int32_t ToLine(size_t pIndex, size_t v0Index, size_t v1Index) const
|
|
{
|
|
// The expression tree has 13 nodes consisting of 6 input
|
|
// leaves and 7 compute nodes.
|
|
|
|
// Use interval arithmetic to determine the sign if possible.
|
|
Vector2<T> const& inP = (pIndex != invalid ? mVertices[pIndex] : mQueryPoint);
|
|
Vector2<T> const& inV0 = mVertices[v0Index];
|
|
Vector2<T> const& inV1 = mVertices[v1Index];
|
|
|
|
auto x0 = SWInterval<T>::Sub(inP[0], inV0[0]);
|
|
auto y0 = SWInterval<T>::Sub(inP[1], inV0[1]);
|
|
auto x1 = SWInterval<T>::Sub(inV1[0], inV0[0]);
|
|
auto y1 = SWInterval<T>::Sub(inV1[1], inV0[1]);
|
|
auto x0y1 = x0 * y1;
|
|
auto x1y0 = x1 * y0;
|
|
auto det = x0y1 - x1y0;
|
|
|
|
T constexpr zero = 0;
|
|
if (det[0] > zero)
|
|
{
|
|
return +1;
|
|
}
|
|
else if (det[1] < zero)
|
|
{
|
|
return -1;
|
|
}
|
|
|
|
// The exact sign of the determinant is not known, so compute
|
|
// the determinant using rational arithmetic.
|
|
|
|
// Name the nodes of the expression tree.
|
|
auto const& irP = (pIndex != invalid ? mIRVertices[pIndex] : mIRQueryPoint);
|
|
Vector2<InputRational> const& irV0 = mIRVertices[v0Index];
|
|
Vector2<InputRational> const& irV1 = mIRVertices[v1Index];
|
|
|
|
auto const& crP0 = Copy(irP[0], mCRPool[0]);
|
|
auto const& crP1 = Copy(irP[1], mCRPool[1]);
|
|
auto const& crV00 = Copy(irV0[0], mCRPool[2]);
|
|
auto const& crV01 = Copy(irV0[1], mCRPool[3]);
|
|
auto const& crV10 = Copy(irV1[0], mCRPool[4]);
|
|
auto const& crV11 = Copy(irV1[1], mCRPool[5]);
|
|
auto& crX0 = mCRPool[6];
|
|
auto& crY0 = mCRPool[7];
|
|
auto& crX1 = mCRPool[8];
|
|
auto& crY1 = mCRPool[9];
|
|
auto& crX0Y1 = mCRPool[10];
|
|
auto& crX1Y0 = mCRPool[11];
|
|
auto& crDet = mCRPool[12];
|
|
|
|
// Evaluate the expression tree.
|
|
crX0 = crP0 - crV00;
|
|
crY0 = crP1 - crV01;
|
|
crX1 = crV10 - crV00;
|
|
crY1 = crV11 - crV01;
|
|
crX0Y1 = crX0 * crY1;
|
|
crX1Y0 = crX1 * crY0;
|
|
crDet = crX0Y1 - crX1Y0;
|
|
return crDet.GetSign();
|
|
}
|
|
|
|
// For a triangle with counterclockwise vertices V0, V1 and V2, operator()
|
|
// returns
|
|
// +1, P outside triangle
|
|
// -1, P inside triangle
|
|
// 0, P on triangle
|
|
int32_t ToTriangle(size_t pIndex, size_t v0Index, size_t v1Index, size_t v2Index) const
|
|
{
|
|
int32_t sign0 = ToLine(pIndex, v1Index, v2Index);
|
|
if (sign0 > 0)
|
|
{
|
|
return +1;
|
|
}
|
|
|
|
int32_t sign1 = ToLine(pIndex, v0Index, v2Index);
|
|
if (sign1 < 0)
|
|
{
|
|
return +1;
|
|
}
|
|
|
|
int32_t sign2 = ToLine(pIndex, v0Index, v1Index);
|
|
if (sign2 > 0)
|
|
{
|
|
return +1;
|
|
}
|
|
|
|
return ((sign0 && sign1 && sign2) ? -1 : 0);
|
|
}
|
|
|
|
// For a triangle with counterclockwise vertices V0, V1 and V2 and a
|
|
// query point P, ToCircumcircle returns
|
|
// +1, P outside circumcircle of triangle
|
|
// -1, P inside circumcircle of triangle
|
|
// 0, P on circumcircle of triangle
|
|
int32_t ToCircumcircle(size_t pIndex, size_t v0Index, size_t v1Index, size_t v2Index) const
|
|
{
|
|
// The expression tree has 43 nodes consisting of 8 input
|
|
// leaves and 35 compute nodes.
|
|
|
|
// Use interval arithmetic to determine the sign if possible.
|
|
auto const& inP = mVertices[pIndex];
|
|
Vector2<T> const& inV0 = mVertices[v0Index];
|
|
Vector2<T> const& inV1 = mVertices[v1Index];
|
|
Vector2<T> const& inV2 = mVertices[v2Index];
|
|
|
|
auto x0 = SWInterval<T>::Sub(inV0[0], inP[0]);
|
|
auto y0 = SWInterval<T>::Sub(inV0[1], inP[1]);
|
|
auto s00 = SWInterval<T>::Add(inV0[0], inP[0]);
|
|
auto s01 = SWInterval<T>::Add(inV0[1], inP[1]);
|
|
auto x1 = SWInterval<T>::Sub(inV1[0], inP[0]);
|
|
auto y1 = SWInterval<T>::Sub(inV1[1], inP[1]);
|
|
auto s10 = SWInterval<T>::Add(inV1[0], inP[0]);
|
|
auto s11 = SWInterval<T>::Add(inV1[1], inP[1]);
|
|
auto x2 = SWInterval<T>::Sub(inV2[0], inP[0]);
|
|
auto y2 = SWInterval<T>::Sub(inV2[1], inP[1]);
|
|
auto s20 = SWInterval<T>::Add(inV2[0], inP[0]);
|
|
auto s21 = SWInterval<T>::Add(inV2[1], inP[1]);
|
|
auto t00 = s00 * x0;
|
|
auto t01 = s01 * y0;
|
|
auto t10 = s10 * x1;
|
|
auto t11 = s11 * y1;
|
|
auto t20 = s20 * x2;
|
|
auto t21 = s21 * y2;
|
|
auto z0 = t00 + t01;
|
|
auto z1 = t10 + t11;
|
|
auto z2 = t20 + t21;
|
|
auto y0z1 = y0 * z1;
|
|
auto y0z2 = y0 * z2;
|
|
auto y1z0 = y1 * z0;
|
|
auto y1z2 = y1 * z2;
|
|
auto y2z0 = y2 * z0;
|
|
auto y2z1 = y2 * z1;
|
|
auto c0 = y1z2 - y2z1;
|
|
auto c1 = y2z0 - y0z2;
|
|
auto c2 = y0z1 - y1z0;
|
|
auto x0c0 = x0 * c0;
|
|
auto x1c1 = x1 * c1;
|
|
auto x2c2 = x2 * c2;
|
|
auto det = x0c0 + x1c1 + x2c2;
|
|
|
|
T constexpr zero = 0;
|
|
if (det[0] > zero)
|
|
{
|
|
return -1;
|
|
}
|
|
else if (det[1] < zero)
|
|
{
|
|
return +1;
|
|
}
|
|
|
|
// The exact sign of the determinant is not known, so compute
|
|
// the determinant using rational arithmetic.
|
|
|
|
// Name the nodes of the expression tree.
|
|
auto const& irP = mIRVertices[pIndex];
|
|
Vector2<InputRational> const& irV0 = mIRVertices[v0Index];
|
|
Vector2<InputRational> const& irV1 = mIRVertices[v1Index];
|
|
Vector2<InputRational> const& irV2 = mIRVertices[v2Index];
|
|
|
|
auto const& crP0 = Copy(irP[0], mCRPool[0]);
|
|
auto const& crP1 = Copy(irP[1], mCRPool[1]);
|
|
auto const& crV00 = Copy(irV0[0], mCRPool[2]);
|
|
auto const& crV01 = Copy(irV0[1], mCRPool[3]);
|
|
auto const& crV10 = Copy(irV1[0], mCRPool[4]);
|
|
auto const& crV11 = Copy(irV1[1], mCRPool[5]);
|
|
auto const& crV20 = Copy(irV2[0], mCRPool[6]);
|
|
auto const& crV21 = Copy(irV2[1], mCRPool[7]);
|
|
|
|
auto& crX0 = mCRPool[8];
|
|
auto& crY0 = mCRPool[9];
|
|
auto& crS00 = mCRPool[10];
|
|
auto& crS01 = mCRPool[11];
|
|
auto& crT00 = mCRPool[12];
|
|
auto& crT01 = mCRPool[13];
|
|
auto& crZ0 = mCRPool[14];
|
|
|
|
auto& crX1 = mCRPool[15];
|
|
auto& crY1 = mCRPool[16];
|
|
auto& crS10 = mCRPool[17];
|
|
auto& crS11 = mCRPool[18];
|
|
auto& crT10 = mCRPool[19];
|
|
auto& crT11 = mCRPool[20];
|
|
auto& crZ1 = mCRPool[21];
|
|
|
|
auto& crX2 = mCRPool[22];
|
|
auto& crY2 = mCRPool[23];
|
|
auto& crS20 = mCRPool[24];
|
|
auto& crS21 = mCRPool[25];
|
|
auto& crT20 = mCRPool[26];
|
|
auto& crT21 = mCRPool[27];
|
|
auto& crZ2 = mCRPool[28];
|
|
|
|
auto& crY0Z1 = mCRPool[29];
|
|
auto& crY0Z2 = mCRPool[30];
|
|
auto& crY1Z0 = mCRPool[31];
|
|
auto& crY1Z2 = mCRPool[32];
|
|
auto& crY2Z0 = mCRPool[33];
|
|
auto& crY2Z1 = mCRPool[34];
|
|
|
|
auto& crC0 = mCRPool[35];
|
|
auto& crC1 = mCRPool[36];
|
|
auto& crC2 = mCRPool[37];
|
|
auto& crX0C0 = mCRPool[38];
|
|
auto& crX1C1 = mCRPool[39];
|
|
auto& crX2C2 = mCRPool[40];
|
|
auto& crTerm = mCRPool[41];
|
|
auto& crDet = mCRPool[42];
|
|
|
|
// Evaluate the expression tree.
|
|
crX0 = crV00 - crP0;
|
|
crY0 = crV01 - crP1;
|
|
crS00 = crV00 + crP0;
|
|
crS01 = crV01 + crP1;
|
|
crT00 = crS00 * crX0;
|
|
crT01 = crS01 * crY0;
|
|
crZ0 = crT00 + crT01;
|
|
|
|
crX1 = crV10 - crP0;
|
|
crY1 = crV11 - crP1;
|
|
crS10 = crV10 + crP0;
|
|
crS11 = crV11 + crP1;
|
|
crT10 = crS10 * crX1;
|
|
crT11 = crS11 * crY1;
|
|
crZ1 = crT10 + crT11;
|
|
|
|
crX2 = crV20 - crP0;
|
|
crY2 = crV21 - crP1;
|
|
crS20 = crV20 + crP0;
|
|
crS21 = crV21 + crP1;
|
|
crT20 = crS20 * crX2;
|
|
crT21 = crS21 * crY2;
|
|
crZ2 = crT20 + crT21;
|
|
|
|
crY0Z1 = crY0 * crZ1;
|
|
crY0Z2 = crY0 * crZ2;
|
|
crY1Z0 = crY1 * crZ0;
|
|
crY1Z2 = crY1 * crZ2;
|
|
crY2Z0 = crY2 * crZ0;
|
|
crY2Z1 = crY2 * crZ1;
|
|
|
|
crC0 = crY1Z2 - crY2Z1;
|
|
crC1 = crY2Z0 - crY0Z2;
|
|
crC2 = crY0Z1 - crY1Z0;
|
|
crX0C0 = crX0 * crC0;
|
|
crX1C1 = crX1 * crC1;
|
|
crX2C2 = crX2 * crC2;
|
|
crTerm = crX0C0 + crX1C1;
|
|
crDet = crTerm + crX2C2;
|
|
return -crDet.GetSign();
|
|
}
|
|
|
|
private:
|
|
// Support for triangulating the removal polygon.
|
|
|
|
// Let Vc be a vertex in the removal polygon. If Vc is not an ear, its
|
|
// weight is +infinity. If Vc is an ear, let Vp be its predecessor and
|
|
// let Vn be its successor when traversing counterclockwise. Let P be
|
|
// the removal point. The weight is
|
|
// W = -H(Vp, Vc, Vn, P) / D(Vp, Vc, Vn) = WNumer / WDenom
|
|
// where
|
|
// +- -+
|
|
// | Vp.x Vc.x Vn.x | +- -+
|
|
// D = det | Vp.y Vc.y Vn.y | = det | Vc.x-Vp.x Vn.x-Vp.x |
|
|
// | 1 1 1 | | Vc.y-Vp.y Vn.y-Vp.y |
|
|
// +- -+ +- -+
|
|
// and
|
|
// +- -+
|
|
// | Vp.x Vc.x Vn.x P.x |
|
|
// H = det | Vp.y Vc.y Vn.y P.y |
|
|
// | |Vp|^2 |Vc|^2 |Vn|^2 |P|^2 |
|
|
// | 1 1 1 1 |
|
|
// +- -+
|
|
//
|
|
// +- -+
|
|
// | Vc.x-Vp.x Vn.x-Vp.x P.x-Vp.x |
|
|
// = -det | Vc.y-Vp.y Vn.y-Vp.y P.y-Vp.y |
|
|
// | |Vc|^2-|Vp|^2 |Vn|^2-|Vp|^2 |P|^2-|Vp|^2 |
|
|
// +- -+
|
|
//
|
|
// To use BSNumber-based rationals, the weight is a ratio stored as a
|
|
// pair (WNumer, WDenom) with WDenom > 0. The comparison of weights is
|
|
// WN0/WD0 < WN1/WD1, implemented as WN0*WD1 < WN1*WD0. Additionally,
|
|
// a Boolean flag isFinite is used to distinguish between a finite
|
|
// ratio (for convex vertices) and a weight that is infinite (for
|
|
// reflex vertices).
|
|
class RPWeight
|
|
{
|
|
public:
|
|
enum class Type
|
|
{
|
|
finite,
|
|
infinite,
|
|
unmodifiable
|
|
};
|
|
|
|
RPWeight(Type inType = Type::unmodifiable)
|
|
:
|
|
numerator(0),
|
|
denominator(inType == Type::finite ? 1 : 0),
|
|
type(inType)
|
|
{
|
|
}
|
|
|
|
bool operator<(RPWeight const& other) const
|
|
{
|
|
// finite < infinite < unmodifiable
|
|
|
|
if (type == Type::finite)
|
|
{
|
|
if (other.type == Type::finite)
|
|
{
|
|
ComputeRational lhs = numerator * other.denominator;
|
|
ComputeRational rhs = other.numerator * denominator;
|
|
return lhs < rhs;
|
|
}
|
|
else // other.type is infinite or unmodifiable
|
|
{
|
|
return true;
|
|
}
|
|
}
|
|
else if (type == Type::infinite)
|
|
{
|
|
if (other.type == Type::finite)
|
|
{
|
|
return false;
|
|
}
|
|
else // other.type is infinite or unmodifiable
|
|
{
|
|
return other.type == Type::unmodifiable;
|
|
}
|
|
}
|
|
else // type is unmodifiable
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
|
|
bool operator<=(RPWeight const& other) const
|
|
{
|
|
return !(other < *this);
|
|
}
|
|
|
|
// The finite weight is numerator/denominator with a nonnegative
|
|
// numerator and a positive denominator. If the weight is
|
|
// infinite, the numerator and denominator values are invalid.
|
|
ComputeRational numerator, denominator;
|
|
Type type;
|
|
};
|
|
|
|
class RPVertex
|
|
{
|
|
public:
|
|
RPVertex()
|
|
:
|
|
vIndex(invalid),
|
|
isConvex(false),
|
|
iPrev(invalid),
|
|
iNext(invalid),
|
|
record(nullptr)
|
|
{
|
|
}
|
|
|
|
// The index relative to IncrementalDelaunay2 mVertices[].
|
|
size_t vIndex;
|
|
|
|
// A vertex is either convex or reflex.
|
|
bool isConvex;
|
|
|
|
// Vertex indices for the polygon. These are indices relative to
|
|
// RPPolygon mVertices[].
|
|
size_t iPrev, iNext;
|
|
|
|
// Support for the priority queue of ears.
|
|
typename MinHeap<size_t, RPWeight>::Record* record;
|
|
};
|
|
|
|
class RPPolygon
|
|
{
|
|
public:
|
|
RPPolygon(std::vector<int32_t> const& polygon,
|
|
std::function<int32_t(size_t, size_t, size_t)> const& ToLine)
|
|
:
|
|
mNumActive(polygon.size()),
|
|
mVertices(polygon.size())
|
|
{
|
|
// Create a circular list of the polygon vertices for dynamic
|
|
// removal of vertices.
|
|
size_t const numVertices = mVertices.size();
|
|
for (size_t i = 0; i < numVertices; ++i)
|
|
{
|
|
RPVertex& vertex = mVertices[i];
|
|
vertex.vIndex = static_cast<size_t>(polygon[i]);
|
|
vertex.iPrev = (i > 0 ? i - 1 : numVertices - 1);
|
|
vertex.iNext = (i < numVertices - 1 ? i + 1 : 0);
|
|
}
|
|
|
|
// Create a linear list of the polygon convex vertices and a
|
|
// linear list of the polygon reflex vertices, both used for
|
|
// dynamic removal of vertices.
|
|
for (size_t i = 0; i < numVertices; ++i)
|
|
{
|
|
// Determine whether the vertex is convex or reflex.
|
|
size_t vPrev = invalid, vCurr = invalid, vNext = invalid;
|
|
GetTriangle(i, vPrev, vCurr, vNext);
|
|
mVertices[i].isConvex = (ToLine(vPrev, vCurr, vNext) < 0);
|
|
}
|
|
}
|
|
|
|
inline RPVertex const& Vertex(size_t i) const
|
|
{
|
|
return mVertices[i];
|
|
}
|
|
|
|
inline RPVertex& Vertex(size_t i)
|
|
{
|
|
return mVertices[i];
|
|
}
|
|
|
|
void GetTriangle(size_t i, size_t& vPrev, size_t& vCurr, size_t& vNext) const
|
|
{
|
|
RPVertex const& vertex = mVertices[i];
|
|
vCurr = vertex.vIndex;
|
|
vPrev = mVertices[vertex.iPrev].vIndex;
|
|
vNext = mVertices[vertex.iNext].vIndex;
|
|
}
|
|
|
|
void Classify(size_t i,
|
|
std::function<int32_t(size_t, size_t, size_t)> const& ToLine)
|
|
{
|
|
size_t vPrev = invalid, vCurr = invalid, vNext = invalid;
|
|
GetTriangle(i, vPrev, vCurr, vNext);
|
|
mVertices[i].isConvex = (ToLine(vPrev, vCurr, vNext) < 0);
|
|
}
|
|
|
|
inline size_t GetNumActive() const
|
|
{
|
|
return mNumActive;
|
|
}
|
|
|
|
size_t GetActive() const
|
|
{
|
|
for (size_t i = 0; i < mVertices.size(); ++i)
|
|
{
|
|
if (mVertices[i].iPrev != invalid)
|
|
{
|
|
return i;
|
|
}
|
|
}
|
|
|
|
LogError(
|
|
"Expecting to find an active vertex.");
|
|
}
|
|
|
|
inline void Remove(size_t i)
|
|
{
|
|
// Remove the vertex from the polygon.
|
|
RPVertex& vertex = mVertices[i];
|
|
size_t iPrev = vertex.iPrev;
|
|
size_t iNext = vertex.iNext;
|
|
mVertices[iPrev].iNext = iNext;
|
|
mVertices[iNext].iPrev = iPrev;
|
|
|
|
vertex.vIndex = invalid;
|
|
vertex.isConvex = false;
|
|
vertex.iPrev = invalid;
|
|
vertex.iNext = invalid;
|
|
vertex.record = nullptr;
|
|
|
|
--mNumActive;
|
|
}
|
|
|
|
private:
|
|
size_t mNumActive;
|
|
std::vector<RPVertex> mVertices;
|
|
};
|
|
|
|
RPWeight ComputeWeight(size_t iConvexIndex, int32_t vRemovalIndex,
|
|
RPPolygon& rpPolygon)
|
|
{
|
|
// Get the triangle <VP,VC,VN> with convex vertex VC.
|
|
size_t vPrev = invalid, vCurr = invalid, vNext = invalid;
|
|
rpPolygon.GetTriangle(iConvexIndex, vPrev, vCurr, vNext);
|
|
|
|
auto const& irVP = mIRVertices[vPrev];
|
|
auto const& irVC = mIRVertices[vCurr];
|
|
auto const& irVN = mIRVertices[vNext];
|
|
auto const& irPR = mIRVertices[vRemovalIndex];
|
|
|
|
CRVector VP, VC, VN, PR;
|
|
Copy(irVP[0], VP[0]);
|
|
Copy(irVP[1], VP[1]);
|
|
Copy(irVC[0], VC[0]);
|
|
Copy(irVC[1], VC[1]);
|
|
Copy(irVN[0], VN[0]);
|
|
Copy(irVN[1], VN[1]);
|
|
Copy(irPR[0], PR[0]);
|
|
Copy(irPR[1], PR[1]);
|
|
|
|
auto subVCVP = VC - VP;
|
|
auto subVNVP = VN - VP;
|
|
auto subPRVP = PR - VP;
|
|
auto addVCVP = VC + VP;
|
|
auto addVNVP = VN + VP;
|
|
auto addPRVP = PR + VP;
|
|
auto c20 = DotPerp(subVNVP, subPRVP);
|
|
auto c21 = DotPerp(subPRVP, subVCVP);
|
|
auto c22 = DotPerp(subVCVP, subVNVP);
|
|
auto a20 = Dot(subVCVP, addVCVP);
|
|
auto a21 = Dot(subVNVP, addVNVP);
|
|
auto a22 = Dot(subPRVP, addPRVP);
|
|
|
|
RPWeight weight(RPWeight::Type::finite);
|
|
weight.numerator = a20 * c20 + a21 * c21 + a22 * c22;
|
|
weight.numerator.SetSign(-weight.numerator.GetSign());
|
|
weight.denominator = std::move(c22);
|
|
if (weight.denominator.GetSign() < 0)
|
|
{
|
|
weight.numerator.SetSign(-weight.numerator.GetSign());
|
|
weight.denominator.SetSign(1);
|
|
}
|
|
return weight;
|
|
}
|
|
|
|
void DoEarClipping(MinHeap<size_t, RPWeight>& earHeap,
|
|
std::function<RPWeight(size_t)> const& WeightFunction,
|
|
RPPolygon& rpPolygon)
|
|
{
|
|
// Remove the finite-weight vertices from the priority queue,
|
|
// one at a time.
|
|
size_t i;
|
|
RPWeight weight{};
|
|
while (earHeap.GetNumElements() >= 3)
|
|
{
|
|
// Get the ear of minimum weight. The vertex at index i must
|
|
// be convex.
|
|
(void)earHeap.GetMinimum(i, weight);
|
|
if (weight.type != RPWeight::Type::finite)
|
|
{
|
|
break;
|
|
}
|
|
earHeap.Remove(i, weight);
|
|
|
|
// Get the triangle associated with the ear.
|
|
size_t vPrev = invalid, vCurr = invalid, vNext = invalid;
|
|
rpPolygon.GetTriangle(i, vPrev, vCurr, vNext);
|
|
|
|
// Insert the triangle into the graph.
|
|
auto inserted = mGraph.Insert(
|
|
static_cast<int32_t>(vPrev),
|
|
static_cast<int32_t>(vCurr),
|
|
static_cast<int32_t>(vNext));
|
|
LogAssert(
|
|
inserted != nullptr,
|
|
"Unexpected insertion failure.");
|
|
if (earHeap.GetNumElements() < 3)
|
|
{
|
|
earHeap.Reset(0);
|
|
break;
|
|
}
|
|
|
|
// Remove the vertex from the polygon. The previous and next
|
|
// neighbor indices are required to update the adjacent
|
|
// vertices after the removal.
|
|
RPVertex const& vertex = rpPolygon.Vertex(i);
|
|
size_t iPrev = vertex.iPrev;
|
|
size_t iNext = vertex.iNext;
|
|
rpPolygon.Remove(i);
|
|
|
|
// Removal of the ear can cause an adjacent vertex to become
|
|
// an ear or to stop being an ear.
|
|
bool wasConvex, nowConvex;
|
|
|
|
RPVertex& vertexP = rpPolygon.Vertex(iPrev);
|
|
wasConvex = vertexP.isConvex;
|
|
rpPolygon.Classify(iPrev, mToLineWrapper);
|
|
nowConvex = vertexP.isConvex;
|
|
if (wasConvex)
|
|
{
|
|
// The 'vertex' is convex. If 'vertexP' was convex, it
|
|
// cannot become reflex after the ear is clipped.
|
|
LogAssert(
|
|
nowConvex,
|
|
"Unexpected condition.");
|
|
|
|
if (vertexP.record->value.type != RPWeight::Type::unmodifiable)
|
|
{
|
|
weight = WeightFunction(iPrev);
|
|
earHeap.Update(vertexP.record, weight);
|
|
}
|
|
}
|
|
else // 'vertexP' was reflex
|
|
{
|
|
if (nowConvex)
|
|
{
|
|
if (vertexP.record->value.type != RPWeight::Type::unmodifiable)
|
|
{
|
|
weight = WeightFunction(iPrev);
|
|
earHeap.Update(vertexP.record, weight);
|
|
}
|
|
}
|
|
}
|
|
|
|
RPVertex& vertexN = rpPolygon.Vertex(iNext);
|
|
wasConvex = vertexN.isConvex;
|
|
rpPolygon.Classify(iNext, mToLineWrapper);
|
|
nowConvex = vertexN.isConvex;
|
|
if (wasConvex)
|
|
{
|
|
// The 'vertex' is convex. If 'vertexN' was convex, it
|
|
// cannot become reflex after the ear is clipped.
|
|
LogAssert(
|
|
nowConvex,
|
|
"Unexpected condition.");
|
|
|
|
if (vertexN.record->value.type != RPWeight::Type::unmodifiable)
|
|
{
|
|
weight = WeightFunction(iNext);
|
|
earHeap.Update(vertexN.record, weight);
|
|
}
|
|
}
|
|
else // 'vertexN' was reflex
|
|
{
|
|
if (nowConvex)
|
|
{
|
|
if (vertexN.record->value.type != RPWeight::Type::unmodifiable)
|
|
{
|
|
weight = WeightFunction(iNext);
|
|
earHeap.Update(vertexN.record, weight);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void DeleteRemovalPolygon(int32_t vRemovalIndex,
|
|
std::unordered_set<Triangle*> const& adjacents,
|
|
std::vector<int32_t>& polygon)
|
|
{
|
|
// Get the edges of the removal polygon. The polygon is star
|
|
// shaped relative to the removal position.
|
|
std::map<int32_t, int32_t> edges;
|
|
for (auto const& adj : adjacents)
|
|
{
|
|
size_t i;
|
|
for (i = 0; i < 3; ++i)
|
|
{
|
|
if (vRemovalIndex == adj->V[i])
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
LogAssert(
|
|
i < 3,
|
|
"Unexpected condition.");
|
|
|
|
int32_t opposite1 = adj->V[(i + 1) % 3];
|
|
int32_t opposite2 = adj->V[(i + 2) % 3];
|
|
edges.insert(std::make_pair(opposite1, opposite2));
|
|
}
|
|
|
|
// Remove the triangles.
|
|
for (auto const& edge : edges)
|
|
{
|
|
bool removed = mGraph.Remove(vRemovalIndex, edge.first, edge.second);
|
|
LogAssert(
|
|
removed,
|
|
"Unexpected removal failure.");
|
|
}
|
|
|
|
// Create the removal polygon; its vertices are counterclockwise
|
|
// ordered.
|
|
polygon.reserve(edges.size());
|
|
polygon.clear();
|
|
int32_t vStart = edges.begin()->first;
|
|
int32_t vCurr = edges.begin()->second;
|
|
polygon.push_back(vStart);
|
|
while (vCurr != vStart)
|
|
{
|
|
polygon.push_back(vCurr);
|
|
auto eIter = edges.find(vCurr);
|
|
LogAssert(
|
|
eIter != edges.end(),
|
|
"Unexpected condition.");
|
|
|
|
vCurr = eIter->second;
|
|
}
|
|
}
|
|
|
|
void RetriangulateInteriorRemovalPolygon(int32_t vRemovalIndex,
|
|
std::vector<int32_t> const& polygon)
|
|
{
|
|
// Create a representation of 'polygon' that can be processed
|
|
// using a priority queue.
|
|
RPPolygon rpPolygon(polygon, mToLineWrapper);
|
|
|
|
auto WeightFunction = [this, &rpPolygon, vRemovalIndex](size_t i)
|
|
{
|
|
return ComputeWeight(i, vRemovalIndex, rpPolygon);
|
|
};
|
|
|
|
// Create a priority queue of vertices. Convex vertices have a
|
|
// finite and positive weight. Reflex vertices have a weight of
|
|
// +infinity.
|
|
MinHeap<size_t, RPWeight> earHeap(static_cast<int32_t>(polygon.size()));
|
|
RPWeight const posInfinity(RPWeight::Type::infinite);
|
|
RPWeight weight{};
|
|
for (size_t i = 0; i < polygon.size(); ++i)
|
|
{
|
|
RPVertex& vertex = rpPolygon.Vertex(i);
|
|
if (vertex.isConvex)
|
|
{
|
|
weight = WeightFunction(i);
|
|
}
|
|
else
|
|
{
|
|
weight = posInfinity;
|
|
}
|
|
vertex.record = earHeap.Insert(i, weight);
|
|
}
|
|
|
|
// Remove the finite-weight vertices from the priority queue,
|
|
// one at a time.
|
|
DoEarClipping(earHeap, WeightFunction, rpPolygon);
|
|
LogAssert(
|
|
earHeap.GetNumElements() == 0,
|
|
"Expecting the hole to be completely filled.");
|
|
}
|
|
|
|
void RetriangulateBoundaryRemovalPolygon(int32_t vRemovalIndex,
|
|
std::vector<int32_t> const& polygon)
|
|
{
|
|
size_t const numPolygon = polygon.size();
|
|
if (numPolygon >= 3)
|
|
{
|
|
// Create a representation of 'polygon' that can be processed
|
|
// using a priority queue.
|
|
RPPolygon rpPolygon(polygon, mToLineWrapper);
|
|
|
|
auto WeightFunction = [this, &rpPolygon, vRemovalIndex](size_t i)
|
|
{
|
|
return ComputeWeight(i, vRemovalIndex, rpPolygon);
|
|
};
|
|
|
|
auto ZeroWeightFunction = [](size_t)
|
|
{
|
|
return RPWeight(RPWeight::Type::finite);
|
|
};
|
|
|
|
// Create a priority queue of vertices. The removal index
|
|
// (polygon[0] = vRemovalIndex) and its two vertex neighbors
|
|
// have a weight of +infinity. Of the other vertices, convex
|
|
// vertices have a finite and positive weight and reflex
|
|
// vertices have a weight of +infinity.
|
|
MinHeap<size_t, RPWeight> earHeap(static_cast<int32_t>(polygon.size()));
|
|
RPWeight const rigid(RPWeight::Type::unmodifiable);
|
|
RPWeight const posInfinity(RPWeight::Type::infinite);
|
|
|
|
size_t iPrev = numPolygon - 2, iCurr = iPrev + 1, iNext = 0;
|
|
for (; iNext < numPolygon; iPrev = iCurr, iCurr = iNext, ++iNext)
|
|
{
|
|
RPVertex& vertexPrev = rpPolygon.Vertex(iPrev);
|
|
RPVertex& vertexCurr = rpPolygon.Vertex(iCurr);
|
|
RPVertex& vertexNext = rpPolygon.Vertex(iNext);
|
|
if (IsSupervertex(vertexPrev.vIndex) ||
|
|
IsSupervertex(vertexCurr.vIndex) ||
|
|
IsSupervertex(vertexNext.vIndex))
|
|
{
|
|
vertexCurr.record = earHeap.Insert(iCurr, rigid);
|
|
}
|
|
else
|
|
{
|
|
if (vertexCurr.isConvex)
|
|
{
|
|
vertexCurr.record = earHeap.Insert(iCurr,
|
|
WeightFunction(iCurr));
|
|
}
|
|
else
|
|
{
|
|
vertexCurr.record = earHeap.Insert(iCurr, posInfinity);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Remove the finite-weight vertices from the priority queue,
|
|
// one at a time. This process fills in the subpolygon of the
|
|
// removal polygon that is contained by the Delaunay
|
|
// triangulation.
|
|
DoEarClipping(earHeap, WeightFunction, rpPolygon);
|
|
|
|
// Get the subpolygon of the removal polygon that is external
|
|
// to the Delaunay triangulation.
|
|
size_t numExternal = rpPolygon.GetNumActive();
|
|
std::vector<size_t> external(numExternal);
|
|
iCurr = rpPolygon.GetActive();
|
|
for (size_t i = 0; i < numExternal; ++i)
|
|
{
|
|
external[i] = iCurr;
|
|
rpPolygon.Classify(iCurr, mToLineWrapper);
|
|
iCurr = rpPolygon.Vertex(iCurr).iNext;
|
|
}
|
|
|
|
earHeap.Reset(static_cast<int32_t>(numExternal));
|
|
for (size_t i = 0; i < numExternal; ++i)
|
|
{
|
|
size_t index = external[i];
|
|
RPVertex& vertex = rpPolygon.Vertex(index);
|
|
if (IsSupervertex(vertex.vIndex))
|
|
{
|
|
vertex.record = earHeap.Insert(index, rigid);
|
|
}
|
|
else
|
|
{
|
|
if (vertex.isConvex)
|
|
{
|
|
vertex.record = earHeap.Insert(index,
|
|
ZeroWeightFunction(index));
|
|
}
|
|
else
|
|
{
|
|
vertex.record = earHeap.Insert(index, posInfinity);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Remove the finite-weight vertices from the priority queue,
|
|
// one at a time. This process fills in a portion or all of
|
|
// the subpolygon of the removal polygon that is external to
|
|
// the Delaunay triangulation.
|
|
DoEarClipping(earHeap, ZeroWeightFunction, rpPolygon);
|
|
if (earHeap.GetNumElements() == 0)
|
|
{
|
|
// The external polygon contained only 1 supervertex.
|
|
return;
|
|
}
|
|
|
|
// The remaining external polygon is a triangle fan with
|
|
// 2 or 3 supervertices.
|
|
numExternal = rpPolygon.GetNumActive();
|
|
external.resize(numExternal);
|
|
iCurr = rpPolygon.GetActive();
|
|
for (size_t i = 0; i < numExternal; ++i)
|
|
{
|
|
external[i] = iCurr;
|
|
rpPolygon.Classify(iCurr, mToLineWrapper);
|
|
iCurr = rpPolygon.Vertex(iCurr).iNext;
|
|
}
|
|
|
|
earHeap.Reset(static_cast<int32_t>(numExternal));
|
|
iPrev = numExternal - 2;
|
|
iCurr = iPrev + 1;
|
|
iNext = 0;
|
|
for (; iNext < numExternal; iPrev = iCurr, iCurr = iNext, ++iNext)
|
|
{
|
|
size_t index = external[iCurr];
|
|
RPVertex& vertexPrev = rpPolygon.Vertex(external[iPrev]);
|
|
RPVertex& vertexCurr = rpPolygon.Vertex(index);
|
|
RPVertex& vertexNext = rpPolygon.Vertex(external[iNext]);
|
|
if (IsSupervertex(vertexCurr.vIndex))
|
|
{
|
|
if (IsDelaunayVertex(vertexPrev.vIndex) ||
|
|
IsDelaunayVertex(vertexNext.vIndex))
|
|
{
|
|
LogAssert(
|
|
vertexCurr.isConvex,
|
|
"Unexpected condition.");
|
|
|
|
vertexCurr.record = earHeap.Insert(index,
|
|
ZeroWeightFunction(index));
|
|
}
|
|
else
|
|
{
|
|
vertexCurr.record = earHeap.Insert(index, rigid);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
vertexCurr.record = earHeap.Insert(index, posInfinity);
|
|
}
|
|
}
|
|
|
|
// Remove the finite-weight vertices from the priority queue,
|
|
// one at a time. This process fills in the triangle fan of
|
|
// the subpolygon of the removal polygon that is external to
|
|
// the Delaunay triangulation.
|
|
DoEarClipping(earHeap, ZeroWeightFunction, rpPolygon);
|
|
LogAssert(
|
|
earHeap.GetNumElements() == 0,
|
|
"Expecting the hole to be completely filled.");
|
|
}
|
|
else // numPolygon == 2
|
|
{
|
|
int32_t vOtherIndex;
|
|
if (polygon[0] == vRemovalIndex)
|
|
{
|
|
vOtherIndex = polygon[1];
|
|
}
|
|
else
|
|
{
|
|
vOtherIndex = polygon[0];
|
|
}
|
|
|
|
mGraph.Clear();
|
|
for (int32_t i0 = 2, i1 = 0; i1 < 3; i0 = i1++)
|
|
{
|
|
auto inserted = mGraph.Insert(vOtherIndex, i0, i1);
|
|
LogAssert(
|
|
inserted != nullptr,
|
|
"Unexpected insertion failure.");
|
|
}
|
|
}
|
|
}
|
|
|
|
private:
|
|
// Support for queries associated with the mesh of Delaunay triangles.
|
|
mutable std::vector<std::array<size_t, 3>> mTriangles;
|
|
mutable std::vector<std::array<size_t, 3>> mAdjacencies;
|
|
mutable bool mTrianglesAndAdjacenciesNeedUpdate;
|
|
mutable Vector2<T> mQueryPoint;
|
|
mutable IRVector mIRQueryPoint;
|
|
|
|
void UpdateTrianglesAndAdjacencies() const
|
|
{
|
|
// Assign integer values to the triangles.
|
|
auto const& tmap = mGraph.GetTriangles();
|
|
if (tmap.size() == 0)
|
|
{
|
|
return;
|
|
}
|
|
|
|
std::unordered_map<Triangle*, size_t> permute;
|
|
permute[nullptr] = invalid;
|
|
size_t numTriangles = 0;
|
|
for (auto const& element : tmap)
|
|
{
|
|
if (IsDelaunayVertex(element.first.V[0]) &&
|
|
IsDelaunayVertex(element.first.V[1]) &&
|
|
IsDelaunayVertex(element.first.V[2]))
|
|
{
|
|
permute[element.second.get()] = numTriangles++;
|
|
}
|
|
else
|
|
{
|
|
permute[element.second.get()] = invalid;
|
|
}
|
|
}
|
|
|
|
mTriangles.resize(numTriangles);
|
|
mAdjacencies.resize(numTriangles);
|
|
size_t t = 0;
|
|
for (auto const& element : tmap)
|
|
{
|
|
auto const& tri = element.second;
|
|
if (permute[element.second.get()] != invalid)
|
|
{
|
|
for (size_t j = 0; j < 3; ++j)
|
|
{
|
|
mTriangles[t][j] = tri->V[j];
|
|
mAdjacencies[t][j] = permute[tri->T[j]];
|
|
}
|
|
++t;
|
|
}
|
|
}
|
|
}
|
|
};
|
|
}
|
|
|