GteDelaunay2.h

Go to the documentation of this file.

// David Eberly, Geometric Tools, Redmond WA 98052
// Copyright (c) 1998-2017
// Distributed under the Boost Software License, Version 1.0.
// http://www.boost.org/LICENSE_1_0.txt
// http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
// File Version: 3.0.0 (2016/06/19)

#pragma once

#include <LowLevel/GteLogger.h>
#include <Mathematics/GteETManifoldMesh.h>
#include <Mathematics/GtePrimalQuery2.h>
#include <Mathematics/GteLine.h>
#include <vector>

// Delaunay triangulation of points (intrinsic dimensionality 2).
//   VQ = number of vertices
//   V  = array of vertices
//   TQ = number of triangles
//   I  = Array of 3-tuples of indices into V that represent the triangles
//        (3*TQ total elements).  Access via GetIndices(*).
//   A  = Array of 3-tuples of indices into I that represent the adjacent
//        triangles (3*TQ total elements).  Access via GetAdjacencies(*).
// The i-th triangle has vertices
//   vertex[0] = V[I[3*i+0]]
//   vertex[1] = V[I[3*i+1]]
//   vertex[2] = V[I[3*i+2]]
// and edge index pairs
//   edge[0] = <I[3*i+0],I[3*i+1]>
//   edge[1] = <I[3*i+1],I[3*i+2]>
//   edge[2] = <I[3*i+2],I[3*i+0]>
// The triangles adjacent to these edges have indices
//   adjacent[0] = A[3*i+0] is the triangle sharing edge[0]
//   adjacent[1] = A[3*i+1] is the triangle sharing edge[1]
//   adjacent[2] = A[3*i+2] is the triangle sharing edge[2]
// If there is no adjacent triangle, the A[*] value is set to -1.  The
// triangle adjacent to edge[j] has vertices
//   adjvertex[0] = V[I[3*adjacent[j]+0]]
//   adjvertex[1] = V[I[3*adjacent[j]+1]]
//   adjvertex[2] = V[I[3*adjacent[j]+2]]
// The only way to ensure a correct result for the input vertices (assumed to
// be exact) is to choose ComputeType for exact rational arithmetic.  You may
// use BSNumber.  No divisions are performed in this computation, so you do
// not have to use BSRational.
//
// The worst-case choices of N for Real of type BSNumber or BSRational with
// integer storage UIntegerFP32<N> are listed in the next table.  The numerical
// computations are encapsulated in PrimalQuery2<Real>::ToLine and
// PrimalQuery2<Real>::ToCircumcircle, the latter query the dominant one in
// determining N.  We recommend using only BSNumber, because no divisions are
// performed in the convex-hull computations.
//
//    input type | compute type | N
//    -----------+--------------+------
//    float      | BSNumber     |    35
//    double     | BSNumber     |   263
//    float      | BSRational   | 12302
//    double     | BSRational   | 92827

namespace gte
{

template <typename InputType, typename ComputeType>
class Delaunay2
{
public:
    // The class is a functor to support computing the Delaunay triangulation
    // of multiple data sets using the same class object.
    virtual ~Delaunay2();
    Delaunay2();

    // The input is the array of vertices whose Delaunay triangulation is
    // required.  The epsilon value is used to determine the intrinsic
    // dimensionality of the vertices (d = 0, 1, or 2).  When epsilon is
    // positive, the determination is fuzzy--vertices approximately the same
    // point, approximately on a line, or planar.  The return value is 'true'
    // if and only if the hull construction is successful.
    bool operator()(int numVertices, Vector2<InputType> const* vertices, InputType epsilon);

    // Dimensional information.  If GetDimension() returns 1, the points lie
    // on a line P+t*D (fuzzy comparison when epsilon > 0).  You can sort
    // these if you need a polyline output by projecting onto the line each
    // vertex X = P+t*D, where t = Dot(D,X-P).
    inline InputType GetEpsilon() const;
    inline int GetDimension() const;
    inline Line2<InputType> const& GetLine() const;

    // Member access.
    inline int GetNumVertices() const;
    inline int GetNumUniqueVertices() const;
    inline int GetNumTriangles() const;
    inline Vector2<InputType> const* GetVertices() const;
    inline PrimalQuery2<ComputeType> const& GetQuery() const;
    inline ETManifoldMesh const& GetGraph() const;
    inline std::vector<int> const& GetIndices() const;
    inline std::vector<int> const& GetAdjacencies() const;

