planar-graph.h

Go to the documentation of this file.

 
/*
 * 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 {
 
template<typename EdgeLabel>
concept EdgeLabelConcept = requires(EdgeLabel el, EdgeLabel const &other) {
    el.onReverse();
    el.onMergeWith(other);
    el.onDetach();
    el = other;
};
 
template<typename EdgeLabel>
requires EdgeLabelConcept<EdgeLabel>
class PlanarGraph
{
public:
 
    struct Incidence
    {
        using Sign = bool;
        inline static Sign const START = false;
        inline static Sign const END   = true;
 
        double azimuth; 
        unsigned index; 
        Sign sign;      
        bool invalid = false; 
 
        Incidence(unsigned edge_index, double departure_azimuth, Sign which_end)
            : azimuth{departure_azimuth}
            , index{edge_index}
            , sign{which_end}
        {
        }
        ~Incidence() = default;
 
        inline bool operator<(Incidence const &other) const { return azimuth < other.azimuth; }
 
        inline bool operator<(double angle) const { return azimuth < angle; }
 
        inline bool operator==(Incidence const &other) const
        {
            return index == other.index && sign == other.sign;
        }
    };
    using IncIt = std::list<Incidence>::iterator;
    using IncConstIt = std::list<Incidence>::const_iterator;
 
    class Vertex
    {
    private:
        Point const _position;            
        std::list<Incidence> _incidences; 
        unsigned mutable _flags = 0;      
 
        inline static Point::LexLess<X> const _cmp; 
 
    public:
        Vertex(Point const &pos)
            : _position{pos}
        {
        }
 
        Point const &point() const { return _position; }
 
        auto const &getIncidences() const { return _incidences; }
 
        bool operator<(Vertex const &other) const { return _cmp(_position, other._position); }
 
        unsigned flags() const { return _flags; } 
        void setFlags(unsigned flags) const { _flags = flags; } 
 
        IncConstIt cyclicNextIncidence(IncConstIt it) const { return cyclic_next(it, _incidences); }
 
        IncConstIt cyclicPrevIncidence(IncConstIt it) const { return cyclic_prior(it, _incidences); }
 
        Incidence *cyclicNextIncidence(Incidence *from)
        {
            return &(*cyclicNextIncidence(_incidencePtr2It(from)));
        }
        Incidence *cyclicPrevIncidence(Incidence *from)
        {
            return &(*cyclicPrevIncidence(_incidencePtr2It(from)));
        }
 
    private:
        IncIt cyclicNextIncidence(IncIt it) { return cyclic_next(it, _incidences); }
        IncIt cyclicPrevIncidence(IncIt it) { return cyclic_prior(it, _incidences); }
 
        Incidence &_addIncidence(unsigned edge_index, double azimuth, 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));
        }
 
        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 = std::list<Vertex>::iterator;
 
    struct Edge
    {
        Vertex *start = nullptr, *end = nullptr; 
        Path path; 
        bool detached = false; 
        bool inserted_as_detached = false; 
        EdgeLabel mutable label; 
 
        Edge(Path &&movein_path, EdgeLabel &&movein_label)
            : path{movein_path}
            , label{movein_label}
        {
        }
 
        void detach()
        {
            detached = true;
            label.onDetach();
        }
    };
    using EdgeIterator = std::vector<Edge>::iterator;
    using EdgeConstIterator = std::vector<Edge>::const_iterator;
 
private:
    double const _precision; 
    std::list<Vertex> _vertices; 
    std::vector<Edge> _edges; 
    std::vector< std::pair<Vertex *, Incidence *> > _junk; 
    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; });
    }
 
    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);
 
    void reserveEdgeCapacity(size_t capacity) { _edges.reserve(capacity); }
 
    unsigned insertEdge(Path &&path, EdgeLabel &&edge = EdgeLabel());
    unsigned insertDetached(Path &&path, EdgeLabel &&edge = EdgeLabel());
 
    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));
    }
 
    std::pair<Vertex *, Incidence *>
    getIncidence(unsigned edge_index, 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)};
    }
 
    inline Incidence const &nextIncidence(VertexIterator const &vertex, IncConstIt const &incidence,
                                          bool clockwise = false) const
    {
        return clockwise ? *(vertex->_cyclicPrevIncidence(incidence))
                         : *(vertex->_cyclicNextIncidence(incidence));
    }
 
    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));
    }
 
    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);
    }
 
    Incidence *nextIncidence(VertexIterator const &vertex, double azimuth,
                             bool clockwise = false) const;
 
    Path getOutgoingPath(Incidence const *incidence) const
    {
        return incidence ? _getPathImpl(incidence, Incidence::START) : Path();
    }
 
    Path getIncomingPath(Incidence const *incidence) const
    {
        return incidence ? _getPathImpl(incidence, Incidence::END) : Path();
    }
 
private:
    inline Path _getPathImpl(Incidence const *incidence, Incidence::Sign origin) const
    {
        return (incidence->sign == origin) ? _edges[incidence->index].path
                                           : _edges[incidence->index].path.reversed();
    }
 
    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 ===
 
    inline static double _getAzimuth(Point const &vec) { return vec.isZero() ? 0.0 : atan2(vec); }
 
    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);
    }
 
    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);
};
 
template<typename EL>
requires EdgeLabelConcept<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)));
}
 
template<typename EdgeLabel>
requires EdgeLabelConcept<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;
}
 
template<typename EdgeLabel>
requires EdgeLabelConcept<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;
}
 
template<typename EdgeLabel>
requires EdgeLabelConcept<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();
}
 
template<typename EdgeLabel>
requires EdgeLabelConcept<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;
}
 
template<typename EdgeLabel>
requires EdgeLabelConcept<EdgeLabel>
void PlanarGraph<EdgeLabel>::_regularizeVertex(Vertex &vertex,
                                               double angle_precision, bool deloop)
{
    auto &incidences = vertex._incidences;
 
    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.
 
    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; 
        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);
    }
}
 
template<typename EL>
requires EdgeLabelConcept<EL>
void PlanarGraph<EL>::_reglueTangentFan(Vertex &vertex,
                                        IncIt const &first,
                                        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;
            }
        }
    }
}
 
template<typename EL>
requires EdgeLabelConcept<EL>
bool PlanarGraph<EL>::_compareAndReglue(Vertex &vertex,
                                        Incidence *first,
                                        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.
}
 
template<typename EL>
requires EdgeLabelConcept<EL>
bool PlanarGraph<EL>::_reglueTeardrop(Vertex &vertex,
                                      Incidence *first,
                                      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;
}
 
template<typename EL>
requires EdgeLabelConcept<EL>
void PlanarGraph<EL>::_reglueLasso(Vertex &vertex,
                                   Incidence *first,
                                   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);
}
 
template<typename EL>
requires EdgeLabelConcept<EL>
void PlanarGraph<EL>::_mergeCoincidingEdges(Incidence *first,
                                            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();
}
 
template<typename EL>
requires EdgeLabelConcept<EL>
void PlanarGraph<EL>::_mergeShorterLonger(Vertex &vertex,
                                          Incidence *shorter,
                                          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;
    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);
}
 
template<typename EL>
requires EdgeLabelConcept<EL>
void PlanarGraph<EL>::_mergeWyeConfiguration(Vertex &vertex,
                                             Incidence *first,
                                             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());
 
    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());
        }
    };
 
    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;
 
    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>
requires EdgeLabelConcept<EL>
PlanarGraph<EL>::Incidence*
PlanarGraph<EL>::nextIncidence(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;
}
 
template<typename EL>
requires EdgeLabelConcept<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
}
 
template<typename EL>
requires EdgeLabelConcept<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 :