From 345434afe5610bc24feeb60854699d42233d22b3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Wed, 3 Oct 2018 15:56:18 +0200
Subject: [PATCH 01/35] Added a comment to clarify some code

---
 Mesh_3/include/CGAL/Mesh_3/Mesher_level.h | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesher_level.h b/Mesh_3/include/CGAL/Mesh_3/Mesher_level.h
index dc629bd5b20..2b5d5506068 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Mesher_level.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Mesher_level.h
@@ -1042,6 +1042,11 @@ public:
     {
       // Lock the element area on the grid
       Element element = derivd.extract_element_from_container_value(ce);
+
+      // This is safe to do with the concurrent compact container because even if the element `ce`
+      // gets removed from the TDS at this point, it is not actually deleted in the cells container and
+      // `ce->vertex(0-3)` still points to a vertex of the vertices container whose `.point()`
+      // can be safely accessed (even if that vertex has itself also been removed from the TDS).
       bool locked = derivd.try_lock_element(element, FIRST_GRID_LOCK_RADIUS);
 
       if( locked )

From e2b06830307fa7516421d5b198c5a10c1cf98264 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Wed, 3 Oct 2018 15:57:40 +0200
Subject: [PATCH 02/35] Fixed not checking 'hint' validity in the parallel
 insertion of a range of WPs

---
 .../include/CGAL/Regular_triangulation_3.h    | 77 +++++++++++++++++--
 1 file changed, 71 insertions(+), 6 deletions(-)

diff --git a/Triangulation_3/include/CGAL/Regular_triangulation_3.h b/Triangulation_3/include/CGAL/Regular_triangulation_3.h
index 8692371ebb2..0f0d27406c3 100644
--- a/Triangulation_3/include/CGAL/Regular_triangulation_3.h
+++ b/Triangulation_3/include/CGAL/Regular_triangulation_3.h
@@ -1324,14 +1324,48 @@ namespace CGAL {
 #endif
 
       Vertex_handle &hint = m_tls_hint.local();
+
       for( size_t i_point = r.begin() ; i_point != r.end() ; ++i_point)
       {
         bool success = false;
         const Weighted_point &p = m_points[i_point];
         while(!success)
         {
-          if (m_rt.try_lock_vertex(hint) && m_rt.try_lock_point(p))
+          // The 'hint' is unsafe to use immediately because we are in a regular triangulation,
+          // and the insertion of a (weighted) point in another thread might have hidden (deleted)
+          // the hint.
+#define CGAL_RT3_IGNORE_TDS_CONCEPT
+#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
+          if(!m_rt.tds().vertices().is_used(hint))
+#else
+          if(!m_rt.tds().is_vertex(hint))
+#endif
           {
+            hint = m_rt.finite_vertices_begin();
+            continue;
+          }
+
+          // We need to make sure that while are locking the position P1 := hint->point(), 'hint'
+          // does not get its position changed to P2 != P1.
+          const Weighted_point hint_point_mem = hint->point();
+
+          if(m_rt.try_lock_point(hint_point_mem) && m_rt.try_lock_point(p))
+          {
+            // Make sure that the hint is still valid (so that we can safely take hint->cell()) and
+            // that its position hasn't changed to ensure that we will start the locate from where
+            // we have locked.
+#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
+            if(!m_rt.tds().vertices().is_used(hint) ||
+#else
+            if(!m_rt.tds().is_vertex(hint) ||
+#endif
+               hint->point() != hint_point_mem)
+            {
+              hint = m_rt.finite_vertices_begin();
+              m_rt.unlock_all_elements();
+              continue;
+            }
+
             bool could_lock_zone;
             Locate_type lt;
             int li, lj;
@@ -1418,14 +1452,46 @@ namespace CGAL {
         const Weighted_point &p = m_points[i_point];
         while (!success)
         {
-          if (m_rt.try_lock_vertex(hint) && m_rt.try_lock_point(p))
+          // The 'hint' is unsafe to use immediately because we are in a regular triangulation,
+          // and the insertion of a (weighted) point in another thread might have hidden (deleted)
+          // the hint.
+#define CGAL_RT3_IGNORE_TDS_CONCEPT
+#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
+          if(!m_rt.tds().vertices().is_used(hint))
+#else
+          if(!m_rt.tds().is_vertex(hint))
+#endif
           {
+            hint = m_rt.finite_vertices_begin();
+            continue;
+          }
+
+          // We need to make sure that while are locking the position P1 := hint->point(), 'hint'
+          // does not get its position changed to P2 != P1.
+          const Weighted_point hint_point_mem = hint->point();
+
+          if(m_rt.try_lock_point(hint_point_mem) && m_rt.try_lock_point(p))
+          {
+            // Make sure that the hint is still valid (so that we can safely take hint->cell()) and
+            // that its position hasn't changed to ensure that we will start the locate from where
+            // we have locked.
+#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
+            if(!m_rt.tds().vertices().is_used(hint) ||
+#else
+            if(!m_rt.tds().is_vertex(hint) ||
+#endif
+               hint->point() != hint_point_mem)
+            {
+              hint = m_rt.finite_vertices_begin();
+              m_rt.unlock_all_elements();
+              continue;
+            }
+
             bool could_lock_zone;
             Locate_type lt;
             int li, lj;
 
-            Cell_handle c = m_rt.locate(p, lt, li, lj, hint->cell(),
-              &could_lock_zone);
+            Cell_handle c = m_rt.locate(p, lt, li, lj, hint->cell(), &could_lock_zone);
             Vertex_handle v;
             if (could_lock_zone)
               v = m_rt.insert(p, lt, c, li, lj, &could_lock_zone);
@@ -1504,8 +1570,7 @@ namespace CGAL {
         bool could_lock_zone, needs_to_be_done_sequentially;
         do
         {
-          needs_to_be_done_sequentially =
-            !m_rt.remove(v, &could_lock_zone);
+          needs_to_be_done_sequentially = !m_rt.remove(v, &could_lock_zone);
           m_rt.unlock_all_elements();
         } while (!could_lock_zone);
 

From e1122cef5c33b128ca6e2adbf2316679565c1d4b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Wed, 3 Oct 2018 15:58:59 +0200
Subject: [PATCH 03/35] Fixed deadlocks in parallel removal of a range of
 weighted points

---
 .../include/CGAL/Regular_triangulation_3.h    | 88 +++++++++++++++----
 1 file changed, 72 insertions(+), 16 deletions(-)

diff --git a/Triangulation_3/include/CGAL/Regular_triangulation_3.h b/Triangulation_3/include/CGAL/Regular_triangulation_3.h
index 0f0d27406c3..3eba767d11c 100644
--- a/Triangulation_3/include/CGAL/Regular_triangulation_3.h
+++ b/Triangulation_3/include/CGAL/Regular_triangulation_3.h
@@ -2453,40 +2453,96 @@ namespace CGAL {
 
   template < class Gt, class Tds, class Lds >
   bool
-    Regular_triangulation_3<Gt,Tds,Lds>::
-    remove(Vertex_handle v, bool *could_lock_zone)
+  Regular_triangulation_3<Gt,Tds,Lds>::
+  remove(Vertex_handle v, bool *could_lock_zone)
   {
     bool removed = true;
 
     // Locking vertex v...
-    if (!this->try_lock_vertex(v))
+    if(!this->try_lock_vertex(v))
     {
       *could_lock_zone = false;
     }
     else
     {
-      Vertex_handle hint = (v->cell()->vertex(0) == v ?
-        v->cell()->vertex(1) : v->cell()->vertex(0));
+      // Check that the vertex hasn't be deleted from the TDS while we were locking it
+#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
+      if(!tds().vertices().is_used(v))
+#else
+      if(!tds().is_vertex(v))
+#endif
+        return true; // vertex is already gone from the TDS, nothing to do
+
+      Vertex_handle hint = v->cell()->vertex(0) == v ? v->cell()->vertex(1) : v->cell()->vertex(0);
 
       Self tmp;
       Vertex_remover<Self> remover(tmp);
       removed = Tr_Base::remove(v, remover, could_lock_zone);
 
-      if (*could_lock_zone && removed)
+      if(*could_lock_zone && removed)
       {
-        // Re-insert the points that v was hiding.
-        for (typename Vertex_remover<Self>::Hidden_points_iterator
-          hi = remover.hidden_points_begin();
-          hi != remover.hidden_points_end(); ++hi)
+        // The vertex has been removed, try to re-insert the points that 'v' was hiding
+
+        // Start by unlocking the area of the removed vertex to avoid deadlocks
+        this->unlock_all_elements();
+
+        for(typename Vertex_remover<Self>::Hidden_points_iterator
+              hi = remover.hidden_points_begin();
+              hi != remover.hidden_points_end(); ++hi)
         {
-          bool could_lock_zone = false;
-          Vertex_handle hv;
-          while (!could_lock_zone)
+          const Weighted_point& wp = *hi;
+
+          // try to lock the positions of the hint and the hidden point
+          bool success = false;
+          while(!success)
           {
-            hv = insert (*hi, hint, &could_lock_zone);
+            // The 'hint' is unsafe to use immediately because we are in a regular triangulation,
+            // and the insertion of a (weighted) point in another thread might have hidden (deleted)
+            // the hint.
+#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
+            if(!tds().vertices().is_used(hint))
+#else
+            if(!tds().is_vertex(hint))
+#endif
+            {
+              hint = finite_vertices_begin();
+              continue;
+            }
+
+            // We need to make sure that while are locking the position P1 := hint->point(), 'hint'
+            // does not get its position changed to P2 != P1.
+            const Weighted_point hint_point_mem = hint->point();
+
+            if(this->try_lock_point(hint_point_mem) && this->try_lock_point(wp))
+            {
+              // Make sure that the hint is still valid (so that we can safely take hint->cell()) and
+              // that its position hasn't changed to ensure that we will start the locate from where
+              // we have locked.
+#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
+              if(!tds().vertices().is_used(hint) ||
+#else
+              if(!tds().is_vertex(hint) ||
+#endif
+                 hint->point() != hint_point_mem)
+              {
+                hint = finite_vertices_begin();
+                this->unlock_all_elements();
+                continue;
+              }
+
+              Vertex_handle hv = insert(wp, hint, could_lock_zone);
+
+              if(*could_lock_zone)
+              {
+                success = true;
+                if(hv != Vertex_handle())
+                  hint = hv;
+              }
+            }
+
+            // This unlocks everything in all cases: partial lock failure, failed insertion, successful insertion
+            this->unlock_all_elements();
           }
-          if (hv != Vertex_handle())
-            hint = hv;
         }
         CGAL_triangulation_expensive_postcondition (is_valid());
       }

From 8a7b11f34411ca20cb2ba4fd18ef3ef3520f5021 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Wed, 3 Oct 2018 16:05:18 +0200
Subject: [PATCH 04/35] Use actual weighted points in parallel regular
 triangulation tests

---
 .../CGAL/_test_cls_parallel_triangulation_3.h | 68 +++++++++++++++----
 .../test/Triangulation_3/test_regular_3.cpp   |  5 +-
 2 files changed, 57 insertions(+), 16 deletions(-)

diff --git a/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h b/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h
index 59abb7b6da2..80e5dbb9181 100644
--- a/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h
+++ b/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h
@@ -27,42 +27,82 @@
 
 #include <boost/mpl/if.hpp>
 
+template <typename Parallel_triangulation,
+          typename WeightedTag = typename Parallel_triangulation::Weighted_tag>
+struct PTR_random_pts_generator
+{
+  typedef typename Parallel_triangulation::Point       Point;
+
+  void operator()(const int num, CGAL::Random& rnd, std::vector<Point>& points) const
+  {
+    CGAL::Random_points_in_cube_3<Point> gen(1., rnd);
+
+    points.reserve(num);
+    for(int i=0; i!=num; ++i)
+      points.push_back(*gen++);
+  }
+};
+
+template <typename Parallel_triangulation>
+struct PTR_random_pts_generator<Parallel_triangulation, CGAL::Tag_true /*weighted triangulation*/>
+{
+  typedef typename Parallel_triangulation::Bare_point      Bare_point;
+  typedef typename Parallel_triangulation::Weighted_point  Weighted_point;
+
+  void operator()(const int num, CGAL::Random& rnd, std::vector<Weighted_point>& points) const
+  {
+    CGAL::Random_points_in_cube_3<Bare_point> gen(1., rnd);
+
+    points.reserve(num);
+    for(int i=0; i!=num; ++i)
+      points.push_back(Weighted_point(*gen++, rnd.get_double(-1., 1.)));
+  }
+};
+
 template <class Parallel_triangulation>
 void
 _test_cls_parallel_triangulation_3(const Parallel_triangulation &)
 {
-  const int NUM_INSERTED_POINTS = 5000;
-
   typedef Parallel_triangulation                                Cls;
 
   typedef typename Cls::Vertex_handle                           Vertex_handle;
   typedef typename Cls::Point                                   Point;
 
-  CGAL::Random_points_in_cube_3<Point> rnd(1.);
+  typedef std::vector<Point>                                    Point_container;
+
+  CGAL::Random rnd;
+  std::cout << "Seed: " << rnd.get_seed() << std::endl;
+
+  const int num_insert = 5000;
+  Point_container points;
+  PTR_random_pts_generator<Parallel_triangulation> points_gen;
+  points_gen(num_insert, rnd, points);
 
-  // Construction from a vector of points
-  std::vector<Point> points;
-  points.reserve(NUM_INSERTED_POINTS);
-  for (int i = 0; i != NUM_INSERTED_POINTS; ++i)
-    points.push_back(*rnd++);
-  
   // Construct the locking data-structure, using the bounding-box of the points
-  typename Cls::Lock_data_structure locking_ds(
-    CGAL::Bbox_3(-1., -1., -1., 1., 1., 1.), 50);
+  typename Cls::Lock_data_structure locking_ds(CGAL::Bbox_3(-1., -1., -1., 1., 1., 1.), 50);
+
   // Contruct the triangulation in parallel
   std::cout << "Construction and parallel insertion" << std::endl;
   Cls tr(points.begin(), points.end(), &locking_ds);
+  std::cout << "Triangulation has " << tr.number_of_vertices() << " vertices" << std::endl;
+  assert(tr.is_valid());
 
   std::cout << "Parallel removal" << std::endl;
-  // Remove the first 100,000 vertices
   std::vector<Vertex_handle> vertices_to_remove;
   typename Cls::Finite_vertices_iterator vit = tr.finite_vertices_begin();
-  for (int i = 0 ; i < NUM_INSERTED_POINTS/10 ; ++i)
+
+  const int num_remove = tr.number_of_vertices() / 10;
+  std::cout << "Removing " << num_remove << " from " << tr.number_of_vertices() << " vertices" << std::endl;
+
+  for(int i=0 ; i<num_remove; ++i) {
     vertices_to_remove.push_back(vit++);
+}
+
   // Parallel remove
   tr.remove(vertices_to_remove.begin(), vertices_to_remove.end());
-  
+  std::cout << "Now, " << tr.number_of_vertices() << " vertices are left" << std::endl;
   assert(tr.is_valid());
+
   tr.clear();
   assert(tr.is_valid());
   assert(tr.dimension()==-1);
diff --git a/Triangulation_3/test/Triangulation_3/test_regular_3.cpp b/Triangulation_3/test/Triangulation_3/test_regular_3.cpp
index af6257bd002..60d0423b248 100644
--- a/Triangulation_3/test/Triangulation_3/test_regular_3.cpp
+++ b/Triangulation_3/test/Triangulation_3/test_regular_3.cpp
@@ -448,9 +448,10 @@ int main()
     CGAL::Parallel_tag >	                            Tds_parallel;
   typedef CGAL::Regular_triangulation_3<
     traits, Tds_parallel, Lock_ds>                    RT_parallel;
-  // The following test won't do things in parallel since it doesn't provide
-  // a lock data structure
+
+  // The following test won't do things in parallel since it doesn't provide a lock data structure
   test_RT<RT_parallel>();
+
   // This test performs parallel operations
   _test_cls_parallel_triangulation_3( RT_parallel() );
 #endif

From 8500add791c314c7399568583e87142d0ffc5551 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Wed, 3 Oct 2018 16:31:42 +0200
Subject: [PATCH 05/35] There is no 'try'!

---
 Triangulation_3/include/CGAL/Regular_triangulation_3.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Triangulation_3/include/CGAL/Regular_triangulation_3.h b/Triangulation_3/include/CGAL/Regular_triangulation_3.h
index 3eba767d11c..e8cc970b0fa 100644
--- a/Triangulation_3/include/CGAL/Regular_triangulation_3.h
+++ b/Triangulation_3/include/CGAL/Regular_triangulation_3.h
@@ -2481,7 +2481,7 @@ namespace CGAL {
 
       if(*could_lock_zone && removed)
       {
-        // The vertex has been removed, try to re-insert the points that 'v' was hiding
+        // The vertex has been removed, re-insert the points that 'v' was hiding
 
         // Start by unlocking the area of the removed vertex to avoid deadlocks
         this->unlock_all_elements();

From aa1e0acbeef802a1ae8cbdfeaec2f0808211d4d5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Wed, 10 Oct 2018 13:08:14 +0200
Subject: [PATCH 06/35] Replaced macro with a legal vertex handle validity
 checker

---
 .../include/CGAL/Regular_triangulation_3.h    | 68 +++++++++----------
 1 file changed, 31 insertions(+), 37 deletions(-)

diff --git a/Triangulation_3/include/CGAL/Regular_triangulation_3.h b/Triangulation_3/include/CGAL/Regular_triangulation_3.h
index e8cc970b0fa..8c01c687981 100644
--- a/Triangulation_3/include/CGAL/Regular_triangulation_3.h
+++ b/Triangulation_3/include/CGAL/Regular_triangulation_3.h
@@ -1291,6 +1291,25 @@ namespace CGAL {
       }
     };
 
+    // In parallel operations, we need to be able to check the health of the 'hint' vertex handle,
+    // which might be invalided by other threads. One way to do that is the 'is_vertex()' function
+    // of the TDS, but it runs in O(sqrt(n)) complexity. When we are using our TDS, we can use
+    // a lower level function from the compact container, which runs in constant time.
+    template <typename TDS_>
+    struct Vertex_validity_checker
+    {
+      bool operator()(const typename TDS_::Vertex_handle vh_, const TDS_& tds_) { return tds_.is_vertex(vh_); }
+    };
+
+    template <typename T1, typename T2, typename T3>
+    struct Vertex_validity_checker<CGAL::Triangulation_data_structure_3<T1, T2, T3> >
+    {
+      typedef CGAL::Triangulation_data_structure_3<T1, T2, T3>   TDS_;
+
+      bool operator()(const typename TDS_::Vertex_handle vh_, const TDS_& tds_) {
+        return tds_.vertices().is_used(vh_); }
+    };
+
   // Functor for parallel insert(begin, end) function
   template <typename RT>
   class Insert_point
@@ -1324,6 +1343,7 @@ namespace CGAL {
 #endif
 
       Vertex_handle &hint = m_tls_hint.local();
+      Vertex_validity_checker<typename RT::Triangulation_data_structure> vertex_validity_check;
 
       for( size_t i_point = r.begin() ; i_point != r.end() ; ++i_point)
       {
@@ -1334,12 +1354,7 @@ namespace CGAL {
           // The 'hint' is unsafe to use immediately because we are in a regular triangulation,
           // and the insertion of a (weighted) point in another thread might have hidden (deleted)
           // the hint.
-#define CGAL_RT3_IGNORE_TDS_CONCEPT
-#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
-          if(!m_rt.tds().vertices().is_used(hint))
-#else
-          if(!m_rt.tds().is_vertex(hint))
-#endif
+          if(!vertex_validity_check(hint, m_rt.tds()))
           {
             hint = m_rt.finite_vertices_begin();
             continue;
@@ -1354,11 +1369,7 @@ namespace CGAL {
             // Make sure that the hint is still valid (so that we can safely take hint->cell()) and
             // that its position hasn't changed to ensure that we will start the locate from where
             // we have locked.
-#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
-            if(!m_rt.tds().vertices().is_used(hint) ||
-#else
-            if(!m_rt.tds().is_vertex(hint) ||
-#endif
+            if(!vertex_validity_check(hint, m_rt.tds()) ||
                hint->point() != hint_point_mem)
             {
               hint = m_rt.finite_vertices_begin();
@@ -1445,6 +1456,8 @@ namespace CGAL {
 #endif
 
       Vertex_handle &hint = m_tls_hint.local();
+      Vertex_validity_checker<typename RT::Triangulation_data_structure> vertex_validity_check;
+
       for (size_t i_idx = r.begin() ; i_idx != r.end() ; ++i_idx)
       {
         bool success = false;
@@ -1455,12 +1468,7 @@ namespace CGAL {
           // The 'hint' is unsafe to use immediately because we are in a regular triangulation,
           // and the insertion of a (weighted) point in another thread might have hidden (deleted)
           // the hint.
-#define CGAL_RT3_IGNORE_TDS_CONCEPT
-#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
-          if(!m_rt.tds().vertices().is_used(hint))
-#else
-          if(!m_rt.tds().is_vertex(hint))
-#endif
+          if(!vertex_validity_check(hint, m_rt.tds()))
           {
             hint = m_rt.finite_vertices_begin();
             continue;
@@ -1475,11 +1483,7 @@ namespace CGAL {
             // Make sure that the hint is still valid (so that we can safely take hint->cell()) and
             // that its position hasn't changed to ensure that we will start the locate from where
             // we have locked.
-#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
-            if(!m_rt.tds().vertices().is_used(hint) ||
-#else
-            if(!m_rt.tds().is_vertex(hint) ||
-#endif
+            if(!vertex_validity_check(hint, m_rt.tds()) ||
                hint->point() != hint_point_mem)
             {
               hint = m_rt.finite_vertices_begin();
@@ -2465,12 +2469,10 @@ namespace CGAL {
     }
     else
     {
+      Vertex_validity_checker<Tds> vertex_validity_check;
+
       // Check that the vertex hasn't be deleted from the TDS while we were locking it
-#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
-      if(!tds().vertices().is_used(v))
-#else
-      if(!tds().is_vertex(v))
-#endif
+      if(!vertex_validity_check(v, tds()))
         return true; // vertex is already gone from the TDS, nothing to do
 
       Vertex_handle hint = v->cell()->vertex(0) == v ? v->cell()->vertex(1) : v->cell()->vertex(0);
@@ -2499,11 +2501,7 @@ namespace CGAL {
             // The 'hint' is unsafe to use immediately because we are in a regular triangulation,
             // and the insertion of a (weighted) point in another thread might have hidden (deleted)
             // the hint.
-#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
-            if(!tds().vertices().is_used(hint))
-#else
-            if(!tds().is_vertex(hint))
-#endif
+            if(!vertex_validity_check(hint, tds()))
             {
               hint = finite_vertices_begin();
               continue;
@@ -2518,11 +2516,7 @@ namespace CGAL {
               // Make sure that the hint is still valid (so that we can safely take hint->cell()) and
               // that its position hasn't changed to ensure that we will start the locate from where
               // we have locked.
-#ifdef CGAL_RT3_IGNORE_TDS_CONCEPT
-              if(!tds().vertices().is_used(hint) ||
-#else
-              if(!tds().is_vertex(hint) ||
-#endif
+              if(!vertex_validity_check(hint, tds()) ||
                  hint->point() != hint_point_mem)
               {
                 hint = finite_vertices_begin();

From 7e75a8a243b629b6fbfbf8c563afe76cdf6d1299 Mon Sep 17 00:00:00 2001
From: Mael <mael.rouxel.labbe@geometryfactory.com>
Date: Wed, 17 Oct 2018 09:24:46 +0200
Subject: [PATCH 07/35] Fixed indentation (Thanks for the suggestion @MaelRL !)

---
 .../include/CGAL/_test_cls_parallel_triangulation_3.h           | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h b/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h
index 80e5dbb9181..2b1314a4066 100644
--- a/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h
+++ b/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h
@@ -96,7 +96,7 @@ _test_cls_parallel_triangulation_3(const Parallel_triangulation &)
 
   for(int i=0 ; i<num_remove; ++i) {
     vertices_to_remove.push_back(vit++);
-}
+  }
 
   // Parallel remove
   tr.remove(vertices_to_remove.begin(), vertices_to_remove.end());

From bb19c965f9e39ca067fa64c3a0aa4c0e59d32b9c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Thu, 11 Oct 2018 10:05:37 +0200
Subject: [PATCH 08/35] Moved helper struct outside of TBB macros to fix
 sequential compilation

---
 .../include/CGAL/Regular_triangulation_3.h    | 38 +++++++++----------
 1 file changed, 19 insertions(+), 19 deletions(-)

diff --git a/Triangulation_3/include/CGAL/Regular_triangulation_3.h b/Triangulation_3/include/CGAL/Regular_triangulation_3.h
index 8c01c687981..36ea66b4765 100644
--- a/Triangulation_3/include/CGAL/Regular_triangulation_3.h
+++ b/Triangulation_3/include/CGAL/Regular_triangulation_3.h
@@ -763,6 +763,25 @@ namespace CGAL {
       return std::copy(vertices.begin(), vertices.end(), res);
     }
 
+    // In parallel operations, we need to be able to check the health of the 'hint' vertex handle,
+    // which might be invalided by other threads. One way to do that is the 'is_vertex()' function
+    // of the TDS, but it runs in O(sqrt(n)) complexity. When we are using our TDS, we can use
+    // a lower level function from the compact container, which runs in constant time.
+    template <typename TDS_>
+    struct Vertex_validity_checker
+    {
+      bool operator()(const typename TDS_::Vertex_handle vh_, const TDS_& tds_) { return tds_.is_vertex(vh_); }
+    };
+
+    template <typename T1, typename T2, typename T3>
+    struct Vertex_validity_checker<CGAL::Triangulation_data_structure_3<T1, T2, T3> >
+    {
+      typedef CGAL::Triangulation_data_structure_3<T1, T2, T3>   TDS_;
+
+      bool operator()(const typename TDS_::Vertex_handle vh_, const TDS_& tds_) {
+        return tds_.vertices().is_used(vh_); }
+    };
+
     void remove (Vertex_handle v);
     // Concurrency-safe
     // See Triangulation_3::remove for more information
@@ -1291,25 +1310,6 @@ namespace CGAL {
       }
     };
 
-    // In parallel operations, we need to be able to check the health of the 'hint' vertex handle,
-    // which might be invalided by other threads. One way to do that is the 'is_vertex()' function
-    // of the TDS, but it runs in O(sqrt(n)) complexity. When we are using our TDS, we can use
-    // a lower level function from the compact container, which runs in constant time.
-    template <typename TDS_>
-    struct Vertex_validity_checker
-    {
-      bool operator()(const typename TDS_::Vertex_handle vh_, const TDS_& tds_) { return tds_.is_vertex(vh_); }
-    };
-
-    template <typename T1, typename T2, typename T3>
-    struct Vertex_validity_checker<CGAL::Triangulation_data_structure_3<T1, T2, T3> >
-    {
-      typedef CGAL::Triangulation_data_structure_3<T1, T2, T3>   TDS_;
-
-      bool operator()(const typename TDS_::Vertex_handle vh_, const TDS_& tds_) {
-        return tds_.vertices().is_used(vh_); }
-    };
-
   // Functor for parallel insert(begin, end) function
   template <typename RT>
   class Insert_point

From 82e458da80106dadbd7eae6b230fc699ed9b6172 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Wed, 14 Nov 2018 07:57:07 +0100
Subject: [PATCH 09/35] Replaced TDS template specialization with tag-based
 detection

---
 .../CGAL/Triangulation_data_structure_3.h     |  6 ++++-
 .../include/CGAL/Regular_triangulation_3.h    | 25 ++++++++++++++-----
 2 files changed, 24 insertions(+), 7 deletions(-)

diff --git a/TDS_3/include/CGAL/Triangulation_data_structure_3.h b/TDS_3/include/CGAL/Triangulation_data_structure_3.h
index 09981b25aa3..54b0a7cb213 100644
--- a/TDS_3/include/CGAL/Triangulation_data_structure_3.h
+++ b/TDS_3/include/CGAL/Triangulation_data_structure_3.h
@@ -77,9 +77,13 @@ class Triangulation_data_structure_3
   typedef Triangulation_data_structure_3<Vb, Cb, Concurrency_tag_> Tds;
 
 public:
-
   typedef Concurrency_tag_            Concurrency_tag;
 
+  // This tag is used in the parallel operations of RT_3 to access some functions
+  // of the TDS (tds.vertices().is_used(Vertex_handle)) that are much more efficient
+  // than what is exposed by the TDS concept (tds.is_vertex(Vertex_handle)).
+  typedef CGAL::Tag_true              Is_CGAL_TDS_3;
+
   // Tools to change the Vertex and Cell types of the TDS.
   template < typename Vb2 >
   struct Rebind_vertex {
diff --git a/Triangulation_3/include/CGAL/Regular_triangulation_3.h b/Triangulation_3/include/CGAL/Regular_triangulation_3.h
index 36ea66b4765..a351e78ae33 100644
--- a/Triangulation_3/include/CGAL/Regular_triangulation_3.h
+++ b/Triangulation_3/include/CGAL/Regular_triangulation_3.h
@@ -767,19 +767,32 @@ namespace CGAL {
     // which might be invalided by other threads. One way to do that is the 'is_vertex()' function
     // of the TDS, but it runs in O(sqrt(n)) complexity. When we are using our TDS, we can use
     // a lower level function from the compact container, which runs in constant time.
+    BOOST_MPL_HAS_XXX_TRAIT_DEF(Is_CGAL_TDS_3)
+
+    template <typename TDS_,
+              bool has_TDS_tag = has_Is_CGAL_TDS_3<TDS_>::value>
+    struct Is_CGAL_TDS_3 : public CGAL::Tag_false
+    { };
+
     template <typename TDS_>
+    struct Is_CGAL_TDS_3<TDS_, true> : public CGAL::Boolean_tag<TDS_::Is_CGAL_TDS_3::value>
+    { };
+
+    template <typename TDS_,
+              bool is_CGAL_TDS_3 = Is_CGAL_TDS_3<TDS_>::value>
     struct Vertex_validity_checker
     {
-      bool operator()(const typename TDS_::Vertex_handle vh_, const TDS_& tds_) { return tds_.is_vertex(vh_); }
+      bool operator()(const typename TDS_::Vertex_handle vh_, const TDS_& tds_) {
+        return tds_.is_vertex(vh_);
+      }
     };
 
-    template <typename T1, typename T2, typename T3>
-    struct Vertex_validity_checker<CGAL::Triangulation_data_structure_3<T1, T2, T3> >
+    template <typename TDS_>
+    struct Vertex_validity_checker<TDS_, true /* is_CGAL_TDS_3 */>
     {
-      typedef CGAL::Triangulation_data_structure_3<T1, T2, T3>   TDS_;
-
       bool operator()(const typename TDS_::Vertex_handle vh_, const TDS_& tds_) {
-        return tds_.vertices().is_used(vh_); }
+        return tds_.vertices().is_used(vh_);
+      }
     };
 
     void remove (Vertex_handle v);

From ab67b0cc34b6dfddaee7e817a1d67661bdfb757f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Mon, 19 Nov 2018 16:00:29 +0100
Subject: [PATCH 10/35] Fixed conversion warning

---
 .../include/CGAL/_test_cls_parallel_triangulation_3.h         | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h b/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h
index 2b1314a4066..facbd70f2d7 100644
--- a/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h
+++ b/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h
@@ -91,10 +91,10 @@ _test_cls_parallel_triangulation_3(const Parallel_triangulation &)
   std::vector<Vertex_handle> vertices_to_remove;
   typename Cls::Finite_vertices_iterator vit = tr.finite_vertices_begin();
 
-  const int num_remove = tr.number_of_vertices() / 10;
+  const std::size_t num_remove = tr.number_of_vertices() / 10;
   std::cout << "Removing " << num_remove << " from " << tr.number_of_vertices() << " vertices" << std::endl;
 
-  for(int i=0 ; i<num_remove; ++i) {
+  for(std::size_t i=0 ; i<num_remove; ++i) {
     vertices_to_remove.push_back(vit++);
   }
 

From 6d6c2f371736e5ee8f0732d65e16068f9cf253b3 Mon Sep 17 00:00:00 2001
From: Efi Fogel <efifogel@gmail.com>
Date: Fri, 30 Nov 2018 10:09:31 +0100
Subject: [PATCH 11/35] Fixed initializer with workaround for the absense of
 c++11 auto feature--added facets properly.

---
 .../Arr_polyhedral_sgm.h                      | 19 +++++++++++++------
 1 file changed, 13 insertions(+), 6 deletions(-)

diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h b/Arrangement_on_surface_2/include/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h
index b2be5687c10..603837c1ebf 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h
@@ -242,6 +242,14 @@ private:
     /*! Set the marked-face index */
     void set_marked_facet_index(size_type id) {m_marked_facet_index = id;}
 
+    /*! Add all facets. */
+    template <typename Iterator, typename Builder>
+    void add_facets(Iterator begin, Iterator end, Builder& B)
+    {
+      for (Iterator it = begin; it != end; ++it) B.add_vertex_to_facet(*it);
+      B.end_facet();
+    }
+
     /*! builds the polyhedron */
     void operator()(HDS& hds)
     {
@@ -262,12 +270,11 @@ private:
       for (CoordIndexIter it = m_indices_begin; it != m_indices_end; ++it) {
         Polyhedron_facet_handle fh = B.begin_facet();
         if (counter == m_marked_facet_index) fh->set_marked(true);
-        //! \todo EF: when upgrading to C++11 change the type of the following
-        // iterator to auto. Better yet, use for (auto blah : foo).
-        for (std::vector<size_t>::const_iterator iit = it->begin();
-             iit != it->end(); ++iit)
-          B.add_vertex_to_facet(*iit);
-        B.end_facet();
+        //! \todo EF: when upgrading to C++11 enable the following code and
+        // remove add_facets().
+        // for (const auto& facet : *it) B.add_vertex_to_facet(facet);
+        // B.end_facet();
+        add_facets(it->begin(), it->end(), B);
         ++counter;
       }
       B.end_surface();

From 6a5857cf6781641d05cb2cef175b39a430018c1b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= <sebastien.loriot@cgal.org>
Date: Mon, 10 Dec 2018 07:22:17 +0100
Subject: [PATCH 12/35] add missing :

---
 .../Straight_skeleton_2/CGAL/Straight_skeleton_converter_2.h    | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Straight_skeleton_2/doc/Straight_skeleton_2/CGAL/Straight_skeleton_converter_2.h b/Straight_skeleton_2/doc/Straight_skeleton_2/CGAL/Straight_skeleton_converter_2.h
index 0497ab07c75..0c26fe5caac 100644
--- a/Straight_skeleton_2/doc/Straight_skeleton_2/CGAL/Straight_skeleton_converter_2.h
+++ b/Straight_skeleton_2/doc/Straight_skeleton_2/CGAL/Straight_skeleton_converter_2.h
@@ -100,7 +100,7 @@ boost::shared_ptr<Target_skeleton> operator()( Source_skeleton const& s) const;
 \tparam SrcSs type of the source straight skeleton
 \tparam TgtSs type of the target straight skeleton
 \tparam NTConverter  a function object that must 
-provide `TgtSs:Traits::FT operator()(SrcSs::Traits::FT n)` that converts `n` to an 
+provide `TgtSs::Traits::FT operator()(SrcSs::Traits::FT n)` that converts `n` to an 
 `TgtSs::Traits::FT` which has the same value. The default value of this parameter is `NT_converter<SrcSs::Traits::FT, TgtSs::Traits::FT>`. 
 
 \cgalModels `StraightSkeletonItemsConverter_2`

From a6c1632f82e40db6626aa95da59e6d59a42e8375 Mon Sep 17 00:00:00 2001
From: Andreas Fabri <andreas.fabri@geometryfactory.com>
Date: Thu, 13 Dec 2018 09:51:47 +0100
Subject: [PATCH 13/35] No need to disable warnings

---
 .../CMakeLists.txt                              | 17 -----------------
 1 file changed, 17 deletions(-)

diff --git a/Surface_mesh_parameterization/test/Surface_mesh_parameterization/CMakeLists.txt b/Surface_mesh_parameterization/test/Surface_mesh_parameterization/CMakeLists.txt
index 9455b197844..e325bf4aafa 100644
--- a/Surface_mesh_parameterization/test/Surface_mesh_parameterization/CMakeLists.txt
+++ b/Surface_mesh_parameterization/test/Surface_mesh_parameterization/CMakeLists.txt
@@ -14,24 +14,7 @@ find_package(CGAL QUIET)
 
 if ( CGAL_FOUND )
 
-  # VisualC++ optimization for applications dealing with large data
-  if (MSVC)
 
-    # Allow Windows applications to use up to 3GB of RAM
-    SET (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} /LARGEADDRESSAWARE")
-
-    # Turn off stupid VC++ warnings
-    SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4267 /wd4311 /wd4800 /wd4503 /wd4244 /wd4345 /wd4996 /wd4396 /wd4018")
-    
-    # Print new compilation options
-    message( STATUS "USING DEBUG CXXFLAGS   = '${CMAKE_CXX_FLAGS} ${CMAKE_CXX_FLAGS_DEBUG}'" )
-    message( STATUS "USING DEBUG EXEFLAGS   = '${CMAKE_EXE_LINKER_FLAGS} ${CMAKE_EXE_LINKER_FLAGS_DEBUG}'" )
-    message( STATUS "USING RELEASE CXXFLAGS = '${CMAKE_CXX_FLAGS} ${CMAKE_CXX_FLAGS_RELEASE}'" )
-    message( STATUS "USING RELEASE EXEFLAGS = '${CMAKE_EXE_LINKER_FLAGS} ${CMAKE_EXE_LINKER_FLAGS_RELEASE}'" )
-  endif(MSVC)
-
-
-  #find Eigen 
   find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater)
   if(EIGEN3_FOUND)
     include( ${EIGEN3_USE_FILE} )

From 858d10070b079fa9c6818b0b6bd57d12c6995939 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Mon, 17 Dec 2018 13:07:18 +0100
Subject: [PATCH 14/35] Improve PLY element handling and thus fix bug when
 reading unknown element

---
 Point_set_3/include/CGAL/Point_set_3/IO.h     | 219 ++++++------
 .../include/CGAL/IO/read_ply_points.h         | 317 +++++++++---------
 Polyhedron_IO/include/CGAL/IO/PLY_reader.h    | 254 ++++++++------
 3 files changed, 429 insertions(+), 361 deletions(-)

diff --git a/Point_set_3/include/CGAL/Point_set_3/IO.h b/Point_set_3/include/CGAL/Point_set_3/IO.h
index 93117799195..287cf45e728 100644
--- a/Point_set_3/include/CGAL/Point_set_3/IO.h
+++ b/Point_set_3/include/CGAL/Point_set_3/IO.h
@@ -135,7 +135,7 @@ private:
   struct Abstract_ply_property_to_point_set_property
   {
     virtual ~Abstract_ply_property_to_point_set_property() { }
-    virtual void assign (PLY_reader& reader, typename Point_set::Index index) = 0;
+    virtual void assign (PLY_element& element, typename Point_set::Index index) = 0;
   };
 
   template <typename Type>
@@ -154,10 +154,10 @@ private:
       m_pmap = ps.push_property_map (m_map);
     }
     
-    virtual void assign (PLY_reader& reader, typename Point_set::Index index)
+    virtual void assign (PLY_element& element, typename Point_set::Index index)
     {
       Type t;
-      reader.assign (t, m_name.c_str());
+      element.assign (t, m_name.c_str());
       put(m_pmap, index, t);
     }
   };
@@ -178,123 +178,121 @@ public:
       delete m_properties[i];
   }
 
-  void instantiate_properties  (PLY_reader& reader)
+  void instantiate_properties  (PLY_element& element)
   {
-    const std::vector<PLY_read_number*>& readers
-      = reader.readers();
-
     bool has_normal[3] = { false, false, false };
     
-    for (std::size_t i = 0; i < readers.size(); ++ i)
+    for (std::size_t j = 0; j < element.number_of_properties(); ++ j)
+    {
+      internal::PLY::PLY_read_number* property = element.property(j);
+        
+      const std::string& name = property->name();
+      if (name == "x" ||
+          name == "y" ||
+          name == "z")
       {
-        const std::string& name = readers[i]->name();
-        if (name == "x" ||
-            name == "y" ||
-            name == "z")
-          {
-            if (dynamic_cast<PLY_read_typed_number<float>*>(readers[i]))
-              m_use_floats = true;
-            continue;
-          }
-        if (name == "nx")
-          {
-            has_normal[0] = true;
-            continue;
-          }
-        if (name == "ny")
-          {
-            has_normal[1] = true;
-            continue;
-          }
-        if (name == "nz")
-          {
-            has_normal[2] = true;
-            continue;
-          }
-
-        if (dynamic_cast<PLY_read_typed_number<boost::int8_t>*>(readers[i]))
-          {
-            m_properties.push_back
-              (new PLY_property_to_point_set_property<boost::int8_t>(m_point_set,
-                                                                     name));
-          }
-        else if (dynamic_cast<PLY_read_typed_number<boost::uint8_t>*>(readers[i]))
-          {
-            m_properties.push_back
-              (new PLY_property_to_point_set_property<boost::uint8_t>(m_point_set,
-                                                                      name));
-          }
-        else if (dynamic_cast<PLY_read_typed_number<boost::int16_t>*>(readers[i]))
-          {
-            m_properties.push_back
-              (new PLY_property_to_point_set_property<boost::int16_t>(m_point_set,
-                                                                      name));
-          }
-        else if (dynamic_cast<PLY_read_typed_number<boost::uint16_t>*>(readers[i]))
-          {
-            m_properties.push_back
-              (new PLY_property_to_point_set_property<boost::uint16_t>(m_point_set,
-                                                                       name));
-          }
-        else if (dynamic_cast<PLY_read_typed_number<boost::int32_t>*>(readers[i]))
-          {
-            m_properties.push_back
-              (new PLY_property_to_point_set_property<boost::int32_t>(m_point_set,
-                                                                      name));
-          }
-        else if (dynamic_cast<PLY_read_typed_number<boost::uint32_t>*>(readers[i]))
-          {
-            m_properties.push_back
-              (new PLY_property_to_point_set_property<boost::uint32_t>(m_point_set,
-                                                                       name));
-          }
-        else if (dynamic_cast<PLY_read_typed_number<float>*>(readers[i]))
-          {
-            m_properties.push_back
-              (new PLY_property_to_point_set_property<float>(m_point_set,
-                                                             name));
-          }
-        else if (dynamic_cast<PLY_read_typed_number<double>*>(readers[i]))
-          {
-            m_properties.push_back
-              (new PLY_property_to_point_set_property<double>(m_point_set,
-                                                             name));
-          }
+        if (dynamic_cast<PLY_read_typed_number<float>*>(property))
+          m_use_floats = true;
+        continue;
+      }
+      if (name == "nx")
+      {
+        has_normal[0] = true;
+        continue;
+      }
+      if (name == "ny")
+      {
+        has_normal[1] = true;
+        continue;
+      }
+      if (name == "nz")
+      {
+        has_normal[2] = true;
+        continue;
       }
 
+      if (dynamic_cast<PLY_read_typed_number<boost::int8_t>*>(property))
+      {
+        m_properties.push_back
+          (new PLY_property_to_point_set_property<boost::int8_t>(m_point_set,
+                                                                 name));
+      }
+      else if (dynamic_cast<PLY_read_typed_number<boost::uint8_t>*>(property))
+      {
+        m_properties.push_back
+          (new PLY_property_to_point_set_property<boost::uint8_t>(m_point_set,
+                                                                  name));
+      }
+      else if (dynamic_cast<PLY_read_typed_number<boost::int16_t>*>(property))
+      {
+        m_properties.push_back
+          (new PLY_property_to_point_set_property<boost::int16_t>(m_point_set,
+                                                                  name));
+      }
+      else if (dynamic_cast<PLY_read_typed_number<boost::uint16_t>*>(property))
+      {
+        m_properties.push_back
+          (new PLY_property_to_point_set_property<boost::uint16_t>(m_point_set,
+                                                                   name));
+      }
+      else if (dynamic_cast<PLY_read_typed_number<boost::int32_t>*>(property))
+      {
+        m_properties.push_back
+          (new PLY_property_to_point_set_property<boost::int32_t>(m_point_set,
+                                                                  name));
+      }
+      else if (dynamic_cast<PLY_read_typed_number<boost::uint32_t>*>(property))
+      {
+        m_properties.push_back
+          (new PLY_property_to_point_set_property<boost::uint32_t>(m_point_set,
+                                                                   name));
+      }
+      else if (dynamic_cast<PLY_read_typed_number<float>*>(property))
+      {
+        m_properties.push_back
+          (new PLY_property_to_point_set_property<float>(m_point_set,
+                                                         name));
+      }
+      else if (dynamic_cast<PLY_read_typed_number<double>*>(property))
+      {
+        m_properties.push_back
+          (new PLY_property_to_point_set_property<double>(m_point_set,
+                                                          name));
+      }
+    }
     if (has_normal[0] && has_normal[1] && has_normal[2])
       m_point_set.add_normal_map();
   }
   
-  void process_line (PLY_reader& reader)
+  void process_line (PLY_element& element)
   {
     m_point_set.insert();
     
     if (m_use_floats)
-      process_line<float>(reader);
+      process_line<float>(element);
     else
-      process_line<double>(reader);
+      process_line<double>(element);
 
     for (std::size_t i = 0; i < m_properties.size(); ++ i)
-      m_properties[i]->assign (reader, *(m_point_set.end() - 1));
+      m_properties[i]->assign (element, *(m_point_set.end() - 1));
   }
 
   template <typename FT>
-  void process_line (PLY_reader& reader)
+  void process_line (PLY_element& element)
   {
     FT x = (FT)0.,y = (FT)0., z = (FT)0.,
       nx = (FT)0., ny = (FT)0., nz = (FT)0.;
-    reader.assign (x, "x");
-    reader.assign (y, "y");
-    reader.assign (z, "z");
+    element.assign (x, "x");
+    element.assign (y, "y");
+    element.assign (z, "z");
     Point point (x, y, z);
     m_point_set.point(*(m_point_set.end() - 1)) = point;
 
     if (m_point_set.has_normal_map())
       {
-        reader.assign (nx, "nx");
-        reader.assign (ny, "ny");
-        reader.assign (nz, "nz");
+        element.assign (nx, "nx");
+        element.assign (ny, "ny");
+        element.assign (nz, "nz");
         Vector normal (nx, ny, nz);
         m_point_set.normal(*(m_point_set.end() - 1)) = normal;
       }
@@ -401,24 +399,33 @@ read_ply_point_set(
   if (comments != NULL)
     *comments = reader.comments();
 
-  filler.instantiate_properties (reader);
+  for (std::size_t i = 0; i < reader.number_of_elements(); ++ i)
+  {
+    internal::PLY::PLY_element& element = reader.element(i);
 
-  point_set.reserve (reader.m_nb_points);
-  
-  std::size_t points_read = 0;
-  
-  while (!(stream.eof()) && points_read < reader.m_nb_points)
+    if (element.name() == "vertex" || element.name() == "vertices")
     {
-      for (std::size_t i = 0; i < reader.readers().size (); ++ i)
-        reader.readers()[i]->get (stream);
-
-      filler.process_line (reader);
-      
-      ++ points_read;
+      point_set.reserve (element.number_of_items());
+      filler.instantiate_properties (element);
     }
-  // Skip remaining lines
+  
+    for (std::size_t j = 0; j < element.number_of_items(); ++ j)
+    {
+      for (std::size_t k = 0; k < element.number_of_properties(); ++ k)
+      {
+        internal::PLY::PLY_read_number* property = element.property(k);
+        property->get (stream);
 
-  return (points_read == reader.m_nb_points);
+        if (stream.eof())
+          return false;
+      }
+
+      if (element.name() == "vertex" || element.name() == "vertices")
+        filler.process_line (element);
+    }
+  }
+
+  return !stream.bad();
 }
   
 /*!
diff --git a/Point_set_processing_3/include/CGAL/IO/read_ply_points.h b/Point_set_processing_3/include/CGAL/IO/read_ply_points.h
index d96d61c76b1..11409af4abb 100644
--- a/Point_set_processing_3/include/CGAL/IO/read_ply_points.h
+++ b/Point_set_processing_3/include/CGAL/IO/read_ply_points.h
@@ -48,12 +48,12 @@
 
 #define TRY_TO_GENERATE_PROPERTY(STD_TYPE, T_TYPE, TYPE)          \
   if (type == STD_TYPE  || type == T_TYPE)                              \
-    m_properties->push_back (new PLY_read_typed_number< TYPE > (name, format))
+    m_elements.back().add_property (new PLY_read_typed_number< TYPE > (name, format))
 
 #define TRY_TO_GENERATE_SIZED_LIST_PROPERTY(STD_SIZE_TYPE, T_SIZE_TYPE, SIZE_TYPE, STD_INDEX_TYPE, T_INDEX_TYPE, INDEX_TYPE) \
   if ((size_type == STD_SIZE_TYPE  || size_type == T_SIZE_TYPE) &&      \
       (index_type == STD_INDEX_TYPE || index_type == T_INDEX_TYPE))     \
-    m_properties->push_back (new PLY_read_typed_list_with_typed_size< SIZE_TYPE , INDEX_TYPE > (name, format))
+    m_elements.back().add_property (new PLY_read_typed_list_with_typed_size< SIZE_TYPE , INDEX_TYPE > (name, format))
 
 #define TRY_TO_GENERATE_LIST_PROPERTY(STD_INDEX_TYPE, T_INDEX_TYPE, INDEX_TYPE) \
   TRY_TO_GENERATE_SIZED_LIST_PROPERTY("uchar", "uint8", boost::uint8_t, STD_INDEX_TYPE, T_INDEX_TYPE, INDEX_TYPE); \
@@ -317,24 +317,145 @@ namespace internal {
     }
   };
 
+  class PLY_element
+  {
+    std::string m_name;
+    std::size_t m_number;
+    
+    std::vector<PLY_read_number*> m_properties;
+  public:
+
+    PLY_element (const std::string& name, std::size_t number)
+      : m_name (name), m_number (number)
+    { }
+
+    PLY_element (const PLY_element& other)
+      : m_name (other.m_name), m_number (other.m_number), m_properties (other.m_properties)
+    {
+      const_cast<PLY_element&>(other).m_properties.clear();
+    }
+
+    PLY_element& operator= (const PLY_element& other)
+    {
+      m_name = other.m_name;
+      m_number = other.m_number;
+      m_properties = other.m_properties;
+      const_cast<PLY_element&>(other).m_properties.clear();
+      return *this;
+    }
+    
+    ~PLY_element()
+    {
+      for (std::size_t i = 0; i < m_properties.size(); ++ i)
+        delete m_properties[i];
+    }
+
+    const std::string& name() const { return m_name; }
+    std::size_t number_of_items() const { return m_number; }
+    std::size_t number_of_properties() const { return m_properties.size(); }
+
+    PLY_read_number* property (std::size_t idx) { return m_properties[idx]; }
+    
+    void add_property (PLY_read_number* read_number)
+    {
+      m_properties.push_back (read_number);
+    }
+
+    template <typename Type>
+    bool has_property (const char* tag)
+    {
+      return has_property (tag, Type());
+    }
+    template <typename Type>
+    bool has_property (const char* tag, const std::vector<Type>&)
+    {
+      for (std::size_t i = 0; i < number_of_properties(); ++ i)
+        if (m_properties[i]->name () == tag)
+          return (dynamic_cast<PLY_read_typed_list<Type>*>(m_properties[i]) != NULL);
+      return false;
+    }
+    
+    template <typename Type>
+    bool has_property (const char* tag, Type)
+    {
+      for (std::size_t i = 0; i < number_of_properties(); ++ i)
+        if (m_properties[i]->name () == tag)
+          return (dynamic_cast<PLY_read_typed_number<Type>*>(m_properties[i]) != NULL);
+      return false;
+    }
+    bool has_property (const char* tag, double)
+    {
+      for (std::size_t i = 0; i < number_of_properties(); ++ i)
+        if (m_properties[i]->name () == tag)
+          return (dynamic_cast<PLY_read_typed_number<double>*>(m_properties[i]) != NULL
+                  || dynamic_cast<PLY_read_typed_number<float>*>(m_properties[i]) != NULL);
+
+      return false;
+    }
+
+    template <typename Type>
+    void assign (Type& t, const char* tag)
+    {
+      for (std::size_t i = 0; i < number_of_properties (); ++ i)
+        if (m_properties[i]->name () == tag)
+        {
+          PLY_read_typed_number<Type>*
+            property = dynamic_cast<PLY_read_typed_number<Type>*>(m_properties[i]);
+          CGAL_assertion (property != NULL);
+          t = property->buffer();
+          return;
+        }
+    }
+
+    template <typename Type>
+    void assign (std::vector<Type>& t, const char* tag)
+    {
+      for (std::size_t i = 0; i < number_of_properties (); ++ i)
+        if (m_properties[i]->name () == tag)
+        {
+          PLY_read_typed_list<Type>*
+            property = dynamic_cast<PLY_read_typed_list<Type>*>(m_properties[i]);
+          CGAL_assertion (property != NULL);
+          t = property->buffer();
+          return;
+        }
+    }
+
+    void assign (double& t, const char* tag)
+    {
+      for (std::size_t i = 0; i < number_of_properties (); ++ i)
+        if (m_properties[i]->name () == tag)
+        {
+          PLY_read_typed_number<double>*
+            property_double = dynamic_cast<PLY_read_typed_number<double>*>(m_properties[i]);
+          if (property_double == NULL)
+          {
+            PLY_read_typed_number<float>*
+              property_float = dynamic_cast<PLY_read_typed_number<float>*>(m_properties[i]);
+            CGAL_assertion (property_float != NULL);
+            t = property_float->buffer();
+          }
+          else
+            t = property_double->buffer();
+          
+          return;
+        }
+    }
+  
+  };
 
   class PLY_reader
   {
-    std::vector<PLY_read_number*> m_point_properties;
-    std::vector<PLY_read_number*> m_face_properties;
-    std::vector<PLY_read_number*>* m_properties;
+    std::vector<PLY_element> m_elements;
     std::string m_comments;
 
   public:
-    std::size_t m_nb_points;
-    std::size_t m_nb_faces;
+    PLY_reader () { }
 
-    PLY_reader () : m_properties (&m_point_properties), m_nb_points (0), m_nb_faces(0)  { }
-
-    const std::vector<PLY_read_number*>& readers() const { return *m_properties; }
-    void read_faces()
+    std::size_t number_of_elements() const { return m_elements.size(); }
+    PLY_element& element (std::size_t idx)
     {
-      m_properties = &m_face_properties;
+      return m_elements[idx];
     }
 
     const std::string& comments() const { return m_comments; }
@@ -349,9 +470,6 @@ namespace internal {
       std::string line;
       std::istringstream iss;
 
-      // Check the order of the properties of the point set
-      bool reading_vertices = false, reading_faces = false;
-  
       while (getline (stream,line))
         {
           iss.clear();
@@ -401,9 +519,6 @@ namespace internal {
 
               if (keyword == "property")
                 {
-                  if (!reading_vertices && !reading_faces)
-                    continue;
-
                   std::string type, name;
                   if (!(iss >> type >> name))
                     {
@@ -455,9 +570,6 @@ namespace internal {
                       m_comments += "\n";
                     }
                 }
-              // When end_header is reached, stop loop and begin reading points
-              else if (keyword == "end_header")
-                break;
               else if (keyword == "element")
                 {
                   std::string type;
@@ -468,120 +580,18 @@ namespace internal {
                       return false;
                     }
                 
-                  if (type == "vertex")
-                    {
-                      m_nb_points = number;
-                      reading_vertices = true;
-                      reading_faces = false;
-                    }
-                
-                  else if (type == "face")
-                    {
-                      m_nb_faces = number;
-                      reading_faces = true;
-                      reading_vertices = false;
-                      m_properties = &m_face_properties;
-                    }
-                  else
-                    {
-                      reading_faces = false;
-                      reading_vertices = false;
-                      continue;
-                    }
+                  m_elements.push_back (PLY_element(type, number));
                 }
-            
+              // When end_header is reached, stop loop and begin reading points
+              else if (keyword == "end_header")
+                break;
             }
         }
-      m_properties = &m_point_properties;
       return true;
     }
 
     ~PLY_reader ()
     {
-      for (std::size_t i = 0; i < m_point_properties.size (); ++ i)
-        delete m_point_properties[i];
-      m_point_properties.clear();
-    }
-  
-    template <typename Type>
-    bool does_tag_exist (const char* tag)
-    {
-      return does_tag_exist (tag, Type());
-    }
-
-    template <typename Type>
-    void assign (Type& t, const char* tag)
-    {
-      for (std::size_t i = 0; i < m_properties->size (); ++ i)
-        if ((*m_properties)[i]->name () == tag)
-          {
-            PLY_read_typed_number<Type>*
-              reader = dynamic_cast<PLY_read_typed_number<Type>*>((*m_properties)[i]);
-            CGAL_assertion (reader != NULL);
-            t = reader->buffer();
-            return;
-          }
-    }
-
-    template <typename Type>
-    void assign (std::vector<Type>& t, const char* tag)
-    {
-      for (std::size_t i = 0; i < m_properties->size (); ++ i)
-        if ((*m_properties)[i]->name () == tag)
-          {
-            PLY_read_typed_list<Type>*
-              reader = dynamic_cast<PLY_read_typed_list<Type>*>((*m_properties)[i]);
-            CGAL_assertion (reader != NULL);
-            t = reader->buffer();
-            return;
-          }
-    }
-
-    template <typename Type>
-    bool does_tag_exist (const char* tag, const std::vector<Type>&)
-    {
-      for (std::size_t i = 0; i < m_properties->size (); ++ i)
-        if ((*m_properties)[i]->name () == tag)
-          return (dynamic_cast<PLY_read_typed_list<Type>*>((*m_properties)[i]) != NULL);
-      return false;
-    }
-    
-    template <typename Type>
-    bool does_tag_exist (const char* tag, Type)
-    {
-      for (std::size_t i = 0; i < m_properties->size (); ++ i)
-        if ((*m_properties)[i]->name () == tag)
-          return (dynamic_cast<PLY_read_typed_number<Type>*>((*m_properties)[i]) != NULL);
-      return false;
-    }
-    bool does_tag_exist (const char* tag, double)
-    {
-      for (std::size_t i = 0; i < m_properties->size (); ++ i)
-        if ((*m_properties)[i]->name () == tag)
-          return (dynamic_cast<PLY_read_typed_number<double>*>((*m_properties)[i]) != NULL
-                  || dynamic_cast<PLY_read_typed_number<float>*>((*m_properties)[i]) != NULL);
-
-      return false;
-    }
-    void assign (double& t, const char* tag)
-    {
-      for (std::size_t i = 0; i < m_properties->size (); ++ i)
-        if ((*m_properties)[i]->name () == tag)
-          {
-            PLY_read_typed_number<double>*
-              reader_double = dynamic_cast<PLY_read_typed_number<double>*>((*m_properties)[i]);
-            if (reader_double == NULL)
-              {
-                PLY_read_typed_number<float>*
-                  reader_float = dynamic_cast<PLY_read_typed_number<float>*>((*m_properties)[i]);
-                CGAL_assertion (reader_float != NULL);
-                t = reader_float->buffer();
-              }
-            else
-              t = reader_double->buffer();
-          
-            return;
-          }
     }
   
   };
@@ -640,12 +650,12 @@ namespace internal {
             typename PropertyMap,
             typename Constructor,
             typename ... T>
-  void process_properties (PLY_reader& reader, OutputValueType& new_element,
+  void process_properties (PLY_element& element, OutputValueType& new_element,
                            std::tuple<PropertyMap, Constructor, PLY_property<T>...>&& current)
   {
     typedef typename PropertyMap::value_type PmapValueType;
     std::tuple<T...> values;
-    Filler<sizeof...(T)-1>::fill(reader, values, current);
+    Filler<sizeof...(T)-1>::fill(element, values, current);
     PmapValueType new_value = call_functor<PmapValueType>(std::get<1>(current), values);
     put (std::get<0>(current), new_element, new_value);
   }
@@ -656,42 +666,42 @@ namespace internal {
             typename ... T,
             typename NextPropertyBinder,
             typename ... PropertyMapBinders>
-  void process_properties (PLY_reader& reader, OutputValueType& new_element,
+  void process_properties (PLY_element& element, OutputValueType& new_element,
                            std::tuple<PropertyMap, Constructor, PLY_property<T>...>&& current,
                            NextPropertyBinder&& next,
                            PropertyMapBinders&& ... properties)
   {
     typedef typename PropertyMap::value_type PmapValueType;
     std::tuple<T...> values;
-    Filler<sizeof...(T)-1>::fill(reader, values, current);
+    Filler<sizeof...(T)-1>::fill(element, values, current);
     PmapValueType new_value = call_functor<PmapValueType>(std::get<1>(current), values);
     put (std::get<0>(current), new_element, new_value);
   
-    process_properties (reader, new_element, std::forward<NextPropertyBinder>(next),
+    process_properties (element, new_element, std::forward<NextPropertyBinder>(next),
                         std::forward<PropertyMapBinders>(properties)...);
   }
 
 
   template <typename OutputValueType, typename PropertyMap, typename T>
-  void process_properties (PLY_reader& reader, OutputValueType& new_element,
+  void process_properties (PLY_element& element, OutputValueType& new_element,
                            std::pair<PropertyMap, PLY_property<T> >&& current)
   {
     T new_value = T();
-    reader.assign (new_value, current.second.name);
+    element.assign (new_value, current.second.name);
     put (current.first, new_element, new_value);
   }
 
   template <typename OutputValueType, typename PropertyMap, typename T,
             typename NextPropertyBinder, typename ... PropertyMapBinders>
-  void process_properties (PLY_reader& reader, OutputValueType& new_element,
+  void process_properties (PLY_element& element, OutputValueType& new_element,
                            std::pair<PropertyMap, PLY_property<T> >&& current,
                            NextPropertyBinder&& next,
                            PropertyMapBinders&& ... properties)
   {
     T new_value = T();
-    reader.assign (new_value, current.second.name);
+    element.assign (new_value, current.second.name);
     put (current.first, new_element, new_value);
-    process_properties (reader, new_element, std::forward<NextPropertyBinder>(next),
+    process_properties (element, new_element, std::forward<NextPropertyBinder>(next),
                         std::forward<PropertyMapBinders>(properties)...);
   }
 
@@ -757,24 +767,31 @@ bool read_ply_points_with_properties (std::istream& stream,
   if (!(reader.init (stream)))
     return false;
   
-  std::size_t points_read = 0;
-  
-  while (!(stream.eof()) && points_read < reader.m_nb_points)
+  for (std::size_t i = 0; i < reader.number_of_elements(); ++ i)
+  {
+    internal::PLY::PLY_element& element = reader.element(i);
+
+    for (std::size_t j = 0; j < element.number_of_items(); ++ j)
     {
-      for (std::size_t i = 0; i < reader.readers().size (); ++ i)
-        reader.readers()[i]->get (stream);
+      for (std::size_t k = 0; k < element.number_of_properties(); ++ k)
+      {
+        internal::PLY::PLY_read_number* property = element.property(k);
+        property->get (stream);
 
-      OutputValueType new_element;
+        if (stream.eof())
+          return false;
+      }
 
-      internal::PLY::process_properties (reader, new_element, std::forward<PropertyHandler>(properties)...);
-
-      *(output ++) = new_element;
-      
-      ++ points_read;
+      if (element.name() == "vertex" || element.name() == "vertices")
+      {
+        OutputValueType new_element;
+        internal::PLY::process_properties (element, new_element, std::forward<PropertyHandler>(properties)...);
+        *(output ++) = new_element;
+      }
     }
-  // Skip remaining lines
+  }
 
-  return (points_read == reader.m_nb_points);
+  return true;
 }
 
 /// \cond SKIP_IN_MANUAL
diff --git a/Polyhedron_IO/include/CGAL/IO/PLY_reader.h b/Polyhedron_IO/include/CGAL/IO/PLY_reader.h
index 371114bd353..36b8a16fef1 100644
--- a/Polyhedron_IO/include/CGAL/IO/PLY_reader.h
+++ b/Polyhedron_IO/include/CGAL/IO/PLY_reader.h
@@ -29,36 +29,40 @@ namespace CGAL{
     template <typename Integer, class Polygon_3, class Color_rgb>
     bool
     read_PLY_faces (std::istream& in,
-                    PLY::PLY_reader& reader,
+                    internal::PLY::PLY_element& element,
                     std::vector< Polygon_3 >& polygons,
                     std::vector< Color_rgb >& fcolors,
                     const char* vertex_indices_tag)
     {
-      std::size_t faces_read = 0;
-
       bool has_colors = false;
       std::string rtag = "r", gtag = "g", btag = "b";
-      if ((reader.does_tag_exist<boost::uint8_t>("red") || reader.does_tag_exist<boost::uint8_t>("r")) &&
-          (reader.does_tag_exist<boost::uint8_t>("green") || reader.does_tag_exist<boost::uint8_t>("g")) &&
-          (reader.does_tag_exist<boost::uint8_t>("blue") || reader.does_tag_exist<boost::uint8_t>("b")))
+      if ((element.has_property<boost::uint8_t>("red") || element.has_property<boost::uint8_t>("r")) &&
+          (element.has_property<boost::uint8_t>("green") || element.has_property<boost::uint8_t>("g")) &&
+          (element.has_property<boost::uint8_t>("blue") || element.has_property<boost::uint8_t>("b")))
       {
         has_colors = true;
-        if (reader.does_tag_exist<boost::uint8_t>("red"))
+        if (element.has_property<boost::uint8_t>("red"))
         {
           rtag = "red"; gtag = "green"; btag = "blue";
         }
       }
       
-      while (!(in.eof()) && faces_read < reader.m_nb_faces)
+      for (std::size_t j = 0; j < element.number_of_items(); ++ j)
       {
-        for (std::size_t i = 0; i < reader.readers().size (); ++ i)
-          reader.readers()[i]->get (in);
+        for (std::size_t k = 0; k < element.number_of_properties(); ++ k)
+        {
+          internal::PLY::PLY_read_number* property = element.property(k);
+          property->get (in);
+
+          if (in.eof())
+            return false;
+        }
 
         cpp11::tuple<std::vector<Integer>, boost::uint8_t, boost::uint8_t, boost::uint8_t> new_face; 
 
         if (has_colors)
         {
-          PLY::process_properties (reader, new_face,
+          PLY::process_properties (element, new_face,
                                    std::make_pair (CGAL::make_nth_of_tuple_property_map<0>(new_face),
                                                    PLY_property<std::vector<Integer> >(vertex_indices_tag)),
                                    std::make_pair (CGAL::make_nth_of_tuple_property_map<1>(new_face),
@@ -71,15 +75,13 @@ namespace CGAL{
           fcolors.push_back (Color_rgb (get<1>(new_face), get<2>(new_face), get<3>(new_face)));
         }
         else
-          PLY::process_properties (reader, new_face,
+          PLY::process_properties (element, new_face,
                                    std::make_pair (CGAL::make_nth_of_tuple_property_map<0>(new_face),
                                                    PLY_property<std::vector<Integer> >(vertex_indices_tag)));
 
         polygons.push_back (Polygon_3(get<0>(new_face).size()));
         for (std::size_t i = 0; i < get<0>(new_face).size(); ++ i)
           polygons.back()[i] = std::size_t(get<0>(new_face)[i]);
-
-        ++ faces_read;
       }
 
       return !in.bad();
@@ -106,45 +108,66 @@ namespace CGAL{
     if (!(reader.init (in)))
       return false;
   
-    std::size_t points_read = 0;
-  
-    while (!(in.eof()) && points_read < reader.m_nb_points)
+    for (std::size_t i = 0; i < reader.number_of_elements(); ++ i)
     {
-      for (std::size_t i = 0; i < reader.readers().size (); ++ i)
-        reader.readers()[i]->get (in);
+      internal::PLY::PLY_element& element = reader.element(i);
 
-      Point_3 new_vertex;
+      if (element.name() == "vertex" || element.name() == "vertices")
+      {
+        for (std::size_t j = 0; j < element.number_of_items(); ++ j)
+        {
+          for (std::size_t k = 0; k < element.number_of_properties(); ++ k)
+          {
+            internal::PLY::PLY_read_number* property = element.property(k);
+            property->get (in);
 
-      internal::PLY::process_properties (reader, new_vertex,
-                                         make_ply_point_reader (CGAL::Identity_property_map<Point_3>()));
-                                         
+            if (in.eof())
+              return false;
+          }
 
-      points.push_back (new_vertex);
+          Point_3 new_vertex;
+
+          internal::PLY::process_properties (element, new_vertex,
+                                             make_ply_point_reader (CGAL::Identity_property_map<Point_3>()));
       
-      ++ points_read;
+          points.push_back (get<0>(new_vertex));
+        }
+      }
+      else if (element.name() == "face" || element.name() == "faces")
+      {
+        std::vector<CGAL::Color> dummy;
+
+        if (element.has_property<std::vector<boost::int32_t> > ("vertex_indices"))
+          internal::read_PLY_faces<boost::int32_t> (in, element, polygons, dummy, "vertex_indices");
+        else if (element.has_property<std::vector<boost::uint32_t> > ("vertex_indices"))
+          internal::read_PLY_faces<boost::uint32_t> (in, element, polygons, dummy, "vertex_indices");
+        else if (element.has_property<std::vector<boost::int32_t> > ("vertex_index"))
+          internal::read_PLY_faces<boost::int32_t> (in, element, polygons, dummy, "vertex_index");
+        else if (element.has_property<std::vector<boost::uint32_t> > ("vertex_index"))
+          internal::read_PLY_faces<boost::uint32_t> (in, element, polygons, dummy, "vertex_index");
+        else
+        {
+          std::cerr << "Error: can't find vertex indices in PLY input" << std::endl;
+          return false;
+        }
+      }
+      else // Read other elements and ignore
+      {
+        for (std::size_t j = 0; j < element.number_of_items(); ++ j)
+        {
+          for (std::size_t k = 0; k < element.number_of_properties(); ++ k)
+          {
+            internal::PLY::PLY_read_number* property = element.property(k);
+            property->get (in);
+
+            if (in.eof())
+              return false;
+          }
+        }
+      }
     }
 
-    if (points_read != reader.m_nb_points)
-      return false;
-
-    std::vector<CGAL::Color> dummy;
-
-    reader.read_faces();
-
-    if (reader.does_tag_exist<std::vector<boost::int32_t> > ("vertex_indices"))
-      return internal::read_PLY_faces<boost::int32_t> (in, reader, polygons, dummy, "vertex_indices");
-
-    if (reader.does_tag_exist<std::vector<boost::uint32_t> > ("vertex_indices"))
-      return internal::read_PLY_faces<boost::uint32_t> (in, reader, polygons, dummy, "vertex_indices");
-
-    if (reader.does_tag_exist<std::vector<boost::int32_t> > ("vertex_index"))
-      return internal::read_PLY_faces<boost::int32_t> (in, reader, polygons, dummy, "vertex_index");
-
-    if (reader.does_tag_exist<std::vector<boost::uint32_t> > ("vertex_index"))
-      return internal::read_PLY_faces<boost::uint32_t> (in, reader, polygons, dummy, "vertex_index");
-
-    std::cerr << "Error: can't find vertex indices in PLY input" << std::endl;
-    return false;
+    return !in.bad();
   }
 
   template <class Point_3, class Polygon_3, class Color_rgb>
@@ -161,75 +184,96 @@ namespace CGAL{
       std::cerr << "Error: cannot open file" << std::endl;
       return false;
     }
-
     internal::PLY::PLY_reader reader;
   
     if (!(reader.init (in)))
       return false;
-  
-    std::size_t points_read = 0;
 
-    bool has_colors = false;
-    std::string rtag = "r", gtag = "g", btag = "b";
-    if ((reader.does_tag_exist<boost::uint8_t>("red") || reader.does_tag_exist<boost::uint8_t>("r")) &&
-        (reader.does_tag_exist<boost::uint8_t>("green") || reader.does_tag_exist<boost::uint8_t>("g")) &&
-        (reader.does_tag_exist<boost::uint8_t>("blue") || reader.does_tag_exist<boost::uint8_t>("b")))
+    for (std::size_t i = 0; i < reader.number_of_elements(); ++ i)
     {
-      has_colors = true;
-      if (reader.does_tag_exist<boost::uint8_t>("red"))
+      internal::PLY::PLY_element& element = reader.element(i);
+
+      if (element.name() == "vertex" || element.name() == "vertices")
       {
-        rtag = "red"; gtag = "green"; btag = "blue";
+        bool has_colors = false;
+        std::string rtag = "r", gtag = "g", btag = "b";
+        if ((element.has_property<boost::uint8_t>("red") || element.has_property<boost::uint8_t>("r")) &&
+            (element.has_property<boost::uint8_t>("green") || element.has_property<boost::uint8_t>("g")) &&
+            (element.has_property<boost::uint8_t>("blue") || element.has_property<boost::uint8_t>("b")))
+        {
+          has_colors = true;
+          if (element.has_property<boost::uint8_t>("red"))
+          {
+            rtag = "red"; gtag = "green"; btag = "blue";
+          }
+        }
+
+        for (std::size_t j = 0; j < element.number_of_items(); ++ j)
+        {
+          for (std::size_t k = 0; k < element.number_of_properties(); ++ k)
+          {
+            internal::PLY::PLY_read_number* property = element.property(k);
+            property->get (in);
+
+            if (in.eof())
+              return false;
+          }
+
+          cpp11::tuple<Point_3, boost::uint8_t, boost::uint8_t, boost::uint8_t> new_vertex;
+
+          if (has_colors)
+          {
+            internal::PLY::process_properties (element, new_vertex,
+                                               make_ply_point_reader (CGAL::make_nth_of_tuple_property_map<0>(new_vertex)),
+                                               std::make_pair (CGAL::make_nth_of_tuple_property_map<1>(new_vertex),
+                                                               PLY_property<boost::uint8_t>(rtag.c_str())),
+                                               std::make_pair (CGAL::make_nth_of_tuple_property_map<2>(new_vertex),
+                                                               PLY_property<boost::uint8_t>(gtag.c_str())),
+                                               std::make_pair (CGAL::make_nth_of_tuple_property_map<3>(new_vertex),
+                                                               PLY_property<boost::uint8_t>(btag.c_str())));
+
+            vcolors.push_back (Color_rgb (get<1>(new_vertex), get<2>(new_vertex), get<3>(new_vertex)));
+          }
+          else
+            internal::PLY::process_properties (element, new_vertex,
+                                               make_ply_point_reader (CGAL::make_nth_of_tuple_property_map<0>(new_vertex)));
+
+          points.push_back (get<0>(new_vertex));
+        }
+      }
+      else if (element.name() == "face" || element.name() == "faces")
+      {
+        if (element.has_property<std::vector<boost::int32_t> > ("vertex_indices"))
+          internal::read_PLY_faces<boost::int32_t> (in, element, polygons, fcolors, "vertex_indices");
+        else if (element.has_property<std::vector<boost::uint32_t> > ("vertex_indices"))
+          internal::read_PLY_faces<boost::uint32_t> (in, element, polygons, fcolors, "vertex_indices");
+        else if (element.has_property<std::vector<boost::int32_t> > ("vertex_index"))
+          internal::read_PLY_faces<boost::int32_t> (in, element, polygons, fcolors, "vertex_index");
+        else if (element.has_property<std::vector<boost::uint32_t> > ("vertex_index"))
+          internal::read_PLY_faces<boost::uint32_t> (in, element, polygons, fcolors, "vertex_index");
+        else
+        {
+          std::cerr << "Error: can't find vertex indices in PLY input" << std::endl;
+          return false;
+        }
+      }
+      else // Read other elements and ignore
+      {
+        for (std::size_t j = 0; j < element.number_of_items(); ++ j)
+        {
+          for (std::size_t k = 0; k < element.number_of_properties(); ++ k)
+          {
+            internal::PLY::PLY_read_number* property = element.property(k);
+            property->get (in);
+
+            if (in.eof())
+              return false;
+          }
+        }
       }
     }
 
-    while (!(in.eof()) && points_read < reader.m_nb_points)
-    {
-      for (std::size_t i = 0; i < reader.readers().size (); ++ i)
-        reader.readers()[i]->get (in);
-
-      cpp11::tuple<Point_3, boost::uint8_t, boost::uint8_t, boost::uint8_t> new_vertex;
-
-      if (has_colors)
-      {
-        internal::PLY::process_properties (reader, new_vertex,
-                                           make_ply_point_reader (CGAL::make_nth_of_tuple_property_map<0>(new_vertex)),
-                                           std::make_pair (CGAL::make_nth_of_tuple_property_map<1>(new_vertex),
-                                                           PLY_property<boost::uint8_t>(rtag.c_str())),
-                                           std::make_pair (CGAL::make_nth_of_tuple_property_map<2>(new_vertex),
-                                                           PLY_property<boost::uint8_t>(gtag.c_str())),
-                                           std::make_pair (CGAL::make_nth_of_tuple_property_map<3>(new_vertex),
-                                                           PLY_property<boost::uint8_t>(btag.c_str())));
-
-        vcolors.push_back (Color_rgb (get<1>(new_vertex), get<2>(new_vertex), get<3>(new_vertex)));
-      }
-      else
-        internal::PLY::process_properties (reader, new_vertex,
-                                           make_ply_point_reader (CGAL::make_nth_of_tuple_property_map<0>(new_vertex)));
-      
-      points.push_back (get<0>(new_vertex));
-      
-      ++ points_read;
-    }
-
-    if (points_read != reader.m_nb_points)
-      return false;
-
-    reader.read_faces();
-    
-    if (reader.does_tag_exist<std::vector<boost::int32_t> > ("vertex_indices"))
-      return internal::read_PLY_faces<boost::int32_t> (in, reader, polygons, fcolors, "vertex_indices");
-
-    if (reader.does_tag_exist<std::vector<boost::uint32_t> > ("vertex_indices"))
-      return internal::read_PLY_faces<boost::uint32_t> (in, reader, polygons, fcolors, "vertex_indices");
-
-    if (reader.does_tag_exist<std::vector<boost::int32_t> > ("vertex_index"))
-      return internal::read_PLY_faces<boost::int32_t> (in, reader, polygons, fcolors, "vertex_index");
-
-    if (reader.does_tag_exist<std::vector<boost::uint32_t> > ("vertex_index"))
-      return internal::read_PLY_faces<boost::uint32_t> (in, reader, polygons, fcolors, "vertex_index");
-
-    std::cerr << "Error: can't find vertex indices in PLY input" << std::endl;
-    return false;
+    return !in.bad();
   }
 
 

From 1cfcb1ba5b57d4989866d4a6749ec311f05186d2 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Tue, 18 Dec 2018 15:28:02 +0100
Subject: [PATCH 15/35] Fix stream status handling everywhere PLY reader is
 used

---
 Point_set_3/include/CGAL/Point_set_3/IO.h     |  8 ++++---
 .../include/CGAL/IO/read_ply_points.h         |  5 ++++-
 Polyhedron_IO/include/CGAL/IO/PLY_reader.h    | 22 ++++++++++++-------
 3 files changed, 23 insertions(+), 12 deletions(-)

diff --git a/Point_set_3/include/CGAL/Point_set_3/IO.h b/Point_set_3/include/CGAL/Point_set_3/IO.h
index 287cf45e728..6d6e3acb2aa 100644
--- a/Point_set_3/include/CGAL/Point_set_3/IO.h
+++ b/Point_set_3/include/CGAL/Point_set_3/IO.h
@@ -394,7 +394,10 @@ read_ply_point_set(
   internal::PLY::Point_set_3_filler<Point, Vector> filler(point_set);
   
   if (!(reader.init (stream)))
+  {
+    stream.setstate(std::ios::failbit);
     return false;
+  }
 
   if (comments != NULL)
     *comments = reader.comments();
@@ -415,8 +418,7 @@ read_ply_point_set(
       {
         internal::PLY::PLY_read_number* property = element.property(k);
         property->get (stream);
-
-        if (stream.eof())
+        if (stream.fail())
           return false;
       }
 
@@ -425,7 +427,7 @@ read_ply_point_set(
     }
   }
 
-  return !stream.bad();
+  return true;
 }
   
 /*!
diff --git a/Point_set_processing_3/include/CGAL/IO/read_ply_points.h b/Point_set_processing_3/include/CGAL/IO/read_ply_points.h
index 11409af4abb..a7cc94683ca 100644
--- a/Point_set_processing_3/include/CGAL/IO/read_ply_points.h
+++ b/Point_set_processing_3/include/CGAL/IO/read_ply_points.h
@@ -765,7 +765,10 @@ bool read_ply_points_with_properties (std::istream& stream,
   internal::PLY::PLY_reader reader;
   
   if (!(reader.init (stream)))
+  {
+    stream.setstate(std::ios::failbit);
     return false;
+  }
   
   for (std::size_t i = 0; i < reader.number_of_elements(); ++ i)
   {
@@ -778,7 +781,7 @@ bool read_ply_points_with_properties (std::istream& stream,
         internal::PLY::PLY_read_number* property = element.property(k);
         property->get (stream);
 
-        if (stream.eof())
+        if (stream.fail())
           return false;
       }
 
diff --git a/Polyhedron_IO/include/CGAL/IO/PLY_reader.h b/Polyhedron_IO/include/CGAL/IO/PLY_reader.h
index 36b8a16fef1..2f7522c36fd 100644
--- a/Polyhedron_IO/include/CGAL/IO/PLY_reader.h
+++ b/Polyhedron_IO/include/CGAL/IO/PLY_reader.h
@@ -54,7 +54,7 @@ namespace CGAL{
           internal::PLY::PLY_read_number* property = element.property(k);
           property->get (in);
 
-          if (in.eof())
+          if (in.fail())
             return false;
         }
 
@@ -84,7 +84,7 @@ namespace CGAL{
           polygons.back()[i] = std::size_t(get<0>(new_face)[i]);
       }
 
-      return !in.bad();
+      return true;
     }
 
   }
@@ -106,7 +106,10 @@ namespace CGAL{
     internal::PLY::PLY_reader reader;
   
     if (!(reader.init (in)))
+    {
+      in.setstate(std::ios::failbit);
       return false;
+    }
   
     for (std::size_t i = 0; i < reader.number_of_elements(); ++ i)
     {
@@ -121,7 +124,7 @@ namespace CGAL{
             internal::PLY::PLY_read_number* property = element.property(k);
             property->get (in);
 
-            if (in.eof())
+            if (in.fail())
               return false;
           }
 
@@ -160,14 +163,14 @@ namespace CGAL{
             internal::PLY::PLY_read_number* property = element.property(k);
             property->get (in);
 
-            if (in.eof())
+            if (in.fail())
               return false;
           }
         }
       }
     }
 
-    return !in.bad();
+    return true;
   }
 
   template <class Point_3, class Polygon_3, class Color_rgb>
@@ -187,7 +190,10 @@ namespace CGAL{
     internal::PLY::PLY_reader reader;
   
     if (!(reader.init (in)))
+    {
+      in.setstate(std::ios::failbit);
       return false;
+    }
 
     for (std::size_t i = 0; i < reader.number_of_elements(); ++ i)
     {
@@ -215,7 +221,7 @@ namespace CGAL{
             internal::PLY::PLY_read_number* property = element.property(k);
             property->get (in);
 
-            if (in.eof())
+            if (in.fail())
               return false;
           }
 
@@ -266,14 +272,14 @@ namespace CGAL{
             internal::PLY::PLY_read_number* property = element.property(k);
             property->get (in);
 
-            if (in.eof())
+            if (in.fail())
               return false;
           }
         }
       }
     }
 
-    return !in.bad();
+    return true;
   }
 
 

From 0f0f03e08bb048eb47efeee0201d41d5a4648188 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= <sebastien.loriot@cgal.org>
Date: Wed, 19 Dec 2018 16:05:08 +0100
Subject: [PATCH 16/35] handle empty meshes

---
 .../include/CGAL/Polygon_mesh_processing/bbox.h          | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/bbox.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/bbox.h
index 033364bac9b..a82985d40a4 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/bbox.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/bbox.h
@@ -76,13 +76,12 @@ namespace CGAL {
       GT gt = choose_param(get_param(np, internal_np::geom_traits), GT());
       typename GT::Construct_bbox_3 get_bbox = gt.construct_bbox_3_object();
 
-      typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
+      typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
 
-      halfedge_descriptor h0 = *(halfedges(pmesh).first);
-      CGAL::Bbox_3 bb = get(vpm, target(h0, pmesh)).bbox();
-      BOOST_FOREACH(halfedge_descriptor h, halfedges(pmesh))
+      CGAL::Bbox_3 bb;
+      BOOST_FOREACH(vertex_descriptor v, vertices(pmesh))
       {
-        bb += get_bbox( get(vpm, target(h, pmesh)) );
+        bb += get_bbox( get(vpm, v) );
       }
       return bb;
     }

From 725e773bdd0e0a43e113f6078db05151844feff7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= <sebastien.loriot@cgal.org>
Date: Thu, 20 Dec 2018 09:31:52 +0100
Subject: [PATCH 17/35] add nef to polygon soup

---
 .../convert_nef_polyhedron_to_polygon_mesh.h  | 41 +++++++++++--------
 1 file changed, 25 insertions(+), 16 deletions(-)

diff --git a/Nef_3/include/CGAL/boost/graph/convert_nef_polyhedron_to_polygon_mesh.h b/Nef_3/include/CGAL/boost/graph/convert_nef_polyhedron_to_polygon_mesh.h
index 369db23902d..7393881de00 100644
--- a/Nef_3/include/CGAL/boost/graph/convert_nef_polyhedron_to_polygon_mesh.h
+++ b/Nef_3/include/CGAL/boost/graph/convert_nef_polyhedron_to_polygon_mesh.h
@@ -343,36 +343,45 @@ void collect_polygon_mesh_info(
 
 } //end of namespace nef_to_pm
 
-
-template <class Nef_polyhedron, class Polygon_mesh>
-void convert_nef_polyhedron_to_polygon_mesh(const Nef_polyhedron& nef, Polygon_mesh& pm, bool triangulate_all_faces = false)
+template <class Output_kernel, class Nef_polyhedron>
+void convert_nef_polyhedron_to_polygon_soup(const Nef_polyhedron& nef,
+                                                  std::vector<typename Output_kernel::Point_3>& points,
+                                                  std::vector< std::vector<std::size_t> >& polygons,
+                                                  bool triangulate_all_faces = false)
 {
   typedef typename Nef_polyhedron::Point_3 Point_3;
-  typedef typename boost::property_traits<typename boost::property_map<Polygon_mesh, vertex_point_t>::type>::value_type PM_Point;
-  typedef typename Kernel_traits<PM_Point>::Kernel PM_Kernel;
   typedef typename Kernel_traits<Point_3>::Kernel Nef_Kernel;
-  typedef Cartesian_converter<Nef_Kernel, PM_Kernel> Converter;
-
+  typedef Cartesian_converter<Nef_Kernel, Output_kernel> Converter;
+  typedef typename Output_kernel::Point_3 Out_point;
   typename Nef_polyhedron::Volume_const_iterator vol_it = nef.volumes_begin(),
                                                  vol_end = nef.volumes_end();
   if ( Nef_polyhedron::Infi_box::extended_kernel() ) ++vol_it; // skip Infi_box
   CGAL_assertion ( vol_it != vol_end );
   ++vol_it; // skip unbounded volume
 
+  Converter to_output;
+  for (;vol_it!=vol_end;++vol_it)
+    nef_to_pm::collect_polygon_mesh_info<Out_point>(points,
+                                                    polygons,
+                                                    nef,
+                                                    vol_it->shells_begin(),
+                                                    to_output,
+                                                    triangulate_all_faces);
+}
+
+template <class Nef_polyhedron, class Polygon_mesh>
+void convert_nef_polyhedron_to_polygon_mesh(const Nef_polyhedron& nef, Polygon_mesh& pm, bool triangulate_all_faces = false)
+{
+  typedef typename boost::property_traits<typename boost::property_map<Polygon_mesh, vertex_point_t>::type>::value_type PM_Point;
+  typedef typename Kernel_traits<PM_Point>::Kernel PM_Kernel;
+
   std::vector<PM_Point> points;
   std::vector< std::vector<std::size_t> > polygons;
-  Converter to_inexact;
-  for (;vol_it!=vol_end;++vol_it)
-    nef_to_pm::collect_polygon_mesh_info<PM_Point>(points,
-                                                  polygons,
-                                                  nef,
-                                                  vol_it->shells_begin(),
-                                                  to_inexact,
-                                                  triangulate_all_faces);
-
+  convert_nef_polyhedron_to_polygon_soup<PM_Kernel>(nef, points, polygons, triangulate_all_faces);
   Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, polygons, pm);
 }
 
+
 } //end of namespace CGAL
 
 

From 2c6e9c2b70b5975827ca1b8cdf6f5a8091aa5992 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Thu, 20 Dec 2018 15:57:51 +0100
Subject: [PATCH 18/35] Fixed memory leaks in partition code

---
 .../CGAL/boost/graph/METIS/partition_dual_graph.h | 15 ++++++++++++---
 .../CGAL/boost/graph/METIS/partition_graph.h      | 15 ++++++++++++---
 2 files changed, 24 insertions(+), 6 deletions(-)

diff --git a/BGL/include/CGAL/boost/graph/METIS/partition_dual_graph.h b/BGL/include/CGAL/boost/graph/METIS/partition_dual_graph.h
index 8bfeb002ae1..7e2a15870a7 100644
--- a/BGL/include/CGAL/boost/graph/METIS/partition_dual_graph.h
+++ b/BGL/include/CGAL/boost/graph/METIS/partition_dual_graph.h
@@ -36,12 +36,15 @@
 #include <boost/graph/graph_traits.hpp>
 #include <boost/tuple/tuple.hpp>
 
+#include <cstdlib>
+
 namespace CGAL {
 
 namespace METIS {
 
 template<typename TriangleMesh, typename METIS_options, typename NamedParameters>
-void partition_dual_graph(const TriangleMesh& tm, int nparts,
+void partition_dual_graph(const TriangleMesh& tm,
+                          int nparts,
                           METIS_options options, // options array
                           const NamedParameters& np)
 {
@@ -93,11 +96,11 @@ void partition_dual_graph(const TriangleMesh& tm, int nparts,
   idx_t objval;
 
   // partition info for the nodes
-  idx_t* npart = (idx_t*) calloc(nn, sizeof(idx_t));
+  idx_t* npart = (idx_t*) calloc(num_vertices(tm), sizeof(idx_t));
   CGAL_assertion(npart != NULL);
 
   // partition info for the elements
-  idx_t* epart = (idx_t*) calloc(ne, sizeof(idx_t));
+  idx_t* epart = (idx_t*) calloc(num_faces(tm), sizeof(idx_t));
   CGAL_assertion(epart != NULL);
 
   // do not support Fortran-style arrays
@@ -118,6 +121,12 @@ void partition_dual_graph(const TriangleMesh& tm, int nparts,
   Output_face_partition_ids fo;
   vo(tm, indices, npart, get_param(np, internal_np::vertex_partition_id));
   fo(tm, epart, get_param(np, internal_np::face_partition_id));
+
+  delete[] eptr;
+  delete[] eind;
+
+  std::free(npart);
+  std::free(epart);
 }
 
 template<typename TriangleMesh, typename NamedParameters>
diff --git a/BGL/include/CGAL/boost/graph/METIS/partition_graph.h b/BGL/include/CGAL/boost/graph/METIS/partition_graph.h
index 65176d37131..d0f3e56b759 100644
--- a/BGL/include/CGAL/boost/graph/METIS/partition_graph.h
+++ b/BGL/include/CGAL/boost/graph/METIS/partition_graph.h
@@ -34,6 +34,8 @@
 #include <boost/graph/graph_traits.hpp>
 #include <boost/tuple/tuple.hpp>
 
+#include <cstdlib>
+
 namespace CGAL {
 
 namespace METIS {
@@ -76,7 +78,8 @@ struct Output_face_partition_ids
 };
 
 template<typename TriangleMesh, typename METIS_options, typename NamedParameters>
-void partition_graph(const TriangleMesh& tm, int nparts,
+void partition_graph(const TriangleMesh& tm,
+                     int nparts,
                      METIS_options options, // pointer to the options array
                      const NamedParameters& np)
 {
@@ -125,11 +128,11 @@ void partition_graph(const TriangleMesh& tm, int nparts,
   idx_t objval;
 
   // partition info for the nodes
-  idx_t* npart = (idx_t*) calloc(nn, sizeof(idx_t));
+  idx_t* npart = (idx_t*) calloc(num_vertices(tm), sizeof(idx_t));
   CGAL_assertion(npart != NULL);
 
   // partition info for the elements
-  idx_t* epart = (idx_t*) calloc(ne, sizeof(idx_t));
+  idx_t* epart = (idx_t*) calloc(num_faces(tm), sizeof(idx_t));
   CGAL_assertion(epart != NULL);
 
   // do not support Fortran-style arrays
@@ -150,6 +153,12 @@ void partition_graph(const TriangleMesh& tm, int nparts,
   Output_face_partition_ids fo;
   vo(tm, indices, npart, get_param(np, internal_np::vertex_partition_id));
   fo(tm, epart, get_param(np, internal_np::face_partition_id));
+
+  delete[] eptr;
+  delete[] eind;
+
+  std::free(npart);
+  std::free(epart);
 }
 
 template<typename TriangleMesh, typename NamedParameters>

From fe5c3a77f97b8f3ee1a644ebb2aa74ac2d7a26b3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= <sebastien.loriot@cgal.org>
Date: Fri, 21 Dec 2018 09:51:58 +0100
Subject: [PATCH 19/35] swap the edge too

---
 .../internal/Isotropic_remeshing/remesh_impl.h                   | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h
index 08d90b8f234..0e0596fc1fb 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h
@@ -686,6 +686,7 @@ namespace internal {
           if (is_on_border(he) || is_on_mesh(he))
           {
             he = opposite(he, mesh_); //he now is PATCH_BORDER
+            e = edge(he, mesh_);
             CGAL_assertion(is_on_patch_border(he));
           }
         }//end if(not on PATCH)

From 04b0a2c48a9cf19cb9b7dc6bc0e82443c88468b0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= <sebastien.loriot@cgal.org>
Date: Mon, 31 Dec 2018 08:15:28 +0100
Subject: [PATCH 20/35] doc clean-up

improve doxygen likeliness
---
 QP_solver/doc/QP_solver/CGAL/QP_functions.h | 115 +++++---------------
 QP_solver/doc/QP_solver/CGAL/QP_models.h    |  30 ++---
 QP_solver/doc/QP_solver/CGAL/QP_options.h   |   3 +-
 QP_solver/doc/QP_solver/CGAL/QP_solution.h  |   6 +-
 4 files changed, 39 insertions(+), 115 deletions(-)

diff --git a/QP_solver/doc/QP_solver/CGAL/QP_functions.h b/QP_solver/doc/QP_solver/CGAL/QP_functions.h
index 1c8b614dd2b..f06b7f6f9a4 100644
--- a/QP_solver/doc/QP_solver/CGAL/QP_functions.h
+++ b/QP_solver/doc/QP_solver/CGAL/QP_functions.h
@@ -4,6 +4,8 @@ namespace CGAL {
 /// @{
 
 /*!
+\tparam LinearProgram a model of `LinearProgram`.
+
 This function writes a linear program to an output stream (in 
 `MPSFormat`). The time 
 complexity is \f$ \Theta (mn)\f$, even if the program is very sparse. 
@@ -12,18 +14,9 @@ It writes the linear program `lp` to `out` in
 `MPSFormat`. The name of the program will be the one provided 
 by `problem_name`.
 
-Requirements 
--------------- 
-
+\cgalHeading{Requirements}
 Output operators are defined for all entry types of `lp`. 
 
-Example 
--------------- 
-
-\ref QP_solver/print_first_lp.cpp 
-
-\sa The concept `LinearProgram` 
-
 */
 template <LinearProgram>
 void print_linear_program 
@@ -32,6 +25,8 @@ const std::string& problem_name = std::string("MY_MPS"));
 
 
 /*!
+\tparam NonnegativeLinearProgram a model of `NonnegativeLinearProgram`
+
 This function writes a nonnegative linear program to an output stream 
 (in `MPSFormat`). The time 
 complexity is \f$ \Theta (mn)\f$, even if the program is very sparse. 
@@ -40,18 +35,9 @@ Writes the nonnegative linear program `lp` to `out` in
 `MPSFormat`. The name of the program will be the one provided 
 by `problem_name`.
 
-Requirements 
--------------- 
-
+\cgalHeading{Requirements}
 Output operators are defined for all entry types of `lp`. 
 
-Example 
--------------- 
-
-\ref QP_solver/print_first_nonnegative_lp.cpp 
-
-\sa The concept `NonnegativeLinearProgram` 
-
 */
 template <NonnegativeLinearProgram>
 void print_nonnegative_linear_program 
@@ -60,6 +46,8 @@ const std::string& problem_name = std::string("MY_MPS"));
 
 
 /*!
+\tparam NonnegativeQuadraticProgram a model of `NonnegativeQuadraticProgram`
+
 This function writes a nonnegative quadratic program 
 to an output stream (in `MPSFormat`). The time 
 complexity is \f$ \Theta (n^2 + mn)\f$, even if the program is very sparse. 
@@ -68,18 +56,9 @@ Writes the nonnegative quadratic program `qp` to `out` in
 `MPSFormat`. The name of the program will be the one provided 
 by `problem_name`.
 
-Requirements 
--------------- 
-
+\cgalHeading{Requirements}
 Output operators are defined for all entry types of `qp`. 
 
-Example 
--------------- 
-
-\ref QP_solver/print_first_nonnegative_qp.cpp 
-
-\sa The concept `NonnegativeQuadraticProgram` 
-
 */
 template <NonnegativeQuadraticProgram>
 void print_nonnegative_quadratic_program 
@@ -88,6 +67,9 @@ const std::string& problem_name = std::string("MY_MPS"));
 
 
 /*!
+\tparam QuadraticProgram a model of `QuadraticProgram`
+
+
 This function writes a quadratic program to an output stream (in 
 `MPSFormat`). The time complexity is \f$ \Theta (n^2 + mn)\f$, even 
 if the program is very sparse. 
@@ -95,17 +77,9 @@ if the program is very sparse.
 Writes the quadratic program `qp` to `out` in `MPSFormat`.
 The name of the program will be the one provided by `problem_name`.
 
-Requirements 
--------------- 
-
+\cgalHeading{Requirements}
 Output operators are defined for all entry types of `qp`. 
 
-Example 
--------------- 
-
-\ref QP_solver/print_first_qp.cpp 
-
-\sa The concept `QuadraticProgram` 
 */
 template <QuadraticProgram>
 void print_quadratic_program 
@@ -118,10 +92,10 @@ This function solves a linear program, using some exact
 Integral Domain `ET` for its computations. Various 
 options may be provided, see `Quadratic_program_options`. 
 
-Requirements 
--------------- 
+\tparam LinearProgram a model of `LinearProgram` such as `Quadratic_program<NT>`, 
+`Quadratic_program_from_mps<NT>`, or `Linear_program_from_iterators<A_it, B_it, R_it, FL_it, L_it, FU_it, U_it, C_it>`
 
-`ET` is a model of the concepts `IntegralDomain` and 
+\tparam ET a model of the concepts `IntegralDomain` and 
 `RealEmbeddable`; it must 
 be an exact type, and all entries of `qp` are convertible to 
 `ET`. 
@@ -144,15 +118,6 @@ factor. For maximum efficiency, it is advisable to define the macros
 \returns the solution of the linear program `lp`, solved
 with exact number type `ET`.
 
-Example 
--------------- 
-
-\ref QP_solver/first_lp.cpp 
-
-\sa `Quadratic_program<NT>` 
-\sa `Quadratic_program_from_mps<NT>` 
-\sa `Linear_program_from_iterators<A_it, B_it, R_it, FL_it, L_it, FU_it, U_it, C_it>` 
-
 */
 template <LinearProgram, ET>
 Quadratic_program_solution<ET> solve_linear_program 
@@ -165,10 +130,10 @@ This function solves a nonnegative linear program, using some exact
 Integral Domain `ET` for its computations. Various 
 options may be provided, see `Quadratic_program_options`. 
 
-Requirements 
--------------- 
+\tparam NonnegativeQuadraticProgram a model of `NonnegativeQuadraticProgram` such as
+`Quadratic_program<NT>`, `Quadratic_program_from_mps<NT>`, or `Nonnegative_quadratic_program_from_iterators<A_it, B_it, R_it, D_it, C_it>`.
 
-`ET` is a model of the concepts `IntegralDomain` and 
+\tparam ET is a model of the concepts `IntegralDomain` and 
 `RealEmbeddable`; it must 
 be an exact type, and all entries of `qp` are convertible to 
 `ET`. 
@@ -191,16 +156,6 @@ factor. For maximum efficiency, it is advisable to define the macros
 \returns the solution of the nonnegative linear program `lp`, solved
 with exact number type `ET`.
 
-Example 
--------------- 
-
-\ref QP_solver/first_nonnegative_lp.cpp 
-
-The models of \ref NonnegativeLinearProgram\:
-\sa `Quadratic_program<NT>` 
-\sa `Quadratic_program_from_mps<NT>` 
-\sa `Nonnegative_linear_program_from_iterators<A_it, B_it, R_it, C_it>` 
-
 */
 template <NonnegativeLinearProgram, ET>
 Quadratic_program_solution<ET> solve_nonnegative_linear_program 
@@ -213,10 +168,10 @@ This function solves a nonnegative quadratic program, using some exact
 Integral Domain `ET` for its computations. Various 
 options may be provided, see `Quadratic_program_options`. 
 
-Requirements 
--------------- 
+\tparam NonnegativeQuadraticProgram a model of `NonnegativeQuadraticProgram` such as
+`Quadratic_program<NT>`, `Quadratic_program_from_mps<NT>`, or `Nonnegative_quadratic_program_from_iterators<A_it, B_it, R_it, D_it, C_it>`.
 
-`ET` is a model of the concepts `IntegralDomain` and 
+\tparam ET a model of the concepts `IntegralDomain` and 
 `RealEmbeddable`; it must 
 be an exact type, and all entries of `qp` are convertible to 
 `ET`. 
@@ -239,15 +194,6 @@ factor. For maximum efficiency, it is advisable to define the macros
 \returns the solution of the nonnegative quadratic program `qp`, solved
 with exact number type `ET`.
 
-Example 
--------------- 
-
-\ref QP_solver/first_nonnegative_qp.cpp 
-
-The models of \ref ::NonnegativeQuadraticProgram\:
-\sa `Quadratic_program<NT>` 
-\sa `Quadratic_program_from_mps<NT>` 
-\sa `Nonnegative_quadratic_program_from_iterators<A_it, B_it, R_it, D_it, C_it>` 
 */
 template <NonnegativeQuadraticProgram, ET>
 Quadratic_program_solution<ET> solve_nonnegative_quadratic_program 
@@ -260,10 +206,11 @@ This function solves a quadratic program, using some exact
 Integral Domain `ET` for its computations. Various 
 options may be provided, see `Quadratic_program_options`. 
 
-Requirements 
--------------- 
 
-`ET` is a model of the concepts `IntegralDomain` and 
+\tparam NonnegativeQuadraticProgram a model of `NonnegativeQuadraticProgram` such as
+`Quadratic_program<NT>`, `Quadratic_program_from_mps<NT>`, or `Nonnegative_quadratic_program_from_iterators<A_it, B_it, R_it, D_it, C_it>`.
+
+\tparam ET is a model of the concepts `IntegralDomain` and 
 `RealEmbeddable`; it must 
 be an exact type, and all entries of `qp` are convertible to 
 `ET`. 
@@ -286,16 +233,6 @@ factor. For maximum efficiency, it is advisable to define the macros
 \returns the solution of the quadratic program `qp`, solved
 with exact number type `ET`.
 
-Example 
--------------- 
-
-\ref QP_solver/first_qp.cpp 
-
-The models of \ref QuadraticProgram\:
-\sa `Quadratic_program<NT>` 
-\sa `Quadratic_program_from_mps<NT>` 
-\sa `Quadratic_program_from_iterators<A_it, B_it, R_it, FL_it, L_it, FU_it, U_it, D_it, C_it>` 
-
 */
 template <QuadraticProgram, ET>
 Quadratic_program_solution<ET> solve_quadratic_program 
diff --git a/QP_solver/doc/QP_solver/CGAL/QP_models.h b/QP_solver/doc/QP_solver/CGAL/QP_models.h
index 1ede1029905..28748e7f103 100644
--- a/QP_solver/doc/QP_solver/CGAL/QP_models.h
+++ b/QP_solver/doc/QP_solver/CGAL/QP_models.h
@@ -57,8 +57,7 @@ namespace CGAL {
   \cgalModels `QuadraticProgram`
   \cgalModels `LinearProgram`
 
-  Example 
-  -------------- 
+  \cgalHeading{Example} 
 
   \ref QP_solver/first_lp_from_iterators.cpp 
 
@@ -113,8 +112,7 @@ public:
 
   \returns an instance of `Linear_program_from_iterators<A_it, B_it, R_it, FL_it, L_it, FU_it, U_it, C_it>`, constructed from the given iterators.
 
-  Example 
-  -------------- 
+  \cgalHeading{Example} 
 
   The following example demonstrates the typical usage of makers 
   with the simpler function `make_nonnegative_linear_program_from_iterators()`. 
@@ -160,8 +158,7 @@ make_linear_program_from_iterators (
   \returns an instance of `Nonnegative_linear_program_from_iterators<A_it, B_it, R_it, C_it>`, constructed from the given iterators.
 
 
-  Example 
-  -------------- 
+  \cgalHeading{Example}
 
   `QP_solver/solve_convex_hull_containment_lp2.h` 
 
@@ -198,8 +195,7 @@ make_nonnegative_linear_program_from_iterators (
   `Nonnegative_quadratic_program_from_iterators<A_it, B_it, R_it, D_it,C_it>`, 
   constructed from the given iterators.
 
-  Example 
-  -------------- 
+  \cgalHeading{Example} 
 
   The following example demonstrates the typical usage of makers 
   with the simpler function `make_nonnegative_linear_program_from_iterators()`. 
@@ -240,8 +236,7 @@ make_nonnegative_quadratic_program_from_iterators (
 
   \returns an instance of `Quadratic_program_from_iterators<A_it, B_it, R_it, FL_it, L_it, FU_it, U_it, D_it, C_it>`, constructed from the given iterators.
 
-  Example 
-  -------------- 
+  \cgalHeading{Example}
 
   The following example demonstrates the typical usage of makers 
   with the simpler function `make_nonnegative_linear_program_from_iterators()`. 
@@ -330,8 +325,7 @@ make_quadratic_program_from_iterators (
   \cgalModels `NonnegativeQuadraticProgram`
   \cgalModels `NonnegativeLinearProgram`
 
-  Example 
-  -------------- 
+  \cgalHeading{Example}
 
   \ref QP_solver/first_nonnegative_lp_from_iterators.cpp 
 
@@ -422,8 +416,7 @@ public:
   \cgalModels `QuadraticProgram`
   \cgalModels `NonnegativeQuadraticProgram`
 
-  Example 
-  -------------- 
+  \cgalHeading{Example}
 
   \ref QP_solver/first_nonnegative_qp_from_iterators.cpp 
 
@@ -522,8 +515,7 @@ public:
 
   \cgalModels `QuadraticProgram`
 
-  Example 
-  -------------- 
+  \cgalHeading{Example}
 
   \ref QP_solver/first_qp_from_iterators.cpp 
 
@@ -630,8 +622,7 @@ public:
   \cgalModels `NonnegativeQuadraticProgram`
   \cgalModels `NonnegativeLinearProgram`
 
-  Example 
-  -------------- 
+  \cgalHeading{Example} 
 
   \ref QP_solver/first_qp_from_mps.cpp 
 
@@ -854,8 +845,7 @@ namespace CGAL {
   \cgalModels `NonnegativeQuadraticProgram`
   \cgalModels `NonnegativeLinearProgram`
 
-  Example 
-  -------------- 
+  \cgalHeading{Example} 
 
   \ref QP_solver/first_qp.cpp 
 
diff --git a/QP_solver/doc/QP_solver/CGAL/QP_options.h b/QP_solver/doc/QP_solver/CGAL/QP_options.h
index 658363321e9..389f4dc81c2 100644
--- a/QP_solver/doc/QP_solver/CGAL/QP_options.h
+++ b/QP_solver/doc/QP_solver/CGAL/QP_options.h
@@ -15,8 +15,7 @@ options referring to
 </OL> 
 The idea is that this list grows in the future. 
 
-Operations 
--------------- 
+\cgalHeading{Operations}
 
 Here we just have set/get pairs for any option type. 
 
diff --git a/QP_solver/doc/QP_solver/CGAL/QP_solution.h b/QP_solver/doc/QP_solver/CGAL/QP_solution.h
index b3b9776dbfe..7130cadaf62 100644
--- a/QP_solver/doc/QP_solver/CGAL/QP_solution.h
+++ b/QP_solver/doc/QP_solver/CGAL/QP_solution.h
@@ -36,13 +36,11 @@ the four functions
 `solve_nonnegative_quadratic_program`, and 
 `solve_nonnegative_linear_program`. 
 
-Example 
--------------- 
+\cgalHeading{Example}
 
 \ref QP_solver/first_qp.cpp 
 
-Terminology 
--------------- 
+\cgalHeading{Terminology}
 
 If there is no \f$ \qpx\f$ that satisfies all the (in)equalities, 
 the program is called <I>infeasible</I>, otherwise, it is <I>feasible</I>, 

From ecc89af2b2eaeb6d3f275b6f20d7a0e2cede9dec Mon Sep 17 00:00:00 2001
From: Andreas Fabri <andreas.fabri@geometryfactory.com>
Date: Sun, 30 Dec 2018 14:57:54 +0100
Subject: [PATCH 21/35] Fix Surface_mesh::join()

---
 .../include/CGAL/Surface_mesh/Surface_mesh.h  | 79 +++++++++++++------
 1 file changed, 54 insertions(+), 25 deletions(-)

diff --git a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
index 45a7b529eec..8e9392bec14 100644
--- a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
+++ b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
@@ -1106,42 +1106,40 @@ public:
     }
     size_type inf_value = (std::numeric_limits<size_type>::max)();
     if(other.vertices_freelist_ != inf_value){
-      if(vertices_freelist_ != inf_value){
-        Vertex_index vi(nv+other.vertices_freelist_);
-        Halfedge_index inf((std::numeric_limits<size_type>::max)());
-        while(vconn_[vi].halfedge_ != inf){
-          Vertex_index corrected_vi = Vertex_index(size_type(vconn_[vi].halfedge_)+nv-nh);
-          vconn_[vi].halfedge_ = Halfedge_index(corrected_vi);
-          vi = corrected_vi;
-        }
-        vconn_[vi].halfedge_ = Halfedge_index(vertices_freelist_);
+      Vertex_index vi(nv+other.vertices_freelist_);
+      Halfedge_index inf((std::numeric_limits<size_type>::max)());
+      while(vconn_[vi].halfedge_ != inf){
+        Vertex_index corrected_vi = Vertex_index(size_type(vconn_[vi].halfedge_)+nv-nh);
+        vconn_[vi].halfedge_ = Halfedge_index(corrected_vi);
+        vi = corrected_vi;
       }
+      vconn_[vi].halfedge_ = Halfedge_index(vertices_freelist_);
+      
       vertices_freelist_ = nv + other.vertices_freelist_; 
     }
     if(other.faces_freelist_ != inf_value){
-      if(faces_freelist_ != inf_value){
-        Face_index fi(nf+other.faces_freelist_);
-        Halfedge_index inf((std::numeric_limits<size_type>::max)());
-        while(fconn_[fi].halfedge_ != inf){
-          Face_index corrected_fi = Face_index(size_type(fconn_[fi].halfedge_)+nf-nh);
-          fconn_[fi].halfedge_ = Halfedge_index(corrected_fi);
-          fi = corrected_fi;
-        }
-        fconn_[fi].halfedge_ = Halfedge_index(faces_freelist_);
+      Face_index fi(nf+other.faces_freelist_);
+      Halfedge_index inf((std::numeric_limits<size_type>::max)());
+      while(fconn_[fi].halfedge_ != inf){
+        Face_index corrected_fi = Face_index(size_type(fconn_[fi].halfedge_)+nf-nh);
+        fconn_[fi].halfedge_ = Halfedge_index(corrected_fi);
+        fi = corrected_fi;
       }
+      fconn_[fi].halfedge_ = Halfedge_index(faces_freelist_);
+
       faces_freelist_ = nf + other.faces_freelist_; 
     }
+    
     if(other.edges_freelist_ != inf_value){
-      if(edges_freelist_ != inf_value){
-        Halfedge_index hi(nh+other.edges_freelist_);
-        Halfedge_index inf((std::numeric_limits<size_type>::max)());
-        while(hconn_[hi].next_halfedge_ != inf){
-          hi = hconn_[hi].next_halfedge_;
-        }
-        hconn_[hi].next_halfedge_ = Halfedge_index(edges_freelist_);
+      Halfedge_index hi(nh+other.edges_freelist_);
+      Halfedge_index inf((std::numeric_limits<size_type>::max)());
+      while(hconn_[hi].next_halfedge_ != inf){
+        hi = hconn_[hi].next_halfedge_;
       }
+      hconn_[hi].next_halfedge_ = Halfedge_index(edges_freelist_);
       edges_freelist_ = nh + other.edges_freelist_; 
     }
+
     garbage_ = garbage_ || other.garbage_;
     removed_vertices_ += other.removed_vertices_;
     removed_edges_ += other.removed_edges_;
@@ -1378,6 +1376,37 @@ public:
         if(!valid && verbose){
           std::cerr << "#faces: iterated: " << fcount << " vs number_of_faces(): " << number_of_faces()<< std::endl;
         }
+
+        size_type inf = (std::numeric_limits<size_type>::max)();
+        size_type vfl = vertices_freelist_;
+        int rv = 0;
+        while(vfl != inf){
+          vfl = (size_type)vconn_[Vertex_index(vfl)].halfedge_;
+          rv++;
+        }
+        assert( rv == removed_vertices_ );
+        valid = valid && ( rv == removed_vertices_ );
+
+
+        size_type efl = edges_freelist_;
+        int re = 0;
+        while(efl != inf){
+          efl = (size_type)hconn_[Halfedge_index(efl)].next_halfedge_;
+          re++;
+        }
+        assert( re == removed_edges_ );
+        valid = valid && ( re == removed_edges_ );
+
+        size_type ffl = faces_freelist_;
+        int rf = 0;
+        while(ffl != inf){
+          ffl = (size_type)fconn_[Face_index(ffl)].halfedge_;
+          rf++;
+        }
+        assert( rf == removed_faces_ );
+        valid = valid && ( rf == removed_faces_ );
+
+  
         return valid;
     }
 

From 3bfc09de3fdb5a91d9ff1383ca0dc14b9f64b822 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= <sebastien.loriot@cgal.org>
Date: Mon, 31 Dec 2018 09:34:17 +0100
Subject: [PATCH 22/35] add comments

---
 .../include/CGAL/Surface_mesh/Surface_mesh.h  | 23 +++++++++++++++----
 1 file changed, 18 insertions(+), 5 deletions(-)

diff --git a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
index 8e9392bec14..573d7d9af91 100644
--- a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
+++ b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
@@ -1067,28 +1067,33 @@ public:
   
   bool join(const Surface_mesh& other)
   {
+    // increase capacity
     const size_type nv = num_vertices(), nh = num_halfedges(), nf = num_faces();
     resize(num_vertices()+  other.num_vertices(),
             num_edges()+  other.num_edges(),
             num_faces()+  other.num_faces());
 
+    // append properties in the free space created by resize
     vprops_.transfer(other.vprops_);
     hprops_.transfer(other.hprops_);
     fprops_.transfer(other.fprops_);
     eprops_.transfer(other.eprops_);
 
+    // translate halfedge index in vertex -> halfedge
     for(size_type i = nv; i < nv+other.num_vertices(); i++){
       Vertex_index vi(i);
       if(vconn_[vi].halfedge_ != null_halfedge()){
         vconn_[vi].halfedge_ = Halfedge_index(size_type(vconn_[vi].halfedge_)+nh);
       }
     }
+    // translate halfedge index in face -> halfedge
     for(size_type i = nf; i < nf+other.num_faces(); i++){
       Face_index fi(i);
       if(fconn_[fi].halfedge_ != null_halfedge()){
         fconn_[fi].halfedge_ = Halfedge_index(size_type(fconn_[fi].halfedge_)+nh);
       }
     }
+    // translate indices in halfedge -> face, halfedge -> target, halfedge -> prev, and halfedge -> next
     for(size_type i = nh; i < nh+other.num_halfedges(); i++){
       Halfedge_index hi(i);
       if(hconn_[hi].face_ != null_face()){
@@ -1105,41 +1110,50 @@ public:
       }
     }
     size_type inf_value = (std::numeric_limits<size_type>::max)();
+
+    // merge vertex free list
     if(other.vertices_freelist_ != inf_value){
       Vertex_index vi(nv+other.vertices_freelist_);
       Halfedge_index inf((std::numeric_limits<size_type>::max)());
+      // correct the indices in the linked list of free vertices copied (due to vconn_ translation)
       while(vconn_[vi].halfedge_ != inf){
         Vertex_index corrected_vi = Vertex_index(size_type(vconn_[vi].halfedge_)+nv-nh);
         vconn_[vi].halfedge_ = Halfedge_index(corrected_vi);
         vi = corrected_vi;
       }
+      // append the vertex free linked list of `this` to the copy of `other`
       vconn_[vi].halfedge_ = Halfedge_index(vertices_freelist_);
-      
+      // update the begin of the vertex free linked list
       vertices_freelist_ = nv + other.vertices_freelist_; 
     }
+    // merge face free list
     if(other.faces_freelist_ != inf_value){
       Face_index fi(nf+other.faces_freelist_);
       Halfedge_index inf((std::numeric_limits<size_type>::max)());
+      // correct the indices in the linked list of free faces copied (due to fconn_ translation)
       while(fconn_[fi].halfedge_ != inf){
         Face_index corrected_fi = Face_index(size_type(fconn_[fi].halfedge_)+nf-nh);
         fconn_[fi].halfedge_ = Halfedge_index(corrected_fi);
         fi = corrected_fi;
       }
+      // append the face free linked list of `this` to the copy of `other`
       fconn_[fi].halfedge_ = Halfedge_index(faces_freelist_);
-
+      // update the begin of the face free linked list
       faces_freelist_ = nf + other.faces_freelist_; 
     }
-    
+    // merge edge free list
     if(other.edges_freelist_ != inf_value){
       Halfedge_index hi(nh+other.edges_freelist_);
       Halfedge_index inf((std::numeric_limits<size_type>::max)());
       while(hconn_[hi].next_halfedge_ != inf){
         hi = hconn_[hi].next_halfedge_;
       }
+      // append the halfedge free linked list of `this` to the copy of `other`
       hconn_[hi].next_halfedge_ = Halfedge_index(edges_freelist_);
+      // update the begin of the halfedge free linked list
       edges_freelist_ = nh + other.edges_freelist_; 
     }
-
+    // update garbage infos
     garbage_ = garbage_ || other.garbage_;
     removed_vertices_ += other.removed_vertices_;
     removed_edges_ += other.removed_edges_;
@@ -1406,7 +1420,6 @@ public:
         assert( rf == removed_faces_ );
         valid = valid && ( rf == removed_faces_ );
 
-  
         return valid;
     }
 

From dc56732a6f6b38ab775d9ae3e69911962ce05f89 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= <sebastien.loriot@cgal.org>
Date: Wed, 2 Jan 2019 08:46:40 +0100
Subject: [PATCH 23/35] fix typo

---
 .../doc/Arrangement_on_surface_2/Arrangement_on_surface_2.txt   | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Arrangement_on_surface_2/doc/Arrangement_on_surface_2/Arrangement_on_surface_2.txt b/Arrangement_on_surface_2/doc/Arrangement_on_surface_2/Arrangement_on_surface_2.txt
index ad8d8120b55..dc4d957e931 100644
--- a/Arrangement_on_surface_2/doc/Arrangement_on_surface_2/Arrangement_on_surface_2.txt
+++ b/Arrangement_on_surface_2/doc/Arrangement_on_surface_2/Arrangement_on_surface_2.txt
@@ -2681,7 +2681,7 @@ manner whenever possible. Thus, it resorts to exact computations only when
 the approximate computation fails to produce an unambiguous result.
 Note that most arrangement vertices are therefore associated with approximated
 points. You cannot access the coordinates of such points and obtain them as
-algebraic numbers, and only access to the approximate coordinates in possible.
+algebraic numbers, and only access to the approximate coordinates is possible.
 See the Reference Manual for the exact interface of the `Point_2`,
 `Curve_2` and `X_monotone_curve_2` defined by the traits class.
 

From 5a1b51bdd05898a1e1718b7f70747055cbb6a79b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= <sebastien.loriot@cgal.org>
Date: Wed, 2 Jan 2019 10:11:28 +0100
Subject: [PATCH 24/35] fix warning

---
 Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
index 573d7d9af91..39437aff178 100644
--- a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
+++ b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
@@ -1393,31 +1393,28 @@ public:
 
         size_type inf = (std::numeric_limits<size_type>::max)();
         size_type vfl = vertices_freelist_;
-        int rv = 0;
+        size_type rv = 0;
         while(vfl != inf){
           vfl = (size_type)vconn_[Vertex_index(vfl)].halfedge_;
           rv++;
         }
-        assert( rv == removed_vertices_ );
         valid = valid && ( rv == removed_vertices_ );
 
 
         size_type efl = edges_freelist_;
-        int re = 0;
+        size_type re = 0;
         while(efl != inf){
           efl = (size_type)hconn_[Halfedge_index(efl)].next_halfedge_;
           re++;
         }
-        assert( re == removed_edges_ );
         valid = valid && ( re == removed_edges_ );
 
         size_type ffl = faces_freelist_;
-        int rf = 0;
+        size_type rf = 0;
         while(ffl != inf){
           ffl = (size_type)fconn_[Face_index(ffl)].halfedge_;
           rf++;
         }
-        assert( rf == removed_faces_ );
         valid = valid && ( rf == removed_faces_ );
 
         return valid;

From ace09977af4da7aab21f10c94a56209607966cab Mon Sep 17 00:00:00 2001
From: Efi Fogel <efifogel@gmail.com>
Date: Wed, 2 Jan 2019 12:22:58 +0200
Subject: [PATCH 25/35] Cleaned up the code

---
 .../Arr_polyhedral_sgm.h                           | 14 ++++++--------
 1 file changed, 6 insertions(+), 8 deletions(-)

diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h b/Arrangement_on_surface_2/include/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h
index 603837c1ebf..8a57cd28808 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h
@@ -242,13 +242,10 @@ private:
     /*! Set the marked-face index */
     void set_marked_facet_index(size_type id) {m_marked_facet_index = id;}
 
-    /*! Add all facets. */
+    /*! Add vertices to the current facet. */
     template <typename Iterator, typename Builder>
-    void add_facets(Iterator begin, Iterator end, Builder& B)
-    {
-      for (Iterator it = begin; it != end; ++it) B.add_vertex_to_facet(*it);
-      B.end_facet();
-    }
+    void add_vertices_to_facet(Iterator begin, Iterator end, Builder& B)
+    { for (Iterator it = begin; it != end; ++it) B.add_vertex_to_facet(*it); }
 
     /*! builds the polyhedron */
     void operator()(HDS& hds)
@@ -271,10 +268,11 @@ private:
         Polyhedron_facet_handle fh = B.begin_facet();
         if (counter == m_marked_facet_index) fh->set_marked(true);
         //! \todo EF: when upgrading to C++11 enable the following code and
-        // remove add_facets().
+        // remove add_vertices_to_facet().
         // for (const auto& facet : *it) B.add_vertex_to_facet(facet);
         // B.end_facet();
-        add_facets(it->begin(), it->end(), B);
+        add_vertices_to_facet(it->begin(), it->end(), B);
+        B.end_facet();
         ++counter;
       }
       B.end_surface();

From 0396dd11182d83ab43e0ef27ab70068086388b56 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 19 Sep 2018 13:49:25 +0200
Subject: [PATCH 26/35] Patch to improve normal orientation using user-defined
 seed points

---
 .../CGAL/boost/graph/named_params_helper.h    |  23 +++
 .../CGAL/boost/graph/parameters_interface.h   |   1 +
 .../include/CGAL/mst_orient_normals.h         | 142 ++++++++++++++----
 3 files changed, 138 insertions(+), 28 deletions(-)

diff --git a/BGL/include/CGAL/boost/graph/named_params_helper.h b/BGL/include/CGAL/boost/graph/named_params_helper.h
index cbf5a2fedce..e585269eba0 100644
--- a/BGL/include/CGAL/boost/graph/named_params_helper.h
+++ b/BGL/include/CGAL/boost/graph/named_params_helper.h
@@ -405,6 +405,29 @@ namespace CGAL {
         > ::type  type;
     };
 
+    template<typename PointRange, typename NamedParameters>
+    class GetIsConstrainedMap
+    {
+      struct DummyConstrainedMap
+      {
+        typedef typename std::iterator_traits<typename PointRange::iterator>::value_type key_type;
+        typedef bool value_type;
+        typedef value_type reference;
+        typedef boost::readable_property_map_tag category;
+
+        typedef DummyConstrainedMap Self;
+        friend reference get(const Self&, const key_type&) { return false; }
+      };
+
+    public:
+      typedef DummyConstrainedMap NoMap;
+      typedef typename boost::lookup_named_param_def <
+        internal_np::point_is_constrained_t,
+        NamedParameters,
+        DummyConstrainedMap //default
+        > ::type  type;
+    };
+
   } // namespace Point_set_processing_3
   
   template<typename NamedParameters, typename DefaultSolver>
diff --git a/BGL/include/CGAL/boost/graph/parameters_interface.h b/BGL/include/CGAL/boost/graph/parameters_interface.h
index d2dd2208420..9e38eec688d 100644
--- a/BGL/include/CGAL/boost/graph/parameters_interface.h
+++ b/BGL/include/CGAL/boost/graph/parameters_interface.h
@@ -117,3 +117,4 @@ CGAL_add_named_parameter(plane_t, plane_map, plane_map)
 CGAL_add_named_parameter(plane_index_t, plane_index_map, plane_index_map)
 CGAL_add_named_parameter(select_percentage_t, select_percentage, select_percentage)
 CGAL_add_named_parameter(require_uniform_sampling_t, require_uniform_sampling, require_uniform_sampling)
+CGAL_add_named_parameter(point_is_constrained_t, point_is_constrained, point_is_constrained_map)
diff --git a/Point_set_processing_3/include/CGAL/mst_orient_normals.h b/Point_set_processing_3/include/CGAL/mst_orient_normals.h
index 4935cd7d208..7f9789d5e2b 100644
--- a/Point_set_processing_3/include/CGAL/mst_orient_normals.h
+++ b/Point_set_processing_3/include/CGAL/mst_orient_normals.h
@@ -121,6 +121,34 @@ public:
     const NormalMap m_normal_map;
 };
 
+template <typename ForwardIterator>
+class Default_constrained_map
+{
+public:
+
+  typedef boost::readable_property_map_tag     category;
+  typedef typename ForwardIterator::value_type key_type;
+  typedef bool                                 value_type;
+  typedef value_type                           reference;
+
+private:
+  
+  ForwardIterator m_source_point;
+
+public:
+
+  Default_constrained_map () { }
+  Default_constrained_map (ForwardIterator source_point)
+    : m_source_point (source_point) { }
+
+  /// Free function to access the map elements.
+  friend inline
+  reference get(const Default_constrained_map& map, key_type p)
+  {
+    return (p == *map.m_source_point);
+  }
+
+};
 
 /// Helper class: Propagate_normal_orientation
 ///
@@ -145,10 +173,12 @@ struct Propagate_normal_orientation
   : public boost::base_visitor< Propagate_normal_orientation<ForwardIterator, NormalMap, Kernel> >
 {
     typedef internal::MST_graph<ForwardIterator, NormalMap, Kernel> MST_graph;
+    typedef typename MST_graph::vertex_descriptor vertex_descriptor;
     typedef boost::on_examine_edge event_filter;
 
-    Propagate_normal_orientation(double angle_max = CGAL_PI/2.) ///< max angle to propagate the normal orientation (radians)
-    : m_angle_max(angle_max)
+    Propagate_normal_orientation(vertex_descriptor source,
+                               double angle_max = CGAL_PI/2.) ///< max angle to propagate the normal orientation (radians)
+      : m_source(source), m_angle_max(angle_max)
     {
         // Precondition: 0 < angle_max <= PI/2
         CGAL_point_set_processing_precondition(0 < angle_max && angle_max <= CGAL_PI/2.);
@@ -158,16 +188,28 @@ struct Propagate_normal_orientation
     void operator()(Edge& edge, const MST_graph& mst_graph)
     {
         typedef typename boost::property_traits<NormalMap>::reference Vector_ref;
-        typedef typename MST_graph::vertex_descriptor vertex_descriptor;
+
+        // Gets source
+        vertex_descriptor source_vertex = source(edge, mst_graph);
+
+        // Gets target
+        vertex_descriptor target_vertex = target(edge, mst_graph);
+        bool& target_normal_is_oriented = ((MST_graph&)mst_graph)[target_vertex].is_oriented;
+        
+        // special case if vertex is source vertex (and thus has no related point/normal)
+        if (source_vertex == m_source)
+        {
+          target_normal_is_oriented = true;
+          return;
+        }
 
         // Gets source normal
-        vertex_descriptor source_vertex = source(edge, mst_graph);
         Vector_ref source_normal = get(mst_graph.m_normal_map, *(mst_graph[source_vertex].input_point) );
         const bool source_normal_is_oriented = mst_graph[source_vertex].is_oriented;
-        // Gets target normal
-        vertex_descriptor target_vertex = target(edge, mst_graph);
+
+        // Gets target
         Vector_ref target_normal = get( mst_graph.m_normal_map, *(mst_graph[target_vertex].input_point) );
-        bool& target_normal_is_oriented = ((MST_graph&)mst_graph)[target_vertex].is_oriented;
+        
         if ( ! target_normal_is_oriented )
         {
           //             ->                        ->
@@ -188,6 +230,7 @@ struct Propagate_normal_orientation
 // Data
 // Implementation note: boost::breadth_first_search() makes copies of this object => data must be constant or shared.
 private:
+    vertex_descriptor m_source;
     const double m_angle_max; ///< max angle to propagate the normal orientation (radians).
 };
 
@@ -264,6 +307,7 @@ template <typename ForwardIterator,
           typename PointMap,
           typename NormalMap,
           typename IndexMap,
+          typename ConstrainedMap,
           typename Kernel
 >
 Riemannian_graph<ForwardIterator>
@@ -273,6 +317,7 @@ create_riemannian_graph(
     PointMap point_map, ///< property map: value_type of ForwardIterator -> Point_3
     NormalMap normal_map, ///< property map: value_type of ForwardIterator -> Vector_3
     IndexMap index_map, ///< property map ForwardIterator -> index
+    ConstrainedMap constrained_map, ///< property map ForwardIterator -> bool
     unsigned int k, ///< number of neighbors
     const Kernel& /*kernel*/) ///< geometric traits.
 {
@@ -337,6 +382,12 @@ create_riemannian_graph(
         CGAL_point_set_processing_assertion(v == get(index_map,it));
         riemannian_graph[v].input_point = it;
     }
+
+    // add source vertex (virtual, does not correspond to a point)
+    typename Riemannian_graph::vertex_descriptor v = add_vertex(riemannian_graph);
+    std::size_t source_point_index = num_input_points;
+    CGAL_point_set_processing_assertion(v == source_point_index);
+
     //
     // add edges
     Riemannian_graph_weight_map riemannian_graph_weight_map = get(boost::edge_weight, riemannian_graph);
@@ -384,6 +435,20 @@ create_riemannian_graph(
 
             search_iterator++;
         }
+
+        // Check if point is source
+        if (get(constrained_map, *it))
+        {
+          typename boost::graph_traits<Riemannian_graph>::edge_descriptor e;
+          bool inserted;
+          boost::tie(e, inserted) = add_edge(vertex(it_index, riemannian_graph),
+                                             vertex(source_point_index, riemannian_graph),
+                                             riemannian_graph);
+          CGAL_point_set_processing_assertion(inserted);
+
+          riemannian_graph_weight_map[e] = 0.;
+        }
+            
     }
 
     return riemannian_graph;
@@ -418,8 +483,7 @@ create_mst_graph(
     IndexMap index_map, ///< property map ForwardIterator -> index
     unsigned int k, ///< number of neighbors
     const Kernel& kernel, ///< geometric traits.
-    const Riemannian_graph<ForwardIterator>& riemannian_graph, ///< graph connecting each vertex to its knn
-    ForwardIterator source_point) ///< source point (with an oriented normal)
+    const Riemannian_graph<ForwardIterator>& riemannian_graph) ///< graph connecting each vertex to its knn
 {
     // prevents warnings
     CGAL_USE(point_map);
@@ -440,16 +504,17 @@ create_mst_graph(
     CGAL_point_set_processing_precondition(first != beyond);
 
     // Number of input points
-    const std::size_t num_input_points = num_vertices(riemannian_graph);
+    const std::size_t num_input_points = num_vertices(riemannian_graph) - 1;
 
     std::size_t memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE("  %ld Mb allocated\n", memory>>20);
     CGAL_TRACE("  Calls boost::prim_minimum_spanning_tree()\n");
 
     // Computes Minimum Spanning Tree.
-    std::size_t source_point_index = get(index_map, source_point);
+    std::size_t source_point_index = num_input_points;
+    
     Riemannian_graph_weight_map riemannian_graph_weight_map = get(boost::edge_weight, riemannian_graph);
     typedef std::vector<typename Riemannian_graph::vertex_descriptor> PredecessorMap;
-    PredecessorMap predecessor(num_input_points);
+    PredecessorMap predecessor(num_input_points + 1);
     boost::prim_minimum_spanning_tree(riemannian_graph, &predecessor[0],
                                       weight_map( riemannian_graph_weight_map )
                                      .root_vertex( vertex(source_point_index, riemannian_graph) ));
@@ -472,8 +537,13 @@ create_mst_graph(
         typename MST_graph::vertex_descriptor v = add_vertex(mst_graph);
         CGAL_point_set_processing_assertion(v == get(index_map,it));
         mst_graph[v].input_point = it;
-        mst_graph[v].is_oriented = (it == source_point);
+        mst_graph[v].is_oriented = false;
     }
+    
+    typename MST_graph::vertex_descriptor v = add_vertex(mst_graph);
+    CGAL_point_set_processing_assertion(v == source_point_index);
+    mst_graph[v].is_oriented = true;
+    
     // add edges
     for (std::size_t i=0; i < predecessor.size(); i++) // add edges
     {
@@ -525,6 +595,11 @@ create_mst_graph(
      If this parameter is omitted, `CGAL::Identity_property_map<geom_traits::Point_3>` is used.\cgalParamEnd
      \cgalParamBegin{normal_map} a model of `ReadWritePropertyMap` with value type
      `geom_traits::Vector_3`.\cgalParamEnd
+     \cgalParamBegin{point_is_constrained_map} a model of `ReadablePropertyMap` with value type
+     `bool`. Points with a `true` value will be used as seed points: their normal will be considered as already
+     oriented, it won't be altered and it will be propagated to its neighbors. If this parameter is omitted, 
+     the heighest point (heightest Z coordinate) will be used as the unique seed with an upward oriented
+     normal\cgalParamEnd
      \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
    \cgalNamedParamsEnd
 
@@ -545,6 +620,7 @@ mst_orient_normals(
     typedef typename Point_set_processing_3::GetPointMap<PointRange, NamedParameters>::type PointMap;
     typedef typename Point_set_processing_3::GetNormalMap<PointRange, NamedParameters>::type NormalMap;
     typedef typename Point_set_processing_3::GetK<PointRange, NamedParameters>::Kernel Kernel;
+    typedef typename Point_set_processing_3::GetIsConstrainedMap<PointRange, NamedParameters>::type ConstrainedMap;
 
     CGAL_static_assertion_msg(!(boost::is_same<NormalMap,
                                 typename Point_set_processing_3::GetNormalMap<PointRange, NamedParameters>::NoMap>::value),
@@ -552,6 +628,7 @@ mst_orient_normals(
 
     PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap());
     NormalMap normal_map = choose_param(get_param(np, internal_np::normal_map), NormalMap());
+    ConstrainedMap constrained_map = choose_param(get_param(np, internal_np::point_is_constrained), ConstrainedMap());
     Kernel kernel;
 
   // Bring private stuff to scope
@@ -584,37 +661,46 @@ mst_orient_normals(
     // and get() requires a lookup in the map.
     IndexMap index_map(points.begin(), points.end());
 
-    // Orients the normal of the point with maximum Z towards +Z axis.
-    typename PointRange::iterator source_point
-      = mst_find_source(points.begin(), points.end(),
-                        point_map, normal_map,
-                        kernel);
-
     // Iterates over input points and creates Riemannian Graph:
     // - vertices are numbered like the input points index.
     // - vertices are empty.
     // - we add the edge (i, j) if either vertex i is in the k-neighborhood of vertex j,
     //   or vertex j is in the k-neighborhood of vertex i.
-    Riemannian_graph riemannian_graph
-      = create_riemannian_graph(points.begin(), points.end(),
-                                point_map, normal_map, index_map,
-                                k,
-                                kernel);
+    Riemannian_graph riemannian_graph;
+
+    if (boost::is_same<ConstrainedMap,
+        typename CGAL::Point_set_processing_3::GetIsConstrainedMap<PointRange, NamedParameters>::NoMap>::value)
+      riemannian_graph = create_riemannian_graph(points.begin(), points.end(),
+                                                 point_map, normal_map, index_map,
+                                                 Default_constrained_map<typename PointRange::iterator>
+                                                 (mst_find_source(points.begin(), points.end(),
+                                                                  point_map, normal_map,
+                                                                  kernel)),
+                                                 k,
+                                                 kernel);
+    else
+      riemannian_graph = create_riemannian_graph(points.begin(), points.end(),
+                                                 point_map, normal_map, index_map,
+                                                 constrained_map,
+                                                 k,
+                                                 kernel);
 
     // Creates a Minimum Spanning Tree starting at source_point
     MST_graph mst_graph = create_mst_graph(points.begin(), points.end(),
                                            point_map, normal_map, index_map,
                                            k,
                                            kernel,
-                                           riemannian_graph,
-                                           source_point);
+                                           riemannian_graph);
 
     memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE("  %ld Mb allocated\n", memory>>20);
     CGAL_TRACE("  Calls boost::breadth_first_search()\n");
 
+    const std::size_t num_input_points = distance(points.begin(), points.end());
+    std::size_t source_point_index = num_input_points;
+    
     // Traverse the point set along the MST to propagate source_point's orientation
-    Propagate_normal_orientation<typename PointRange::iterator, NormalMap, Kernel> orienter;
-    std::size_t source_point_index = get(index_map, source_point);
+    Propagate_normal_orientation<typename PointRange::iterator, NormalMap, Kernel> orienter(source_point_index);
+
     boost::breadth_first_search(mst_graph,
                                 vertex(source_point_index, mst_graph), // source
                                 visitor(boost::make_bfs_visitor(orienter)));

From 8d1a0eea9ac3d73f9eb9586ee128b6a59399b5a7 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 19 Sep 2018 12:36:02 +0200
Subject: [PATCH 27/35] Specialize QMultiInputDialog to handle correctly radio
 buttons

---
 .../Polyhedron/include/QMultipleInputDialog.h    | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/Polyhedron/demo/Polyhedron/include/QMultipleInputDialog.h b/Polyhedron/demo/Polyhedron/include/QMultipleInputDialog.h
index 17147651f29..4825f59a6c0 100644
--- a/Polyhedron/demo/Polyhedron/include/QMultipleInputDialog.h
+++ b/Polyhedron/demo/Polyhedron/include/QMultipleInputDialog.h
@@ -9,6 +9,7 @@
 #include <QDialog>
 #include <QFormLayout>
 #include <QDialogButtonBox>
+#include <QRadioButton>
 
 class QMultipleInputDialog
 {
@@ -29,8 +30,19 @@ public:
   template <typename QObjectType>
   QObjectType* add (const char* name)
   {
-    QObjectType* out = new QObjectType (dialog);
-    form->addRow (QString(name), out);
+    QObjectType* out = NULL;
+        
+    if (boost::is_same<QObjectType, QRadioButton>::value)
+    {
+      out = dynamic_cast<QObjectType*>(new QRadioButton (QString(name), dialog));
+      form->addRow (out);
+    }
+    else
+    {
+      out = new QObjectType (dialog);
+      form->addRow (QString(name), out);
+    }
+    
     return out;
   }
 

From c67bc374adfb9ea250296ecfb8251363ab251253 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 19 Sep 2018 12:36:36 +0200
Subject: [PATCH 28/35] Orient normal plugin with selected seed + separation
 from estimation

---
 .../Point_set_normal_estimation_plugin.cpp    | 169 +++++++++++++-----
 .../Point_set_normal_estimation_plugin.ui     | 101 +----------
 2 files changed, 124 insertions(+), 146 deletions(-)

diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp
index ab53c6c887b..f0cd446b369 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp
@@ -18,6 +18,8 @@
 #include <QInputDialog>
 #include <QMessageBox>
 
+#include <QMultipleInputDialog.h>
+
 #include "run_with_qprogressdialog.h"
 
 #include "ui_Point_set_normal_estimation_plugin.h"
@@ -77,6 +79,25 @@ struct Jet_estimate_normals_functor
   }
 };
 
+struct Vector_to_pmap
+{
+  typedef boost::readable_property_map_tag     category;
+  typedef typename Point_set::Index            key_type;
+  typedef bool                                 value_type;
+  typedef value_type                           reference;
+  
+  std::vector<bool>* vec;
+
+  Vector_to_pmap (std::vector<bool>* vec = NULL) : vec (vec) { }
+
+  friend inline
+  reference get(const Vector_to_pmap& map, key_type p)
+  {
+    return (*map.vec)[p];
+  }
+
+};
+
 using namespace CGAL::Three;
 
 class Polyhedron_demo_point_set_normal_estimation_plugin :
@@ -88,6 +109,7 @@ class Polyhedron_demo_point_set_normal_estimation_plugin :
   Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0")
 
   QAction* actionNormalEstimation;
+  QAction* actionNormalOrientation;
   QAction* actionNormalInversion;
 
 public:
@@ -99,6 +121,10 @@ public:
     actionNormalEstimation->setObjectName("actionNormalEstimation");
     actionNormalEstimation->setProperty("subMenuName","Point Set Processing");
 
+    actionNormalOrientation = new QAction(tr("Normal Orientation"), mainWindow);
+    actionNormalOrientation->setObjectName("actionNormalOrientation");
+    actionNormalOrientation->setProperty("subMenuName","Point Set Processing");
+
     actionNormalInversion = new QAction(tr("Inverse Normal Orientations"), mainWindow);
     actionNormalInversion->setObjectName("actionNormalInversion");
     actionNormalInversion->setProperty("subMenuName","Point Set Processing");
@@ -106,7 +132,7 @@ public:
   }
 
   QList<QAction*> actions() const {
-    return QList<QAction*>() << actionNormalEstimation << actionNormalInversion;
+    return QList<QAction*>() << actionNormalEstimation << actionNormalOrientation << actionNormalInversion;
   }
 
   bool applicable(QAction* action) const {
@@ -124,6 +150,7 @@ public:
 
 public Q_SLOTS:
   void on_actionNormalEstimation_triggered();
+  void on_actionNormalOrientation_triggered();
   void on_actionNormalInversion_triggered();
 
 }; // end PS_demo_smoothing_plugin
@@ -142,7 +169,6 @@ class Point_set_demo_normal_estimation_dialog : public QDialog, private Ui::Norm
   unsigned int convolution_neighbors() const { return m_convolution_neighbors->value(); }
   double convolution_radius() const { return m_convolution_radius->value(); }
   double offset_radius() const { return m_offset_radius->value(); }
-  int orient_neighbors() const { return m_orient_neighbors->value(); }
 
   unsigned int method () const
   {
@@ -151,10 +177,6 @@ class Point_set_demo_normal_estimation_dialog : public QDialog, private Ui::Norm
     if (buttonVCM->isChecked ())       return 2;
     return -1;
   }
-  bool orient () const
-  {
-    return buttonOrient->isChecked();
-  }
   bool use_convolution_radius () const
   {
     return buttonRadius->isChecked();
@@ -275,59 +297,112 @@ void Polyhedron_demo_point_set_normal_estimation_plugin::on_actionNormalEstimati
                                                   << (memory>>20) << " Mb allocated"
                                                   << std::endl;
     }
+
     item->resetMenu();
     item->setRenderingMode(ShadedPoints);
 
+    // Updates scene
+    item->invalidateOpenGLBuffers();
+    scene->itemChanged(index);
+
+    QApplication::restoreOverrideCursor();
+  }
+#endif // !CGAL_DISABLE_NORMAL_ESTIMATION_PLUGIN
+}
+
+void Polyhedron_demo_point_set_normal_estimation_plugin::on_actionNormalOrientation_triggered()
+{
+  const CGAL::Three::Scene_interface::Item_id index = scene->mainSelectionIndex();
+
+  Scene_points_with_normal_item* item =
+    qobject_cast<Scene_points_with_normal_item*>(scene->item(index));
+
+  if(item)
+  {
+    // Gets point set
+    Point_set* points = item->point_set();
+    if(points == NULL)
+        return;
+
+    // Gets options
+    QMultipleInputDialog dialog ("Normal Orientation", mw);
+    QSpinBox* neighborhood = dialog.add<QSpinBox> ("Neighborhood Size: ");
+    neighborhood->setRange (1, 10000000);
+    neighborhood->setValue (18);
+
+    QRadioButton* use_seed_points = NULL;
+    QRadioButton* orient_selection = NULL;
+
+    if (points->nb_selected_points() != 0)
+    {
+      use_seed_points = dialog.add<QRadioButton> ("Use selection as seed points and orient the unselected points");
+      use_seed_points->setChecked(true);
+      orient_selection = dialog.add<QRadioButton> ("Orient selection");
+      orient_selection->setChecked(false);
+    }
+    
+    if(!dialog.exec())
+      return;
+      
+    QApplication::setOverrideCursor(Qt::BusyCursor);
+    QApplication::processEvents();
+
+    // First point to delete
+    Point_set::iterator first_unoriented_point = points->end();
+
     //***************************************
     // normal orientation
     //***************************************
 
-    if (dialog.orient ())
-      {
-        CGAL::Timer task_timer; task_timer.start();
-        std::cerr << "Orient normals with a Minimum Spanning Tree (k=" << dialog.orient_neighbors() << ")...\n";
+    CGAL::Timer task_timer; task_timer.start();
+    std::cerr << "Orient normals with a Minimum Spanning Tree (k=" << neighborhood->value() << ")...\n";
 
-        // Tries to orient normals
-        first_unoriented_point =
-          CGAL::mst_orient_normals(points->all_or_selection_if_not_empty(),
-                                   dialog.orient_neighbors(),
-                                   points->parameters());
+    // Tries to orient normals
+    if (points->nb_selected_points() != 0 && use_seed_points->isChecked())
+    {
+      std::vector<bool> constrained_map (points->size(), false);
 
-        std::size_t nb_unoriented_normals = std::distance(first_unoriented_point, points->end());
-        std::size_t memory = CGAL::Memory_sizer().virtual_size();
-        std::cerr << "Orient normals: " << nb_unoriented_normals << " point(s) with an unoriented normal are selected ("
-                  << task_timer.time() << " seconds, "
-                  << (memory>>20) << " Mb allocated)"
-                  << std::endl;
-
-        // Selects points with an unoriented normal
-        points->set_first_selected (first_unoriented_point);
-
-        // Updates scene
-        item->invalidateOpenGLBuffers();
-        scene->itemChanged(index);
-
-        QApplication::restoreOverrideCursor();
-
-        // Warns user
-        if (nb_unoriented_normals > 0)
-          {
-            QMessageBox::information(NULL,
-                                     tr("Points with an unoriented normal"),
-                                     tr("%1 point(s) with an unoriented normal are selected.\nPlease orient them or remove them before running Poisson reconstruction.")
-                                     .arg(nb_unoriented_normals));
-          }
-      }
+      for (typename Point_set::iterator it = points->first_selected(); it != points->end(); ++ it)
+        constrained_map[*it] = true;
+              
+      first_unoriented_point =
+        CGAL::mst_orient_normals(*points,
+                                 std::size_t(neighborhood->value()),
+                                 points->parameters().
+                                 point_is_constrained_map(Vector_to_pmap(&constrained_map)));
+    }
     else
-      {
-        // Updates scene
-        item->invalidateOpenGLBuffers();
-        scene->itemChanged(index);
+      first_unoriented_point =
+        CGAL::mst_orient_normals(points->all_or_selection_if_not_empty(),
+                                 std::size_t(neighborhood->value()),
+                                 points->parameters());
 
-        QApplication::restoreOverrideCursor();
-      }
+    std::size_t nb_unoriented_normals = std::distance(first_unoriented_point, points->end());
+    std::size_t memory = CGAL::Memory_sizer().virtual_size();
+    std::cerr << "Orient normals: " << nb_unoriented_normals << " point(s) with an unoriented normal are selected ("
+              << task_timer.time() << " seconds, "
+              << (memory>>20) << " Mb allocated)"
+              << std::endl;
+
+    // Selects points with an unoriented normal
+    points->set_first_selected (first_unoriented_point);
+
+    // Updates scene
+    item->invalidateOpenGLBuffers();
+    scene->itemChanged(index);
+
+    QApplication::restoreOverrideCursor();
+
+    // Warns user
+    if (nb_unoriented_normals > 0)
+    {
+      QMessageBox::information(NULL,
+                               tr("Points with an unoriented normal"),
+                               tr("%1 point(s) with an unoriented normal are selected.\nPlease orient them or remove them before running Poisson reconstruction.")
+                               .arg(nb_unoriented_normals));
+    }
   }
-#endif // !CGAL_DISABLE_NORMAL_ESTIMATION_PLUGIN
 }
 
+
 #include "Point_set_normal_estimation_plugin.moc"
diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.ui b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.ui
index 905ba695a9d..8ceaeaea887 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.ui
+++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.ui
@@ -6,8 +6,8 @@
    <rect>
     <x>0</x>
     <y>0</y>
-    <width>400</width>
-    <height>473</height>
+    <width>301</width>
+    <height>372</height>
    </rect>
   </property>
   <property name="windowTitle">
@@ -201,74 +201,6 @@
      </property>
     </spacer>
    </item>
-   <item>
-    <widget class="Line" name="line">
-     <property name="orientation">
-      <enum>Qt::Horizontal</enum>
-     </property>
-    </widget>
-   </item>
-   <item>
-    <spacer name="verticalSpacer_2">
-     <property name="orientation">
-      <enum>Qt::Vertical</enum>
-     </property>
-     <property name="sizeHint" stdset="0">
-      <size>
-       <width>20</width>
-       <height>40</height>
-      </size>
-     </property>
-    </spacer>
-   </item>
-   <item>
-    <widget class="QCheckBox" name="buttonOrient">
-     <property name="text">
-      <string>Orient Estimated Normals</string>
-     </property>
-     <property name="checked">
-      <bool>true</bool>
-     </property>
-    </widget>
-   </item>
-   <item>
-    <widget class="QFrame" name="frame_4">
-     <property name="frameShape">
-      <enum>QFrame::StyledPanel</enum>
-     </property>
-     <property name="frameShadow">
-      <enum>QFrame::Raised</enum>
-     </property>
-     <layout class="QFormLayout" name="formLayout_5">
-      <property name="fieldGrowthPolicy">
-       <enum>QFormLayout::AllNonFixedFieldsGrow</enum>
-      </property>
-      <item row="0" column="0">
-       <widget class="QLabel" name="neighborhoodSizeLabel">
-        <property name="text">
-         <string>Neighborhood Size</string>
-        </property>
-       </widget>
-      </item>
-      <item row="0" column="1">
-       <widget class="QSpinBox" name="m_orient_neighbors">
-        <property name="suffix">
-         <string> Nearest Neighbors</string>
-        </property>
-        <property name="minimum">
-         <number>6</number>
-        </property>
-        <property name="maximum">
-         <number>100000000</number>
-        </property>
-        <property name="value">
-         <number>18</number>
-        </property>
-       </widget>
-      </item>
-     </layout>
-    </widget>
-   </item>
    <item>
     <widget class="QDialogButtonBox" name="buttonBox">
      <property name="orientation">
@@ -279,19 +211,6 @@
      </property>
     </widget>
    </item>
-   <item>
-    <spacer name="verticalSpacer_3">
-     <property name="orientation">
-      <enum>Qt::Vertical</enum>
-     </property>
-     <property name="sizeHint" stdset="0">
-      <size>
-       <width>20</width>
-       <height>40</height>
-      </size>
-     </property>
-    </spacer>
-   </item>
   </layout>
  </widget>
  <resources/>
@@ -344,22 +263,6 @@
     </hint>
    </hints>
   </connection>
-  <connection>
-   <sender>buttonOrient</sender>
-   <signal>toggled(bool)</signal>
-   <receiver>frame_4</receiver>
-   <slot>setEnabled(bool)</slot>
-   <hints>
-    <hint type="sourcelabel">
-     <x>199</x>
-     <y>316</y>
-    </hint>
-    <hint type="destinationlabel">
-     <x>199</x>
-     <y>355</y>
-    </hint>
-   </hints>
-  </connection>
   <connection>
    <sender>buttonRadius</sender>
    <signal>toggled(bool)</signal>

From f406473fd3c47c5e613d1391f8eebaaa1ae78c45 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 19 Sep 2018 14:11:25 +0200
Subject: [PATCH 29/35] Document new named parameter

---
 .../doc/Point_set_processing_3/NamedParameters.txt       | 9 +++++++++
 Point_set_processing_3/include/CGAL/mst_orient_normals.h | 2 +-
 2 files changed, 10 insertions(+), 1 deletion(-)

diff --git a/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt b/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt
index dcb0f56f626..9d6702cdd6a 100644
--- a/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt
+++ b/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt
@@ -182,6 +182,15 @@ give better result if the distribution of the input points is highly non-uniform
 <b>Default value</b>: `false`
 \cgalNPEnd
 
+\cgalNPBegin{point_is_constrained_map} \anchor PSP_point_is_constrained_map
+is the property map containing information about points being constrained or not.
+Constrained points are left unaltered and are used as seeds in `mst_orient_normals()`.\n
+\b Type: a class model of `ReadablePropertyMap` with
+`PointRange::iterator::value_type` as key type and
+`bool` as value type. \n
+<b>Default value</b>: a property map with only the highest point constrained.
+\cgalNPEnd
+
 \cgalNPTableEnd
 
 */
diff --git a/Point_set_processing_3/include/CGAL/mst_orient_normals.h b/Point_set_processing_3/include/CGAL/mst_orient_normals.h
index 7f9789d5e2b..d0f981d7824 100644
--- a/Point_set_processing_3/include/CGAL/mst_orient_normals.h
+++ b/Point_set_processing_3/include/CGAL/mst_orient_normals.h
@@ -598,7 +598,7 @@ create_mst_graph(
      \cgalParamBegin{point_is_constrained_map} a model of `ReadablePropertyMap` with value type
      `bool`. Points with a `true` value will be used as seed points: their normal will be considered as already
      oriented, it won't be altered and it will be propagated to its neighbors. If this parameter is omitted, 
-     the heighest point (heightest Z coordinate) will be used as the unique seed with an upward oriented
+     the highest point (hightest Z coordinate) will be used as the unique seed with an upward oriented
      normal\cgalParamEnd
      \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
    \cgalNamedParamsEnd

From ad8ebeaa015ab05fb79af3cf03a1b8b8703e28bb Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 19 Sep 2018 14:12:40 +0200
Subject: [PATCH 30/35] Fix documentation of named parameters, replacing
 iterators by value_type

---
 .../Point_set_processing_3/NamedParameters.txt   | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt b/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt
index 9d6702cdd6a..ee321dfd95e 100644
--- a/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt
+++ b/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt
@@ -31,17 +31,17 @@ typename CGAL::Kernel_traits<point_map::value_type>::Kernel
 \cgalNPEnd
 
 \cgalNPBegin{point_map} \anchor PSP_point_map
- is the property map containing the points associated to the iterators of the point range `points`.\n
+ is the property map containing the points associated to the elements of the point range `points`.\n
 \b Type: a class model of `ReadablePropertyMap` with
-`PointRange::iterator` as key type and
+`PointRange::iterator::value_type` as key type and
 `geom_traits::Point_3` as value type. \n
 <b>Default value</b>: \code CGAL::Identity_property_map<geom_traits::Point_3>\endcode
 \cgalNPEnd
 
 \cgalNPBegin{normal_map} \anchor PSP_normal_map
- is the property map containing the normal vectors associated to the iterators of the point range `points`.\n
+ is the property map containing the normal vectors associated to the elements of the point range `points`.\n
 \b Type: a class model of `ReadablePropertyMap` with
-`PointRange::iterator` as key type and
+`PointRange::iterator::value_type` as key type and
 `geom_traits::Vector_3` as value type. \n
 No default value.
 \cgalNPEnd
@@ -67,9 +67,9 @@ Note that when a callback is run on a parallelized algorithm with `CGAL::Paralle
 \cgalNPEnd
 
 \cgalNPBegin{query_point_map} \anchor PSP_query_point_map
- is the property map containing the points associated to the iterators of the point range `queries`.\n
+ is the property map containing the points associated to the elements of the point range `queries`.\n
 \b Type: a class model of `ReadablePropertyMap` with
-`PointRange::iterator` as key type and
+`PointRange::iterator::value_type` as key type and
 `geom_traits::Point_3` as value type. \n
 <b>Default value</b>: \code CGAL::Identity_property_map<geom_traits::Point_3>\endcode
 \cgalNPEnd
@@ -146,9 +146,9 @@ multiple of a tolerance `epsilon` used to connect simplices.
 \cgalNPEnd
 
 \cgalNPBegin{plane_map} \anchor PSP_plane_map
- is the property map containing the planes associated to the iterators of the plane range `planes`.\n
+ is the property map containing the planes associated to the elements of the plane range `planes`.\n
 \b Type: a class model of `ReadablePropertyMap` with
-`PlaneRange::iterator` as key type and
+`PlaneRange::iterator::value_type` as key type and
 `geom_traits::Plane_3` as value type. \n
 <b>Default value</b>: \code CGAL::Identity_property_map<geom_traits::Plane_3>\endcode
 \cgalNPEnd

From 7656b167027a7520f455139fee4f7289b545f8e9 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 4 Oct 2018 11:52:12 +0200
Subject: [PATCH 31/35] Fix typo hightest->highest

---
 Point_set_processing_3/include/CGAL/mst_orient_normals.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Point_set_processing_3/include/CGAL/mst_orient_normals.h b/Point_set_processing_3/include/CGAL/mst_orient_normals.h
index d0f981d7824..b5da8304180 100644
--- a/Point_set_processing_3/include/CGAL/mst_orient_normals.h
+++ b/Point_set_processing_3/include/CGAL/mst_orient_normals.h
@@ -598,7 +598,7 @@ create_mst_graph(
      \cgalParamBegin{point_is_constrained_map} a model of `ReadablePropertyMap` with value type
      `bool`. Points with a `true` value will be used as seed points: their normal will be considered as already
      oriented, it won't be altered and it will be propagated to its neighbors. If this parameter is omitted, 
-     the highest point (hightest Z coordinate) will be used as the unique seed with an upward oriented
+     the highest point (highest Z coordinate) will be used as the unique seed with an upward oriented
      normal\cgalParamEnd
      \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
    \cgalNamedParamsEnd

From 52ba7cac4661d8de1f7bec6e2ac907c035665a45 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= <sebastien.loriot@cgal.org>
Date: Wed, 2 Jan 2019 10:15:53 +0100
Subject: [PATCH 32/35] remove extra typename

---
 .../Plugins/Point_set/Point_set_normal_estimation_plugin.cpp  | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp
index f0cd446b369..6ae9d2ba4c7 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp
@@ -82,7 +82,7 @@ struct Jet_estimate_normals_functor
 struct Vector_to_pmap
 {
   typedef boost::readable_property_map_tag     category;
-  typedef typename Point_set::Index            key_type;
+  typedef Point_set::Index                     key_type;
   typedef bool                                 value_type;
   typedef value_type                           reference;
   
@@ -362,7 +362,7 @@ void Polyhedron_demo_point_set_normal_estimation_plugin::on_actionNormalOrientat
     {
       std::vector<bool> constrained_map (points->size(), false);
 
-      for (typename Point_set::iterator it = points->first_selected(); it != points->end(); ++ it)
+      for (Point_set::iterator it = points->first_selected(); it != points->end(); ++ it)
         constrained_map[*it] = true;
               
       first_unoriented_point =

From 0bdce6935d34d30714837371b62fd2fb3f098bfe Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 2 Jan 2019 11:35:33 +0100
Subject: [PATCH 33/35] Remove useless variable

---
 .../Plugins/Point_set/Point_set_normal_estimation_plugin.cpp   | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp
index 6ae9d2ba4c7..5193fec725f 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp
@@ -231,9 +231,6 @@ void Polyhedron_demo_point_set_normal_estimation_plugin::on_actionNormalEstimati
     QApplication::setOverrideCursor(Qt::BusyCursor);
     QApplication::processEvents();
 
-    // First point to delete
-    Point_set::iterator first_unoriented_point = points->end();
-
     //***************************************
     // normal estimation
     //***************************************

From 294ebaf8f5155d9fdd79b0cacf5f5340bf5d87c4 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 3 Jan 2019 14:02:33 +0100
Subject: [PATCH 34/35] Remove unused variable and useless assertion

---
 Point_set_processing_3/include/CGAL/mst_orient_normals.h | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/Point_set_processing_3/include/CGAL/mst_orient_normals.h b/Point_set_processing_3/include/CGAL/mst_orient_normals.h
index b5da8304180..72cd4e3b48f 100644
--- a/Point_set_processing_3/include/CGAL/mst_orient_normals.h
+++ b/Point_set_processing_3/include/CGAL/mst_orient_normals.h
@@ -384,9 +384,8 @@ create_riemannian_graph(
     }
 
     // add source vertex (virtual, does not correspond to a point)
-    typename Riemannian_graph::vertex_descriptor v = add_vertex(riemannian_graph);
+    add_vertex(riemannian_graph);
     std::size_t source_point_index = num_input_points;
-    CGAL_point_set_processing_assertion(v == source_point_index);
 
     //
     // add edges

From 3b246294bfe908efdbac57f2feb3760dfb844b70 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Fri, 4 Jan 2019 14:39:29 +0100
Subject: [PATCH 35/35] Update CHANGES.md

---
 Installation/CHANGES.md | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md
index 86dd48044b6..9f2c4f38d50 100644
--- a/Installation/CHANGES.md
+++ b/Installation/CHANGES.md
@@ -32,6 +32,14 @@ Release date: March 2019
 -   Added the class `CGAL::Rigid_triangle_mesh_collision_detection` to detect intersections between meshes
     and volumes undergoing affine transformations.
 
+### Point Set Processing
+
+-   `CGAL::mst_orient_normals()` can now be called with a set of user-selected
+    seed points that are known to be already oriented. A new optional named
+    parameter `point_is_constrained_map` is added for this purpose. The
+    original behavior (using one unique and automatically selected seed) is
+    kept if this parameter is not used.
+
 ### 3D Fast Intersection and Distance Computation
 
 -   The primitives `AABB_face_graph_triangle_primitive` and