    // If 'vertices' has no duplicates, GetDuplicates()[i] = i for all i.
    // If vertices[i] is the first occurrence of a vertex and if vertices[j]
    // is found later, then GetDuplicates()[j] = i.
    inline std::vector<int> const& GetDuplicates() const;

    // Locate those triangle edges that do not share other triangles.  The
    // returned array has hull.size() = 2*numEdges, each pair representing an
    // edge.  The edges are not ordered, but the pair of vertices for an edge
    // is ordered so that they conform to a counterclockwise traversal of the
    // hull.  The return value is 'true' iff the dimension is 2.
    bool GetHull(std::vector<int>& hull) const;

    // Get the vertex indices for triangle i.  The function returns 'true'
    // when the dimension is 2 and i is a valid triangle index, in which case
    // the vertices are valid; otherwise, the function returns 'false' and the
    // vertices are invalid.
    bool GetIndices(int i, std::array<int, 3>& indices) const;

    // Get the indices for triangles adjacent to triangle i.  The function
    // returns 'true' when the dimension is 2 and if i is a valid triangle
    // index, in which case the adjacencies are valid; otherwise, the function
    // returns 'false' and the adjacencies are invalid.
    bool GetAdjacencies(int i, std::array<int, 3>& adjacencies) const;

    // Support for searching the triangulation for a triangle that contains
    // a point.  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, -1 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 those (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.
    struct SearchInfo
    {
        int initialTriangle;
        int numPath;
        std::vector<int> path;
        int finalTriangle;
        std::array<int, 3> finalV;
    };
    int GetContainingTriangle(Vector2<InputType> const& p, SearchInfo& info) const;

protected:
    // Support for incremental Delaunay triangulation.
    typedef ETManifoldMesh::Triangle Triangle;
    bool GetContainingTriangle(int i, std::shared_ptr<Triangle>& tri) const;
    bool GetAndRemoveInsertionPolygon(int i, std::set<std::shared_ptr<Triangle>>& candidates,
        std::set<EdgeKey<true>>& boundary);
    bool Update(int i);

    // The epsilon value is used for fuzzy determination of intrinsic
    // dimensionality.  If the dimension is 0 or 1, the constructor returns
    // early.  The caller is responsible for retrieving the dimension and
    // taking an alternate path should the dimension be smaller than 2.
    // If the dimension is 0, the caller may as well treat all vertices[]
    // as a single point, say, vertices[0].  If the dimension is 1, the
    // caller can query for the approximating line and project vertices[]
    // onto it for further processing.
    InputType mEpsilon;
    int mDimension;
    Line2<InputType> mLine;

    // The array of vertices used for geometric queries.  If you want to be
    // certain of a correct result, choose ComputeType to be BSNumber.
    std::vector<Vector2<ComputeType>> mComputeVertices;
    PrimalQuery2<ComputeType> mQuery;

    // The graph information.
    int mNumVertices;
    int mNumUniqueVertices;
    int mNumTriangles;
    Vector2<InputType> const* mVertices;
    ETManifoldMesh mGraph;
    std::vector<int> mIndices;
    std::vector<int> mAdjacencies;

    // If a vertex occurs multiple times in the 'vertices' input to the
    // constructor, the first processed occurrence of that vertex has an
    // index stored in this array.  If there are no duplicates, then
    // mDuplicates[i] = i for all i.

    struct ProcessedVertex
    {
        ProcessedVertex();
        ProcessedVertex(Vector2<InputType> const& inVertex, int inLocation);
        bool operator<(ProcessedVertex const& v) const;

        Vector2<InputType> vertex;
        int location;
    };

    std::vector<int> mDuplicates;

