diff options
Diffstat (limited to 'src/2geom/planar-graph.h')
| -rw-r--r-- | src/2geom/planar-graph.h | 1252 |
1 files changed, 1252 insertions, 0 deletions
diff --git a/src/2geom/planar-graph.h b/src/2geom/planar-graph.h new file mode 100644 index 0000000..fb5f1ac --- /dev/null +++ b/src/2geom/planar-graph.h @@ -0,0 +1,1252 @@ +/** @file PlanarGraph – a graph geometrically embedded in the plane. + */ +/* + * Authors: + * Rafał Siejakowski <rs@rs-math.net> + * + * Copyright 2022 the Authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + +// WARNING: This is a private header. Do not include it directly. + +#ifndef LIB2GEOM_SEEN_PLANAR_GRAPH_H +#define LIB2GEOM_SEEN_PLANAR_GRAPH_H + +#include <algorithm> +#include <iterator> +#include <list> + +#include <2geom/angle.h> +#include <2geom/coord.h> +#include <2geom/line.h> +#include <2geom/point.h> +#include <2geom/path.h> +#include <2geom/path-intersection.h> +#include <2geom/utils.h> + +namespace Geom { + +/** + * \class PlanarGraph + * \brief Planar graph - a graph geometrically embedded in the plane. + * + * A planar graph is composed of vertices with assigned locations (as points in the plane) + * and of edges (arcs), which are imagined as non-intersecting paths in the plane connecting + * the vertices. The edges can hold user-supplied labels (e.g., weights) which support event + * callbacks for when the graph is reconfigured, allowing the labels to be updated accordingly. + * + * \tparam EdgeLabel A user-supplied type; an object of this type will be attached to each + * edge of the planar graph (e.g., a "weight" of the edge). The type must + * satisfy requirements described further below. + * + * In order to construct a planar graph, you should specify the clumping precision (passed as + * a constructor argument) and then use the method insertEdge() to add edges to the graph, as + * many times as necessary. The graph will automatically figure out the locations of the + * vertices based on the endpoints of the inserted edges. Vertices will be combined into one + * when they are positioned within the distance specified as the clumping threshold, and the + * inserted edges will be attached to them accordingly. It is also possible to insert paths + * (typically, closed) not attached to any vertices, using the method insertDetached(). + * + * After the edges are inserted, the graph is in a potentially degenerate state, where several + * edges may exactly coincide in part or in full. If this is not desired, you can regularize + * the graph by calling regularize(). During the regularization process, any overlapping edges + * are combined into one. Partially overlapping edges are first split into overlapping and + * non-overlapping portions, after which the overlapping portions are combined. If the edges + * or their parts overlap but run in opposite directions, one of them will be reversed before + * being merged with the other one. The overlaps are detected using the precision setting passed + * as the clumping precision in the constructor argument. + * + * Note however that the regularization procedure does NOT detect transverse intersections + * between the edge paths: if such intersections are not desired, the user must pass non-\ + * intersecting paths to the insertEdge() method (the paths may still have common endpoints, + * which is fine: that's how common vertices are created). + * + * The insertion of new edges invalidates the regularized status, which you can check at any + * time by calling isRegularized(). + * + * The vertices stored by the graph are sorted by increasing X-coordinate, and if they have + * equal X-coordinates, by increasing Y-coordinate. Even before regularization, incidences of + * edges to each vertex are sorted by increasing azimuth of the outgoing tangent (departure + * heading, but in radians, in the interval \f$(-\pi, \pi]\f$). After regularization, the edges + * around each vertex are guaranteed to be sorted counterclockwise (when the Y-axis points up) + * by where they end up going eventually, even if they're tangent at the vertex and therefore + * have equal or nearly equal departure azimuths. + * + * \note + * Requirements on the \c EdgeLabel template parameter type. + * In order for the template to instantiate correctly, the following must be satisfied: + * \li The \c EdgeLabel type provides a method \c onReverse() which gets called whenever + * the orientation of the labeled edge is reversed. This is useful when implementing + * a directed graph, since the label can keep track of the logical direction. + * \li The \c EdgeLabel type provides a method \c onMergeWith(EdgeLabel const&) which gets + * called when the labeled edge is combined with a geometrically identical (coinciding) + * edge (both combined edges having the same orientations). The label of the edge merged + * with the current one is provided as an argument to the method. This is useful when + * implementing a graph with weights: for example, when two edges are merged, you may + * want to combine their weights in some way. + * \li There is a method \c onDetach() called when the edge is removed from the graph. The + * edge objects are never destroyed but may be disconnected from the graph when they're no + * longer needed; this allows the user to put the labels of such edges in a "dead" state. + * \li The \c EdgeLabel objects must be copy-constructible and copy-assignable. This is + * because when an edge is subdivided into two, the new edges replacing it get decorated + * with copies of the original edge's label. + */ +template<typename EdgeLabel> +#if __cplusplus >= 202002L +requires requires(EdgeLabel el, EdgeLabel const &other) { + el.onReverse(); + el.onMergeWith(other); + el.onDetach(); + el = other; +} +#endif +class PlanarGraph +{ +public: + + /** Represents the joint between an edge and a vertex. */ + struct Incidence + { + using Sign = bool; + inline static Sign const START = false; + inline static Sign const END = true; + + double azimuth; ///< Angle of the edge's departure. + unsigned index; ///< Index of the edge in the parent graph. + Sign sign; ///< Whether this is the start or end of the edge. + bool invalid = false; ///< Whether this incidence has been marked for deletion. + + Incidence(unsigned edge_index, double departure_azimuth, Sign which_end) + : azimuth{departure_azimuth} + , index{edge_index} + , sign{which_end} + { + } + ~Incidence() = default; + + /// Compare incidences based on their azimuths in radians. + inline bool operator<(Incidence const &other) const { return azimuth < other.azimuth; } + + /// Compare the azimuth of an incidence with the given angle. + inline bool operator<(double angle) const { return azimuth < angle; } + + /// Check equality (only tests edges and their ends). + inline bool operator==(Incidence const &other) const + { + return index == other.index && sign == other.sign; + } + }; + using IncIt = typename std::list<Incidence>::iterator; + using IncConstIt = typename std::list<Incidence>::const_iterator; + + /** Represents the vertex of a planar graph. */ + class Vertex + { + private: + Point const _position; ///< Geometric position of the vertex. + std::list<Incidence> _incidences; ///< List of incidences of edges to this vertex. + unsigned mutable _flags = 0; ///< User-settable flags. + + inline static Point::LexLess<X> const _cmp; ///< Point sorting function object. + + public: + Vertex(Point const &pos) + : _position{pos} + { + } + + /** Get the geometric position of the vertex. */ + Point const &point() const { return _position; } + + /** Get the list of incidences to the vertex. */ + auto const &getIncidences() const { return _incidences; } + + /** Compare vertices based on their coordinates (lexicographically). */ + bool operator<(Vertex const &other) const { return _cmp(_position, other._position); } + + unsigned flags() const { return _flags; } ///< Get the user flags. + void setFlags(unsigned flags) const { _flags = flags; } ///< Set the user flags. + + /** Get the cyclic-next incidence after the passed one, in the CCW direction. */ + IncConstIt cyclicNextIncidence(IncConstIt it) const { return cyclic_next(it, _incidences); } + + /** Get the cyclic-next incidence after the passed one, in the CW direction. */ + IncConstIt cyclicPrevIncidence(IncConstIt it) const { return cyclic_prior(it, _incidences); } + + /** The same but with pointers. */ + Incidence *cyclicNextIncidence(Incidence *from) + { + return &(*cyclicNextIncidence(_incidencePtr2It(from))); + } + Incidence *cyclicPrevIncidence(Incidence *from) + { + return &(*cyclicPrevIncidence(_incidencePtr2It(from))); + } + + private: + /** Same as above, but not const (so only for private use). */ + IncIt cyclicNextIncidence(IncIt it) { return cyclic_next(it, _incidences); } + IncIt cyclicPrevIncidence(IncIt it) { return cyclic_prior(it, _incidences); } + + /** Insert an incidence; for internal use by the PlanarGraph class. */ + Incidence &_addIncidence(unsigned edge_index, double azimuth, typename Incidence::Sign sign) + { + auto where = std::find_if(_incidences.begin(), _incidences.end(), [=](auto &inc) -> bool { + return inc.azimuth >= azimuth; + }); + return *(_incidences.emplace(where, edge_index, azimuth, sign)); + } + + /** Return a valid iterator to an incidence passed by pointer; + * if the pointer is invalid, return a start iterator. */ + IncIt _incidencePtr2It(Incidence *pointer) + { + auto it = std::find_if(_incidences.begin(), _incidences.end(), + [=](Incidence const &i) -> bool { return &i == pointer; }); + return (it == _incidences.end()) ? _incidences.begin() : it; + } + + friend class PlanarGraph<EdgeLabel>; + }; + using VertexIterator = typename std::list<Vertex>::iterator; + + /** Represents an edge of the planar graph. */ + struct Edge + { + Vertex *start = nullptr, *end = nullptr; ///< Start and end vertices. + Path path; ///< The path associated to the edge. + bool detached = false; ///< Whether the edge is detached from the graph. + bool inserted_as_detached = false; ///< Whether the edge was inserted as detached. + EdgeLabel mutable label; ///< The user-supplied label of the edge. + + /** Construct an edge with a given label. */ + Edge(Path &&movein_path, EdgeLabel &&movein_label) + : path{movein_path} + , label{movein_label} + { + } + + /** Detach the edge from the graph. */ + void detach() + { + detached = true; + label.onDetach(); + } + }; + using EdgeIterator = typename std::vector<Edge>::iterator; + using EdgeConstIterator = typename std::vector<Edge>::const_iterator; + +private: + double const _precision; ///< Numerical epsilon for vertex clumping. + std::list<Vertex> _vertices; ///< Vertices of the graph. + std::vector<Edge> _edges; ///< Edges of the graph. + std::vector< std::pair<Vertex *, Incidence *> > _junk; ///< Incidences that should be purged. + bool _regularized = true; // An empty graph is (trivially) regularized. + +public: + PlanarGraph(Coord precision = EPSILON) + : _precision{precision} + { + } + + std::list<Vertex> const &getVertices() const { return _vertices; } + std::vector<Edge> const &getEdges() const { return _edges; } + Edge const &getEdge(size_t index) const { return _edges.at(index); } + size_t getEdgeIndex(Edge const &edge) const { return &edge - _edges.data(); } + double getPrecision() const { return _precision; } + size_t numVertices() const { return _vertices.size(); } + size_t numEdges(bool include_detached = true) const + { + if (include_detached) { + return _edges.size(); + } + return std::count_if(_edges.begin(), _edges.end(), + [](Edge const &e) -> bool { return !e.detached; }); + } + + /** Check if the graph has been regularized. */ + bool isRegularized() const { return _regularized; } + + // 0x1p-50 is about twice the distance between M_PI and the next representable double. + void regularize(double angle_precision = 0x1p-50, bool remove_collapsed_loops = true); + + /** Allocate memory to store the specified number of edges. */ + void reserveEdgeCapacity(size_t capacity) { _edges.reserve(capacity); } + + unsigned insertEdge(Path &&path, EdgeLabel &&edge = EdgeLabel()); + unsigned insertDetached(Path &&path, EdgeLabel &&edge = EdgeLabel()); + + /** Edge insertion with a copy of the path. */ + unsigned insertEdge(Path const &path, EdgeLabel &&edge = EdgeLabel()) + { + return insertEdge(Path(path), std::forward<EdgeLabel>(edge)); + } + unsigned insertDetached(Path const &path, EdgeLabel &&edge = EdgeLabel()) + { + return insertDetached(Path(path), std::forward<EdgeLabel>(edge)); + } + + /** \brief Find the incidence at the specified endpoint of the edge. + * + * \param edge_index The index of the edge whose incidence we wish to query. + * \param sign Which end of the edge do we want an incidence of. + * \return A pair consisting of pointers to the vertex and the incidence. + * If not found, both pointers will be null. + */ + std::pair<Vertex *, Incidence *> + getIncidence(unsigned edge_index, typename Incidence::Sign sign) const + { + if (edge_index >= _edges.size() || _edges[edge_index].detached) { + return {nullptr, nullptr}; + } + Vertex *vertex = (sign == Incidence::START) ? _edges[edge_index].start + : _edges[edge_index].end; + if (!vertex) { + return {nullptr, nullptr}; + } + auto it = std::find(vertex->_incidences.begin(), vertex->_incidences.end(), + Incidence(edge_index, 42, sign)); // azimuth doesn't matter. + if (it == vertex->_incidences.end()) { + return {nullptr, nullptr}; + } + return {vertex, &(*it)}; + } + + /** + * \brief Go clockwise or counterclockwise around the vertex and find the next incidence. + * The notions of "clockwise"/"counterclockwise" correspond to the y-axis pointing up. + * + * \param vertex The vertex around which to orbit. + * \param incidence The incidence from which to start traversal. + * \param clockwise Whether to go clockwise instead of (default) counterclockwise. + * \return The next incidence encountered going in the specified direction. + */ + inline Incidence const &nextIncidence(VertexIterator const &vertex, IncConstIt const &incidence, + bool clockwise = false) const + { + return clockwise ? *(vertex->_cyclicPrevIncidence(incidence)) + : *(vertex->_cyclicNextIncidence(incidence)); + } + + /** As above, but taking references instead of iterators. */ + inline Incidence const &nextIncidence(Vertex const &vertex, Incidence const &incidence, + bool clockwise = false) const + { + IncConstIt it = std::find(vertex._incidences.begin(), vertex._incidences.end(), incidence); + if (it == vertex._incidences.end()) { + return incidence; + } + return clockwise ? *(vertex.cyclicPrevIncidence(it)) + : *(vertex.cyclicNextIncidence(it)); + } + + /** As above, but return an iterator to a const incidence. */ + inline IncConstIt nextIncidenceIt(Vertex const &vertex, Incidence const &incidence, + bool clockwise = false) const + { + IncConstIt it = std::find(vertex._incidences.begin(), vertex._incidences.end(), incidence); + if (it == vertex._incidences.end()) { + return vertex._incidences.begin(); + } + return clockwise ? vertex.cyclicPrevIncidence(it) + : vertex.cyclicNextIncidence(it); + } + inline IncConstIt nextIncidenceIt(Vertex const &vertex, IncConstIt const &incidence, + bool clockwise = false) const + { + return clockwise ? vertex.cyclicPrevIncidence(incidence) + : vertex.cyclicNextIncidence(incidence); + } + + /** As above, but start at the prescribed departure azimuth (in radians). + * + * \return A pointer to the incidence emanating from the vertex at or immediately after + * the specified azimuth, when going around the vertex in the specified direction. + * If the vertex has no incidences, return value is nullptr. + */ + Incidence *nextIncidence(VertexIterator const &vertex, double azimuth, + bool clockwise = false) const; + + /** Get the incident path, always oriented away from the vertex. */ + Path getOutgoingPath(Incidence const *incidence) const + { + return incidence ? _getPathImpl(incidence, Incidence::START) : Path(); + } + + /** Get the incident path, always oriented towards the vertex. */ + Path getIncomingPath(Incidence const *incidence) const + { + return incidence ? _getPathImpl(incidence, Incidence::END) : Path(); + } + +private: + inline Path _getPathImpl(Incidence const *incidence, typename Incidence::Sign origin) const + { + return (incidence->sign == origin) ? _edges[incidence->index].path + : _edges[incidence->index].path.reversed(); + } + + /** Earmark an incidence for future deletion. */ + inline void _throwAway(Vertex *vertex, Incidence *incidence) + { + if (!vertex || !incidence) { + return; + } + incidence->invalid = true; + _junk.emplace_back(vertex, incidence); + } + + // Topological reconfiguration functions; see their definitions for documentation. + bool _compareAndReglue(Vertex &vertex, Incidence *first, Incidence *second, bool deloop); + Vertex *_ensureVertexAt(Point const &position); + void _mergeCoincidingEdges(Incidence *first, Incidence *second); + void _mergeShorterLonger(Vertex &vertex, Incidence *shorter, Incidence *longer, + PathTime const &time_on_longer); + void _mergeWyeConfiguration(Vertex &vertex, Incidence *first, Incidence *second, + PathIntersection const &split); + void _purgeJunkIncidences(); + void _reglueLasso(Vertex &vertex, Incidence *first, Incidence *second, + PathIntersection const &split); + bool _reglueTeardrop(Vertex &vertex, Incidence *first, Incidence *second, bool deloop); + void _reglueTangentFan(Vertex &vertex, IncIt const &first, IncIt const &last, bool deloop); + void _regularizeVertex(Vertex &vertex, double angle_precision, bool deloop); + + // === Static stuff === + + /** Return the angle between the vector and the positive X axis or 0 if undefined. */ + inline static double _getAzimuth(Point const &vec) { return vec.isZero() ? 0.0 : atan2(vec); } + + /** Return path time corresponding to the same point on the reversed path. */ + inline static PathTime _reversePathTime(PathTime const &time, Path const &path) + { + int new_index = path.size() - time.curve_index - 1; + Coord new_time = 1.0 - time.t; + if (new_index < 0) { + new_index = 0; + new_time = 0; + } + return PathTime(new_index, new_time); + } + + /** Return path time at the end of the path. */ + inline static PathTime _pathEnd(Path const &path) { return PathTime(path.size() - 1, 1.0); } + inline static auto const PATH_START = PathTime(0, 0); + +public: + static double closedPathArea(Path const &path); + static bool deviatesLeft(Path const &first, Path const &second); +}; + +/** + * \brief Insert a new vertex or reuse an existing one. + * + * Ensures that there is a vertex at or near the specified position + * (within the distance of _precision). + * + * \param pos The desired geometric position of the new vertex. + * \return A pointer to the inserted vertex or a pre-existing vertex near the + * desired position. + */ +template<typename EL> +typename PlanarGraph<EL>::Vertex *PlanarGraph<EL>::_ensureVertexAt(Point const &pos) +{ + auto const insert_at_front = [&, this]() -> Vertex* { + _vertices.emplace_front(pos); + return &(_vertices.front()); + }; + + if (_vertices.empty()) { + return insert_at_front(); + } + + // TODO: Use a heap? + auto it = std::find_if(_vertices.begin(), _vertices.end(), [&](Vertex const &v) -> bool { + return Vertex::_cmp(pos, v._position); // existing vertex position > pos. + }); + + if (it != _vertices.end()) { + if (are_near(pos, it->_position, _precision)) { + return &(*it); // Reuse existing vertex. + } + if (it == _vertices.begin()) { + return insert_at_front(); + } + } + // Look at the previous element, reuse if near, insert before `it` otherwise. + return &(*(are_near(pos, std::prev(it)->_position, _precision) ? std::prev(it) + : _vertices.emplace(it, pos))); +} + +/** + * \brief Move-insert a new labeled edge into the planar graph. + * + * \param path The geometric path of the edge. + * \param label Optionally, the label (extra user data) associated to this edge. + * If absent, a default-constructed label will be used. + * \return The index of the inserted edge. + */ +template<typename EdgeLabel> +unsigned PlanarGraph<EdgeLabel>::insertEdge(Path &&path, EdgeLabel &&label) +{ + unsigned edge_index = _edges.size(); + auto &inserted = _edges.emplace_back(std::forward<Path>(path), + std::forward<EdgeLabel>(label)); + + // Calculate the outgoing azimuths at both endpoints. + double const start_azimuth = _getAzimuth(inserted.path.initialUnitTangent()); + double const end_azimuth = _getAzimuth(-inserted.path.finalUnitTangent()); + + // Get the endpoints into the graph. + auto start = _ensureVertexAt(inserted.path.initialPoint()); + auto end = _ensureVertexAt(inserted.path.finalPoint()); + + // Inform the edge about its endpoints. + inserted.start = start; + inserted.end = end; + + // Add incidences at the endpoints. + start->_addIncidence(edge_index, start_azimuth, Incidence::START); + end->_addIncidence(edge_index, end_azimuth, Incidence::END); + + _regularized = false; + return edge_index; +} + +/** + * \brief Move-insert a new labeled edge but do not connect it to the graph. + * + * Although the graph will hold the path data of an edge inserted in this way, the edge + * will not be connected to any vertex. This can be used to store information about closed + * paths (loops) in the instance, without having to specify starting points for them. + * + * \param path The geometric path of the edge. + * \param label Optionally, the label (extra user data) associated to this edge; if absent, + * the label will be default-constructed. + * \return The index of the inserted edge. + */ +template<typename EdgeLabel> +unsigned PlanarGraph<EdgeLabel>::insertDetached(Path &&path, EdgeLabel &&label) +{ + unsigned edge_index = _edges.size(); + auto &inserted = _edges.emplace_back(std::forward<Path>(path), + std::forward<EdgeLabel>(label)); + inserted.detached = true; + inserted.inserted_as_detached = true; + return edge_index; +} + +/** Remove incidences previously marked as junk. */ +template<typename EdgeLabel> +void PlanarGraph<EdgeLabel>::_purgeJunkIncidences() +{ + for (auto &[vertex, incidence] : _junk) { + Incidence *to_remove = incidence; + auto it = std::find_if(vertex->_incidences.begin(), vertex->_incidences.end(), + [=](Incidence const &inc) -> bool { return &inc == to_remove; }); + if (it != vertex->_incidences.end()) { + vertex->_incidences.erase(it); + } + } + _junk.clear(); +} + +/** + * \brief Merge overlapping edges or their portions, adding vertices if necessary. + * + * \param angle_precision The numerical epsilon for radian angle comparisons. + * \param remove_collapsed_loops Whether to detach edges with both ends incident to the same + * vertex (loops) when these loops don't enclose any area. + * + * This function performs the following operations: + * \li Edges that are tangent at a vertex but don't otherwise overlap are sorted correctly + * in the counterclockwise cyclic order around the vertex. + * \li Degenerate loops which don't enclose any area are removed if the argument is true. + * \li Edges that coincide completely are reversed if needed and merged into one. + * \li Edges that coincide partially are split into overlapping and non-overlapping portions. + * Any overlapping portions are oriented consistently and then merged. + * \li As a sub-case of the above, any non-degenerate loop with an initial self-everlap + * (a "lasso") is replaced with a shorter non-overlapping loop and a simple path leading + * to it. + */ +template<typename EdgeLabel> +void PlanarGraph<EdgeLabel>::regularize(double angle_precision, bool remove_collapsed_loops) +{ + for (auto it = _vertices.begin(); it != _vertices.end(); ++it) { + // Note: the list of vertices may grow during the execution of this loop, + // so don't replace it with a range-for (which stores the end iterator). + // We want the loop to carry on going over the elements it inserted. + if (it->_incidences.size() < 2) { + continue; + } + _regularizeVertex(*it, angle_precision, remove_collapsed_loops); + } + _purgeJunkIncidences(); + _regularized = true; +} + +/** + * \brief Analyze and regularize all edges emanating from a given vertex. + * + * This function goes through the list of incidences at the vertex (roughly sorted by + * azimuth, i.e., departure heading in radians), picking out runs of mutually tangent + * edges and calling _reglueTangentFan() on each run. The algorithm is quite complicated + * because the incidences have to be treated as a cyclic container and a run of mutually + * tangent edges may straddle the "end" of the list, including the possibility that the + * entire list is a single such run. + * + * \param vertex The vertex whose incidences should be analyzed. + * \param angle_precision The numerical epsilon for radian angle comparisons. + * \param deloop Whether loops that don't enclose any area should be detached. + */ +template<typename EdgeLabel> +void PlanarGraph<EdgeLabel>::_regularizeVertex(typename PlanarGraph<EdgeLabel>::Vertex &vertex, + double angle_precision, bool deloop) +{ + auto &incidences = vertex._incidences; + + /// Compare two polar angles in the interval [-π, π] modulo 2π to within angle_precision: + auto const angles_equal = [=](double az1, double az2) -> bool { + static double const twopi = 2.0 * M_PI; + return are_near(std::fmod(az1 + twopi, twopi), std::fmod(az2 + twopi, twopi), + angle_precision); + }; + + IncIt run_begin; // First element in the last discovered run of equal azimuths. + + /// Find and reglue runs of nearly identical azimuths in the specified range. + auto const process_runs = [&](IncIt begin, IncIt end) -> bool + { + double current_azimuth = 42; // Invalid radian angle. + bool in_a_run = false; + + for (auto it = begin; it != end; ++it) { + bool const equal = angles_equal(it->azimuth, current_azimuth); + if (equal && !in_a_run) { + run_begin = std::prev(it); // Save to enclosing scope. + in_a_run = true; + } else if (!equal && in_a_run) { + _reglueTangentFan(vertex, run_begin, std::prev(it), deloop); + in_a_run = false; + } + current_azimuth = it->azimuth; + } + return in_a_run; + }; + + double const last_azimuth = incidences.back().azimuth; + + if (angles_equal(incidences.front().azimuth, last_azimuth)) { + // The cyclic list contains a run of equal azimuths which straddles the "end". + // This means that we must skip the part of this run on the "begin" side on the + // first pass and handle it once we've traversed the remainder of the list. + + bool processed = false; ///< Whether we've cleared the straddling run. + double previous_azimuth = last_azimuth; + IncIt straddling_run_last; + + for (auto it = incidences.begin(); it != incidences.end(); ++it) { + if (!angles_equal(it->azimuth, previous_azimuth)) { + straddling_run_last = std::prev(it); + process_runs(it, incidences.end()); + processed = true; + break; + } + previous_azimuth = it->azimuth; + } + if (processed) { + // Find the first element of the straddling run. + auto it = std::prev(incidences.end()); + while (angles_equal(it->azimuth, last_azimuth)) { + --it; + } + ++it; // Now we're at the start of the straddling run. + _reglueTangentFan(vertex, it, straddling_run_last, deloop); + } else { + // We never encountered anything outside of the straddling run: reglue everything. + _reglueTangentFan(vertex, incidences.begin(), std::prev(incidences.end()), deloop); + } + } else if (process_runs(incidences.begin(), incidences.end())) { + // Our run got rudely interrupted by the end of the container; reglue till the end. + _reglueTangentFan(vertex, run_begin, std::prev(incidences.end()), deloop); + } +} + +/** + * \brief Regularize a fan of mutually tangent edges emanating from a vertex. + * + * This function compares the tangent edges pairwise and ensures that the sequence of their + * incidences to the vertex ends up being sorted by the ultimate direction in which the + * emanating edges fan out, in the counterclockwise order. + * + * If a partial or complete overlap between edges is detected, these edges are reglued. + * + * \param vertex The vertex from which the fan emanates. + * \param first An iterator pointing to the first incidence in the fan. + * \param last An iterator pointing to the last incidence in the fan. + * NOTE: This iterator must point to the actual last incidence, not "past" it. + * The reason is that we're iterating over a cyclic collection, so there + * isn't really a meaningful end. + * \param deloop Whether loops that don't enclose any area should be detached. + */ +template<typename EL> +void PlanarGraph<EL>::_reglueTangentFan(typename PlanarGraph<EL>::Vertex &vertex, + typename PlanarGraph<EL>::IncIt const &first, + typename PlanarGraph<EL>::IncIt const &last, bool deloop) +{ + // Search all pairs (triangular pattern), skipping invalid incidences. + for (auto it = first; it != last; it = vertex.cyclicNextIncidence(it)) { + if (it->invalid) { + continue; + } + for (auto is = vertex.cyclicNextIncidence(it); true; is = vertex.cyclicNextIncidence(is)) { + if (!is->invalid && _compareAndReglue(vertex, &(*it), &(*is), deloop)) { + // Swap the incidences, effectively implementing "bubble sort". + std::swap(*it, *is); + } + if (is == last) { + break; + } + } + } +} + +/** + * \brief Compare a pair of edges emanating from the same vertex in the same direction. + * + * If the edges overlap in part or in full, they get reglued, which means that the topology + * of the graph may get modified. Otherwise, if the detailed comparison shows that the edges + * aren't correctly ordered around the vertex (because the second edge deviates to the right + * instead of to the left of the first, when looking away from the vertex), then the function + * will return true, signalling that the incidences should be swapped. + * + * \param vertex The vertex where the mutually tangent paths meet. + * \param first The incidence appearing as the first one in the provisional cyclic order. + * \param second The incidence appearing as the second one in the provisional cyclic order. + * \param deloop Whether to detach collapsed loops (backtracks) which don't enclose any area. + * \return Whether the incidences should be swapped. + */ +template<typename EL> +bool PlanarGraph<EL>::_compareAndReglue(typename PlanarGraph<EL>::Vertex &vertex, + typename PlanarGraph<EL>::Incidence *first, + typename PlanarGraph<EL>::Incidence *second, bool deloop) +{ + if (first->index == second->index) { + return _reglueTeardrop(vertex, first, second, deloop); + } + + // Get paths corresponding to the edges but travelling away from the vertex. + auto first_path_out = getOutgoingPath(first); + auto second_path_out = getOutgoingPath(second); + auto split = parting_point(first_path_out, second_path_out, _precision); + + if (are_near(split.point(), vertex.point(), _precision)) { + // Paths deviate immediately, so no gluing is needed. The incidences should + // be swapped if the first edge path departs to the left of the second one. + return deviatesLeft(first_path_out, second_path_out); + } + + // Determine the nature of the initial overlap between the paths. + bool const till_end_of_1st = are_near(split.point(), first_path_out.finalPoint(), _precision); + bool const till_end_of_2nd = are_near(split.point(), second_path_out.finalPoint(), _precision); + + if (till_end_of_1st && till_end_of_2nd) { // Paths coincide completely. + _mergeCoincidingEdges(first, second); + } else if (till_end_of_1st) { + // Paths coincide until the end of the 1st one, which however isn't the end of the + // 2nd one; for example, the first one could be the vertical riser of the letter L + // whereas the second one – the entire letter stroke. + _mergeShorterLonger(vertex, first, second, split.second); + } else if (till_end_of_2nd) { + // The same but with with the second edge shorter than the first one. + _mergeShorterLonger(vertex, second, first, split.first); + } else { // A Y-shaped split. + _mergeWyeConfiguration(vertex, first, second, split); + } + return false; // We've glued so no need to swap anything. +} + +/** + * \brief Analyze a loop path a with self-tangency at the attachment point (a teardrop). + * + * The following steps are taken: + * \li If the loop encloses zero area and \c deloop is true, the loop is detached. + * \li If the two arms of the loop split out immediately, the loop is left alone and we + * only check whether the incidences should be swapped. + * \li If the loop overlaps itself near the vertex, resembling a lasso, we split it into + * a shorter simple path and a smaller loop attached to the end of the shorter path. + * + * \param vertex The vertex at which the teardrop originates. + * \param first The first incidence of the loop to the vertex. + * \param second The second incidence of the loop to the vertex. + * \param deloop Whether the loop should be removed if it doesn't enclose any area + * (i.e., the path exactly backtracks on itself). + * \return Whether the two incidences of the loop to the vertex should be swapped. + */ +template<typename EL> +bool PlanarGraph<EL>::_reglueTeardrop(typename PlanarGraph<EL>::Vertex &vertex, + typename PlanarGraph<EL>::Incidence *first, + typename PlanarGraph<EL>::Incidence *second, bool deloop) +{ + // Calculate the area enclosed by the teardrop. + // The convention is that the unit circle (cos(t), sint(t)), t from 0 to 2pi, + // encloses an area of +pi. + auto &edge = _edges[first->index]; + Path loop = edge.path; loop.close(); + double signed_area = closedPathArea(loop); + + if (deloop && are_near(signed_area, 0.0, _precision)) { + edge.detach(); + _throwAway(&vertex, first); + _throwAway(&vertex, second); + return false; + } + + auto split = parting_point(loop, loop.reversed(), _precision); + if (are_near(split.point(), vertex.point(), _precision)) { + // The loop spreads out immediately. We simply check if the incidences should be swapped. + // We want them to be ordered such that the signed area encircled by the path going out + // at the first incidence and coming back at the second (with this orientation) is > 0. + return (first->sign == Incidence::START) ^ (signed_area > 0.0); + } + + // The loop encloses a nonzero area, but the two branches don't separate at the starting + // point. Instead, they travel together for a while before they split like a lasso. + _reglueLasso(vertex, first, second, split); + return false; +} + +/** + * \brief Reglue a lasso-shaped loop, separating it into the "free rope" and the "hoop". + * + * The lasso is an edge looping back to the same vertex, where the closed path encloses + * a non-zero area, but its two branches don't separate at the starting point. Instead, + * they travel together for a while (forming the doubled-up "free rope") before they + * split like a lasso. This function cuts the lasso at the split point: + * \code{.unparsed} + * ____ ____ + * / \ / \ + * VERTEX =====< | ==> VERTEX ------ NEW + NEW < | + * \____/ (lasso) (rope) \____/ (hoop) + * + * \endcode + * + * \param vertex A reference to the vertex where the lasso is attached. + * \param first The first incidence of the lasso to the vertex. + * \param second The second incidence of the lasso to the vertex. + * \param split The point where the free rope of the lasso ends and the hoop begins. + */ +template<typename EL> +void PlanarGraph<EL>::_reglueLasso(typename PlanarGraph<EL>::Vertex &vertex, + typename PlanarGraph<EL>::Incidence *first, + typename PlanarGraph<EL>::Incidence *second, + PathIntersection const &split) +{ + unsigned lasso = first->index; + + // Create the "free rope" path. + auto rope = _edges[lasso].path.portion(PATH_START, split.first); + rope.setInitial(vertex.point()); + rope.setFinal(split.point()); + double const rope_final_backward_azimuth = _getAzimuth(-rope.finalUnitTangent()); + + // Compute the new label of the rope edge. + auto oriented_as_loop = _edges[lasso].label; + auto reversed = oriented_as_loop; reversed.onReverse(); + oriented_as_loop.onMergeWith(reversed); + + // Insert the rope and its endpoint. + unsigned const rope_index = _edges.size(); + auto &rope_edge = _edges.emplace_back(std::move(rope), std::move(oriented_as_loop)); + auto const new_split_vertex = _ensureVertexAt(split.point()); + + // Reuse lasso's first incidence as the incidence to the rope (azimuth can stay). + first->index = rope_index; + first->sign = Incidence::START; + + // Connect the rope to the newly created split vertex. + new_split_vertex->_addIncidence(rope_index, rope_final_backward_azimuth, Incidence::END); + rope_edge.start = &vertex; + rope_edge.end = new_split_vertex; + + // Insert the hoop + auto hoop = _edges[lasso].path.portion(split.first, + _reversePathTime(split.second, _edges[lasso].path)); + hoop.setInitial(split.point()); + hoop.setFinal(split.point()); + insertEdge(std::move(hoop), EL(_edges[lasso].label)); + + // Detach the original lasso edge and mark the second incidence for cleanup. + _edges[lasso].detach(); + _throwAway(&vertex, second); +} + +/** + * \brief Completely coallesce two fully overlapping edges. + * + * In practice, the first edge stays and the second one gets detached from the graph. + * + * \param first An iterator to the first edge's incidence to a common vertex. + * \param second An iterator to the second edge's incidence to a common vertex. + */ +template<typename EL> +void PlanarGraph<EL>::_mergeCoincidingEdges(typename PlanarGraph<EL>::Incidence *first, + typename PlanarGraph<EL>::Incidence *second) +{ + auto &surviver = _edges[first->index]; + auto &casualty = _edges[second->index]; + + auto other_label = casualty.label; + if (first->sign != second->sign) { // Logically reverse the label before merging. + other_label.onReverse(); + } + surviver.label.onMergeWith(other_label); + + // Mark both incidences of the second edge as junk and detach it. + auto [start_vertex, start_inc] = getIncidence(second->index, Incidence::START); + _throwAway(start_vertex, start_inc); + auto [end_vertex, end_inc] = getIncidence(second->index, Incidence::END); + _throwAway(end_vertex, end_inc); + casualty.detach(); +} + +/** + * \brief Merge a longer edge with a shorter edge that overlaps it. + * + * In practice, the shorter edge remains unchanged and the longer one is trimmed to + * become just the part extending past the shorter one. + * + * \param vertex The vertex where the overlap starts. + * \param shorter The incidence of the shorter edge to the common vertex. + * \param longer The incidence of the longer edge to the common vertex. + * \param time_on_longer The PathTime on the longer edge at which it passes through + * the endpoint of the shorter edge. + */ +template<typename EL> +void PlanarGraph<EL>::_mergeShorterLonger(typename PlanarGraph<EL>::Vertex &vertex, + typename PlanarGraph<EL>::Incidence *shorter, + typename PlanarGraph<EL>::Incidence *longer, + PathTime const &time_on_longer) +{ + auto &shorter_edge = _edges[shorter->index]; + auto &longer_edge = _edges[longer->index]; + + // Get the vertices at the far ends of both edges. + auto shorter_far_end = (shorter->sign == Incidence::START) ? shorter_edge.end + : shorter_edge.start; + /// Whether the longer edge heads out of the vertex. + bool const longer_out = (longer->sign == Incidence::START); + auto longer_far_end = longer_out ? longer_edge.end : longer_edge.start; + + // Copy the longer edge's label and merge with that of the shorter. + auto longer_label = longer_edge.label; + if (shorter->sign != longer->sign) { + longer_label.onReverse(); + } + shorter_edge.label.onMergeWith(longer_label); + + // Create the trimmed path (longer minus shorter). + Path trimmed; + double trimmed_departure_azimuth; + if (longer_out) { + trimmed = longer_edge.path.portion(time_on_longer, _pathEnd(longer_edge.path)); + longer_edge.start = shorter_far_end; + trimmed.setInitial(shorter_far_end->point()); + trimmed.setFinal(longer_far_end->point()); + trimmed_departure_azimuth = _getAzimuth(trimmed.initialUnitTangent()); + } else { + trimmed = longer_edge.path.portion(PATH_START, _reversePathTime(time_on_longer, + longer_edge.path)); + longer_edge.end = shorter_far_end; + trimmed.setInitial(longer_far_end->point()); + trimmed.setFinal(shorter_far_end->point()); + trimmed_departure_azimuth = _getAzimuth(-trimmed.finalUnitTangent()); + } + + // Set the trimmed path as the new path of the longer edge and set up the incidences: + longer_edge.path = std::move(trimmed); + shorter_far_end->_addIncidence(longer->index, trimmed_departure_azimuth, longer->sign); + + // Throw away the old start incidence of the longer edge. + _throwAway(&vertex, longer); +} + +/** + * \brief Merge a pair of partially overlapping edges, producing a Y-split at a new vertex. + * + * This topological modification is performed by inserting a new vertex at the three-way + * point (where the two paths separate) and clipping the original edges to that point. + * In this way, the original edges become the "arms" of the Y-shape. In addition, a new + * edge is inserted, forming the "stem" of the Y. + * + * \param vertex The vertex from which the partially overlapping edges originate (bottom of Y). + * \param first The incidence to the first edge (whose path is the stem and one arm of the Y). + * \param second The incidence to the second edge (stem and the other arm of the Y). + * \param fork The splitting point of the two paths. + */ +template<typename EL> +void PlanarGraph<EL>::_mergeWyeConfiguration(typename PlanarGraph<EL>::Vertex &vertex, + typename PlanarGraph<EL>::Incidence *first, + typename PlanarGraph<EL>::Incidence *second, + PathIntersection const &fork) +{ + bool const first_is_out = (first->sign == Incidence::START); + bool const second_is_out = (second->sign == Incidence::START); + + auto &first_edge = _edges[first->index]; + auto &second_edge = _edges[second->index]; + + // Calculate the path forming the stem of the Y: + auto stem_path = getOutgoingPath(first).portion(PATH_START, fork.first); + stem_path.setInitial(vertex.point()); + stem_path.setFinal(fork.point()); + + /// A closure to clip the path of an original edge to the fork point. + auto const clip_to_fork = [&](PathTime const &t, Edge &e, bool out) { + if (out) { // Trim from time to end + e.path = e.path.portion(t, _pathEnd(e.path)); + e.path.setInitial(fork.point()); + } else { // Trim from reverse-end to reverse-time + e.path = e.path.portion(PATH_START, _reversePathTime(t, e.path)); + e.path.setFinal(fork.point()); + } + }; + + /// A closure to find the departing azimuth of an edge at the fork point. + auto const departing_azimuth = [&](Edge const &e, bool out) -> double { + return _getAzimuth((out) ? e.path.initialUnitTangent() + : -e.path.finalUnitTangent()); + }; + + // Clip the paths obtaining the arms of the Y. + clip_to_fork(fork.first, first_edge, first_is_out); + clip_to_fork(fork.second, second_edge, second_is_out); + + // Create the fork vertex and set up its incidences. + auto const fork_vertex = _ensureVertexAt(fork.point()); + fork_vertex->_addIncidence(first->index, departing_azimuth(first_edge, first_is_out), + first->sign); + fork_vertex->_addIncidence(second->index, departing_azimuth(second_edge, second_is_out), + second->sign); + + // Repoint the ends of the edges that were clipped + (first_is_out ? first_edge.start : first_edge.end) = fork_vertex; + (second_is_out ? second_edge.start : second_edge.end) = fork_vertex; + + /// A closure to get a consistently oriented label of an edge. + auto upwards_oriented_label = [&](Edge const &e, bool out) -> EL { + auto label = e.label; + if (!out) { + label.onReverse(); + } + return label; + }; + + auto stem_label = upwards_oriented_label(first_edge, first_is_out); + stem_label.onMergeWith(upwards_oriented_label(second_edge, second_is_out)); + auto stem_departure_from_fork = _getAzimuth(-stem_path.finalUnitTangent()); + + // Insert the stem of the Y-configuration. + unsigned const stem_index = _edges.size(); + auto &stem_edge = _edges.emplace_back(std::move(stem_path), std::move(stem_label)); + stem_edge.start = &vertex; + stem_edge.end = fork_vertex; + + // Set up the incidences. + fork_vertex->_addIncidence(stem_index, stem_departure_from_fork, Incidence::END); + first->index = stem_index; + first->sign = Incidence::START; + _throwAway(&vertex, second); +} + +template<typename EL> +typename PlanarGraph<EL>::Incidence* +PlanarGraph<EL>::nextIncidence(typename PlanarGraph<EL>::VertexIterator const &vertex, + double azimuth, bool clockwise) const +{ + auto &incidences = vertex._incidences; + Incidence *result = nullptr; + + if (incidences.empty()) { + return result; + } + // Normalize azimuth to the interval [-pi; pi]. + auto angle = Angle(azimuth); + + if (clockwise) { // Go backwards and find a lower bound + auto it = std::find_if(incidences.rbegin(), incidences.rend(), [=](auto inc) -> bool { + return inc.azimuth <= angle; + }); + if (it == incidences.rend()) { + // azimuth is lower than the azimuths of all incidences; + // going clockwise we wrap back to the highest azimuth (last incidence). + return &incidences.back(); + } + result = &(*it); + } else { + auto it = std::find_if(incidences.begin(), incidences.end, [=](auto inc) -> bool { + return inc.azimuth >= angle; + }); + if (it == incidences.end()) { + // azimuth is higher than the azimuths of all incidences; + // going counterclockwise we wrap back to the lowest azimuth. + return &incidences.front(); + } + result = &(*it); + } + return result; +} + +/** Return the signed area enclosed by a closed path. */ +template<typename EL> +double PlanarGraph<EL>::closedPathArea(Path const &path) +{ + double area; + Point _; + centroid(path.toPwSb(), _, area); + return -area; // Our convention is that the Y-axis points up +} + +/** \brief Determine whether the first path deviates to the left of the second. + * + * The two paths are assumed to have identical or nearly identical starting points + * but not an overlapping initial portion. The concept of "left" is based on the + * y-axis pointing up. + * + * \param first The first path. + * \param second The second path. + * + * \return True if the first path deviates towards the left of the second; + * False if the first path deviates towards the right of the second. + */ +template<typename EL> +bool PlanarGraph<EL>::deviatesLeft(Path const &first, Path const &second) +{ + auto start = middle_point(first.initialPoint(), second.initialPoint()); + auto tangent_between = middle_point(first.initialUnitTangent(), second.initialUnitTangent()); + if (tangent_between.isZero()) { + return false; + } + auto tangent_line = Line::from_origin_and_vector(start, tangent_between); + + // Find the first non-degenerate curves on both paths + std::unique_ptr<Curve> c[2]; + auto const find_first_nondegen = [](std::unique_ptr<Curve> &pointer, Path const &path) { + for (auto const &c : path) { + if (!c.isDegenerate()) { + pointer.reset(c.duplicate()); + return; + } + } + }; + + find_first_nondegen(c[0], first); + find_first_nondegen(c[1], second); + if (!c[0] || !c[1]) { + return false; + } + + // Find the bounding boxes + Rect const bounding_boxes[] { + c[0]->boundsExact(), + c[1]->boundsExact() + }; + + // For a bounding box, find the corner that goes the furthest in the direction of the + // tangent vector. + auto const furthest_corner = [&](Rect const &r) -> unsigned { + Coord max_dot = dot(r.corner(0) - start, tangent_between); + unsigned result = 0; + for (unsigned i = 1; i < 4; i++) { + auto current_dot = dot(r.corner(i), tangent_between); + if (current_dot > max_dot) { + max_dot = current_dot; + result = i; + } else if (current_dot == max_dot) { + // Disambiguate based on proximity to the tangent line. + auto const offset = start + tangent_between; + if (distance(offset, r.corner(i)) < distance(offset, r.corner(result))) { + result = i; + } + } + } + return result; + }; + + // Calculate the corner points overlooking the "rift" between the paths. + Point corner_points[2]; + for (size_t i : {0, 1}) { + corner_points[i] = bounding_boxes[i].corner(furthest_corner(bounding_boxes[i])); + } + + // Find a vantage point from which we can best observe the splitting paths. + Point vantage_point; + bool found = false; + if (corner_points[0] != corner_points[1]) { + auto line_connecting_corners = Line(corner_points[0], corner_points[1]); + auto xing = line_connecting_corners.intersect(tangent_line); + if (!xing.empty()) { + vantage_point = xing[0].point(); + found = true; + } + } + if (!found) { + vantage_point = tangent_line.pointAt(tangent_line.timeAtProjection(corner_points[0])); + } + + // Move to twice as far in the direction of the vantage point. + vantage_point += vantage_point - start; + + // Find the points on both curves that are nearest to the vantage point. + Coord nearest[2]; + for (size_t i : {0, 1}) { + nearest[i] = c[i]->nearestTime(vantage_point); + } + + // Clip to the nearest points and examine the closed contour. + Path closed_contour(start); + closed_contour.setStitching(true); + closed_contour.append(c[0]->portion(0, nearest[0])); + closed_contour = closed_contour.reversed(); + closed_contour.setStitching(true); + closed_contour.append(c[1]->portion(0, nearest[1])); + closed_contour.close(); + return !path_direction(closed_contour); // Reverse to match the convention that y-axis is up. +} + +} // namespace Geom + +#endif // LIB2GEOM_SEEN_PLANAR_GRAPH_H + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : |