From 90fccd526c74b3a3fe0615cc996cb2d1c068c233 Mon Sep 17 00:00:00 2001 From: Tinko Bartels Date: Tue, 13 Aug 2019 00:22:43 +0200 Subject: [PATCH] Add triangulation model and delaunay triangulation algorithm. --- extensions/example/Jamfile | 1 + extensions/example/triangulation/Jamfile.v2 | 11 + .../triangulation/triangulation_example.cpp | 56 ++ extensions/test/triangulation/Jamfile | 1 + .../test/triangulation/triangulation.cpp | 45 + .../algorithms/delaunay_triangulation.hpp | 67 ++ .../geometries/triangulation.hpp | 814 ++++++++++++++++++ .../geometries/voronoi_adaptor.hpp | 132 +++ .../cartesian/accelerated_shull.hpp | 63 ++ .../cartesian/detail/accelerated_shull.hpp | 448 ++++++++++ .../strategies/delaunay_triangulation.hpp | 41 + .../triangulation/triangulation.hpp | 18 + 12 files changed, 1697 insertions(+) create mode 100644 extensions/example/triangulation/Jamfile.v2 create mode 100644 extensions/example/triangulation/triangulation_example.cpp create mode 100644 extensions/test/triangulation/triangulation.cpp create mode 100644 include/boost/geometry/extensions/triangulation/algorithms/delaunay_triangulation.hpp create mode 100644 include/boost/geometry/extensions/triangulation/geometries/triangulation.hpp create mode 100644 include/boost/geometry/extensions/triangulation/geometries/voronoi_adaptor.hpp create mode 100644 include/boost/geometry/extensions/triangulation/strategies/cartesian/accelerated_shull.hpp create mode 100644 include/boost/geometry/extensions/triangulation/strategies/cartesian/detail/accelerated_shull.hpp create mode 100644 include/boost/geometry/extensions/triangulation/strategies/delaunay_triangulation.hpp create mode 100644 include/boost/geometry/extensions/triangulation/triangulation.hpp diff --git a/extensions/example/Jamfile b/extensions/example/Jamfile index 8110621689..1ae69cce26 100644 --- a/extensions/example/Jamfile +++ b/extensions/example/Jamfile @@ -13,3 +13,4 @@ project boost-geometry-examples-extensions ; build-project gis ; +build-project triangulation ; diff --git a/extensions/example/triangulation/Jamfile.v2 b/extensions/example/triangulation/Jamfile.v2 new file mode 100644 index 0000000000..80fcf24de2 --- /dev/null +++ b/extensions/example/triangulation/Jamfile.v2 @@ -0,0 +1,11 @@ +# Boost.Geometry (aka GGL, Generic Geometry Library) +# +# Use, modification and distribution is subject to the Boost Software License, +# Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +# http://www.boost.org/LICENSE_1_0.txt) + +project boost-geometry-example-extensions-triangulation + : # requirements + ; + +exe triangulation_example : triangulation_example.cpp ; diff --git a/extensions/example/triangulation/triangulation_example.cpp b/extensions/example/triangulation/triangulation_example.cpp new file mode 100644 index 0000000000..9b984fb18e --- /dev/null +++ b/extensions/example/triangulation/triangulation_example.cpp @@ -0,0 +1,56 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2019 Tinko Bartels, Berlin, Germany. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +// Triangulation Example - computes a delaunay triangulation and draws the +// triangulation and the corresponding finite voronoi edges + +#include +#include +#include + +#include +#include +#include + +const int samples = 100; + +int main() +{ + namespace bg = boost::geometry; + typedef bg::model::point point; + typedef bg::model::triangulation triangulation; + + std::default_random_engine gen(1); + std::uniform_real_distribution<> dist(0.0, 1.0); + std::vector in; + for (int i = 0 ; i < samples ; ++i) + { + in.push_back(point(dist(gen), dist(gen))); + } + triangulation out(samples); + bg::delaunay_triangulation(in, out); + + std::ofstream svg("triangulation.svg"); + bg::svg_mapper mapper(svg, 720, 720); + mapper.add(point(0, 0)); + mapper.add(point(1, 1)); + + for (auto const& f : bg::face_range(out)) + { + mapper.map(f, "fill-opacity:0.3;fill:rgb(102,102,201);stroke:rgb(51,51,152);"); + } + + typedef bg::model::voronoi_face_view voronoi_face_view; + for (auto it = out.vertices_begin() ; it != out.vertices_end() ; ++it) + { + voronoi_face_view vfv(out, it); + mapper.map(voronoi_face_view(out, it), "opacity:1.0;fill:none;stroke:rgb(255,0,0);stroke-width:1"); + } + + return 0; +} diff --git a/extensions/test/triangulation/Jamfile b/extensions/test/triangulation/Jamfile index ad0b333c80..1f141984d7 100644 --- a/extensions/test/triangulation/Jamfile +++ b/extensions/test/triangulation/Jamfile @@ -8,5 +8,6 @@ test-suite boost-geometry-extensions-triangulation : [ run side_robust.cpp ] [ run in_circle_robust.cpp ] + [ run triangulation.cpp ] ; diff --git a/extensions/test/triangulation/triangulation.cpp b/extensions/test/triangulation/triangulation.cpp new file mode 100644 index 0000000000..6bafb48ade --- /dev/null +++ b/extensions/test/triangulation/triangulation.cpp @@ -0,0 +1,45 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// Unit Test + +// Copyright (c) 2019 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2019 program. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#include + +#include + +#include +#include + +namespace bg = boost::geometry; + +template +void test_all() +{ + typedef bg::model::triangulation triangulation; + std::vector