    // 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.  TODO: Use the "opposite edge" to be consistent with that of
    // TetrahedronKey.  The "opposite" idea extends easily to higher
    // dimensions.
    std::array<std::array<int, 2>, 3> mIndex;
};


template <typename InputType, typename ComputeType>
Delaunay2<InputType, ComputeType>::~Delaunay2()
{
}

template <typename InputType, typename ComputeType>
Delaunay2<InputType, ComputeType>::Delaunay2()
    :
    mEpsilon((InputType)0),
    mDimension(0),
    mLine(Vector2<InputType>::Zero(), Vector2<InputType>::Zero()),
    mNumVertices(0),
    mNumUniqueVertices(0),
    mNumTriangles(0),
    mVertices(nullptr)
{
    // INVESTIGATE.  If the initialization of mIndex is placed in the
    // constructor initializer list, MSVS 2012 generates an internal
    // compiler error.
    mIndex = { { { 0, 1 }, { 1, 2 }, { 2, 0 } } };
}

template <typename InputType, typename ComputeType>
bool Delaunay2<InputType, ComputeType>::operator()(int numVertices,
    Vector2<InputType> const* vertices, InputType epsilon)
{
    mEpsilon = std::max(epsilon, (InputType)0);
    mDimension = 0;
    mLine.origin = Vector2<InputType>::Zero();
    mLine.direction = Vector2<InputType>::Zero();
    mNumVertices = numVertices;
    mNumUniqueVertices = 0;
    mNumTriangles = 0;
    mVertices = vertices;
    mGraph.Clear();
    mIndices.clear();
    mAdjacencies.clear();
    mDuplicates.resize(std::max(numVertices, 3));

    int i, j;
    if (mNumVertices < 3)
    {
        // Delaunay2 should be called with at least three points.
        return false;
    }

    IntrinsicsVector2<InputType> info(mNumVertices, vertices, mEpsilon);
    if (info.dimension == 0)
    {
        // mDimension is 0; mGraph, mIndices, and mAdjacencies are empty
        return false;
    }

    if (info.dimension == 1)
    {
        // The set is (nearly) collinear.
        mDimension = 1;
        mLine = Line2<InputType>(info.origin, info.direction[0]);
        return false;
    }

    mDimension = 2;

    // Compute the vertices for the queries.
    mComputeVertices.resize(mNumVertices);
    mQuery.Set(mNumVertices, &mComputeVertices[0]);
    for (i = 0; i < mNumVertices; ++i)
    {
        for (j = 0; j < 2; ++j)
        {
            mComputeVertices[i][j] = vertices[i][j];
        }
    }

    // Insert the (nondegenerate) triangle constructed by the call to
    // GetInformation.  This is necessary for the circumcircle-visibility
    // algorithm to work correctly.
    if (!info.extremeCCW)
    {
        std::swap(info.extreme[1], info.extreme[2]);
    }
    if (!mGraph.Insert(info.extreme[0], info.extreme[1], info.extreme[2]))
    {
        return false;
    }

    // Incrementally update the triangulation.  The set of processed points
    // is maintained to eliminate duplicates, either in the original input
    // points or in the points obtained by snap rounding.
    std::set<ProcessedVertex> processed;
    for (i = 0; i < 3; ++i)
    {
        j = info.extreme[i];
        processed.insert(ProcessedVertex(vertices[j], j));
        mDuplicates[j] = j;
    }
    for (i = 0; i < mNumVertices; ++i)
    {
        ProcessedVertex v(vertices[i], i);
        auto iter = processed.find(v);
        if (iter == processed.end())
        {
            if (!Update(i))
            {
                // A failure can occur if ComputeType is not an exact
                // arithmetic type.
                return false;
            }
            processed.insert(v);
            mDuplicates[i] = i;
        }
        else
        {
            mDuplicates[i] = iter->location;
        }
    }
    mNumUniqueVertices = static_cast<int>(processed.size());

    // Assign integer values to the triangles for use by the caller.
    std::map<std::shared_ptr<Triangle>, int> permute;
    i = -1;
    permute[nullptr] = i++;
    for (auto const& element : mGraph.GetTriangles())
    {
        permute[element.second] = i++;
    }

    // Put Delaunay triangles into an array (vertices and adjacency info).
    mNumTriangles = static_cast<int>(mGraph.GetTriangles().size());
    int numindices = 3 * mNumTriangles;
    if (numindices > 0)
    {
        mIndices.resize(numindices);
        mAdjacencies.resize(numindices);
        i = 0;
        for (auto const& element : mGraph.GetTriangles())
        {
            std::shared_ptr<Triangle> tri = element.second;
            for (j = 0; j < 3; ++j, ++i)
            {
                mIndices[i] = tri->V[j];
                mAdjacencies[i] = permute[tri->T[j].lock()];
            }
        }
    }

    return true;
}

template <typename InputType, typename ComputeType> inline
InputType Delaunay2<InputType, ComputeType>::GetEpsilon() const
{
    return mEpsilon;
}

template <typename InputType, typename ComputeType> inline
int Delaunay2<InputType, ComputeType>::GetDimension() const
{
    return mDimension;
}

template <typename InputType, typename ComputeType> inline
Line2<InputType> const& Delaunay2<InputType, ComputeType>::GetLine() const
{
    return mLine;
}

template <typename InputType, typename ComputeType> inline
int Delaunay2<InputType, ComputeType>::GetNumVertices() const
{
    return mNumVertices;
}

template <typename InputType, typename ComputeType> inline
int Delaunay2<InputType, ComputeType>::GetNumUniqueVertices() const
{
    return mNumUniqueVertices;
}

template <typename InputType, typename ComputeType> inline
int Delaunay2<InputType, ComputeType>::GetNumTriangles() const
{
    return mNumTriangles;
}

template <typename InputType, typename ComputeType> inline
Vector2<InputType> const* Delaunay2<InputType, ComputeType>::GetVertices() const
{
    return mVertices;
}

template <typename InputType, typename ComputeType> inline
PrimalQuery2<ComputeType> const& Delaunay2<InputType, ComputeType>::GetQuery() const
{
    return mQuery;
}

template <typename InputType, typename ComputeType> inline
ETManifoldMesh const& Delaunay2<InputType, ComputeType>::GetGraph() const
{
    return mGraph;
}

template <typename InputType, typename ComputeType> inline
std::vector<int> const& Delaunay2<InputType, ComputeType>::GetIndices() const
{
    return mIndices;
}

template <typename InputType, typename ComputeType> inline
std::vector<int> const& Delaunay2<InputType, ComputeType>::GetAdjacencies() const
{
    return mAdjacencies;
}

template <typename InputType, typename ComputeType> inline
std::vector<int> const& Delaunay2<InputType, ComputeType>::GetDuplicates() const
{
    return mDuplicates;
}

template <typename InputType, typename ComputeType>
bool Delaunay2<InputType, ComputeType>::GetHull(std::vector<int>& hull) const
{
    if (mDimension == 2)
    {
        // Count the number of edges that are not shared by two triangles.
        int numEdges = 0;
        for (auto adj : mAdjacencies)
        {
            if (adj == -1)
            {
                ++numEdges;
            }
        }

        if (numEdges > 0)
        {
            // Enumerate the edges.
            hull.resize(2 * numEdges);
            int current = 0, i = 0;
            for (auto adj : mAdjacencies)
            {
                if (adj == -1)
                {
                    int tri = i / 3, j = i % 3;
                    hull[current++] = mIndices[3 * tri + j];
                    hull[current++] = mIndices[3 * tri + ((j + 1) % 3)];
                }
                ++i;
            }
            return true;
        }
        else
        {
            LogError("Unexpected.  There must be at least one triangle.");
        }
    }
    else
    {
        LogError("The dimension must be 2.");
    }
    return false;
}

template <typename InputType, typename ComputeType>
bool Delaunay2<InputType, ComputeType>::GetIndices(int i, std::array<int, 3>& indices) const
{
    if (mDimension == 2)
    {
        int numTriangles = static_cast<int>(mIndices.size() / 3);
        if (0 <= i && i < numTriangles)
        {
            indices[0] = mIndices[3 * i];
            indices[1] = mIndices[3 * i + 1];
            indices[2] = mIndices[3 * i + 2];
            return true;
        }
    }
    else
    {
        LogError("The dimension must be 2.");
    }
    return false;
}

template <typename InputType, typename ComputeType>
bool Delaunay2<InputType, ComputeType>::GetAdjacencies(int i, std::array<int, 3>& adjacencies) const
{
    if (mDimension == 2)
    {
        int numTriangles = static_cast<int>(mIndices.size() / 3);
        if (0 <= i && i < numTriangles)
        {
            adjacencies[0] = mAdjacencies[3 * i];
            adjacencies[1] = mAdjacencies[3 * i + 1];
            adjacencies[2] = mAdjacencies[3 * i + 2];
            return true;
        }
    }
    else
    {
        LogError("The dimension must be 2.");
    }
    return false;
}

template <typename InputType, typename ComputeType>
int Delaunay2<InputType, ComputeType>::GetContainingTriangle(
    Vector2<InputType> const& p, SearchInfo& info) const
{
    if (mDimension == 2)
    {
        Vector2<ComputeType> test{ p[0], p[1] };

        int numTriangles = static_cast<int>(mIndices.size() / 3);
        info.path.resize(numTriangles);
        info.numPath = 0;
        int triangle;
        if (0 <= info.initialTriangle && info.initialTriangle < numTriangles)
        {
            triangle = info.initialTriangle;
        }
        else
        {
            info.initialTriangle = 0;
            triangle = 0;
        }

        // Use triangle edges as binary separating lines.
        for (int i = 0; i < numTriangles; ++i)
        {
            int ibase = 3 * triangle;
            int const* v = &mIndices[ibase];

            info.path[info.numPath++] = triangle;
            info.finalTriangle = triangle;
            info.finalV[0] = v[0];
            info.finalV[1] = v[1];
            info.finalV[2] = v[2];

            if (mQuery.ToLine(test, v[0], v[1]) > 0)
            {
                triangle = mAdjacencies[ibase];
                if (triangle == -1)
                {
                    info.finalV[0] = v[0];
                    info.finalV[1] = v[1];
                    info.finalV[2] = v[2];
                    return -1;
                }
                continue;
            }

            if (mQuery.ToLine(test, v[1], v[2]) > 0)
            {
                triangle = mAdjacencies[ibase + 1];
                if (triangle == -1)
                {
                    info.finalV[0] = v[1];
                    info.finalV[1] = v[2];
                    info.finalV[2] = v[0];
                    return -1;
                }
                continue;
            }

            if (mQuery.ToLine(test, v[2], v[0]) > 0)
            {
                triangle = mAdjacencies[ibase + 2];
                if (triangle == -1)
                {
                    info.finalV[0] = v[2];
                    info.finalV[1] = v[0];
                    info.finalV[2] = v[1];
                    return -1;
                }
                continue;
            }

            return triangle;
        }
    }
    else
    {
        LogError("The dimension must be 2.");
    }
    return -1;
}

template <typename InputType, typename ComputeType>
bool Delaunay2<InputType, ComputeType>::GetContainingTriangle(int i, std::shared_ptr<Triangle>& tri) const
{
    int numTriangles = static_cast<int>(mGraph.GetTriangles().size());
    for (int t = 0; t < numTriangles; ++t)
    {
        int j;
        for (j = 0; j < 3; ++j)
        {
            int v0 = tri->V[mIndex[j][0]];
            int v1 = tri->V[mIndex[j][1]];
            if (mQuery.ToLine(i, v0, v1) > 0)
            {
                // Point i sees edge <v0,v1> from outside the triangle.
                auto adjTri = tri->T[j].lock();
                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.  TODO:  Once a hull data structure is in place,
                    // return tri->T[j] as the candidate for starting a search
                    // for visible hull edges.
                    return false;
                }
            }

        }

        if (j == 3)
        {
            // The point is inside all four edges, so the point is inside
            // a triangle.
            return true;
        }
    }

