// 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;
                }
            }
        }
    };
}