in; + in.push_back(P(-1, -2)); + in.push_back(P(1, -2)); + in.push_back(P(-2, 0)); + in.push_back(P(2, 0)); + in.push_back(P(1, 2)); + in.push_back(P(-1, 2)); + in.push_back(P(0, 0)); + in.push_back(P(0, 1)); + triangulation t(8); + bg::delaunay_triangulation(in, t); + BOOST_CHECK( t.valid() ); +} + + +int test_main(int, char* []) +{ + test_all >(); + return 0; +} diff --git a/include/boost/geometry/extensions/triangulation/algorithms/delaunay_triangulation.hpp b/include/boost/geometry/extensions/triangulation/algorithms/delaunay_triangulation.hpp new file mode 100644 index 0000000000..d7bc8988da --- /dev/null +++ b/include/boost/geometry/extensions/triangulation/algorithms/delaunay_triangulation.hpp @@ -0,0 +1,67 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2019 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2019 program. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_ALGORITHMS_DELAUNAY_TRIANGULATION_HPP +#define BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_ALGORITHMS_DELAUNAY_TRIANGULATION_HPP + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + + +#include +#include +#include +#include +#include +#include +#include + +namespace boost { namespace geometry +{ + +template +< + typename PointContainer, + typename Triangulation, + typename SideStrategy = strategy::side::side_by_triangle<>, + typename InCircleStrategy = strategy::in_circle::in_circle_robust<>, + typename ConstructionStrategy = + typename strategy::delaunay_triangulation::services::default_strategy + < + typename tag::type, + typename cs_tag::type, + dimension::value + >::type +> +inline void delaunay_triangulation(PointContainer const & in, + Triangulation& out, + bool legalize = true) +{ + ConstructionStrategy::template apply + < + PointContainer, + Triangulation, + SideStrategy, + InCircleStrategy + >(in, out, legalize); +} + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_ALGORITHMS_DELAUNAY_TRIANGULATION_HPP diff --git a/include/boost/geometry/extensions/triangulation/geometries/triangulation.hpp b/include/boost/geometry/extensions/triangulation/geometries/triangulation.hpp new file mode 100644 index 0000000000..509a9eba37 --- /dev/null +++ b/include/boost/geometry/extensions/triangulation/geometries/triangulation.hpp @@ -0,0 +1,814 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2019 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2019 program. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_GEOMETRIES_TRIANGULATION_HPP +#define BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_GEOMETRIES_TRIANGULATION_HPP + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace boost { namespace geometry +{ +template +struct face_range_type {}; + +template +typename face_range_type::type + face_range(Triangulation const& t); + +template +struct triangulation_face {}; + +namespace model +{ + +template +< + typename Point, + bool ClockWise = true, + template class VertexContainer = std::vector, + template class FaceContainer = std::vector, + template class VertexAllocator = std::allocator, + template class FaceAllocator = std::allocator +> +class triangulation; + +template +struct vertex_ref +{ + typename Triangulation::point_type m_p; + typename Triangulation::face_iterator m_f; +}; + +template +struct face_ref +{ + typedef typename Triangulation::point_type point_type; + typedef typename Triangulation::vertex_iterator vertex_iterator; + typedef typename Triangulation::face_iterator face_iterator; + typedef std::array vertex_container; + vertex_container m_v; + std::array m_f; + std::array m_o; +public: + typedef boost::indirect_iterator + const_iterator; + typedef boost::indirect_iterator + iterator; + + const_iterator begin() const { return const_iterator(m_v.begin()); } + const_iterator end() const { return const_iterator(m_v.end()); } +}; + +template +< + typename Value, + template class Container, + template class Allocator +> +struct reserve_if_vector{ + static void apply(Container>& c, std::size_t n) {} +}; + +template class Allocator> +struct reserve_if_vector{ + static void apply(std::vector>& c, + std::size_t n) { c.reserve(n);} +}; + +template +< + typename Point, + bool ClockWise, + template class VertexContainer, + template class FaceContainer, + template class VertexAllocator, + template class FaceAllocator +> +class triangulation +{ +private: + BOOST_CONCEPT_ASSERT( (concepts::Point) ); +public: + typedef Point point_type; + typedef face_ref face_type; + typedef vertex_ref vertex_type; + typedef unsigned short face_vertex_index; + typedef VertexContainer> + vertex_container; + typedef FaceContainer> face_container; + typedef typename face_container::iterator face_iterator; + typedef typename face_container::const_iterator const_face_iterator; + typedef typename vertex_container::iterator vertex_iterator; + typedef typename vertex_container::const_iterator const_vertex_iterator; + typedef typename coordinate_type::type coordinate_type; + typedef typename model::segment segment_type; + + static const order_selector point_order = + ClockWise ? clockwise : counterclockwise; + struct halfedge_index + { + halfedge_index(face_iterator f, face_vertex_index v):m_f(f), m_v(v) {} + face_iterator m_f; + face_vertex_index m_v; + }; + + struct fulledge_index + { + fulledge_index(face_iterator f1, + face_vertex_index v1, + face_iterator f2, + face_vertex_index v2) + : m_f1(f1), + m_f2(f2), + m_v1(v1), + m_v2(v2) {} + fulledge_index(halfedge_index e) + : m_f1(e.m_f), + m_f2(e.m_f->m_f[ e.m_v ]), + m_v1(e.m_v), + m_v2(e.m_f->m_o[ e.m_v ]) {} + face_iterator m_f1, m_f2; + face_vertex_index m_v1, m_v2; + }; + + static face_iterator invalid() { return face_iterator(); } + triangulation(std::size_t points = 3) + { + reserve_if_vector:: + apply(m_vertices, points); + reserve_if_vector:: + apply(m_faces, 2 * points - 5); + } + + triangulation(std::size_t points, std::size_t faces) + { + reserve_if_vector:: + apply(m_vertices, points); + reserve_if_vector:: + apply(m_faces, faces); + } + + template + triangulation(InputIt begin, InputIt end) + { + m_vertices.assign(begin, end); + reserve_if_vector + ::apply(m_faces, 2 * m_vertices.size() - 5); + } + + typename vertex_container::iterator vertices_begin() + { + return m_vertices.begin(); + } + + typename vertex_container::const_iterator vertices_begin() const + { + return m_vertices.cbegin(); + } + + typename vertex_container::iterator vertices_end() + { + return m_vertices.end(); + } + + typename vertex_container::const_iterator vertices_end() const + { + return m_vertices.cend(); + } + + vertex_iterator add_vertex(const Point& p) + { + return m_vertices.insert(m_vertices.end(), vertex_type{p, invalid()}); + } + + face_container const& face_range() const + { + return m_faces; + } + + vertex_container const& vertex_range() const + { + return m_vertices; + } + + typename face_container::const_iterator faces_cbegin() const + { + return m_faces.cbegin(); + } + + typename face_container::iterator faces_begin() + { + return m_faces.begin(); + } + + typename face_container::const_iterator faces_cend() const + { + return m_faces.cend(); + } + + typename face_container::iterator faces_end() + { + return m_faces.end(); + } + + template + void assign_vertices(InputIt begin, InputIt end) + { + m_vertices.assign(begin, end); + m_faces.reserve(2 * m_vertices.size() - 5); + } + + static Point const& face_vertex(face_iterator f, face_vertex_index v) + { + return f -> m_v[ v ] -> m_p; + } + + static segment_type face_segment(halfedge_index e) + { + return segment_type( + face_vertex(e.m_f, (e.m_v == 2 ? 0 : e.m_v + 1)), + face_vertex(e.m_f, (e.m_v == 0 ? 2 : e.m_v - 1)) ); + } + + static Point& vertex(vertex_iterator v) + { + return v -> m_p; + } + + static face_iterator neighbour(face_iterator f, unsigned short v) + { + return f -> m_f[ v ]; + } + + static const_face_iterator neighbour(const_face_iterator f, + unsigned short v) + { + return f -> m_f[ v ]; + } + + static face_vertex_index opposite(face_iterator f, unsigned short v) + { + return f -> m_o[ v ]; + } + + static halfedge_index opposite(halfedge_index const& e) + { + return halfedge_index{ e.m_f -> m_f[ e.m_v ], e.m_f -> m_o[ e.m_v ] }; + } + + static halfedge_index next(halfedge_index const& e) + { + return halfedge_index{ + e.m_f, + static_cast(e.m_v == 2 ? 0 : e.m_v + 1) }; + } + + static halfedge_index prev(halfedge_index const& e) + { + return halfedge_index{ + e.m_f, + static_cast(e.m_v == 0 ? 2 : e.m_v - 1)}; + } + + vertex_iterator boundary_vertex() const + { + return m_boundary_vertex; + } + + static vertex_iterator boundary_next(vertex_iterator v) + { + face_iterator fi = v -> m_f; + if(fi -> m_v[ 0 ] == v ) + return fi -> m_v[ 1 ]; + else if(fi -> m_v[ 1 ] == v) + return fi -> m_v[ 2 ]; + else + return fi -> m_v[ 0 ]; + } + + vertex_iterator boundary_prev(vertex_iterator v) const + { + face_iterator fi = v -> m_f; + unsigned short vi; + if(fi -> m_v[ 0 ] == v) vi = 0; + else if(fi -> m_v[ 1 ] == v) vi = 1; + else vi = 2; + if(m_faces.size() == 1) + return vi == 0 ? fi -> m_v[ 2 ] : fi -> m_v[ vi - 1 ]; + halfedge_index e = next(halfedge_index{fi, vi}); + while( opposite(e).m_f != invalid() ) + { + e = prev(opposite(e)); + } + return e.m_f -> m_v[ e.m_v == 2 ? 0 : e.m_v + 1 ]; + } + + void clear() + { + m_vertices.clear(); + m_faces.clear(); + } + + std::size_t vertices() const + { + return m_vertices.size(); + } + + std::size_t faces() const + { + return m_faces.size(); + } + + static void flip(const halfedge_index& e) + { + face_iterator fi1 = e.m_f; + face_type& f1 = *fi1; + face_iterator fi2 = f1.m_f[ e.m_v ]; + face_type& f2 = *fi2; + unsigned short const& v1 = e.m_v; + unsigned short const v2 = f1.m_o[ v1 ]; + + if( f1.m_v[ v1 == 0 ? 2 : v1 - 1 ] -> m_f == fi1 ) + f1.m_v[ v1 == 0 ? 2 : v1 - 1 ] -> m_f = fi2; + if( f2.m_v[ v2 == 0 ? 2 : v2 - 1 ]-> m_f == fi2) + f2.m_v[ v2 == 0 ? 2 : v2 - 1 ]-> m_f = fi1; + f1.m_v[ v1 == 0 ? 2 : v1 - 1 ] = f2.m_v[ v2 ]; + f2.m_v[ v2 == 0 ? 2 : v2 - 1 ] = f1.m_v[ v1 ]; + + f1.m_f[v1] = f2.m_f[ v2 == 2 ? 0 : v2 + 1 ]; + f1.m_o[v1] = f2.m_o[ v2 == 2 ? 0 : v2 + 1 ]; + if(f1.m_f[v1] != invalid()) + { + f1.m_f[v1] -> m_f[ f2.m_o[ v2 == 2 ? 0 : v2 + 1 ] ] = fi1; + f1.m_f[v1] -> m_o[ f2.m_o[ v2 == 2 ? 0 : v2 + 1 ] ] = v1; + } + f2.m_f[v2] = f1.m_f[ v1 == 2 ? 0 : v1 + 1 ]; + f2.m_o[v2] = f1.m_o[ v1 == 2 ? 0 : v1 + 1 ]; + if(f2.m_f[v2] != invalid()) + { + f2.m_f[v2] -> m_f[ f1.m_o[ v1 == 2 ? 0 : v1 + 1 ] ] = fi2; + f2.m_f[v2] -> m_o[ f1.m_o[ v1 == 2 ? 0 : v1 + 1 ] ] = v2; + } + f1.m_f[ v1 == 2 ? 0 : v1 + 1 ] = fi2; + f1.m_o[ v1 == 2 ? 0 : v1 + 1 ] = v2 == 2 ? 0 : v2 + 1; + f2.m_f[ v2 == 2 ? 0 : v2 + 1 ] = fi1; + f2.m_o[ v2 == 2 ? 0 : v2 + 1 ] = v1 == 2 ? 0 : v1 + 1; + } + + face_iterator add_face_on_boundary(halfedge_index e, vertex_iterator v) + { + const face_iterator f = e.m_f; + const face_vertex_index adj = e.m_v; + f -> m_o[adj] = 0; + m_boundary_vertex = v; + face_iterator pos = m_faces.insert(m_faces.end(), + face_type{ {{v, f->m_v[ adj == 0 ? 2 : adj - 1 ], + f->m_v[ adj == 2 ? 0 : adj + 1 ] }}, + {{f, invalid(), invalid()}}, + {{adj, 4, 4}}}); + f -> m_f[adj] = pos; + v -> m_f = pos; + m_faces.back().m_v[ 2 ]->m_f = pos; + return pos; + } + + fulledge_index next_around_vertex(fulledge_index e) + { + if(e.m_f2 == invalid()) + { + face_vertex_index left_vi = (e.m_v1 == 0 ? 2 : e.m_v1 - 1); + vertex_iterator left_v = e.m_f1->m_v[ left_vi ]; + face_iterator next_f = left_v->m_f; + face_vertex_index next_vi = + (next_f->m_v[ 0 ] == left_v) ? 2 : + ((next_f->m_v[ 1 ] == left_v) ? 0 : 1); + return fulledge_index(invalid(), 4, next_f, next_vi); + } else { + return fulledge_index(halfedge_index(e.m_f2, + e.m_v2 == 0 ? 2 : e.m_v2 - 1)); + } + } + + static fulledge_index begin_vertex_edge(vertex_iterator vi) + { + face_iterator first_f = vi -> m_f; + face_vertex_index next_vi = + (first_f->m_v[ 0 ] == vi) ? 2 : + ((first_f->m_v[ 1 ] == vi) ? 0 : 1); + return fulledge_index(invalid(), 4, first_f, next_vi); + } + + face_iterator add_isolated_face(vertex_iterator v1, + vertex_iterator v2, + vertex_iterator v3) + { + m_boundary_vertex = v1; + face_iterator pos = m_faces.insert( m_faces.end(), + face_type{ + {{ v1, v2, v3 }}, + {{ invalid(), invalid(), invalid() }}, + {{4, 4, 4}} } ); + v1 -> m_f = v2 -> m_f = v3 -> m_f = pos; + return pos; + } + + static halfedge_index face_edge(face_iterator f, face_vertex_index v = 0) + { + return halfedge_index{f, v}; + } + + static void connect(halfedge_index e1, halfedge_index e2) + { + face_iterator f1 = e1.m_f; + face_iterator f2 = e2.m_f; + unsigned short& v1 = e1.m_v; + unsigned short& v2 = e2.m_v; + f1 -> m_f[ v1 ] = f2; + f1 -> m_o[ v1 ] = v2; + f2 -> m_f[ v2 ] = f1; + f2 -> m_o[ v2 ] = v1; + if(f1 -> m_v[ (v1 == 2 ? 0 : v1 + 1) ] -> m_f == f1) + { + f1 -> m_v[ (v1 == 2 ? 0 : v1 + 1) ] -> m_f = f2; + } + if(f2 -> m_v[ (v2 == 2 ? 0 : v2 + 1) ] -> m_f == f2) + { + f2 -> m_v[ (v2 == 2 ? 0 : v2 + 1) ] -> m_f = f1; + } + } + + bool valid() const + { + bool valid = true; + for(const_face_iterator fi = m_faces.begin(); + fi != m_faces.end(); + ++fi) + { + face_type const& f = *fi; + for(unsigned short v = 0 ; v < 3 ; ++v) + { + if(f.m_f[ v ] == invalid()) + continue; + face_vertex_index const& o = f.m_o[ v ]; + valid = valid && (f.m_f[ v ]->m_o[ o ] == v); + if(!valid) return false; + valid = valid && (f.m_f[ v ]->m_f[ o ] == fi); + if(!valid) return false; + valid = valid + && (f.m_v[ (v + 1) % 3 ] == f.m_f[ v ]->m_v[ (o + 2) % 3 ]) + && (f.m_v[ (v + 2) % 3 ] == f.m_f[ v ]->m_v[ (o + 1) % 3 ]); + if(!valid) return false; + if(f.m_o[v] == 4) + { + unsigned short next = (v + 1) % 3; + valid = valid && f.m_v[ next ]->m_f == fi; + if(!valid) return false; + } + } + } + for(const_vertex_iterator vi = m_vertices.cbegin(); + vi != m_vertices.cend(); + ++vi) + { + vertex_type const& v = *vi; + if(v.m_f == invalid()) continue; + + bool found = false; + for(face_vertex_index vj = 0; vj < 3 ; ++vj) + { + found = found || (v.m_f -> m_v[ vj ] == vi); + } + valid = valid && found; + if(!valid) return false; + } + return valid; + } +private: + vertex_container m_vertices; + face_container m_faces; + vertex_iterator m_boundary_vertex; +}; + +template< typename Point > +struct edge_ref +{ + typename triangulation::halfedge_index m_e; + triangulation& m_t; +}; + +template< typename Point > +using triangulation_face_range = typename triangulation::face_container; + +template< typename Point > +using triangulation_vertex_range = + typename triangulation::vertex_container; + +} // namespace model + +#ifndef DOXYGEN_NO_TRAITS_SPECIALIZATIONS +namespace traits +{ + +template +struct tag< model::vertex_ref > +{ typedef point_tag type; }; + +template +struct coordinate_type< model::vertex_ref > +{ + typedef typename coordinate_type + ::type type; +}; + +template +struct dimension< model::vertex_ref > : boost::mpl::int_<2> {}; + +template +struct coordinate_system< model::vertex_ref > +{ + typedef typename coordinate_system + ::type type; +}; + +template +struct access, Dimension> +{ + static typename coordinate_type::type + get(model::vertex_ref const& p) + { + return boost::geometry::get(p.m_p); + } +}; + +template struct tag< model::edge_ref > +{ typedef segment_tag type; }; + +template struct point_type< model::edge_ref > +{ typedef Point type; }; + +template +struct indexed_access, 0, Dimension> +{ + static typename coordinate_type::type + get(model::edge_ref const& p) + { + return get( + p.m_t.face_vertex(p.m_e.m_f, (p.m_e.m_v == 2 ? 0 : p.m_e.m_v + 1)) + ); + } +}; + +template +struct indexed_access, 1, Dimension> +{ + static typename coordinate_type::type + get(model::edge_ref const& p) + { + return get( + p.m_t.face_vertex(p.m_e.m_f, (p.m_e.m_v == 0 ? 2 : p.m_e.m_v - 1))); + } +}; + +template +struct tag> +{ typedef ring_tag type; }; + +template +struct point_order> +{ static const order_selector value = Triangulation::point_order; }; + +template +struct closure> +{ static const closure_selector value = open; }; + +} // namespace traits +#endif // DOXYGEN_NO_TRAITS_SPECIALIZATIONS + +template +struct face_range_type> { + typedef typename model::triangulation_face_range type; +}; + +template +< + typename Point, + bool ClockWise, + template class VertexContainer, + template class FaceContainer, + template class VertexAllocator, + template class FaceAllocator +> +inline typename model::triangulation +< + Point, + ClockWise, + VertexContainer, + FaceContainer, + VertexAllocator, + FaceAllocator +>::vertex_container const& vertex_range( + model::triangulation const& t) +{ return t.vertex_range(); } + +template +< + typename Point, + bool ClockWise, + template class VertexContainer, + template class FaceContainer, + template class VertexAllocator, + template class FaceAllocator +> +inline typename model::triangulation +< + Point, + ClockWise, + VertexContainer, + FaceContainer, + VertexAllocator, + FaceAllocator +>::face_container const& face_range( + model::triangulation const& t) +{ return t.face_range(); } + +template +< + typename Point, + bool ClockWise, + template class VertexContainer, + template class FaceContainer, + template class VertexAllocator, + template class FaceAllocator, + typename FaceIterator, + typename OutputIterator +> +inline void face_adjacent_range( + model::triangulation& t, + FaceIterator fi, + OutputIterator out) +{ + typedef typename model::triangulation triangulation; + typedef typename triangulation::face_iterator face_iterator; + face_iterator const invalid = t.invalid(); + if(fi -> m_f[ 0 ] != invalid ) *out++ = fi->m_f[ 0 ]; + if(fi -> m_f[ 1 ] != invalid ) *out++ = fi->m_f[ 1 ]; + if(fi -> m_f[ 2 ] != invalid ) *out++ = fi->m_f[ 2 ]; +} + +template +< + typename Point, + bool ClockWise, + template class VertexContainer, + template class FaceContainer, + template class VertexAllocator, + template class FaceAllocator, + typename FaceIterator, + typename OutputIterator +> +inline void face_incident_faces( + model::triangulation & t, + FaceIterator fi, + OutputIterator out) +{ + typedef typename model::triangulation triangulation; + typedef typename triangulation::face_iterator face_iterator; + typedef typename triangulation::vertex_iterator vertex_iterator; + typedef typename triangulation::face_vertex_index face_vertex_index; + typedef typename triangulation::face_type face_type; + face_iterator const invalid = t.invalid(); + for(face_vertex_index i = 0; i < 3; ++i) + { + face_iterator n = fi -> m_f[i]; + face_iterator m = fi -> m_f[i == 2 ? 0 : i + 1]; + face_iterator f_prev = fi; + face_vertex_index v_prev = i; + if(n != invalid) + { + f_prev = n; + v_prev = fi -> m_o[i]; + v_prev = (v_prev == 0 ? 2 : v_prev - 1); + } + while(true) + { + face_iterator next = f_prev -> m_f[v_prev]; + if(next == invalid) + { + face_vertex_index j = (i == 2 ? 0 : i + 1); + face_iterator m = fi -> m_f[j]; + if(m == invalid) break; + face_vertex_index prev_vertex_index = (i == 0 ? 2 : i - 1 ); + vertex_iterator const& prev_vertex_it = + fi->m_v[ prev_vertex_index ]; + face_type& first = *prev_vertex_it->m_f; + *out++ = prev_vertex_it->m_f; + f_prev = next = prev_vertex_it->m_f; + if(first.m_v[0] == prev_vertex_it) v_prev = 1; + else if(first.m_v[1] == prev_vertex_it) v_prev = 2; + else v_prev = 0; + continue; + } else { + if(&(*next) == &(*fi)) break; + *out++ = next; + v_prev = f_prev -> m_o[v_prev]; + v_prev = (v_prev == 0 ? 2 : v_prev - 1); + f_prev = next; + } + } + } +} + +template +< + typename Point, + bool ClockWise, + template class VertexContainer, + template class FaceContainer, + template class VertexAllocator, + template class FaceAllocator, + typename OutputIterator +> +inline void vertex_incident_faces( + model::triangulation & t, + typename model::triangulation::vertex_iterator vi, + OutputIterator out) +{ + typedef typename model::triangulation triangulation; + typedef typename triangulation::face_iterator face_iterator; + typedef typename triangulation::fulledge_index fulledge_index; + fulledge_index e = t.begin_vertex_edge(vi); + face_iterator first_face = e.m_f2; + *out++ = first_face; + while(true) + { + e = t.next_around_vertex(e); + if(e.m_f2 == t.invalid() || e.m_f2 == first_face) break; + *out++ = e.m_f2; + } +} + +template +< + typename Point, + bool ClockWise, + template class VertexContainer, + template class FaceContainer, + template class VertexAllocator, + template class FaceAllocator, + typename OutputIterator +> +inline void vertex_incident_vertices( + model::triangulation & t, + typename model::triangulation::vertex_iterator vi, + OutputIterator out) +{ + typedef model::triangulation triangulation_type; + typedef typename triangulation_type::face_iterator face_iterator; + typedef typename triangulation_type::fulledge_index fulledge_index; + fulledge_index e = t.begin_vertex_edge(vi); + face_iterator first_face = e.m_f2; + *out++ = e.m_f2 -> m_v[ e.m_v2 == 0 ? 2 : e.m_v2 - 1 ]; + while(true) + { + e = t.next_around_vertex(e); + if(e.m_f2 == t.invalid() || e.m_f2 == first_face) break; + *out++ = e.m_f2 -> m_v[ e.m_v2 == 0 ? 2 : e.m_v2 - 1 ]; + } +} + +} // namespace geometry + +} // namespace boost + +#endif // BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_GEOMETRIES_TRIANGULATION_HPP diff --git a/include/boost/geometry/extensions/triangulation/geometries/voronoi_adaptor.hpp b/include/boost/geometry/extensions/triangulation/geometries/voronoi_adaptor.hpp new file mode 100644 index 0000000000..07d6318718 --- /dev/null +++ b/include/boost/geometry/extensions/triangulation/geometries/voronoi_adaptor.hpp @@ -0,0 +1,132 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2019 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2019 program. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_GEOMETRIES_VORONOI_ADAPTOR_HPP +#define BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_GEOMETRIES_VORONOI_ADAPTOR_HPP + +#include +#include + +namespace boost { namespace geometry +{ + +namespace model +{ + +template +< + typename Triangulation +> +class voronoi_vertex_view { +public: + typedef typename Triangulation::point_type point_type; + typedef typename Triangulation::face_iterator face_iterator; + voronoi_vertex_view(face_iterator f): + m_p(strategy::delaunay_triangulation::detail::circumcircle_center + < + point_type, + boost::geometry::point_order::value == clockwise + >( *f->begin(), + *(f->begin() + 1), + *(f->begin() + 2))), + m_f(f) {} + point_type const* point() const { return m_p; } + face_iterator base() const { return m_f; } + point_type m_p; +private: + face_iterator m_f; +}; + +template +< + typename Triangulation +> +class voronoi_face_view { +public: + typedef typename Triangulation::point_type point_type; + typedef typename Triangulation::vertex_iterator vertex_iterator; + typedef typename Triangulation::face_iterator face_iterator; + typedef typename std::vector>::const_iterator const_iterator; + typedef typename std::vector>::iterator iterator; + + typedef voronoi_vertex_view voronoi_vertex; + voronoi_face_view(Triangulation& t, vertex_iterator v):m_v(v) { + struct transforming_back_insert { + transforming_back_insert(std::back_insert_iterator> bi): + m_underlying(bi) {}; + std::back_insert_iterator> m_underlying; + transforming_back_insert& operator*() { return *this; } + transforming_back_insert& operator++() { return *this; } + transforming_back_insert& operator++(int) { return *this; } + transforming_back_insert& operator=(face_iterator const& f) + { + m_underlying = voronoi_vertex(f); + return *this; + } + }; + vertex_incident_faces(t, m_v, transforming_back_insert(std::back_inserter(m_voronoi_v))); + } + const_iterator begin() const { return m_voronoi_v.begin(); } + const_iterator end() const { return m_voronoi_v.end(); } +private: + std::vector m_voronoi_v; + vertex_iterator m_v; +}; + +} // namespace model + +#ifndef DOXYGEN_NO_TRAITS_SPECIALIZATIONS +namespace traits +{ + +template +struct dimension< model::voronoi_vertex_view > : boost::mpl::int_<2> {}; + +template +struct tag< model::voronoi_vertex_view > +{ typedef point_tag type; }; + +template +struct coordinate_type< model::voronoi_vertex_view > +{ typedef typename coordinate_type::type type; }; + +template +struct coordinate_system< model::voronoi_vertex_view > +{ typedef typename coordinate_system::type type; }; + +template +struct access, Dimension> +{ + static typename coordinate_type::type + get(model::voronoi_vertex_view const& p) + { + return boost::geometry::get(p.m_p); + } +}; + +template struct tag> +{ typedef ring_tag type; }; + +template struct point_order> +{ static const order_selector value = counterclockwise; }; + +template +struct closure> +{ + static const closure_selector value = open; +}; + +} // namespace traits +#endif // DOXYGEN_NO_TRAITS_SPECIALIZATIONS + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_GEOMETRIES_VORONOI_ADAPTOR_HPP diff --git a/include/boost/geometry/extensions/triangulation/strategies/cartesian/accelerated_shull.hpp b/include/boost/geometry/extensions/triangulation/strategies/cartesian/accelerated_shull.hpp new file mode 100644 index 0000000000..f57d42f41e --- /dev/null +++ b/include/boost/geometry/extensions/triangulation/strategies/cartesian/accelerated_shull.hpp @@ -0,0 +1,63 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2019 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2019 program. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_ACCELERATED_SHULL_HPP +#define BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_ACCELERATED_SHULL_HPP + +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace boost { namespace geometry { namespace strategy +{ + +namespace delaunay_triangulation +{ + +struct accelerated_shull +{ + template + < + typename PointContainer, + typename Triangulation, + typename SideStrategy = strategy::side::side_by_triangle<>, + typename InCircleStrategy = strategy::in_circle::in_circle_robust<> + > + static inline void apply(PointContainer const & in, + Triangulation& out, + bool legalize = true) + { + detail::accelerated_shull::apply< + PointContainer, + Triangulation, + SideStrategy, + InCircleStrategy>(in, out, legalize); + } +}; + +namespace services +{ +template <> +struct default_strategy +{ + typedef accelerated_shull type; +}; +} + +}}}} // namespace boost::geometry::strategy::delaunay:triangulation + +#endif // BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_ACCELERATED_SHULL_HPP diff --git a/include/boost/geometry/extensions/triangulation/strategies/cartesian/detail/accelerated_shull.hpp b/include/boost/geometry/extensions/triangulation/strategies/cartesian/detail/accelerated_shull.hpp new file mode 100644 index 0000000000..cf4494c3f0 --- /dev/null +++ b/include/boost/geometry/extensions/triangulation/strategies/cartesian/detail/accelerated_shull.hpp @@ -0,0 +1,448 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2019 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2019 program. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_DETAIL_ACCELERATED_SHULL_HPP +#define BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_DETAIL_ACCELERATED_SHULL_HPP + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace boost { namespace geometry { namespace strategy +{ + +#ifndef DOXYGEN_NO_DETAIL + +namespace delaunay_triangulation { namespace detail +{ + +template +inline CalculationType comparable_circumcircle_diameter( + const Point& p1, const Point& p2, const Point& p3) +{ + typedef typename coordinate_type::type coordinate_type; + CalculationType comp_area = SideStrategy::template + side_value(p1, p2, p3); + return comparable_distance(p1, p2) + * comparable_distance(p1, p3) + * comparable_distance(p2, p3) + / (comp_area * comp_area); +} + +template +< + typename PointOut, + bool ClockWise, + typename PointIn +> +inline PointOut circumcircle_center(const PointIn& p1, + const PointIn& p2, + const PointIn& p3) +{ + typedef typename coordinate_type::type coordinate_type; + coordinate_type ax = get<0>( p1 ); + coordinate_type ay = get<1>( p1 ); + coordinate_type bx = get<0>( !ClockWise ? p2 : p3 ); + coordinate_type by = get<1>( !ClockWise ? p2 : p3 ); + coordinate_type cx = get<0>( !ClockWise ? p3 : p2 ); + coordinate_type cy = get<1>( !ClockWise ? p3 : p2 ); + coordinate_type d = + 2 * ( ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)); + coordinate_type x = ( (ax * ax + ay * ay) * (by - cy) + + (bx * bx + by * by) * (cy - ay) + + (cx * cx + cy * cy) * (ay - by)) / d; + coordinate_type y = ( (ax * ax + ay * ay) * (cx - bx ) + + (bx * bx + by * by) * (ax - cx) + + (cx * cx + cy * cy) * (bx - ax)) / d; + return boost::geometry::make(x, y); +} + +/* + * The following algorithm is based on s-hull by David Sinclair. More details can be found at + * http://s-hull.org. A more detailed description can be found in the paper: + * http://www.s-hull.org/paper/s_hull.pdf. The segments that correspond to the steps in the paper + * are marked by comments. An attempt has been made to speed up the algorithm by using a binary + * search for the first and last visible edge in step 7. This reduces the number of calls to the + * orientation test (SideStrategy). + */ +struct accelerated_shull +{ + template + < + typename PointContainer, + typename Triangulation, + typename SideStrategy, + typename InCircleStrategy, + typename CalculationType = double + > + static inline void apply(PointContainer const & in, + Triangulation& out, + bool legalize) + { + typedef typename PointContainer::value_type point_type; + typedef typename default_distance_result::type + distance_type; + const bool is_cw = + boost::geometry::point_order::value + == clockwise; + typedef CalculationType ct; + std::vector> points; + points.reserve(in.size()); + //Step 1 & 2 & 3 + point_type zero; + set<0>(zero, 0); + set<1>(zero, 0); + std::transform(std::begin(in), std::end(in), std::back_inserter(points), + [&points, &zero](point_type const& p) + { + return std::tuple(p, + boost::geometry::comparable_distance(p, zero), + 0); + }); + std::sort(std::begin(points), + std::end(points), + [](std::tuple const& p0, + std::tuple const& p1) + { return std::get<1>(p0) < std::get<1>(p1); }); + for(std::tuple& p : points) + { + std::get<1>(p) = boost::geometry::comparable_distance(std::get<0>(p), + std::get<0>(points[ 0 ])); + } + std::sort(std::begin(points) + 1, + std::end(points), + [](std::tuple const& p0, + std::tuple const& p1) + { return std::get<1>(p0) < std::get<1>(p1); }); + //Step 4 + { + distance_type min_circumdiameter = + std::numeric_limits::max(); + std::size_t min_index = 2; + for( std::size_t i = 2 ; + std::get<1>(points[i]) < min_circumdiameter && i < points.size() ; + ++i ) + { + distance_type diam = + comparable_circumcircle_diameter( + std::get<0>(points[ 0 ]), + std::get<0>(points[ 1 ]), + std::get<0>(points[ i ])); + if(diam < min_circumdiameter) + { + min_index = i; + min_circumdiameter = diam; + } + } + std::swap(points[ 2 ], points[ min_index ]); + } + //Step 5 + if( ( !is_cw + && SideStrategy::apply(std::get<0>(points[ 0 ]), + std::get<0>(points[ 1 ]), + std::get<0>(points[ 2 ])) < 0 ) + || ( is_cw + && SideStrategy::apply(std::get<0>(points[ 0 ]), + std::get<0>(points[ 1 ]), + std::get<0>(points[ 2 ])) > 0 )) + { + std::swap(points[ 1 ], points[ 2 ]); + } + const double PI = 3.14159265358979323846; + //Step 6 + + typedef typename Triangulation::vertex_iterator vertex_iterator; + typedef typename Triangulation::face_iterator face_iterator; + vertex_iterator const p1 = out.add_vertex(std::get<0>(points[ 0 ])); + vertex_iterator const p2 = out.add_vertex(std::get<0>(points[ 1 ])); + vertex_iterator const p3 = out.add_vertex(std::get<0>(points[ 2 ])); + face_iterator const seed_face = out.add_isolated_face(p1, p2, p3); + point_type cen; + set<0>(cen, + ( get<0>(std::get<0>(points[ 0 ])) + + get<0>(std::get<0>(points[ 1 ])) + + get<0>(std::get<0>(points[ 2 ])) ) / 3); + set<1>(cen, + ( get<1>(std::get<0>(points[ 0 ])) + + get<1>(std::get<0>(points[ 1 ])) + + get<1>(std::get<0>(points[ 2 ])) ) / 3); + + point_type C = circumcircle_center( + std::get<0>(points[ 0 ]), + std::get<0>(points[ 1 ]), + std::get<0>(points[ 2 ])); + { + std::for_each(std::begin(points), std::end(points), + [&C, &cen, &PI](std::tuple& p){ + std::get<1>(p) = comparable_distance(std::get<0>(p), C); + std::get<2>(p) = + ! is_cw ? + std::atan2( + get<1>(std::get<0>(p)) - get<1>(cen), + get<0>(std::get<0>(p)) - get<0>(cen)) + PI + : 2 * PI - (std::atan2( + get<1>(std::get<0>(p)) - get<1>(cen), + get<0>(std::get<0>(p)) - get<0>(cen)) + PI); + }); + std::sort(std::begin(points) + 3, + std::end(points), + [](std::tuple const& p0, + std::tuple const& p1) + { return std::get<1>(p0) < std::get<1>(p1); }); + } + //Step 7 & 8 + typedef typename Triangulation::halfedge_index halfedge_index; + typedef typename Triangulation::point_type out_point_type; + typedef typename Triangulation::segment_type segment_type; + { + halfedge_index e1 = out.face_edge(seed_face); + halfedge_index e2 = out.next(e1); + halfedge_index e3 = out.next(e2); + typedef std::pair + convex_hull_edge; + typedef std::vector convex_hull_container; + typedef typename convex_hull_container::iterator convex_hull_iterator; + convex_hull_container convex_hull { + std::make_pair(e1, ct(std::get<2>(points[1]))), + std::make_pair(e2, ct(std::get<2>(points[2]))), + std::make_pair(e3, ct(std::get<2>(points[0]))) }; + for (std::size_t i = 3; i < points.size(); ++i) + { + vertex_iterator new_vertex = + out.add_vertex(std::get<0>(points[ i ])); + out_point_type const& p = out.vertex(new_vertex); + auto is_visible = [&out, &p](convex_hull_edge const& bep) + { + const halfedge_index& be = std::get<0>(bep); + const segment_type s = out.face_segment(be); + const out_point_type p1 = s.first; + const out_point_type p2 = s.second; + auto det = SideStrategy::apply(p1,p2,p); + bool result = (!is_cw && det < 0) || (is_cw && det > 0); + return result; + }; + const ct& ref_angle = std::get<1>(*convex_hull.begin()); + ct angle = std::get<2>(points[ i ]) - ref_angle; + if(angle < 0) angle += 2 * PI; + ct opposite = angle - PI; + if(opposite < 0) opposite += 2 * PI; + + auto pred = [&ref_angle, &PI] + (convex_hull_edge const& be, ct const& a) + { + ct ba = std::get<1>(be) - ref_angle; + if(ba < 0) ba += 2*PI; + return ba < a; + }; + convex_hull_iterator vis_edge = std::lower_bound( + convex_hull.begin(), convex_hull.end(), + angle, pred); + if( vis_edge == convex_hull.begin() ) + vis_edge = convex_hull.end() - 1; + else --vis_edge; + convex_hull_iterator invis_edge = std::lower_bound( + convex_hull.begin(), convex_hull.end(), + opposite, pred); + if( invis_edge == convex_hull.begin() ) + invis_edge = convex_hull.end() - 1; + else --invis_edge; + if(!is_visible(*vis_edge)) + { + invis_edge = vis_edge; + for(auto it = convex_hull.begin(); + it != convex_hull.end(); + ++it) + { + if(is_visible(*it)) + { + vis_edge = it; + break; + } + } + } + auto vis_small = [&is_visible] + (convex_hull_edge const& be, int const&) + { + if(is_visible(be)) return true; + else return false; + }; + auto vis_large = [&is_visible] + (convex_hull_edge const& be, int const&) + { + if(is_visible(be)) return false; + else return true; + }; + convex_hull_iterator first_visible, last_visible, + fv2 = convex_hull.end(); + if(vis_edge < invis_edge) + { + if(is_visible(*convex_hull.begin())) + { + first_visible = convex_hull.begin(); + last_visible = std::lower_bound(vis_edge, invis_edge, + 0, vis_small); + fv2 = std::lower_bound(invis_edge, convex_hull.end(), + 0, vis_large); + } + else + { + first_visible = std::lower_bound(convex_hull.begin(), + vis_edge, 0, vis_large); + last_visible = std::lower_bound(vis_edge, + convex_hull.end(), 0, vis_small); + } + } + else + { + if(is_visible(*convex_hull.begin())) + { + first_visible = convex_hull.begin(); + last_visible = std::lower_bound( + convex_hull.begin()+1, invis_edge, + 0, vis_small); + fv2 = std::lower_bound(invis_edge, vis_edge, + 0, vis_large); + } + else + { + first_visible = std::lower_bound(invis_edge, vis_edge, + 0, vis_large); + last_visible = std::lower_bound(vis_edge, + convex_hull.end(), 0, vis_small); + } + } + ct first_visible_angle = std::get<1>(*first_visible); + const bool begin_visible = first_visible == std::begin(convex_hull); + face_iterator prev, fnf2; + bool looped = false; + ct fv2_angle; + if(begin_visible && is_visible(convex_hull.back())) + { + looped = true; + fv2_angle = std::get<1>(*fv2); + fnf2 = out.add_face_on_boundary(std::get<0>(*fv2), new_vertex); + prev = fnf2; + for(convex_hull_iterator it = fv2 + 1 ; + it != std::end(convex_hull) ; + ++it) + { + face_iterator next = + out.add_face_on_boundary(std::get<0>(*it), new_vertex); + out.connect(out.next(out.face_edge(next)), + out.prev(out.face_edge(prev))); + prev = next; + } + convex_hull.erase(fv2, std::end(convex_hull)); + } + const face_iterator fnf = + out.add_face_on_boundary(std::get<0>(*first_visible), + new_vertex); + if(looped) + { + out.connect(out.next(out.face_edge(fnf)), + out.prev(out.face_edge(prev))); + } + prev = fnf; + for (convex_hull_iterator it = first_visible + 1; + it != last_visible; + ++it) + { + face_iterator next = + out.add_face_on_boundary(std::get<0>(*it), new_vertex); + out.connect(out.next(out.face_edge(next)), + out.prev(out.face_edge(prev))); + prev = next; + } + convex_hull_iterator ip = + convex_hull.erase(first_visible, last_visible); + ip = convex_hull.insert(ip, + std::make_pair(out.prev(out.face_edge(prev)), + ct(std::get<2>(points[i])))); + if(looped) + { + convex_hull.insert(ip, + std::make_pair(out.next(out.face_edge(fnf2)), + ct(fv2_angle))); + } + else + { + convex_hull.insert(ip, + std::make_pair(out.next(out.face_edge(fnf)), + ct(first_visible_angle))); + } + } + } + //Step 9 + if(legalize) + { + typedef typename Triangulation::face_vertex_index face_vertex_index; + std::vector L; + for(face_iterator i = out.faces_begin(); i != out.faces_end(); ++i) + { + for(face_vertex_index j = 0; j<2 ; ++j) + if( out.neighbour(i, j)!=Triangulation::invalid() + && &(*i) > &(*out.neighbour(i, j))) + { + L.push_back(halfedge_index(i,j)); + } + } + auto edge_legal = [&out](halfedge_index const& e) + { + out_point_type const& p1 = out.face_vertex(e.m_f, 0); + out_point_type const& p2 = out.face_vertex(e.m_f, 1); + out_point_type const& p3 = out.face_vertex(e.m_f, 2); + out_point_type const& p = out.face_vertex( + out.neighbour(e.m_f, e.m_v), + out.opposite(e.m_f, e.m_v) + ); + return !is_cw ? + InCircleStrategy::apply(p1, p2, p3, p) <= 0 + || InCircleStrategy::apply(p1, p2, p, p3) > 0 + : InCircleStrategy::apply(p1, p3, p2, p) <= 0 + || InCircleStrategy::apply(p1, p3, p, p2) > 0; + }; + while( !L.empty() ) + { + halfedge_index e = L.back(); + L.pop_back(); + face_iterator const f1 = e.m_f; + face_vertex_index const v1 = e.m_v; + face_iterator const f2 = out.neighbour(f1, v1); + if(f2 == face_iterator()) continue; + face_vertex_index const v2 = out.opposite(f1, v1); + if( !edge_legal(e) ) + { + L.emplace_back(f1, v1); + L.emplace_back(f2, v2); + L.emplace_back(f1, v1 == 0 ? 2 : v1 - 1); + L.emplace_back(f2, v2 == 0 ? 2 : v2 - 1); + out.flip(e); + } + } + } + } +}; + +}} // namespace delaunay_triangulation::detail + +#endif // DOXYGEN_NO_DETAIL + +}}} // namespace boost::geometry::strategy + +#endif // BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_DETAIL_ACCELERATED_SHULL_HPP diff --git a/include/boost/geometry/extensions/triangulation/strategies/delaunay_triangulation.hpp b/include/boost/geometry/extensions/triangulation/strategies/delaunay_triangulation.hpp new file mode 100644 index 0000000000..31e57a9a57 --- /dev/null +++ b/include/boost/geometry/extensions/triangulation/strategies/delaunay_triangulation.hpp @@ -0,0 +1,41 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2019 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2019 program. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_DELAUNAY_TRIANGULATION_HPP +#define BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_DELAUNAY_TRIANGULATION_HPP + +#include + +namespace boost { namespace geometry { namespace strategy +{ + +namespace delaunay_triangulation { namespace services +{ + +/*! +\brief Traits class binding a default delaunay triangulation construction strategy to a coordinate system +\ingroup util +\tparam Tag tag of geometry +\tparam CSTag tag of coordinate system +*/ +template +struct default_strategy +{ + BOOST_MPL_ASSERT_MSG + ( + false, NOT_IMPLEMENTED_FOR_THIS_TYPE + , (types) + ); +}; + +}}}}} // namesspace boost::geometry::strategy::delaunay_triangulation::services + +#endif // BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_DELAUNAY_TRIANGULATION_HPP diff --git a/include/boost/geometry/extensions/triangulation/triangulation.hpp b/include/boost/geometry/extensions/triangulation/triangulation.hpp new file mode 100644 index 0000000000..fcd924b208 --- /dev/null +++ b/include/boost/geometry/extensions/triangulation/triangulation.hpp @@ -0,0 +1,18 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2019 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2019 program. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_HPP +#define BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_HPP + +#include +#include + +#endif // BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_HPP