    LogError("Unexpected termination of GetContainingTriangle.");
    return false;
}

template <typename InputType, typename ComputeType>
bool Delaunay2<InputType, ComputeType>::GetAndRemoveInsertionPolygon(int i,
    std::set<std::shared_ptr<Triangle>>& candidates, std::set<EdgeKey<true>>& boundary)
{
    // Locate the triangles that make up the insertion polygon.
    ETManifoldMesh polygon;
    while (candidates.size() > 0)
    {
        std::shared_ptr<Triangle> tri = *candidates.begin();
        candidates.erase(candidates.begin());

        for (int j = 0; j < 3; ++j)
        {
            auto adj = tri->T[j].lock();
            if (adj && candidates.find(adj) == candidates.end())
            {
                int a0 = adj->V[0];
                int a1 = adj->V[1];
                int a2 = adj->V[2];
                if (mQuery.ToCircumcircle(i, a0, a1, a2) <= 0)
                {
                    // Point i is in the circumcircle.
                    candidates.insert(adj);
                }
            }
        }

        if (!polygon.Insert(tri->V[0], tri->V[1], tri->V[2]))
        {
            return false;
        }
        if (!mGraph.Remove(tri->V[0], tri->V[1], tri->V[2]))
        {
            return false;
        }
    }

    // Get the boundary edges of the insertion polygon.
    for (auto const& element : polygon.GetTriangles())
    {
        std::shared_ptr<Triangle> tri = element.second;
        for (int j = 0; j < 3; ++j)
        {
            if (!tri->T[j].lock())
            {
                boundary.insert(EdgeKey<true>(tri->V[mIndex[j][0]], tri->V[mIndex[j][1]]));
            }
        }
    }
    return true;
}

