GteIsPlanarGraph.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 <algorithm>
#include <array>
#include <map>
#include <set>
#include <vector>

// Test whether an undirected graph is planar.  The input positions must be
// unique and the input edges must be unique, so the number of positions is
// at least 2 and the number of edges is at least one.  The elements of the
// edges array must be indices in {0,..,positions.size()-1}.
// 
// A sort-and-sweep algorithm is used to determine edge-edge intersections.
// If none of the intersections occur at edge-interior points, the graph is
// planar.  See Game Physics (2nd edition), Section 6.2.2: Culling with
// Axis-Aligned Bounding Boxes for such an algorithm.  The operator()
// returns 'true' when the graph is planar.  If it returns 'false', the
// 'invalidIntersections' set contains pairs of edges that intersect at an
// edge-interior point (that by definition is not a graph vertex).  Each set
// element is a pair of indices into the 'edges' array; the indices are
// ordered so that element[0] < element[1].  The Real type must be chosen
// to guarantee exact computation of edge-edge intersections.

namespace gte
{

template <typename Real>
class IsPlanarGraph
{
public:
    IsPlanarGraph();

    enum
    {
        IPG_IS_PLANAR_GRAPH = 0,
        IPG_INVALID_INPUT_SIZES = 1,
        IPG_DUPLICATED_POSITIONS = 2,
        IPG_DUPLICATED_EDGES = 4,
        IPG_DEGENERATE_EDGES = 8,
        IPG_EDGES_WITH_INVALID_VERTICES = 16,
        IPG_INVALID_INTERSECTIONS = 32
    };

    // The function returns a combination of the IPG_* flags listed in the
    // previous enumeration.  A combined value of 0 indicates the input
    // forms a planar graph.  If the combined value is not zero, you may
    // examine the flags for the failure conditions and use the Get* member
    // accessors to obtain specific information about the failure.  If the
    // positions.size() < 2 or edges.size() == 0, the IPG_INVALID_INPUT_SIZES
    // flag is set.

    int operator()(std::vector<std::array<Real, 2>> const& positions,
        std::vector<std::array<int, 2>> const& edges);

    // A pair of indices (v0,v1) into the positions array is stored as
    // (min(v0,v1), max(v0,v1)).  This supports sorted containers of edges.
    struct OrderedEdge
    {
        OrderedEdge(int v0 = -1, int v1 = -1);
        bool operator<(OrderedEdge const& edge) const;
        std::array<int, 2> V;
    };

    inline std::vector<std::vector<int>> const& GetDuplicatedPositions() const;
    inline std::vector<std::vector<int>> const& GetDuplicatedEdges() const;
    inline std::vector<int> const& GetDegenerateEdges() const;
    inline std::vector<int> const& GetEdgesWithInvalidVertices() const;
    inline std::vector<OrderedEdge> const& GetInvalidIntersections() const;

private:
    class Endpoint
    {
    public:
        Real value;     // endpoint value
        int type;       // '0' if interval min, '1' if interval max.
        int index;      // index of interval containing this endpoint

        // Comparison operator for sorting.
        bool operator<(Endpoint const& endpoint) const;
    };

    int IsValidTopology(std::vector<std::array<Real, 2>> const& positions,
        std::vector<std::array<int, 2>> const& edges);

    void ComputeOverlappingRectangles(std::vector<std::array<Real, 2>> const& positions,
        std::vector<std::array<int, 2>> const& edges,
        std::set<OrderedEdge>& overlappingRectangles) const;

    bool InvalidSegmentIntersection(
        std::array<Real, 2> const& p0, std::array<Real, 2> const& p1,
        std::array<Real, 2> const& q0, std::array<Real, 2> const& q1) const;

