diff --git a/src/CGAL_typedefs.h b/src/CGAL_typedefs.h
new file mode 100644
index 0000000000000000000000000000000000000000..7c88a59410f3c3d43f564847774f0ab739bbbf0a
--- /dev/null
+++ b/src/CGAL_typedefs.h
@@ -0,0 +1,97 @@
+// compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#ifndef CGAL_TYPEDEFS_H
+#define CGAL_TYPEDEFS_H
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Linear_cell_complex_for_combinatorial_map.h>
+#include <CGAL/Compact_container.h>
+#include <CGAL/Real_timer.h>
+#include <CGAL/Timer.h>
+#include <vector>
+
+template<typename LCC>
+class Inner_point;
+
+///////////////////////////////////////////////////////////////////////////////
+typedef CGAL::Exact_predicates_inexact_constructions_kernel EPICK;
+typedef CGAL::Exact_predicates_exact_constructions_kernel EPECK;
+typedef EPECK::FT ECoord;
+typedef EPECK::Point_2 EPoint;
+typedef EPECK::Segment_2 ESegment;
+typedef EPECK::Iso_rectangle_2 ERectangle;
+typedef EPICK::FT ICoord;
+typedef EPICK::Point_2 IPoint;
+typedef EPICK::Segment_2 ISegment;
+typedef EPICK::Iso_rectangle_2 IRectangle;
+///////////////////////////////////////////////////////////////////////////////
+// typedef CGAL::Timer My_timer;
+typedef CGAL::Real_timer My_timer;
+///////////////////////////////////////////////////////////////////////////////
+template<typename Ref>
+struct Info_for_face
+{
+ using Face_attribute_handle=
+ typename Ref::template Attribute_handle<2>::type;
+
+ Info_for_face() : m_finite(true), m_father(nullptr)
+ {}
+
+ bool m_finite;
+ std::vector<Face_attribute_handle> m_son;
+ Face_attribute_handle m_father;
+};
+////////////////////////////////////////////////////////////////////////////////
+template<typename Ref>
+struct Info_for_vertex
+{
+ using Dart_handle=typename Ref::Dart_handle;
+
+ Info_for_vertex(Dart_handle dh=nullptr) : m_dart_above(dh)
+ {}
+
+ Dart_handle m_dart_above;
+};
+////////////////////////////////////////////////////////////////////////////////
+struct Myitem
+{
+ template<class Ref>
+ struct Dart_wrapper
+ {
+ typedef CGAL::Cell_attribute_with_point<Ref, Info_for_vertex<Ref>> Vertex_attribute;
+ typedef CGAL::Cell_attribute<Ref, Info_for_face<Ref>> Face_attribute;
+ typedef std::tuple<Vertex_attribute,void,Face_attribute> Attributes;
+ };
+};
+typedef CGAL::Linear_cell_complex_traits
+<2, CGAL::Exact_predicates_exact_constructions_kernel> Traits;
+ typedef CGAL::Linear_cell_complex_for_combinatorial_map
+ <2,2,Traits,Myitem,CGAL_ALLOCATOR(int),
+ CGAL::Combinatorial_map_base,
+ CGAL::CMap_linear_cell_complex_storage_1
+ <2, 2, Traits, Myitem, CGAL_ALLOCATOR(int), CGAL::Tag_true>
+> LCC_2;
+///////////////////////////////////////////////////////////////////////////////
+typedef LCC_2::Dart_handle Dart_handle;
+typedef LCC_2::Dart_const_handle Dart_const_handle;
+///////////////////////////////////////////////////////////////////////////////
+
+#endif // CGAL_TYPEDEFS_H
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c5e704758aea7cfef3b53cb99d1155039dce23c8
--- /dev/null
+++ b/src/CMakeLists.txt
@@ -0,0 +1,81 @@
+project(Parallel_arrangement)
+
+cmake_minimum_required(VERSION 3.1)
+if("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}" VERSION_GREATER 2.6)
+ if("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}.${CMAKE_PATCH_VERSION}" VERSION_GREATER 2.8.3)
+ cmake_policy(VERSION 2.8.4)
+ else()
+ cmake_policy(VERSION 2.6)
+ endif()
+endif()
+
+cmake_policy(SET CMP0064 NEW)
+
+set( CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true )
+if ( COMMAND cmake_policy )
+ cmake_policy( SET CMP0003 NEW )
+endif()
+
+# ADD_DEFINITIONS("-Wall -Wextra")
+# ADD_DEFINITIONS("-O3")
+# -fsanitize=thread
+
+set(CMAKE_CXX_STANDARD 17)
+find_package(Boost COMPONENTS system filesystem REQUIRED)
+
+# ##########################################################
+# Find some packages
+
+find_package(CGAL COMPONENTS Qt5 REQUIRED) # CGAL and QT5 component
+if(CGAL_Qt5_FOUND)
+ add_definitions(-DCGAL_USE_BASIC_VIEWER -DQT_NO_KEYWORDS)
+endif()
+
+# ##########################################################
+## For profilling with gprof
+# add_definitions("-pg")
+# SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -pg")
+# add_definitions("-D_GLIBCXX_DEBUG")
+add_definitions("-g") # for vtune (perf evaluator)
+
+# ##########################################################
+include_directories("${CMAKE_CURRENT_SOURCE_DIR}/")
+
+set(HEADER_FILES
+ CGAL_typedefs.h
+ draw_shds.h
+ Export_xfig.hh
+ Hds.h
+ memory.h
+ My_arr_segment_traits_2.h
+ My_construction_subcurve.h
+ My_event.h
+ My_surface_sweep.h
+ My_visitor.h
+ Partial_hds.h
+ SHds.h
+ SHds_to_lcc.h
+ SHds_to_cgal_arrangement.h
+ Segment_readers.h
+ Print_txt.h)
+
+set(SOURCE_FILES
+ parallel-arrangement-bench.cpp
+ parallel-arrangement.cpp
+ compute-arrangement-streamed.cpp)
+
+foreach(cppfile ${SOURCE_FILES})
+ create_single_source_cgal_program("${cppfile}" "${HEADER_FILES}")
+endforeach()
+# ##########################################################
+# Link
+if (CGAL_Qt5_FOUND)
+ target_link_libraries(parallel-arrangement PUBLIC CGAL::CGAL_Qt5)
+endif()
+
+find_package(TBB REQUIRED)
+include(CGAL_TBB_support)
+target_link_libraries(parallel-arrangement PUBLIC pthread CGAL::TBB_support ${Boost_FILESYSTEM_LIBRARY} ${Boost_SYSTEM_LIBRARY})
+target_link_libraries(parallel-arrangement-bench PUBLIC pthread CGAL::TBB_support ${Boost_FILESYSTEM_LIBRARY} ${Boost_SYSTEM_LIBRARY})
+target_link_libraries(compute-arrangement-streamed PUBLIC pthread CGAL::TBB_support ${Boost_FILESYSTEM_LIBRARY} ${Boost_SYSTEM_LIBRARY})
+# ##########################################################
diff --git a/src/Export_xfig.hh b/src/Export_xfig.hh
new file mode 100644
index 0000000000000000000000000000000000000000..43a34ec23fef540297df01ce15d19d58ae24e07d
--- /dev/null
+++ b/src/Export_xfig.hh
@@ -0,0 +1,205 @@
+// compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#ifndef EXPORT_XFIG_H
+#define EXPORT_XFIG_H
+
+#include <cmath>
+#include <string>
+#include <fstream>
+
+class Export_xfig
+{
+ constexpr static unsigned int MAX_PROF=999; // max depth for xfig
+
+public:
+ // xsize is the length of the bbox in x-coordinates
+ // We resize the scene so that the x-length is xfig_x_size.
+ Export_xfig(const std::string& filename, double xsize,
+ double xfig_x_size=20000.,
+ bool invert_y=false) :
+ m_nb_couls(0),
+ m_invert_y(invert_y),
+ m_x_factor(xfig_x_size/xsize),
+ m_x_shift(0),
+ m_y_shift(0)
+ {
+ m_ofstream.open(filename);
+ if (!m_ofstream)
+ { std::cerr<<"ERROR opening "<<filename<<std::endl; }
+
+ sauveEntete();
+ }
+
+ ~Export_xfig()
+ { m_ofstream.close(); }
+
+ void set_x_shift(int x)
+ { m_x_shift=x; }
+ void set_y_shift(int y)
+ { m_y_shift=y; }
+
+ uint8_t add_color(const CGAL::Color& c)
+ {
+ m_ofstream<<std::setfill('0');
+ m_ofstream<<"0 "<<32+m_nb_couls<<" #"
+ <<std::hex<<std::uppercase<<std::setw(2)
+ <<static_cast<int>(c.red())
+ <<static_cast<int>(c.green())
+ <<static_cast<int>(c.blue())<<std::dec<<std::endl;
+ ++m_nb_couls;
+ return m_nb_couls-1;
+ }
+
+ /* TODO void debutComposante( ofstream &os,
+ const CVertex & min, const CVertex & max )
+ {
+ m_ofstream<<"6 "
+ <<int(RINT(min.getX()+decalageX))<<" "
+ <<int(RINT(min.getY()+decalageY))<<" "
+ <<int(RINT(max.getX()+decalageX))<<" "
+ <<int(RINT(max.getY()+decalageY))<<endl;
+ }*/
+
+ // Fin d'une composante : m_ofstream<<"-6\n";
+
+
+ void point(double x, double y, uint16_t rayon,
+ u_int8_t coul, uint16_t larg, uint16_t prof)
+ {
+ m_ofstream<<"1 4 0 "<<static_cast<int>(larg)<<" "
+ <<static_cast<int>(32+coul)<<" "
+ <<static_cast<int>(32+coul)<<" "
+ <<static_cast<int>(prof);
+ m_ofstream <<" 0 20 0.000 1 0.0000 ";
+
+
+ m_ofstream<<point_to_string(x, y)<<" " // Centre x et y
+ <<static_cast<int>(rayon)<<" "<<static_cast<int>(rayon)<<" " //Rayon en x et en y
+ <<point_to_string(x, y)<<" " // Premier point == centre
+ <<point_to_string(x+rayon, y)<<" " // deuxieme point (sur le cercle)
+ <<std::endl;
+ }
+
+ void segment(double x1, double y1, double x2, double y2,
+ u_int8_t coul, uint16_t larg, bool arrow, uint16_t prof)
+ {
+ int ix1=p_to_x(x1);
+ int iy1=p_to_y(y1);
+ int ix2=p_to_x(x2);
+ int iy2=p_to_y(y2);
+
+ if (ix1!=ix2 || iy1!=iy2)
+ {
+ m_ofstream<<"2 1 0 "<<static_cast<int>(larg)<<" "
+ <<static_cast<int>(32+coul)<<" "
+ <<static_cast<int>(32+coul)<<" "
+ <<static_cast<int>(prof)
+ <<" 0 -1 0.000 1 0 7 ";
+ if (arrow)
+ {
+ m_ofstream<<"1 0 2\n"
+ <<" 2 0 1.00 60.00 120.00\n";
+ }
+ else
+ { m_ofstream<<"0 0 2\n"; }
+
+ m_ofstream<<ix1<<" "<<iy1<<" "<<ix2<<" "<<iy2<<std::endl;
+ }
+ }
+
+ void rectangle(double x1, double y1, double x2, double y2,
+ u_int8_t coul, uint16_t larg, bool arrow, uint16_t prof)
+ {
+ segment(x1, y1, x2, y1, coul, larg, arrow, prof);
+ segment(x2, y1, x2, y2, coul, larg, arrow, prof);
+ segment(x2, y2, x1, y2, coul, larg, arrow, prof);
+ segment(x1, y2, x1, y1, coul, larg, arrow, prof);
+ }
+
+protected:
+
+ int transfoMmPixel(double val) const
+ {
+ // Xfig est en 1200 dpi, et donc le nombre de pixel par mm vaut 1200/25.4
+ // (sachant que 1inch==25.4mm)
+ // double dpmX=(1200/25.4);
+ return std::round(val*m_x_factor); //(val*dpmX)/COEF_DE_REDUCTION);
+ }
+
+ int p_to_x(double val) const
+ { return transfoMmPixel(val)+m_x_shift; }
+ int p_to_y(double val) const
+ {
+ int res=transfoMmPixel(val)+m_y_shift;
+ if (m_invert_y) { res=-res; }
+ return res;
+ }
+
+ std::string value_to_string(double val)
+ { return std::to_string(transfoMmPixel(val)); }
+
+ std::string point_to_string(double x, double y)
+ {
+ int ix=p_to_x(x);
+ int iy=p_to_y(y);
+ return std::string(std::to_string(ix)+" "+std::to_string(iy));
+ }
+
+ void sauveEntete()
+ {
+ m_ofstream<<"#FIG 3.2"<<std::endl
+ <<"Portrait"<<std::endl
+ <<"Metric"<<std::endl
+ <<"A4"<<std::endl
+ <<"100.00"<<std::endl
+ <<"Single"<<std::endl
+ <<"-2"<<std::endl
+ <<"1200 2"<<std::endl;
+ }
+
+
+/* void CControlerGMap::sauveFace(const std::vector<Point2d>& pts,
+ int coulFill, int larg, int AProf )
+{
+ assert(pts.size()>=3); // 3 => polygone à 2 arêtes...
+
+
+ os<<"2 1 0 " << larg << " "
+ << 32+COUL_DART << " " << 32+coulFill << " "
+ << (1+AProf) <<" 0 20 0.000 1 0 7 0 0 "
+ << nbPts << endl;
+
+ for (int i=0; i<nbPts; ++i)
+ {
+ os << int(RINT(p[i].getX()+decalageX)) << " "
+ << int(RINT(p[i].getY()+decalageY)) << endl;
+ }
+}
+*/
+
+protected:
+ std::ofstream m_ofstream;
+ uint8_t m_nb_couls;
+ bool m_invert_y;
+ double m_x_factor;
+ int m_x_shift, m_y_shift;
+};
+////////////////////////////////////////////////////////////////////////////////
+#endif // EXPORT_XFIG_H
diff --git a/src/Hds.h b/src/Hds.h
new file mode 100644
index 0000000000000000000000000000000000000000..30a6c4c36da285888d9c9865866f5d979b10490f
--- /dev/null
+++ b/src/Hds.h
@@ -0,0 +1,525 @@
+// compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#ifndef HDS_H
+#define HDS_H
+
+#include <bitset>
+
+#include "CGAL_typedefs.h"
+#include "Print_txt.h"
+#include "Export_xfig.hh"
+
+///////////////////////////////////////////////////////////////////////////////
+class Halfedge;
+class Vertex;
+class Face;
+class HDS;
+class SHds;
+////////////////////////////////////////////////////////////////////////////////
+typedef CGAL::Compact_container<Halfedge>::iterator Halfedge_handle;
+typedef CGAL::Compact_container<Halfedge>::const_iterator Halfedge_const_handle;
+typedef CGAL::Compact_container<Vertex>::iterator Vertex_handle;
+typedef CGAL::Compact_container<Face>::iterator Face_handle;
+////////////////////////////////////////////////////////////////////////////////
+class Halfedge
+{
+public:
+ friend class HDS;
+ friend class Vertex;
+ friend class SHds;
+
+ typedef uint32_t size_type;
+ static const size_type NB_MARKS = 32;
+
+ Halfedge(Halfedge_handle opposite=nullptr,
+ Vertex_handle vertex=nullptr,
+ int edge_id=0) :
+ m_opposite(opposite),
+ m_prev_around_vertex(nullptr),
+ m_next_around_vertex(nullptr),
+ m_vertex(vertex),
+ m_face(nullptr),
+ m_edge_id(edge_id)
+ {}
+
+ void set( Halfedge_handle opposite,
+ Halfedge_handle prev_around_vertex,
+ Halfedge_handle next_around_vertex,
+ Vertex_handle vertex,
+ int edge_id)
+ {
+ m_opposite=opposite;
+ m_prev_around_vertex=prev_around_vertex;
+ m_next_around_vertex=next_around_vertex;
+ m_vertex=vertex;
+ m_edge_id=edge_id;
+ }
+
+ Halfedge_handle opposite() const
+ { return m_opposite; }
+ Halfedge_handle prev_around_vertex() const
+ { return m_prev_around_vertex; }
+ Halfedge_handle next_around_vertex() const
+ { return m_next_around_vertex; }
+ Vertex_handle vertex() const
+ { return m_vertex; }
+ Face_handle face() const
+ { return m_face; }
+ int edge_id() const
+ { return m_edge_id; }
+
+ bool get_mark(size_type amark) const
+ {
+ CGAL_assertion(amark>=0 && amark<NB_MARKS);
+ return m_marks[amark];
+ }
+ void set_mark(size_type amark, bool avalue) const
+ {
+ CGAL_assertion(amark>=0 && amark<NB_MARKS);
+ m_marks.set(amark, avalue);
+ }
+ void flip_mark(size_type amark) const
+ {
+ CGAL_assertion(amark>=0 && amark<NB_MARKS);
+ m_marks.flip(amark);
+ }
+ std::bitset<NB_MARKS> get_marks() const
+ { return m_marks; }
+ void set_marks(const std::bitset<NB_MARKS>& amarks) const
+ { m_marks=amarks; }
+
+ void* for_compact_container() const
+ { return m_opposite.for_compact_container(); }
+ void for_compact_container(void* p)
+ { m_opposite.for_compact_container(p);}
+
+protected:
+ Halfedge_handle m_opposite;
+ Halfedge_handle m_prev_around_vertex;
+ Halfedge_handle m_next_around_vertex;
+ Vertex_handle m_vertex;
+ Face_handle m_face;
+ int m_edge_id; // >0 if halfedge from left to right or vertical and bottom to top, <0 otherwise
+ // two halfedges of the same edge have the same value and opposite sign
+ mutable std::bitset<NB_MARKS> m_marks;
+};
+////////////////////////////////////////////////////////////////////////////////
+class Global_halfedge
+{
+public:
+ friend class SHds;
+
+ Global_halfedge(std::size_t strip_id=0,
+ Halfedge_handle hh=nullptr):
+ m_strip_id(strip_id),
+ m_hh(hh)
+ {}
+
+ void set(std::size_t strip_id=0, Halfedge_handle hh=nullptr)
+ { m_strip_id=strip_id; m_hh=hh; }
+
+ void set_halfedge(Halfedge_handle hh)
+ { m_hh=hh; }
+
+ std::size_t strip_id() const
+ { return m_strip_id; }
+
+ const Halfedge_handle& halfedge() const
+ { return m_hh; }
+ Halfedge_handle& halfedge()
+ { return m_hh; }
+
+ void go_to_left_strip()
+ { assert(m_strip_id>0); --m_strip_id; }
+ void go_to_right_strip()
+ { ++m_strip_id; }
+
+ bool operator==(const Global_halfedge& other) const
+ { return m_hh==other.m_hh; } // Handle are unique, thus we don't need to compare m_strip_id
+ bool operator!=(const Global_halfedge& other) const
+ { return !operator==(other); }
+
+ bool operator<(const Global_halfedge& other) const
+ { return m_strip_id<other.m_strip_id ||
+ (m_strip_id==other.m_strip_id && m_hh<other.m_hh); }
+ bool operator>(const Global_halfedge& other) const
+ { return other.operator<(*this); }
+ bool operator>=(const Global_halfedge& other) const
+ { return !operator<(other); }
+ bool operator<=(const Global_halfedge& other) const
+ { return !operator>(other); }
+
+protected:
+ std::size_t m_strip_id;
+ Halfedge_handle m_hh;
+};
+////////////////////////////////////////////////////////////////////////////////
+class Vertex
+{
+public:
+ friend class HDS;
+ friend class Halfedge;
+
+ Vertex() :
+ m_first_halfedge(nullptr),
+ m_above_halfedge(nullptr)
+ {}
+
+ Vertex(const EPoint& p) :
+ m_point(p),
+ m_first_halfedge(nullptr),
+ m_above_halfedge(nullptr)
+ {}
+
+ void set(const EPoint& p, Halfedge_handle first_halfedge,
+ Halfedge_handle above_halfedge)
+ {
+ m_point=p;
+ m_first_halfedge=first_halfedge;
+ m_above_halfedge=above_halfedge;
+ }
+
+ const EPoint& point() const
+ { return m_point; }
+
+ Halfedge_handle first_halfedge() const
+ { return m_first_halfedge; }
+
+ Halfedge_handle above_halfedge() const
+ { return m_above_halfedge; }
+ void set_above_halfedge(Halfedge_handle dh)
+ { m_above_halfedge=dh; }
+
+ bool no_left_edge() const
+ {
+ return m_first_halfedge->m_edge_id>0; // If the first halfedge is from left to right
+ // or vertical and bottom to up, there is no left edge around the vertex
+ }
+
+ void serialize_point(std::ostream& os) const
+ {
+ CGAL::exact(m_point);
+ os<<std::setprecision(17)<<m_point;
+ }
+
+ void* for_compact_container() const
+ { return m_first_halfedge.for_compact_container(); }
+ void for_compact_container(void* p)
+ { m_first_halfedge.for_compact_container(p);}
+
+protected:
+ EPoint m_point;
+ Halfedge_handle m_first_halfedge;
+ Halfedge_handle m_above_halfedge; // Halfedge handle of the edge above this vertex (nullptr if no such edge)
+};
+////////////////////////////////////////////////////////////////////////////////
+class Face
+{
+public:
+ friend class HDS;
+ friend class Halfedge;
+ friend class SHds;
+
+ Face() :
+ m_strip_id(0),
+ m_first_halfedge(nullptr),
+ m_father(nullptr),
+ m_finite(true)
+ {}
+
+ void set(std::size_t strip_id,
+ Halfedge_handle first_halfedge,
+ Face_handle father,
+ bool finite)
+ {
+ m_strip_id=strip_id;
+ m_first_halfedge=first_halfedge;
+ m_father=father;
+ m_finite=finite;
+ }
+
+ std::size_t strip_id() const
+ { return m_strip_id; }
+
+ Halfedge_handle first_halfedge() const
+ { return m_first_halfedge; }
+
+ Face_handle father() const
+ { return m_father; }
+
+ bool finite() const
+ { return m_finite; }
+
+ const std::vector<Face_handle>& sons() const
+ { return m_son; }
+
+ void* for_compact_container() const
+ { return m_first_halfedge.for_compact_container(); }
+ void for_compact_container(void* p)
+ { m_first_halfedge.for_compact_container(p);}
+
+protected:
+ std::size_t m_strip_id; // strip id of the first halfedge
+ Halfedge_handle m_first_halfedge;
+ Face_handle m_father;
+ std::vector<Face_handle> m_son;
+ bool m_finite;
+};
+////////////////////////////////////////////////////////////////////////////////
+class HDS
+{
+public:
+ typedef uint32_t size_type;
+
+ typedef CGAL::Compact_container<Halfedge> Halfedges_container;
+ typedef CGAL::Compact_container<Vertex> Vertices_container;
+
+ void clear()
+ {
+ m_halfedges.clear();
+ m_vertices.clear();
+ }
+
+ std::size_t number_of_halfedges() const
+ { return m_halfedges.size(); }
+ std::size_t number_of_vertices() const
+ { return m_vertices.size(); }
+
+ Halfedges_container& halfedges()
+ { return m_halfedges; }
+ Vertices_container& vertices()
+ { return m_vertices; }
+
+ const Halfedges_container& halfedges() const
+ { return m_halfedges; }
+ const Vertices_container& vertices() const
+ { return m_vertices; }
+
+ Halfedge_handle opposite(Halfedge_handle hh) const
+ { return hh->m_opposite; }
+ Halfedge_handle prev_around_vertex(Halfedge_handle hh) const
+ { return hh->m_prev_around_vertex; }
+ Halfedge_handle next_around_vertex(Halfedge_handle hh) const
+ { return hh->m_next_around_vertex; }
+ Vertex_handle vertex(Halfedge_handle hh) const
+ { return hh->m_vertex; }
+ Face_handle face(Halfedge_handle hh) const
+ { return hh->m_face; }
+ const EPoint& point(Vertex_handle vh) const
+ { return vh->point(); }
+ const EPoint& point(Halfedge_handle hh) const
+ { return point(vertex(hh)); }
+ int edge_id(Halfedge_handle hh) const
+ { return hh->m_edge_id; }
+
+ void move_to_opposite(Halfedge_handle& hh) const
+ { hh=opposite(hh); }
+ void move_to_prev_around_vertex(Halfedge_handle& hh) const
+ { hh=prev_around_vertex(hh); }
+ void move_to_next_around_vertex(Halfedge_handle& hh) const
+ { hh=next_around_vertex(hh); }
+
+ // @return true iff halfedge hh is from left to right
+ // (or vertical and from bottom to top)
+ bool left_to_right(Halfedge_handle hh) const
+ { return edge_id(hh)>0; }
+ // @return true iff halfedge hh is from right to left
+ bool right_to_left(Halfedge_handle hh) const
+ { return edge_id(hh)<0; }
+
+ Vertex_handle create_vertex(const EPoint& p)
+ { return m_vertices.emplace(p); }
+ void set_vertex(Halfedge_handle hh, Vertex_handle vh)
+ { hh->m_vertex=vh; }
+ void set_vertex(Halfedge_handle hh, const EPoint& p)
+ { set_vertex(hh, create_vertex(p)); }
+
+ Halfedge_handle create_edge(Vertex_handle vh, int segmentid)
+ {
+ assert(segmentid>0);
+
+ Halfedge_handle hh1=m_halfedges.emplace(nullptr, vh, segmentid);
+ Halfedge_handle hh2=m_halfedges.emplace(hh1, nullptr, -segmentid);
+ hh1->m_opposite=hh2;
+
+ return hh1;
+ }
+
+ void erase_halfedge(Halfedge_handle hh)
+ { m_halfedges.erase(hh); }
+
+ void erase_edge(Halfedge_handle hh)
+ {
+ erase_halfedge(hh->m_opposite);
+ erase_halfedge(hh);
+ }
+
+ Halfedge_handle add_segment_around_vertex(Halfedge_handle hh,
+ Vertex_handle vh,
+ Halfedge_handle prevhh)
+ {
+ if (prevhh==nullptr)
+ {
+ assert(vh->first_halfedge()==nullptr);
+ vh->m_first_halfedge=hh;
+ }
+ else
+ {
+ assert(vh->first_halfedge()!=nullptr);
+ prevhh->m_next_around_vertex=hh;
+ hh->m_prev_around_vertex=prevhh;
+ }
+ return hh;
+ }
+
+ void end_of_vertex(Vertex_handle vh, Halfedge_handle prevhh)
+ {
+ assert(vh->first_halfedge()!=nullptr);
+ {
+ prevhh->m_next_around_vertex=vh->first_halfedge();
+ vh->first_halfedge()->m_prev_around_vertex=prevhh;
+ }
+ }
+
+ friend std::ostream& operator<<(std::ostream& os,
+ const HDS& hds)
+ {
+ os<<hds.halfedges().size()<<" "<<hds.vertices().size()<<std::endl;
+ for (auto it=hds.vertices().begin(), itend=hds.vertices().end(); it!=itend; ++it)
+ {
+ assert(it->first_halfedge()!=nullptr);
+ it->serialize_point(os);
+ os<<" "<<hds.halfedges().index(it->first_halfedge())
+ <<" "<<(it->above_halfedge()==nullptr?'-':'+');
+ if (it->above_halfedge()!=nullptr)
+ { os<<" "<<hds.halfedges().index(it->above_halfedge()); }
+ os<<std::endl;
+ }
+
+ os<<std::endl;
+ for (auto it=hds.halfedges().begin(), itend=hds.halfedges().end(); it!=itend; ++it)
+ {
+ assert(it->opposite()!=nullptr);
+ assert(it->vertex()==nullptr || it->prev_around_vertex()!=nullptr);
+ assert(it->vertex()==nullptr || it->next_around_vertex()!=nullptr);
+
+ os<<" "<<hds.halfedges().index(it->opposite())
+ <<" "<<(it->vertex()!=nullptr?'+':'-');
+ if (it->vertex()!=nullptr)
+ {
+ os<<" "<<hds.halfedges().index(it->prev_around_vertex())
+ <<" "<<hds.halfedges().index(it->next_around_vertex())
+ <<" "<<hds.vertices().index(it->vertex());
+ }
+ os<<" "<<it->edge_id()<<std::endl;
+ }
+ return os;
+ }
+
+ void load(std::istream& is,
+ std::vector<Halfedge_handle>& halfedges,
+ std::vector<Vertex_handle>& vertices)
+ {
+ clear();
+ halfedges.clear();
+ vertices.clear();
+
+ std::size_t nb_halfedges, nb_vertices;
+ is>>nb_halfedges>>nb_vertices;
+
+ halfedges.reserve(nb_halfedges);
+ vertices.reserve(nb_vertices);
+
+ for (std::size_t i=0; i<nb_vertices; ++i)
+ { vertices.push_back(m_vertices.emplace()); }
+ for (std::size_t i=0; i<nb_halfedges; ++i)
+ { halfedges.push_back(m_halfedges.emplace()); }
+
+ EPoint p;
+ std::size_t ind1, ind2, ind3, ind4;
+ int int1;
+ char c;
+
+ for (std::size_t i=0; i<nb_vertices; ++i)
+ {
+ is>>p>>ind1>>c;
+ if (c=='+') is>>ind2;
+ vertices[i]->set(p, halfedges[ind1], (c=='+'?halfedges[ind2]:nullptr));
+ }
+ for (std::size_t i=0; i<nb_halfedges; ++i)
+ {
+ is>>ind1>>c;
+ if (c=='+') is>>ind2>>ind3>>ind4;
+ is>>int1;
+ halfedges[i]->set(halfedges[ind1],
+ (c=='+'?halfedges[ind2]:nullptr),
+ (c=='+'?halfedges[ind3]:nullptr),
+ (c=='+'?vertices[ind4]:nullptr), int1);
+ }
+ }
+
+ friend std::istream& operator>>(std::istream& is, HDS& hds)
+ {
+ std::vector<Halfedge_handle> halfedges;
+ std::vector<Vertex_handle> vertices;
+ hds.load(is, halfedges, vertices);
+ return is;
+ }
+
+ void export_to_xfig(Export_xfig& export_xfig,
+ double dx, double dy,
+ uint8_t coul_segments,
+ uint8_t coul_disks,
+ uint16_t width_segments,
+ uint16_t radius_disks) const
+ {
+ if (width_segments>0)
+ {
+ for (auto it=halfedges().begin(), itend=halfedges().end(); it!=itend; ++it)
+ {
+ if (it->vertex()!=nullptr && it->opposite()->vertex()!=nullptr &&
+ it<it->opposite())
+ {
+ export_xfig.segment(CGAL::to_double(it->vertex()->point().x())+dx,
+ CGAL::to_double(it->vertex()->point().y())+dy,
+ CGAL::to_double(it->opposite()->vertex()->point().x())+dx,
+ CGAL::to_double(it->opposite()->vertex()->point().y())+dy,
+ coul_segments, width_segments, false, 90);
+ }
+ }
+ }
+
+ if (radius_disks>0)
+ {
+ for (auto it=vertices().begin(), itend=vertices().end(); it!=itend; ++it)
+ {
+ export_xfig.point(CGAL::to_double(it->point().x())+dx,
+ CGAL::to_double(it->point().y())+dy,
+ radius_disks, coul_disks, false, 70);
+ }
+ }
+ }
+
+protected:
+ CGAL::Compact_container<Halfedge> m_halfedges;
+ CGAL::Compact_container<Vertex> m_vertices;
+};
+////////////////////////////////////////////////////////////////////////////////
+#endif // HDS_H
diff --git a/src/My_arr_segment_traits_2.h b/src/My_arr_segment_traits_2.h
new file mode 100644
index 0000000000000000000000000000000000000000..7fed9307e7fc7c55ed6c9e0a75c9f07b13a839b5
--- /dev/null
+++ b/src/My_arr_segment_traits_2.h
@@ -0,0 +1,1556 @@
+// Copyright (c) 2006,2007,2009,2010,2011 Tel-Aviv University (Israel).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s): Ron Wein <wein@post.tau.ac.il>
+// Efi Fogel <efif@gmail.com>
+// Waqar Khan <wkhan@mpi-inf.mpg.de>
+// Simon Giraudot <simon.giraudot@geometryfactory.com>
+
+#ifndef MY_ARR_SEGMENT_TRAITS_2_H
+#define MY_ARR_SEGMENT_TRAITS_2_H
+
+#include <CGAL/disable_warnings.h>
+
+/*! \file
+ * The segment traits-class for the arrangement package.
+ */
+
+#include <fstream>
+
+#include <boost/variant.hpp>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/tags.h>
+#include <CGAL/intersections.h>
+#include <CGAL/Arr_tags.h>
+#include <CGAL/Arr_enums.h>
+#include <CGAL/Arr_geometry_traits/Segment_assertions.h>
+
+/*! \class A traits class for maintaining an arrangement of segments, avoiding
+ * cascading of computations as much as possible.
+ *
+ * The class is derived from the parameterized kernel to extend the traits
+ * with all the types and operations supported by the kernel. This makes it
+ * possible to use the traits class for data structures that extend the
+ * Arrangement_2 type and require objects and operations supported by the
+ * kernel, but not defined in this derived class.
+ */
+
+template <typename Kernel_ = CGAL::Exact_predicates_exact_constructions_kernel>
+class My_arr_segment_2;
+
+template <typename Kernel_ = CGAL::Exact_predicates_exact_constructions_kernel>
+class My_arr_segment_traits_2 : public Kernel_ {
+ friend class My_arr_segment_2<Kernel_>;
+
+public:
+ typedef Kernel_ Kernel;
+ typedef typename Kernel::FT FT;
+
+ typedef typename CGAL::Algebraic_structure_traits<FT>::Is_exact
+ Has_exact_division;
+
+ // Category tags:
+ typedef CGAL::Tag_true Has_left_category;
+ typedef CGAL::Tag_true Has_merge_category;
+ typedef CGAL::Tag_false Has_do_intersect_category;
+
+ typedef CGAL::Arr_oblivious_side_tag Left_side_category;
+ typedef CGAL::Arr_oblivious_side_tag Bottom_side_category;
+ typedef CGAL::Arr_oblivious_side_tag Top_side_category;
+ typedef CGAL::Arr_oblivious_side_tag Right_side_category;
+
+ typedef typename Kernel::Line_2 Line_2;
+ typedef CGAL::Segment_assertions<My_arr_segment_traits_2<Kernel> >
+ Segment_assertions;
+
+ using Comparison_result=CGAL::Comparison_result;
+#define SMALLER CGAL::SMALLER
+#define EQUAL CGAL::EQUAL
+#define LARGER CGAL::LARGER
+ // using Bbox_2=CGAL::Bbox_2;
+
+ /*! \class Representation of a segment with cached data.
+ */
+ class _Segment_cached_2 {
+ public:
+ typedef typename Kernel::Line_2 Line_2;
+ typedef typename Kernel::Segment_2 Segment_2;
+ typedef typename Kernel::Point_2 Point_2;
+
+ protected:
+ Line_2 m_l; // the line that supports the segment.
+ Point_2 m_ps; // the source point of the segment.
+ Point_2 m_pt; // the target point of the segment.
+ bool m_is_directed_right; // is (lexicographically) directed left to right.
+ bool m_is_vert; // is this a vertical segment.
+ bool m_is_degen; // is the segment degenerate (a single point).
+
+ public:
+
+ /// \name Creation
+ //@{
+
+ /*! Construct default. */
+ _Segment_cached_2();
+
+ /*! Construct a segment from a Kernel segment.
+ * \param seg the segment.
+ * \pre the segment is not degenerate.
+ */
+ _Segment_cached_2(const Segment_2& seg);
+
+ /*! Construct a segment from two endpoints.
+ * \param source the source point.
+ * \param target the target point.
+ * \param `source` and `target` are not equal.
+ */
+ _Segment_cached_2(const Point_2& source, const Point_2& target);
+
+ /*! Construct a segment from two endpoints on a supporting line.
+ * \param line the supporting line.
+ * \param source the source point.
+ * \param target the target point.
+ * \pre `source` and `target` are not equal and both lie on `line`.
+ */
+ _Segment_cached_2(const Line_2& line,
+ const Point_2& source, const Point_2& target);
+
+ /*! Construct a segment from all fields.
+ * \param line the supporting line.
+ * \param source the source point.
+ * \param target the target point.
+ * \param is_directed_right is (lexicographically) directed left to right.
+ * \param is_vert is the segment vertical.
+ * \param is_degen is the segment degenerate (a single point).
+ */
+ _Segment_cached_2(const Line_2& line,
+ const Point_2& source, const Point_2& target,
+ bool is_directed_right, bool is_vert, bool is_degen);
+
+ /*! Assign.
+ * \param seg the source segment to copy from
+ * \pre the segment is not degenerate.
+ */
+ const _Segment_cached_2& operator=(const Segment_2& seg);
+
+ //@}
+
+ /// \name Accessors
+ //@{
+
+ /*! Obtain the supporting line.
+ * \return the supporting line.
+ */
+ const Line_2& line() const;
+
+ /*! Obtain the segment source.
+ * \return the segment source.
+ */
+ const Point_2& source() const;
+
+ /*! Obtain the segment target.
+ * \return the segment target.
+ */
+ const Point_2& target() const;
+
+ /*! Determine whether the curve is vertical.
+ * \return a Boolean flag indicating whether the curve is vertical.
+ */
+ bool is_vertical() const;
+
+ /*! Determine whether the curve is degenerate.
+ * return a Boolean flag indicating whether the curve is degenerate.
+ */
+ bool is_degenerate() const;
+
+ /*! Determine whether the curve is lexicographically directed from left to
+ * right.
+ * \return a Boolean flag indicating whether the curve is lexicographically
+ * directed from left to right.
+ */
+ bool is_directed_right() const;
+
+ /*! Obtain the (lexicographically) left endpoint.
+ * \return the (lexicographically) left endpoint.
+ */
+ const Point_2& left() const;
+
+ /*! Obtain the (lexicographically) right endpoint.
+ * \return the (lexicographically) right endpoint.
+ */
+ const Point_2& right() const;
+
+ //@}
+
+ /// \name Modifiers
+ //@{
+
+ /*! Set the (lexicographically) left endpoint.
+ * \param p the point to set.
+ * \pre p lies on the supporting line to the left of the right endpoint.
+ */
+ void set_left(const Point_2& p);
+
+ /*! Set the (lexicographically) right endpoint.
+ * \param p the point to set.
+ * \pre p lies on the supporting line to the right of the left endpoint.
+ */
+ void set_right(const Point_2& p);
+
+ //@}
+
+ /// \name Deprecated
+ //@{
+
+ /*! Determine whether the given point is in the x-range of the segment.
+ * \param p the query point.
+ * \return (true) is in the x-range of the segment; (false) if it is not.
+ */
+ CGAL_DEPRECATED bool is_in_x_range(const Point_2& p) const;
+
+ /*! Determine whether the given point is in the y-range of the segment.
+ * \param p the query point.
+ * \return (true) is in the y-range of the segment; (false) if it is not.
+ */
+ CGAL_DEPRECATED bool is_in_y_range(const Point_2& p) const;
+
+ //@}
+ };
+
+public:
+ // Traits objects
+ typedef typename Kernel::Point_2 Point_2;
+ typedef My_arr_segment_2<Kernel> X_monotone_curve_2;
+ typedef My_arr_segment_2<Kernel> Curve_2;
+ typedef unsigned int Multiplicity;
+
+public:
+ /*! Construct default. */
+ My_arr_segment_traits_2() {}
+
+ /// \name Basic functor definitions.
+ //@{
+
+ class Compare_x_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ //! The traits (in case it has state).
+ const Traits& m_traits;
+
+ /*! Constructor
+ * \param traits the traits (in case it has state)
+ */
+ Compare_x_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ public:
+ /*! Compare the x-coordinates of two points.
+ * \param p1 the first point.
+ * \param p2 the second point.
+ * \return LARGER if x(p1) > x(p2);
+ * SMALLER if x(p1) < x(p2);
+ * EQUAL if x(p1) = x(p2).
+ */
+ Comparison_result operator()(const Point_2& p1, const Point_2& p2) const
+ {
+ const Kernel& kernel = m_traits;
+ return (kernel.compare_x_2_object()(p1, p2));
+ }
+ };
+
+ /*! Obtain a Compare_x_2 functor object. */
+ Compare_x_2 compare_x_2_object() const { return Compare_x_2(*this); }
+
+ class Compare_xy_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ /*! The traits (in case it has state) */
+ const Traits& m_traits;
+
+ /*! Constructor
+ * \param traits the traits (in case it has state)
+ */
+ Compare_xy_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ public:
+ /*! Compare two points lexicographically: by x, then by y.
+ * \param p1 the first point.
+ * \param p2 the second point.
+ * \return LARGER if x(p1) > x(p2), or if x(p1) = x(p2) and y(p1) > y(p2);
+ * SMALLER if x(p1) < x(p2), or if x(p1) = x(p2) and y(p1) < y(p2);
+ * EQUAL if the two points are equal.
+ */
+ Comparison_result operator()(const Point_2& p1, const Point_2& p2) const
+ {
+ const Kernel& kernel = m_traits;
+ return (kernel.compare_xy_2_object()(p1, p2));
+ }
+ };
+
+ /*! Obtain a Compare_xy_2 functor object. */
+ Compare_xy_2 compare_xy_2_object() const { return Compare_xy_2(*this); }
+
+ class Construct_min_vertex_2 {
+ public:
+ /*! Obtain the left endpoint of the x-monotone curve (segment).
+ * \param cv the curve.
+ * \return the left endpoint.
+ */
+ const Point_2& operator()(const X_monotone_curve_2& cv) const
+ { return (cv.left()); }
+ };
+
+ /*! Obtain a Construct_min_vertex_2 functor object. */
+ Construct_min_vertex_2 construct_min_vertex_2_object() const
+ { return Construct_min_vertex_2(); }
+
+ class Construct_max_vertex_2 {
+ public:
+ /*! Obtain the right endpoint of the x-monotone curve (segment).
+ * \param cv the curve.
+ * \return the right endpoint.
+ */
+ const Point_2& operator()(const X_monotone_curve_2& cv) const
+ { return (cv.right()); }
+ };
+
+ /*! Obtain a Construct_max_vertex_2 functor object. */
+ Construct_max_vertex_2 construct_max_vertex_2_object() const
+ { return Construct_max_vertex_2(); }
+
+ class Is_vertical_2 {
+ public:
+ /*! Check whether the given x-monotone curve is a vertical segment.
+ * \param cv the curve.
+ * \return (true) if the curve is a vertical segment; (false) otherwise.
+ */
+ bool operator()(const X_monotone_curve_2& cv) const
+ { return (cv.is_vertical()); }
+ };
+
+ /*! Obtain an Is_vertical_2 functor object. */
+ Is_vertical_2 is_vertical_2_object () const { return Is_vertical_2(); }
+
+ class Compare_y_at_x_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ /*! the traits (in case it has state) */
+ const Traits& m_traits;
+
+ /*! Constructor
+ * \param traits the traits (in case it has state)
+ */
+ Compare_y_at_x_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ public:
+ /*! Return the location of the given point with respect to the input curve.
+ * \param cv the curve.
+ * \param p the point.
+ * \pre p is in the x-range of cv.
+ * \return SMALLER if y(p) < cv(x(p)), i.e. the point is below the curve;
+ * LARGER if y(p) > cv(x(p)), i.e. the point is above the curve;
+ * EQUAL if p lies on the curve.
+ */
+ Comparison_result operator()(const Point_2& p,
+ const X_monotone_curve_2& cv) const
+ {
+ CGAL_precondition(m_traits.is_in_x_range_2_object()(cv, p));
+
+ const Kernel& kernel = m_traits;
+
+ if (! cv.is_vertical()) {
+ // Compare p with the segment supporting line.
+ CGAL_assertion_code(auto cmp_x = kernel.compare_x_2_object());
+ CGAL_assertion(cmp_x(cv.left(), cv.right()) == SMALLER);
+ return kernel.orientation_2_object()(cv.left(), cv.right(), p);
+ }
+
+ // Compare with the vertical segment endpoints.
+ auto compare_y = kernel.compare_y_2_object();
+ Comparison_result res1 = compare_y(p, cv.left());
+ Comparison_result res2 = compare_y(p, cv.right());
+ return (res1 == res2) ? res1 : EQUAL;
+ }
+ };
+
+ /*! Obtain a Compare_y_at_x_2 functor object. */
+ Compare_y_at_x_2 compare_y_at_x_2_object() const
+ { return Compare_y_at_x_2(*this); }
+
+ class Compare_y_at_x_left_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ /*! The traits (in case it has state) */
+ const Traits& m_traits;
+
+ /*! Constructor
+ * \param traits the traits (in case it has state)
+ */
+ Compare_y_at_x_left_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ public:
+ /*! Compare the y value of two x-monotone curves immediately to the left
+ * of their intersection point.
+ * \param cv1 the first curve.
+ * \param cv2 the second curve.
+ * \param p the intersection point.
+ * \pre the point p lies on both curves, and both of them must be also be
+ * defined (lexicographically) to its left.
+ * \return the relative position of cv1 with respect to cv2 immediately to
+ * the left of p: SMALLER, LARGER or EQUAL.
+ */
+ Comparison_result operator()(const X_monotone_curve_2& cv1,
+ const X_monotone_curve_2& cv2,
+ const Point_2& CGAL_assertion_code(p)) const
+ {
+ const Kernel& kernel = m_traits;
+
+ // Make sure that p lies on both curves, and that both are defined to its
+ // left (so their left endpoint is lexicographically smaller than p).
+ CGAL_precondition_code(auto compare_xy = kernel.compare_xy_2_object());
+
+ CGAL_precondition((m_traits.compare_y_at_x_2_object()(p, cv1) == EQUAL) &&
+ (m_traits.compare_y_at_x_2_object()(p, cv2) == EQUAL));
+
+ CGAL_precondition(compare_xy(cv1.left(), p) == SMALLER &&
+ compare_xy(cv2.left(), p) == SMALLER);
+
+ // Compare the slopes of the two segments to determine their relative
+ // position immediately to the left of q.
+ // Notice we use the supporting lines in order to compare the slopes,
+ // and that we swap the order of the curves in order to obtain the
+ // correct result to the left of p.
+ return (kernel.compare_slope_2_object()(cv2.line(), cv1.line()));
+ }
+ };
+
+ /*! Obtain a Compare_y_at_x_left_2 functor object. */
+ Compare_y_at_x_left_2 compare_y_at_x_left_2_object() const
+ { return Compare_y_at_x_left_2(*this); }
+
+ class Compare_y_at_x_right_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ /*! The traits (in case it has state) */
+ const Traits& m_traits;
+
+ /*! Constructor
+ * \param traits the traits (in case it has state)
+ */
+ Compare_y_at_x_right_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ public:
+ /*! Compare the y value of two x-monotone curves immediately to the right
+ * of their intersection point.
+ * \param cv1 the first curve.
+ * \param cv2 the second curve.
+ * \param p the intersection point.
+ * \pre the point p lies on both curves, and both of them must be also be
+ * defined (lexicographically) to its right.
+ * \return the relative position of cv1 with respect to cv2 immediately to
+ * the right of p: SMALLER, LARGER or EQUAL.
+ */
+ Comparison_result operator()(const X_monotone_curve_2& cv1,
+ const X_monotone_curve_2& cv2,
+ const Point_2& CGAL_assertion_code(p)) const
+ {
+ const Kernel& kernel = m_traits;
+
+ // Make sure that p lies on both curves, and that both are defined to its
+ // right (so their right endpoint is lexicographically larger than p).
+ CGAL_precondition_code(auto compare_xy = kernel.compare_xy_2_object());
+
+ CGAL_precondition((m_traits.compare_y_at_x_2_object()(p, cv1) == EQUAL) &&
+ (m_traits.compare_y_at_x_2_object()(p, cv2) == EQUAL));
+
+ CGAL_precondition(compare_xy(cv1.right(), p) == LARGER &&
+ compare_xy(cv2.right(), p) == LARGER);
+
+ // Compare the slopes of the two segments to determine their relative
+ // position immediately to the left of q.
+ // Notice we use the supporting lines in order to compare the slopes.
+ return (kernel.compare_slope_2_object()(cv1.line(), cv2.line()));
+ }
+ };
+
+ /*! Obtain a Compare_y_at_x_right_2 functor object. */
+ Compare_y_at_x_right_2 compare_y_at_x_right_2_object() const
+ { return Compare_y_at_x_right_2(*this); }
+
+ class Equal_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ /*! The traits (in case it has state) */
+ const Traits& m_traits;
+
+ /*! Constructor
+ * \param traits the traits (in case it has state)
+ */
+ Equal_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ public:
+ /*! Check whether the two x-monotone curves are the same (have the same
+ * graph).
+ * \param cv1 the first curve.
+ * \param cv2 the second curve.
+ * \return (true) if the two curves are the same; (false) otherwise.
+ */
+ bool operator()(const X_monotone_curve_2& cv1,
+ const X_monotone_curve_2& cv2) const
+ {
+ const Kernel& kernel = m_traits;
+ typename Kernel::Equal_2 equal = kernel.equal_2_object();
+
+ return (equal(cv1.left(), cv2.left()) &&
+ equal(cv1.right(), cv2.right()));
+ }
+
+ /*! Determine whether the two points are the same.
+ * \param p1 the first point.
+ * \param p2 the second point.
+ * \return (true) if the two point are the same; (false) otherwise.
+ */
+ bool operator()(const Point_2& p1, const Point_2& p2) const
+ {
+ const Kernel& kernel = m_traits;
+ return (kernel.equal_2_object()(p1, p2));
+ }
+ };
+
+ /*! Obtain an Equal_2 functor object. */
+ Equal_2 equal_2_object() const { return Equal_2(*this); }
+
+ //@}
+
+ //! \name Intersections, subdivisions, and mergings
+ //@{
+
+ /*! \class Make_x_monotone_2
+ * A functor for subdividing a curve into x-monotone curves.
+ */
+ class Make_x_monotone_2 {
+ public:
+ /*! Subdivide a given curve into x-monotone subcurves and insert them into
+ * a given output iterator. As segments are always x_monotone a single
+ * object is inserted.
+ * \param cv the curve.
+ * \param oi the output iterator for the result. Its dereference type is a
+ * variant that wraps a \c Point_2 or an \c X_monotone_curve_2
+ * objects.
+ * \return the past-the-end output iterator.
+ */
+ template <typename OutputIterator>
+ OutputIterator operator()(const Curve_2& cv, OutputIterator oi) const
+ {
+ // Wrap the segment with a variant.
+ typedef boost::variant<Point_2, X_monotone_curve_2>
+ Make_x_monotone_result;
+ *oi++ = Make_x_monotone_result(cv);
+ return oi;
+ }
+ };
+
+ /*! Obtain a Make_x_monotone_2 functor object. */
+ Make_x_monotone_2 make_x_monotone_2_object() const
+ { return Make_x_monotone_2(); }
+
+ class Split_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ /*! The traits (in case it has state) */
+ const Traits& m_traits;
+
+ /*! Constructor
+ * \param traits the traits (in case it has state)
+ */
+ Split_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ public:
+ /*! Split a given x-monotone curve at a given point into two sub-curves.
+ * \param cv the curve to split
+ * \param p the split point.
+ * \param c1 Output: the left resulting subcurve (p is its right endpoint).
+ * \param c2 Output: the right resulting subcurve (p is its left endpoint).
+ * \pre p lies on cv but is not one of its endpoints.
+ */
+ void operator()(const X_monotone_curve_2& cv, const Point_2& p,
+ X_monotone_curve_2& c1, X_monotone_curve_2& c2) const
+ {
+ // Make sure that p lies on the interior of the curve.
+ CGAL_precondition_code(const Kernel& kernel = m_traits;
+ auto compare_xy = kernel.compare_xy_2_object());
+
+ CGAL_precondition((m_traits.compare_y_at_x_2_object()(p, cv) == EQUAL) &&
+ compare_xy(cv.left(), p) == SMALLER &&
+ compare_xy(cv.right(), p) == LARGER);
+
+ // Perform the split.
+ c1 = cv;
+ c1.set_right(p);
+
+ c2 = cv;
+ c2.set_left(p);
+ }
+ };
+
+ /*! Obtain a Split_2 functor object. */
+ Split_2 split_2_object() const { return Split_2(*this); }
+
+ class Intersect_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ /*! The traits (in case it has state) */
+ const Traits& m_traits;
+
+ /*! Construct
+ * \param traits the traits (in case it has state)
+ */
+ Intersect_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ // Specialized do_intersect with many tests skipped because at
+ // this point, we already know which point is left / right for
+ // both segments
+ bool do_intersect(const Point_2& A1, const Point_2& A2,
+ const Point_2& B1, const Point_2& B2) const
+ {
+ const Kernel& kernel = m_traits;
+ auto compare_xy = kernel.compare_xy_2_object();
+ namespace interx = CGAL::Intersections::internal;
+
+ switch(make_certain(compare_xy(A1,B1))) {
+ case SMALLER:
+ switch(make_certain(compare_xy(A2,B1))) {
+ case SMALLER: return false;
+ case EQUAL: return true;
+ default: // LARGER
+ switch(make_certain(compare_xy(A2,B2))) {
+ case SMALLER:
+ return interx::seg_seg_do_intersect_crossing(A1,A2,B1,B2, kernel);
+ case EQUAL: return true;
+ default: // LARGER
+ return interx::seg_seg_do_intersect_contained(A1,A2,B1,B2, kernel);
+ }
+ }
+ case EQUAL: return true;
+ default: // LARGER
+ switch(make_certain(compare_xy(B2,A1))) {
+ case SMALLER: return false;
+ case EQUAL: return true;
+ default: // LARGER
+ switch(make_certain(compare_xy(B2,A2))) {
+ case SMALLER:
+ return interx::seg_seg_do_intersect_crossing(B1,B2,A1,A2, kernel);
+ case EQUAL: return true;
+ default: // LARGER
+ return interx::seg_seg_do_intersect_contained(B1,B2,A1,A2, kernel);
+ }
+ }
+ }
+ CGAL_error(); // never reached
+ return false;
+ }
+
+ /*! Determine whether the bounding boxes of two segments overlap
+ */
+ bool do_bboxes_overlap(const X_monotone_curve_2& cv1,
+ const X_monotone_curve_2& cv2) const
+ {
+ const Kernel& kernel = m_traits;
+ auto construct_bbox = kernel.construct_bbox_2_object();
+ auto bbox1 = construct_bbox(cv1.source()) + construct_bbox(cv1.target());
+ auto bbox2 = construct_bbox(cv2.source()) + construct_bbox(cv2.target());
+ return CGAL::do_overlap(bbox1, bbox2);
+ }
+
+ public:
+ /*! Find the intersections of the two given curves and insert them into the
+ * given output iterator. As two segments may intersect only once, only a
+ * single intersection will be contained in the iterator.
+ * \param cv1 the first curve.
+ * \param cv2 the second curve.
+ * \param oi the output iterator.
+ * \return the past-the-end iterator.
+ */
+ template <typename OutputIterator>
+ OutputIterator operator()(const X_monotone_curve_2& cv1,
+ const X_monotone_curve_2& cv2,
+ OutputIterator oi) const
+ {
+ typedef std::pair<Point_2, Multiplicity> Intersection_point;
+ typedef boost::variant<Intersection_point, X_monotone_curve_2>
+ Intersection_result;
+
+ // Early ending with Bbox overlapping test
+ if (! do_bboxes_overlap(cv1, cv2)) return oi;
+
+ // Early ending with specialized do_intersect
+ const Kernel& kernel = m_traits;
+ if (! do_intersect(cv1.left(), cv1.right(), cv2.left(), cv2.right()))
+ return oi;
+
+ // An intersection is guaranteed.
+
+ // Intersect the two supporting lines.
+ auto res = kernel.intersect_2_object()(cv1.line(), cv2.line());
+ CGAL_assertion(bool(res));
+
+ // Check if we have a single intersection point.
+ const Point_2* ip = boost::get<Point_2>(&*res);
+ if (ip != nullptr) {
+ CGAL_assertion(cv1.is_vertical() ?
+ m_traits.is_in_y_range_2_object()(cv1, *ip) :
+ m_traits.is_in_x_range_2_object()(cv1, *ip));
+ CGAL_assertion(cv2.is_vertical() ?
+ m_traits.is_in_y_range_2_object()(cv2, *ip) :
+ m_traits.is_in_x_range_2_object()(cv2, *ip));
+ Intersection_point ip_mult(*ip, 1);
+ *oi++ = Intersection_result(ip_mult);
+ return oi;
+ }
+
+ // In this case, the two supporting lines overlap.
+ // The overlapping segment is therefore [p_l,p_r], where p_l is the
+ // rightmost of the two left endpoints and p_r is the leftmost of the
+ // two right endpoints.
+ auto compare_xy = kernel.compare_xy_2_object();
+ const Point_2& p_l = (compare_xy(cv1.left(), cv2.left()) == SMALLER) ?
+ cv2.left() : cv1.left();
+ const Point_2& p_r = (compare_xy(cv1.right(), cv2.right()) == SMALLER) ?
+ cv1.right() : cv2.right();
+
+ // Examine the resulting segment.
+ const Comparison_result cmp_res = compare_xy(p_l, p_r);
+ if (cmp_res == EQUAL) {
+ // The two segment have the same supporting line, but they just share
+ // a common endpoint. Thus we have an intersection point, but we leave
+ // the multiplicity of this point undefined.
+ Intersection_point ip_mult(p_r, 0);
+ *oi++ = Intersection_result(ip_mult);
+ return oi;
+ }
+
+ CGAL_assertion(cmp_res == SMALLER);
+ // We have discovered an overlapping segment:
+ if (cv1.is_directed_right() == cv2.is_directed_right()) {
+ // cv1 and cv2 have the same directions, maintain this direction
+ // in the overlap segment
+ if (cv1.is_directed_right()) {
+ X_monotone_curve_2 overlap_seg(cv1.line(), p_l, p_r, cv1.m_original_segment_id);
+ *oi++ = Intersection_result(overlap_seg);
+ return oi;
+ }
+ X_monotone_curve_2 overlap_seg(cv1.line(), p_r, p_l, cv1.m_original_segment_id);
+ *oi++ = Intersection_result(overlap_seg);
+ return oi;
+ }
+ // cv1 and cv2 have opposite directions, the overlap segment
+ // will be directed from left to right
+ X_monotone_curve_2 overlap_seg(cv1.line(), p_l, p_r, cv1.m_original_segment_id);
+ *oi++ = Intersection_result(overlap_seg);
+ return oi;
+ }
+ };
+
+ /*! Obtain an Intersect_2 functor object. */
+ Intersect_2 intersect_2_object() const { return Intersect_2(*this); }
+
+ class Are_mergeable_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ /*! The traits (in case it has state) */
+ const Traits& m_traits;
+
+ /*! Constructor
+ * \param traits the traits (in case it has state)
+ */
+ Are_mergeable_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ public:
+ /*! Check whether it is possible to merge two given x-monotone curves.
+ * \param cv1 the first curve.
+ * \param cv2 the second curve.
+ * \return (true) if the two curves are mergeable, that is, if they are
+ * supported by the same line; (false) otherwise.
+ * \pre cv1 and cv2 share a common endpoint.
+ */
+ bool operator()(const X_monotone_curve_2& cv1,
+ const X_monotone_curve_2& cv2) const
+ {
+ const Kernel& kernel = m_traits;
+ typename Kernel::Equal_2 equal = kernel.equal_2_object();
+ if (! equal(cv1.right(), cv2.left()) &&
+ ! equal(cv2.right(), cv1.left()))
+ return false;
+
+ // Check whether the two curves have the same supporting line.
+ return (equal(cv1.line(), cv2.line()) ||
+ equal(cv1.line(),
+ kernel.construct_opposite_line_2_object()(cv2.line())));
+ }
+ };
+
+ /*! Obtain an Are_mergeable_2 functor object. */
+ Are_mergeable_2 are_mergeable_2_object() const
+ { return Are_mergeable_2(*this); }
+
+ /*! \class Merge_2
+ * A functor that merges two x-monotone arcs into one.
+ */
+ class Merge_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ /*! The traits (in case it has state) */
+ const Traits& m_traits;
+
+ /*! Constructor
+ * \param traits the traits (in case it has state)
+ */
+ Merge_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ public:
+ /*! Merge two given x-monotone curves into a single curve (segment).
+ * \param cv1 the first curve.
+ * \param cv2 the second curve.
+ * \param c Output: the merged curve.
+ * \pre the two curves are mergeable.
+ */
+ void operator()(const X_monotone_curve_2& cv1,
+ const X_monotone_curve_2& cv2,
+ X_monotone_curve_2& c) const
+ {
+ CGAL_precondition(m_traits.are_mergeable_2_object()(cv1, cv2));
+
+ const Kernel& kernel = m_traits;
+ auto equal = kernel.equal_2_object();
+
+ // Check which curve extends to the right of the other.
+ if (equal(cv1.right(), cv2.left())) {
+ // cv2 extends cv1 to the right.
+ c = cv1;
+ c.set_right(cv2.right());
+ return;
+ }
+
+ CGAL_precondition(equal(cv2.right(), cv1.left()));
+
+ // cv1 extends cv2 to the right.
+ c = cv2;
+ c.set_right(cv1.right());
+ }
+ };
+
+ /*! Obtain a Merge_2 functor object. */
+ Merge_2 merge_2_object() const { return Merge_2(*this); }
+ //@}
+
+ /// \name Functor definitions for the landmarks point-location strategy.
+ //@{
+ typedef double Approximate_number_type;
+
+ class Approximate_2 {
+ public:
+ /*! Obtain an approximation of a point coordinate.
+ * \param p the exact point.
+ * \param i the coordinate index (either 0 or 1).
+ * \pre i is either 0 or 1.
+ * \return An approximation of p's x-coordinate (if i == 0), or an
+ * approximation of p's y-coordinate (if i == 1).
+ */
+ Approximate_number_type operator()(const Point_2& p, int i) const
+ {
+ CGAL_precondition((i == 0) || (i == 1));
+ return (i == 0) ? (CGAL::to_double(p.x())) : (CGAL::to_double(p.y()));
+ }
+ };
+
+ /*! Obtain an Approximate_2 functor object. */
+ Approximate_2 approximate_2_object() const { return Approximate_2(); }
+
+ class Construct_x_monotone_curve_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ //! The traits (in case it has state).
+ const Traits& m_traits;
+
+ /*! Constructor
+ * \param traits the traits (in case it has state)
+ */
+ Construct_x_monotone_curve_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ public:
+ typedef typename Kernel::Segment_2 Segment_2;
+
+ /*! Obtain an x-monotone curve connecting two given endpoints.
+ * \param source the first point.
+ * \param target the second point.
+ * \pre `source` and `target` must not be equal.
+ * \return a segment connecting `source` and `target`.
+ */
+ X_monotone_curve_2 operator()(const Point_2& source,
+ const Point_2& target) const
+ {
+ const Kernel& kernel = m_traits;
+ auto line = kernel.construct_line_2_object()(source, target);
+ Comparison_result res = kernel.compare_xy_2_object()(source, target);
+ auto is_degen = (res == EQUAL);
+ auto is_directed_right = (res == SMALLER);
+ CGAL_precondition_msg(! is_degen,
+ "Cannot construct a degenerate segment.");
+ auto is_vert = kernel.is_vertical_2_object()(line);
+ return X_monotone_curve_2(line, source, target,
+ is_directed_right, is_vert, is_degen);
+ }
+
+ /*! Obtain an \f$x\f$-monotone curve given a Kernel segment.
+ * \param seg the segment.
+ * \return the \f$x\f$-monotone curve.
+ * \pre the segment is not degenerate.
+ * \return a segment that is the same as `seg`..
+ */
+ X_monotone_curve_2 operator()(const Segment_2& seg) const
+ {
+ const Kernel& kernel = m_traits;
+ auto line = kernel.construct_line_2_object()(seg);
+ auto vertex_ctr = kernel.construct_vertex_2_object();
+ auto source = vertex_ctr(seg, 0);
+ auto target = vertex_ctr(seg, 1);
+ Comparison_result res = kernel.compare_xy_2_object()(source, target);
+ auto is_degen = (res == EQUAL);
+ auto is_directed_right = (res == SMALLER);
+ CGAL_precondition_msg(! is_degen,
+ "Cannot construct a degenerate segment.");
+ auto is_vert = kernel.is_vertical_2_object()(seg);
+ return X_monotone_curve_2(line, source, target,
+ is_directed_right, is_vert, is_degen);
+ }
+
+ /*! Obtain an \f$x\f$-monotone curve given two endpoints and the supporting
+ * line.
+ * \param line the supporting line.
+ * \param the source point.
+ * \param the target point.
+ * \pre `ps` and `pt` are not equal and both lie on `line`.
+ */
+ X_monotone_curve_2 operator()(const Line_2& line,
+ const Point_2& source,
+ const Point_2& target) const
+ {
+ const Kernel& kernel = m_traits;
+ CGAL_precondition
+ (Segment_assertions::_assert_is_point_on(source, line,
+ Has_exact_division()) &&
+ Segment_assertions::_assert_is_point_on(target, line,
+ Has_exact_division()));
+ auto is_vert = kernel.is_vertical_2_object()(line);
+ Comparison_result res = kernel.compare_xy_2_object()(source, target);
+ auto is_degen = (res == EQUAL);
+ auto is_directed_right = (res == SMALLER);
+ CGAL_precondition_msg(! is_degen,
+ "Cannot construct a degenerate segment.");
+ return X_monotone_curve_2(line, source, target,
+ is_directed_right, is_vert, is_degen);
+ }
+ };
+
+ /*! Obtain a Construct_x_monotone_curve_2 functor object. */
+ Construct_x_monotone_curve_2 construct_x_monotone_curve_2_object() const
+ { return Construct_x_monotone_curve_2(*this); }
+ //@}
+
+ /// \name Functor definitions for the Boolean set-operation traits.
+ //@{
+
+ class Trim_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ /*! The traits (in case it has state). */
+ const Traits& m_traits;
+
+ /*! Constructor
+ * \param traits the traits (in case it has state)
+ */
+ Trim_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ /*! Obtain a trimmed version of a line.
+ * \param xseg the x-monotone segment.
+ * \param src the new start endpoint.
+ * \param tgt the new end endpoint.
+ * \return the trimmed x-monotone segment.
+ * \pre src != tgt
+ * \pre both points must lie on segment
+ */
+ public:
+ X_monotone_curve_2 operator()(const X_monotone_curve_2& xcv,
+ const Point_2& src,
+ const Point_2& tgt)const
+ {
+ CGAL_precondition_code(Equal_2 equal = m_traits.equal_2_object());
+ CGAL_precondition_code(Compare_y_at_x_2 compare_y_at_x =
+ m_traits.compare_y_at_x_2_object());
+ Compare_x_2 compare_x_2 = m_traits.compare_x_2_object();
+
+ // check whether source and taget are two distinct points and they lie
+ // on the line.
+ CGAL_precondition(!equal(src, tgt));
+ CGAL_precondition(compare_y_at_x(src, xcv) == EQUAL);
+ CGAL_precondition(compare_y_at_x(tgt, xcv) == EQUAL);
+
+ // exchange src and tgt IF they do not conform with the direction
+ X_monotone_curve_2 trimmed_segment;
+ if (xcv.is_directed_right() && compare_x_2(src, tgt) == LARGER)
+ trimmed_segment = X_monotone_curve_2(tgt, src);
+ else if (! xcv.is_directed_right() && (compare_x_2(src, tgt) == SMALLER))
+ trimmed_segment = X_monotone_curve_2(tgt, src);
+ else trimmed_segment = X_monotone_curve_2(src, tgt);
+ return trimmed_segment;
+ }
+ };
+
+ /*! Obtain a Trim_2 functor object */
+ Trim_2 trim_2_object() const { return Trim_2(*this); }
+
+ class Compare_endpoints_xy_2 {
+ public:
+ /*! Compare the endpoints of an $x$-monotone curve lexicographically.
+ * (assuming the curve has a designated source and target points).
+ * \param cv the curve.
+ * \return SMALLER if the curve is directed right;
+ * LARGER if the curve is directed left.
+ */
+ Comparison_result operator()(const X_monotone_curve_2& cv) const
+ { return (cv.is_directed_right()) ? (SMALLER) : (LARGER); }
+ };
+
+ /*! Obtain a Compare_endpoints_xy_2 functor object. */
+ Compare_endpoints_xy_2 compare_endpoints_xy_2_object() const
+ { return Compare_endpoints_xy_2(); }
+
+ class Construct_opposite_2 {
+ public:
+ /*! Construct an opposite x-monotone (with swapped source and target).
+ * \param cv the curve.
+ * \return the opposite curve.
+ */
+ X_monotone_curve_2 operator()(const X_monotone_curve_2& cv) const
+ { return (cv.flip()); }
+ };
+
+ /*! Obtain a Construct_opposite_2 functor object. */
+ Construct_opposite_2 construct_opposite_2_object() const
+ { return Construct_opposite_2(); }
+ //@}
+
+ /// \name Utilities.
+ //@{
+
+ class Is_in_x_range_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ //! The traits (in case it has state).
+ const Traits& m_traits;
+
+ /*! Construct
+ * \param traits the traits (in case it has state)
+ */
+ Is_in_x_range_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ public:
+ /*! Determine whether a given point is in the \f$x\f$-range of a given
+ * segment.
+ * \param cv the segment.
+ * \param p the point.
+ * \return true if p is in the \f$x\f$-range of cv; false otherwise.
+ */
+ bool operator()(const X_monotone_curve_2& cv, const Point_2& p) const
+ {
+ const Kernel& kernel = m_traits;
+ auto compare_x = kernel.compare_x_2_object();
+ Comparison_result res1 = compare_x(p, cv.left());
+
+ if (res1 == SMALLER) return false;
+ else if (res1 == EQUAL) return true;
+
+ Comparison_result res2 = compare_x(p, cv.right());
+ return (res2 != LARGER);
+ }
+ };
+
+ /*! Obtain an Is_in_x_range_2 functor object */
+ Is_in_x_range_2 is_in_x_range_2_object() const
+ { return Is_in_x_range_2(*this); }
+
+ class Is_in_y_range_2 {
+ protected:
+ typedef My_arr_segment_traits_2<Kernel> Traits;
+
+ //! The traits (in case it has state).
+ const Traits& m_traits;
+
+ /*! Construct
+ * \param traits the traits (in case it has state)
+ */
+ Is_in_y_range_2(const Traits& traits) : m_traits(traits) {}
+
+ friend class My_arr_segment_traits_2<Kernel>;
+
+ public:
+ /*! Determine whether a given point is in the \f$y\f$-range of a given
+ * segment.
+ * \param cv the segment.
+ * \param p the point.
+ * \return true if p is in the \f$y\f$-range of cv; false otherwise.
+ */
+ bool operator()(const X_monotone_curve_2& cv, const Point_2& p) const
+ {
+ const Kernel& kernel = m_traits;
+ auto compare_y = kernel.compare_y_2_object();
+ Comparison_result res1 = compare_y(p, cv.left());
+
+ if (res1 == SMALLER) return false;
+ else if (res1 == EQUAL) return true;
+
+ Comparison_result res2 = compare_y(p, cv.right());
+ return (res2 != LARGER);
+ }
+ };
+
+ /*! Obtain an Is_in_y_range_2 functor object */
+ Is_in_y_range_2 is_in_y_range_2_object() const
+ { return Is_in_y_range_2(*this); }
+
+ //@}
+};
+
+// Creation
+
+//! \brief constructs default.
+template <typename Kernel>
+My_arr_segment_traits_2<Kernel>::_Segment_cached_2::_Segment_cached_2() :
+ m_is_directed_right(false),
+ m_is_vert(false),
+ m_is_degen(true)
+{}
+
+//! \brief constructs a segment from a Kernel segment.
+template <typename Kernel>
+My_arr_segment_traits_2<Kernel>::
+_Segment_cached_2::_Segment_cached_2(const Segment_2& seg)
+{
+ Kernel kernel;
+ auto vertex_ctr = kernel.construct_vertex_2_object();
+
+ m_ps = vertex_ctr(seg, 0);
+ m_pt = vertex_ctr(seg, 1);
+
+ Comparison_result res = kernel.compare_xy_2_object()(m_ps, m_pt);
+ m_is_degen = (res == EQUAL);
+ m_is_directed_right = (res == SMALLER);
+
+ CGAL_precondition_msg(! m_is_degen, "Cannot construct a degenerate segment.");
+
+ m_l = kernel.construct_line_2_object()(seg);
+ m_is_vert = kernel.is_vertical_2_object()(seg);
+}
+
+//! \brief Constructs a segment from two endpoints.
+template <typename Kernel>
+My_arr_segment_traits_2<Kernel>::
+_Segment_cached_2::_Segment_cached_2(const Point_2& source,
+ const Point_2& target) :
+ m_ps(source),
+ m_pt(target)
+{
+ Kernel kernel;
+
+ Comparison_result res = kernel.compare_xy_2_object()(m_ps, m_pt);
+ m_is_degen = (res == EQUAL);
+ m_is_directed_right = (res == SMALLER);
+
+ CGAL_precondition_msg(! m_is_degen, "Cannot construct a degenerate segment.");
+
+ m_l = kernel.construct_line_2_object()(source, target);
+ m_is_vert = kernel.is_vertical_2_object()(m_l);
+}
+
+//! \brief constructs a segment from two endpoints on a supporting line.
+template <typename Kernel>
+My_arr_segment_traits_2<Kernel>::
+_Segment_cached_2::_Segment_cached_2(const Line_2& line,
+ const Point_2& source,
+ const Point_2& target) :
+ m_l(line),
+ m_ps(source),
+ m_pt(target)
+{
+ Kernel kernel;
+
+ CGAL_precondition
+ (Segment_assertions::_assert_is_point_on(source, m_l,
+ Has_exact_division()) &&
+ Segment_assertions::_assert_is_point_on(target, m_l,
+ Has_exact_division()));
+
+ m_is_vert = kernel.is_vertical_2_object()(m_l);
+
+ Comparison_result res = kernel.compare_xy_2_object()(m_ps, m_pt);
+ m_is_degen = (res == EQUAL);
+ m_is_directed_right = (res == SMALLER);
+
+ CGAL_precondition_msg(! m_is_degen, "Cannot construct a degenerate segment.");
+}
+
+//! \brief constructs a segment from all fields.
+template <typename Kernel>
+My_arr_segment_traits_2<Kernel>::_Segment_cached_2::
+_Segment_cached_2(const Line_2& line,
+ const Point_2& source, const Point_2& target,
+ bool is_directed_right, bool is_vert, bool is_degen) :
+ m_l(line),
+ m_ps(source),
+ m_pt(target),
+ m_is_directed_right(is_directed_right),
+ m_is_vert(is_vert),
+ m_is_degen(is_degen)
+{}
+
+//! \brief assigns.
+template <typename Kernel>
+const typename My_arr_segment_traits_2<Kernel>::_Segment_cached_2&
+My_arr_segment_traits_2<Kernel>::_Segment_cached_2::operator=(const Segment_2& seg)
+{
+ Kernel kernel;
+ auto vertex_ctr = kernel.construct_vertex_2_object();
+
+ m_ps = vertex_ctr(seg, 0);
+ m_pt = vertex_ctr(seg, 1);
+
+ Comparison_result res = kernel.compare_xy_2_object()(m_ps, m_pt);
+ m_is_degen = (res == EQUAL);
+ m_is_directed_right = (res == SMALLER);
+
+ CGAL_precondition_msg(! m_is_degen, "Cannot construct a degenerate segment.");
+
+ m_l = kernel.construct_line_2_object()(seg);
+ m_is_vert = kernel.is_vertical_2_object()(seg);
+
+ return (*this);
+}
+
+// Accessors
+
+//! \brief obtains the supporting line.
+template <typename Kernel>
+const typename Kernel::Line_2&
+My_arr_segment_traits_2<Kernel>::_Segment_cached_2::line() const { return m_l; }
+
+//! \brief determines whether the curve is vertical.
+template <typename Kernel>
+bool My_arr_segment_traits_2<Kernel>::_Segment_cached_2::is_vertical() const
+{ return m_is_vert; }
+
+//! \brief determines whether the curve is degenerate.
+template <typename Kernel>
+bool My_arr_segment_traits_2<Kernel>::_Segment_cached_2::is_degenerate() const
+{ return m_is_degen; }
+
+/*! \brief determines whether the curve is lexicographically directed from
+ * left to right
+ */
+template <typename Kernel>
+bool My_arr_segment_traits_2<Kernel>::_Segment_cached_2::is_directed_right() const
+{ return m_is_directed_right; }
+
+//! \brief obtain the segment source.
+template <typename Kernel>
+const typename Kernel::Point_2&
+My_arr_segment_traits_2<Kernel>::_Segment_cached_2::source() const { return m_ps; }
+
+//! \brief obtains the segment target.
+template <typename Kernel>
+const typename Kernel::Point_2&
+My_arr_segment_traits_2<Kernel>::_Segment_cached_2::target() const { return m_pt; }
+
+//! \brief obtains the (lexicographically) left endpoint.
+template <typename Kernel>
+const typename Kernel::Point_2&
+My_arr_segment_traits_2<Kernel>::_Segment_cached_2::left() const
+{ return (m_is_directed_right ? m_ps : m_pt); }
+
+//! \brief obtains the (lexicographically) right endpoint.
+template <typename Kernel>
+const typename Kernel::Point_2&
+My_arr_segment_traits_2<Kernel>::_Segment_cached_2::right() const
+{ return (m_is_directed_right ? m_pt : m_ps); }
+
+// Modifiers
+
+//! \brief sets the (lexicographically) left endpoint.
+template <typename Kernel>
+void My_arr_segment_traits_2<Kernel>::_Segment_cached_2::set_left(const Point_2& p)
+{
+ CGAL_precondition(! m_is_degen);
+ CGAL_precondition_code(Kernel kernel);
+ CGAL_precondition
+ (Segment_assertions::_assert_is_point_on(p, m_l, Has_exact_division()) &&
+ (kernel.compare_xy_2_object()(p, right()) == SMALLER));
+
+ if (m_is_directed_right) m_ps = p;
+ else m_pt = p;
+}
+
+//! \brief sets the (lexicographically) right endpoint.
+template <typename Kernel>
+void My_arr_segment_traits_2<Kernel>::_Segment_cached_2::set_right(const Point_2& p)
+{
+ CGAL_precondition(! m_is_degen);
+ CGAL_precondition_code(Kernel kernel);
+ CGAL_precondition
+ (Segment_assertions::_assert_is_point_on(p, m_l, Has_exact_division()) &&
+ (kernel.compare_xy_2_object()(p, left()) == LARGER));
+
+ if (m_is_directed_right) m_pt = p;
+ else m_ps = p;
+}
+
+//! \brief determines whether the given point is in the x-range of the segment.
+template <typename Kernel>
+bool My_arr_segment_traits_2<Kernel>::_Segment_cached_2::
+is_in_x_range(const Point_2& p) const
+{
+ Kernel kernel;
+ typename Kernel::Compare_x_2 compare_x = kernel.compare_x_2_object();
+ const Comparison_result res1 = compare_x(p, left());
+
+ if (res1 == SMALLER) return false;
+ else if (res1 == EQUAL) return true;
+
+ const Comparison_result res2 = compare_x(p, right());
+ return (res2 != LARGER);
+}
+
+//! \brief determines whether the given point is in the y-range of the segment.
+template <typename Kernel>
+bool My_arr_segment_traits_2<Kernel>::_Segment_cached_2::
+is_in_y_range(const Point_2& p) const
+{
+ Kernel kernel;
+ typename Kernel::Compare_y_2 compare_y = kernel.compare_y_2_object();
+ const Comparison_result res1 = compare_y(p, left());
+
+ if (res1 == SMALLER) return false;
+ else if (res1 == EQUAL) return true;
+
+ const Comparison_result res2 = compare_y(p, right());
+ return (res2 != LARGER);
+}
+
+/*! \class A representation of a segment, as used by the My_arr_segment_traits_2
+ * traits-class.
+ */
+template <typename Kernel_>
+class My_arr_segment_2 : public My_arr_segment_traits_2<Kernel_>::_Segment_cached_2 {
+ typedef Kernel_ Kernel;
+
+ typedef typename My_arr_segment_traits_2<Kernel>::_Segment_cached_2 Base;
+ typedef typename Kernel::Segment_2 Segment_2;
+ typedef typename Kernel::Point_2 Point_2;
+ typedef typename Kernel::Line_2 Line_2;
+
+public:
+ /*! Construct default. */
+ My_arr_segment_2();
+
+ /*! Construct a segment from a "kernel" segment.
+ * \param seg the segment.
+ * \pre the segment is not degenerate.
+ */
+ My_arr_segment_2(const Segment_2& seg, int original_segment_id);
+
+ /*! Construct a segment from two endpoints.
+ * \param source the source point.
+ * \param target the target point.
+ * \pre `source` and `target` are not equal.
+ */
+ My_arr_segment_2(const Point_2& source, const Point_2& target,
+ int original_segment_id);
+
+ /*! Construct a segment from a line and two endpoints.
+ * \param line the supporting line.
+ * \param source the source point.
+ * \param target the target point.
+ * \pre Both `source` and `target` must be on `line`.
+ * \pre `source` and `target` are not equal.
+ */
+ My_arr_segment_2(const Line_2& line,
+ const Point_2& source, const Point_2& target,
+ int original_segment_id);
+
+ /*! Construct a segment from all fields.
+ * \param line the supporting line.
+ * \param source the source point.
+ * \param target the target point.
+ * \param is_directed_right is (lexicographically) directed left to right.
+ * \param is_vert is the segment vertical.
+ * \param is_degen is the segment degenerate (a single point).
+ */
+ My_arr_segment_2(const Line_2& line,
+ const Point_2& source, const Point_2& target,
+ bool is_directed_right, bool is_vert, bool is_degen,
+ int original_segment_id);
+
+ int original_segment_id() const
+ { return m_original_segment_id; }
+
+ /*! Cast to a segment.
+ */
+ operator Segment_2() const;
+
+ /*! Flip the segment (swap its source and target).
+ */
+ My_arr_segment_2 flip() const;
+
+ /*! Create a bounding box for the segment.
+ */
+ CGAL::Bbox_2 bbox() const
+ {
+ Kernel kernel;
+ auto construct_bbox = kernel.construct_bbox_2_object();
+ return construct_bbox(this->m_ps) + construct_bbox(this->m_pt);
+ }
+
+ int m_original_segment_id;
+};
+
+//! \brief constructs default.
+template <typename Kernel>
+My_arr_segment_2<Kernel>::My_arr_segment_2() : Base() { }
+
+//! \brief constructs from a "kernel" segment.
+template <typename Kernel>
+My_arr_segment_2<Kernel>::My_arr_segment_2(const Segment_2& seg,
+ int original_segment_id) :
+ Base(seg),
+ m_original_segment_id(original_segment_id)
+{}
+
+//! \brief constructs a segment from two end-points.
+template <typename Kernel>
+My_arr_segment_2<Kernel>::My_arr_segment_2(const Point_2& source,
+ const Point_2& target,
+ int original_segment_id) :
+ Base(source, target),
+ m_original_segment_id(original_segment_id)
+{}
+
+//! \brief constructs a segment from a line and two end-points.
+template <typename Kernel>
+My_arr_segment_2<Kernel>::My_arr_segment_2(const Line_2& line,
+ const Point_2& source,
+ const Point_2& target,
+ int original_segment_id) :
+ Base(line,source, target),
+ m_original_segment_id(original_segment_id)
+{}
+
+//! \brief constructs a segment from all fields.
+template <typename Kernel>
+My_arr_segment_2<Kernel>::My_arr_segment_2(const Line_2& line,
+ const Point_2& source,
+ const Point_2& target,
+ bool is_directed_right,
+ bool is_vert, bool is_degen,
+ int original_segment_id) :
+ Base(line, source, target, is_directed_right, is_vert, is_degen),
+ m_original_segment_id(original_segment_id)
+{}
+
+//! \brief casts to a segment.
+template <typename Kernel>
+My_arr_segment_2<Kernel>::operator typename Kernel::Segment_2() const
+{
+ Kernel kernel;
+ auto seg_ctr = kernel.construct_segment_2_object();
+ return seg_ctr(this->source(), this->target());
+}
+
+//! \brief flips the segment (swap its source and target).
+template <typename Kernel>
+My_arr_segment_2<Kernel> My_arr_segment_2<Kernel>::flip() const
+{
+ return My_arr_segment_2(this->line(), this->target(), this->source(),
+ ! (this->is_directed_right()), this->is_vertical(),
+ this->is_degenerate(), m_original_segment_id);
+}
+
+/*! Exporter for the segment class used by the traits-class.
+ */
+template <typename Kernel, typename OutputStream>
+OutputStream& operator<<(OutputStream& os, const My_arr_segment_2<Kernel>& seg)
+{
+ os << static_cast<typename Kernel::Segment_2>(seg);
+ return (os);
+}
+
+/*! Importer for the segment class used by the traits-class.
+ */
+template <typename Kernel, typename InputStream>
+InputStream& operator>>(InputStream& is, My_arr_segment_2<Kernel>& seg)
+{
+ typename Kernel::Segment_2 kernel_seg;
+ is >> kernel_seg;
+ seg = kernel_seg;
+ return is;
+}
+
+#undef SMALLER
+#undef EQUAL
+#undef LARGER
+
+#include <CGAL/enable_warnings.h>
+
+#endif
diff --git a/src/My_construction_subcurve.h b/src/My_construction_subcurve.h
new file mode 100644
index 0000000000000000000000000000000000000000..9faba448bddd31570a49a55e7b2dbfb7daca47ca
--- /dev/null
+++ b/src/My_construction_subcurve.h
@@ -0,0 +1,184 @@
+// Copyright (c) 2006,2007,2009,2010,2011 Tel-Aviv University (Israel).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s) : Tali Zvi <talizvi@post.tau.ac.il>
+// Baruch Zukerman <baruchzu@post.tau.ac.il>
+
+#ifndef MY_CONSTRUCTION_SUBCURVE_H
+#define MY_CONSTRUCTION_SUBCURVE_H
+
+/*! \file
+ *
+ * Definition of the Arr_construction_subcurve class-template, which is an
+ * extended curve type, referred to as Subcurve, used by the surface-sweep
+ * framework.
+ *
+ * The surface-sweep framework is implemented as a template that is
+ * parameterized, among the other, by the Subcurve and Event types. That is,
+ * instance types of Subcurve and Event must be available when the
+ * surface-sweep template is instantiated.
+ *
+ * Arr_construction_subcurve derives from an instance of the Default_subcurve
+ * class template. The user is allowed to introduce new types that derive from
+ * an instance of the Arr_construction_subcurve class template. However, some of
+ * the fields of this template depends on the Subcurve type. We use the
+ * curiously recurring template pattern (CRTP) idiom to force the correct
+ * matching of these types.
+ */
+
+#include <CGAL/Surface_sweep_2/Default_subcurve.h>
+#include <CGAL/Default.h>
+#include "Hds.h"
+
+namespace Ss2 = CGAL::Surface_sweep_2;
+
+/*! \class Arr_construction_subcurve_base
+ *
+ * This is the base class of the Arr_construction_subcurve class template used
+ * by the (CRTP) idiom.
+ * \tparam GeometryTraits_2 the geometry traits.
+ * \tparam Event_ the event type.
+ * \tparam Allocator_ a type of an element that is used to acquire/release
+ * memory for elements of the event queue and the status
+ * structure, and to construct/destroy the elements in that
+ * memory. The type must meet the requirements of Allocator.
+ * \tparam Subcurve_ the subcurve actual type.
+ *
+ * The information contained in this class last:
+ * - ishe event that was handled on the curve.
+ * - The index for a subcurve that may represent a hole
+ * - Indices of all halfedge below the curve that may represent a hole.
+ */
+template <typename GeometryTraits_2, typename Event_, typename Allocator_,
+ typename Subcurve_>
+class My_construction_subcurve_base :
+ public CGAL::Surface_sweep_2::Default_subcurve_base<GeometryTraits_2, Event_, Allocator_, Subcurve_>
+{
+public:
+ typedef GeometryTraits_2 Geometry_traits_2;
+ typedef Subcurve_ Subcurve;
+ typedef Event_ Event;
+ typedef Allocator_ Allocator;
+
+private:
+ typedef Geometry_traits_2 Gt2;
+ typedef CGAL::Surface_sweep_2::Default_subcurve_base<Gt2, Event, Allocator, Subcurve>
+ Base;
+
+public:
+ typedef typename Gt2::X_monotone_curve_2 X_monotone_curve_2;
+ typedef Event* Event_ptr;
+
+ /*! Construct default. */
+ My_construction_subcurve_base() :
+ Base(),
+ m_last_event(nullptr),
+ m_original_segment_id(std::numeric_limits<int>::max()),
+ m_hh(nullptr)
+ {}
+
+ /*! Constructor from an x-monotone curve. */
+ My_construction_subcurve_base(X_monotone_curve_2& curve) :
+ Base(curve),
+ m_last_event(nullptr),
+ m_original_segment_id(curve.original_segment_id()),
+ m_hh(nullptr)
+ {}
+
+protected:
+ Event_ptr m_last_event; // The last event that was handled on the curve.
+ int m_original_segment_id;
+ Halfedge_handle m_hh; // Halfedge_handle of the edge between m_last_event, and lying
+ // on this curve
+
+public:
+ /*! Initialize the curve. */
+ void init(const X_monotone_curve_2& curve)
+ {
+ Base::init(curve);
+ m_original_segment_id=curve.original_segment_id();
+ }
+
+ /*! Set the event associated with the left end of the subcurve. */
+ template <typename SweepEvent>
+ void set_left_event(SweepEvent* left)
+ {
+ Base::set_left_event(left);
+ set_last_event(left);
+ }
+
+ /*! Set the last event on the subcurve. */
+ void set_last_event(Event_ptr e) { m_last_event = e; }
+
+ /*! Obtain the last event. */
+ Event_ptr last_event() const { return m_last_event; }
+
+ int original_segment_id() const
+ { return m_original_segment_id; }
+ void set_original_segment_id(int id)
+ { m_original_segment_id=id; }
+
+ Halfedge_handle halfedge() const
+ { return m_hh; }
+ void set_halfedge(Halfedge_handle hh)
+ { m_hh=hh; }
+};
+
+/*! \class Arr_construction_subcurve
+ *
+ * This class that holds information about a curve that is added to the
+ * arrangement. In addition to the information that is contained in
+ * Surface_sweep_subcurve, when an arrangement is constructed, a pointer to the
+ * last handled event on the curve is stored (in the base class). This
+ * information is used to retrieve hints when a subcurve of this curve is
+ * inserted into the planar map.
+ *
+ * Inherits from `Surface_sweep_subcurve`
+ * \sa `Surface_sweep_subcurve`
+ */
+template <typename GeometryTraits_2, typename Event_,
+ typename Allocator_ = CGAL_ALLOCATOR(int),
+ typename Subcurve_ = CGAL::Default>
+class My_construction_subcurve :
+ public My_construction_subcurve_base<
+ GeometryTraits_2, Event_, Allocator_,
+ typename CGAL::Default::Get<Subcurve_,
+ My_construction_subcurve<GeometryTraits_2, Event_,
+ Allocator_,
+ Subcurve_> >::type>
+{
+public:
+ typedef GeometryTraits_2 Geometry_traits_2;
+ typedef Event_ Event;
+ typedef Allocator_ Allocator;
+
+private:
+ typedef Geometry_traits_2 Gt2;
+ typedef My_construction_subcurve<Gt2, Event, Allocator,
+ Subcurve_> Self;
+ typedef typename CGAL::Default::Get<Subcurve_, Self>::type Subcurve;
+ typedef My_construction_subcurve_base<Gt2, Event, Allocator,
+ Subcurve> Base;
+
+public:
+ typedef typename Gt2::X_monotone_curve_2 X_monotone_curve_2;
+
+ typedef typename Base::Event_ptr Event_ptr;
+
+ /*! Construct deafult. */
+ My_construction_subcurve() {}
+
+ /*! Constructor from an x-monotone curve. */
+ My_construction_subcurve(X_monotone_curve_2& curve) :
+ Base(curve)
+ {}
+};
+
+
+#endif
diff --git a/src/My_event.h b/src/My_event.h
new file mode 100644
index 0000000000000000000000000000000000000000..631f46aa46ee2636753e59dc8738ea03f08e5381
--- /dev/null
+++ b/src/My_event.h
@@ -0,0 +1,120 @@
+// Copyright (c) 2006,2007,2009,2010,2011 Tel-Aviv University (Israel).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s) : Tali Zvi <talizvi@post.tau.ac.il>
+// Baruch Zukerman <baruchzu@post.tau.ac.il>
+// Efi Fogel <efif@post.tau.ac.il>
+
+#ifndef CGAL_MY_EVENT_H
+#define CGAL_MY_EVENT_H
+
+/*! \file
+ *
+ * Definition of the Arr_construction_event_base class-template.
+ */
+
+#include <vector>
+
+#include <CGAL/Surface_sweep_2/Default_event_base.h>
+#include <CGAL/assertions.h>
+#include "Hds.h"
+#include "My_construction_subcurve.h"
+
+/*! \class Arr_construction_event_base
+ *
+ * This template represents an event used by the surface-sweep framework. It
+ * inherits either from `Default_event_base` (the default) or
+ * 'No_overlap_event_base' depending on whether the curve may overlap or not.
+ * It stores the data associated with an event in addition to the information
+ * stored in its base class. When constructing an arrangement, additional
+ * information is stored, in order to expedite the insertion of curves into the
+ * arrangement.
+ *
+ * The additional infomation contains the following:
+ * - among the left curves of the event, we keep the highest halfedge that
+ * was inserted into the arrangement at any given time and when there are no
+ * left curves, we keep the highest halfedge that was inseted to the right.
+ *
+ * \tparam GeometryTraits_2 the geometry traits.
+ * \tparam Allocator_ a type of an element that is used to acquire/release
+ * memory for elements of the event queue and the status
+ * structure, and to construct/destroy the elements in that
+ * memory. The type must meet the requirements of Allocator.
+ * \tparam SurfaceSweepEvent a template, an instance of which is used as the
+ * base class.
+ *
+ * \sa `Default_event_base`
+ * \sa `No_overlap_event_base`
+ */
+ template <typename GeometryTraits_2, typename Subcurve_,
+ template <typename, typename>
+ class SurfaceSweepEvent = CGAL::Surface_sweep_2::Default_event_base>
+class My_event_base :
+ public SurfaceSweepEvent<GeometryTraits_2, Subcurve_>
+{
+public:
+ typedef GeometryTraits_2 Geometry_traits_2;
+ typedef Subcurve_ Subcurve;
+
+private:
+ typedef Geometry_traits_2 Gt2;
+ typedef SurfaceSweepEvent<Gt2, Subcurve> Base;
+ typedef My_event_base<Gt2, Subcurve, SurfaceSweepEvent> Self;
+
+public:
+ typedef typename Gt2::X_monotone_curve_2 X_monotone_curve_2;
+ typedef typename Gt2::Point_2 Point_2;
+
+ typedef typename Base::Subcurve_container Subcurve_container;
+ typedef typename Base::Subcurve_iterator Subcurve_iterator;
+ typedef typename Base::Subcurve_reverse_iterator Subcurve_reverse_iterator;
+
+protected:
+ Vertex_handle m_vertex; // The vertex handle associated with this event.
+ Halfedge_handle m_above_halfedge; // The halfedge handle of the subcurve just above this event.
+ bool m_before_strip;
+public:
+ /*! Default constructor. */
+ My_event_base(): m_vertex(nullptr), m_above_halfedge(nullptr), m_before_strip(false)
+ {}
+
+ /*! Destructor */
+ ~My_event_base() {}
+
+ /*! Set the vertex handle. */
+ void set_vertex_handle(Vertex_handle vh) { m_vertex=vh; }
+
+ /*! Get the vertex handle. */
+ Vertex_handle vertex_handle() const { return m_vertex; }
+
+ void set_above_halfedge(Halfedge_handle hh)
+ { m_above_halfedge=hh; }
+ Halfedge_handle above_halfedge() const
+ { return m_above_halfedge; }
+
+ bool before_strip() const
+ { return m_before_strip; }
+ void set_before_strip()
+ { m_before_strip=true; }
+};
+
+template <typename GeometryTraits_2,
+ typename Allocator_ = CGAL_ALLOCATOR(int)>
+class My_event :
+ public My_event_base<GeometryTraits_2,
+ My_construction_subcurve<GeometryTraits_2,
+ My_event<GeometryTraits_2, Allocator_>,
+ Allocator_> >
+{
+public:
+ /*! Construct default. */
+ My_event() {}
+};
+
+#endif
diff --git a/src/My_surface_sweep.h b/src/My_surface_sweep.h
new file mode 100644
index 0000000000000000000000000000000000000000..835e0d089a524dc215a362e719f09de205f748b1
--- /dev/null
+++ b/src/My_surface_sweep.h
@@ -0,0 +1,115 @@
+// compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#ifndef MY_SURFACE_SWEEP_H
+#define MY_SURFACE_SWEEP_H
+
+#include <CGAL/Surface_sweep_2.h>
+
+template <typename Visitor_>
+class My_surface_sweep_2:
+ CGAL::Surface_sweep_2::Surface_sweep_2<Visitor_>
+{
+public:
+ typedef My_surface_sweep_2<Visitor_> Self;
+ using Base=CGAL::Surface_sweep_2::Surface_sweep_2<Visitor_>;
+
+ using Geometry_traits_2=typename Base::Geometry_traits_2;
+ using Point_2=typename Base::Point_2;
+ using Event_queue_iterator=typename Base::Event_queue_iterator;
+ using Event_subcurve_iterator=typename Base::Event_subcurve_iterator;
+ using Subcurve=typename Base::Subcurve;
+
+ My_surface_sweep_2(Visitor_* visitor) : Base(visitor) {}
+ My_surface_sweep_2(const Geometry_traits_2* traits, Visitor_* visitor) :
+ Base(traits, visitor) {}
+
+ template <typename CurveInputIterator>
+ void sweep_until(CurveInputIterator curves_begin,
+ CurveInputIterator curves_end,
+ const Point_2& p)
+ {
+ this->m_visitor->before_sweep();
+ this->_init_sweep(curves_begin, curves_end);
+ //m_visitor->after_init();
+ this->_sweep_until(p);
+ this->_complete_sweep();
+ this->m_visitor->after_sweep();
+ }
+
+protected:
+ void _sweep_until(const Point_2& p)
+ {
+ // Looping over the events in the queue.
+ Event_queue_iterator eventIter = this->m_queue->begin();
+
+ std::size_t number_of_strips=this->m_visitor->number_of_strips();
+ std::size_t strip_id=this->m_visitor->strip_id();
+ bool before_strip=true;
+ bool after_strip=false;
+
+ if (strip_id==0)
+ {
+ before_strip=false;
+ this->m_visitor->enter_in_strip();
+ }
+
+ while (eventIter != this->m_queue->end()) {
+ // Get the next event from the queue.
+ this->m_currentEvent = *eventIter;
+
+ if (number_of_strips>1)
+ {
+ if (before_strip)
+ {
+ before_strip=(CGAL::compare_x(this->m_currentEvent->point(),
+ this->m_visitor->pmin())==CGAL::SMALLER);
+ if (!before_strip)
+ { this->m_visitor->enter_in_strip(); }
+ }
+ if (!after_strip && strip_id<number_of_strips-1)
+ {
+ after_strip=(CGAL::compare_x(this->m_currentEvent->point(),
+ this->m_visitor->pmax())!=CGAL::SMALLER);
+ if (after_strip) { this->m_visitor->end_of_strip(); }
+ }
+ }
+
+ if (!after_strip)
+ {
+ this->_handle_left_curves();
+ this->_handle_right_curves();
+ }
+
+ if (this->m_visitor->after_handle_event(this->m_currentEvent,
+ this->m_status_line_insert_hint,
+ this->m_is_event_on_above))
+ {
+ this->deallocate_event(this->m_currentEvent);
+ }
+
+ this->m_queue->erase(eventIter);
+ eventIter = this->m_queue->begin();
+ }
+
+ this->m_statusLine.clear();
+ }
+};
+
+#endif // MY_SURFACE_SWEEP_H
diff --git a/src/My_visitor.h b/src/My_visitor.h
new file mode 100644
index 0000000000000000000000000000000000000000..dbad32e1230aa77ce31786b3080fdb7278de360f
--- /dev/null
+++ b/src/My_visitor.h
@@ -0,0 +1,319 @@
+// Copyright (c) 2005,2006,2007,2009,2010,2011 Tel-Aviv University (Israel).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s): Baruch Zukerman <baruchzu@post.tau.ac.il>
+// Efi Fogel <efif@post.tau.ac.il>
+// (based on old version by Tali Zvi)
+
+#ifndef CGAL_MY_VISITOR_H
+#define CGAL_MY_VISITOR_H
+
+/*! \file
+ *
+ * Definition of a surface-sweep visitor that reports all maximal x-monotone
+ * non-intersecting x-monotone curves induced by a set of input curves.
+ */
+
+#include <vector>
+
+#include <CGAL/Surface_sweep_2/Default_visitor.h>
+#include <CGAL/Surface_sweep_2/Surface_sweep_2_utils.h>
+#include "Hds.h"
+#include "My_surface_sweep.h"
+#include "My_event.h"
+#include "My_construction_subcurve.h"
+#include "Partial_hds.h"
+
+/*! \class Subcurves_visitor
+ *
+ * A simple sweep-line visitor that reports all maximal x-monotone
+ * non-intersecting x-monotone curves induced by a set of input curves. Used by
+ * compute_subcurves().
+ */
+ template <typename GeometryTraits_2,
+ typename Allocator_ = CGAL_ALLOCATOR(int)>
+class My_visitor :
+ public CGAL::Surface_sweep_2::Default_visitor
+ <My_visitor<GeometryTraits_2, Allocator_>,
+ GeometryTraits_2,
+ Allocator_,
+ My_event<GeometryTraits_2, Allocator_>,
+ My_construction_subcurve<GeometryTraits_2,
+ My_event<GeometryTraits_2, Allocator_>,
+ Allocator_>
+ >
+{
+public:
+ typedef GeometryTraits_2 Geometry_traits_2;
+ typedef Allocator_ Allocator;
+
+private:
+ typedef Geometry_traits_2 Gt2;
+ typedef My_visitor<Gt2, Allocator> Self;
+ typedef CGAL::Surface_sweep_2::Default_visitor
+ <Self, Gt2, Allocator,
+ My_event<GeometryTraits_2, Allocator_>,
+ My_construction_subcurve<GeometryTraits_2,
+ My_event<GeometryTraits_2, Allocator_>,
+ Allocator_>> Base;
+
+public:
+ typedef typename Gt2::X_monotone_curve_2 X_monotone_curve_2;
+ typedef typename Gt2::Point_2 Point_2;
+
+ typedef typename Base::Event Event;
+ typedef typename Base::Subcurve Subcurve;
+
+ typedef typename Subcurve::Status_line_iterator Status_line_iterator;
+
+protected:
+ // Data members:
+ Partial_hds& m_partial_hds;
+ bool m_before_strip, m_after_strip;
+ std::size_t m_nb_intersections;
+
+public:
+ My_visitor(Partial_hds& a_partial_hds) :
+ m_partial_hds(a_partial_hds),
+ m_before_strip(true),
+ m_after_strip(false),
+ m_nb_intersections(0)
+ {}
+
+ template<typename Container>
+ void sweep(const Container& all_segments, bool crop=false, std::size_t id=1)
+ {
+ m_partial_hds.m_nb_input_segments=0;
+ std::vector<X_monotone_curve_2> curves_vec;
+ if (m_partial_hds.number_of_strips()==1)
+ { curves_vec.reserve(all_segments.size()); }
+ else
+ {
+ curves_vec.reserve((2*all_segments.size())/(m_partial_hds.number_of_strips()));
+ }
+
+ bool concerned=false, cropx1, cropx2;
+ EPECK::Line_2 line;
+ double val;
+ ECoord x1, y1;
+ ECoord x2, y2;
+ auto it=all_segments.begin();
+ for (; it!=all_segments.end(); ++id, ++it)
+ {
+ if (it->source()!=it->target())
+ {
+ concerned=false;
+
+ if (it->source().x()<m_partial_hds.m_minx)
+ {
+ if (it->target().x()>=m_partial_hds.m_minx)
+ { concerned=true; }
+ }
+ else
+ {
+ if (it->source().x()<m_partial_hds.m_maxx)
+ { concerned=true; }
+ }
+
+ if (concerned)
+ {
+ if (crop)
+ {
+ if (it->source().x()<m_partial_hds.m_minx) { cropx1=true; }
+ else { cropx1=false; }
+ if (it->target().x()>m_partial_hds.m_maxx) { cropx2=true; }
+ else { cropx2=false; }
+
+ if (cropx1 || cropx2)
+ {
+ line=EPECK::Line_2(EPoint(it->source().x(), it->source().y()),
+ EPoint(it->target().x(), it->target().y()));
+ if (cropx1)
+ {
+ val=m_partial_hds.m_minx-
+ ((m_partial_hds.m_maxx-m_partial_hds.m_minx)/10000.);
+ if (val<it->source().x())
+ { cropx1=false; }
+ else
+ { x1=val; y1=line.y_at_x(x1); } // Exact coordinates
+ }
+ if (cropx2)
+ { x2=m_partial_hds.m_maxx; y2=line.y_at_x(x2); }
+ }
+
+ if (!cropx1)
+ { x1=it->source().x(); y1=it->source().y(); }
+ if (!cropx2)
+ { x2=it->target().x(); y2=it->target().y(); }
+ curves_vec.push_back(X_monotone_curve_2(ESegment(EPoint(x1,y1),
+ EPoint(x2, y2)),
+ id));
+ }
+ else
+ {
+ curves_vec.push_back(X_monotone_curve_2
+ (ESegment(EPoint(it->source().x(), it->source().y()),
+ EPoint(it->target().x(), it->target().y())),
+ id));
+ }
+ ++m_partial_hds.m_nb_input_segments;
+ }
+ }
+ }
+
+ // 2) Perform the sweep.
+ My_surface_sweep_2<Self>* sl = reinterpret_cast<My_surface_sweep_2<Self>*>
+ (this->surface_sweep());
+ sl->sweep_until(curves_vec.begin(), curves_vec.end(),
+ m_partial_hds.m_max);
+ }
+
+ void enter_in_strip()
+ { m_before_strip=false; }
+
+ void end_of_strip()
+ { m_after_strip=true; }
+
+ HDS& hds()
+ { return m_partial_hds.hds(); }
+
+ const EPoint& pmin() const
+ { return m_partial_hds.pmin(); }
+ const EPoint& pmax() const
+ { return m_partial_hds.pmax(); }
+
+ std::size_t strip_id() const
+ { return m_partial_hds.strip_id(); }
+ std::size_t number_of_strips() const
+ { return m_partial_hds.number_of_strips(); }
+
+ std::size_t number_of_intersections() const
+ { return m_nb_intersections; }
+
+ // Before to process the given event in the sweep line
+ void before_handle_event(Event* event)
+ {
+ assert(event->is_closed());
+ assert(event->vertex_handle()==nullptr);
+
+ if (m_before_strip)
+ { event->set_before_strip(); return; }
+
+ if (m_after_strip) return;
+
+ assert(m_partial_hds.point_in_strip(event->point()));
+ event->set_vertex_handle(hds().create_vertex(event->point()));
+
+ /* print_txt_with_endl("[strip", m_partial_hds.strip_id(),
+ "] Create inner point for [",event->point()); */
+ }
+
+ // After having proceed the given even in the sweep line
+ bool after_handle_event(Event* event,
+ Status_line_iterator iter,
+ bool /*is_event_on_above*/)
+ {
+ // print_txt_with_endl("after_handle_event ", event, " ", event->point());
+ Halfedge_handle prevhh=nullptr;
+ if (event->vertex_handle()!=nullptr)
+ { // Point in strip
+ typename Event::Subcurve_reverse_iterator subcurve_rit;
+ for (subcurve_rit=event->left_curves_rbegin();
+ subcurve_rit!=event->left_curves_rend(); ++subcurve_rit)
+ {
+ hds().set_vertex(hds().opposite((*subcurve_rit)->halfedge()),
+ event->vertex_handle());
+ prevhh=hds().add_segment_around_vertex
+ (hds().opposite((*subcurve_rit)->halfedge()),
+ event->vertex_handle(), prevhh);
+ if ((*subcurve_rit)->last_event()->vertex_handle()==nullptr)
+ { m_partial_hds.add_left_external_halfedge((*subcurve_rit)->halfedge()); }
+ }
+
+ if (event->is_intersection() || event->is_weak_intersection())
+ { ++m_nb_intersections; }
+ }
+ else if (m_after_strip)
+ {
+ typename Event::Subcurve_reverse_iterator subcurve_rit;
+ for (subcurve_rit=event->left_curves_rbegin();
+ subcurve_rit!=event->left_curves_rend(); ++subcurve_rit)
+ {
+ if ((*subcurve_rit)->last_event()->vertex_handle()!=nullptr ||
+ (*subcurve_rit)->last_event()->before_strip())
+ {
+ /* print_txt_with_endl("right_segment ",
+ (*subcurve_rit)->last_event()->point(),
+ " -> ", event->point()); */
+ m_partial_hds.add_right_external_halfedge
+ (hds().opposite((*subcurve_rit)->halfedge()));
+
+ if ((*subcurve_rit)->last_event()->before_strip())
+ { m_partial_hds.add_left_external_halfedge((*subcurve_rit)->halfedge()); }
+ }
+ }
+
+ typename Event::Subcurve_iterator subcurve_it;
+ for (subcurve_it=event->right_curves_begin();
+ subcurve_it!=event->right_curves_end(); ++subcurve_it)
+ { (*subcurve_it)->set_last_event(event); }
+ }
+ else
+ {
+ assert(m_before_strip);
+ typename Event::Subcurve_reverse_iterator subcurve_rit;
+ for (subcurve_rit=event->left_curves_rbegin();
+ subcurve_rit!=event->left_curves_rend(); ++subcurve_rit)
+ {
+ // Case of an edge entirely before the strip: ignore
+ hds().erase_edge((*subcurve_rit)->halfedge());
+ (*subcurve_rit)->set_halfedge(nullptr);
+ }
+ }
+
+ if (!m_after_strip)
+ {
+ typename Event::Subcurve_iterator subcurve_it;
+ for (subcurve_it=event->right_curves_begin();
+ subcurve_it!=event->right_curves_end(); ++subcurve_it)
+ {
+ (*subcurve_it)->set_halfedge
+ (hds().create_edge
+ (event->vertex_handle(), (*subcurve_it)->original_segment_id()));
+ if (event->vertex_handle()!=nullptr)
+ {
+ prevhh=hds().add_segment_around_vertex
+ ((*subcurve_it)->halfedge(), event->vertex_handle(), prevhh);
+ }
+ (*subcurve_it)->set_last_event(event);
+ }
+ }
+
+ if (event->vertex_handle()!=nullptr)
+ {
+ hds().end_of_vertex(event->vertex_handle(), prevhh);
+ // Iter is on the status line just after event
+ Halfedge_handle fatherdart=(iter!=this->status_line_end()?
+ (*iter)->halfedge():nullptr);
+ m_partial_hds.set_above_halfedge(event->vertex_handle(), fatherdart);
+ }
+
+ return (event->number_of_right_curves() == 0);
+ }
+
+ void found_overlap(Subcurve* sc1,
+ Subcurve* sc2,
+ Subcurve* ov_sc)
+ {
+ ov_sc->set_original_segment_id(std::min(sc1->original_segment_id(),
+ sc2->original_segment_id()));
+ }
+};
+
+#endif
diff --git a/src/Partial_hds.h b/src/Partial_hds.h
new file mode 100644
index 0000000000000000000000000000000000000000..8f7c6675b77abb7aea80dd68e3b4dc82267f02e1
--- /dev/null
+++ b/src/Partial_hds.h
@@ -0,0 +1,336 @@
+// compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#ifndef PARTIAL_HDS_H
+#define PARTIAL_HDS_H
+
+#include <iostream>
+#include <thread>
+#include <vector>
+#include <tuple>
+#include <stack>
+#include <map>
+
+#include "Print_txt.h"
+#include "Hds.h"
+#include "Export_xfig.hh"
+
+template <typename GeometryTraits_2, typename Allocator_>
+class My_visitor;
+template <typename Kernel_>
+class My_arr_segment_traits_2;
+class SHds;
+
+////////////////////////////////////////////////////////////////////////////////
+class Partial_hds
+{
+public:
+ using Self=Partial_hds;
+
+ template <typename GeometryTraits_2, typename Allocator_>
+ friend class My_visitor;
+
+ friend class SHds;
+ friend class MySimpleSHDSViewerQt;
+
+ Partial_hds(double x1=0., double y1=0., double x2=0., double y2=0.,
+ std::size_t strip_id=0, std::size_t nb_strips=0):
+ m_minx(x1),
+ m_maxx(x2),
+ m_min(x1, y1),
+ m_max(x2, y2),
+ m_strip_id(strip_id),
+ m_nb_strips(nb_strips),
+ m_nb_input_segments(0)
+ {
+ /*print_txt_with_endl("[Strip ", m_strip_id, "]: ", x1, ", ",
+ y1, ", ", x2, ", ", y2); */
+ }
+
+ std::size_t strip_id() const
+ { return m_strip_id; }
+ std::size_t number_of_strips() const
+ { return m_nb_strips; }
+
+ ECoord xmin() const
+ { return m_min.x(); }
+ ECoord ymin() const
+ { return m_min.y(); }
+ ECoord xmax() const
+ { return m_max.x(); }
+ ECoord ymax() const
+ { return m_max.y(); }
+
+ const EPoint& pmin() const
+ { return m_min; }
+ const EPoint& pmax() const
+ { return m_max; }
+
+ HDS& hds()
+ { return m_hds; }
+ const HDS& hds() const
+ { return m_hds; }
+
+ std::size_t number_of_left_externals() const
+ { return m_left_external_halfedges.size(); }
+ std::size_t number_of_right_externals() const
+ { return m_right_external_halfedges.size(); }
+ std::size_t number_of_externals() const
+ { return number_of_left_externals()+number_of_right_externals(); }
+ std::size_t number_of_non_externals() const
+ { return hds().number_of_halfedges()-number_of_externals(); }
+ std::size_t number_of_input_segments() const
+ { return m_nb_input_segments; }
+
+ void clear()
+ {
+ m_minx=m_maxx=0;
+ m_min=m_max=CGAL::ORIGIN;
+ m_strip_id=m_nb_strips=0;
+ m_hds.clear();
+ m_left_external_halfedges.clear();
+ m_right_external_halfedges.clear();
+ }
+
+ bool point_in_strip(const EPoint& p) const
+ {
+ // in strip means p.x >= m_min.x && p.x < m_max.x
+ return (CGAL::compare_x(p, m_min)!=CGAL::SMALLER &&
+ CGAL::compare_x(p, m_max)==CGAL::SMALLER);
+ }
+
+ void add_left_external_halfedge(Halfedge_handle hh)
+ {
+ /*print_txt_with_endl("[strip ", strip_id(), "] add_left_external_halfedge ",
+ hds().edge_id(hh), " ptr=", &*hh);*/
+
+ assert(is_left_external(hh));
+ assert(m_left_external_halfedges.count(hds().edge_id(hh))==0);
+ assert(hds().left_to_right(hh));
+ m_left_external_halfedges.insert(std::make_pair(hds().edge_id(hh), hh));
+ }
+
+ void add_right_external_halfedge(Halfedge_handle hh)
+ {
+ /*print_txt_with_endl("[strip ", strip_id(), "] add_right_external_halfedge ",
+ -hds().edge_id(hh), " ptr=", &*hh);*/
+
+ assert(is_right_external(hh));
+ assert(m_right_external_halfedges.count(-hds().edge_id(hh))==0);
+ assert(hds().right_to_left(hh));
+
+ m_right_external_halfedges.insert(std::make_pair(-hds().edge_id(hh), hh));
+ }
+
+ void set_above_halfedge(Vertex_handle vh, Halfedge_handle hh)
+ { vh->set_above_halfedge(hh); }
+
+ // @return true iff hh is external, i.e. not linked with an inner point in
+ // this strip.
+ bool is_external(Halfedge_handle hh) const
+ { return hds().vertex(hh)==nullptr; }
+
+ // @return true iff hh is left external, i.e. external and intersect the
+ // left side of the strip (and thus from left to right).
+ bool is_left_external(Halfedge_handle hh) const
+ { return is_external(hh) && hds().left_to_right(hh); }
+
+ // @return true iff hh is right external, i.e. external and intersect the
+ // right side of the strip (and thus from right to left).
+ bool is_right_external(Halfedge_handle hh) const
+ { return is_external(hh) && hds().right_to_left(hh); }
+
+ bool source_in_strip(Halfedge_handle hh) const
+ { return !is_external(hh); }
+ bool target_in_strip(Halfedge_handle hh) const
+ { return !is_external(hds().opposite(hh)); }
+
+ bool is_valid_external_halfedges()
+ {
+ bool res=true;
+
+ if (m_strip_id==0 && !m_left_external_halfedges.empty())
+ {
+ print_txt_with_endl("[strip 0] : ERROR m_left_external_halfedges "
+ "is not empty.");
+ }
+ if (m_strip_id==m_nb_strips-1 && !m_right_external_halfedges.empty())
+ {
+ print_txt_with_endl("[strip ", strip_id(),
+ "] (last strip) : ERROR, m_right_external_halfedges "
+ "is not empty.");
+ }
+
+ for (auto it=hds().halfedges().begin(), itend=hds().halfedges().end();
+ it!=itend; ++it)
+ {
+ if (is_left_external(it))
+ {
+ auto resfind=m_left_external_halfedges.find(hds().edge_id(it));
+ if (resfind==m_left_external_halfedges.end())
+ { print_txt_with_endl("[strip ", strip_id(),
+ "] : ERROR, halfedge ",
+ hds().edge_id(it), " is left external but ",
+ "is not in array m_left_external_halfedges.");
+ res=false;
+ }
+ else
+ {
+ if (resfind->second!=it)
+ { print_txt_with_endl("[strip ", strip_id(),
+ "] : ERROR, halfedge ",
+ hds().edge_id(it), " is left external but ",
+ "its associated halfedge in array ",
+ "m_left_external_halfedges is different ",
+ hds().edge_id(resfind->second),
+ " PTR ", &*(resfind->second), " AND ",
+ &*(it));
+ res=false;
+ }
+ }
+ }
+ if (is_right_external(it))
+ {
+ auto resfind=m_right_external_halfedges.find(-hds().edge_id(it));
+ if (resfind==m_right_external_halfedges.end())
+ { print_txt_with_endl("[strip ", strip_id(),
+ "] : ERROR, halfedge ",
+ hds().edge_id(it), " is right external but ",
+ "is not in array m_right_external_halfedges.");
+ res=false;
+ }
+ else
+ {
+ if (resfind->second!=it)
+ { print_txt_with_endl("[strip ", strip_id(),
+ "] : ERROR, halfedge ",
+ hds().edge_id(it), " is right external but ",
+ "its associated halfedge in array ",
+ "m_right_external_halfedges is different ",
+ hds().edge_id(resfind->second));
+ res=false;
+ }
+ }
+ }
+ }
+ return res;
+ }
+
+ friend std::ostream& operator<<(std::ostream& os,
+ const Partial_hds& la)
+ {
+ CGAL::exact(la.m_min); CGAL::exact(la.m_max);
+ os<<std::setprecision(17)
+ <<la.m_minx<<" "<<la.m_maxx<<" "<<la.m_min<<" "<<la.m_max<<" "
+ <<la.m_strip_id<<" "<<la.m_nb_strips<<" "<<la.m_nb_input_segments
+ <<std::endl<<std::endl;
+ os<<la.m_hds<<std::endl;
+
+ os<<la.m_left_external_halfedges.size()<<" "
+ <<la.m_right_external_halfedges.size()<<std::endl;
+
+ for (auto it=la.m_left_external_halfedges.begin(),
+ itend=la.m_left_external_halfedges.end(); it!=itend; ++it)
+ { os<<it->first<<" "<<la.hds().halfedges().index(it->second)<<std::endl; }
+ for (auto it=la.m_right_external_halfedges.begin(),
+ itend=la.m_right_external_halfedges.end(); it!=itend; ++it)
+ { os<<it->first<<" "<<la.hds().halfedges().index(it->second)<<std::endl; }
+
+ return os;
+ }
+
+ friend std::istream& operator>>(std::istream& is,
+ Partial_hds& la)
+ {
+ la.clear();
+
+ is>>la.m_minx>>la.m_maxx>>la.m_min>>la.m_max
+ >>la.m_strip_id>>la.m_nb_strips>>la.m_nb_input_segments;
+
+ std::vector<Halfedge_handle> halfedges;
+ std::vector<Vertex_handle> vertices;
+ la.m_hds.load(is, halfedges, vertices);
+
+ std::size_t nb_left, nb_right;
+ std::size_t ind1;
+ int int1;
+
+ is>>nb_left>>nb_right;
+ for (std::size_t i=0; i<nb_left; ++i)
+ {
+ is>>int1>>ind1;
+ la.m_left_external_halfedges[int1]=halfedges[ind1];
+ }
+ for (std::size_t i=0; i<nb_right; ++i)
+ {
+ is>>int1>>ind1;
+ la.m_right_external_halfedges[int1]=halfedges[ind1];
+ }
+ return is;
+ }
+
+ void display_info() const
+ {
+ print_txt_with_endl("[strip ", m_strip_id, "]: (", m_minx, " -> ", m_maxx,
+ ") #input-segments ", m_nb_input_segments,
+ " #left-external=", m_left_external_halfedges.size(),
+ " #right-external=", m_right_external_halfedges.size(),
+ " #halfedges=", hds().number_of_halfedges(),
+ " #vertices=", hds().number_of_vertices());
+ }
+
+ void export_to_xfig(Export_xfig& export_xfig,
+ double dx, double dy,
+ uint8_t coul_segments,
+ uint8_t coul_disks,
+ uint8_t coul_borders,
+ uint16_t width_segments,
+ uint16_t radius_disks,
+ uint16_t width_borders) const
+ {
+ // Border of the strip
+ if (width_borders>0)
+ {
+ export_xfig.rectangle(CGAL::to_double(m_min.x())+dx,
+ CGAL::to_double(m_min.y())+dy,
+ CGAL::to_double(m_max.x())+dx,
+ CGAL::to_double(m_max.y())+dy,
+ coul_borders, width_borders, false, 200);
+ }
+
+ // HDS
+ hds().export_to_xfig(export_xfig, dx, dy,
+ coul_segments, coul_disks,
+ width_segments, radius_disks);
+ }
+
+protected:
+ double m_minx, m_maxx;
+ EPoint m_min, m_max;
+ std::size_t m_strip_id;
+ std::size_t m_nb_strips;
+ std::size_t m_nb_input_segments;
+
+ HDS m_hds;
+ std::unordered_map<int, Halfedge_handle> m_left_external_halfedges;
+ std::unordered_map<int, Halfedge_handle> m_right_external_halfedges;
+};
+////////////////////////////////////////////////////////////////////////////////
+
+#endif // PARTIAL_HDS_H
diff --git a/src/Print_txt.h b/src/Print_txt.h
new file mode 100644
index 0000000000000000000000000000000000000000..2844c34c8fab5e456fbdd0b1b35372dd7e8d8aef
--- /dev/null
+++ b/src/Print_txt.h
@@ -0,0 +1,62 @@
+// compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#ifndef PRINT_TXT_H
+#define PRINT_TXT_H
+
+#include <sstream>
+#include <utility>
+#include <chrono>
+
+///////////////////////////////////////////////////////////////////////////////
+template<typename T>
+void put_one_value(std::ostringstream& txt, const T& v)
+{ txt<<v; }
+//----------------------------------------------------------------------------
+template<>
+void put_one_value(std::ostringstream& txt,
+ const std::chrono::duration<double>& v)
+{
+ std::chrono::seconds s=std::chrono::duration_cast<std::chrono::seconds>(v);
+ txt<<v.count();
+}
+///////////////////////////////////////////////////////////////////////////////
+template <typename Arg, typename... Args>
+void print_txt(Arg&& arg, Args&&... args)
+{
+ std::ostringstream txt;
+ put_one_value(txt, std::forward<Arg>(arg));
+ using expander = int[];
+ (void)expander{0, (void(put_one_value(txt, std::forward<Args>(args))),0)...};
+ std::cout<<txt.str()<<std::flush;
+}
+///////////////////////////////////////////////////////////////////////////////
+template <typename Arg, typename... Args>
+void print_txt_with_endl(Arg&& arg, Args&&... args)
+{
+ std::ostringstream txt;
+ put_one_value(txt, std::forward<Arg>(arg));
+ using expander = int[];
+ (void)expander{0, (void(put_one_value(txt, std::forward<Args>(args))),0)...};
+ txt<<std::endl;
+ std::cout<<txt.str()<<std::flush;
+}
+///////////////////////////////////////////////////////////////////////////////
+
+#endif // PRINT_TXT_H
diff --git a/src/SHds.h b/src/SHds.h
new file mode 100644
index 0000000000000000000000000000000000000000..90229ca773bd884c56bcafdc8b0fac489e21384d
--- /dev/null
+++ b/src/SHds.h
@@ -0,0 +1,1195 @@
+// compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#ifndef SHDS_H
+#define SHDS_H
+
+#include <iostream>
+#include <thread>
+#include <vector>
+#include <mutex>
+#include <boost/filesystem.hpp>
+
+#include <CGAL/Compact_container.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include "My_construction_subcurve.h"
+#include "My_arr_segment_traits_2.h"
+#include "My_surface_sweep.h"
+#include "My_visitor.h"
+#include "memory.h"
+#include "CGAL_typedefs.h"
+#include "Partial_hds.h"
+#include "Hds.h"
+
+////////////////////////////////////////////////////////////////////////////////
+class SHds
+{
+public:
+ using Self=SHds;
+
+ typedef uint32_t size_type;
+ static const size_type NB_MARKS=Halfedge::NB_MARKS;
+ static const size_type INVALID_MARK=NB_MARKS;
+ class Exception_no_more_available_mark {};
+
+ SHds():
+ m_xmin(0), m_ymin(0), m_xmax(0), m_ymax(0),
+ m_nb_strips(0),
+ m_number_segments(0),
+ m_number_of_halfedges(0),
+ m_nb_finite_faces(0)
+ { init_all_marks(); }
+
+ void init_all_marks()
+ {
+ mnb_used_marks = 0;
+ mmask_marks.reset();
+ for (size_type i=0; i<NB_MARKS; ++i)
+ {
+ mfree_marks_stack[i] = i;
+ mindex_marks[i] = i;
+ mnb_marked_halfedges[i] = 0;
+ mnb_times_reserved_marks[i] = 0;
+ }
+ }
+
+ SHds(const std::string& filename,
+ std::size_t nb_strips,
+ std::size_t nb_threads,
+ bool optimized_strips,
+ std::vector<double>& times,
+ bool crop=false,
+ std::vector<double>* xpos_cuts=nullptr):
+ m_xmin(0), m_ymin(0), m_xmax(0), m_ymax(0),
+ m_nb_strips(nb_strips),
+ m_number_segments(0),
+ m_number_of_halfedges(0),
+ m_nb_finite_faces(0)
+ {
+ init_all_marks();
+
+ if (times.empty())
+ { times.resize(4, 0); }
+
+ // times[0] : global time
+ // times[1] : load and dispatch
+ // times[2] : compute all partial arrangements
+ // times[3] : compute faces
+
+ My_timer global_timer; global_timer.start();
+ My_timer local_timer;
+ local_timer.start();
+
+ std::vector<ISegment> all_segments;
+ std::vector<double> xpos_strips;
+
+ if (xpos_cuts!=nullptr)
+ { compute_strips_v0(filename, all_segments, xpos_strips, xpos_cuts);}
+ else if (!optimized_strips)
+ { compute_strips_v1(filename, all_segments, xpos_strips); }
+ else
+ { compute_strips_v2(filename, all_segments, xpos_strips); }
+
+ std::cout<<"strips xpos=(";
+ for (std::size_t i=0; i<xpos_strips.size(); ++i)
+ { std::cout<<xpos_strips[i]<<" "; }
+ std::cout<<") #size="<<xpos_strips.size()<<std::endl;
+
+ // Test if definition of strips is valid.
+ bool strip_error=(xpos_strips.size()!=number_of_strips()+1);
+ for (std::size_t i=1; i<nb_strips; ++i)
+ {
+ if (xpos_strips[i]<=m_xmin || xpos_strips[i]>=m_xmax ||
+ xpos_strips[i]<=xpos_strips[i-1])
+ { strip_error=true; }
+ }
+ if (xpos_strips.front()!=m_xmin || xpos_strips.back()!=m_xmax)
+ { strip_error=true; }
+ if (strip_error)
+ {
+ std::cerr<<"[ERROR] in definition of strips: #nbstrips="<<number_of_strips()
+ <<std::endl;
+ exit(EXIT_FAILURE);
+ }
+
+ local_timer.stop();
+ times[1]+=local_timer.time();
+
+ // Now compute partial hds, possibly in parallel
+ m_tab_partial_hds.resize(number_of_strips());
+ std::vector<std::thread> tab_threads;
+ tab_threads.reserve(nb_threads);
+
+ std::mutex queue_mutex;
+ std::queue<std::size_t> strips_to_process;
+ for (std::size_t i=0; i<nb_strips; ++i)
+ { strips_to_process.push(i); }
+
+ local_timer.reset();
+ local_timer.start();
+ for (std::size_t i=0; i<nb_threads; ++i)
+ {
+ tab_threads.push_back
+ (std::thread
+ ([i, nb_threads, &queue_mutex, &strips_to_process, &all_segments,
+ &xpos_strips, crop,this]()
+ {
+ std::size_t nb_intersections=0;
+ std::ostringstream ss;
+ My_timer timer; timer.start();
+ std::size_t strip_id;
+ queue_mutex.lock();
+ while(!strips_to_process.empty())
+ {
+ strip_id=strips_to_process.front();
+ strips_to_process.pop();
+ queue_mutex.unlock();
+ ss<<strip_id<<" ";
+
+ m_tab_partial_hds[strip_id]=new Partial_hds
+ (xpos_strips[strip_id], m_ymin, xpos_strips[strip_id+1], m_ymax,
+ strip_id, m_nb_strips);
+
+ typedef My_visitor<My_arr_segment_traits_2<EPECK>> Visitor;
+ My_arr_segment_traits_2<EPECK> m_traits;
+ Visitor visitor(*m_tab_partial_hds[strip_id]);
+ My_surface_sweep_2<Visitor> surface_sweep(&m_traits, &visitor);
+ visitor.sweep(all_segments, crop);
+ nb_intersections+=visitor.number_of_intersections();
+ assert(m_tab_partial_hds[strip_id]->is_valid_external_halfedges());
+ queue_mutex.lock();
+ }
+ queue_mutex.unlock();
+ timer.stop();
+ print_txt(" [Thread ", i, "]: strips=", ss.str(), "; #intersections=",
+ nb_intersections, " time=",timer.time());
+ }
+ ));
+ }
+ for (std::size_t i=0; i<tab_threads.size(); ++i)
+ { tab_threads[i].join(); }
+ std::cout<<std::endl;
+
+ local_timer.stop();
+ times[2]+=local_timer.time();
+
+ local_timer.reset();
+ local_timer.start();
+ build_faces_array();
+ local_timer.stop();
+ times[3]+=local_timer.time();
+
+ global_timer.stop();
+ times[0]+=global_timer.time();
+ }
+
+ ~SHds()
+ {
+ for (auto it=m_tab_partial_hds.begin(),
+ itend=m_tab_partial_hds.end(); it!=itend; ++it)
+ { delete *it; *it=nullptr; }
+ }
+
+ std::size_t number_of_strips() const
+ { return m_nb_strips; }
+ std::size_t number_of_faces() const
+ { return m_faces.size(); }
+ std::size_t number_of_finite_faces() const
+ { return m_nb_finite_faces; }
+
+ double xmin() const
+ { return m_xmin; }
+ double ymin() const
+ { return m_ymin; }
+ double xmax() const
+ { return m_xmax; }
+ double ymax() const
+ { return m_ymax; }
+
+ HDS& hds(std::size_t nb)
+ {
+ assert(nb<number_of_strips());
+ return m_tab_partial_hds[nb]->hds();
+ }
+ const HDS& hds(std::size_t nb) const
+ {
+ assert(nb<number_of_strips());
+ return m_tab_partial_hds[nb]->hds();
+ }
+
+ HDS& hds(const Global_halfedge& phe)
+ { return hds(phe.strip_id()); }
+ const HDS& hds(const Global_halfedge& phe) const
+ { return hds(phe.strip_id()); }
+
+ Partial_hds& partial_hds(std::size_t nb)
+ {
+ assert(nb<number_of_strips());
+ return *(m_tab_partial_hds[nb]);
+ }
+ const Partial_hds& partial_hds(std::size_t nb) const
+ {
+ assert(nb<number_of_strips());
+ return *(m_tab_partial_hds[nb]);
+ }
+
+ Partial_hds& partial_hds(const Global_halfedge& phe)
+ { return partial_hds(phe.strip_id()); }
+ const Partial_hds& partial_hds(const Global_halfedge& phe) const
+ { return partial_hds(phe.strip_id()); }
+
+ void clear()
+ {
+ m_xmin=m_ymin=m_xmax=m_ymax=0;
+ m_nb_strips=m_number_segments=m_number_of_halfedges=0;
+ for (auto it=m_tab_partial_hds.begin(),
+ itend=m_tab_partial_hds.end(); it!=itend; ++it)
+ { delete *it; *it=nullptr; }
+ m_faces.clear();
+ init_all_marks();
+ }
+
+ void read_segments_and_compute_bbox(const std::string& filename,
+ std::vector<ISegment>& all_segments)
+ {
+ all_segments.clear();
+ all_segments.reserve(number_of_lines(filename));
+
+ CGAL::Bbox_2 bbox;
+ m_number_segments=0;
+ Read_from_file<EPICK> reader(filename);
+ while (reader.ok())
+ {
+ EPICK::Segment_2 s=reader.next_segment();
+ all_segments.push_back(s);
+ if (m_number_segments==0) { bbox=s.bbox(); }
+ else { bbox+=s.bbox(); }
+ ++m_number_segments;
+ }
+ if (m_number_segments==0)
+ {
+ std::cout<<"Empty set of segments for file "<<filename<<", abort."<<std::endl;
+ exit(EXIT_FAILURE);
+ }
+
+ double deltax=(bbox.xmax()-bbox.xmin())/10000.; // 0.01 %
+ double deltay=(bbox.ymax()-bbox.ymin())/10000.; // 0.01 %
+ m_xmin=bbox.xmin()-deltax;
+ m_ymin=bbox.ymin()-deltay;
+ m_xmax=bbox.xmax()+deltax;
+ m_ymax=bbox.ymax()+deltay;
+ }
+
+ void compute_strips_v0(const std::string& filename,
+ std::vector<ISegment>& all_segments,
+ std::vector<double>& xpos_strips,
+ std::vector<double>* xpos_cuts)
+ {
+ if (xpos_cuts==nullptr || xpos_cuts->size()!=number_of_strips()-1)
+ {
+ std::cout<<"Error in compute_strips_v0, switch to default strip computation."<<std::endl;
+ compute_strips_v1(filename, all_segments, xpos_strips);
+ return;
+ }
+
+ read_segments_and_compute_bbox(filename, all_segments);
+
+ xpos_strips.clear();
+ xpos_strips.resize(number_of_strips()+1, 0);
+
+ xpos_strips[0]=m_xmin;
+ for (std::size_t i=0; i<xpos_cuts->size(); ++i)
+ { xpos_strips[i+1]=xpos_cuts->at(i); }
+ xpos_strips[number_of_strips()]=m_xmax;
+ }
+
+ void compute_strips_v1(const std::string& filename,
+ std::vector<ISegment>& all_segments,
+ std::vector<double>& xpos_strips)
+ {
+ read_segments_and_compute_bbox(filename, all_segments);
+
+ xpos_strips.clear();
+ xpos_strips.resize(number_of_strips()+1, 0);
+
+ xpos_strips[0]=m_xmin;
+ for (std::size_t i=1; i<number_of_strips(); ++i)
+ { xpos_strips[i]=m_xmin+((m_xmax-m_xmin)*i)/number_of_strips(); }
+ xpos_strips[number_of_strips()]=m_xmax;
+ }
+
+ void compute_strips_v2(const std::string& filename,
+ std::vector<ISegment>& all_segments,
+ std::vector<double>& xpos_strips)
+ {
+ constexpr std::size_t array_size=1000;
+
+ read_segments_and_compute_bbox(filename, all_segments);
+
+ // Now compute histogram of segments (only x coordinates)
+ std::vector<std::size_t> nbsegments(array_size,0);
+ std::size_t sum=0;
+ for (std::size_t i=0; i<all_segments.size(); ++i)
+ {
+ ISegment& s=all_segments[i];
+ std::size_t imin=floor((array_size*(s.source().x()-m_xmin))/(m_xmax-m_xmin));
+ std::size_t imax=floor((array_size*(s.target().x()-m_xmin))/(m_xmax-m_xmin));
+ for (; imin<=imax; ++imin)
+ {
+ assert(imin<array_size);
+ ++nbsegments[imin]; ++sum;
+ }
+ }
+
+ // Now we need to dispatch sum parts of segment into number_of_threads strips.
+ std::size_t cursum=0;
+ xpos_strips.clear();
+ xpos_strips.reserve(number_of_strips()+1);
+ xpos_strips.push_back(m_xmin);
+ for (std::size_t i=0; i<array_size && xpos_strips.size()<number_of_strips(); ++i)
+ {
+ if (cursum+nbsegments[i]>sum/number_of_strips())
+ {
+ xpos_strips.push_back(m_xmin+i*(m_xmax-m_xmin)/array_size);
+ cursum=0;
+ }
+ cursum+=nbsegments[i];
+ }
+ xpos_strips.push_back(m_xmax);
+ assert(xpos_strips.size()==number_of_strips()+1);
+ }
+
+ // phe is a left_to_right external halfedge.
+ // Move it in the halfedge with the same id in the left strip.
+ void move_to_same_halfedge_in_left_strip(Global_halfedge& phe) const
+ {
+ assert(partial_hds(phe).is_external(phe.halfedge()));
+ assert(hds(phe).left_to_right(phe.halfedge()));
+ phe.go_to_left_strip();
+ auto find_external=partial_hds(phe).m_right_external_halfedges.find(hds(phe).edge_id(phe.halfedge()));
+ assert(find_external!=partial_hds(phe).m_right_external_halfedges.end());
+ phe.halfedge()=hds(phe).opposite(find_external->second);
+ }
+
+ // phe is a right_to_left external halfedge.
+ // Move it in the halfedge with the same id in the right strip.
+ void move_to_same_halfedge_in_right_strip(Global_halfedge& phe) const
+ {
+ assert(partial_hds(phe).is_external(phe.halfedge()));
+ assert(hds(phe).right_to_left(phe.halfedge()));
+ phe.go_to_right_strip();
+ auto find_external=partial_hds(phe).m_left_external_halfedges.find(-hds(phe).edge_id(phe.halfedge()));
+ assert(find_external!=partial_hds(phe).m_left_external_halfedges.end());
+ phe.halfedge()=hds(phe).opposite(find_external->second);
+ }
+
+ void move_to_same_halfedge_in_adjacent_strip(Global_halfedge& phe) const
+ {
+ assert(partial_hds(phe).is_external(phe.halfedge()));
+ if (hds(phe).left_to_right(phe.halfedge()))
+ { move_to_same_halfedge_in_left_strip(phe); }
+ else
+ { move_to_same_halfedge_in_right_strip(phe); }
+ }
+
+ // phe is a parallel halfedge. If this is a external halfedge, move to the
+ // corresponding non external halfedge.
+ void move_to_non_external(Global_halfedge& phe) const
+ {
+ if (partial_hds(phe).is_external(phe.halfedge()))
+ {
+ if (hds(phe).left_to_right(phe.halfedge()))
+ {
+ do
+ { move_to_same_halfedge_in_left_strip(phe); }
+ while(partial_hds(phe).is_external(phe.halfedge()));
+ }
+ else
+ {
+ do
+ { move_to_same_halfedge_in_right_strip(phe); }
+ while(partial_hds(phe).is_external(phe.halfedge()));
+ }
+ }
+ }
+
+ // The following move_to_xxx methods traverse the shds jumping over external
+ // halfedges. This allows to see the stripped hds like a normal hds.
+ void move_to_previous(Global_halfedge& phe) const
+ {
+ assert(!partial_hds(phe).is_external(phe.halfedge()));
+ hds(phe.strip_id()).move_to_prev_around_vertex(phe.halfedge());
+ move_to_opposite(phe);
+ }
+ void move_to_next(Global_halfedge& phe) const
+ {
+ assert(!partial_hds(phe).is_external(phe.halfedge()));
+ move_to_opposite(phe);
+ hds(phe.strip_id()).move_to_next_around_vertex(phe.halfedge());
+ }
+ void move_to_opposite(Global_halfedge& phe) const
+ {
+ assert(!partial_hds(phe).is_external(phe.halfedge()));
+ hds(phe.strip_id()).move_to_opposite(phe.halfedge());
+ move_to_non_external(phe);
+ }
+
+ Global_halfedge next(const Global_halfedge& phe) const
+ {
+ Global_halfedge res(phe);
+ move_to_next(res);
+ return res;
+ }
+
+ Global_halfedge previous(const Global_halfedge& phe) const
+ {
+ Global_halfedge res(phe);
+ move_to_previous(res);
+ return res;
+ }
+
+ Global_halfedge opposite(const Global_halfedge& phe) const
+ {
+ Global_halfedge res(phe);
+ move_to_opposite(res);
+ return res;
+ }
+
+ // The following local_move_to_xxx methods traverse the hds by considering
+ // also external halfedges (contrary to methods above).
+ void local_move_to_previous(Global_halfedge& phe) const
+ {
+ if (partial_hds(phe).is_external(phe.halfedge()))
+ { move_to_same_halfedge_in_adjacent_strip(phe); }
+ else
+ {
+ hds(phe.strip_id()).move_to_prev_around_vertex(phe.halfedge());
+ hds(phe.strip_id()).move_to_opposite(phe.halfedge());
+ }
+ }
+ void local_move_to_next(Global_halfedge& phe) const
+ {
+ hds(phe.strip_id()).move_to_opposite(phe.halfedge());
+ if(partial_hds(phe).is_external(phe.halfedge()))
+ {
+ move_to_same_halfedge_in_adjacent_strip(phe);
+ hds(phe.strip_id()).move_to_opposite(phe.halfedge());
+ }
+ else
+ {
+ hds(phe.strip_id()).move_to_next_around_vertex(phe.halfedge());
+ }
+ }
+ void local_move_to_opposite(Global_halfedge& phe) const
+ {
+ phe.m_hh=hds(phe.strip_id()).opposite(phe.halfedge());
+ }
+ Global_halfedge local_next(const Global_halfedge& phe) const
+ {
+ Global_halfedge res(phe);
+ local_move_to_next(res);
+ return res;
+ }
+
+ Global_halfedge local_previous(const Global_halfedge& phe) const
+ {
+ Global_halfedge res(phe);
+ local_move_to_previous(res);
+ return res;
+ }
+
+ Global_halfedge local_opposite(const Global_halfedge& phe) const
+ {
+ Global_halfedge res(phe);
+ local_move_to_opposite(res);
+ return res;
+ }
+
+ const EPoint& point(const Global_halfedge& phe) const
+ { return hds(phe.strip_id()).point(phe.halfedge()); }
+
+ Vertex_handle vertex(Halfedge_handle hh)
+ { return hh->m_vertex; }
+ Vertex_handle vertex(const Global_halfedge& phh)
+ { return vertex(phh.halfedge()); }
+
+ Face_handle create_face()
+ { return m_faces.emplace(); }
+ Face_handle face(Halfedge_handle hh)
+ { return hh->m_face; }
+ Face_handle face(const Global_halfedge& phh)
+ { return face(phh.halfedge()); }
+ void set_face(Halfedge_handle hh, Face_handle fh)
+ { hh->m_face=fh; }
+ const CGAL::Compact_container<Face>& faces() const
+ { return m_faces; }
+
+ std::size_t number_of_halfedges() const
+ {
+ std::size_t res=0;
+ for (std::size_t i=0; i<number_of_strips(); ++i)
+ { res+=hds(i).number_of_halfedges(); }
+ return res;
+ }
+
+ std::size_t number_of_non_external_halfedges() const
+ {
+ std::size_t res=0;
+ for (std::size_t i=0; i<number_of_strips(); ++i)
+ { res+=Partial_hds(i).number_of_non_externals(); }
+ return res;
+ }
+
+ std::size_t number_of_external_halfedges() const
+ {
+ std::size_t res=0;
+ for (std::size_t i=0; i<number_of_strips(); ++i)
+ { res+=m_tab_partial_hds[i]->number_of_externals(); }
+ return res;
+ }
+
+ /** Count the number of used marks.
+ * @return the number of used marks.
+ */
+ size_type number_of_used_marks() const
+ { return mnb_used_marks; }
+
+ /** Test if a given mark is reserved.
+ * @return true iff the mark is reserved (ie in used).
+ */
+ bool is_reserved(size_type amark) const
+ {
+ CGAL_assertion(amark<NB_MARKS);
+ return (mnb_times_reserved_marks[amark]!=0);
+ }
+
+ /** Count the number of marked halfedges for a given mark.
+ * @param amark the mark index.
+ * @return the number of marked halfedges for amark.
+ */
+ size_type number_of_marked_halfedges(size_type amark) const
+ {
+ CGAL_assertion( is_reserved(amark) );
+ return mnb_marked_halfedges[amark];
+ }
+
+ /** Count the number of unmarked halfedges for a given mark.
+ * @param amark the mark index.
+ * @return the number of unmarked halfedges for amark.
+ */
+ size_type number_of_unmarked_halfedges(size_type amark) const
+ {
+ CGAL_assertion( is_reserved(amark) );
+ return number_of_halfedges()-number_of_marked_halfedges(amark);
+ }
+
+ /** Test if all the halfedges are unmarked for a given mark.
+ * @param amark the mark index.
+ * @return true iff all the halfedges are unmarked for amark.
+ */
+ bool is_whole_hds_unmarked(size_type amark) const
+ { return number_of_marked_halfedges(amark)==0; }
+
+ /** Test if all the halfedges are marked for a given mark.
+ * @param amark the mark index.
+ * @return true iff all the halfedges are marked for amark.
+ */
+ bool is_whole_hds_marked(size_type amark) const
+ { return number_of_marked_halfedges(amark)==number_of_halfedges(); }
+
+ /** Reserve a new mark.
+ * Get a new free mark and return its index.
+ * All the halfedges are unmarked for this mark.
+ * @return the index of the new mark.
+ * @pre mnb_used_marks < NB_MARKS
+ */
+ size_type get_new_mark() const
+ {
+ if (mnb_used_marks == NB_MARKS)
+ {
+ std::cerr << "Not enough Boolean marks: "
+ "increase NB_MARKS." << std::endl;
+ std::cerr << " (exception launched)" << std::endl;
+ throw Exception_no_more_available_mark();
+ }
+
+ size_type m = mfree_marks_stack[mnb_used_marks];
+ mused_marks_stack[mnb_used_marks] = m;
+
+ mindex_marks[m] = mnb_used_marks;
+ mnb_times_reserved_marks[m]=1;
+
+ ++mnb_used_marks;
+ CGAL_assertion(is_whole_hds_unmarked(m));
+
+ return m;
+ }
+
+ /** Increase the number of times a mark is reserved.
+ * @param amark the mark to share.
+ */
+ void share_a_mark(size_type amark) const
+ {
+ CGAL_assertion( is_reserved(amark) );
+ ++mnb_times_reserved_marks[amark];
+ }
+
+ /** @return the number of times a mark is reserved.
+ * @param amark the mark to share.
+ */
+ size_type get_number_of_times_mark_reserved(size_type amark) const
+ {
+ CGAL_assertion( amark<NB_MARKS );
+ return mnb_times_reserved_marks[amark];
+ }
+
+ /** Negate the mark of all the halfedges for a given mark.
+ * After this call, all the marked halfedges become unmarked and all the
+ * unmarked halfedges become marked (in constant time operation).
+ * @param amark the mark index
+ */
+ void negate_mark(size_type amark) const
+ {
+ CGAL_assertion( is_reserved(amark) );
+
+ mnb_marked_halfedges[amark] = number_of_halfedges()-mnb_marked_halfedges[amark];
+
+ mmask_marks.flip(amark);
+ }
+
+ /// Return the mark value of he a given mark number.
+ bool get_he_mark(Halfedge_const_handle he, size_type amark) const
+ { return he->get_mark(amark); }
+
+ /// Set the mark of a given mark number to a given value.
+ void set_he_mark(Halfedge_const_handle he, size_type amark, bool avalue) const
+ { he->set_mark(amark, avalue); }
+
+ /// Flip the mark of a given mark number to a given value.
+ void flip_he_mark(Halfedge_const_handle he, size_type amark) const
+ { he->flip_mark(amark); }
+
+ std::bitset<NB_MARKS> get_he_marks(Halfedge_const_handle he) const
+ { return he->get_marks(); }
+ void set_he_marks(Halfedge_const_handle he,
+ const std::bitset<NB_MARKS>& amarks) const
+ { he->set_marks(amarks); }
+
+ /** Test if a given halfedge is marked for a given mark.
+ * @param he the halfedge to test.
+ * @param amark the given mark.
+ * @return true iff halfedge is marked for the mark amark.
+ */
+ bool is_marked(Halfedge_const_handle he, size_type amark) const
+ {
+ CGAL_assertion( is_reserved(amark) );
+
+ return get_he_mark(he, amark)!=mmask_marks[amark];
+ }
+
+ /** Set the mark of a given halfedge to a state (on or off).
+ * @param ahalfedge the halfedge.
+ * @param amark the given mark.
+ * @param astate the state of the mark (on or off).
+ */
+ void set_mark_to(Halfedge_const_handle he, size_type amark,
+ bool astate) const
+ {
+ CGAL_assertion(is_reserved(amark));
+
+ if (is_marked(he, amark)!=astate)
+ {
+ if (astate) ++mnb_marked_halfedges[amark];
+ else --mnb_marked_halfedges[amark];
+
+ flip_he_mark(he, amark);
+ }
+ }
+
+ /** Mark the given halfedge.
+ * @param he the halfedge.
+ * @param amark the given mark.
+ */
+ void mark(Halfedge_const_handle he, size_type amark) const
+ {
+ CGAL_assertion(is_reserved(amark));
+
+ if (is_marked(he, amark)) return;
+
+ ++mnb_marked_halfedges[amark];
+ flip_he_mark(he, amark);
+ }
+
+ /** Unmark the given halfedge.
+ * @param he the halfedge.
+ * @param amark the given mark.
+ */
+ void unmark(Halfedge_const_handle he, size_type amark) const
+ {
+ CGAL_assertion(is_reserved(amark));
+
+ if (!is_marked(he, amark)) return;
+
+ --mnb_marked_halfedges[amark];
+ flip_he_mark(he, amark);
+ }
+
+ bool is_marked(const Global_halfedge& phe, size_type amark) const
+ { return is_marked(phe.halfedge(), amark); }
+ void mark(const Global_halfedge& phe, size_type amark) const
+ { mark(phe.halfedge(), amark); }
+ void unmark(const Global_halfedge& phe, size_type amark) const
+ { unmark(phe.halfedge(), amark); }
+
+ /** Unmark all the halfedges of the hds for a given mark.
+ * If all the halfedges are marked or unmarked, this operation takes O(1)
+ * operations, otherwise it traverses all the halfedges of the map.
+ * @param amark the given mark.
+ */
+ void unmark_all(size_type amark) const
+ {
+ CGAL_assertion( is_reserved(amark) );
+
+ if ( is_whole_hds_marked(amark) )
+ { negate_mark(amark); }
+ else if ( !is_whole_hds_unmarked(amark) )
+ {
+ for (std::size_t i=0; i<number_of_strips(); ++i)
+ {
+ for (auto it(hds(i).halfedges().begin()), itend(hds(i).halfedges().end());
+ it!=itend; ++it)
+ unmark(it, amark);
+ }
+ }
+ CGAL_assertion(is_whole_hds_unmarked(amark));
+ }
+
+ /** Free a given mark, previously calling unmark_all.
+ * @param amark the given mark.
+ */
+ void free_mark(size_type amark) const
+ {
+ CGAL_assertion( is_reserved(amark) );
+
+ if ( mnb_times_reserved_marks[amark]>1 )
+ {
+ --mnb_times_reserved_marks[amark];
+ return;
+ }
+
+ unmark_all(amark);
+
+ // 1) We remove amark from the array mused_marks_stack by
+ // replacing it with the last mark in this array.
+ mused_marks_stack[mindex_marks[amark]] =
+ mused_marks_stack[--mnb_used_marks];
+ mindex_marks[mused_marks_stack[mnb_used_marks]] =
+ mindex_marks[amark];
+
+ // 2) We add amark in the array mfree_marks_stack and update its index.
+ mfree_marks_stack[ mnb_used_marks ] = amark;
+ mindex_marks[amark] = mnb_used_marks;
+
+ mnb_times_reserved_marks[amark]=0;
+ }
+
+ // Build the faces of the given lcc (and the father/son relations between holes)
+ // Return the number of finite faces.
+ std::size_t build_faces_array()
+ {
+ // display_info();
+
+ Face_handle m_root=create_face();
+ m_root->m_finite=false; // Infinite face
+
+ std::size_t finite_faces=0;
+ auto treated=get_new_mark();
+ Face_handle newf, father;
+ Global_halfedge dh, initdh, odh;
+ Halfedge_handle firsthalfedge, fatherhalfedge;
+ std::stack<Global_halfedge> totreat;
+
+ // Iterate through all inner points, ordered from left to right
+ for (std::size_t i=0; i<number_of_strips(); ++i)
+ {
+ for (auto it=hds(i).vertices().begin(),
+ itend=hds(i).vertices().end(); it!=itend; ++it)
+ {
+ if (!Partial_hds(i).is_external(it->first_halfedge()) &&
+ !is_marked(it->first_halfedge(), treated))
+ { // first_halfedge belongs necessarily to one infinite face.
+ firsthalfedge=it->first_halfedge();
+ newf=create_face();
+
+ // std::size_t nbhalfedges=0;
+ initdh.set(i, firsthalfedge);
+ dh=initdh;
+ do
+ {
+ // ++nbhalfedges;
+ mark(dh, treated);
+ set_face(dh.halfedge(), newf);
+ odh=local_opposite(dh);
+ if (!is_marked(odh, treated))
+ { totreat.push(odh); }
+ local_move_to_next(dh);
+ }
+ while(dh!=initdh);
+
+ fatherhalfedge=it->above_halfedge();
+ father=nullptr;
+ if (fatherhalfedge==nullptr)
+ { father=m_root; }
+ else
+ {
+ assert(face(fatherhalfedge)!=nullptr);
+
+ if (face(fatherhalfedge)->m_finite)
+ { // Father is the halfedge of the finite face containing this hole
+ father=face(fatherhalfedge);
+ }
+ else
+ { // Father is a halfedge of another hole, inside the same finite face
+ // This hole was already considered before.
+ assert(face(fatherhalfedge)->m_father!=nullptr);
+ father=face(fatherhalfedge)->m_father;
+ }
+ }
+ assert(father==m_root || father->m_finite);
+ newf->set(i, firsthalfedge, father, false); // Infinite face
+ father->m_son.push_back(newf);
+
+ // Now we iterate through the cc containing it->first_halfedge().
+ // Each new face is necessarily a finite one.
+ while(!totreat.empty())
+ {
+ initdh=totreat.top(); totreat.pop();
+ if (!is_marked(initdh, treated))
+ {
+ newf=create_face();
+ newf->set(initdh.strip_id(), initdh.halfedge(), nullptr, true); // Finite face
+ ++finite_faces;
+
+ // nbhalfedges=0;
+ dh=initdh;
+ do
+ {
+ // ++nbhalfedges;
+ mark(dh, treated);
+ set_face(dh.halfedge(), newf);
+ odh=local_opposite(dh);
+ if (!is_marked(odh, treated))
+ { totreat.push(odh); }
+ local_move_to_next(dh);
+ }
+ while(dh!=initdh);
+ }
+ }
+ }
+ }
+ }
+ assert(is_whole_hds_marked(treated));
+ free_mark(treated);
+
+ m_nb_finite_faces=finite_faces;
+ return finite_faces;
+ }
+
+ std::size_t traverse_all_faces() const
+ {
+ std::size_t nb=0;
+ for (auto itf=faces().begin(), itfend=faces().end(); itf!=itfend; ++itf)
+ {
+ if (itf->finite() || itf==faces().begin())
+ {
+ if (itf!=faces().begin())
+ { // Because the root face does not have any first_halfedge
+ Global_halfedge dhinit(itf->strip_id(), itf->first_halfedge());
+ Global_halfedge dh=dhinit;
+ do
+ {
+ ++nb;
+ move_to_next(dh);
+ }
+ while(dh!=dhinit);
+ }
+ for (auto its=itf->sons().begin(), itsend=itf->sons().end();
+ its!=itsend; ++its)
+ {
+ Global_halfedge dhinit((*its)->strip_id(), (*its)->first_halfedge());
+ Global_halfedge dh=dhinit;
+ do
+ {
+ ++nb;
+ move_to_next(dh);
+ }
+ while(dh!=dhinit);
+ }
+ }
+ }
+ return nb;
+ }
+
+ void display_info() const
+ {
+ std::size_t nbexternal=0, nbleft=0, nbright=0, nbhalfedges=0, nbvertices=0, sumsegments=0;
+
+ std::cout<<"#halfedges, #vertices #input-segments (#left-external, #right-externals) for each strip: "<<std::endl;
+ for (std::size_t i=0; i<number_of_strips(); ++i)
+ {
+ std::cout<<m_tab_partial_hds[i]->hds().number_of_halfedges()
+ <<" "<<m_tab_partial_hds[i]->hds().number_of_vertices()
+ <<" "<<m_tab_partial_hds[i]->number_of_input_segments()
+ <<" ("
+ <<m_tab_partial_hds[i]->number_of_left_externals()
+ <<", "
+ <<m_tab_partial_hds[i]->number_of_right_externals()
+ <<") | ";
+ nbhalfedges+=m_tab_partial_hds[i]->hds().number_of_halfedges();
+ nbvertices+=m_tab_partial_hds[i]->hds().number_of_vertices();
+ nbleft+=m_tab_partial_hds[i]->m_left_external_halfedges.size();
+ nbright+=m_tab_partial_hds[i]->m_right_external_halfedges.size();
+ nbexternal+=m_tab_partial_hds[i]->number_of_externals();
+ sumsegments+=m_tab_partial_hds[i]->number_of_input_segments();
+ }
+
+ std::cout<<std::endl<<"Total: #halfedges="<<nbhalfedges
+ <<", #vertices="<<nbvertices<<", #faces="<<m_faces.size()
+ <<" (#finite-faces="<<m_nb_finite_faces
+ <<"), #nbexternal="<<nbexternal
+ <<" (left="<<nbleft<<", right="<<nbright<<")"<<std::endl;
+ std::cout<<"Nb input segments: "<<m_number_segments
+ <<", Sum segments in strips: "<<sumsegments<<std::endl;
+ std::cout<<"CURRENT RSS: "<<getCurrentRSS()<<" bytes (i.e. "<<
+ double(getCurrentRSS())/double(1024*1024)<<" mega-bytes and "<<
+ double(getCurrentRSS())/double(1024*1024*1024)<<" giga-bytes)"<<std::endl;
+ }
+
+ friend std::ostream& operator<<(std::ostream& os,
+ const SHds& pa)
+ {
+ os<<pa.m_xmin<<" "<<pa.m_ymin<<" "<<pa.m_xmax<<" "<<pa.m_ymax<<" "
+ <<pa.m_nb_strips<<" "<<pa.m_number_segments<<std::endl;
+
+ for (auto it=pa.m_tab_partial_hds.begin(),
+ itend=pa.m_tab_partial_hds.end(); it!=itend; ++it)
+ { os<<**it<<std::endl; }
+
+ return os;
+ }
+
+ friend std::istream& operator>>(std::istream& is,
+ SHds& pa)
+ {
+ pa.clear();
+
+ is>>pa.m_xmin>>pa.m_ymin>>pa.m_xmax>>pa.m_ymax
+ >>pa.m_nb_strips>>pa.m_number_segments;
+
+ pa.m_tab_partial_hds.resize(pa.m_nb_strips);
+ for (std::size_t i=0; i<pa.m_nb_strips; ++i)
+ {
+ pa.m_tab_partial_hds[i]=new Partial_hds;
+ is>>*(pa.m_tab_partial_hds[i]);
+ }
+
+ return is;
+ }
+
+ // Load the streamed files; filename contains a % which will be replaced
+ // by the number of the strip.
+ void load_streamed_files(const std::string& filename)
+ {
+ m_xmin=m_ymin=m_xmax=m_ymax=0.;
+ m_number_segments=0;
+ m_number_of_halfedges=0;
+ m_tab_partial_hds.clear();
+
+ std::size_t i=0;
+ std::string real_file=filename;
+ real_file.replace(real_file.find("%"), 1, std::to_string(i));
+ while(boost::filesystem::exists(boost::filesystem::path(real_file)))
+ {
+ m_tab_partial_hds.push_back(new Partial_hds);
+ std::ifstream is(real_file);
+ is>>*(m_tab_partial_hds.back());
+ is.close();
+
+ if (m_tab_partial_hds.size()==1)
+ {
+ m_xmin=CGAL::to_double(m_tab_partial_hds.back()->xmin());
+ m_ymin=CGAL::to_double(m_tab_partial_hds.back()->ymin());
+ m_xmax=CGAL::to_double(m_tab_partial_hds.back()->xmax());
+ m_ymax=CGAL::to_double(m_tab_partial_hds.back()->ymax());
+ }
+ else
+ {
+ m_xmin=std::min(m_xmin, CGAL::to_double(m_tab_partial_hds.back()->xmin()));
+ m_ymin=std::min(m_ymin, CGAL::to_double(m_tab_partial_hds.back()->ymin()));
+ m_xmax=std::max(m_xmax, CGAL::to_double(m_tab_partial_hds.back()->xmax()));
+ m_ymax=std::max(m_ymax, CGAL::to_double(m_tab_partial_hds.back()->ymax()));
+ }
+ m_number_segments+=(m_tab_partial_hds.back()->hds().number_of_halfedges())/2;
+
+ ++i;
+ real_file=filename;
+ real_file.replace(real_file.find("%"), 1, std::to_string(i));
+ }
+
+ m_nb_strips=m_tab_partial_hds.size();
+
+ for (std::size_t i=0; i<m_tab_partial_hds.size(); ++i)
+ {
+ m_tab_partial_hds[i]->m_strip_id=i;
+ m_tab_partial_hds[i]->m_nb_strips=m_nb_strips;
+ m_tab_partial_hds[i]->m_min=EPoint(m_tab_partial_hds[i]->m_minx,
+ m_ymin);
+ m_tab_partial_hds[i]->m_max=EPoint(m_tab_partial_hds[i]->m_maxx,
+ m_ymax);
+ }
+ }
+
+ // Export the current stripped arrangement into the given filename.
+ // If draw_external draw external edges, otherwise they are not drawn.
+ // if nbstrips!=0, draw only the first nbstrips. Otherwise draw all strips.
+ void export_to_xfig(const std::string& filename,
+ bool invert_y=true, double xfig_x_size=50000,
+ std::size_t nbstrips=0,
+ const CGAL::Color& color_segments=CGAL::Color(40,40,220),
+ const CGAL::Color& color_disks=CGAL::Color(220,40,40),
+ const CGAL::Color& color_externals=CGAL::Color(90,190,110),
+ const CGAL::Color& color_borders=CGAL::Color(150,150,150),
+ uint16_t width_segments=2, uint16_t radius_disks=30,
+ uint16_t width_externals=2, uint16_t width_borders=1,
+ uint16_t space_between_strips=100)
+ {
+ Export_xfig export_xfig(filename, (m_xmax-m_xmin),
+ xfig_x_size,
+ invert_y);
+
+ uint8_t coul_segments=export_xfig.add_color(color_segments);
+ uint8_t coul_disks=export_xfig.add_color(color_disks);
+ uint8_t coul_externals=export_xfig.add_color(color_externals);
+ uint8_t coul_borders=export_xfig.add_color(color_borders);
+
+ double dx=-m_xmin;
+ double dy=-m_ymin;
+
+ std::size_t i=0;
+ for (auto it=m_tab_partial_hds.begin(),
+ itend=m_tab_partial_hds.end();
+ (nbstrips==0 || i<nbstrips) && it!=itend; ++it, ++i)
+ {
+ export_xfig.set_x_shift(i*space_between_strips);
+
+ if (width_externals>0)
+ {
+ double x1=CGAL::to_double((*it)->xmin());
+ double y1, x2, y2;
+ for (auto itc=(*it)->m_left_external_halfedges.begin(),
+ itcend=(*it)->m_left_external_halfedges.end(); itc!=itcend; ++itc)
+ {
+ Global_halfedge h1((*it)->strip_id(), itc->second);
+ move_to_non_external(h1);
+ Global_halfedge h2=opposite(h1);
+ move_to_non_external(h2);
+
+ EPECK::Line_2 line(point(h1), point(h2));
+ y1=CGAL::to_double(line.y_at_x((*it)->xmin()));
+
+ if (h2.strip_id()==(*it)->strip_id())
+ { // Non traversing segment
+ x2=CGAL::to_double(point(h2).x());
+ y2=CGAL::to_double(point(h2).y());
+ }
+ else
+ { // Traversing segment
+ x2=CGAL::to_double((*it)->xmax());
+ y2=CGAL::to_double(line.y_at_x((*it)->xmax()));
+ }
+
+ export_xfig.segment(x1+dx, y1+dy, x2+dx, y2+dy,
+ coul_externals, width_externals, false, 100);
+ }
+
+ x1=CGAL::to_double((*it)->xmax());
+ for (auto itc=(*it)->m_right_external_halfedges.begin(),
+ itcend=(*it)->m_right_external_halfedges.end(); itc!=itcend; ++itc)
+ {
+ Global_halfedge h1((*it)->strip_id(), itc->second);
+ move_to_non_external(h1);
+ Global_halfedge h2=opposite(h1);
+ move_to_non_external(h2);
+
+ EPECK::Line_2 line(point(h1), point(h2));
+ y1=CGAL::to_double(line.y_at_x((*it)->xmax()));
+
+ if (h2.strip_id()==(*it)->strip_id())
+ { // Non traversing segment; traversing segments were alreary drawn above
+ x2=CGAL::to_double(point(h2).x());
+ y2=CGAL::to_double(point(h2).y());
+
+ export_xfig.segment(x1+dx, y1+dy, x2+dx, y2+dy,
+ coul_externals, width_externals, false, 100);
+ }
+ }
+ }
+
+ (*it)->export_to_xfig(export_xfig,
+ dx, dy,
+ coul_segments, coul_disks, coul_borders,
+ width_segments, radius_disks, width_borders);
+ }
+ }
+
+protected:
+ double m_xmin, m_ymin, m_xmax, m_ymax;
+ std::size_t m_nb_strips;
+ std::size_t m_number_segments;
+ std::size_t m_number_of_halfedges;
+ std::size_t m_nb_finite_faces;
+
+ std::vector<Partial_hds*> m_tab_partial_hds;
+ CGAL::Compact_container<Face> m_faces;
+
+ /// Number of times each mark is reserved. 0 if the mark is free.
+ mutable size_type mnb_times_reserved_marks[NB_MARKS];
+
+ /// Mask marks to know the value of unmark halfedge, for each index i.
+ mutable std::bitset<NB_MARKS> mmask_marks;
+
+ /// Number of used marks.
+ mutable size_type mnb_used_marks;
+
+ /// Index of each mark, in mfree_marks_stack or in mfree_marks_stack.
+ mutable size_type mindex_marks[NB_MARKS];
+
+ /// "Stack" of free marks.
+ mutable size_type mfree_marks_stack[NB_MARKS];
+
+ /// "Stack" of used marks.
+ mutable size_type mused_marks_stack[NB_MARKS];
+
+ /// Number of marked halfedges for each used marks.
+ mutable size_type mnb_marked_halfedges[NB_MARKS];
+};
+////////////////////////////////////////////////////////////////////////////////
+
+#endif // SHDS_H
diff --git a/src/SHds_to_cgal_arrangement.h b/src/SHds_to_cgal_arrangement.h
new file mode 100644
index 0000000000000000000000000000000000000000..f2907ee95da2bd1a4e43f6040413c7f0eb9a9237
--- /dev/null
+++ b/src/SHds_to_cgal_arrangement.h
@@ -0,0 +1,54 @@
+ // compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#ifndef SHDS_TO_CGAL_ARRANGEMENT_H
+#define SHDS_TO_CGAL_ARRANGEMENT_H
+
+#include "SHds.h"
+#include <CGAL/Arrangement_2.h>
+
+// Build a CGAL arrangement from a parallel arrangement.
+template<typename Arr>
+void shds_to_cgal_arrangement(SHds& shds, Arr& arr)
+{
+ arr.clear();
+ std::vector<ESegment> all_segments;
+ all_segments.reserve(shds.number_of_non_external_halfedges()/2);
+
+ Global_halfedge hh;
+ for (std::size_t i=0; i<shds.number_of_strips(); ++i)
+ {
+ for (auto it=shds.hds(i).halfedges().begin(),
+ itend=shds.hds(i).halfedges().end(); it!=itend; ++it)
+ {
+ if (!shds.partial_hds(i).is_external(it))
+ { // First halfedge is it
+ hh=shds.opposite(Global_halfedge(i, it)); // Opposite halfedge, a non critical dart
+ if (it<hh.halfedge()) // To add a segment only once
+ { all_segments.push_back(ESegment(shds.hds(i).point(it),
+ shds.hds(hh).point(hh.halfedge()))); }
+ }
+ }
+ }
+
+ CGAL::insert_non_intersecting_curves(arr,
+ all_segments.begin(), all_segments.end());
+}
+
+#endif // SHDS_TO_CGAL_ARRANGEMENT_H
diff --git a/src/SHds_to_lcc.h b/src/SHds_to_lcc.h
new file mode 100644
index 0000000000000000000000000000000000000000..12fbb7e1df01e9d094b348ce299426a4e01e4b3b
--- /dev/null
+++ b/src/SHds_to_lcc.h
@@ -0,0 +1,259 @@
+// compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#ifndef PARALLEL_ARRANGEMENT_TO_LCC_H
+#define PARALLEL_ARRANGEMENT_TO_LCC_H
+
+#include "SHds.h"
+#include "CGAL_typedefs.h"
+#include <tbb/parallel_for.h>
+#include <tbb/blocked_range.h>
+#include <unordered_map>
+
+template<typename LCC>
+class SHDS_to_LCC
+{
+public:
+ typedef std::unordered_map<Halfedge_handle, typename LCC::Dart_handle> Assoc;
+
+ SHDS_to_LCC(SHds& ashds, LCC& alcc): m_shds(ashds), m_lcc(alcc)
+ {}
+
+ std::size_t shds_to_lcc()
+ {
+ m_lcc.clear();
+ create_darts();
+ set_beta_and_vertices_links();
+ return build_faces_array();
+ }
+
+ typename LCC::Dart_handle halfedge_to_dart(Halfedge_handle hh)
+ {
+ return assoc[hh];
+ }
+
+ void create_darts()
+ {
+ for (std::size_t i=0; i<m_shds.number_of_strips(); ++i)
+ {
+ for (auto it=m_shds.hds(i).halfedges().begin(),
+ itend=m_shds.hds(i).halfedges().end(); it!=itend; ++it)
+ {
+ if (!m_shds.partial_hds(i).is_external(it))
+ { assoc[it]=m_lcc.create_dart(); }
+ }
+ }
+ }
+
+ void set_beta_and_vertices_links()
+ {
+ tbb::parallel_for(tbb::blocked_range<std::size_t>(0, m_shds.number_of_strips()),
+ [&](tbb::blocked_range<std::size_t> r)
+ {
+ for (std::size_t i=r.begin(); i<r.end(); ++i)
+ {
+ Global_halfedge hh;
+ for (auto it=m_shds.hds(i).halfedges().begin(),
+ itend=m_shds.hds(i).halfedges().end(); it!=itend; ++it)
+ {
+ if (!m_shds.partial_hds(i).is_external(it))
+ { // First halfedge is it
+ hh=m_shds.opposite(Global_halfedge(i, it)); // Opposite halfedge, non external
+ if (it<hh.halfedge())
+ { m_lcc.template basic_link_beta_for_involution<2>
+ (halfedge_to_dart(it),
+ halfedge_to_dart(hh.halfedge())); }
+ m_lcc.template basic_link_beta_1
+ (halfedge_to_dart(it),
+ halfedge_to_dart(hh.halfedge()->next_around_vertex()));
+ }
+ }
+ }
+ }
+ );
+
+
+ // Last step not in parallel to keep vertices in the same order.
+ typename LCC::Dart_handle above=nullptr;
+ for (std::size_t i=0; i<m_shds.number_of_strips(); ++i)
+ {
+ Global_halfedge hh;
+ for (auto it=m_shds.hds(i).vertices().begin(),
+ itend=m_shds.hds(i).vertices().end(); it!=itend; ++it)
+ {
+ above=nullptr;
+ if (it->above_halfedge()!=nullptr)
+ {
+ hh=Global_halfedge(i, it->above_halfedge());
+ m_shds.move_to_non_external(hh);
+ above=halfedge_to_dart(hh.halfedge());
+ assert(above!=nullptr);
+ }
+ m_lcc.set_vertex_attribute(halfedge_to_dart(it->first_halfedge()),
+ m_lcc.create_vertex_attribute(m_shds.hds(i).point(it), above));
+ }
+ }
+ }
+
+ // Build the faces of the given lcc (and the father/son relations between holes)
+ // Return the number of finite faces. TODO copy the face array of m_shds
+ // instead of re-compute it.
+ std::size_t build_faces_array()
+ {
+ // display_info();
+
+ auto m_root=m_lcc.template create_attribute<2>();
+ m_root->info().m_finite=true; // Infinite face but considered as finite
+ // to ensure the property: each infinite face is included into a finite one.
+
+ std::size_t finite_faces=0;
+ auto treated=m_lcc.get_new_mark();
+ typename LCC::template Attribute_handle<2>::type newf, father;
+ Dart_handle dh, initdh, firstdart, fatherdart;
+ std::stack<Dart_handle> totreat;
+
+ // Iterate through all inner points, ordered from left to right
+ for (std::size_t i=0; i<m_shds.number_of_strips(); ++i)
+ {
+ for (auto it=m_shds.hds(i).vertices().begin(),
+ itend=m_shds.hds(i).vertices().end(); it!=itend; ++it)
+ {
+ if (!m_shds.partial_hds(i).is_external(it->first_halfedge()) &&
+ !m_lcc.is_marked(halfedge_to_dart(it->first_halfedge()),
+ treated))
+ { // first_halfedge belongs necessarily to one infinite face.
+ firstdart=halfedge_to_dart(it->first_halfedge());
+ assert(m_lcc.is_dart_used(firstdart));
+
+ newf=m_lcc.template create_attribute<2>();
+ newf->info().m_finite=false; // Infinite face
+
+ // std::size_t nbdarts=0;
+ initdh=dh=firstdart;
+ do
+ {
+ // ++nbdarts;
+ m_lcc.mark(dh, treated);
+ m_lcc.template set_dart_attribute<2>(dh, newf);
+ if (!m_lcc.is_marked(m_lcc.template beta<2>(dh), treated))
+ { totreat.push(m_lcc.template beta<2>(dh)); }
+ dh=m_lcc.template beta<1>(dh);
+ }
+ while(dh!=initdh);
+
+ /* print_txt_with_endl("[Infinite Face ",
+ m_lcc.template attributes<2>().index(newf), "]: ",
+ "first segment: ",
+ m_lcc.point(firstdart), " , ",
+ m_lcc.point(m_lcc.other_extremity(firstdart)),
+ "; #darts=", nbdarts); */
+ fatherdart=m_lcc.template info<0>(firstdart).m_dart_above;
+ father=nullptr;
+ if (fatherdart==nullptr)
+ { father=m_root; }
+ else
+ {
+ assert(m_lcc.is_dart_used(fatherdart));
+ assert(m_lcc.template attribute<2>(fatherdart)!=nullptr);
+ assert(m_lcc.point(fatherdart)<
+ m_lcc.point(m_lcc.other_extremity(fatherdart)));
+
+ /*print_txt_with_endl("father_of_first_halfedge: ",
+ m_lcc.point(fatherdart), " , ",
+ m_lcc.point(m_lcc.other_extremity(fatherdart))); */
+
+ if (m_lcc.template info<2>(fatherdart).m_finite)
+ { // Father is the dart of the finite face containing this hole
+ father=m_lcc.template attribute<2>(fatherdart);
+ }
+ else
+ { // Father is a dart of another hole, inside the same finite face
+ // This hole was already considered before.
+ assert(m_lcc.template info<2>(fatherdart).m_father!=nullptr);
+ father=m_lcc.template info<2>(fatherdart).m_father;
+ }
+ }
+ /* print_txt_with_endl(" Father of face ",
+ m_lcc.template attributes<2>().index(newf),
+ " -> ", m_lcc.template attributes<2>().index(father));
+ */
+ assert(father->info().m_finite);
+ newf->info().m_father=father;
+ father->info().m_son.push_back(newf);
+
+ // Now we iterate through the cc containing it->first_halfedge().
+ // Each new face is necessarily a finite one.
+ while(!totreat.empty())
+ {
+ initdh=totreat.top(); totreat.pop();
+ if (!m_lcc.is_marked(initdh, treated))
+ {
+ newf=m_lcc.template create_attribute<2>();
+ newf->info().m_finite=true; // Finite face
+ ++finite_faces;
+
+ // nbdarts=0;
+ dh=initdh;
+ do
+ {
+ // ++nbdarts;
+ m_lcc.mark(dh, treated);
+ m_lcc.template set_dart_attribute<2>(dh, newf);
+ if (!m_lcc.is_marked(m_lcc.template beta<2>(dh), treated))
+ { totreat.push(m_lcc.template beta<2>(dh)); }
+ dh=m_lcc.template beta<1>(dh);
+ }
+ while(dh!=initdh);
+
+ /* print_txt_with_endl("[Finite Face ",
+ m_lcc.template attributes<2>().index(newf), "]: ",
+ "father: ",
+ m_lcc.template attributes<2>().index(father),
+ " first segment: ",
+ m_lcc.point(initdh), " , ",
+ m_lcc.point(m_lcc.other_extremity(initdh)),
+ "; #darts=", nbdarts); */
+ }
+ }
+ }
+ }
+ }
+ assert(m_lcc.is_whole_map_marked(treated));
+ m_lcc.free_mark(treated);
+
+ return finite_faces;
+ }
+
+protected:
+ SHds& m_shds;
+ LCC& m_lcc;
+ Assoc assoc;
+};
+
+
+// Build a global lcc from a parallel arrangement.
+// @Return the number of finite faces.
+template<typename LCC>
+std::size_t shds_to_lcc(SHds& ashds, LCC& alcc)
+{
+ SHDS_to_LCC<LCC> shds_to_lcc(ashds, alcc);
+ return shds_to_lcc.shds_to_lcc();
+}
+
+#endif // SHDS_TO_LCC_H
diff --git a/src/Segment_readers.h b/src/Segment_readers.h
new file mode 100644
index 0000000000000000000000000000000000000000..a6a76842ac7cbbbead3bb8f540c9b3de378235e2
--- /dev/null
+++ b/src/Segment_readers.h
@@ -0,0 +1,274 @@
+// compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#ifndef SEGMENT_READERS_H
+#define SEGMENT_READERS_H
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <CGAL/Random.h>
+#include <CGAL/point_generators_2.h>
+
+///////////////////////////////////////////////////////////////////////////////
+/// \brief The Reader class
+template<typename Kernel>
+class Reader
+{
+public:
+ using FT=typename Kernel::FT;
+ using Point=typename Kernel::Point_2;
+ using Segment=typename Kernel::Segment_2;
+
+ Reader() : m_cur_segment(0), m_split(false)
+ {}
+
+ void split(std::size_t splitLow, std::size_t splitHight)
+ {
+ m_split=true;
+ m_splitLow=splitLow;
+ m_splitHight=splitHight;
+ }
+
+ virtual ~Reader(){}
+ virtual bool reader_ok() const =0;
+ virtual Segment get_the_next_segment() =0;
+
+ bool ok() const
+ { return reader_ok() && (!m_split || m_cur_segment<m_splitHight); }
+
+ Segment next_segment()
+ {
+ Segment res;
+ do
+ {
+ res=get_the_next_segment();
+ ++m_cur_segment;
+ }
+ while(reader_ok() &&
+ (m_split &&
+ (m_cur_segment<m_splitLow || m_cur_segment>m_splitHight)));
+ return res;
+ }
+
+protected:
+ std::size_t m_cur_segment;
+ bool m_split;
+ std::size_t m_splitLow, m_splitHight;
+};
+///////////////////////////////////////////////////////////////////////////////
+template<typename Kernel>
+class Read_from_file: public Reader<Kernel>
+{
+public:
+ using Base=Reader<Kernel>;
+ using Segment=typename Base::Segment;
+ using Point=typename Base::Point;
+
+ Read_from_file(const std::string& filename): ifs(filename)
+ {
+ // std::cout<<"Segment file: "<<filename<<std::endl;
+ if (!ifs.is_open())
+ {
+ std::cerr<<"Error opening file '"<<filename<<"'"<<std::endl;
+ }
+ else if (ifs.peek()==EOF)
+ {
+ std::cerr<<"Error, empty file '"<<filename<<"'"<<std::endl;
+ ifs.close();
+ }
+ }
+
+ bool reader_ok() const
+ { return ifs.is_open(); }
+
+ Segment get_the_next_segment()
+ {
+ Segment s;
+ std::string line;
+ double x1, y1, x2, y2;
+ do
+ {
+ std::getline(ifs, line);
+ if (ifs.good())
+ { line.erase(line.find_last_not_of(" \n\r\t")+1); }
+ }
+ while (!ifs.fail() && line[0]=='#'); // jump over comment lines
+
+ if (ifs.good())
+ {
+ std::istringstream is(line);
+ is>>x1>>y1>>x2>>y2;
+ if (!is.fail())
+ {
+ if (ifs.peek()==EOF) { ifs.close(); }
+ if (x1>x2 || (x1==x2 && y1>y2))
+ { std::swap(x1, x2); std::swap(y1, y2); }
+ return Segment(Point(x1, y1), Point(x2, y2));
+ }
+ }
+ std::cerr<<"Error when reading file."<<std::endl;
+ ifs.close();
+ return Segment(Point(0,0),Point(0,0)); // An error occured.
+ }
+
+protected:
+ std::ifstream ifs;
+};
+///////////////////////////////////////////////////////////////////////////////
+template<typename Kernel>
+class Random_segments: public Reader<Kernel>
+{
+public:
+ using Base=Reader<Kernel>;
+ using Segment=typename Base::Segment;
+ using Point=typename Base::Point;
+ using FT=typename Base::FT;
+
+ Random_segments(std::size_t N, CGAL::Random& random,
+ bool size, double L) :
+ m_N(N), m_random(random), m_size(size), m_L(L), m_nb_generated(0)
+ {
+ CGAL_assertion(L>0); // Otherwise, crash of Random_points_in_disc_2 constructor
+ std::cout<<"Random seed: "<<random.get_seed();
+ if (size) std::cout<<" with size<="<<L;
+ std::cout<<std::endl;
+ }
+
+ bool reader_ok() const
+ { return m_nb_generated<m_N; }
+
+ Segment get_the_next_segment()
+ {
+ Segment s;
+ Point p2;
+ typedef CGAL::Creator_uniform_2<FT, Point> Creator;
+ CGAL::Random_points_in_disc_2<Point,Creator> rand_in_disk(m_L, m_random);
+ if (m_nb_generated==0)
+ {
+ m_p1=Point(m_random.get_double(0.0, 100.0),
+ m_random.get_double(0.0, 100.0));
+ }
+ if (!m_size)
+ { p2=Point(m_random.get_double(0.0, 100.0),
+ m_random.get_double(0.0, 100.0)); }
+ else
+ {
+ do
+ {
+ p2=(*rand_in_disk); ++rand_in_disk;
+ p2=Point(m_p1.x()+p2.x(), m_p1.y()+p2.y());
+ }
+ while(p2.x()<0.0 || p2.x()>100.0 || p2.y()<0.0 || p2.y()>100.0);
+ }
+ /* if (random.get_int(0,101)>60 && number_of_segments()>0) // dans 40% des cas
+ { // Get an existing segment
+ std::size_t existing_s=
+ static_cast<std::size_t>(random.get_int
+ (0,static_cast<int>(number_of_segments())));
+ double d1=random.get_double(0., 1.);
+ double d2=random.get_double(0., 1.);
+ typename MyKernel::Vector_2 vect
+ (get_segment(existing_s).source(), get_segment(existing_s).target());
+ p1=Point(get_segment(existing_s).source()+vect*d1);
+ p2=Point(get_segment(existing_s).target()-vect*d2);
+ }
+ else if (random.get_int(0,101)>50 && number_of_segments()>0) // dans 50% des cas restants
+ {
+ std::size_t existing_s=
+ static_cast<std::size_t>(random.get_int
+ (0,static_cast<int>(number_of_segments())));
+ p1=get_segment(existing_s).target();
+ do
+ {
+ p2=(*rand_in_disk); ++rand_in_disk;
+ p2=Point(p1.x()+p2.x(), p1.y()+p2.y());
+ }
+ while(p2.x()<0.0 || p2.x()>100.0 || p2.y()<0.0 || p2.y()>100.0);
+ }*/
+ Segment res;
+ if (m_p1.x()>p2.x() || (m_p1.x()==p2.x() && m_p1.y()>p2.y()))
+ { res=Segment(p2, m_p1); }
+ else { res=Segment(m_p1, p2); }
+
+ ++m_nb_generated;
+ m_p1=p2;
+ return res;
+ // std::cout<<"Segment "<<p1<<" -> "<<p2<<std::endl;
+ }
+
+protected:
+ std::size_t m_N;
+ CGAL::Random m_random;
+ bool m_size;
+ double m_L;
+ std::size_t m_nb_generated;
+ Point m_p1;
+};
+///////////////////////////////////////////////////////////////////////////////
+template<typename Kernel>
+class Vector_of_segments: public Reader<Kernel>
+{
+public:
+ using Base=Reader<Kernel>;
+ using Segment=typename Base::Segment;
+ using Point=typename Base::Point;
+
+ Vector_of_segments(const std::vector<Segment>& all_segments) :
+ m_all_segments(all_segments),
+ m_nb_generated(0)
+ {}
+
+ bool reader_ok() const
+ { return m_nb_generated<m_all_segments.size(); }
+
+ Segment get_the_next_segment()
+ { return m_all_segments[m_nb_generated++]; }
+
+protected:
+ const std::vector<Segment>& m_all_segments;
+ std::size_t m_nb_generated;
+};
+///////////////////////////////////////////////////////////////////////////////
+template<typename Kernel>
+inline
+void fill_segment_array(Reader<Kernel>& reader,
+ std::vector<typename Kernel::Segment_2>& all_segments)
+{
+ all_segments.clear();
+ while (reader.ok())
+ {
+ typename Kernel::Segment_2 s=reader.next_segment();
+ if (s.source()!=s.target())
+ { all_segments.push_back(s); }
+ }
+}
+///////////////////////////////////////////////////////////////////////////////
+inline std::size_t number_of_lines(const std::string& filename)
+{
+ std::ifstream in(filename);
+ std::size_t nb_lignes=0;
+ while(in.ignore(std::numeric_limits<std::streamsize>::max(), '\n'))
+ { ++nb_lignes; }
+ in.close();
+ return nb_lignes;
+}
+///////////////////////////////////////////////////////////////////////////////
+#endif // SEGMENT_READERS_H
diff --git a/src/compute-arrangement-streamed.cpp b/src/compute-arrangement-streamed.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6c14b9c7417f798cff96852dc573e64334614b55
--- /dev/null
+++ b/src/compute-arrangement-streamed.cpp
@@ -0,0 +1,248 @@
+// compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <stack>
+#include <tuple>
+#include <list>
+
+#include <CGAL/Memory_sizer.h>
+#include <CGAL/assertions.h>
+#include <CGAL/Surface_sweep_2.h>
+#include <CGAL/Surface_sweep_2_algorithms.h>
+#include <CGAL/Arr_segment_traits_2.h>
+
+#include "CGAL_typedefs.h"
+#include "Segment_readers.h"
+#include "My_visitor.h"
+#include "SHds.h"
+#include "SHds_to_lcc.h"
+#include "SHds_to_cgal_arrangement.h"
+
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////
+[[ noreturn ]] void usage(int /*argc*/, char** argv)
+{
+ // Name
+ std::cout<<"Name"<<std::endl;
+ std::cout<<" "<<argv[0]<<" - a new method to compute arrangement of segments, streamed version";
+ std::cout<<std::endl<<std::endl;
+ // Synopsis
+ std::cout<<"SYNOPSIS"<<std::endl;
+ std::cout<<" "<<argv[0]<<" [--help|-h|-?] [-t1 filename] [-nbs N]"
+ <<std::endl<<std::endl;
+ // Description
+ std::cout<<"DESCRIPTION"<<std::endl;
+ std::cout<<" "<<"Compute the arrangement of segments given in the filename."
+ <<std::endl<<std::endl;
+ // Options
+ std::cout<<" --help, -h, -?"<<std::endl
+ <<" display this help and exit."
+ <<std::endl<<std::endl;
+ std::cout<<" -t1 filename"
+ <<std::endl
+ <<" load the text segment file, and creates the arrangement: required."
+ <<std::endl<<std::endl;
+ std::cout<<" -nbs N"
+ <<std::endl
+ <<" fix the number of segments to load for one strip."
+ <<std::endl<<std::endl;
+ // Author
+ std::cout<<"AUTHOR"<<std::endl;
+ std::cout<<" "<<"Written by Guillaume Damiand."
+ <<std::endl<<std::endl;
+ // Copyright
+ std::cout<<"COPYRIGHT"<<std::endl;
+ std::cout<<" "<<"Copyright (C) 2020 CNRS and LIRIS' Establishments (France). "
+ <<"License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>"
+ <<std::endl
+ <<" This is free software: you are free to change and redistribute it under certain conditions. "
+ <<"There is NO WARRANTY, to the extent permitted by law."
+ <<std::endl<<std::endl;
+ // Reporting bugs
+ std::cout<<"REPORTING BUGS"<<std::endl;
+ std::cout<<" "<<"Email bug reports to Guillaume Damiand ⟨guillaume.damiand@liris.cnrs.fr⟩."
+ <<std::endl<<std::endl;
+ exit(EXIT_FAILURE);
+}
+////////////////////////////////////////////////////////////////////////////////
+[[ noreturn ]] void error_command_line(int argc, char** argv, const char* msg)
+{
+ std::cout<<"ERROR: "<<msg<<std::endl;
+ usage(argc, argv);
+}
+////////////////////////////////////////////////////////////////////////////////
+void process_command_line(int argc, char** argv,
+ bool& t1,
+ std::string& filename,
+ std::size_t& nbs)
+{
+ t1=false;
+ filename="";
+ nbs=10000;
+
+ bool helprequired=false;
+ std::string arg;
+
+ for (int i=1; i<argc; ++i)
+ {
+ arg=std::string(argv[i]);
+ if (arg==std::string("-h") || arg==std::string("--help") || arg==std::string("-?"))
+ { helprequired=true; }
+ else if (arg==std::string("-t1"))
+ {
+ t1=true;
+ if (argc-1-i<1)
+ { error_command_line(argc, argv, "Error: no filename after -t1 option."); }
+ filename=std::string(argv[++i]);
+ }
+ else if (arg==std::string("-nbs"))
+ {
+ if (argc-1-i<1)
+ { error_command_line(argc, argv, "Error: no number after -nbs option."); }
+ nbs=std::stoi(std::string(argv[++i]));
+ }
+ else if (arg[0]=='-')
+ {
+ std::cout<<"Unknown option "<<arg<<std::endl;
+ helprequired=true;
+ }
+ }
+
+ if (!t1) { std::cout<<"ERROR: missing file name."<<std::endl; }
+ if (helprequired || !t1) { usage(argc, argv); }
+}
+///////////////////////////////////////////////////////////////////////////////
+void process_file(const std::string& filename, std::size_t nbs)
+{
+ My_timer timer, global_timer, savetimer;
+
+ global_timer.reset();
+ global_timer.start();
+
+ std::vector<double> times(4,0);
+
+ typedef My_visitor<My_arr_segment_traits_2<EPECK>> Visitor;
+ My_arr_segment_traits_2<EPECK> m_traits;
+
+ double strip_xmin=0, strip_xmax=0, strip_ymin=0, strip_ymax=0, xmax=0;
+ bool first_strip=true;
+ std::list<ISegment> all_segments;
+ Read_from_file<EPICK> reader(filename);
+ std::size_t strip_id=0;
+ std::size_t id=1; // segment id
+ while (reader.ok())
+ {
+ // 1) Load new segments
+ timer.reset(); timer.start();
+ for (std::size_t i=0; reader.ok() && i<nbs; ++i)
+ {
+ EPICK::Segment_2 s=reader.next_segment();
+ if (first_strip && all_segments.empty())
+ { strip_xmin=s.source().x()-1.; xmax=s.target().x(); }
+ else { xmax=std::max(xmax, s.target().x()); }
+ strip_xmax=s.source().x();
+
+ // Note that we don't need ymin and ymax !
+ if (all_segments.empty())
+ { strip_ymin=s.source().y(); strip_ymax=strip_ymin;}
+ else
+ {
+ if (s.source().y()<strip_ymin) { strip_ymin=s.source().y(); }
+ if (s.source().y()>strip_ymax) { strip_ymax=s.source().y(); }
+ }
+ if (s.target().y()<strip_ymin) { strip_ymin=s.target().y(); }
+ if (s.target().y()>strip_ymax) { strip_ymax=s.target().y(); }
+
+ all_segments.push_back(s);
+ }
+ if (!reader.ok())
+ { strip_xmax=xmax+1; }
+ timer.stop();
+ times[0]+=timer.time();
+
+ // 2) Now all_segments contains all the segments of the current strip;
+ // compute the partial hds of this strip.
+ timer.reset(); timer.start();
+ {
+ Partial_hds phds(strip_xmin, strip_ymin, strip_xmax, strip_ymax,
+ (first_strip?0:(!reader.ok()?2:1)), 3);
+ Visitor visitor(phds);
+ My_surface_sweep_2<Visitor> surface_sweep(&m_traits, &visitor);
+ visitor.sweep(all_segments, false, id);
+ phds.display_info();
+
+ // 3) Save the partial hds
+ savetimer.reset(); savetimer.start();
+ std::ofstream outputfile(std::string("strip-")+std::to_string(strip_id++)+".save");
+ if (outputfile)
+ { outputfile<<phds; }
+ outputfile.close();
+ savetimer.stop();
+ times[2]+=savetimer.time();
+ }
+ timer.stop();
+ times[1]+=(timer.time()-savetimer.time()); // time to compute partial hds without
+ // including the save time; and including the de-allocation of the phds
+ // (reason of the block)
+
+ // 4) Remove all segments no more concerned by the new strip
+ timer.reset(); timer.start();
+ first_strip=false;
+ strip_xmin=strip_xmax; // xmin of the new strip is xmax of the previous strip
+ while (!all_segments.empty() &&
+ (all_segments.begin())->target().x()<strip_xmin)
+ {
+ all_segments.pop_front();
+ ++id; // id is incremented to keep global id the same
+ }
+ timer.stop();
+ times[3]+=timer.time();
+ }
+
+ global_timer.stop();
+
+ print_txt_with_endl("[TIME]: to load segments ", times[0]);
+ print_txt_with_endl("[TIME]: to compute all partial HDS ", times[1]);
+ print_txt_with_endl("[TIME]: to save partial HDS ", times[2]);
+ print_txt_with_endl("[TIME]: to update segments ", times[3]);
+ print_txt_with_endl("[GLOBAL]: time to compute all local arrangements ",
+ global_timer.time());
+ print_txt_with_endl("[GLOBAL]: time to compute all local arrangements without save/load times ",
+ times[1]+times[3]);
+}
+///////////////////////////////////////////////////////////////////////////////
+int main(int argc, char** argv)
+{
+ std::string filename;
+ bool t1;
+ std::size_t nbs;
+
+ process_command_line(argc, argv, t1, filename, nbs);
+ process_file(filename, nbs);
+
+ return EXIT_SUCCESS;
+}
+///////////////////////////////////////////////////////////////////////////////
diff --git a/src/draw_shds.h b/src/draw_shds.h
new file mode 100644
index 0000000000000000000000000000000000000000..186d121e1acaf9e0af3f077bf5c87cd2ae30911a
--- /dev/null
+++ b/src/draw_shds.h
@@ -0,0 +1,475 @@
+// compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#ifndef DRAW_SHDS_H
+#define DRAW_SHDS_H
+
+#include <CGAL/Qt/Basic_viewer_qt.h>
+
+#ifdef CGAL_USE_BASIC_VIEWER
+
+#include <CGAL/Random.h>
+#include <CGAL/Qt/Basic_viewer_qt.h>
+#include "SHds.h"
+
+// Viewer class for sHDS
+class MySimpleSHDSViewerQt : public CGAL::Basic_viewer_qt
+{
+ typedef CGAL::Basic_viewer_qt Base;
+ typedef typename Base::Local_point Point;
+
+public:
+ /// Construct the viewer.
+ /// @param ashds the sHDS to view
+ /// @param title the title of the window
+ /// @param anofaces if true, do not draw faces (faces are not computed; this can be
+ /// usefull for very big object where this time could be long)
+ MySimpleSHDSViewerQt(QWidget* parent,
+ SHds& ashds,
+ const char* title="Basic Viewer for sHDS",
+ bool anofaces=false) :
+ // First draw: vertices; edges, faces; multi-color; inverse normal
+ Base(parent, title, true, true, true, false, false),
+ m_shds(ashds),
+ m_nofaces(anofaces),
+ m_random_face_color(true),
+ m_num_face(0),
+ m_draw_grid(true),
+ m_cur_strip(0),
+ m_color_points(10, 10, 255),
+ m_color_grid(220, 220, 220),
+ m_color_selected_he(10,255,10),
+ m_color_above_links(85,210,180)
+ {
+ compute_elements();
+ }
+
+protected:
+ void draw_grid()
+ {
+ for (std::size_t i=0; i<m_shds.number_of_strips(); ++i)
+ {
+ const Partial_hds& la=m_shds.partial_hds(i);
+ add_segment(EPoint(la.xmin(), la.ymin()),
+ EPoint(la.xmax(), la.ymin()), m_color_grid);
+
+ add_segment(EPoint(la.xmin(), la.ymax()),
+ EPoint(la.xmax(), la.ymax()), m_color_grid);
+
+ add_segment(EPoint(la.xmin(), la.ymin()),
+ EPoint(la.xmin(), la.ymax()), m_color_grid);
+ }
+
+ const Partial_hds& la=m_shds.
+ partial_hds(m_shds.number_of_strips()-1);
+ add_segment(EPoint(la.xmax(), la.ymin()),
+ EPoint(la.xmax(), la.ymax()), m_color_grid);
+ }
+
+ void compute_face(const Global_halfedge& dh)
+ {
+ if (m_nofaces) return;
+ assert (m_shds.face(dh)->finite());
+
+ // We fill only closed faces.
+ Global_halfedge cur=dh;
+ Global_halfedge min=dh;
+ do
+ {
+ if (cur<min) min=cur;
+ m_shds.move_to_next(cur);
+ }
+ while(cur!=dh);
+
+ if (m_random_face_color)
+ {
+ CGAL::Random random(m_shds.hds(min).halfedges().index(min.halfedge()));
+ CGAL::Color c=get_random_color(random);
+ face_begin(c);
+ }
+ else
+ { face_begin(); }
+
+ cur=dh;
+ do
+ {
+ add_point_in_face(m_shds.point(cur));
+ m_shds.move_to_next(cur);
+ }
+ while(cur!=dh);
+
+ for (std::size_t i=0; i<m_shds.face(dh)->sons().size(); ++i)
+ {
+ assert(!m_shds.face(dh)->sons()[i]->finite());
+ add_point_in_face(m_shds.point(dh));
+ cur=Global_halfedge(m_shds.face(dh)->sons()[i]->strip_id(),
+ m_shds.face(dh)->sons()[i]->first_halfedge());
+ min=cur;
+ do
+ {
+ add_point_in_face(m_shds.point(cur));
+ if (m_num_face!=0)
+ { compute_edge(cur); }
+ m_shds.move_to_next(cur);
+ }
+ while(cur!=min);
+ add_point_in_face(m_shds.point(cur));
+ }
+ face_end();
+ }
+
+ void compute_edge(const Global_halfedge& phe)
+ {
+ const auto& p1 = m_shds.point(phe);
+ Global_halfedge d2=m_shds.opposite(phe);
+
+ if (m_cur_she.halfedge()!=nullptr && (phe==m_cur_she || m_cur_she==d2))
+ {
+ add_segment(p1, m_shds.point(d2), m_color_selected_he);
+ return;
+ }
+
+ bool concerned=(m_cur_strip==0?true:false);
+ if (m_cur_strip!=0)
+ {
+ CGAL::Comparison_result c1=CGAL::compare_x(p1, m_shds.partial_hds(m_cur_strip-1).m_min);
+ CGAL::Comparison_result c2=CGAL::compare_x(p1, m_shds.partial_hds(m_cur_strip-1).m_max);
+ if (c1!=CGAL::SMALLER && c2==CGAL::SMALLER) // Point in strip
+ { concerned=true; }
+ else
+ {
+ CGAL::Comparison_result c3=CGAL::compare_x(m_shds.point(d2), m_shds.partial_hds(m_cur_strip-1).m_min);
+ CGAL::Comparison_result c4=CGAL::compare_x(m_shds.point(d2), m_shds.partial_hds(m_cur_strip-1).m_max);
+ if (c3!=CGAL::SMALLER && c4==CGAL::SMALLER) // Point in strip
+ { concerned=true; }
+ else
+ { concerned=(c1==CGAL::SMALLER && c4!=CGAL::SMALLER); }
+ }
+ }
+ if (!concerned) return;
+
+ add_segment(p1, m_shds.point(d2));
+ }
+
+ void compute_vertex(const Global_halfedge& phe)
+ {
+ if (m_cur_she.halfedge()!=nullptr && m_shds.vertex(phe)==m_shds.vertex(m_cur_she))
+ {
+ add_point(m_shds.point(phe), m_color_selected_he);
+
+ Global_halfedge above(phe.strip_id(), m_shds.vertex(phe)->above_halfedge());
+ if (above.halfedge()!=nullptr)
+ {
+ m_shds.move_to_non_external(above);
+ add_segment(m_shds.point(phe), m_shds.point(above), m_color_above_links);
+ add_segment(m_shds.point(phe), m_shds.point(m_shds.opposite(above)),
+ m_color_above_links);
+ }
+ return;
+ }
+
+ bool concerned=(m_cur_strip==0?true:false);
+ if (m_cur_strip!=0)
+ {
+ CGAL::Comparison_result c1=CGAL::compare_x
+ (m_shds.point(phe), m_shds.partial_hds(m_cur_strip-1).m_min);
+ CGAL::Comparison_result c2=CGAL::compare_x
+ (m_shds.point(phe), m_shds.partial_hds(m_cur_strip-1).m_max);
+ if (c1!=CGAL::SMALLER && c2==CGAL::SMALLER) // Point in strip
+ { concerned=true; }
+ }
+ if (!concerned) return;
+
+ add_point(m_shds.point(phe));
+ }
+
+ void compute_elements()
+ {
+ clear();
+
+ if (m_draw_grid)
+ { draw_grid(); }
+
+ typename SHds::size_type markedges = m_shds.get_new_mark();
+ typename SHds::size_type markvertices = m_shds.get_new_mark();
+
+ std::size_t n=0;
+ auto it=m_shds.faces().begin(), itend=m_shds.faces().end();
+ for (++it; it!=itend; ++it ) // We need to jump over the first face
+ {
+ if (m_num_face==0 || it->finite())
+ {
+ ++n;
+ if (m_cur_she.halfedge()==nullptr && (m_num_face==0 || m_num_face==n))
+ { m_cur_she.set(it->strip_id(), it->first_halfedge()); }
+ if (m_num_face==0 || m_num_face==n)
+ {
+ assert(it->first_halfedge()!=nullptr);
+ Global_halfedge itf(it->strip_id(), it->first_halfedge());
+ Global_halfedge itfend=itf;
+
+ assert(!m_shds.partial_hds(itf).is_external(itf.halfedge()));
+
+ if (it->finite())
+ { compute_face(itf); }
+
+ do
+ {
+ if (!m_shds.is_marked(itf, markvertices))
+ {
+ compute_vertex(itf);
+ Global_halfedge itv=itf;
+ do
+ {
+ m_shds.mark(itv, markvertices);
+ itv.set_halfedge(m_shds.hds(itv).next_around_vertex(itv.halfedge()));
+ }
+ while(itv!=itf);
+ }
+
+ if (!m_shds.is_marked(itf, markedges))
+ {
+ compute_edge(itf);
+ m_shds.mark(itf, markedges);
+ m_shds.mark(m_shds.opposite(itf), markedges);
+ }
+
+ m_shds.move_to_next(itf);
+ }
+ while(itf!=itfend);
+ }
+ }
+ }
+
+ for (std::size_t i=0; i<m_shds.number_of_strips(); ++i)
+ {
+ for (auto it=m_shds.hds(i).halfedges().begin(),
+ itend=m_shds.hds(i).halfedges().end(); it!=itend; ++it)
+ {
+ m_shds.unmark(it, markvertices);
+ m_shds.unmark(it, markedges);
+ }
+ }
+
+ m_shds.free_mark(markedges);
+ m_shds.free_mark(markvertices);
+ }
+
+ virtual void init()
+ {
+ Base::init();
+
+ setKeyDescription(::Qt::Key_B, "Toggles display grid borders");
+
+ setKeyDescription(::Qt::Key_D, "Move current halfedge by next");
+ setKeyDescription(::Qt::ControlModifier+::Qt::Key_D, "Move current halfedge by previous");
+ setKeyDescription(::Qt::ShiftModifier+::Qt::Key_D, "Move current halfedge by opposite");
+
+ setKeyDescription(::Qt::Key_Z+::Qt::ControlModifier, "Draw all faces");
+ setKeyDescription(::Qt::ControlModifier+::Qt::Key_Left, "Decrement current face drawn");
+ setKeyDescription(::Qt::ControlModifier+::Qt::Key_Right, "Increment current face drawn");
+
+ setKeyDescription(::Qt::Key_T, "Increment current strip drawn");
+ setKeyDescription(::Qt::ShiftModifier+::Qt::Key_T, "Decrement current strip drawn");
+ setKeyDescription(::Qt::Key_T+::Qt::ControlModifier, "Draw all strips");
+
+ setMouseBindingDescription(::Qt::ControlModifier, ::Qt::LeftButton, "Display position of the point");
+ }
+
+ virtual void keyPressEvent(QKeyEvent *e)
+ {
+ const ::Qt::KeyboardModifiers modifiers = e->modifiers();
+ if ((e->key()==::Qt::Key_B) && (modifiers==::Qt::NoButton))
+ {
+ m_draw_grid=!m_draw_grid;
+ displayMessage(QString("Draw grid=%1.").arg(m_draw_grid?"true":"false"));
+ compute_elements();
+ redraw();
+ }
+ else if ((e->key()==::Qt::Key_D) && (modifiers==::Qt::NoButton))
+ {
+ m_shds.move_to_next(m_cur_she);
+ displayMessage(QString("Cur halfedge move by next."));
+ std::cout<<"Dart: "<<m_shds.point(m_cur_she)<<", "
+ <<m_shds.point(m_shds.opposite(m_cur_she))<<std::endl;
+ compute_elements();
+ redraw();
+ }
+ else if ((e->key()==::Qt::Key_D) && (modifiers==::Qt::ControlModifier))
+ {
+ m_shds.move_to_previous(m_cur_she);
+ displayMessage(QString("Cur halfedge move by previous."));
+ std::cout<<"Dart: "<<m_shds.point(m_cur_she)<<", "
+ <<m_shds.point(m_shds.opposite(m_cur_she))<<std::endl;
+ compute_elements();
+ redraw();
+ }
+ else if ((e->key()==::Qt::Key_D) && (modifiers==::Qt::ShiftModifier))
+ {
+ m_shds.move_to_opposite(m_cur_she);
+ std::cout<<"Dart: "<<m_shds.point(m_cur_she)<<", "
+ <<m_shds.point(m_shds.opposite(m_cur_she))<<std::endl;
+ displayMessage(QString("Cur halfedge move by opposite."));
+ compute_elements();
+ redraw();
+ }
+ else if ((e->key()==::Qt::Key_Z) && (modifiers==::Qt::ControlModifier))
+ {
+ m_num_face=0;
+ if (m_cur_she.halfedge()==nullptr && !m_shds.faces().empty())
+ {
+ m_cur_she.set(0, m_shds.hds(0).halfedges().begin());
+ m_cur_she.set(m_shds.faces().begin()->strip_id(),
+ m_shds.faces().begin()->first_halfedge());
+ }
+ displayMessage(QString("Draw all faces."));
+ compute_elements();
+ redraw();
+ }
+ else if ((e->key()==::Qt::Key_Right) && (modifiers==::Qt::ControlModifier))
+ {
+ ++m_num_face; m_cur_she.set_halfedge(nullptr);
+ if (m_num_face==m_shds.number_of_finite_faces()+1)
+ {
+ m_num_face=0;
+ displayMessage(QString("Draw all faces."));
+ }
+ else
+ {
+ displayMessage(QString("Draw face num %1.").arg(m_num_face));
+ }
+ compute_elements();
+ redraw();
+ }
+ else if ((e->key()==::Qt::Key_Left) && (modifiers==::Qt::ControlModifier))
+ {
+ m_cur_she.set_halfedge(nullptr);
+ if (m_num_face==0)
+ { m_num_face=m_shds.number_of_finite_faces(); }
+ else { --m_num_face; }
+
+ if (m_num_face==0)
+ { displayMessage(QString("Draw all faces.")); }
+ else
+ { displayMessage(QString("Draw face num %1.").arg(m_num_face)); }
+ compute_elements();
+ redraw();
+ }
+ else if ((e->key()==::Qt::Key_T) && (modifiers==::Qt::NoButton))
+ {
+ ++m_cur_strip;
+ if (m_cur_strip==m_shds.number_of_strips()+1)
+ {
+ m_cur_strip=0;
+ displayMessage(QString("Draw all strips."));
+ }
+ else
+ { displayMessage(QString("Draw strip num %1.").arg(m_cur_strip-1)); }
+ compute_elements();
+ redraw();
+ }
+ else if ((e->key()==::Qt::Key_T) && (modifiers==::Qt::ShiftModifier))
+ {
+ if (m_cur_strip==0)
+ { m_cur_strip=m_shds.number_of_strips(); }
+ else { --m_cur_strip; }
+
+ if (m_cur_strip==0)
+ { displayMessage(QString("Draw all strips.")); }
+ else
+ { displayMessage(QString("Draw strip num %1.").arg(m_cur_strip-1)); }
+ compute_elements();
+ redraw();
+ }
+ else if ((e->key()==::Qt::Key_T) && (modifiers==::Qt::ControlModifier))
+ {
+ m_cur_strip=0;
+ displayMessage(QString("Draw all strips."));
+ compute_elements();
+ redraw();
+ }
+ else
+ { Base::keyPressEvent(e); } // Call the base method to process others/classicals key
+
+ // Call: * compute_elements() if the model changed, followed by
+ // * redraw() if some viewing parameters changed that implies some
+ // modifications of the buffers
+ // (eg. type of normal, color/mono)
+ // * update() just to update the drawing
+ }
+
+ void mousePressEvent (QMouseEvent * m)
+ {
+ if (m->modifiers()==Qt::ControlModifier)
+ {
+ CGAL::qglviewer::Vec globalp=camera()->unprojectedCoordinatesOf
+ (CGAL::qglviewer::Vec(m->x(), m->y(), 0));
+ displayMessage(QString("Point (%1, %2).").arg(globalp.x).arg(globalp.y));
+ }
+ Base::mousePressEvent(m);
+ }
+
+protected:
+ SHds& m_shds;
+ bool m_nofaces;
+ bool m_random_face_color;
+ std::size_t m_num_face;
+ Global_halfedge m_cur_she;
+ bool m_draw_grid;
+ std::size_t m_cur_strip;
+
+ CGAL::Color m_color_points;
+ CGAL::Color m_color_grid;
+ CGAL::Color m_color_selected_he;
+ CGAL::Color m_color_above_links;
+};
+
+void draw_shds(SHds& ashds,
+ const char* title="Basic Viewer for sHDS",
+ bool nofill=false)
+{
+#if defined(CGAL_TEST_SUITE)
+ bool cgal_test_suite=true;
+#else
+ bool cgal_test_suite=qEnvironmentVariableIsSet("CGAL_TEST_SUITE");
+#endif
+
+ if (!cgal_test_suite)
+ {
+ int argc=1;
+ const char* argv[2]={"hdsviewer","\0"};
+ QApplication app(argc,const_cast<char**>(argv));
+ MySimpleSHDSViewerQt mainwindow(app.activeWindow(), ashds, title, nofill);
+ mainwindow.show();
+ app.exec();
+ }
+}
+
+#else // CGAL_USE_BASIC_VIEWER
+
+void draw_shds(SHds& /*ashds*/,
+ const char* /*title*/="Basic Viewer for sHDS",
+ bool /*nofill*/=false)
+{
+}
+
+
+#endif // CGAL_USE_BASIC_VIEWER
+
+#endif // DRAW_SHDS_H
diff --git a/src/memory.h b/src/memory.h
new file mode 100644
index 0000000000000000000000000000000000000000..1c6715252d38c50e0b22f978b5b7a71022f531cb
--- /dev/null
+++ b/src/memory.h
@@ -0,0 +1,127 @@
+/*
+ * Author: David Robert Nadeau
+ * Site: http://NadeauSoftware.com/
+ * License: Creative Commons Attribution 3.0 Unported License
+ * http://creativecommons.org/licenses/by/3.0/deed.en_US
+ */
+
+#ifndef MEMORY_H
+#define MEMORY_H
+
+#if defined(_WIN32)
+#include <windows.h>
+#include <psapi.h>
+
+#elif defined(__unix__) || defined(__unix) || defined(unix) || (defined(__APPLE__) && defined(__MACH__))
+#include <unistd.h>
+#include <sys/resource.h>
+
+#if defined(__APPLE__) && defined(__MACH__)
+#include <mach/mach.h>
+
+#elif (defined(_AIX) || defined(__TOS__AIX__)) || (defined(__sun__) || defined(__sun) || defined(sun) && (defined(__SVR4) || defined(__svr4__)))
+#include <fcntl.h>
+#include <procfs.h>
+
+#elif defined(__linux__) || defined(__linux) || defined(linux) || defined(__gnu_linux__)
+#include <stdio.h>
+
+#endif
+
+#else
+#error "Cannot define getPeakRSS( ) or getCurrentRSS( ) for an unknown OS."
+#endif
+
+
+
+
+
+/**
+ * Returns the peak (maximum so far) resident set size (physical
+ * memory use) measured in bytes, or zero if the value cannot be
+ * determined on this OS.
+ */
+size_t getPeakRSS( )
+{
+#if defined(_WIN32)
+ /* Windows -------------------------------------------------- */
+ PROCESS_MEMORY_COUNTERS info;
+ GetProcessMemoryInfo( GetCurrentProcess( ), &info, sizeof(info) );
+ return (size_t)info.PeakWorkingSetSize;
+
+#elif (defined(_AIX) || defined(__TOS__AIX__)) || (defined(__sun__) || defined(__sun) || defined(sun) && (defined(__SVR4) || defined(__svr4__)))
+ /* AIX and Solaris ------------------------------------------ */
+ struct psinfo psinfo;
+ int fd = -1;
+ if ( (fd = open( "/proc/self/psinfo", O_RDONLY )) == -1 )
+ return (size_t)0L; /* Can't open? */
+ if ( read( fd, &psinfo, sizeof(psinfo) ) != sizeof(psinfo) )
+ {
+ close( fd );
+ return (size_t)0L; /* Can't read? */
+ }
+ close( fd );
+ return (size_t)(psinfo.pr_rssize * 1024L);
+
+#elif defined(__unix__) || defined(__unix) || defined(unix) || (defined(__APPLE__) && defined(__MACH__))
+ /* BSD, Linux, and OSX -------------------------------------- */
+ struct rusage rusage;
+ getrusage( RUSAGE_SELF, &rusage );
+#if defined(__APPLE__) && defined(__MACH__)
+ return (size_t)rusage.ru_maxrss;
+#else
+ return (size_t)(rusage.ru_maxrss * 1024L);
+#endif
+
+#else
+ /* Unknown OS ----------------------------------------------- */
+ return (size_t)0L; /* Unsupported. */
+#endif
+}
+
+
+
+
+
+/**
+ * Returns the current resident set size (physical memory use) measured
+ * in bytes, or zero if the value cannot be determined on this OS.
+ */
+size_t getCurrentRSS( )
+{
+#if defined(_WIN32)
+ /* Windows -------------------------------------------------- */
+ PROCESS_MEMORY_COUNTERS info;
+ GetProcessMemoryInfo( GetCurrentProcess( ), &info, sizeof(info) );
+ return (size_t)info.WorkingSetSize;
+
+#elif defined(__APPLE__) && defined(__MACH__)
+ /* OSX ------------------------------------------------------ */
+ struct mach_task_basic_info info;
+ mach_msg_type_number_t infoCount = MACH_TASK_BASIC_INFO_COUNT;
+ if ( task_info( mach_task_self( ), MACH_TASK_BASIC_INFO,
+ (task_info_t)&info, &infoCount ) != KERN_SUCCESS )
+ return (size_t)0L; /* Can't access? */
+ return (size_t)info.resident_size;
+
+#elif defined(__linux__) || defined(__linux) || defined(linux) || defined(__gnu_linux__)
+ /* Linux ---------------------------------------------------- */
+ long rss = 0L;
+ FILE* fp = NULL;
+ if ( (fp = fopen( "/proc/self/statm", "r" )) == NULL )
+ return (size_t)0L; /* Can't open? */
+ if ( fscanf( fp, "%*s%ld", &rss ) != 1 )
+ {
+ fclose( fp );
+ return (size_t)0L; /* Can't read? */
+ }
+ fclose( fp );
+ return (size_t)rss * (size_t)sysconf( _SC_PAGESIZE);
+
+#else
+ /* AIX, BSD, Solaris, and Unknown OS ------------------------ */
+ return (size_t)0L; /* Unsupported. */
+#endif
+}
+
+#endif // MEMORY_H
diff --git a/src/parallel-arrangement-bench.cpp b/src/parallel-arrangement-bench.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..164746b6da5a25916bd28c870a828b0d9236edaf
--- /dev/null
+++ b/src/parallel-arrangement-bench.cpp
@@ -0,0 +1,282 @@
+// compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <stack>
+#include <tuple>
+#include <thread>
+#include <boost/filesystem.hpp>
+
+#include <CGAL/Memory_sizer.h>
+#include <CGAL/assertions.h>
+#include <CGAL/Surface_sweep_2.h>
+#include <CGAL/Surface_sweep_2_algorithms.h>
+#include <CGAL/Arr_segment_traits_2.h>
+#include <CGAL/Arrangement_2.h>
+#include <CGAL/Arr_consolidated_curve_data_traits_2.h>
+#include <CGAL/Arr_non_caching_segment_basic_traits_2.h>
+
+#include "CGAL_typedefs.h"
+#include "Segment_readers.h"
+#include "My_visitor.h"
+#include "SHds.h"
+#include "SHds_to_lcc.h"
+
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////
+[[ noreturn ]] void usage(int /*argc*/, char** argv)
+{
+ // Name
+ std::cout<<"Name"<<std::endl;
+ std::cout<<" "<<argv[0]<<" - a new method to compute arrangement of segments";
+ std::cout<<std::endl<<std::endl;
+ // Synopsis
+ std::cout<<"SYNOPSIS"<<std::endl;
+ std::cout<<" "<<argv[0]<<" [--help|-h|-?] "
+ <<"[-t1 filename] [-nbs S] [-nbt N] [-crop] [-optimized-strips] [-cgal]"
+ <<std::endl<<std::endl;
+ // Description
+ std::cout<<"DESCRIPTION"<<std::endl;
+ std::cout<<" "<<"Compute the arrangement of segments given in the filename."
+ <<std::endl<<std::endl;
+ // Options
+ std::cout<<" --help, -h, -?"<<std::endl
+ <<" display this help and exit."
+ <<std::endl<<std::endl;
+ std::cout<<" -t1 filename"
+ <<std::endl
+ <<" load the text segment file, and creates the arrangement: required."
+ <<std::endl<<std::endl;
+ std::cout<<" -nbs S"
+ <<std::endl
+ <<" number of strips (1 by default)."
+ <<std::endl<<std::endl;
+ std::cout<<" -nbt N"
+ <<std::endl
+ <<" number of threads (1 by default)."
+ <<std::endl<<std::endl;
+ std::cout<<" -crop"
+ <<std::endl
+ <<" crop the segments to fit in its strip."
+ <<std::endl<<std::endl;
+ std::cout<<" -repeat R"
+ <<std::endl
+ <<" repeat each test R times (1 by default)."
+ <<std::endl<<std::endl;
+ std::cout<<" -optimized-strips"
+ <<std::endl
+ <<" compute the length of each strips dynamically (by default each strip has the same size) (used only if the number of threads is > 1)."
+ <<std::endl<<std::endl;
+ std::cout<<" -cgal"
+ <<std::endl
+ <<" run cgal arrangement (and not our new method)."
+ <<std::endl<<std::endl;
+// Author
+ std::cout<<"AUTHOR"<<std::endl;
+ std::cout<<" "<<"Written by Guillaume Damiand."
+ <<std::endl<<std::endl;
+ // Copyright
+ std::cout<<"COPYRIGHT"<<std::endl;
+ std::cout<<" "<<"Copyright (C) 2020 CNRS and LIRIS' Establishments (France). "
+ <<"License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>"
+ <<std::endl
+ <<" This is free software: you are free to change and redistribute it under certain conditions. "
+ <<"There is NO WARRANTY, to the extent permitted by law."
+ <<std::endl<<std::endl;
+ // Reporting bugs
+ std::cout<<"REPORTING BUGS"<<std::endl;
+ std::cout<<" "<<"Email bug reports to Guillaume Damiand ⟨guillaume.damiand@liris.cnrs.fr⟩."
+ <<std::endl<<std::endl;
+ exit(EXIT_FAILURE);
+}
+////////////////////////////////////////////////////////////////////////////////
+[[ noreturn ]] void error_command_line(int argc, char** argv, const char* msg)
+{
+ std::cout<<"ERROR: "<<msg<<std::endl;
+ usage(argc, argv);
+}
+////////////////////////////////////////////////////////////////////////////////
+void process_command_line(int argc, char** argv,
+ bool& t1,
+ std::string& filename,
+ std::size_t& nbstrips,
+ std::size_t& nbthreads,
+ std::size_t& repeat,
+ bool &optimized_strips,
+ bool &cgal,
+ bool &crop)
+{
+ t1=false;
+ filename="";
+ nbstrips=1;
+ nbthreads=1;
+ optimized_strips=false;
+ cgal=false;
+ repeat=1;
+ crop=false;
+
+ bool helprequired=false;
+ std::string arg;
+
+ for (int i=1; i<argc; ++i)
+ {
+ arg=std::string(argv[i]);
+ if (arg==std::string("-h") || arg==std::string("--help") || arg==std::string("-?"))
+ { helprequired=true; }
+ else if (arg==std::string("-t1"))
+ {
+ t1=true;
+ if (argc-1-i<1)
+ { error_command_line(argc, argv, "Error: no filename after -t1 option."); }
+ filename=std::string(argv[++i]);
+ }
+ else if (arg==std::string("-nbs"))
+ {
+ if (argc-1-i<1)
+ { error_command_line(argc, argv, "Error: no number after -nbs option."); }
+ nbstrips=std::stoi(std::string(argv[++i]));
+ }
+ else if (arg==std::string("-nbt"))
+ {
+ if (argc-1-i<1)
+ { error_command_line(argc, argv, "Error: no number after -nbt option."); }
+ nbthreads=std::stoi(std::string(argv[++i]));
+ }
+ else if (arg==std::string("-repeat"))
+ {
+ if (argc-1-i<1)
+ { error_command_line(argc, argv, "Error: no number after -repeat option."); }
+ repeat=std::stoi(std::string(argv[++i]));
+ }
+ else if (arg==std::string("-optimized-strips"))
+ { optimized_strips=true; }
+ else if (arg==std::string("-cgal"))
+ { cgal=true; }
+ else if (arg==std::string("-crop"))
+ { crop=true; }
+ else if (arg[0]=='-')
+ {
+ std::cout<<"Unknown option "<<arg<<std::endl;
+ helprequired=true;
+ }
+ }
+
+ if (nbthreads==0) { nbthreads=1; }
+ if (repeat==0) { repeat=1; }
+ if (!t1) { std::cout<<"ERROR: missing txt file name."<<std::endl; }
+ if (helprequired || !t1) { usage(argc, argv); }
+}
+///////////////////////////////////////////////////////////////////////////////
+int main(int argc, char** argv)
+{
+ std::string filename;
+ bool t1;
+ std::size_t nbstrips, nbthreads, repeat;
+ bool optimized_strips;
+ bool cgal;
+ bool crop;
+
+ process_command_line(argc, argv, t1, filename, nbstrips, nbthreads, repeat,
+ optimized_strips, cgal, crop);
+
+ std::ofstream ftimes("times.dat", std::ios_base::app);
+ std::ofstream fdetailedtimes("detailed-times.dat", std::ios_base::app);
+ std::ofstream fcells("nbcells.dat", std::ios_base::app);
+ ftimes.setf( std::ios::fixed, std:: ios::floatfield );
+ ftimes.precision(2);
+ fdetailedtimes.setf( std::ios::fixed, std:: ios::floatfield );
+ fdetailedtimes.precision(2);
+
+ boost::filesystem::path p(filename);
+
+ if (cgal)
+ {
+ ftimes<<p.stem()<<'\t';
+ fdetailedtimes<<p.stem()<<'\t';
+ fcells<<"####################"<<p.stem()<<std::endl;
+ }
+
+ if (cgal)
+ {
+ typedef CGAL::Arr_segment_traits_2<EPECK> Traits_2;
+ typedef CGAL::Arrangement_2<Traits_2> Arrangement_2;
+ std::vector<ESegment> all_segments;
+ all_segments.reserve(number_of_lines(filename));
+ Read_from_file<EPECK> reader(filename);
+ fill_segment_array<EPECK>(reader, all_segments);
+
+ My_timer timer;
+ for (std::size_t i=0; i<repeat; ++i)
+ {
+ timer.start();
+ Arrangement_2 arr;
+ CGAL::insert(arr, all_segments.begin(), all_segments.end());
+ timer.stop();
+
+ if (i==0)
+ {
+ fcells<<arr.number_of_halfedges()<<'\t'
+ <<arr.number_of_vertices()<<'\t'
+ <<arr.number_of_edges()<<'\t'
+ <<arr.number_of_faces()-arr.number_of_unbounded_faces()<<'\t'
+ <<arr.number_of_faces()<<std::endl;
+ }
+ }
+
+ ftimes<<timer.time()/repeat<<'\t';
+ }
+ else
+ {
+ std::vector<double> times;
+ // times[0] : global time
+ // times[1] : load and dispatch
+ // times[2] : compute all local arrangements
+ // times[3] : compute faces
+
+ for (std::size_t i=0; i<repeat; ++i)
+ {
+ SHds shds(filename, nbstrips, nbthreads, optimized_strips, times, crop);
+ if (i==0)
+ {
+ LCC_2 lcc;
+ std::size_t finite_faces=shds_to_lcc(shds, lcc); // TODO NO MORE NEED !!
+ std::vector<unsigned int> nbcells=lcc.count_all_cells();
+
+ fcells<<shds.number_of_halfedges()<<'\t'
+ <<nbcells[0]<<'\t'<<nbcells[1]<<'\t'
+ <<finite_faces<<'\t'<<nbcells[2]
+ <<'\t'<<shds.number_of_external_halfedges()<<std::endl;
+ }
+ }
+
+ ftimes<<(times[0]-times[1])/repeat<<'\t';
+ fdetailedtimes<<times[1]/repeat<<'\t'
+ <<times[2]/repeat<<'\t'
+ <<times[3]/repeat<<'\t';
+ }
+
+ return EXIT_SUCCESS;
+}
+///////////////////////////////////////////////////////////////////////////////
diff --git a/src/parallel-arrangement.cpp b/src/parallel-arrangement.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ae895ffb62f24990a9709b5cb4bf868039682a7f
--- /dev/null
+++ b/src/parallel-arrangement.cpp
@@ -0,0 +1,590 @@
+// compute_arrangement: a new method to compute arrangement of segments.
+// Copyright (C) 2020 CNRS and LIRIS' Establishments (France).
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <https://www.gnu.org/licenses/>.
+//
+// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+//
+
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <stack>
+#include <tuple>
+#include <thread>
+#include <boost/filesystem.hpp>
+
+#include <CGAL/assertions.h>
+#include <CGAL/Surface_sweep_2.h>
+#include <CGAL/Surface_sweep_2_algorithms.h>
+#include <CGAL/Arr_segment_traits_2.h>
+#include <CGAL/Arrangement_2.h>
+#include <CGAL/Arr_consolidated_curve_data_traits_2.h>
+#include <CGAL/Arr_non_caching_segment_basic_traits_2.h>
+
+#include "CGAL_typedefs.h"
+#include "Segment_readers.h"
+#include "My_visitor.h"
+#include "SHds.h"
+#include "SHds_to_lcc.h"
+#include "SHds_to_cgal_arrangement.h"
+#include "memory.h"
+#include "draw_shds.h"
+
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////
+[[ noreturn ]] void usage(int /*argc*/, char** argv)
+{
+ // Name
+ std::cout<<"Name"<<std::endl;
+ std::cout<<" "<<argv[0]<<" - a new method to compute arrangement of segments";
+ std::cout<<std::endl<<std::endl;
+ // Synopsis
+ std::cout<<"SYNOPSIS"<<std::endl;
+ std::cout<<" "<<argv[0]<<" [--help|-h|-?] "
+ <<"[-t1 FILENAME] [-t2 FILENAME] [-t3 FILENAME] [-cgal] [-lcc] [-traverse]"
+ <<"[-nbs S] [-nbt N] [-crop] [-optimized-strips] [-xpos X1 X2 ...] "
+ <<"[-export-cgal] [-draw] [-save] [-export-xfig XFIGFILE] "
+ <<"[-export-xfig-sequence] [-xfig-invert-y] [-xfig-x-size S] "
+ <<"[-xfig-width W1 W2 W3 W4 W5] "
+ <<"[-xfig-color R1 G1 B1 R2 G2 B2 R3 G3 B3 R4 G4 B4]"
+ <<std::endl<<std::endl;
+ // Description
+ std::cout<<"DESCRIPTION"<<std::endl;
+ std::cout<<" "<<"Compute the arrangement of segments given in the FILENAME."
+ <<std::endl<<std::endl;
+ // Options
+ std::cout<<" --help, -h, -?"<<std::endl
+ <<" display this help and exit."
+ <<std::endl<<std::endl;
+ std::cout<<" -t1 FILENAME"
+ <<std::endl
+ <<" load the text segment file, and creates the arrangement."
+ <<std::endl<<std::endl;
+ std::cout<<" -t2 FILENAME"
+ <<std::endl
+ <<" load the sHDS from FILENAME."
+ <<std::endl<<std::endl;
+ std::cout<<" -t3 FILENAME"
+ <<std::endl
+ <<" load the sequence of HDS from FILENAME, which must contains % at the position of the number."
+ <<std::endl<<std::endl;
+ std::cout<<" -cgal"
+ <<std::endl
+ <<" use CGAL arrangement (and not our method as default)."
+ <<std::endl<<std::endl;
+ std::cout<<" -lcc"
+ <<std::endl
+ <<" compute a 2D Linear cell complex from the sHDS."
+ <<std::endl<<std::endl;
+ std::cout<<" -traverse"
+ <<std::endl
+ <<" traverse the sHDS and compute time."
+ <<std::endl<<std::endl;
+ std::cout<<" -nbs S"
+ <<std::endl
+ <<" number of strips (1 by default)."
+ <<std::endl<<std::endl;
+ std::cout<<" -nbt N"
+ <<std::endl
+ <<" number of threads (1 by default)."
+ <<std::endl<<std::endl;
+ std::cout<<" -crop"
+ <<std::endl
+ <<" crop the segments to fit in its strip."
+ <<std::endl<<std::endl;
+ std::cout<<" -optimized-strips"
+ <<std::endl
+ <<" compute the length of each strips dynamically (by default each strip has the same size) (used only if the number of strips is > 1) (not used if option -xpos is used)."
+ <<std::endl<<std::endl;
+ std::cout<<" -xpos X1 X2 ..."
+ <<std::endl
+ <<" gives the x positions of the strips (we must have S-1 numbers, B being the number given after -nbs option, and this option must be used after -nbs option)."
+ <<std::endl<<std::endl;
+ std::cout<<" -export-cgal"
+ <<std::endl
+ <<" compute a CGAL arrangement from the sHDS."
+ <<std::endl<<std::endl;
+ std::cout<<" -draw"
+ <<std::endl
+ <<" draw the sHDS."
+ <<std::endl<<std::endl;
+ std::cout<<" -save"
+ <<std::endl
+ <<" save the sHDS."
+ <<std::endl<<std::endl;
+ std::cout<<" -export-xfig XFIGFILE"
+ <<std::endl
+ <<" export the sHDS into the given xfig file."
+ <<std::endl<<std::endl;
+ std::cout<<" -export-xfig-sequence"
+ <<std::endl
+ <<" export the strips as a sequence of xfig files."
+ <<std::endl<<std::endl;
+ std::cout<<" -xfig-invert-y"
+ <<std::endl
+ <<" invert y axis for the export xfig."
+ <<std::endl<<std::endl;
+ std::cout<<" -xfig-x-size S"
+ <<std::endl
+ <<" set the x size of the xfig file to S (default 50000.)."
+ <<std::endl<<std::endl;
+ std::cout<<" -xfig-width W1 W2 W3 W4 W5"
+ <<std::endl
+ <<" fix the size of: segments=W1; disks=W2; criticals=W3; borders=W4; "
+ <<"space between strips=W5 (default 2 30 2 1 100) (if one width is 0, corresponding elements are not drawn)."
+ <<std::endl<<std::endl;
+ std::cout<<" -xfig-color R1 G1 B1 R2 G2 B2 R3 G3 B3 R4 G4 B4"
+ <<std::endl
+ <<" fix the color of: segments=R1 G1 B1; disks=R2 G2 B2; "
+ <<"criticals=R3 G3 B3; borders=R4 G4 B4 (default 40 40 220 220 40 40 90 190 110 150 150 150)."
+ <<std::endl<<std::endl;
+ // Author
+ std::cout<<"AUTHOR"<<std::endl;
+ std::cout<<" "<<"Written by Guillaume Damiand."
+ <<std::endl<<std::endl;
+ // Copyright
+ std::cout<<"COPYRIGHT"<<std::endl;
+ std::cout<<" "<<"Copyright (C) 2020 CNRS and LIRIS' Establishments (France). "
+ <<"License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>"
+ <<std::endl
+ <<" This is free software: you are free to change and redistribute it under certain conditions. "
+ <<"There is NO WARRANTY, to the extent permitted by law."
+ <<std::endl<<std::endl;
+ // Reporting bugs
+ std::cout<<"REPORTING BUGS"<<std::endl;
+ std::cout<<" "<<"Email bug reports to Guillaume Damiand ⟨guillaume.damiand@liris.cnrs.fr⟩."
+ <<std::endl<<std::endl;
+ exit(EXIT_FAILURE);
+}
+////////////////////////////////////////////////////////////////////////////////
+[[ noreturn ]] void error_command_line(int argc, char** argv, const char* msg)
+{
+ std::cout<<"ERROR: "<<msg<<std::endl;
+ usage(argc, argv);
+}
+////////////////////////////////////////////////////////////////////////////////
+void process_command_line(int argc, char** argv,
+ std::string& filename,
+ bool& t1,
+ bool &t2,
+ bool &t3,
+ bool &cgal,
+ bool &lcc,
+ bool &traverse,
+ std::size_t& nbstrips,
+ std::size_t& nbthreads,
+ bool &crop,
+ bool &optimized_strips,
+ std::vector<double>& xpos_cuts,
+ bool &export_cgal,
+ bool &draw,
+ bool &save,
+ std::string& xfig_file,
+ bool &export_xfig_sequence,
+ bool &xfig_invert_y,
+ double &xfig_x_size,
+ uint16_t& width_segments,
+ uint16_t& radius_disks,
+ uint16_t& width_criticals,
+ uint16_t& width_borders,
+ uint16_t& space_between_strips,
+ CGAL::Color& color_segments,
+ CGAL::Color& color_disks,
+ CGAL::Color& color_criticals,
+ CGAL::Color& color_borders)
+ {
+ t1=false;
+ t2=false;
+ t3=false;
+ cgal=false;
+ lcc=false;
+ traverse=false;
+ filename="";
+ nbstrips=1;
+ nbthreads=1;
+ draw=false;
+ optimized_strips=false;
+ export_cgal=false;
+ save=false;
+ crop=false;
+ xfig_file="";
+ export_xfig_sequence=false;
+ xpos_cuts.clear();
+
+ width_segments=2;
+ radius_disks=30;
+ width_criticals=2;
+ width_borders=1;
+ space_between_strips=100;
+
+ xfig_invert_y=false;
+ xfig_x_size=50000;
+
+ color_segments=CGAL::Color(40,40,220);
+ color_disks=CGAL::Color(220,40,40);
+ color_criticals=CGAL::Color(90,190,110),
+ color_borders=CGAL::Color(150,150,150);
+
+ bool helprequired=false;
+ std::string arg;
+
+ for (int i=1; i<argc; ++i)
+ {
+ arg=std::string(argv[i]);
+ if (arg==std::string("-h") || arg==std::string("--help") || arg==std::string("-?"))
+ { helprequired=true; }
+ else if (arg==std::string("-t1"))
+ {
+ t1=true;
+ if (argc-1-i<1)
+ { error_command_line(argc, argv, "Error: no filename after -t1 option."); }
+ filename=std::string(argv[++i]);
+ }
+ else if (arg==std::string("-t2"))
+ {
+ t2=true;
+ if (argc-1-i<1)
+ { error_command_line(argc, argv, "Error: no filename after -t2 option."); }
+ filename=std::string(argv[++i]);
+ }
+ else if (arg==std::string("-t3"))
+ {
+ t3=true;
+ if (argc-1-i<1)
+ { error_command_line(argc, argv, "Error: no filename after -t3 option."); }
+ filename=std::string(argv[++i]);
+ if (filename.find("%")==std::string::npos)
+ { error_command_line(argc, argv, "Error: filename after -t3 option does not contain % character."); }
+ }
+ else if (arg==std::string("-cgal"))
+ { cgal=true; }
+ else if (arg==std::string("-lcc"))
+ { lcc=true; }
+ else if (arg==std::string("-traverse"))
+ { traverse=true; }
+ else if (arg==std::string("-nbs"))
+ {
+ if (argc-1-i<1)
+ { error_command_line(argc, argv, "Error: no number after -nbs option."); }
+ nbstrips=std::stoi(std::string(argv[++i]));
+ }
+ else if (arg==std::string("-nbt"))
+ {
+ if (argc-1-i<1)
+ { error_command_line(argc, argv, "Error: no number after -nbt option."); }
+ nbthreads=std::stoi(std::string(argv[++i]));
+ }
+ else if (arg==std::string("-crop"))
+ { crop=true; }
+ else if (arg==std::string("-optimized-strips"))
+ { optimized_strips=true; }
+ else if (arg==std::string("-xpos"))
+ {
+ if (nbstrips<=1)
+ { error_command_line(argc, argv, "Error: -xpos option can only be used after-nbs option, and with a number of strips >1."); }
+ if (argc-1-i<nbstrips-1)
+ { error_command_line(argc, argv, "Error: we need (numberofstrips-1) numbers -xpos option."); }
+
+ xpos_cuts.resize(nbstrips-1);
+ for (std::size_t j=0; j<nbstrips-1; ++j)
+ { xpos_cuts[j]=std::stod(std::string(argv[++i])); }
+ }
+ else if (arg==std::string("-export-cgal"))
+ { export_cgal=true; }
+ else if (arg==std::string("-draw"))
+ { draw=true; }
+ else if (arg==std::string("-save"))
+ { save=true; }
+ else if (arg==std::string("-export-xfig"))
+ {
+ if (argc-1-i<1)
+ { error_command_line(argc, argv, "Error: no number after -export-xfig option."); }
+ xfig_file=std::string(argv[i+1]);
+ }
+ else if (arg==std::string("-export-xfig-sequence"))
+ { export_xfig_sequence=true; }
+ else if (arg==std::string("-xfig-invert-y"))
+ { xfig_invert_y=true; }
+ else if (arg==std::string("-xfig-x-size"))
+ {
+ if (argc-1-i<1)
+ { error_command_line(argc, argv, "Error: no number after -xfig-x-size option."); }
+ xfig_x_size=std::stod(std::string(argv[++i]));
+ }
+ else if (arg==std::string("-xfig-width"))
+ {
+ if (argc-1-i<5)
+ { error_command_line(argc, argv, "Error: we need 5 numbers after -xfig-width option."); }
+ width_segments=std::stoi(std::string(argv[++i]));
+ radius_disks=std::stod(std::string(argv[++i]));
+ width_criticals=std::stoi(std::string(argv[++i]));
+ width_borders=std::stoi(std::string(argv[++i]));
+ space_between_strips=std::stoi(std::string(argv[++i]));
+ }
+ else if (arg==std::string("-xfig-color"))
+ {
+ if (argc-1-i<12)
+ { error_command_line(argc, argv, "Error: we need 12 numbers after -xfig-color option."); }
+ color_segments=CGAL::Color(std::stoi(std::string(argv[++i])),
+ std::stoi(std::string(argv[++i])),
+ std::stoi(std::string(argv[++i])));
+ color_disks=CGAL::Color(std::stoi(std::string(argv[++i])),
+ std::stoi(std::string(argv[++i])),
+ std::stoi(std::string(argv[++i])));
+ color_criticals=CGAL::Color(std::stoi(std::string(argv[++i])),
+ std::stoi(std::string(argv[++i])),
+ std::stoi(std::string(argv[++i])));
+ color_borders=CGAL::Color(std::stoi(std::string(argv[++i])),
+ std::stoi(std::string(argv[++i])),
+ std::stoi(std::string(argv[++i])));
+ }
+ else if (arg[0]=='-')
+ {
+ std::cout<<"Unknown option "<<arg<<std::endl;
+ helprequired=true;
+ }
+ }
+
+ if (nbthreads==0) { nbthreads=1; }
+ if (!t1 && !t2 && !t3) { std::cout<<"ERROR: missing file name (either with -t1 or -t2 or -t3)."<<std::endl; }
+ if (cgal && (t2 || t3)) { std::cout<<"ERROR: -t2 and -t3 option cannot be used with -cgal."<<std::endl; }
+ if (helprequired || (!t1 && !t2 && !t3)) { usage(argc, argv); }
+}
+///////////////////////////////////////////////////////////////////////////////
+void run_cgal_method(const std::string& filename)
+{
+ typedef CGAL::Arr_segment_traits_2<EPECK> Traits_2;
+ typedef CGAL::Arrangement_2<Traits_2> Arrangement_2;
+
+ My_timer timer;
+ timer.start();
+ Arrangement_2 arr;
+
+ My_timer load_timer;
+ load_timer.start();
+
+ std::vector<ESegment> all_segments;
+ all_segments.reserve(number_of_lines(filename));
+ Read_from_file<EPECK> reader(filename);
+ fill_segment_array<EPECK>(reader, all_segments);
+
+ load_timer.stop();
+
+ CGAL::insert(arr, all_segments.begin(), all_segments.end());
+ timer.stop();
+
+ print_txt_with_endl("Nb input segments: ", all_segments.size());
+ print_txt_with_endl("#halfedges=", arr.number_of_halfedges(), ", #0-cells=",
+ arr.number_of_vertices(), ", #1-cells=",
+ arr.number_of_edges(), ", #2-cells=",
+ arr.number_of_faces(), ", #finite-2-cells=",
+ arr.number_of_faces()-arr.number_of_unbounded_faces());
+
+ std::cout<<"CURRENT RSS: "<<getCurrentRSS()<<" bytes (i.e. "<<
+ double(getCurrentRSS())/double(1024*1024)<<" mega-bytes and "<<
+ double(getCurrentRSS())/double(1024*1024*1024)<<" giga-bytes)"<<std::endl;
+
+ print_txt_with_endl("[TIME]: to load segments ", load_timer.time());
+ print_txt_with_endl("[GLOBAL]: time to build cgal arrangement ", timer.time());
+}
+///////////////////////////////////////////////////////////////////////////////
+int main(int argc, char** argv)
+{
+ std::string filename;
+ bool t1;
+ bool t2;
+ bool t3;
+ bool cgal;
+ bool lcc;
+ bool traverse;
+ std::size_t nbstrips;
+ std::size_t nbthreads;
+ bool crop;
+ bool optimized_strips;
+ std::vector<double> xpos_cuts;
+ bool export_cgal;
+ bool draw;
+ bool save;
+ std::string xfig_file;
+ bool export_xfig_sequence;
+ bool xfig_invert_y;
+ double xfig_x_size;
+ uint16_t width_segments;
+ uint16_t radius_disks;
+ uint16_t width_criticals;
+ uint16_t width_borders;
+ uint16_t space_between_strips;
+ CGAL::Color color_segments; //(svg::Color::Transparent);
+ CGAL::Color color_disks;
+ CGAL::Color color_criticals;
+ CGAL::Color color_borders;
+
+ std::vector<double> times;
+ process_command_line(argc, argv, filename, t1, t2, t3, cgal, lcc, traverse,
+ nbstrips, nbthreads, crop, optimized_strips,
+ xpos_cuts, export_cgal, draw,
+ save, xfig_file, export_xfig_sequence, xfig_invert_y,
+ xfig_x_size, width_segments, radius_disks,
+ width_criticals, width_borders, space_between_strips,
+ color_segments, color_disks, color_criticals,
+ color_borders);
+
+ if (cgal)
+ {
+ run_cgal_method(filename);
+ exit(EXIT_SUCCESS);
+ }
+
+ SHds *shds=nullptr;
+ if (t1)
+ {
+ shds=new SHds(filename, nbstrips, nbthreads,
+ optimized_strips, times, crop,
+ (xpos_cuts.size()>0?&xpos_cuts:nullptr));
+ }
+ else if (t2)
+ {
+ My_timer timer; timer.start();
+ shds=new SHds;
+ std::ifstream inputfile(filename);
+ if (inputfile)
+ { inputfile>>*shds; }
+ inputfile.close();
+ timer.stop();
+ times.push_back(timer.time());
+
+ timer.reset(); timer.start();
+ shds->build_faces_array(); // For now faces are not saved, thus recomputed here
+ timer.stop();
+ times.push_back(timer.time());
+ }
+ else if (t3)
+ {
+ My_timer timer; timer.start();
+ shds=new SHds;
+ shds->load_streamed_files(filename);
+ timer.stop();
+ times.push_back(timer.time());
+
+ timer.reset(); timer.start();
+ shds->build_faces_array(); // For now faces are not saved, thus recomputed here
+ timer.stop();
+ times.push_back(timer.time());
+ }
+
+ shds->display_info();
+
+ // times[0] : global time
+ // times[1] : load and dispatch
+ // times[2] : compute all partial HDS
+ // times[3] : compute faces
+
+ if (t1)
+ {
+ print_txt_with_endl("[TIME]: to load and dispatch segments ", times[1]);
+ print_txt_with_endl("[TIME]: to compute all partial HDS ", times[2]);
+ print_txt_with_endl("[TIME]: to compute faces ", times[3]);
+ print_txt_with_endl("[GLOBAL]: time to compute the sHDS ",
+ times[0]);
+
+ if (save)
+ {
+ std::ofstream outputfile(filename+"output");
+ if (outputfile)
+ { outputfile<<*shds; }
+ outputfile.close();
+ }
+ }
+ else
+ {
+ print_txt_with_endl("[TIME]: to load the sHDS ", times[0]);
+ print_txt_with_endl("[TIME]: to compute faces ", times[1]);
+ }
+
+ if (traverse)
+ {
+ My_timer local_timer;
+ local_timer.start();
+ std::size_t nb=0;
+ for (int i=0; i<10; ++i)
+ { nb+=shds->traverse_all_faces(); }
+ local_timer.stop();
+ print_txt_with_endl("[TIME]: to traverse the sHDS 10 times ", local_timer.time(),
+ " (",nb/10," halfedges)");
+ }
+
+ if (xfig_file!="")
+ {
+ if (!export_xfig_sequence)
+ {
+ shds->export_to_xfig(xfig_file, !xfig_invert_y, xfig_x_size, 0,
+ color_segments, color_disks, color_criticals, color_borders,
+ width_segments, radius_disks, width_criticals,
+ width_borders, space_between_strips);
+ }
+ else
+ {
+ for (std::size_t i=0; i<shds->number_of_strips(); ++i)
+ {
+ boost::filesystem::path p1=xfig_file;
+ boost::filesystem::path p2=std::to_string(i)+std::string(".fig");
+ boost::filesystem::path p=p1.replace_extension(p2);
+
+ shds->export_to_xfig(p.string(), !xfig_invert_y, xfig_x_size, i+1,
+ color_segments, color_disks, color_criticals, color_borders,
+ width_segments, radius_disks, width_criticals,
+ width_borders, space_between_strips);
+ }
+ }
+ }
+
+ if (lcc)
+ {
+ My_timer lcc_timer;
+ lcc_timer.start();
+ LCC_2 lcc;
+ std::size_t nb_finite_faces=shds_to_lcc(*shds, lcc);
+ lcc_timer.stop();
+
+ lcc.display_characteristics(std::cout)<<", #finite-2-cells="<<nb_finite_faces
+ <<(lcc.is_valid()?", valid.":", NOT VALID.")<<std::endl;
+ print_txt_with_endl("[TIME]: to compute a LCC from the sHDS ", lcc_timer.time());
+ }
+
+ if (export_cgal)
+ {
+ typedef CGAL::Arr_segment_traits_2<EPECK> Traits_2;
+ typedef CGAL::Arrangement_2<Traits_2> Arrangement_2;
+ My_timer arr_timer;
+ arr_timer.start();
+ Arrangement_2 arr;
+ shds_to_cgal_arrangement(*shds, arr);
+ arr_timer.stop();
+ print_txt_with_endl("[TIME]: to export the sHDS into cgal arrangement ", arr_timer.time());
+ print_txt_with_endl("[GLOBAL]: time to build cgal arrangement from the sHDS ",
+ times[0]+arr_timer.time());
+ }
+
+ if (draw)
+ { draw_shds(*shds, "Viewer", false); }
+
+ delete shds;
+ return EXIT_SUCCESS;
+}
+///////////////////////////////////////////////////////////////////////////////