template <typename InputType, typename ComputeType>
bool Delaunay2<InputType, ComputeType>::Update(int i)
{
    // The return value of mGraph.Insert(...) is nullptr if there was a
    // failure to insert.  The Update function will return 'false' when
    // the insertion fails.

    auto const& tmap = mGraph.GetTriangles();
    std::shared_ptr<Triangle> tri = tmap.begin()->second;
    if (GetContainingTriangle(i, 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 i.
        std::set<std::shared_ptr<Triangle>> candidates;
        candidates.insert(tri);

        // Get the boundary of the insertion polygon C that contains the
        // triangles whose circumcircles contain point i.  C contains the
        // point i.
        std::set<EdgeKey<true>> boundary;
        if (!GetAndRemoveInsertionPolygon(i, candidates, boundary))
        {
            return false;
        }

        // The insertion polygon consists of the triangles formed by
        // point i and the faces of C.
        for (auto const& key : boundary)
        {
            int v0 = key.V[0];
            int v1 = key.V[1];
            if (mQuery.ToLine(i, v0, v1) < 0)
            {
                if (!mGraph.Insert(i, v0, v1))
                {
                    return false;
                }
            }
            // else:  Point i is on an edge of 'tri', so the
            // subdivision has degenerate triangles.  Ignore these.
        }
    }
    else
    {
        // The point is outside the convex hull.  The insertion polygon
        // is formed by point i and any triangles in the current
        // triangulation whose circumcircles contain point i.

        // Locate the convex hull of the triangles.  TODO:  Maintain a hull
        // data structure that is updated incrementally.
        std::set<EdgeKey<true>> hull;
        for (auto const& element : tmap)
        {
            std::shared_ptr<Triangle> t = element.second;
            for (int j = 0; j < 3; ++j)
            {
                if (!t->T[j].lock())
                {
                    hull.insert(EdgeKey<true>(t->V[mIndex[j][0]], t->V[mIndex[j][1]]));
                }
            }
        }

        // TODO:  Until the hull change, for now just iterate over all the
        // hull edges and use the ones visible to point i to locate the
        // insertion polygon.
        auto const& emap = mGraph.GetEdges();
        std::set<std::shared_ptr<Triangle>> candidates;
        std::set<EdgeKey<true>> visible;
        for (auto const& key : hull)
        {
            int v0 = key.V[0];
            int v1 = key.V[1];
            if (mQuery.ToLine(i, v0, v1) > 0)
            {
                auto iter = emap.find(EdgeKey<false>(v0, v1));
                if (iter != emap.end() && iter->second->T[1].lock() == nullptr)
                {
                    auto adj = iter->second->T[0].lock();
                    if (adj && candidates.find(adj) == candidates.end())
                    {
                        int a0 = adj->V[0];
                        int a1 = adj->V[1];
                        int a2 = adj->V[2];
                        if (mQuery.ToCircumcircle(i, a0, a1, a2) <= 0)
                        {
                            // Point i is in the circumcircle.
                            candidates.insert(adj);
                        }
                        else
                        {
                            // Point i is not in the circumcircle but the hull
                            // edge is visible.
                            visible.insert(key);
                        }
                    }
                }
                else
                {
                    LogError("Unexpected condition (ComputeType not exact?)");
                    return false;
                }
            }
        }

        // Get the boundary of the insertion subpolygon C that contains the
        // triangles whose circumcircles contain point i.
        std::set<EdgeKey<true>> boundary;
        if (!GetAndRemoveInsertionPolygon(i, candidates, boundary))
        {
            return false;
        }

        // The insertion polygon P consists of the triangles formed by point i
        // and the back edges of C *and* the visible edges of mGraph-C.
        for (auto const& key : boundary)
        {
            int v0 = key.V[0];
            int v1 = key.V[1];
            if (mQuery.ToLine(i, v0, v1) < 0)
            {
                // This is a back edge of the boundary.
                if (!mGraph.Insert(i, v0, v1))
                {
                    return false;
                }
            }
        }
        for (auto const& key : visible)
        {
            if (!mGraph.Insert(i, key.V[1], key.V[0]))
            {
                return false;
            }
        }
    }

    return true;
}

template <typename InputType, typename ComputeType>
Delaunay2<InputType, ComputeType>::ProcessedVertex::ProcessedVertex()
{
}

template <typename InputType, typename ComputeType>
Delaunay2<InputType, ComputeType>::ProcessedVertex::ProcessedVertex(
    Vector2<InputType> const& inVertex, int inLocation)
    :
    vertex(inVertex),
    location(inLocation)
{
}

template <typename InputType, typename ComputeType>
bool Delaunay2<InputType, ComputeType>::ProcessedVertex::operator<(
    ProcessedVertex const& v) const
{
    return vertex < v.vertex;
}


}