    std::vector<std::vector<int>> mDuplicatedPositions;
    std::vector<std::vector<int>> mDuplicatedEdges;
    std::vector<int> mDegenerateEdges;
    std::vector<int> mEdgesWithInvalidVertices;
    std::vector<OrderedEdge> mInvalidIntersections;
    Real mZero, mOne;
};


template <typename Real>
IsPlanarGraph<Real>::IsPlanarGraph()
    :
    mZero(0),
    mOne(1)
{
}

template <typename Real>
int IsPlanarGraph<Real>::operator()(
    std::vector<std::array<Real, 2>> const& positions,
    std::vector<std::array<int, 2>> const& edges)
{
    mDuplicatedPositions.clear();
    mDuplicatedEdges.clear();
    mDegenerateEdges.clear();
    mEdgesWithInvalidVertices.clear();
    mInvalidIntersections.clear();

    int flags = IsValidTopology(positions, edges);
    if (flags == IPG_INVALID_INPUT_SIZES)
    {
        return flags;
    }

    std::set<OrderedEdge> overlappingRectangles;
    ComputeOverlappingRectangles(positions, edges, overlappingRectangles);
    for (auto key : overlappingRectangles)
    {
        // Get the endpoints of the line segments for the edges whose bounding
        // rectangles overlapped.  Determine whether the line segments
        // intersect.  If they do, determine how they intersect.
        std::array<int, 2> e0 = edges[key.V[0]];
        std::array<int, 2> e1 = edges[key.V[1]];
        std::array<Real, 2> const& p0 = positions[e0[0]];
        std::array<Real, 2> const& p1 = positions[e0[1]];
        std::array<Real, 2> const& q0 = positions[e1[0]];
        std::array<Real, 2> const& q1 = positions[e1[1]];
        if (InvalidSegmentIntersection(p0, p1, q0, q1))
        {
            mInvalidIntersections.push_back(key);
        }
    }

    if (mInvalidIntersections.size() > 0)
    {
        flags |= IPG_INVALID_INTERSECTIONS;
    }

    return flags;
}

template <typename Real>
inline std::vector<std::vector<int>> const&
IsPlanarGraph<Real>::GetDuplicatedPositions() const
{
    return mDuplicatedPositions;
}

template <typename Real>
inline std::vector<std::vector<int>> const&
IsPlanarGraph<Real>::GetDuplicatedEdges() const
{
    return mDuplicatedEdges;
}

template <typename Real>
inline std::vector<int> const&
IsPlanarGraph<Real>::GetDegenerateEdges() const
{
    return mDegenerateEdges;
}

template <typename Real>
inline std::vector<int> const&
IsPlanarGraph<Real>::GetEdgesWithInvalidVertices() const
{
    return mEdgesWithInvalidVertices;
}

template <typename Real>
inline std::vector<typename IsPlanarGraph<Real>::OrderedEdge> const&
IsPlanarGraph<Real>::GetInvalidIntersections() const
{
    return mInvalidIntersections;
}

template <typename Real>
int IsPlanarGraph<Real>::IsValidTopology(std::vector<std::array<Real, 2>> const& positions,
    std::vector<std::array<int, 2>> const& edges)
{
    int const numPositions = static_cast<int>(positions.size());
    int const numEdges = static_cast<int>(edges.size());
    if (numPositions < 2 || numEdges == 0)
    {
        // The graph must have at least one edge.
        return IPG_INVALID_INPUT_SIZES;
    }

    // The positions must be unique.
    int flags = IPG_IS_PLANAR_GRAPH;
    std::map<std::array<Real, 2>, std::vector<int>> uniquePositions;
    for (int i = 0; i < numPositions; ++i)
    {
        std::array<Real, 2> p = positions[i];
        auto iter = uniquePositions.find(p);
        if (iter == uniquePositions.end())
        {
            std::vector<int> indices;
            indices.push_back(i);
            uniquePositions.insert(std::make_pair(p, indices));
        }
        else
        {
            iter->second.push_back(i);
        }
    }
    if (uniquePositions.size() < positions.size())
    {
        // At least two positions are duplicated.
        for (auto const& element : uniquePositions)
        {
            if (element.second.size() > 1)
            {
                mDuplicatedPositions.push_back(element.second);
            }
        }
        flags |= IPG_DUPLICATED_POSITIONS;
    }

    // The edges must be unique.
    std::map<OrderedEdge, std::vector<int>> uniqueEdges;
    for (int i = 0; i < numEdges; ++i)
    {
        OrderedEdge key(edges[i][0], edges[i][1]);
        auto iter = uniqueEdges.find(key);
        if (iter == uniqueEdges.end())
        {
            std::vector<int> indices;
            indices.push_back(i);
            uniqueEdges.insert(std::make_pair(key, indices));
        }
        else
        {
            iter->second.push_back(i);
        }
    }
    if (uniqueEdges.size() < edges.size())
    {
        // At least two edges are duplicated, possibly even a pair of
        // edges (v0,v1) and (v1,v0) which is not allowed because the
        // graph is undirected.
        for (auto const& element : uniqueEdges)
        {
            if (element.second.size() > 1)
            {
                mDuplicatedEdges.push_back(element.second);
            }
        }
        flags |= IPG_DUPLICATED_EDGES;
    }

    // The edges are represented as pairs of indices into the 'positions'
    // array.  The indices for a single edge must be different (no edges
    // allowed from a vertex to itself) and all indices must be valid.
    // At the same time, keep track of unique edges.
    for (int i = 0; i < numEdges; ++i)
    {
        std::array<int, 2> e = edges[i];
        if (e[0] == e[1])
        {
            // The edge is degenerate, originating and terminating at the
            // same vertex.
            mDegenerateEdges.push_back(i);
            flags |= IPG_DEGENERATE_EDGES;
        }

        if (e[0] < 0 || e[0] >= numPositions || e[1] < 0 || e[1] >= numPositions)
        {
            // The edge has an index that references a non-existant vertex.
            mEdgesWithInvalidVertices.push_back(i);
            flags |= IPG_EDGES_WITH_INVALID_VERTICES;
        }
    }

    return flags;
}

template <typename Real>
void IsPlanarGraph<Real>::ComputeOverlappingRectangles(std::vector<std::array<Real, 2>> const& positions,
    std::vector<std::array<int, 2>> const& edges, std::set<OrderedEdge>& overlappingRectangles) const
{
    // Compute axis-aligned bounding rectangles for the edges.
    int const numEdges = static_cast<int>(edges.size());
    std::vector<std::array<Real, 2>> emin(numEdges);
    std::vector<std::array<Real, 2>> emax(numEdges);
    for (int i = 0; i < numEdges; ++i)
    {
        std::array<int, 2> e = edges[i];
        std::array<Real, 2> const& p0 = positions[e[0]];
        std::array<Real, 2> const& p1 = positions[e[1]];

        for (int j = 0; j < 2; ++j)
        {
            if (p0[j] < p1[j])
            {
                emin[i][j] = p0[j];
                emax[i][j] = p1[j];
            }
            else
            {
                emin[i][j] = p1[j];
                emax[i][j] = p0[j];
            }
        }
    }

    // Get the rectangle endpoints.
    int const numEndpoints = 2 * numEdges;
    std::vector<Endpoint> xEndpoints(numEndpoints);
    std::vector<Endpoint> yEndpoints(numEndpoints);
    for (int i = 0, j = 0; i < numEdges; ++i)
    {
        xEndpoints[j].type = 0;
        xEndpoints[j].value = emin[i][0];
        xEndpoints[j].index = i;
        yEndpoints[j].type = 0;
        yEndpoints[j].value = emin[i][1];
        yEndpoints[j].index = i;
        ++j;

        xEndpoints[j].type = 1;
        xEndpoints[j].value = emax[i][0];
        xEndpoints[j].index = i;
        yEndpoints[j].type = 1;
        yEndpoints[j].value = emax[i][1];
        yEndpoints[j].index = i;
        ++j;
    }

    // Sort the rectangle endpoints.
    std::sort(xEndpoints.begin(), xEndpoints.end());
    std::sort(yEndpoints.begin(), yEndpoints.end());

    // Sweep through the endpoints to determine overlapping x-intervals.  Use
    // an active set of rectangles to reduce the complexity of the algorithm.
    std::set<int> active;
    for (int i = 0; i < numEndpoints; ++i)
    {
        Endpoint const& endpoint = xEndpoints[i];
        int index = endpoint.index;
        if (endpoint.type == 0)  // an interval 'begin' value
        {
            // In the 1D problem, the current interval overlaps with all the
            // active intervals.  In 2D this we also need to check for
            // y-overlap.
            for (auto activeIndex : active)
            {
                // Rectangles activeIndex and index overlap in the
                // x-dimension.  Test for overlap in the y-dimension.
                std::array<Real, 2> const& r0min = emin[activeIndex];
                std::array<Real, 2> const& r0max = emax[activeIndex];
                std::array<Real, 2> const& r1min = emin[index];
                std::array<Real, 2> const& r1max = emax[index];
                if (r0max[1] >= r1min[1] && r0min[1] <= r1max[1])
                {
                    if (activeIndex < index)
                    {
                        overlappingRectangles.insert(OrderedEdge(activeIndex, index));
                    }
                    else
                    {
                        overlappingRectangles.insert(OrderedEdge(index, activeIndex));
                    }
                }
            }
            active.insert(index);
        }
        else  // an interval 'end' value
        {
            active.erase(index);
        }
    }
}

template <typename Real>
bool IsPlanarGraph<Real>::InvalidSegmentIntersection(
    std::array<Real, 2> const& p0, std::array<Real, 2> const& p1,
    std::array<Real, 2> const& q0, std::array<Real, 2> const& q1) const
{
    // We must solve the two linear equations
    //   p0 + t0 * (p1 - p0) = q0 + t1 * (q1 - q0)
    // for the unknown variables t0 and t1.  These may be written as
    //   t0 * (p1 - p0) - t1 * (q1 - q0) = q0 - p0
    // If denom = Dot(p1 - p0, Perp(q1 - q0)) is not zero, then
    //   t0 = Dot(q0 - p0, Perp(q1 - q0)) / denom = numer0 / denom
    //   t1 = Dot(q0 - p0, Perp(p1 - p0)) / denom = numer1 / denom
    // For an invalid intersection, we need (t0,t1) with:
    // ((0 < t0 < 1) and (0 <= t1 <= 1)) or ((0 <= t0 <= 1) and (0 < t1 < 1).
    // If denom is zero, then

    std::array<Real, 2> p1mp0, q1mq0, q0mp0;
    for (int j = 0; j < 2; ++j)
    {
        p1mp0[j] = p1[j] - p0[j];
        q1mq0[j] = q1[j] - q0[j];
        q0mp0[j] = q0[j] - p0[j];
    }

    Real denom = p1mp0[0] * q1mq0[1] - p1mp0[1] * q1mq0[0];
    Real numer0 = q0mp0[0] * q1mq0[1] - q0mp0[1] * q1mq0[0];
    Real numer1 = q0mp0[0] * p1mp0[1] - q0mp0[1] * p1mp0[0];

    if (denom != mZero)
    {
        // The lines of the segments are not parallel.
        if (denom > mZero)
        {
            if (mZero <= numer0 && numer0 <= denom && mZero <= numer1 && numer1 <= denom)
            {
                // The segments intersect.
                return (numer0 != mZero && numer0 != denom) || (numer1 != mZero && numer1 != denom);
            }
            else
            {
                return false;
            }
        }
        else  // denom < mZero
        {
            if (mZero >= numer0 && numer0 >= denom && mZero >= numer1 && numer1 >= denom)
            {
                // The segments intersect.
                return (numer0 != mZero && numer0 != denom) || (numer1 != mZero && numer1 != denom);
            }
            else
            {
                return false;
            }
        }
    }
    else
    {
        // The lines of the segments are parallel.
        if (numer0 != mZero || numer1 != mZero)
        {
            // The lines of the segments are separated.
            return false;
        }
        else
        {
            // The segments lie on the same line.  Compute the parameter
            // intervals for the segments in terms of the t0-parameter and
            // determine their overlap (if any).
            std::array<Real, 2> q1mp0;
            for (int j = 0; j < 2; ++j)
            {
                q1mp0[j] = q1[j] - p0[j];
            }
            Real sqrLenP1mP0 = p1mp0[0] * p1mp0[0] + p1mp0[1] * p1mp0[1];
            Real value0 = q0mp0[0] * p1mp0[0] + q0mp0[1] * p1mp0[1];
            Real value1 = q1mp0[0] * p1mp0[0] + q1mp0[1] * p1mp0[1];
            if ((value0 >= sqrLenP1mP0 && value1 >= sqrLenP1mP0)
                || (value0 <= mZero && value1 <= mZero))
            {
                // If the segments intersect, they must do so at endpoints of
                // the segments.
                return false;
            }
            else
            {
                // The segments overlap in a positive-length interval.
                return true;
            }
        }
    }
}

template <typename Real>
bool IsPlanarGraph<Real>::Endpoint::operator<(Endpoint const& endpoint) const
{
    if (value < endpoint.value)
    {
        return true;
    }
    if (value > endpoint.value)
    {
        return false;
    }
    return type < endpoint.type;
}

template <typename Real>
IsPlanarGraph<Real>::OrderedEdge::OrderedEdge(int v0, int v1)
{
    if (v0 < v1)
    {
        // v0 is minimum
        V[0] = v0;
        V[1] = v1;
    }
    else
    {
        // v1 is minimum
        V[0] = v1;
        V[1] = v0;
    }
}

template <typename Real>
bool IsPlanarGraph<Real>::OrderedEdge::operator<(OrderedEdge const& edge) const
{
    // Lexicographical ordering used by std::array<int,2>.
    return V < edge.V;
}

}