diff --git a/.github/workflows/cmake-all.yml b/.github/workflows/cmake-all.yml
new file mode 100644
index 00000000000..bfe53b55a0f
--- /dev/null
+++ b/.github/workflows/cmake-all.yml
@@ -0,0 +1,17 @@
+name: CMake all CGAL
+
+on: [push]
+
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+    
+    steps:
+    - uses: actions/checkout@v1
+      with:
+        fetch-depth: 0
+    - name: install dependencies
+      run: sudo apt-get install -y libboost-dev libboost-program-options-dev libmpfr-dev libeigen3-dev
+    - name: configure all
+      run: mkdir build && cd build && CXX=clang++ cmake -DWITH_examples=ON -DWITH_tests=ON -DWITH_demos=ON ..
diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Arr_naive_point_location_impl.h b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Arr_naive_point_location_impl.h
index ab194ac8361..afdd2654f16 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Arr_naive_point_location_impl.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Arr_naive_point_location_impl.h
@@ -91,7 +91,7 @@ Arr_naive_point_location<Arrangement>::locate(const Point_2& p) const
       else if (! fh->is_unbounded() && fh->number_of_outer_ccbs() > 0)
       {
         // As we have already some other containing face, one face must be
-        // contained inside the other. Two check that, we select a
+        // contained inside the other. To check that, we select a
         // representative vertex of inner_f and check whether it is contained
         // in our current face.
 
diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_active_edge.h b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_active_edge.h
index abb8fbb0942..766e5f951f1 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_active_edge.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_active_edge.h
@@ -98,6 +98,7 @@ public:
   //friend class declarations:
 
   friend class Trapezoidal_decomposition_2<Traits>;
+  friend struct internal::Non_recursive_td_map_item_destructor<Traits>;
   
 #ifdef CGAL_PM_FRIEND_CLASS
 #if defined(__SUNPRO_CC) || defined(__PGI) || defined(__INTEL_COMPILER)
@@ -114,7 +115,10 @@ public:
   friend class In_face_iterator;
 #endif
 #endif
-  
+  ~Td_active_edge(){
+    if (this->refs()==1)
+      internal::Non_recursive_td_map_item_destructor<Traits>(*this);
+  }
    /*! \class
    * Inner class Data derived from Rep class
    */
@@ -281,6 +285,10 @@ public:
   /*! Access DAG node. */
   Dag_node* dag_node() const            {return ptr()->p_node; } //m_dag_node;}
   
+  bool is_last_reference() const
+  {
+    return this->refs()==1;
+  }
   
   //@}
   
diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_active_trapezoid.h b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_active_trapezoid.h
index 4ef7826fc53..5bec4cd66ba 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_active_trapezoid.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_active_trapezoid.h
@@ -95,6 +95,7 @@ public:
   //friend class declarations:
 
   friend class Trapezoidal_decomposition_2<Traits>;
+  friend struct internal::Non_recursive_td_map_item_destructor<Traits>;
   
 #ifdef CGAL_PM_FRIEND_CLASS
 #if defined(__SUNPRO_CC) || defined(__PGI) || defined(__INTEL_COMPILER)
@@ -154,7 +155,14 @@ public:
   
   Data* ptr() const { return (Data*)(PTR); }
 	
+public:
 
+  ~Td_active_trapezoid(){
+    if (this->refs()==1)
+      internal::Non_recursive_td_map_item_destructor<Traits>(*this);
+  }
+
+private:
   
 
 #ifndef CGAL_TD_DEBUG
@@ -414,6 +422,11 @@ public:
   /*! Access DAG node. */
   Dag_node* dag_node() const            {return ptr()->p_node; }
   
+  bool is_last_reference() const
+  {
+    return this->refs()==1;
+  }
+
   void clear_neighbors()
   {
     set_lb(Td_map_item(0));
@@ -421,6 +434,11 @@ public:
     set_rb(Td_map_item(0));
     set_rt(Td_map_item(0));
   }
+
+  void non_recursive_clear_neighbors()
+  {
+    internal::Non_recursive_td_map_item_destructor<Traits>(*this);
+  }
   
   //@}
   
diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_active_vertex.h b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_active_vertex.h
index 47d4f78c5b0..086a485c848 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_active_vertex.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_active_vertex.h
@@ -129,7 +129,10 @@ public:
     Data(Vertex_const_handle _v, Halfedge_const_handle _cw_he,
          Dag_node* _p_node) :
       v(_v), cw_he(_cw_he), p_node(_p_node)
-    {}
+    {
+      CGAL_assertion( _cw_he==Halfedge_const_handle()
+                     || _cw_he->source() == v );
+    }
 
     ~Data() {}
 
@@ -176,9 +179,8 @@ public:
    */
   inline void set_cw_he(Halfedge_const_handle he)
   {
-    ptr()->cw_he = ((cw_he() != Traits::empty_he_handle()) &&
-                    (cw_he()->direction() != he->direction())) ?
-      he->twin() : he;
+    ptr()->cw_he = he->twin()->source()==ptr()->v ? he->twin() : he;
+    CGAL_assertion( ptr()->v == ptr()->cw_he->source() );
   }
 
   /*! Reset the first he going clockwise starting at 12 o'clock.
diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_dag_node.h b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_dag_node.h
index 75e20cd15c2..ac3b81a7494 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_dag_node.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Td_dag_node.h
@@ -68,8 +68,8 @@ public:
 
   //bool operator!() const {  return PTR == 0;  } //MICHAL: maybe use ptr(), and also can change to is_null or something similar
   bool is_null() const { return PTR == 0; }
-protected:
   Rep * ptr() const { return (Rep*) PTR; }
+protected:
   //Rep *& ptr() { return (Rep*) PTR; }
   void set_ptr(Rep* rep) { PTR = rep; }
 
@@ -126,7 +126,7 @@ protected:
     public:
       void operator()(Td_active_trapezoid& t) const
       {
-        t.clear_neighbors();
+        t.non_recursive_clear_neighbors();
       }
 
       template < typename Tp >
@@ -209,7 +209,49 @@ public:
   }
 
   //d'tor
-  ~Td_dag_node() { }
+  ~Td_dag_node()
+  {
+    // The following code is used to avoid recursive calls to ~Handle
+    // Node being holding elements of type Td_dag_node_handle
+    if (!this->is_null() && this->refs()==1)
+    {
+      std::vector<Td_dag_node_handle> children;
+      if (!node()->m_left_child.is_null())
+      {
+        children.push_back(node()->m_left_child);
+        node()->m_left_child.reset();
+      }
+      if (!node()->m_right_child.is_null())
+      {
+        children.push_back(node()->m_right_child);
+        node()->m_right_child.reset();
+      }
+
+      while(!children.empty())
+      {
+        Td_dag_node_handle child = children.back();
+        children.pop_back();
+        Node* child_node = (Node*) child.ptr();
+
+        if (child_node != NULL)
+        {
+          if (child.refs()==1)
+          {
+            if (!child_node->m_left_child.is_null())
+            {
+              children.push_back(child_node->m_left_child);
+              child_node->m_left_child.reset();
+            }
+            if( !child_node->m_right_child.is_null() )
+            {
+              children.push_back(child_node->m_right_child);
+              child_node->m_right_child.reset();
+            }
+          }
+        }
+      }
+    }
+  }
 
   //information retrieval
 
diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Trapezoidal_decomposition_2.h b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Trapezoidal_decomposition_2.h
index 3dc5e4c26d5..48ac5e49304 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Trapezoidal_decomposition_2.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Trapezoidal_decomposition_2.h
@@ -49,6 +49,111 @@
 
 namespace CGAL {
 
+namespace internal{
+
+// struct used to avoid recursive deletion of elements of
+// Td_map_item. Td_active_edge and Td_active_edge_item are
+// both refering to elements of the same type creating
+// recursive call to ~Handle() if we let the regular
+// calls of destructors. Here elements are copied in
+// a vector and the true deletion is done when the vector
+// is cleared.
+template <class Traits>
+struct Non_recursive_td_map_item_destructor
+{
+  typedef typename Traits::Td_map_item Td_map_item;
+  typedef typename Traits::Td_active_trapezoid Td_active_trapezoid;
+  typedef typename Traits::Td_active_edge Td_active_edge;
+
+
+  struct Child_visitor
+  {
+    typedef void result_type;
+    std::vector<Td_map_item>& m_queue;
+
+    Child_visitor(std::vector<Td_map_item>& queue)
+      : m_queue(queue)
+    {}
+
+    void operator()(Td_active_trapezoid& item)
+    {
+      if (item.is_last_reference())
+        m_queue.push_back(item);
+    }
+
+    void operator()(Td_active_edge& item)
+    {
+      if (item.is_last_reference())
+        m_queue.push_back(item);
+    }
+
+    template <class T>
+    void operator()(T&) {} // nothing to do for the other types of the variant
+  };
+
+  struct Item_visitor
+  {
+    typedef void result_type;
+    Child_visitor& m_child_visitor;
+
+    Item_visitor(Child_visitor& child_visitor)
+      : m_child_visitor(child_visitor)
+    {}
+
+    void operator()(Td_active_trapezoid& item)
+    {
+      boost::apply_visitor(m_child_visitor, item.lb());
+      boost::apply_visitor(m_child_visitor, item.lt());
+      boost::apply_visitor(m_child_visitor, item.rb());
+      boost::apply_visitor(m_child_visitor, item.rt());
+      item.clear_neighbors();
+    }
+
+    void operator()(Td_active_edge& item)
+    {
+      boost::apply_visitor(m_child_visitor, item.next());
+      item.set_next(Td_map_item(0));
+    }
+
+    template <class T>
+    void operator()(T&) {} // nothing to do for the other types of the variant
+  };
+
+  std::vector<Td_map_item> queue;
+  Child_visitor child_visitor;
+  Item_visitor item_visitor;
+
+  Non_recursive_td_map_item_destructor(Td_active_trapezoid& item)
+    : child_visitor(queue)
+    , item_visitor(child_visitor)
+  {
+    item_visitor(item);
+
+    while (!queue.empty())
+    {
+      Td_map_item item = queue.back();
+      queue.pop_back();
+      boost::apply_visitor(item_visitor, item);
+    }
+  }
+
+  Non_recursive_td_map_item_destructor(Td_active_edge& item)
+    : child_visitor(queue)
+    , item_visitor(child_visitor)
+  {
+    item_visitor(item);
+
+    while (!queue.empty())
+    {
+      Td_map_item item = queue.back();
+      queue.pop_back();
+      boost::apply_visitor(item_visitor, item);
+    }
+  }
+};
+
+} // internal
+
 /*! \class Trapezoidal_decomposition_2
  * parameters    Traits
  * Description   Implementation for a planar trapezoidal map also known as
@@ -1447,7 +1552,9 @@ public:
       {
         if (do_rebuild && not_within_limits())
         {
+#ifdef CGAL_TD_DEBUG
           std::cout << "starting over after " << number_of_curves() << std::flush;
+#endif
           start_over = true;
           clear();
           break;
@@ -1853,14 +1960,14 @@ public:
     ds->filter(representatives, Td_active_edge_item(*traits));
 
 #ifndef CGAL_TD_DEBUG
-    CGAL_warning(sz == representatives.size());
+    CGAL_warning(sz <= representatives.size());
 #else
 
     unsigned long rep = representatives.size();
-    if (sz != rep) {
+    if (sz > rep) {
       std::cerr << "\nnumber_of_curves()=" << sz;
       std::cerr << "\nrepresentatives.size()=" << rep;
-      CGAL_assertion(number_of_curves()==representatives.size());
+      CGAL_assertion(number_of_curves()<=representatives.size());
     }
 #endif
 
@@ -1873,7 +1980,15 @@ public:
       }
     }
     if (! container.empty()) {
-      CGAL::cpp98::random_shuffle(container.begin(),container.end());
+      if (sz != representatives.size())
+      {
+        std::sort(container.begin(),container.end());
+        typename Halfedge_container::iterator last = std::unique(container.begin(), container.end());
+        container.erase(last, container.end());
+      }
+      CGAL_assertion(sz==container.size());
+
+      // CGAL::cpp98::random_shuffle(container.begin(),container.end()); // already done in insert()
     }
     return sz;
   }
diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Trapezoidal_decomposition_2_impl.h b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Trapezoidal_decomposition_2_impl.h
index 715a9e74d7e..99c4d32faef 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Trapezoidal_decomposition_2_impl.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_point_location/Trapezoidal_decomposition_2_impl.h
@@ -229,7 +229,7 @@ deactivate_trapezoid(Dag_node& trpz_node, Dag_node* active_node) const
   CGAL_precondition(traits->is_td_trapezoid(trpz_node.get_data()));
   if (Td_active_trapezoid* trap =
       boost::get<Td_active_trapezoid>(&trpz_node.get_data()))
-    trap->clear_neighbors();
+    trap->non_recursive_clear_neighbors();
   trpz_node.set_data(Td_inactive_trapezoid());
   if (active_node) trpz_node.set_left_child(*active_node);
 }
@@ -510,10 +510,10 @@ update_vtx_cw_he_after_remove(Halfedge_const_handle old_he,
 
   Halfedge_const_handle cw_he(boost::apply_visitor(cw_he_visitor(), vtx_item));
   if ((old_he == cw_he) || (old_he->twin() == cw_he)) {
-    Halfedge_const_handle new_he = cw_he->twin()->prev();
-    if (new_he != old_he)
+    Halfedge_const_handle new_he = cw_he->twin()->next();
+    if (new_he != cw_he)
       boost::apply_visitor(set_cw_he_visitor(new_he), vtx_item);
-    else boost::apply_visitor(reset_cw_he_visitor(), vtx_item);
+    else boost::apply_visitor(reset_cw_he_visitor(), vtx_item); // dangling edge removed
   }
 }
 
diff --git a/CGAL_ImageIO/include/CGAL/read_vtk_image_data.h b/CGAL_ImageIO/include/CGAL/read_vtk_image_data.h
index 1b54ef99ea0..b7a058af5e2 100644
--- a/CGAL_ImageIO/include/CGAL/read_vtk_image_data.h
+++ b/CGAL_ImageIO/include/CGAL/read_vtk_image_data.h
@@ -69,9 +69,9 @@ read_vtk_image_data(vtkImageData* vtk_image, Image_3::Own owning = Image_3::OWN_
   image->ydim = dims[1];
   image->zdim = dims[2];
   image->vdim = 1;
-  image->vx = spacing[0];
-  image->vy = spacing[1];
-  image->vz = spacing[2];
+  image->vx = (spacing[0] == 0) ? 1 : spacing[0];
+  image->vy = (spacing[1] == 0) ? 1 : spacing[1];
+  image->vz = (spacing[2] == 0) ? 1 : spacing[2];
   image->tx = static_cast<float>(offset[0]);
   image->ty = static_cast<float>(offset[1]);
   image->tz = static_cast<float>(offset[2]);
diff --git a/Classification/include/CGAL/Classification/OpenCV/Random_forest_classifier.h b/Classification/include/CGAL/Classification/OpenCV/Random_forest_classifier.h
index 2b657e52cb7..2024c686972 100644
--- a/Classification/include/CGAL/Classification/OpenCV/Random_forest_classifier.h
+++ b/Classification/include/CGAL/Classification/OpenCV/Random_forest_classifier.h
@@ -25,14 +25,24 @@
 
 #include <CGAL/Classification/Feature_set.h>
 #include <CGAL/Classification/Label_set.h>
-#if (CV_MAJOR_VERSION < 3)
-#include <cv.h>
-#include <ml.h>
+
+#include <opencv2/opencv.hpp>
+
+//In opencv version 2.X the first digit is named EPOCH,
+//until version 3.0 where EPOCH disappears and it becomes MAJOR. Hence this
+//weird condition
+#ifdef CV_VERSION_EPOCH
+  #if  CV_VERSION_MAJOR == 4 && CV_VERSION_MINOR>= 11
+    #include <opencv2/ml.hpp>
+  #else
+    #include <opencv2/ml/ml.hpp>
+  #endif
 #else
-#include <opencv/cv.h>
-#include <opencv/ml.h>
+  #include <opencv2/ml.hpp>
 #endif
 
+
+
 namespace CGAL {
 
 namespace Classification {
diff --git a/Cone_spanners_2/doc/Cone_spanners_2/PackageDescription.txt b/Cone_spanners_2/doc/Cone_spanners_2/PackageDescription.txt
index 6ffa6f32a52..9dfe45083f8 100644
--- a/Cone_spanners_2/doc/Cone_spanners_2/PackageDescription.txt
+++ b/Cone_spanners_2/doc/Cone_spanners_2/PackageDescription.txt
@@ -7,7 +7,7 @@
 \cgalPkgPicture{Logo-ConeSpanners.png}
 
 \cgalPkgSummaryBegin
-\cgalPkgAuthors{Weisheng Si, Quincy Tse and Frédérk Paradis}
+\cgalPkgAuthors{Weisheng Si, Quincy Tse and Frédérik Paradis}
 \cgalPkgDesc{This package provides functors for constructing two kinds of cone-based spanners:
 Yao graph and Theta graph, given a set of vertices on the plane and the directions of cone boundaries.
 Both exact and inexact constructions are supported.
diff --git a/Documentation/doc/Documentation/Installation.txt b/Documentation/doc/Documentation/Installation.txt
index 8dbe782d0db..541befc0636 100644
--- a/Documentation/doc/Documentation/Installation.txt
+++ b/Documentation/doc/Documentation/Installation.txt
@@ -140,6 +140,8 @@ supporting <a href="https://isocpp.org/wiki/faq/cpp14">C++14</a> or later.
 
 It may work for older versions of the above listed compilers.
 
+\attention Recent versions of CMake are needed for recent versions of MS Visual C++. Please refer to CMake's documentation for further information.
+
 
 \section secconfigwithcmake Configuring CGAL with CMake 
 
diff --git a/Documentation/doc/resources/1.8.13/menu_version.js b/Documentation/doc/resources/1.8.13/menu_version.js
index fd75226e1e9..1f8a8eef476 100644
--- a/Documentation/doc/resources/1.8.13/menu_version.js
+++ b/Documentation/doc/resources/1.8.13/menu_version.js
@@ -3,22 +3,23 @@
 
   var url_re =  /(cgal\.geometryfactory\.com\/CGAL\/doc\/|doc\.cgal\.org\/)(master|latest|(\d\.\d+|\d\.\d+\.\d+))\//;
   var url_local =  /.*\/doc_output\//;
-  var current_version_local = '4.14-beta1'
+  var current_version_local = '5.0-beta2'
   var all_versions = [
-    'master',
-    'latest',
-    '4.14',
-    '4.13.1',
-    '4.12.2',
-    '4.11.3',
-    '4.10.2',
-    '4.9.1',
-    '4.8.2',
-    '4.7',
-    '4.6.3',
-    '4.5.2',
-    '4.4',
-    '4.3'
+      'master',
+      'latest',
+      '5.0',
+      '4.14.1',
+      '4.13.2',
+      '4.12.2',
+      '4.11.3',
+      '4.10.2',
+      '4.9.1',
+      '4.8.2',
+      '4.7',
+      '4.6.3',
+      '4.5.2',
+      '4.4',
+      '4.3'
   ];
 
   function build_select(current_version) {
diff --git a/Documentation/doc/resources/1.8.13/stylesheet.css b/Documentation/doc/resources/1.8.13/stylesheet.css
index 4f1ab9195b4..41df79908ef 100644
--- a/Documentation/doc/resources/1.8.13/stylesheet.css
+++ b/Documentation/doc/resources/1.8.13/stylesheet.css
@@ -498,7 +498,7 @@ table.memberdecls {
         white-space: nowrap;
 }
 
-.memItemRight {
+.memItemRight, .memTemplItemRight {
 	width: 100%;
 }
 
diff --git a/Documentation/doc/resources/1.8.14/menu_version.js b/Documentation/doc/resources/1.8.14/menu_version.js
index fd75226e1e9..1f8a8eef476 100644
--- a/Documentation/doc/resources/1.8.14/menu_version.js
+++ b/Documentation/doc/resources/1.8.14/menu_version.js
@@ -3,22 +3,23 @@
 
   var url_re =  /(cgal\.geometryfactory\.com\/CGAL\/doc\/|doc\.cgal\.org\/)(master|latest|(\d\.\d+|\d\.\d+\.\d+))\//;
   var url_local =  /.*\/doc_output\//;
-  var current_version_local = '4.14-beta1'
+  var current_version_local = '5.0-beta2'
   var all_versions = [
-    'master',
-    'latest',
-    '4.14',
-    '4.13.1',
-    '4.12.2',
-    '4.11.3',
-    '4.10.2',
-    '4.9.1',
-    '4.8.2',
-    '4.7',
-    '4.6.3',
-    '4.5.2',
-    '4.4',
-    '4.3'
+      'master',
+      'latest',
+      '5.0',
+      '4.14.1',
+      '4.13.2',
+      '4.12.2',
+      '4.11.3',
+      '4.10.2',
+      '4.9.1',
+      '4.8.2',
+      '4.7',
+      '4.6.3',
+      '4.5.2',
+      '4.4',
+      '4.3'
   ];
 
   function build_select(current_version) {
diff --git a/Documentation/doc/resources/1.8.14/stylesheet.css b/Documentation/doc/resources/1.8.14/stylesheet.css
index 4f1ab9195b4..41df79908ef 100644
--- a/Documentation/doc/resources/1.8.14/stylesheet.css
+++ b/Documentation/doc/resources/1.8.14/stylesheet.css
@@ -498,7 +498,7 @@ table.memberdecls {
         white-space: nowrap;
 }
 
-.memItemRight {
+.memItemRight, .memTemplItemRight {
 	width: 100%;
 }
 
diff --git a/Documentation/doc/resources/1.8.4/menu_version.js b/Documentation/doc/resources/1.8.4/menu_version.js
index fd75226e1e9..1f8a8eef476 100644
--- a/Documentation/doc/resources/1.8.4/menu_version.js
+++ b/Documentation/doc/resources/1.8.4/menu_version.js
@@ -3,22 +3,23 @@
 
   var url_re =  /(cgal\.geometryfactory\.com\/CGAL\/doc\/|doc\.cgal\.org\/)(master|latest|(\d\.\d+|\d\.\d+\.\d+))\//;
   var url_local =  /.*\/doc_output\//;
-  var current_version_local = '4.14-beta1'
+  var current_version_local = '5.0-beta2'
   var all_versions = [
-    'master',
-    'latest',
-    '4.14',
-    '4.13.1',
-    '4.12.2',
-    '4.11.3',
-    '4.10.2',
-    '4.9.1',
-    '4.8.2',
-    '4.7',
-    '4.6.3',
-    '4.5.2',
-    '4.4',
-    '4.3'
+      'master',
+      'latest',
+      '5.0',
+      '4.14.1',
+      '4.13.2',
+      '4.12.2',
+      '4.11.3',
+      '4.10.2',
+      '4.9.1',
+      '4.8.2',
+      '4.7',
+      '4.6.3',
+      '4.5.2',
+      '4.4',
+      '4.3'
   ];
 
   function build_select(current_version) {
diff --git a/Documentation/doc/resources/1.8.4/stylesheet.css b/Documentation/doc/resources/1.8.4/stylesheet.css
index dc14a84c031..688b81b2557 100644
--- a/Documentation/doc/resources/1.8.4/stylesheet.css
+++ b/Documentation/doc/resources/1.8.4/stylesheet.css
@@ -466,7 +466,7 @@ table.memberdecls {
         white-space: nowrap;
 }
 
-.memItemRight {
+.memItemRight, .memTemplItemRight {
 	width: 100%;
 }
 
diff --git a/Geomview/package_info/Geomview/dependencies b/Geomview/package_info/Geomview/dependencies
index cbbfc58b29b..7c8f7ab1d75 100644
--- a/Geomview/package_info/Geomview/dependencies
+++ b/Geomview/package_info/Geomview/dependencies
@@ -1,5 +1,6 @@
 Algebraic_foundations
 Circulator
+Geomview
 Installation
 Interval_support
 Kernel_23
diff --git a/GraphicsView/include/CGAL/Buffer_for_vao.h b/GraphicsView/include/CGAL/Buffer_for_vao.h
index beeb5907486..3fc91ec4f10 100644
--- a/GraphicsView/include/CGAL/Buffer_for_vao.h
+++ b/GraphicsView/include/CGAL/Buffer_for_vao.h
@@ -39,9 +39,6 @@
 
 namespace CGAL
 {
-  typedef CGAL::Exact_predicates_inexact_constructions_kernel Local_kernel;
-  typedef Local_kernel::Point_3  Local_point;
-  typedef Local_kernel::Vector_3 Local_vector;
 
 //------------------------------------------------------------------------------
 namespace internal
@@ -58,90 +55,77 @@ namespace internal
                n.z()+((p.x()-q.x())*(p.y()+q.y())));
     }
 
-  inline
-  Local_vector compute_normal_of_face(const std::vector<Local_point>& points)
+  template <class Point, class Vector>
+  Vector compute_normal_of_face(const std::vector<Point>& points)
   {
-    Local_vector normal(CGAL::NULL_VECTOR);
+    Vector normal(CGAL::NULL_VECTOR);
     unsigned int nb = 0;
     for (std::size_t i=0; i<points.size(); ++i)
     {
-      newell_single_step_3(points[i], points[(i+1)%points.size()], normal);
+      internal::newell_single_step_3(points[i], points[(i+1)%points.size()],
+                                     normal);
       ++nb;
     }
     
     assert(nb>0);
-    return (Local_kernel::Construct_scaled_vector_3()(normal, 1.0/nb));
+    return (typename Kernel_traits<Vector>::Kernel::Construct_scaled_vector_3()
+            (normal, 1.0/nb));
   }
 
   ////////////////////////////////////////////////////////////////
   // Structs to transform any CGAL point/vector into a Local_point/Local_vector
-  template<typename K>
+  template<typename K, typename Local_kernel>
   struct Geom_utils
   {
-    static Local_point get_local_point(const typename K::Point_2& p)
+    static typename Local_kernel::Point_3 get_local_point(const typename K::Point_2& p)
     {
       CGAL::Cartesian_converter<K, Local_kernel> converter;
       return Local_point(converter(p.x()), 0, converter(p.y()));
     }
-    static Local_point get_local_point(const typename K::Weighted_point_2& p)
+    static typename Local_kernel::Point_3 get_local_point(const typename K::Weighted_point_2& p)
     {
       typename K::Point_2 lp(p);
-      return Geom_utils<K>::get_local_point(lp);
+      return Geom_utils<K, Local_kernel>::get_local_point(lp);
     }
-    static Local_point get_local_point(const typename K::Point_3& p)
+    static typename Local_kernel::Point_3 get_local_point(const typename K::Point_3& p)
     {
       CGAL::Cartesian_converter<K, Local_kernel> converter;
       return converter(p);
     }
-    static Local_point get_local_point(const typename K::Weighted_point_3& p)
+    static typename Local_kernel::Point_3 get_local_point(const typename K::Weighted_point_3& p)
     {
       typename K::Point_3 lp(p);
-      return Geom_utils<K>::get_local_point(lp);
+      return Geom_utils<K, Local_kernel>::get_local_point(lp);
     }
-    static Local_vector get_local_vector(const typename K::Vector_2& v)
+    static typename Local_kernel::Vector_3 get_local_vector(const typename K::Vector_2& v)
     {
       CGAL::Cartesian_converter<K, Local_kernel> converter;
       return Local_vector(converter(v.x()), 0, converter(v.y()));
     }
-    static Local_vector get_local_vector(const typename K::Vector_3& v)
+    static typename Local_kernel::Vector_3 get_local_vector(const typename K::Vector_3& v)
     {
       CGAL::Cartesian_converter<K, Local_kernel> converter;
       return converter(v);
     }
   };
 
-  // Specialization for Local_kernel, because there is no need of convertion here.
-  template<>
-  struct Geom_utils<Local_kernel>
+  // Specialization when K==Local_kernel, because there is no need of convertion here.
+  template<typename Local_kernel>
+  struct Geom_utils<Local_kernel, Local_kernel>
   {
-    static Local_point get_local_point(const Local_kernel::Point_2& p)
-    { return Local_point(p.x(), 0, p.y()); }
-    static Local_point get_local_point(const Local_kernel::Weighted_point_2& p)
-    { return Local_point(p.point().x(), 0, p.point().y());}
-    static const Local_point & get_local_point(const Local_kernel::Point_3& p)
+    static typename Local_kernel::Point_3 get_local_point(const typename Local_kernel::Point_2& p)
+    { return typename Local_kernel::Point_3(p.x(), 0, p.y()); }
+    static typename Local_kernel::Point_3 get_local_point(const typename Local_kernel::Weighted_point_2& p)
+    { return typename Local_kernel::Point_3(p.point().x(), 0, p.point().y());}
+    static const typename Local_kernel::Point_3 & get_local_point(const typename Local_kernel::Point_3& p)
     { return p; }
-    static Local_point get_local_point(const Local_kernel::Weighted_point_3& p)
-    { return Local_point(p);}
-    static Local_vector get_local_vector(const Local_kernel::Vector_2& v)
-    { return Local_vector(v.x(), 0, v.y()); }
-    static const Local_vector& get_local_vector(const Local_kernel::Vector_3& v)
+    static typename Local_kernel::Point_3 get_local_point(const typename Local_kernel::Weighted_point_3& p)
+    { return typename Local_kernel::Point_3(p);}
+    static typename Local_kernel::Vector_3 get_local_vector(const typename Local_kernel::Vector_2& v)
+    { return typename Local_kernel::Vector_3(v.x(), 0, v.y()); }
+    static const typename Local_kernel::Vector_3& get_local_vector(const typename Local_kernel::Vector_3& v)
     { return v; }
-  };    
-
-  ////////////////////////////////////////////////////////////////
-  // Global function to simplify function calls.
-  template<typename KPoint>
-  Local_point get_local_point(const KPoint& p)
-  {
-    return Geom_utils<typename CGAL::Kernel_traits<KPoint>::Kernel>::
-      get_local_point(p);
-  }
-  template<typename KVector>
-  Local_vector get_local_vector(const KVector& v)
-  {
-    return Geom_utils<typename CGAL::Kernel_traits<KVector>::Kernel>::
-      get_local_vector(v);
-  }
+  };
 } // End namespace internal
 
 //------------------------------------------------------------------------------
@@ -149,6 +133,10 @@ template<typename BufferType=float, typename IndexType=std::size_t>
 class Buffer_for_vao
 {
 public:
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Local_kernel;
+  typedef Local_kernel::Point_3  Local_point;
+  typedef Local_kernel::Vector_3 Local_vector;
+
   Buffer_for_vao(std::vector<BufferType>* pos=nullptr,
                  std::vector<IndexType>* indices=nullptr,
                  CGAL::Bbox_3* bbox=nullptr,
@@ -220,7 +208,7 @@ public:
   {
     if (!has_position()) return (std::size_t)-1;
     
-    Local_point p=internal::get_local_point(kp);
+    Local_point p=get_local_point(kp);
     add_point_in_buffer(p, *m_pos_buffer);
 
     if (m_bb!=nullptr)
@@ -293,7 +281,7 @@ public:
   template<typename KNormal>
   void face_begin(const KNormal& kv)
   {
-    m_normal_of_face=internal::get_local_vector(kv);
+    m_normal_of_face=get_local_vector(kv);
     face_begin_internal(false, true);
   }
 
@@ -302,7 +290,7 @@ public:
   void face_begin(const CGAL::Color& c, const KNormal& kv)
   {
     m_color_of_face=c;
-    m_normal_of_face=internal::get_local_vector(kv);
+    m_normal_of_face=get_local_vector(kv);
     face_begin_internal(true, true);
   }
 
@@ -314,7 +302,7 @@ public:
   {
     if (!is_a_face_started()) return false;
 
-    Local_point p=internal::get_local_point(kp);
+    Local_point p=get_local_point(kp);
     if (m_points_of_face.empty() || m_points_of_face.back()!=p) // TODO test if the distance between prev point and kp is smaller than an epsilon (?? not sure ??)
     {
       m_points_of_face.push_back(p);
@@ -331,7 +319,7 @@ public:
   {
     if (add_point_in_face(kp))
     {
-      m_vertex_normals_for_face.push_back(internal::get_local_vector(p_normal));
+      m_vertex_normals_for_face.push_back(get_local_vector(p_normal));
       return true;
     }
     return false;
@@ -388,7 +376,8 @@ public:
     }
     
     Local_vector normal=(m_started_face_has_normal?m_normal_of_face:
-                         internal::compute_normal_of_face(m_points_of_face));
+                         internal::compute_normal_of_face
+                         <Local_point, Local_vector>(m_points_of_face));
 
     if (m_points_of_face.size()==3)
     { triangular_face_end_internal(normal); } // Triangle: no need to triangulate
@@ -411,7 +400,7 @@ public:
   template<typename KPoint>
   static void add_point_in_buffer(const KPoint& kp, std::vector<float>& buffer)
   {
-    Local_point p=internal::get_local_point(kp);
+    Local_point p=get_local_point(kp);
     buffer.push_back(p.x());
     buffer.push_back(p.y());
     buffer.push_back(p.z());
@@ -421,7 +410,7 @@ public:
   template<typename KVector>
   static void add_normal_in_buffer(const KVector& kv, std::vector<float>& buffer)
   {
-    Local_vector n=internal::get_local_vector(kv);
+    Local_vector n=get_local_vector(kv);
     buffer.push_back(n.x());
     buffer.push_back(n.y());
     buffer.push_back(n.z());
@@ -796,6 +785,21 @@ protected:
     { add_normal_in_buffer(kv, *m_gouraud_normal_buffer); }
   }
 
+protected:
+  // Shortcuts to simplify function calls.
+  template<typename KPoint>
+  static Local_point get_local_point(const KPoint& p)
+  {
+    return internal::Geom_utils<typename CGAL::Kernel_traits<KPoint>::Kernel, Local_kernel>::
+      get_local_point(p);
+  }
+  template<typename KVector>
+  static Local_vector get_local_vector(const KVector& v)
+  {
+    return internal::Geom_utils<typename CGAL::Kernel_traits<KVector>::Kernel, Local_kernel>::
+      get_local_vector(v);
+  }
+
 protected:
   // Types usefull for triangulation
   struct Vertex_info
diff --git a/GraphicsView/include/CGAL/Qt/Basic_viewer_qt.h b/GraphicsView/include/CGAL/Qt/Basic_viewer_qt.h
index 8b23594cd7c..29e615e26c1 100644
--- a/GraphicsView/include/CGAL/Qt/Basic_viewer_qt.h
+++ b/GraphicsView/include/CGAL/Qt/Basic_viewer_qt.h
@@ -446,6 +446,20 @@ public:
   }
 
 protected:
+  // Shortcuts to simplify function calls.
+  template<typename KPoint>
+  static Local_point get_local_point(const KPoint& p)
+  {
+    return internal::Geom_utils<typename CGAL::Kernel_traits<KPoint>::Kernel, Local_kernel>::
+      get_local_point(p);
+  }
+  template<typename KVector>
+  static Local_vector get_local_vector(const KVector& v)
+  {
+    return internal::Geom_utils<typename CGAL::Kernel_traits<KVector>::Kernel, Local_kernel>::
+      get_local_vector(v);
+  }
+
   void compile_shaders()
   {
     rendering_program_face.removeAllShaders();
diff --git a/GraphicsView/include/CGAL/Qt/GraphicsItem.h b/GraphicsView/include/CGAL/Qt/GraphicsItem.h
index 7c24add2bef..939831ce432 100644
--- a/GraphicsView/include/CGAL/Qt/GraphicsItem.h
+++ b/GraphicsView/include/CGAL/Qt/GraphicsItem.h
@@ -43,9 +43,7 @@ namespace Qt {
 class CGAL_QT_EXPORT GraphicsItem : public QObject, public QGraphicsItem {
 
   Q_OBJECT
-#ifdef CGAL_HEADER_ONLY
   Q_INTERFACES(QGraphicsItem)
-#endif
 
 public Q_SLOTS:
 
diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md
index b0d3bf73941..51e45d1c998 100644
--- a/Installation/CHANGES.md
+++ b/Installation/CHANGES.md
@@ -1,144 +1,180 @@
 Release History
 ===============
 
-Release 5.0
+[Release 5.0](https://github.com/CGAL/cgal/releases/tag/releases%2FCGAL-5.0)
 -----------
 
-Release date: September 2019
+Release date: October 2019
 
-### Polygonal Surface Reconstruction (new package)
+### General changes
 
--   This package provides a method for piecewise planar object reconstruction from point clouds. 
-    The method takes as input an unordered point set sampled from a piecewise planar object 
-    and outputs a compact and watertight surface mesh interpolating the input point set. 
-    The method assumes that all necessary major planes are provided (or can be extracted from 
-    the input point set using the shape detection method described in Point Set Shape Detection, 
-    or any other alternative methods).The method can handle arbitrary piecewise planar objects 
-    and is capable of recovering sharp features and is robust to noise and outliers.
-    
-### Solver Interface
+- CGAL 5.0 is the first release of CGAL that requires a C++ compiler
+  with the support of C++14 or later. The new list of supported
+  compilers is:
+  - Visual C++ 14.0 (from Visual Studio 2015 Update 3) or later,
+  - Gnu g++ 6.3 or later (on Linux or MacOS),
+  - LLVM Clang version 8.0 or later (on Linux or MacOS), and
+  - Apple Clang compiler versions 7.0.2 and 10.0.1 (on MacOS).
+- Since CGAL 4.9, CGAL can be used as a header-only library, with
+  dependencies. Since CGAL 5.0, that is now the default, unless
+  specified differently in the (optional) CMake configuration.
 
- - Added concepts and models for solving Mixed Integer Programming (MIP) problems with or without constraints. 
-      
-## 2D Triangulations
--   Added range types and functions that return ranges, for example
-    for all vertices, which enables to use C++11 for-loops.
--   **Breaking change**: Removed the functions `CGAL::Constrained_triangulation_plus_2::
-    vertices_in_constraint_{begin/end}(Vertex_handle va, Vertex_handle vb) const;`,
-    and `CGAL::Constrained_triangulation_plus_2::remove_constraint((Vertex_handle va, Vertex_handle vb)`,
-    that is a pair of vertex handles is no longer a key for a polyline constraint.
-    Users must use a version prior to 5.0 if they need this functionality.
--   **Breaking change**: Removed the deprecated classes `CGAL::Regular_triangulation_euclidean_traits_2`, `CGAL::Regular_triangulation_filtered_traits_2`.   Users must use a version prior to 5.0 if they need these classes.
--   **Breaking change**: The constructor and the `insert()` function of `CGAL::Triangulation_2` which takes
-    a range of points as argument no longer performs a `spatial_sort()` of the points.
--   Add constructor and `insert()` function to `CGAL::Triangulation_2` that takes a range of points with info.  
- - **Breaking change**: The graph traits enabling CGAL's 2D triangulations to be used as a parameter
-    for any graph-based algorithm of CGAL (or boost) have been improved to fully model the `FaceGraph` concept.
-    In addition, only the finite simplicies (those not incident to the infinite vertex) of the 2D triangulations
-    are now visibile through this scope. The complete triangulation can still be accessed as a graph,
-    by using the graph traits of the underlying triangulation data structure (usually,
-    `CGAL::Triangulation_data_structure_2`).
- -  Introduced a new face base class, `Triangulation_face_base_with_id_2` which enables storing
-    user-defined integer IDs in the face of any 2D triangulation, a precondition to use some
-    BGL algorithms.
-### 3D Triangulations
--   Added range types and functions that return ranges, for example
-    for all vertices, which enables to use C++11 for-loops.
--   **Breaking change**: The constructor and the `insert()` function of `CGAL::Triangulation_3` which takes
-    a range of points as argument no longer performs a `spatial_sort()` of the points.
--   Add constructor and `insert()` function to `CGAL::Triangulation_3` that takes a range of points with info.  
 
-### Surface Mesh
- -   New functions to read and write using the PLY format,
-     `CGAL::read_ply()` and `CGAL::write_ply()`, allowing to save and
-     load additional property maps of the surface mesh.
+### [Polygonal Surface Reconstruction](https://doc.cgal.org/5.0/Manual/packages.html#PkgPolygonalSurfaceReconstruction) (new package)
 
-### 3D Point Set
- -   The PLY IO functions now take an additional optional parameter to
-     read/write comments from/in the PLY header.
+ -   This package provides a method for piecewise planar object reconstruction from point clouds. 
+     The method takes as input an unordered point set sampled from a piecewise planar object 
+     and outputs a compact and watertight surface mesh interpolating the input point set. 
+     The method assumes that all necessary major planes are provided (or can be extracted from 
+     the input point set using the shape detection method described in Point Set Shape Detection, 
+     or any other alternative methods).The method can handle arbitrary piecewise planar objects 
+     and is capable of recovering sharp features and is robust to noise and outliers. See also 
+     the associated [blog entry](https://www.cgal.org/2019/08/05/Polygonal_surface_reconstruction/).
 
-### Point Set Processing
- -   **Breaking change**: the old API using iterators and overloads
-     for optional parameters is now removed (it was deprecated since
-     CGAL 4.12). The current (and now only) API uses ranges and Named
-     Parameters.
- -   Added the possibility to use the named parameter
-     `neighbor_radius` to use spherical neighbor queries instead of
-     K-nearest neighbors queries for the following functions:
-     `CGAL::bilateral_smooth_point_set()`,
-     `CGAL::jet_estimate_normals()`, `CGAL::jet_smooth_point_set()`,
-     `CGAL::mst_orient_normals()`, `CGAL::pca_estimate_normals()` and
-     `CGAL::remove_outliers()`.
+### [Shape Detection](https://doc.cgal.org/5.0/Manual/packages.html#PkgShapeDetection) (major changes)
+ -   **Breaking change:** The concept `ShapeDetectionTraits` has been renamed to [`EfficientRANSACTraits`](https://doc.cgal.org/5.0/Shape_detection/classEfficientRANSACTraits.html).
+ -   **Breaking change:** The `Shape_detection_3` namespace has been renamed to [`Shape_detection`](https://doc.cgal.org/5.0/Shape_detection/annotated.html).
+ -   Added a new, generic implementation of region growing. This enables for example applying region growing to inputs such as 2D and 3D point sets,
+     or models of the [`FaceGraph`](https://doc.cgal.org/5.0/BGL/classFaceGraph.html) concept. Learn more about this new algorithm with this [blog entry](https://www.cgal.org/2019/07/30/Shape_detection/).
 
-### Polygon Mesh Processing
--    Added the function `CGAL::Polygon_mesh_processing::non_manifold_vertices()`,
-     which can be used to collect all the non-manifold vertices (i.e. pinched vertices,
-     or vertices appearing in multiple umbrellas) of a mesh.
- -   Introduced a wide range of new functions related to location of queries on a triangle mesh,
-     such as `CGAL::Polygon_mesh_processing::locate(Point, Mesh)`, . The location of a point on a triangle mesh
-     is expressed as the pair of a face and the barycentric coordinates of the point in this face,
-     enabling robust manipulation of locations (for example, intersections of two 3D segments living
-     within the same face).
- -   Added the function `CGAL::Polygon_mesh_processing::centroid()`, which computes
-     the centroid of a closed triangle mesh.
- -   Added the functions `CGAL::Polygon_mesh_processing::stitch_boundary_cycle()` and
-     `CGAL::Polygon_mesh_processing::stitch_boundary_cycles()`, which can be used
-     to try and merge together geometrically compatible but combinatorially different halfedges
+### [dD Geometry Kernel](https://doc.cgal.org/5.0/Manual/packages.html#PkgKernelD)
+ -   A new exact kernel, [`Epeck_d`](https://doc.cgal.org/5.0/Kernel_d/structCGAL_1_1Epeck__d.html), is now available.
+
+### [2D and 3D Linear Geometry Kernel](https://doc.cgal.org/5.0/Manual/packages.html#PkgKernel23)
+ -   Added a new concept, [`ComputeApproximateAngle_3`](https://doc.cgal.org/5.0/Kernel_23/classKernel_1_1ComputeApproximateAngle__3.html),
+     to the 3D Kernel concepts to compute the approximate angle between two 3D vectors. Corresponding functors
+     in the model ([`Compute_approximate_angle_3`](https://doc.cgal.org/5.0/Kernel_23/classKernel.html#a183c9ac358a4ccddc04e680f8ed16c0b)) 
+     and free function ([`approximate_angle`](https://doc.cgal.org/5.0/Kernel_23/group__approximate__angle__grp.html))
+     have also been added.
+ -   The following objects are now hashable and thus trivially usable 
+     with [`std::unordered_set`](https://en.cppreference.com/w/cpp/container/unordered_set) 
+     and [`std::unordered_map`](https://en.cppreference.com/w/cpp/header/unordered_map):
+     `CGAL::Aff_transformation_2`, `CGAL::Aff_transformation_3`,
+     `CGAL::Bbox_2`, `CGAL::Bbox_3`, `CGAL::Circle_2`,
+     `CGAL::Iso_cuboid_3`, `CGAL::Iso_rectangle_2`, `CGAL::Point_2`,
+     `CGAL::Point_3`, `CGAL::Segment_2`, `CGAL::Segment_3`,
+     `CGAL::Sphere_3`, `CGAL::Vector_2`, `CGAL::Vector_3`,
+     `CGAL::Weighted_point_2` and `CGAL::Weighted_point_3`.
+
+### [Polygon Mesh Processing](https://doc.cgal.org/latest/Manual/packages.html#PkgPolygonMeshProcessing)
+ -   Introduced a [wide range of new functions](https://doc.cgal.org/5.0/Polygon_mesh_processing/index.html#title36) 
+     related to location of queries on a triangle mesh,
+     such as [`CGAL::Polygon_mesh_processing::locate(Point, Mesh)`](https://doc.cgal.org/5.0/Polygon_mesh_processing/group__PMP__locate__grp.html#gada09bd8740ba69ead9deca597d53cf15).
+     The location of a point on a triangle mesh is expressed as the pair of a face and the barycentric
+     coordinates of the point in this face, enabling robust manipulation of locations 
+     (for example, intersections of two 3D segments living within the same face).
+ -   Added the mesh smoothing function [`smooth_mesh()`](https://doc.cgal.org/5.0/Polygon_mesh_processing/group__PMP__meshing__grp.html#gaa0551d546f6ab2cd9402bea12d8332a3),
+     which can be used to improve the quality of triangle elements based on various geometric characteristics.
+ -   Added the shape smoothing function [`smooth_shape()`](https://doc.cgal.org/5.0/Polygon_mesh_processing/group__PMP__meshing__grp.html#gaaa083ec78bcecf351e04d1bbf460b4a2),
+     which can be used to smooth the surface of a triangle mesh, using the mean curvature flow to perform noise removal.
+     (See also the new entry in the [User Manual](https://doc.cgal.org/5.0/Polygon_mesh_processing/index.html#title8))
+ -   Added the function [`CGAL::Polygon_mesh_processing::centroid()`](https://doc.cgal.org/5.0/Polygon_mesh_processing/group__measure__grp.html#ga6da5119ce2c50729fda11a90ae7fb9ba),
+     which computes the centroid of a closed triangle mesh.
+ -   Added the functions [`CGAL::Polygon_mesh_processing::stitch_boundary_cycle()`](https://doc.cgal.org/5.0/Polygon_mesh_processing/group__PMP__repairing__grp.html#ga9c12c4878c08a117b3733bb45f1a34cf)
+     and [`CGAL::Polygon_mesh_processing::stitch_boundary_cycles()`](https://doc.cgal.org/5.0/Polygon_mesh_processing/group__PMP__repairing__grp.html#ga24d5ae37f62064b3fc576ba48a4ccc63),
+     which can be used to try and merge together geometrically compatible but combinatorially different halfedges
      that belong to the same boundary cycle.
--    It is now possible to pass a face-size property map to `CGAL::Polygon_mesh_processing::keep_large_connected_components()`
-     and `CGAL::Polygon_mesh_processing::keep_largest_connected_components()`, enabling users to define
+ -   It is now possible to pass a face-size property map to [`CGAL::Polygon_mesh_processing::keep_large_connected_components()`](https://doc.cgal.org/5.0/Polygon_mesh_processing/group__keep__connected__components__grp.html#ga48e7b3e6922ee78cf8ce801e3e325d9a)
+     and [`CGAL::Polygon_mesh_processing::keep_largest_connected_components()`](https://doc.cgal.org/5.0/Polygon_mesh_processing/group__keep__connected__components__grp.html#ga68c6c29dfc6a26a6a2f8befe6944f19d), enabling users to define
      how the size of a face is computed (the size of the connected component is the sum of the sizes of its faces).
      If no property map is passed, the behavior is unchanged to previous versions: the size
      of a connected component is the number of faces it contains.
- -   Added the mesh smoothing function `smooth_mesh()`, which can be used to
-     improve the quality of triangle elements based on various geometric characteristics.
- -   Added the shape smoothing function `smooth_shape()`, which can be used to
-     smooth the surface of a triangle mesh, using the mean curvature flow to perform noise removal.
+ -   Added the function [`CGAL::Polygon_mesh_processing::non_manifold_vertices()`](https://doc.cgal.org/5.0/Polygon_mesh_processing/group__PMP__repairing__grp.html#ga36098d2415efd0604b7b996163bc22db),
+     which can be used to collect all the non-manifold vertices (i.e. pinched vertices,
+     or vertices appearing in multiple umbrellas) of a mesh.
 
-### IO Streams
- -   **Breaking change:** The API of `CGAL::Color` has been cleaned up.
+### [3D Point Set](https://doc.cgal.org/5.0/Manual/packages.html#PkgPointSet3)
+ -   The [PLY IO functions](https://doc.cgal.org/5.0/Point_set_3/group__PkgPointSet3IO.html) now take an additional optional parameter to 
+     read/write comments from/in the PLY header.
 
-### Shape Detection
- -   Added a new generic implementation of region growing.
- -   New region growing can be launched on points in 2D and 3D and on a face graph.
- -   **Breaking change:** ShapeDetectionTraits is renamed to EfficientRANSACTraits.
- -   **Breaking change:** Shape_detection_3 namespace is renamed to Shape_detection.
+### [Point Set Processing](https://doc.cgal.org/latest/Manual/packages.html#PkgPointSetProcessing3)
+ -   **Breaking change**: the API using iterators and overloads for optional parameters (deprecated since
+     CGAL 4.12) has been removed. The current (and now only) API uses ranges and Named Parameters.
+ -   Added the possibility to use the named parameter
+     [`neighbor_radius`](https://doc.cgal.org/5.0/Point_set_processing_3/group__psp__namedparameters.html#PSP_neighbor_radius) 
+     to use spherical neighbor queries instead of K-nearest neighbors queries for the following functions:
+     [`CGAL::bilateral_smooth_point_set()`](https://doc.cgal.org/5.0/Point_set_processing_3/group__PkgPointSetProcessing3Algorithms.html#ga4f82723e2f0bb33f3677e29e0208a256),
+     [`CGAL::jet_estimate_normals()`](https://doc.cgal.org/5.0/Point_set_processing_3/group__PkgPointSetProcessing3Algorithms.html#ga0cd0f87de690d4edf82740e856efa491),
+     [`CGAL::jet_smooth_point_set()`](https://doc.cgal.org/5.0/Point_set_processing_3/group__PkgPointSetProcessing3Algorithms.html#ga549402c0a8a8b6b71875181e93961521),
+     [`CGAL::mst_orient_normals()`](https://doc.cgal.org/5.0/Point_set_processing_3/group__PkgPointSetProcessing3Algorithms.html#ga50c98d5c5ae5535bce6f32eddbd03f33),
+     [`CGAL::pca_estimate_normals()`](https://doc.cgal.org/5.0/Point_set_processing_3/group__PkgPointSetProcessing3Algorithms.html#ga8c642da96a025ab32445aeb6cc219b0b) and
+     [`CGAL::remove_outliers()`](https://doc.cgal.org/5.0/Point_set_processing_3/group__PkgPointSetProcessing3Algorithms.html#gafd0b5a21ec5042e4bca09cb43f1847f9).
 
-### 3D Boolean Operations on Nef Polyhedra
- -   Added a function to convert a Nef_polyhedron_3 to a polygon soup: `CGAL::convert_nef_polyhedron_to_polygon_soup()`
+### [2D Triangulations](https://doc.cgal.org/5.0/Manual/packages.html#PkgTriangulation2)
+ -   **Breaking change**: Removed the deprecated functions `CGAL::Constrained_triangulation_plus_2::
+     vertices_in_constraint_{begin/end}(Vertex_handle va, Vertex_handle vb) const;`,
+     and `CGAL::Constrained_triangulation_plus_2::remove_constraint(Vertex_handle va, Vertex_handle vb)`,
+     that is a pair of vertex handles is no longer a key for a polyline constraint.
+     Users must use a version prior to 5.0 if they need this functionality.
+ -   **Breaking change**: Removed the deprecated classes `CGAL::Regular_triangulation_euclidean_traits_2`,
+     `CGAL::Regular_triangulation_filtered_traits_2`. Users must use a version prior to 5.0 if they need these classes.
+ -   **Breaking change**: The [graph traits](https://doc.cgal.org/5.0/BGL/group__PkgBGLTraits.html) enabling CGAL's 2D triangulations to be used as a parameter
+     for any graph-based algorithm of CGAL (or boost) have been improved to fully model the [`FaceGraph`](https://doc.cgal.org/5.0/BGL/classFaceGraph.html) concept.
+     In addition, only the finite simplicies (those not incident to the infinite vertex) of the 2D triangulations
+     are now visibile through this scope. The complete triangulation can still be accessed as a graph,
+     by using the graph traits of the underlying triangulation data structure (usually,
+     [`CGAL::Triangulation_data_structure_2`](https://doc.cgal.org/5.0/TDS_2/classCGAL_1_1Triangulation__data__structure__2.html)).
+ -   **Breaking change**: The `insert()` function
+     of
+     [`CGAL::Triangulation_2`](https://doc.cgal.org/latest/Triangulation_2/classCGAL_1_1Triangulation__2.html)
+     which takes a range of points as argument is now guaranteed to
+     insert the points following the order of `InputIterator`.  Note
+     that this change only affects the base class `Triangulation_2`
+     and not any derived class, such as `Delaunay_triangulation_2`.
+-   Added a new [constructor](https://doc.cgal.org/5.0/Triangulation_2/classCGAL_1_1Triangulation__2.html#a6cfa7d3aaa375a25d217858b49e2eb07=)
+     and [`insert()`](https://doc.cgal.org/5.0/Triangulation_2/classCGAL_1_1Triangulation__2.html#ac5e9bc8adef80dc01a0b31c2d0234545)
+     function to [`CGAL::Triangulation_2`](https://doc.cgal.org/5.0/Triangulation_2/classCGAL_1_1Triangulation__2.html)
+     that takes a range of points with info.
+ -   Introduced a new face base class, [`Triangulation_face_base_with_id_2`](https://doc.cgal.org/5.0/BGL/classCGAL_1_1Triangulation__face__base__with__id__2.html)
+     which enables storing user-defined integer IDs in the face of any 2D triangulation, a precondition to use some
+     BGL algorithms.
+ -   Added range types and functions that return ranges, for example for all vertices, enabling the use of `C++11` `for`-loops.
+     See [this new example](https://doc.cgal.org/5.0/Triangulation_2/Triangulation_2_2for_loop_2_8cpp-example.html) for a usage demonstration.
 
-### 2D and 3D Linear Geometry Kernel
- - Add `ComputeApproximateAngle_3` in the 2D/3D Kernel concept to compute
-   the approximate dihedral angle between 2 vectors. Corresponding functors
-   in the model (`Compute_approximate_angle_3`) and free function (`approximate_angle`)
-   are also added.
- - The following objects are now hashable and thus trivially usable
-   with `std::unordered_set` and `std::unordered_map`:
-   `CGAL::Aff_transformation_2`, `CGAL::Aff_transformation_3`,
-   `CGAL::Bbox_2`, `CGAL::Bbox_3`, `CGAL::Circle_2`,
-   `CGAL::Iso_cuboid_3`, `CGAL::Iso_rectangle_2`, `CGAL::Point_2`,
-   `CGAL::Point_3`, `CGAL::Segment_2`, `CGAL::Segment_3`,
-   `CGAL::Sphere_3`, `CGAL::Vector_2`, `CGAL::Vector_3`,
-   `CGAL::Weighted_point_2` and `CGAL::Weighted_point_3`.
+### [3D Triangulations](https://doc.cgal.org/5.0/Manual/packages.html#PkgTriangulation3)
+ -   **Breaking change**: The [constructor](https://doc.cgal.org/5.0/Triangulation_3/classCGAL_1_1Triangulation__3.html#a63f67cf6aaadcee14318cf56a36d247a)
+     and the [`insert()`](https://doc.cgal.org/5.0/Triangulation_3/classCGAL_1_1Triangulation__3.html#ad3353128386bbb51f79d0263e7f67337)
+     function of [`CGAL::Triangulation_3`](https://doc.cgal.org/5.0/Triangulation_3/classCGAL_1_1Triangulation__3.html) 
+     which take a range of points as argument are now guaranteed to
+     insert the points following the order of `InputIterator`. Note
+     that this change only affects the base class `Triangulation_3`
+     and not any derived class, such as `Delaunay_triangulation_3`.
+ -   Added constructor and [`insert()`](https://doc.cgal.org/5.0/Triangulation_3/classCGAL_1_1Triangulation__3.html#a8aa85f88733d30aa3ec5385538e13ace)
+     function to `CGAL::Triangulation_3` that takes a range of points with info.
+ -   Added range types and functions that return ranges, for example for all vertices, which enables to use C++11 for-loops.
+     See [this new example](https://doc.cgal.org/5.0/Triangulation_3/Triangulation_3_2for_loop_8cpp-example.html) for a usage demonstration.
 
-### dD Geometry Kernel
-- New exact kernel `Epeck_d`
+### [Surface Mesh](https://doc.cgal.org/5.0/Manual/packages.html#PkgSurfaceMesh)
+ -   Introduced new functions to read and write using the PLY format,
+     [`CGAL::read_ply()`](https://doc.cgal.org/5.0/Surface_mesh/group__PkgSurface__mesh.html#ga42f6ad486ddab74e13d3dc53f511c343) 
+     and [`CGAL::write_ply()`](https://doc.cgal.org/5.0/Surface_mesh/group__PkgSurface__mesh.html#ga77bbb79d449c981895eedb6c3c23bd14), 
+     enabling users to save and load additional property maps of the surface mesh.
+
+###  [CGAL and Solvers](https://doc.cgal.org/5.0/Manual/packages.html#PkgSolverInterface)
+ -   Added [concepts](https://doc.cgal.org/5.0/Solver_interface/group__PkgSolverInterfaceConcepts.html)
+     and [models](https://doc.cgal.org/5.0/Solver_interface/group__PkgSolverInterfaceRef.html) 
+     for solving Mixed Integer Programming (MIP) problems with or without constraints. 
+
+### [3D Boolean Operations on Nef Polyhedra](https://doc.cgal.org/5.0/Manual/packages.html#PkgNef3)
+ -   Added a function to convert a Nef_polyhedron_3 to a polygon soup: [`CGAL::convert_nef_polyhedron_to_polygon_soup()`](https://doc.cgal.org/5.0/Nef_3/group__PkgNef3IOFunctions.html#ga28a9eb4da0cd6153f0c16f7f9eaf6665)
+
+### [IO Streams](https://doc.cgal.org/5.0/Manual/packages.html#PkgStreamSupport)
+- **Breaking change:** The API of [`CGAL::Color`](https://doc.cgal.org/5.0/Stream_support/classCGAL_1_1Color.html) has been cleaned up.
+- Added new functions to support some parts of the WKT file format:
+    * [`CGAL::read_WKT()`](https://doc.cgal.org/5.0/Stream_support/group__PkgStreamSupportRef.html#gad2872abfe6fcf17d705d38567fdd6248)
+    * [`CGAL::read_point_WKT()`](https://doc.cgal.org/5.0/Stream_support/group__PkgStreamSupportRef.html#gadbd2705b183e467507abd2f167446eba)
+    * [`CGAL::read_multi_point_WKT()`](https://doc.cgal.org/5.0/Stream_support/group__PkgStreamSupportRef.html#ga4fb72e49a1fd385bbed35ea20297aa8d)
+    * [`CGAL::read_linestring_WKT()`](https://doc.cgal.org/5.0/Stream_support/group__PkgStreamSupportRef.html#gaaa236308b9da5dbf217ef281fdb55de4)
+    * [`CGAL::read_multi_linestring_WKT()`](https://doc.cgal.org/5.0/Stream_support/group__PkgStreamSupportRef.html#gad6046c7f9d36512b8a014be82c1e2220)
+    * [`CGAL::read_polygon_WKT()`](https://doc.cgal.org/5.0/Stream_support/group__PkgStreamSupportRef.html#gaa36ccd3ac4b3fe3e3fd8a76715c56b9a)
+    * [`CGAL::read_multi_polygon_WKT()`](https://doc.cgal.org/5.0/Stream_support/group__PkgStreamSupportRef.html#ga4ceaa71b9cb3b3f7984bed19afff6fc6)
+    * [`CGAL::write_point_WKT()`](https://doc.cgal.org/5.0/Stream_support/group__PkgStreamSupportRef.html#gab1a2d277b43c218bf128a2056eb53ced)
+    * [`CGAL::write_polygon_WKT()`](https://doc.cgal.org/5.0/Stream_support/group__PkgStreamSupportRef.html#gab5365a4726893aa4f51739ede63f5a09)
+    * [`CGAL::write_linestring_WKT()`](https://doc.cgal.org/5.0/Stream_support/group__PkgStreamSupportRef.html#gaa37ed77d1a01567b93c872a48198efa6)
+    * [`CGAL::write_multi_point_WKT()`](https://doc.cgal.org/5.0/Stream_support/group__PkgStreamSupportRef.html#ga98de4b4e5cccb370febe5daf66bb582d)
+    * [`CGAL::write_multi_polygon_WKT()`](https://doc.cgal.org/5.0/Stream_support/group__PkgStreamSupportRef.html#ga4ded40ab50f57e0b410640e28964935e)
+    * [`CGAL::write_multi_linestring_WKT()`](https://doc.cgal.org/5.0/Stream_support/group__PkgStreamSupportRef.html#ga219987f7a9c0b871c1733aa0c38f26b3)
 
-### IO Streams
--   Added new functions to support some parts of the WKT file format:
-    - `CGAL::read_point_WKT()`
-    - `CGAL::read_multi_point_WKT()`
-    - `CGAL::read_linestring_WKT()`
-    - `CGAL::read_multi_linestring_WKT()`
-    - `CGAL::read_polygon_WKT()`
-    - `CGAL::read_multi_polygon_WKT()`
-    - `CGAL::write_point_WKT()`
-    - `CGAL::write_polygon_WKT()`
-    - `CGAL::write_linestring_WKT()`
-    - `CGAL::write_multi_point_WKT()`
-    - `CGAL::write_multi_polygon_WKT()`
-    - `CGAL::write_multi_linestring_WKT()`
-    - `CGAL:read_WKT()`
 
 Release 4.14
 ------------
diff --git a/Installation/CMakeLists.txt b/Installation/CMakeLists.txt
index f16e99d4a1a..c149c941c4e 100644
--- a/Installation/CMakeLists.txt
+++ b/Installation/CMakeLists.txt
@@ -925,7 +925,10 @@ if(NOT CGAL_HEADER_ONLY)
     ${CMAKE_BINARY_DIR}/config/CGALConfig.cmake
     DESTINATION ${CGAL_INSTALL_CMAKE_DIR} )
 else()
+  configure_file(${CMAKE_CURRENT_LIST_DIR}/lib/cmake/CGAL/CGALConfig-installation-dirs.cmake.in
+    ${CMAKE_BINARY_DIR}/config/CGALConfig-installation-dirs.cmake)
   install(FILES
+    ${CMAKE_BINARY_DIR}/config/CGALConfig-installation-dirs.cmake
     ${CMAKE_CURRENT_LIST_DIR}/lib/cmake/CGAL/CGALConfig.cmake
     DESTINATION ${CGAL_INSTALL_CMAKE_DIR} )
 endif()
diff --git a/Installation/cmake/modules/CGAL_CreateSingleSourceCGALProgram.cmake b/Installation/cmake/modules/CGAL_CreateSingleSourceCGALProgram.cmake
index b6b7c44bd52..18d7f418ce7 100644
--- a/Installation/cmake/modules/CGAL_CreateSingleSourceCGALProgram.cmake
+++ b/Installation/cmake/modules/CGAL_CreateSingleSourceCGALProgram.cmake
@@ -69,6 +69,11 @@ function(create_single_source_cgal_program firstfile )
     add_to_cached_list( CGAL_EXECUTABLE_TARGETS ${exe_name} )
 
     target_link_libraries(${exe_name} PRIVATE CGAL::CGAL)
+    foreach(comp ${CGAL_REQUESTED_COMPONENTS})
+      if(TARGET CGAL::CGAL_${comp})
+        target_link_libraries(${exe_name} PRIVATE CGAL::CGAL_${comp})
+      endif()
+    endforeach()
     if(CGAL_3RD_PARTY_LIBRARIES)
       target_link_libraries(${exe_name} PRIVATE ${CGAL_3RD_PARTY_LIBRARIES})
     endif()
diff --git a/Installation/cmake/modules/CGAL_GeneratorSpecificSettings.cmake b/Installation/cmake/modules/CGAL_GeneratorSpecificSettings.cmake
index 9aaed8dd84f..366378ea74e 100644
--- a/Installation/cmake/modules/CGAL_GeneratorSpecificSettings.cmake
+++ b/Installation/cmake/modules/CGAL_GeneratorSpecificSettings.cmake
@@ -8,9 +8,9 @@ if ( NOT CGAL_GENERATOR_SPECIFIC_SETTINGS_FILE_INCLUDED )
     set(CGAL_AUTO_LINK_ENABLED TRUE)
   endif()
 
-  if     ( MSVC15 )
-    set(CGAL_TOOLSET "vc150")
-    message( STATUS "Using VC15 compiler." )
+  if ( MSVC_TOOLSET_VERSION )
+    set(CGAL_TOOLSET "vc${MSVC_TOOLSET_VERSION}")
+    message( STATUS "Using VC toolset ${MSVC_TOOLSET_VERSION}." )
   elseif ( MSVC14 )
     set(CGAL_TOOLSET "vc140")
     message( STATUS "Using VC14 compiler." )
diff --git a/Installation/cmake/modules/CGAL_Qt5_moc_and_resource_files.cmake b/Installation/cmake/modules/CGAL_Qt5_moc_and_resource_files.cmake
index 63609720936..45f531f7f93 100644
--- a/Installation/cmake/modules/CGAL_Qt5_moc_and_resource_files.cmake
+++ b/Installation/cmake/modules/CGAL_Qt5_moc_and_resource_files.cmake
@@ -3,7 +3,7 @@ if(CGAL_Qt5_moc_and_resource_files_included)
 endif()
 set(CGAL_Qt5_moc_and_resource_files_included TRUE)
 
-if(NOT CGAL_HEADER_ONLY)
+if(NOT CGAL_HEADER_ONLY AND CGAL_BUILDING_LIBS)
   qt5_wrap_cpp(_CGAL_Qt5_MOC_FILES_private
     ${CGAL_GRAPHICSVIEW_PACKAGE_DIR}/include/CGAL/Qt/GraphicsViewNavigation.h
     ${CGAL_GRAPHICSVIEW_PACKAGE_DIR}/include/CGAL/Qt/DemosMainWindow.h
@@ -16,6 +16,7 @@ if(NOT CGAL_HEADER_ONLY)
     ${CGAL_GRAPHICSVIEW_PACKAGE_DIR}/include/CGAL/Qt/manipulatedFrame.h
     ${CGAL_GRAPHICSVIEW_PACKAGE_DIR}/include/CGAL/Qt/qglviewer.h
     ${CGAL_GRAPHICSVIEW_PACKAGE_DIR}/include/CGAL/Qt/image_interface.h
+    TARGET CGAL_Qt5
     )
 endif()#CGAL_HEADER_ONLY 
 
diff --git a/Installation/cmake/modules/CGAL_SetupCGALDependencies.cmake b/Installation/cmake/modules/CGAL_SetupCGALDependencies.cmake
index 4900c35e593..f0a2f7c3ea3 100644
--- a/Installation/cmake/modules/CGAL_SetupCGALDependencies.cmake
+++ b/Installation/cmake/modules/CGAL_SetupCGALDependencies.cmake
@@ -96,8 +96,8 @@ function(CGAL_setup_CGAL_dependencies target)
     use_CGAL_LEDA_support(${target} ${keyword})
   endif()
   
-  if (CGAL_HEADER_ONLY)
-    target_compile_definitions(${target} ${keyword} CGAL_HEADER_ONLY=1)
+  if (NOT CGAL_HEADER_ONLY)
+    target_compile_definitions(${target} ${keyword} CGAL_NOT_HEADER_ONLY=1)
   endif()
   if (RUNNING_CGAL_AUTO_TEST OR CGAL_TEST_SUITE)
     target_compile_definitions(${target} ${keyword} CGAL_TEST_SUITE=1)
diff --git a/Installation/cmake/modules/FindTBB.cmake b/Installation/cmake/modules/FindTBB.cmake
index b1a3d884c65..23caa1cbf9a 100644
--- a/Installation/cmake/modules/FindTBB.cmake
+++ b/Installation/cmake/modules/FindTBB.cmake
@@ -192,6 +192,11 @@ endmacro()
 # Get path, convert backslashes as ${ENV_${var}}
 getenv_path(TBB_ROOT)
 
+if(NOT ENV_TBB_ROOT)
+  getenv_path(TBBROOT)
+  set(ENV_TBB_ROOT ${ENV_TBBROOT})
+endif()
+
 # initialize search paths
 set(TBB_PREFIX_PATH ${TBB_ROOT} ${ENV_TBB_ROOT})
 set(TBB_INC_SEARCH_PATH "")
diff --git a/Installation/cmake/modules/list_of_whitelisted_headers.cmake b/Installation/cmake/modules/list_of_whitelisted_headers.cmake
index 600c0197b53..b6055fee1b9 100644
--- a/Installation/cmake/modules/list_of_whitelisted_headers.cmake
+++ b/Installation/cmake/modules/list_of_whitelisted_headers.cmake
@@ -9,22 +9,27 @@ set(list_of_whitelisted_headers_txt [=[
   CGAL/IO/read_ply_points.h
   CGAL/IO/write_ply_points.h
   CGAL/Surface_mesh_parameterization/internal/shortest_path.h
+  CGAL/Three/Edge_container.h
   CGAL/Three/exceptions.h
-  CGAL/Three/Polyhedron_demo_plugin_interface.h
-  CGAL/Three/Scene_interface.h
-  CGAL/Three/Scene_item_with_properties.h
-  CGAL/Three/Scene_zoomable_item_interface.h
-  CGAL/Three/Viewer_interface.h
+  CGAL/Three/Point_container.h
   CGAL/Three/Polyhedron_demo_io_plugin_interface.h
-  CGAL/Three/Scene_draw_interface.h
-  CGAL/Three/Scene_item_config.h
-  CGAL/Three/Scene_print_item_interface.h
-  CGAL/Three/TextRenderer.h
   CGAL/Three/Polyhedron_demo_plugin_helper.h
+  CGAL/Three/Polyhedron_demo_plugin_interface.h
+  CGAL/Three/Primitive_container.h
+  CGAL/Three/Scene_draw_interface.h
   CGAL/Three/Scene_group_item.h
+  CGAL/Three/Scene_interface.h
+  CGAL/Three/Scene_item_config.h
   CGAL/Three/Scene_item.h
+  CGAL/Three/Scene_item_rendering_helper.h
+  CGAL/Three/Scene_item_with_properties.h
+  CGAL/Three/Scene_print_item_interface.h
   CGAL/Three/Scene_transparent_interface.h
+  CGAL/Three/Scene_zoomable_item_interface.h
+  CGAL/Three/TextRenderer.h
+  CGAL/Three/Triangle_container.h
   CGAL/Three/Viewer_config.h
+  CGAL/Three/Viewer_interface.h
 
 ]=])
 
diff --git a/Installation/include/CGAL/auto_link/auto_link.h b/Installation/include/CGAL/auto_link/auto_link.h
index f09070d95a6..b964a3ba097 100644
--- a/Installation/include/CGAL/auto_link/auto_link.h
+++ b/Installation/include/CGAL/auto_link/auto_link.h
@@ -173,14 +173,21 @@ CGAL_VERSION:        Defined in <CGAL/version.h>
 
 #  elif defined(BOOST_MSVC) && (BOOST_MSVC < 1900)
 
+
      // vc12:
 #    define CGAL_LIB_TOOLSET "vc120"
 
-# elif defined(BOOST_MSVC)
-
+#  elif defined(BOOST_MSVC) && (BOOST_MSVC < 1910)
    // vc14:
 #  define CGAL_LIB_TOOLSET "vc140"
 
+#  elif defined(BOOST_MSVC) && (BOOST_MSVC < 1920)
+// vc141:
+#  define CGAL_LIB_TOOLSET "vc141"
+#  elif defined(BOOST_MSVC)
+// vc142:
+#  define CGAL_LIB_TOOLSET "vc142"
+
 #  elif defined(__BORLANDC__)
 
      // CBuilder 6:
diff --git a/Installation/include/CGAL/config.h b/Installation/include/CGAL/config.h
index 4654728fa22..ead1ddbd949 100644
--- a/Installation/include/CGAL/config.h
+++ b/Installation/include/CGAL/config.h
@@ -30,6 +30,11 @@
 #ifndef CGAL_CONFIG_H
 #define CGAL_CONFIG_H
 
+// CGAL is header-only by default since CGAL-5.0.
+#if !defined(CGAL_HEADER_ONLY) && ! CGAL_NOT_HEADER_ONLY
+#  define CGAL_HEADER_ONLY 1
+#endif
+
 #ifdef CGAL_HEADER_ONLY
 #  define CGAL_NO_AUTOLINK 1
 #endif
@@ -132,7 +137,7 @@
 
 // workaround for the bug https://svn.boost.org/trac10/ticket/12534
 // That bug was introduced in Boost 1.62 and fixed in 1.63.
-#if BOOST_VERSION >= 106200 && BOOSTS_VERSION < 106300
+#if BOOST_VERSION >= 106200 && BOOST_VERSION < 106300
 #  include <boost/container/flat_map.hpp>
 #endif
 
@@ -361,8 +366,6 @@
 #    define CGAL_BIG_ENDIAN
 #  elif BOOST_ENDIAN_LITTLE_BYTE
 #    define CGAL_LITTLE_ENDIAN
-#  else
-#    error Unknown endianness
 #  endif
 #elif defined (__GLIBC__)
 #  include <endian.h>
@@ -370,8 +373,6 @@
 #    define CGAL_LITTLE_ENDIAN
 #  elif (__BYTE_ORDER == __BIG_ENDIAN)
 #    define CGAL_BIG_ENDIAN
-#  else
-#    error Unknown endianness
 #  endif
 #elif defined(__sparc) || defined(__sparc__) \
    || defined(_POWER) || defined(__powerpc__) \
@@ -385,11 +386,19 @@
    || defined(_M_IX86) || defined(_M_IA64) \
    || defined(_M_ALPHA) || defined(_WIN64)
 #  define CGAL_LITTLE_ENDIAN
-#else
-#  error Unknown endianness
 #endif
 
-
+#if ! defined(CGAL_LITTLE_ENDIAN) && ! defined(CGAL_BIG_ENDIAN)
+#  ifdef CGAL_DEFAULT_IS_LITTLE_ENDIAN
+#    if CGAL_DEFAULT_IS_LITTLE_ENDIAN
+#      define CGAL_LITTLE_ENDIAN
+#    else
+#      define CGAL_BIG_ENDIAN
+#    endif
+#  else
+#    error Unknown endianness: Define CGAL_DEFAULT_IS_LITTLE_ENDIAN to 1 for little endian and to 0 for big endian.
+#  endif
+#endif
 // Symbolic constants to tailor inlining. Inlining Policy.
 // =======================================================
 #ifndef CGAL_MEDIUM_INLINE
@@ -550,7 +559,7 @@ using std::max;
 #  define CGAL_NORETURN  __attribute__ ((__noreturn__))
 #elif defined (_MSC_VER)
 #  define CGAL_NORETURN __declspec(noreturn)
-#else  
+#else
 #  define CGAL_NORETURN
 #endif
 
diff --git a/Installation/include/CGAL/version.h b/Installation/include/CGAL/version.h
index 6d4097f591f..497e9b9a27f 100644
--- a/Installation/include/CGAL/version.h
+++ b/Installation/include/CGAL/version.h
@@ -25,11 +25,11 @@
 #ifndef CGAL_VERSION_H
 #define CGAL_VERSION_H
 
-#define CGAL_VERSION 5.0-beta1
-#define CGAL_VERSION_NR 1050000910
+#define CGAL_VERSION 5.0-beta2
+#define CGAL_VERSION_NR 1050000920
 #define CGAL_SVN_REVISION 99999
 #define CGAL_GIT_HASH abcdef
-#define CGAL_RELEASE_DATE 20190812
+#define CGAL_RELEASE_DATE 20190930
 
 #include <CGAL/version_macros.h>
 
diff --git a/Installation/lib/cmake/CGAL/CGALConfig-installation-dirs.cmake.in b/Installation/lib/cmake/CGAL/CGALConfig-installation-dirs.cmake.in
new file mode 100644
index 00000000000..81ca96bb586
--- /dev/null
+++ b/Installation/lib/cmake/CGAL/CGALConfig-installation-dirs.cmake.in
@@ -0,0 +1 @@
+set(CGAL_ROOT @CMAKE_INSTALL_PREFIX@)
diff --git a/Installation/lib/cmake/CGAL/CGALConfig.cmake b/Installation/lib/cmake/CGAL/CGALConfig.cmake
index 8e2894ac205..50f2ed24e1c 100644
--- a/Installation/lib/cmake/CGAL/CGALConfig.cmake
+++ b/Installation/lib/cmake/CGAL/CGALConfig.cmake
@@ -43,15 +43,17 @@ if(BRANCH_BUILD)
     endif()
   endforeach()
 else()
-  set(CGAL_ROOT ${CGAL_CONFIG_DIR})
-  get_filename_component(CGAL_ROOT "${CGAL_ROOT}" DIRECTORY)
-  if(NOT EXISTS ${CGAL_ROOT}/include/CGAL/config.h)
+  include(${CGAL_CONFIG_DIR}/CGALConfig-installation-dirs.cmake OPTIONAL RESULT_VARIABLE _has_installation_dirs)
+  if(NOT _has_installation_dirs)
+    set(CGAL_ROOT ${CGAL_CONFIG_DIR})
     get_filename_component(CGAL_ROOT "${CGAL_ROOT}" DIRECTORY)
+    if(NOT EXISTS ${CGAL_ROOT}/include/CGAL/config.h)
+      get_filename_component(CGAL_ROOT "${CGAL_ROOT}" DIRECTORY)
+    endif()
+    if(NOT EXISTS ${CGAL_ROOT}/include/CGAL/config.h)
+      get_filename_component(CGAL_ROOT "${CGAL_ROOT}" DIRECTORY)
+    endif()
   endif()
-  if(NOT EXISTS ${CGAL_ROOT}/include/CGAL/config.h)
-    get_filename_component(CGAL_ROOT "${CGAL_ROOT}" DIRECTORY)
-  endif()
-
   # not BRANCH_BUILD: it can be an installed CGAL, or the tarball layout
   if(EXISTS ${CGAL_CONFIG_DIR}/CGAL_add_test.cmake)
     # installed CGAL
@@ -95,7 +97,9 @@ endforeach()
 
 set(CGALConfig_all_targets_are_defined TRUE)
 foreach(cgal_lib ${CGAL_LIBRARIES})
-  if(NOT TARGET CGAL::${cgal_lib})
+  if(TARGET CGAL::${cgal_lib})
+    set(${cgal_lib}_FOUND TRUE)
+  else()
     set(CGALConfig_all_targets_are_defined FALSE)
   endif()
 endforeach()
@@ -139,9 +143,6 @@ foreach(cgal_lib ${CGAL_LIBRARIES})
     if(NOT TARGET CGAL::${cgal_lib})
       add_library(CGAL::${cgal_lib} ALIAS ${cgal_lib})
     endif()
-    if(${cgal_lib} STREQUAL CGAL)
-      target_compile_definitions(CGAL INTERFACE CGAL_HEADER_ONLY=1)
-    endif()
     CGAL_setup_target_dependencies(${cgal_lib} INTERFACE)
   endif()
 endforeach()
diff --git a/Installation/src/CGAL/CMakeLists.txt b/Installation/src/CGAL/CMakeLists.txt
index cf661e93fa7..34e2be2081f 100644
--- a/Installation/src/CGAL/CMakeLists.txt
+++ b/Installation/src/CGAL/CMakeLists.txt
@@ -10,8 +10,8 @@ endif()
 
 CGAL_setup_CGAL_dependencies(CGAL ${keyword})
 
-if(CGAL_HEADER_ONLY)
-  target_compile_definitions(CGAL INTERFACE CGAL_HEADER_ONLY=1)
+if(NOT CGAL_HEADER_ONLY)
+  target_compile_definitions(CGAL INTERFACE CGAL_NOT_HEADER_ONLY=1)
 endif()
 
 if(NOT CGAL_DISABLE_GMP)
diff --git a/Linear_cell_complex/include/CGAL/draw_linear_cell_complex.h b/Linear_cell_complex/include/CGAL/draw_linear_cell_complex.h
index 1b750f2dc0b..39dc2dce637 100644
--- a/Linear_cell_complex/include/CGAL/draw_linear_cell_complex.h
+++ b/Linear_cell_complex/include/CGAL/draw_linear_cell_complex.h
@@ -127,29 +127,30 @@ struct DefaultDrawingFunctorLCC
   }
 };
 
-template<class LCC, class Kernel, int dim=LCC::ambient_dimension>
+template<class LCC, class Local_kernel, int dim=LCC::ambient_dimension>
 struct LCC_geom_utils;
 
-template<class LCC, class Kernel>
-struct LCC_geom_utils<LCC, Kernel, 3>
+template<class LCC, class Local_kernel>
+struct LCC_geom_utils<LCC, Local_kernel, 3>
 {
-  static typename Kernel::Vector_3
+  static typename Local_kernel::Vector_3
   get_vertex_normal(const LCC& lcc, typename LCC::Dart_const_handle dh)
   {
-    typename Kernel::Vector_3 n = internal::Geom_utils<typename LCC::Traits>::
+    typename Local_kernel::Vector_3 n = internal::Geom_utils
+      <typename LCC::Traits, Local_kernel>::
       get_local_vector(CGAL::compute_normal_of_cell_0<LCC>(lcc,dh));
     n = n/(CGAL::sqrt(n*n));
     return n;
   }
 };
 
-template<class LCC, class Kernel>
-struct LCC_geom_utils<LCC, Kernel, 2>
+template<class LCC, class Local_kernel>
+struct LCC_geom_utils<LCC, Local_kernel, 2>
 {
-  static typename Kernel::Vector_3
+  static typename Local_kernel::Vector_3
   get_vertex_normal(const LCC&, typename LCC::Dart_const_handle)
   {
-    typename Kernel::Vector_3 n=CGAL::NULL_VECTOR;
+    typename Local_kernel::Vector_3 n=CGAL::NULL_VECTOR;
     return n;
   }
 };
diff --git a/Maintenance/infrastructure/cgal.geometryfactory.com/crontab b/Maintenance/infrastructure/cgal.geometryfactory.com/crontab
index cd4b38dab5c..239636ec6d6 100644
--- a/Maintenance/infrastructure/cgal.geometryfactory.com/crontab
+++ b/Maintenance/infrastructure/cgal.geometryfactory.com/crontab
@@ -23,13 +23,13 @@ LC_CTYPE=en_US.UTF-8
 # The script also updates the manual tools.
 
 # "master" alone
-0 21  *  *  Sun                      cd $HOME/CGAL/create_internal_release && /usr/bin/time $HOME/bin/create_release $HOME/CGAL/branches/master.git --public --do-it || echo ERROR
+0 21  *  *  Sun                      cd $HOME/CGAL/create_internal_release && /usr/bin/time $HOME/bin/create_release $HOME/CGAL/branches/master.git --public --do-it --beta 2 || echo ERROR
 # "integration"
-0 21  *  *  Mon,Tue,Wed,Thu      cd $HOME/CGAL/create_internal_release && /usr/bin/time $HOME/bin/create_release $HOME/CGAL/branches/integration.git $HOME/CGAL/branches/empty-dir --do-it || echo ERROR
+0 21  *  *  Mon,Tue,Wed,Thu,Fri      cd $HOME/CGAL/create_internal_release && /usr/bin/time $HOME/bin/create_release $HOME/CGAL/branches/integration.git $HOME/CGAL/branches/empty-dir --do-it || echo ERROR
 # from branch 4.14
 0 21  *  *  Sat                      cd $HOME/CGAL/create_internal_release-4.14-branch && /usr/bin/time $HOME/bin/create_release $HOME/CGAL/branches/CGAL-4.14-branch.git --public --do-it || echo ERROR
 # from branch 4.13
-0 21  *  *  Fri                      cd $HOME/CGAL/create_internal_release-4.13-branch && /usr/bin/time $HOME/bin/create_release $HOME/CGAL/branches/CGAL-4.13-branch.git --public --do-it || echo ERROR
+#0 21  *  *  Fri                      cd $HOME/CGAL/create_internal_release-4.13-branch && /usr/bin/time $HOME/bin/create_release $HOME/CGAL/branches/CGAL-4.13-branch.git --public --do-it || echo ERROR
 # from branch 4.12
 #0 21  *  *  Sat                      cd $HOME/CGAL/create_internal_release-4.12-branch && /usr/bin/time $HOME/bin/create_release $HOME/CGAL/branches/CGAL-4.12-branch.git --public --do-it || echo ERROR
 # from branch 4.11
diff --git a/Maintenance/public_release/announcement/mailing-beta.eml b/Maintenance/public_release/announcement/mailing-beta.eml
index a3a24b3c777..262f7185e04 100644
--- a/Maintenance/public_release/announcement/mailing-beta.eml
+++ b/Maintenance/public_release/announcement/mailing-beta.eml
@@ -1,52 +1,109 @@
-Subject: CGAL 4.14 Beta 1 Released, Computational Geometry Algorithms Library
+Subject: CGAL 5.0 Beta 1 Released, Computational Geometry Algorithms Library
 Content-Type: text/plain; charset="utf-8"
 Body: 
 
-The CGAL Open Source Project is pleased to announce the release 4.14 Beta 1
+The CGAL Open Source Project is pleased to announce the release 5.0 Beta 1
 of CGAL, the Computational Geometry Algorithms Library.
 
 
-CGAL version 4.14 Beta 1 is a public testing release. It should provide
+CGAL version 5.0 Beta 1 is a public testing release. It should provide
 a solid ground to report bugs that need to be tackled before the
-release of the final version of CGAL 4.14 in September.
+release of the final version of CGAL 5.0 in October.
+
+
+CGAL 5.0 is the first release of CGAL that requires a C++ compiler
+with the support of C++14 or later. The new list of supported
+compilers is:
+
+  -   Visual C++ 14.0 (from Visual Studio 2015 Update 3) or later,
+  -   Gnu g++ 6.3 or later (on Linux or MacOS),
+  -   LLVM Clang version 8.0 or later (on Linux or MacOS), and
+  -   Apple Clang compiler versions 7.0.2 and 10.0.1 (on MacOS).
+
+Since CGAL 4.9, CGAL can be used as a header-only library, with
+dependencies. Since CGAL 5.0, that is now the default, unless
+specified differently in the (optional) CMake configuration.
+
 
 Besides fixes and general enhancement to existing packages, the following
-has changed since CGAL 4.13:
+has changed since CGAL 4.14:
 
 
-### 2D Periodic Hyperbolic Triangulations (new package)
+Polygonal Surface Reconstruction (new package)
 
- -   This package allows the computation of Delaunay triangulations of
-     the Bolza surface.  The Bolza surface is the most symmetric
-     hyperbolic surface of genus 2. Its fundamental domain is the
-     regular hyperbolic octagon with angles π/4 centered at the origin
-     of the Poincaré disk. Triangulations of the Bolza surface can be
-     seen as triangulations of the hyperbolic plane that are periodic
-     in the four directions defined by the sides of this regular
-     octagon.
-
-### 2D Hyperbolic Triangulations (new package)
-
- -   This package allows the computation of Delaunay Triangulations of
-     sets of points in the Poincaré disk, which is one of the
-     conformal models for the hyperbolic plane.
-
-### The Heat Method (new package)
-
--   This package provides an algorithm that solves the single- or
-    multiple-source shortest path problem by returning an
-    approximation of the geodesic distance for all vertices of a
-    triangle mesh to the closest vertex in a given set of source
-    vertices.
-
-### Triangulated Surface Mesh Approximation (new package)
-
--   This package implements the Variational Shape Approximation method
-    to approximate an input surface triangle mesh by a simpler surface
-    triangle mesh.
+-   This package provides a method for piecewise planar object
+    reconstruction from point clouds. The method takes as input an
+    unordered point set sampled from a piecewise planar object and
+    outputs a compact and watertight surface mesh interpolating the
+    input point set. The method assumes that all necessary major planes
+    are provided (or can be extracted from the input point set using the
+    shape detection method described in Point Set Shape Detection, or
+    any other alternative methods).The method can handle arbitrary
+    piecewise planar objects and is capable of recovering sharp features
+    and is robust to noise and outliers. See also the associated blog
+    entry:
+    
+      https://www.cgal.org/2019/08/05/Polygonal_surface_reconstruction/
 
 
-See https://www.cgal.org/2019/03/04/cgal414-beta1/ for a complete list of
+Shape Detection (major changes)
+
+-   BREAKING CHANGE: The concept ShapeDetectionTraits has been renamed
+    to EfficientRANSACTraits.
+-   BREAKING CHANGE: The Shape_detection_3 namespace has been renamed to
+    Shape_detection.
+-   Added a new, generic implementation of region growing. This enables
+    for example applying region growing to inputs such as 2D and 3D
+    point sets, or models of the FaceGraph concept. Learn more about
+    this new algorithm with this blog entry:
+    
+      https://www.cgal.org/2019/07/30/Shape_detection/
+
+
+dD Geometry Kernel
+
+-   A new exact kernel, Epeck_d, is now available.
+
+
+2D and 3D Triangulations
+
+-   BREAKING CHANGE: Several deprecated functions and classes have been
+    removed. See the full list of breaking changes in the release
+    notes.
+
+-   BREAKING CHANGE: The constructor and the insert() function of
+    CGAL::Triangulation_2 or CGAL::Triangulation_3 which take a range
+    of points as argument are now guaranteed to insert the points
+    following the order of InputIterator.  Note that this change only
+    affects the base class CGAL::Triangulation_[23] and not any
+    derived class, such as CGAL::Delaunay_triangulation_[23].
+
+
+Polygon Mesh Processing
+
+-   Introduced a wide range of new functions related to location of
+    queries on a triangle mesh, such as
+    CGAL::Polygon_mesh_processing::locate(Point, Mesh). The location of
+    a point on a triangle mesh is expressed as the pair of a face and
+    the barycentric coordinates of the point in this face, enabling
+    robust manipulation of locations (for example, intersections of two
+    3D segments living within the same face).
+-   Added the mesh smoothing function smooth_mesh(), which can be used
+    to improve the quality of triangle elements based on various
+    geometric characteristics.
+-   Added the shape smoothing function smooth_shape(), which can be used
+    to smooth the surface of a triangle mesh, using the mean curvature
+    flow to perform noise removal.
+
+
+Point Set Processing
+
+-   BREAKING CHANGE: the API using iterators and overloads for optional
+    parameters (deprecated since CGAL 4.12) has been removed. The
+    current (and now only) API uses ranges and Named Parameters.
+
+
+See https://www.cgal.org/2019/09/30/cgal50-beta1/ for a complete list of
 changes.
 
 
diff --git a/Maintenance/release_building/public_release_name b/Maintenance/release_building/public_release_name
index e9cf1a7f937..35f0924c662 100644
--- a/Maintenance/release_building/public_release_name
+++ b/Maintenance/release_building/public_release_name
@@ -1 +1 @@
-CGAL-5.0-beta1
+CGAL-5.0-beta2
diff --git a/Number_types/include/CGAL/Mpzf.h b/Number_types/include/CGAL/Mpzf.h
index d78ccbf5eb4..5a14cf72d95 100644
--- a/Number_types/include/CGAL/Mpzf.h
+++ b/Number_types/include/CGAL/Mpzf.h
@@ -86,10 +86,11 @@
 
 #if defined(BOOST_MSVC)
 #  pragma warning(push)
-#  pragma warning(disable:4146 4244 4267 4800)
+#  pragma warning(disable:4146 4244 4267 4702 4800)
      // warning on - applied on unsigned number
      // conversion with loss of data
      // conversion with loss of data
+     // unreachable code
      // int to bool performance
 #endif
 
diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in
index 82bd44cdbf4..a2b73d0f5d2 100644
--- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in
+++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in
@@ -28,4 +28,6 @@ EXPAND_AS_DEFINED      = CGAL_PMP_NP_TEMPLATE_PARAMETERS \
 EXCLUDE = ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Polygon_mesh_processing/internal
 EXCLUDE_SYMBOLS += experimental
 
-HTML_EXTRA_FILES           =  ${CGAL_PACKAGE_DOC_DIR}/fig/selfintersections.jpg
+HTML_EXTRA_FILES           = ${CGAL_PACKAGE_DOC_DIR}/fig/selfintersections.jpg \
+                             ${CGAL_PACKAGE_DOC_DIR}/fig/mesh_smoothing.png \
+                             ${CGAL_PACKAGE_DOC_DIR}/fig/shape_smoothing.png
diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
index f39fce7d045..eb78d77b85c 100644
--- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
+++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
@@ -70,7 +70,7 @@
 \cgalPkgPicture{hole_filling_ico.png}
  
 \cgalPkgSummaryBegin
-\cgalPkgAuthor{Sébastien Loriot, Jane Tournois, Ilker %O. Yaz}
+\cgalPkgAuthor{Sébastien Loriot, Mael Rouxel-Labbé, Jane Tournois, Ilker %O. Yaz}
 \cgalPkgDesc{This package provides
 a collection of methods and classes for polygon mesh processing,
 ranging from basic operations on simplices, to complex geometry processing algorithms.}
diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt
index 3cf3d142e12..00bf800f868 100644
--- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt
+++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt
@@ -4,7 +4,7 @@ namespace CGAL {
 \anchor Chapter_PolygonMeshProcessing
 
 \cgalAutoToc
-\authors Sébastien Loriot, Jane Tournois, Ilker %O. Yaz
+\authors Sébastien Loriot, Mael Rouxel-Labbé, Jane Tournois, Ilker %O. Yaz
 
 \image html neptun_head.jpg
 \image latex neptun_head.jpg
@@ -874,6 +874,12 @@ and the number of surface patches that are separated by these edges.
 A first version of this package was started by Ilker %O. Yaz and Sébastien Loriot.
 Jane Tournois worked on the finalization of the API, code, and documentation.
 
+A prototype of mesh and shape smoothing was developed during the 2017 edition of the Google Summer of Code
+by Konstantinos Katrioplas, under supervision of Jane Tournois. It was finalized by Mael Rouxel-Labbé and
+integrated in \cgal 5.0.
+
+Functionalities related to mesh and polygon soup reparation have been introduced steadily over multiple versions
+since \cgal 4.10, in joint work between Sébastien Loriot and Mael Rouxel-Labbé.
 
 */
 } /* namespace CGAL */
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 678703009ec..f13f6bffe40 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/bbox.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/bbox.h
@@ -33,6 +33,16 @@
 #include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
 
 
+#ifdef DOXYGEN_RUNNING
+#define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters
+#define CGAL_PMP_NP_CLASS NamedParameters
+#endif
+
+#ifdef DOXYGEN_RUNNING
+#define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters
+#define CGAL_PMP_NP_CLASS NamedParameters
+#endif
+
 namespace CGAL {
 
   namespace Polygon_mesh_processing {
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/extrude.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/extrude.h
index 400f265c4e2..a68c8bbeaf2 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/extrude.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/extrude.h
@@ -34,6 +34,11 @@
 #include <CGAL/boost/graph/Euler_operations.h>
 #include <vector>
 
+#ifdef DOXYGEN_RUNNING
+#define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters
+#define CGAL_PMP_NP_CLASS NamedParameters
+#endif
+
 namespace CGAL {
 namespace Polygon_mesh_processing {
 namespace extrude_impl{
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h
index f13fdb9bbc5..774ca73f713 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h
@@ -69,18 +69,18 @@ namespace internal{
                            const PolygonMesh& pmesh,
                            const NamedParameters& np)
   {
-    using parameters::choose_parameter;
+
     using parameters::get_parameter;
 
     CGAL_assertion(halfedge(v_max, pmesh)!=boost::graph_traits<PolygonMesh>::null_halfedge());
 
     //VertexPointMap
     typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type VPMap;
-    VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
+    VPMap vpmap = CGAL::parameters::choose_parameter(get_parameter(np, internal_np::vertex_point),
                                get_const_property_map(vertex_point, pmesh));
     //Kernel
     typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GT;
-    GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
+    GT gt = CGAL::parameters::choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
 
     //among the incoming edges of `v_max`, find one edge `e` with the minimal slope
     typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
@@ -182,16 +182,16 @@ bool is_outward_oriented(const PolygonMesh& pmesh,
   if (faces(pmesh).first == faces(pmesh).second)
     return true;
 
-  using parameters::choose_parameter;
+
   using parameters::get_parameter;
 
   //VertexPointMap
   typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type VPMap;
-  VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
+  VPMap vpmap = CGAL::parameters::choose_parameter(get_parameter(np, internal_np::vertex_point),
                              get_const_property_map(vertex_point, pmesh));
   //Kernel
   typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GT;
-  GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
+  GT gt = CGAL::parameters::choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
 
   //find the vertex with maximal z coordinate
   internal::Compare_vertex_points_z_3<GT, VPMap> less_z(vpmap, gt);
@@ -482,16 +482,15 @@ void orient(TriangleMesh& tm, const NamedParameters& np)
   CGAL_assertion(is_valid_polygon_mesh(tm));
   CGAL_assertion(is_closed(tm));
 
-  using parameters::choose_parameter;
   using parameters::get_parameter;
 
-  bool orient_outward = choose_parameter(
+  bool orient_outward = CGAL::parameters::choose_parameter(
                           get_parameter(np, internal_np::outward_orientation),true);
 
-  Vpm vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
+  Vpm vpm = CGAL::parameters::choose_parameter(get_parameter(np, internal_np::vertex_point),
                          get_const_property_map(boost::vertex_point, tm));
 
-  Fid_map fid_map = choose_parameter(get_parameter(np, internal_np::face_index),
+  Fid_map fid_map = CGAL::parameters::choose_parameter(get_parameter(np, internal_np::face_index),
                                  get_const_property_map(boost::face_index, tm));
 
   std::vector<std::size_t> face_cc(num_faces(tm), std::size_t(-1));
@@ -589,16 +588,16 @@ void orient_to_bound_a_volume(TriangleMesh& tm,
   if (!is_closed(tm)) return;
   if (!is_triangle_mesh(tm)) return;
 
-  using parameters::choose_parameter;
+
   using parameters::get_parameter;
 
-  bool orient_outward = choose_parameter(
+  bool orient_outward = CGAL::parameters::choose_parameter(
                           get_parameter(np, internal_np::outward_orientation),true);
 
-  Vpm vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
+  Vpm vpm = CGAL::parameters::choose_parameter(get_parameter(np, internal_np::vertex_point),
                          get_const_property_map(boost::vertex_point, tm));
 
-  Fid_map fid_map = choose_parameter(get_parameter(np, internal_np::face_index),
+  Fid_map fid_map = CGAL::parameters::choose_parameter(get_parameter(np, internal_np::face_index),
                                  get_const_property_map(boost::face_index, tm));
 
   std::vector<std::size_t> face_cc(num_faces(tm), std::size_t(-1));
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair.h
index 7c4c1736139..8e8923af87c 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair.h
@@ -62,6 +62,11 @@
 #include <utility>
 #include <vector>
 
+#ifdef DOXYGEN_RUNNING
+#define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters
+#define CGAL_PMP_NP_CLASS NamedParameters
+#endif
+
 namespace CGAL{
 namespace Polygon_mesh_processing {
 namespace debug{
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h
index d9f9b0f42c0..0d50efdc020 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h
@@ -575,7 +575,7 @@ std::size_t merge_duplicate_points_in_polygon_soup(PointRange& points,
 
   const std::size_t removed_points_n = ini_points_n - points.size();
 
-#ifdef CGAL_PMP_REPAIR_POLYGON_SOUP_VERBOSE_PP
+#ifdef CGAL_PMP_REPAIR_POLYGON_SOUP_VERBOSE
   std::cout << "Removed (merged) " << removed_points_n << " duplicate points" << std::endl;
 #endif
 
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_mesh.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_mesh.h
index 5eefc513f33..6db07aae65c 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_mesh.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_mesh.h
@@ -34,6 +34,11 @@
 
 #include <CGAL/property_map.h>
 
+#ifdef DOXYGEN_RUNNING
+#define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters
+#define CGAL_PMP_NP_CLASS NamedParameters
+#endif
+
 namespace CGAL {
 namespace Polygon_mesh_processing {
 
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h
index 740cc875b89..a4ee2971294 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h
@@ -39,6 +39,11 @@
 #include <sstream>
 #endif
 
+#ifdef DOXYGEN_RUNNING
+#define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters
+#define CGAL_PMP_NP_CLASS NamedParameters
+#endif
+
 namespace CGAL {
 namespace Polygon_mesh_processing {
 
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_hole.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_hole.h
index c6ce20869cf..5de2789c1b7 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_hole.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_hole.h
@@ -39,6 +39,11 @@
 
 #include <vector>
 
+#ifdef DOXYGEN_RUNNING
+#define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters
+#define CGAL_PMP_NP_CLASS NamedParameters
+#endif
+
 namespace CGAL {
 
 namespace Polygon_mesh_processing {
diff --git a/Polyhedron/demo/Polyhedron/Plugins/Classification/Point_set_item_classification.h b/Polyhedron/demo/Polyhedron/Plugins/Classification/Point_set_item_classification.h
index 4de6fce22fc..3ddee5310c3 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/Classification/Point_set_item_classification.h
+++ b/Polyhedron/demo/Polyhedron/Plugins/Classification/Point_set_item_classification.h
@@ -331,6 +331,10 @@ class Point_set_item_classification : public Item_classification_base
     std::vector<int> indices (m_points->point_set()->size(), -1);
 
     m_label_probabilities.clear();
+    m_label_probabilities.resize (m_labels.size());
+    for (std::size_t i = 0; i < m_label_probabilities.size(); ++ i)
+      m_label_probabilities[i].resize (m_points->point_set()->size(), -1);
+    
     if (method == 0)
       CGAL::Classification::classify<Concurrency_tag> (*(m_points->point_set()),
                                                        m_labels, classifier,
diff --git a/Polyhedron/demo/Polyhedron/Plugins/IO/OFF_io_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/IO/OFF_io_plugin.cpp
index 76d53b4014b..580023022f7 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/IO/OFF_io_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/IO/OFF_io_plugin.cpp
@@ -2,6 +2,7 @@
 #include "Scene_polygon_soup_item.h"
 #include "Scene_points_with_normal_item.h"
 #include <CGAL/Three/Three.h>
+#include <CGAL/Polygon_mesh_processing/repair.h>
 
 
 #include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
@@ -25,17 +26,17 @@ class Polyhedron_demo_off_plugin :
   Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.IOPluginInterface/1.90" FILE "off_io_plugin.json")
 
 public:
-  bool isDefaultLoader(const Scene_item *item) const 
-  { 
+  bool isDefaultLoader(const Scene_item *item) const
+  {
     if(qobject_cast<const Scene_surface_mesh_item*>(item)
-       || qobject_cast<const Scene_polygon_soup_item*>(item)) 
-      return true; 
+       || qobject_cast<const Scene_polygon_soup_item*>(item))
+      return true;
     return false;
   }
-  bool isDefaultLoader(const QString& name) const 
-  { 
-    if(name == QString("off")) 
-      return true; 
+  bool isDefaultLoader(const QString& name) const
+  {
+    if(name == QString("off"))
+      return true;
     return false;
   }
   QString name() const { return "off_plugin"; }
@@ -44,7 +45,7 @@ public:
   QList<Scene_item*> load(QFileInfo fileinfo, bool& ok, bool add_to_scene=true);
   CGAL::Three::Scene_item* load_off(QFileInfo fileinfo);
   CGAL::Three::Scene_item* load_obj(QFileInfo fileinfo);
-  
+
   bool canSave(const CGAL::Three::Scene_item*);
   bool save(QFileInfo fileinfo,QList<CGAL::Three::Scene_item*>& );
 };
@@ -55,7 +56,7 @@ bool Polyhedron_demo_off_plugin::canLoad(QFileInfo) const {
 
 QList<Scene_item*> Polyhedron_demo_off_plugin::
 load(QFileInfo fileinfo, bool& ok, bool add_to_scene) {
-  
+
   if(fileinfo.size() == 0)
   {
     CGAL::Three::Three::warning( tr("The file you are trying to load is empty."));
@@ -108,10 +109,10 @@ Polyhedron_demo_off_plugin::load_off(QFileInfo fileinfo) {
     std::cerr << "Error! Cannot open file " << (const char*)fileinfo.filePath().toUtf8() << std::endl;
     return NULL;
   }
-  
-  
+
+
   CGAL::File_scanner_OFF scanner( in, false);
-  
+
   // Try to read .off in a point set
   if (scanner.size_of_facets() == 0)
   {
@@ -124,10 +125,10 @@ Polyhedron_demo_off_plugin::load_off(QFileInfo fileinfo) {
       delete item;
       return 0;
     }
-    
+
     return item;
   }
-  
+
   in.seekg(0);
   // Try to read .off in a surface_mesh
   SMesh *surface_mesh = new SMesh();
@@ -181,6 +182,19 @@ Polyhedron_demo_off_plugin::load_off(QFileInfo fileinfo) {
                          tr("%1 isolated vertices found")
                          .arg(item->getNbIsolatedvertices()));
   }
+  typedef boost::function_output_iterator<CGAL::internal::Throw_at_output> OutputIterator;
+  try{
+    CGAL::Polygon_mesh_processing::non_manifold_vertices(*surface_mesh, OutputIterator());
+  }
+  catch( CGAL::internal::Throw_at_output::Throw_at_output_exception& )
+  {
+
+    QApplication::restoreOverrideCursor();
+    QMessageBox::warning((QWidget*)NULL,
+                         tr("Non Manifold Vertices"),
+                         tr("Non-manifold vertices have been found"));
+  }
+
   if(item->isItemMulticolor())
     item->computeItemColorVectorAutomatically(true);
   return item;
@@ -218,9 +232,9 @@ save(QFileInfo fileinfo,QList<CGAL::Three::Scene_item*>& items)
   // This plugin supports point sets, surface_meshes and polygon soups
   const Scene_points_with_normal_item* points_item =
     qobject_cast<const Scene_points_with_normal_item*>(item);
-  const Scene_surface_mesh_item* sm_item = 
+  const Scene_surface_mesh_item* sm_item =
     qobject_cast<const Scene_surface_mesh_item*>(item);
-  const Scene_polygon_soup_item* soup_item = 
+  const Scene_polygon_soup_item* soup_item =
     qobject_cast<const Scene_polygon_soup_item*>(item);
 
   if(!sm_item && !soup_item && !points_item)
diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Volume_plane.h b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Volume_plane.h
index 53796d48b86..bf6e61287a1 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Volume_plane.h
+++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Volume_plane.h
@@ -278,18 +278,37 @@ private:
   QString name(y_tag) const { return tr("Y Slice for %2").arg(name_); }
   QString name(z_tag) const { return tr("Z Slice for %2").arg(name_); }
 
-  double compute_maxDim() const
+  //according to the tag, a,b,c dim change but not the scale. We look for the max dimension of the whole image.
+  //A high scale factor will often go with a low dimesion, to compensate it. So we don't want a max being the
+  //higher scale * the higher dim, hence the tag specialisation.
+//TODO: set the scale factors according to the dimensipon to avoid doing that.
+  double compute_maxDim(x_tag) const
   {
-    double ax((adim_ - 1) * xscale_), ay((adim_ - 1) * yscale_), az((adim_ - 1) * zscale_),
-           bx((bdim_ - 1) * xscale_), by((bdim_ - 1) * yscale_), bz((bdim_ - 1) * zscale_),
-           cx((cdim_ - 1) * xscale_), cy((cdim_ - 1) * yscale_), cz((cdim_ - 1) * zscale_);
+    double max_a((adim_ - 1) * yscale_),
+        max_b((bdim_ - 1) * zscale_),
+        max_c((cdim_ - 1) * xscale_);
 
-    double max_a = (std::max)((std::max)(ax, ay), az);
-    double max_b = (std::max)((std::max)(bx, by), bz);
-    double max_c = (std::max)((std::max)(cx, cy), cz);
     return (std::max)((std::max)(max_a, max_b), max_c);
-
   }
+
+  double compute_maxDim(y_tag) const
+  {
+    double max_a((adim_ - 1) * xscale_),
+        max_b((bdim_ - 1) * zscale_),
+        max_c((cdim_ - 1) * yscale_);
+
+    return (std::max)((std::max)(max_a, max_b), max_c);
+  }
+
+  double compute_maxDim(z_tag) const
+  {
+    double max_a((adim_ - 1) * xscale_),
+        max_b((bdim_ - 1) * yscale_),
+        max_c((cdim_ - 1) * zscale_);
+
+    return (std::max)((std::max)(max_a, max_b), max_c);
+  }
+
   void drawRectangle(x_tag, bool is_loop) const {
 
 
@@ -348,20 +367,23 @@ private:
 
   void drawSpheres(x_tag) const
   {
-      double max_dim = compute_maxDim();
+      double max_dim = compute_maxDim(x_tag());
       sphere_radius = max_dim / 40.0f;
       create_flat_sphere(1.0f, v_spheres, n_spheres, 18);
 
-      c_spheres.push_back(0.0f); c_spheres.push_back((adim_ - 1) * yscale_/2.0f + max_dim/15.0f); c_spheres.push_back(0.0f);
-      c_spheres.push_back(0.0f); c_spheres.push_back((adim_ - 1) * yscale_ ); c_spheres.push_back((bdim_ - 1) * zscale_/2.0f + max_dim/15.0f);
-      c_spheres.push_back(0.0f); c_spheres.push_back((adim_ - 1) * yscale_/2.0f + max_dim/15.0f); c_spheres.push_back((bdim_ - 1 ) * zscale_);
-      c_spheres.push_back(0.0f); c_spheres.push_back(0.0f); c_spheres.push_back((bdim_ - 1) * zscale_/2.0f + max_dim/15.0f);
+      c_spheres.push_back(0.0f); c_spheres.push_back((adim_ - 1) * yscale_/2.0f + 1.1*sphere_radius); c_spheres.push_back(0.0f);
+
+      c_spheres.push_back(0.0f); c_spheres.push_back((adim_ - 1) * yscale_ ); c_spheres.push_back((bdim_ - 1) * zscale_/2.0f + 1.1*sphere_radius);
+
+      c_spheres.push_back(0.0f); c_spheres.push_back((adim_ - 1) * yscale_/2.0f + 1.1*sphere_radius); c_spheres.push_back((bdim_ - 1 ) * zscale_);
+
+      c_spheres.push_back(0.0f); c_spheres.push_back(0.0f); c_spheres.push_back((bdim_ - 1) * zscale_/2.0f + 1.1*sphere_radius);
 
   }
 
   void drawSpheres(y_tag) const
   {
-      double max_dim = compute_maxDim();
+      double max_dim = compute_maxDim(y_tag());
       sphere_radius = max_dim / 40.0f;
       create_flat_sphere(1.0f, v_spheres, n_spheres,18);
 
@@ -373,14 +395,14 @@ private:
 
   void drawSpheres(z_tag) const
   {
-      double max_dim = compute_maxDim();
+      double max_dim = compute_maxDim(z_tag());
       sphere_radius = max_dim / 40.0f;
       create_flat_sphere(1.0f, v_spheres, n_spheres,18);
 
-      c_spheres.push_back(0.0f); c_spheres.push_back((bdim_ - 1) * yscale_/2.0f - max_dim/15.0f); c_spheres.push_back(0.0f);
-      c_spheres.push_back((adim_ - 1) * xscale_/2.0f-max_dim/15.0f); c_spheres.push_back((bdim_ - 1) * yscale_); c_spheres.push_back(0.0f);
-      c_spheres.push_back((adim_ - 1) * xscale_); c_spheres.push_back((bdim_ - 1) * yscale_/2.0f-max_dim/15.0f); c_spheres.push_back(0.0f);
-      c_spheres.push_back((adim_ - 1) * xscale_/2.0f-max_dim/15.0f); c_spheres.push_back(0.0f); c_spheres.push_back(0.0f);
+      c_spheres.push_back(0.0f); c_spheres.push_back((bdim_ - 1) * yscale_/2.0f - 1.1*sphere_radius); c_spheres.push_back(0.0f);
+      c_spheres.push_back((adim_ - 1) * xscale_/2.0f-1.1*sphere_radius); c_spheres.push_back((bdim_ - 1) * yscale_); c_spheres.push_back(0.0f);
+      c_spheres.push_back((adim_ - 1) * xscale_); c_spheres.push_back((bdim_ - 1) * yscale_/2.0f-1.1*sphere_radius); c_spheres.push_back(0.0f);
+      c_spheres.push_back((adim_ - 1) * xscale_/2.0f-1.1*sphere_radius); c_spheres.push_back(0.0f); c_spheres.push_back(0.0f);
   }
 
   CGAL::qglviewer::Constraint* setConstraint(x_tag) {
@@ -544,7 +566,7 @@ void Volume_plane<T>::draw(Viewer_interface *viewer) const {
      bDim() <= 1 ||
      cDim() <=1)
     return;
-  double max_dim = compute_maxDim();
+  double max_dim = compute_maxDim(*this);
   sphere_radius = max_dim/20.0f * sphere_Slider->value()/100.0f;
   getTriangleContainer(1)->getVbo(Tc::Radius)->bind();
   getTriangleContainer(1)->getVao(viewer)->program->setAttributeValue("radius", sphere_radius);
diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Polyhedron_stitching_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Polyhedron_stitching_plugin.cpp
index 011d5e9dc79..79cb4dbec8c 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Polyhedron_stitching_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Polyhedron_stitching_plugin.cpp
@@ -14,25 +14,6 @@
 #include <CGAL/boost/graph/helpers.h>
 #include <boost/graph/filtered_graph.hpp>
 
-template <typename G>
-struct Is_border {
-  const G& g;
-  Is_border(const G& g)
-    : g(g)
-  {}
-
- template <typename Descriptor>
-  bool operator()(const Descriptor& d) const {
-   return is_border(d,g);
-  }
-
-  bool operator()(typename boost::graph_traits<G>::vertex_descriptor d) const {
-    return is_border(d,g) != boost::none;
-  }
-
-};
-
-
 using namespace CGAL::Three;
 class Polyhedron_demo_polyhedron_stitching_plugin :
   public QObject,
@@ -90,30 +71,6 @@ public Q_SLOTS:
 }; // end Polyhedron_demo_polyhedron_stitching_plugin
 
 
-template <typename Poly>
-struct Polyline_visitor
-{
-  Scene_polylines_item* new_item;
-  typename boost::property_map<Poly, CGAL::vertex_point_t>::const_type vpm;
-
-  Polyline_visitor(const Poly& poly, Scene_polylines_item* new_item)
-    : new_item(new_item), vpm(get(CGAL::vertex_point,poly))
-  {}
-
-  void start_new_polyline()
-  {
-    new_item->polylines.push_back( Scene_polylines_item::Polyline() );
-  }
-
-  void add_node(typename boost::graph_traits<Poly>::vertex_descriptor vd)
-  {
-    
-    new_item->polylines.back().push_back(get(vpm,vd));
-  }
-
-  void end_polyline(){}
-};
-
 
 template <typename Item>
 void Polyhedron_demo_polyhedron_stitching_plugin::on_actionDetectBorders_triggered(Scene_interface::Item_id index)
@@ -127,18 +84,16 @@ void Polyhedron_demo_polyhedron_stitching_plugin::on_actionDetectBorders_trigger
 
       FaceGraph* pMesh = item->polyhedron();
       normalize_border(*pMesh);
+      for(auto ed : edges(*pMesh))
+      {
+        if(pMesh->is_border(ed))
+        {
+          new_item->polylines.push_back(Scene_polylines_item::Polyline());
+          new_item->polylines.back().push_back(pMesh->point(pMesh->source(pMesh->halfedge(ed))));
+          new_item->polylines.back().push_back(pMesh->point(pMesh->target(pMesh->halfedge(ed))));
+        }
+      }
 
-
-      typedef boost::filtered_graph<FaceGraph,Is_border<FaceGraph>, Is_border<FaceGraph> > BorderGraph;
-      
-      Is_border<FaceGraph> ib(*pMesh);
-      BorderGraph bg(*pMesh,ib,ib);
-      Polyline_visitor<FaceGraph> polyline_visitor(*pMesh, new_item); 
-      CGAL::split_graph_into_polylines( bg,
-                                        polyline_visitor,
-                                        CGAL::internal::IsTerminalDefault() );
-
-      
       if (new_item->polylines.empty())
         {
           delete new_item;
diff --git a/Polyhedron/demo/Polyhedron/Scene_surface_mesh_item.cpp b/Polyhedron/demo/Polyhedron/Scene_surface_mesh_item.cpp
index 4b80cdaee92..09dd204eff6 100644
--- a/Polyhedron/demo/Polyhedron/Scene_surface_mesh_item.cpp
+++ b/Polyhedron/demo/Polyhedron/Scene_surface_mesh_item.cpp
@@ -260,6 +260,7 @@ struct Scene_surface_mesh_item_priv{
   double volume, area;
   unsigned int number_of_null_length_edges;
   unsigned int number_of_degenerated_faces;
+  bool has_nm_vertices;
   int genus;
   bool self_intersect;
   mutable QSlider* alphaSlider;
@@ -1515,6 +1516,7 @@ invalidate_stats()
 {
   number_of_degenerated_faces = (unsigned int)(-1);
   number_of_null_length_edges = (unsigned int)(-1);
+  has_nm_vertices = false;
   volume = -std::numeric_limits<double>::infinity();
   area = -std::numeric_limits<double>::infinity();
   self_intersect = false;
@@ -1568,11 +1570,30 @@ QString Scene_surface_mesh_item::computeStats(int type)
     }
     faces_aspect_ratio(d->smesh_, min_altitude, min_ar, max_ar, mean_ar);
   }
+  if(type == HAS_NM_VERTICES)
+  {
 
+    d->has_nm_vertices = false;
+    typedef boost::function_output_iterator<CGAL::internal::Throw_at_output> OutputIterator;
+    try{
+      CGAL::Polygon_mesh_processing::non_manifold_vertices(*d->smesh_, OutputIterator());
+    }
+    catch( CGAL::internal::Throw_at_output::Throw_at_output_exception& )
+    {
+      d->has_nm_vertices = true;
+    }
+
+  }
   switch(type)
   {
   case NB_VERTICES:
     return QString::number(num_vertices(*d->smesh_));
+  case HAS_NM_VERTICES:
+  {
+    if(d->has_nm_vertices)
+      return QString("Yes");
+    return QString("No");
+  }
 
   case NB_FACETS:
     return QString::number(num_faces(*d->smesh_));
@@ -1721,14 +1742,15 @@ CGAL::Three::Scene_item::Header_data Scene_surface_mesh_item::header() const
   CGAL::Three::Scene_item::Header_data data;
   //categories
 
-  data.categories.append(std::pair<QString,int>(QString("Properties"),9));
+  data.categories.append(std::pair<QString,int>(QString("Properties"),11));
   data.categories.append(std::pair<QString,int>(QString("Faces"),10));
-  data.categories.append(std::pair<QString,int>(QString("Edges"),7));
-  data.categories.append(std::pair<QString,int>(QString("Angles"),2));
+  data.categories.append(std::pair<QString,int>(QString("Edges"),6));
+  data.categories.append(std::pair<QString,int>(QString("Angles"),3));
 
 
   //titles
   data.titles.append(QString("#Vertices"));
+  data.titles.append(QString("Has Non-manifold Vertices"));
   data.titles.append(QString("#Connected Components"));
   data.titles.append(QString("#Border Edges"));
   data.titles.append(QString("Pure Triangle"));
diff --git a/Polyhedron/demo/Polyhedron/Scene_surface_mesh_item.h b/Polyhedron/demo/Polyhedron/Scene_surface_mesh_item.h
index 0306acce865..ddd39a9f616 100644
--- a/Polyhedron/demo/Polyhedron/Scene_surface_mesh_item.h
+++ b/Polyhedron/demo/Polyhedron/Scene_surface_mesh_item.h
@@ -106,6 +106,7 @@ public:
   //statistics
   enum STATS {
     NB_VERTICES = 0,
+    HAS_NM_VERTICES,
     NB_CONNECTED_COMPOS,
     NB_BORDER_EDGES,
     IS_PURE_TRIANGLE,
diff --git a/Polyhedron/include/CGAL/draw_polyhedron.h b/Polyhedron/include/CGAL/draw_polyhedron.h
index 66cf3a7f3bd..1e005623df8 100644
--- a/Polyhedron/include/CGAL/draw_polyhedron.h
+++ b/Polyhedron/include/CGAL/draw_polyhedron.h
@@ -157,8 +157,8 @@ protected:
     do
     {
       internal::newell_single_step_3
-        (internal::Geom_utils<Kernel>::get_local_point(he->vertex()->point()),
-         internal::Geom_utils<Kernel>::get_local_point(he->next()->vertex()->point()),
+        (this->get_local_point(he->vertex()->point()),
+         this->get_local_point(he->next()->vertex()->point()),
          normal);
       ++nb;
       he=he->next();
diff --git a/Polytope_distance_d/include/CGAL/width_assertions.h b/Polytope_distance_d/include/CGAL/width_assertions.h
index 5df80c2ce3e..1d26635e507 100644
--- a/Polytope_distance_d/include/CGAL/width_assertions.h
+++ b/Polytope_distance_d/include/CGAL/width_assertions.h
@@ -24,7 +24,6 @@
 
 #include <CGAL/license/Polytope_distance_d.h>
 
-
 #ifdef SIMPLIFY
 #define GCD_COMPUTATION 1
 #endif
@@ -98,16 +97,16 @@
  #include<iostream>
 
  #define DEBUGENDL(doit,msg,var)\
- if(doit!=0) std::cout << msg << " " << var << endl;
+ if(doit!=0) std::cout << msg << " " << var << std::endl;
 
  #define DEBUGPRINT(doit,msg,var)\
  if(doit!=0) std::cout << msg << " " << var;
 
  #define DEBUGMSG(doit,msg)\
- if(doit!=0) std::cout << msg << endl;
+ if(doit!=0) std::cout << msg << std::endl;
 
  #define INFOMSG(doit,msg)\
- if(doit!=0) std::cerr<<msg<<endl;
+ if(doit!=0) std::cerr<<msg<<std::endl;
 
 #else
  #define DEBUGENDL(doit,msg,var)
diff --git a/Random_numbers/package_info/Random_numbers/dependencies b/Random_numbers/package_info/Random_numbers/dependencies
index 1598b3438f7..ed142edb24d 100644
--- a/Random_numbers/package_info/Random_numbers/dependencies
+++ b/Random_numbers/package_info/Random_numbers/dependencies
@@ -5,5 +5,6 @@ Kernel_23
 Modular_arithmetic
 Number_types
 Profiling_tools
+Random_numbers
 STL_Extension
 Stream_support
diff --git a/STL_Extension/include/CGAL/Handle.h b/STL_Extension/include/CGAL/Handle.h
index 14304c4a7b8..e01cf446a79 100644
--- a/STL_Extension/include/CGAL/Handle.h
+++ b/STL_Extension/include/CGAL/Handle.h
@@ -75,6 +75,16 @@ class Handle
       return *this;
     }
 
+    void reset()
+    {
+      if (PTR)
+      {
+        if (--PTR->count==0)
+          delete PTR;
+        PTR=0;
+      }
+    }
+
     int
     refs()  const { return PTR->count; }
 
diff --git a/Scripts/developer_scripts/tag_pr_per_release.sh b/Scripts/developer_scripts/tag_pr_per_release.sh
index af182a31c83..17e58a7379d 100644
--- a/Scripts/developer_scripts/tag_pr_per_release.sh
+++ b/Scripts/developer_scripts/tag_pr_per_release.sh
@@ -1,6 +1,9 @@
 #!/bin/bash
 
 # this script requires ghi: https://github.com/stephencelis/ghi
+# See the wiki for how to assign a token to connect without password:
+# https://github.com/stephencelis/ghi/wiki/FAQ
+#
 # example calls within a git repo
 # bash tag_pr_per_release.sh 4.12 4.12.1
 # bash tag_pr_per_release.sh 4.12 4.13
@@ -26,7 +29,7 @@ for i in ${PR_LIST}; do
 done
 
 read -p "Please confirm operation by typing YES? " -n 4 -r
-echo    # (optional) move to a new line
+echo
 if [[ $REPLY =~ ^YES$ ]]; then
 
   for i in ${PR_LIST}; do
diff --git a/Scripts/scripts/cgal_create_CMakeLists b/Scripts/scripts/cgal_create_CMakeLists
index 6e4e7dbfcb9..05324bc4f53 100755
--- a/Scripts/scripts/cgal_create_CMakeLists
+++ b/Scripts/scripts/cgal_create_CMakeLists
@@ -53,6 +53,8 @@
 
 create_cmake_script_with_options()
 {
+  qt4='n'
+
   # parse options file
   if [ -e "$OPTIONS_FILE" ]; then
 
@@ -97,8 +99,6 @@ create_cmake_script_with_options()
 # Created by the script cgal_create_CMakeLists
 # This is the CMake script for compiling a set of CGAL applications.
 
-cmake_minimum_required(VERSION 3.1...3.14)
-
 EOF
   #---------------------------------------------------------------------------
   if [ "$SINGLE_SOURCE" = "n" ]; then
@@ -110,15 +110,22 @@ EOF
   #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
   cat << 'EOF'
 
+
+cmake_minimum_required(VERSION 2.8.11)
+
 # CGAL and its components
 EOF
+if [ -n "$ENABLE_CTEST" ]; then
+  echo "enable_testing()"
+fi
   #---------------------------------------------------------------------------
   # echo "CGAL_COMPONENTS: $CGAL_COMPONENTS"
   if [ ! -z "$CGAL_COMPONENTS" ]; then
     # ensure capitalization
 
-    # CGAL: Core, PDB, ImageIO
+    # CGAL: Core, Qt4, PDB, ImageIO
     CGAL_COMPONENTS=${CGAL_COMPONENTS//[c|C][o|O][r|R][e|E]/Core}
+    CGAL_COMPONENTS=${CGAL_COMPONENTS//[q|Q][t|T]4/Qt4}
     CGAL_COMPONENTS=${CGAL_COMPONENTS//[i|I][m|M][a|A][g|G][e|E][i|I][o|O]/ImageIO}
   
     # external libs
@@ -154,35 +161,120 @@ EOF
     IFS=':'
     for cgal_component in $CGAL_COMPONENTS; do
       COMPONENT=`echo $cgal_component | tr '[:upper:]' '[:lower:]'`
+
+      # for qtmoc
+      if [ "$COMPONENT" = "qt4" ]; then
+        qt4='y'
+      fi
+
     done
     IFS=$OLDIFS
+
   fi
 
-  if [ -n "${CGAL_COMPONENTS}" ]; then
-    echo "find_package( CGAL REQUIRED OPTIONAL_COMPONENTS ${CGAL_COMPONENTS//:/ } )"
-  else
-    echo "find_package( CGAL REQUIRED )"
-  fi
+  echo "find_package( CGAL QUIET COMPONENTS ${CGAL_COMPONENTS//:/ } )"
+
   #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
+  cat << 'EOF'
 
-  if [ ! -z "$BOOST_COMPONENTS" ]; then
-    cat << 'EOF'
+if ( NOT CGAL_FOUND )
+
+  message(STATUS "This project requires the CGAL library, and will not be compiled.")
+  return()  
+
+endif()
+
+EOF
+  #---------------------------------------------------------------------------
+
+  #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
+  cat << 'EOF'
 
 # Boost and its components
 EOF
-    echo "find_package( Boost REQUIRED OPTIONAL_COMPONENTS ${BOOST_COMPONENTS//:/ } )"
-
-  fi # additional Boost components
   #---------------------------------------------------------------------------
 
-  if [ -d include ] ; then
-    cat << 'EOF'
+  if [ ! -z "$BOOST_COMPONENTS" ]; then
+
+    echo "find_package( Boost REQUIRED COMPONENTS ${BOOST_COMPONENTS//:/ } )"
+
+  else
+
+    echo "find_package( Boost REQUIRED )"
+
+  fi # additional Boost components
+
+  #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
+  cat << 'EOF'
+
+if ( NOT Boost_FOUND )
+
+  message(STATUS "This project requires the Boost library, and will not be compiled.")
+
+  return()  
+
+endif()
+EOF
+  #---------------------------------------------------------------------------
+
+
+  #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
+  cat << 'EOF'
 
 # include for local directory
 EOF
+  #---------------------------------------------------------------------------
+
+  if [ -d include ] ; then
     echo 'include_directories( BEFORE include )'
   fi		
 
+  #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
+  cat << 'EOF'
+
+# include for local package
+EOF
+  #---------------------------------------------------------------------------
+
+  # includes for local package
+  if [ -d ../include ] ; then
+    echo 'include_directories( BEFORE ../include )'
+  fi
+
+  if [ ! -z "$PACKAGES" ]; then
+    #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
+    cat << 'EOF'
+
+# include of additional packages
+EOF
+    #-------------------------------------------------------------------------
+  fi
+
+  # Qt4
+  if [ "$qt4" = "y" ]; then
+
+    #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
+    cat << 'EOF'
+
+# Qt4
+set( QT_USE_QTXML     true )
+set( QT_USE_QTMAIN    true )
+set( QT_USE_QTSCRIPT  true )
+set( QT_USE_QTOPENGL  true )
+
+find_package(Qt4)  
+
+if ( NOT QT_FOUND )
+
+  message(STATUS "This project requires the Qt4 library, and will not be compiled.")
+  return()  
+
+endif()
+EOF
+    #-------------------------------------------------------------------------
+
+  fi #qt4
+
   if [ ! -z "$BOOST_COMPONENTS" ]; then
     #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
     cat << 'EOF'
@@ -196,7 +288,7 @@ EOF
     for boost_component in $BOOST_COMPONENTS; do
       BOOST_COMPONENT=`echo $boost_component | tr '[:lower:]' '[:upper:]'`
       echo "add_definitions( \"-DCGAL_USE_BOOST_${BOOST_COMPONENT}\" )"
-      echo "list(APPEND CGAL_3RD_PARTY_LIBRARIES Boost::${boost_component} )"
+      echo "list(APPEND CGAL_3RD_PARTY_LIBRARIES \${Boost_${BOOST_COMPONENT}_LIBRARY} )"
     done
     IFS=$OLDIFS
 
@@ -214,12 +306,33 @@ EOF
     #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
     cat << 'EOF'
 
+
 # Creating entries for all C++ files with "main" routine
 # ##########################################################
+
 EOF
     #-------------------------------------------------------------------------
 
 
+    # add a new line
+    echo
+
+    # Qt4
+    if [ "$qt4" = "y" ]; then
+      #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
+      cat << 'EOF'
+
+if ( CGAL_Qt4_FOUND AND QT_FOUND )
+
+  include( ${QT_USE_FILE} )
+  include_directories( ${QT_INCLUDE_DIR} )  
+
+endif()
+
+EOF
+      #-----------------------------------------------------------------------
+    fi # qt4
+
     for file in `ls *.cc *.cp *.cxx *.cpp *.CPP *.c++ *.C 2> /dev/null | sort` ; do
       # Create an executable for each cpp that  contains a function "main()"
       BASE=`basename $file .cc`
@@ -231,8 +344,32 @@ EOF
       BASE=`basename $BASE .C`
       egrep '\bmain[ \t]*\(' $file >/dev/null 2>&1
       if [ $? -eq 0 ]; then
-        echo "create_single_source_cgal_program( \"$file\" )"
+        if [ "$qt4" = "y" ]; then
+          echo "create_single_source_cgal_program_qt4( \"$file\" )"
+        else
+          echo "create_single_source_cgal_program( \"$file\" )"
+        fi
+        if [ -n "$ENABLE_CTEST" ]; then 
+          if [ -f "$BASE.cin" ] ; then
+            CIN=" < $BASE.cin"
+          else
+            CIN=
+          fi
+          cat <<EOF
+add_test( "$BASE" \${CMAKE_CTEST_COMMAND}
+  --build-and-test "\${CMAKE_CURRENT_SOURCE_DIR}"
+                   "\${CMAKE_CURRENT_BINARY_DIR}"
+  --build-generator "\${CMAKE_GENERATOR}"
+  --build-makeprogram "\${CMAKE_MAKE_PROGRAM}"
+  --build-target $BASE
+  --build-no-clean
+  --build-run-dir "\${CMAKE_CURRENT_SOURCE_DIR}"
+  --test-command sh -c "\${CMAKE_CURRENT_BINARY_DIR}/$BASE$CIN" )
+EOF
+        fi
       fi
+      #add a new line
+      echo 
     done
     
   else #======================================================================
@@ -256,6 +393,73 @@ EOF
       all="$all $file"
     done
 
+    # Qt4
+    if [ "$qt4" = "y" ]; then
+      #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
+      cat << 'EOF'
+
+if ( CGAL_Qt4_FOUND AND QT_FOUND )
+
+  include( ${QT_USE_FILE} )
+  include_directories( ${QT_INCLUDE_DIR} )  
+
+EOF
+      #-----------------------------------------------------------------------
+
+      echo "  # UI files (Qt Designer files)"
+      for file in `ls *.ui 2> /dev/null | sort`; do
+        echo "  qt4_wrap_ui( DT_UI_FILES $file )"
+      done
+      echo
+      echo "  # qrc files (resources files, that contain icons, at least)"
+      for file in `ls *.qrc 2> /dev/null | sort`; do
+        echo "  qt4_add_resources ( DT_RESOURCE_FILES ./$file )"
+      done
+      echo
+      MOC_FILES=""
+      echo "  # use the Qt MOC preprocessor on classes that derives from QObject"
+      for file in `ls include/*.h 2> /dev/null | sort`; do
+        BASE=`basename $file .h`
+        egrep 'Q_OBJECT' $file >/dev/null 2>&1
+        if [ $? -eq 0 ]; then
+          echo "  qt4_generate_moc( include/${BASE}.h ${BASE}.moc )"
+          MOC_FILES="${BASE}.moc $MOC_FILES"
+        fi
+      done
+      for file in `ls *.h 2> /dev/null | sort`; do
+        BASE=`basename $file .h`
+        egrep 'Q_OBJECT' $file >/dev/null 2>&1
+        if [ $? -eq 0 ]; then
+          echo "  qt4_generate_moc( ${BASE}.h ${BASE}.moc )"
+          MOC_FILES="${BASE}.moc $MOC_FILES"
+        fi
+      done
+      for file in `ls *.cc *.cp *.cxx *.cpp *.CPP *.c++ *.C 2> /dev/null | sort` ; do
+        BASE=`basename $file .cc`
+        BASE=`basename $BASE .cp`
+        BASE=`basename $BASE .cxx`
+        BASE=`basename $BASE .cpp`
+        BASE=`basename $BASE .CPP`
+        BASE=`basename $BASE .c++`
+        BASE=`basename $BASE .C`
+        egrep 'Q_OBJECT' $file >/dev/null 2>&1
+        if [ $? -eq 0 ]; then
+          echo "  qt4_generate_moc( ${BASE}.cpp ${BASE}.moc )"
+          MOC_FILES="${BASE}.moc $MOC_FILES"
+        fi
+      done
+
+      #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
+      cat << 'EOF'
+
+endif()
+EOF
+      #-----------------------------------------------------------------------
+
+      all="${all} ${MOC_FILES} \${DT_UI_FILES} \${DT_RESOURCE_FILES}"
+
+    fi # qt4
+
     # no 'cat' here, as variable substitution required
     echo
     echo "add_executable( ${target_name} ${all} )"
@@ -264,11 +468,20 @@ EOF
     echo
     echo "# Link the executable to CGAL and third-party libraries"
     LIBS=""
-    LIBS=$LIBS" \${CGAL_LIBRARIES} \${CGAL_3RD_PARTY_LIBRARIES}"
-    echo "target_link_libraries(${target_name}  $LIBS )"
+    for comp in ${CGAL_COMPONENTS}; do
+        LIBS=$LIBS" CGAL::CGAL_${comp}"
+    done
+    echo "target_link_libraries(${target_name} PRIVATE CGAL::CGAL${LIBS} )"
 
   fi # single source or all files #===========================================
 
+  #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
+  cat << 'EOF'
+EOF
+  #---------------------------------------------------------------------------
+
+echo
+
 }
 
 usage()
@@ -276,9 +489,9 @@ usage()
   echo "Usage: `basename $0` [-s source] [-c cgal-component1:cgal-component2:...] [-b boost-component1:boost-component2:...] [-p] [-o options_file='`pwd`/cgal_cmake_options:$HOME/.cgal_cmake_options_rc'] [-v] [-h]" >&2
   echo >&2
   echo "  -s source If this parameter is given the script will create one single executable for 'source' with all source files; otherwise it creates one executable for each main'ed source." >&2
-  echo "  cgal_componentX - must be a valid cgal component, examples are 'Core','ImageIO'." >&2
+  echo "  cgal_componentX - must be a valid cgal component, examples are 'Core','ImageIO','Qt4' ('benchmark', 'symbolic')." >&2
   echo "  boost_componentX - must be a valid boost component, like 'filesystem', 'program_options'." >&2
-  echo "  -o options_file - file with CGAL_COMPONENT, and BOOST_COMPONENT directives" >&2
+  echo "  -o options_file - file with PACKAGE, DIRECTORY, CGAL_COMPONENT, and BOOST_COMPONENT directives" >&2
   echo "  -v the version" >&2
   echo "  -h this info screen" >&2
   echo  >&2
@@ -363,6 +576,8 @@ while getopts s:c:b:o:phvt OPT; do
           exit 1
         fi
         ;;
+    t)  ENABLE_CTEST='y'
+        ;;
     h)  usage
         exit 0
         ;;
diff --git a/Spatial_searching/include/CGAL/Search_traits_adapter.h b/Spatial_searching/include/CGAL/Search_traits_adapter.h
index 2e2039fc172..55f637ca602 100644
--- a/Spatial_searching/include/CGAL/Search_traits_adapter.h
+++ b/Spatial_searching/include/CGAL/Search_traits_adapter.h
@@ -157,7 +157,8 @@ public:
   public:
     
     No_lvalue_iterator() : point(NULL), idx(0) { }
-    No_lvalue_iterator(const Point& point, std::size_t idx = 0) : point(new Point(point)), idx(idx) { }
+    No_lvalue_iterator(const Point& point) : point(new Point(point)), idx(0) { }
+    No_lvalue_iterator(const Point& point, int) : point(new Point(point)), idx(Base::Dimension::value) { }
     
   private:
     
@@ -279,9 +280,6 @@ class Distance_adapter : public Base_distance {
   PointPropertyMap ppmap;
   typedef typename Base_distance::FT FT;
 
-  CGAL_static_assertion( ( boost::is_same< boost::lvalue_property_map_tag,
-                           typename boost::property_traits<PointPropertyMap>::category
-                         >::value ) );
 public:
     
   Distance_adapter( const PointPropertyMap& ppmap_=PointPropertyMap(),
diff --git a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
index 63e3eeb4a63..daf7f836b56 100644
--- a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
+++ b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
@@ -1758,6 +1758,8 @@ public:
     /// With `check_all_incident_halfedges == false` the function returns `true`, if the incident
     /// halfedge associated to vertex `v` is a border halfedge, or if the vertex is isolated.
     /// \cgalAdvancedEnd
+    /// \attention If the data contained in the `Surface_mesh` is not a 2-manifold, then
+    /// this operation is not guaranteed to return the right result.
   bool is_border(Vertex_index v, bool check_all_incident_halfedges = true) const
     {
         Halfedge_index h(halfedge(v));
diff --git a/Surface_mesh/include/CGAL/draw_surface_mesh.h b/Surface_mesh/include/CGAL/draw_surface_mesh.h
index bfd76350505..029aec8b769 100644
--- a/Surface_mesh/include/CGAL/draw_surface_mesh.h
+++ b/Surface_mesh/include/CGAL/draw_surface_mesh.h
@@ -165,8 +165,8 @@ protected:
     do
     {
       internal::newell_single_step_3
-        (internal::Geom_utils<Kernel>::get_local_point(sm.point(sm.source(he))),
-         internal::Geom_utils<Kernel>::get_local_point(sm.point(sm.target(he))), normal);
+        (this->get_local_point(sm.point(sm.source(he))),
+         this->get_local_point(sm.point(sm.target(he))), normal);
       ++nb;
       he=sm.next(he);
     }
diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/ARAP_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/ARAP_parameterizer_3.h
index e917ec3f89f..fbe9302a683 100644
--- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/ARAP_parameterizer_3.h
+++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/ARAP_parameterizer_3.h
@@ -42,6 +42,9 @@
 
 #if defined(CGAL_EIGEN3_ENABLED)
 #include <CGAL/Eigen_solver_traits.h>
+#ifdef CGAL_SMP_USE_SPARSESUITE_SOLVERS
+#include <Eigen/UmfPackSupport>
+#endif
 #endif
 
 #include <CGAL/assertions.h>
@@ -156,7 +159,14 @@ namespace Surface_mesh_parameterization {
 ///         and `CGAL_EIGEN3_ENABLED` is defined, then an overload of `Eigen_solver_traits`
 ///         is provided as default parameter:
 /// \code
-///   CGAL::Eigen_solver_traits<Eigen::BICGSTAB< Eigen::SparseMatrix<double> > >
+///   CGAL::Eigen_solver_traits<
+///           Eigen::SparseLU<Eigen_sparse_matrix<double>::EigenType> >
+/// \endcode
+///         Moreover, if SparseSuite solvers are available, which is greatly preferable for speed,
+///         then the default parameter is:
+/// \code
+///   CGAL::Eigen_solver_traits<
+///           Eigen::UmfPackLU<Eigen_sparse_matrix<double>::EigenType> >
 /// \endcode
 ///
 /// \sa `CGAL::Surface_mesh_parameterization::Fixed_border_parameterizer_3<TriangleMesh, BorderParameterizer, SolverTraits>`
@@ -175,7 +185,13 @@ public:
   typedef typename Default::Get<
     SolverTraits_,
   #if defined(CGAL_EIGEN3_ENABLED)
-    Eigen_solver_traits< > // defaults to Eigen::BICGSTAB with Eigen_sparse_matrix
+    #ifdef CGAL_SMP_USE_SPARSESUITE_SOLVERS
+      CGAL::Eigen_solver_traits<
+        Eigen::UmfPackLU<Eigen_sparse_matrix<double>::EigenType> >
+    #else
+      CGAL::Eigen_solver_traits<
+        Eigen::SparseLU<Eigen_sparse_matrix<double>::EigenType> >
+    #endif
   #else
     #pragma message("Error: You must either provide 'SolverTraits_' or link CGAL with the Eigen library")
     SolverTraits_ // no parameter provided, and Eigen is not enabled: so don't compile!
@@ -289,19 +305,26 @@ private:
     std::ostringstream out_ss;
     out_ss << filename << iter << ".off" << std::ends;
     std::ofstream out(out_ss.str().c_str());
-    output_uvmap_to_off(mesh, vertices, faces, uvmap, vimap, out);
+    IO::output_uvmap_to_off(mesh, vertices, faces, uvmap, vimap, out);
   }
 
   // Copy the data from two vectors to the UVmap.
   template <typename VertexUVMap,
-            typename VertexIndexMap>
+            typename VertexIndexMap,
+            typename VertexParameterizedMap>
   void assign_solution(const Vector& Xu,
                        const Vector& Xv,
                        const Vertex_set& vertices,
                        VertexUVMap uvmap,
-                       const VertexIndexMap vimap)
+                       const VertexIndexMap vimap,
+                       const VertexParameterizedMap vpmap)
   {
     for(vertex_descriptor vd : vertices) {
+      // The solver might not have managed to exactly constrain the vertex that was marked
+      // as constrained; simply don't update its position.
+      if(get(vpmap, vd))
+        continue;
+
       int index = get(vimap, vd);
       NT u = Xu(index);
       NT v = Xv(index);
@@ -1115,13 +1138,13 @@ private:
       for(vertex_descriptor vd : vertices) {
         if(get(vpmap, vd)) {
           int index = get(vimap, vd);
-          CGAL_postcondition(std::abs(Xu[index] - Bu[index] ) < 1e-10);
-          CGAL_postcondition(std::abs(Xv[index] - Bv[index] ) < 1e-10);
+          CGAL_warning(std::abs(Xu[index] - Bu[index] ) < 1e-7);
+          CGAL_warning(std::abs(Xv[index] - Bv[index] ) < 1e-7);
         }
       }
     )
 
-    assign_solution(Xu, Xv, vertices, uvmap, vimap);
+    assign_solution(Xu, Xv, vertices, uvmap, vimap, vpmap);
     return status;
   }
 
@@ -1339,49 +1362,62 @@ public:
 #endif
 
     // main loop
-    for(unsigned int ite=1; ite<=m_iterations; ++ite)
+    unsigned int ite = 1;
+    for(;;)
     {
       compute_optimal_Lt_matrices(mesh, faces, ctmap, lp, lpmap, uvmap, ltmap);
       status = update_solution(mesh, vertices, ctmap, lp, lpmap, ltmap,
                                                uvmap, vimap, vpmap, A);
 
-      // Output the current situation
+      // Output the current parameterization
 #ifdef CGAL_SMP_ARAP_DEBUG
-      output_uvmap("ARAP_iteration_", ite, mesh, vertices, faces, uvmap vimap);
+      output_uvmap("ARAP_iteration_", ite, mesh, vertices, faces, uvmap, vimap);
 #endif
-      energy_last = energy_this;
-      energy_this = compute_current_energy(mesh, faces, ctmap, lp, lpmap,
-                                                        ltmap, uvmap);
-#ifdef CGAL_PARAMETERIZATION_ARAP_VERBOSE
-      std::cout << "Energy at iteration " << ite << " : " << energy_this << std::endl;
-#endif
-      CGAL_warning(energy_this >= 0);
 
       if(status != OK)
         return status;
 
       // energy based termination
-      if(m_tolerance > 0.0 && ite <= m_iterations) // if tolerance <= 0, don't compute energy
-      {  // also no need compute energy if this iteration is the last iteration
+      if(m_tolerance > 0. && ite <= m_iterations) { // if tolerance <= 0, don't compute energy
+        energy_last = energy_this;
+        energy_this = compute_current_energy(mesh, faces, ctmap, lp, lpmap, ltmap, uvmap);
+
+#ifdef CGAL_PARAMETERIZATION_ARAP_VERBOSE
+        std::cout << "Energy at iteration " << ite << " : " << energy_this << std::endl;
+#endif
+
+        if(energy_this < 0) {
+          // numerical issues can make it so you may get an energy of -1e-17,
+          // but it shouldn't be too wrong
+          CGAL_assertion(energy_this >= - std::numeric_limits<NT>::epsilon());
+          break;
+        }
+
         double energy_diff = std::abs((energy_last - energy_this) / energy_this);
         if(energy_diff < m_tolerance) {
-#ifdef CGAL_PARAMETERIZATION_ARAP_VERBOSE
-          std::cout << "Minimization process ended after: "
-                    << ite + 1 << " iterations. "
-                    << "Energy diff: " << energy_diff << std::endl;
-#endif
           break;
         }
       }
+
+      if(ite >= m_iterations)
+        break;
+      else
+        ++ite;
     }
 
+#ifdef CGAL_PARAMETERIZATION_ARAP_VERBOSE
+    std::cout << "Minimization process ended after: " << ite << " iterations. " << std::endl;
+#endif
+
 #ifdef CGAL_SMP_ARAP_DEBUG
     output_uvmap("ARAP_final_pre_processing.off", mesh, vertices, faces, uvmap, vimap);
 #endif
 
     if(!is_one_to_one_mapping(mesh, faces, uvmap)) {
      // Use post processing to handle flipped elements
-      std::cerr << "Parameterization is not valid; calling post processor" << std::endl;
+#ifdef CGAL_PARAMETERIZATION_ARAP_VERBOSE
+      std::cout << "Parameterization is not valid; calling post processor" << std::endl;
+#endif
       status = post_process(mesh, vertices, faces, bhd, uvmap, vimap);
     }
 
diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/MVC_post_processor_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/MVC_post_processor_3.h
index 1daced092cc..8d0cdba57af 100644
--- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/MVC_post_processor_3.h
+++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/MVC_post_processor_3.h
@@ -291,7 +291,7 @@ private:
   {
     // Build the constrained triangulation
 
-    // Since the border is closed and we are interest in triangles that are outside
+    // Since the border is closed and we are interested in triangles that are outside
     // of the border, we actually only need to insert points on the border
     for(halfedge_descriptor hd : halfedges_around_face(bhd, mesh)) {
       vertex_descriptor s = source(hd, mesh);
diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h
index b4d79c41816..7196f539672 100644
--- a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h
+++ b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h
@@ -126,6 +126,7 @@ private:
   typename GeomTraits::Construct_normal_3              normal_functor;
   typename GeomTraits::Construct_translated_point_3    translated_point_functor;
   typename GeomTraits::Construct_centroid_3            centroid_functor;
+  typename GeomTraits::Collinear_3                     collinear_functor;
 
   Tree tree;
 
@@ -153,11 +154,11 @@ public:
     normal_functor(traits.construct_normal_3_object()),
     translated_point_functor(traits.construct_translated_point_3_object()),
     centroid_functor(traits.construct_centroid_3_object()),
+    collinear_functor(traits.collinear_3_object()),
     tree(create_traits(mesh, vertex_point_map)),
     use_diagonal(use_diagonal) 
   {
     typedef typename boost::property_traits<VertexPointPmap>::reference Point_ref;
-    typename GeomTraits::Collinear_3  collinear = traits.collinear_3_object();
     face_iterator it, end;
     for(it = faces(mesh).begin(), end = faces(mesh).end(); it!=end; it++)
     {
@@ -167,7 +168,7 @@ public:
         Point_ref b(get(vertex_point_map,target(h, mesh)));
         h = next(h, mesh);
         Point_ref c(get(vertex_point_map, target(h, mesh)));
-        bool test = collinear(a,b,c);
+        bool test = collinear_functor(a,b,c);
         if(!test)
           tree.insert(Primitive(it, mesh, vertex_point_map));
     }
@@ -388,10 +389,11 @@ private:
     const Point p2 = get(vertex_point_map,target(next(halfedge(facet,mesh),mesh),mesh));
     const Point p3 = get(vertex_point_map,target(prev(halfedge(facet,mesh),mesh),mesh));
     const Point center  = centroid_functor(p1, p2, p3);
+    if (collinear_functor(p1, p2, p3)) return boost::none;
     Vector normal = normal_functor(p2, p1, p3);
     normal=scale_functor(normal,
                          FT(1.0/std::sqrt(to_double(normal.squared_length()))));
-
+    if (normal!=normal) return boost::none;
     CGAL::internal::SkipPrimitiveFunctor<face_handle>
     skip(facet);
     CGAL::internal::FirstIntersectionVisitor<face_handle>
diff --git a/Surface_mesh_shortest_path/benchmark/Surface_mesh_shortest_path/benchmark_shortest_paths.cpp b/Surface_mesh_shortest_path/benchmark/Surface_mesh_shortest_path/benchmark_shortest_paths.cpp
index e2c8384402f..a37756a682c 100644
--- a/Surface_mesh_shortest_path/benchmark/Surface_mesh_shortest_path/benchmark_shortest_paths.cpp
+++ b/Surface_mesh_shortest_path/benchmark/Surface_mesh_shortest_path/benchmark_shortest_paths.cpp
@@ -152,9 +152,9 @@ typename Traits::Barycentric_coordinates random_coordinates(CGAL::Random& rand)
 {
   typedef typename Traits::FT FT;
   typename Traits::Construct_barycentric_coordinates construct_barycentric_coordinates;
-  FT u = rand.uniform_real(FT(0.0), FT(1.0));
-  FT v = rand.uniform_real(FT(0.0), FT(1.0) - u);
-  return construct_barycentric_coordinates(u, v, FT(1.0) - u - v);
+  FT u = rand.uniform_real(FT(0), FT(1));
+  FT v = rand.uniform_real(FT(0), FT(1) - u);
+  return construct_barycentric_coordinates(u, v, FT(1) - u - v);
 }
 
 template <class Kernel>
diff --git a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h
index a612338a861..9f6d751be34 100644
--- a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h
+++ b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h
@@ -42,6 +42,8 @@
 #include <CGAL/assertions.h>
 #include <CGAL/AABB_tree.h>
 #include <CGAL/Default.h>
+#include <CGAL/enum.h>
+#include <CGAL/number_utils.h>
 
 #include <CGAL/Surface_mesh_shortest_path/barycentric.h>
 #include <CGAL/Surface_mesh_shortest_path/internal/Cone_tree.h>
@@ -169,7 +171,8 @@ public:
   \brief A model of `BidirectionalIterator` to access the source points
 
   \details An iterator becomes invalid if:
-   - the corresponding point is removed (either with `Surface_mesh_shortest_path::remove_source_point()` or `Surface_mesh_shortest_path::remove_all_source_points()`).
+   - the corresponding point is removed (either with `Surface_mesh_shortest_path::remove_source_point()`
+     or `Surface_mesh_shortest_path::remove_all_source_points()`).
    - the structure is re-built (triggered by a shortest path query or a call to `Surface_mesh_shortest_path::build_sequence_tree()`).
    - the structure is cleared (`Surface_mesh_shortest_path::clear()`).
 
@@ -290,10 +293,12 @@ private:
   typedef typename Traits::Point_2 Point_2;
   typedef typename Traits::Vector_2 Vector_2;
 
-  typedef internal::Cone_tree_node<Traits> Cone_tree_node;
-  typedef internal::Cone_expansion_event<Traits> Cone_expansion_event;
+  typedef Surface_mesh_shortest_paths_3::internal::Cone_tree_node<Traits> Cone_tree_node;
+  typedef Surface_mesh_shortest_paths_3::internal::Cone_expansion_event<Traits> Cone_expansion_event;
 
-  typedef std::priority_queue<Cone_expansion_event*, std::vector<Cone_expansion_event*>, internal::Cone_expansion_event_min_priority_queue_comparator<Traits> > Expansion_priqueue;
+  typedef std::priority_queue<Cone_expansion_event*,
+                              std::vector<Cone_expansion_event*>,
+                              Surface_mesh_shortest_paths_3::internal::Cone_expansion_event_min_priority_queue_comparator<Traits> > Expansion_priqueue;
   typedef std::pair<Cone_tree_node*, FT> Node_distance_pair;
 
 private:
@@ -330,8 +335,9 @@ private:
   };
 
 private:
-  Traits m_traits;
+  const Traits& m_traits;
   Triangle_mesh& m_graph;
+
   Vertex_index_map m_vertexIndexMap;
   Halfedge_index_map m_halfedgeIndexMap;
   Face_index_map m_faceIndexMap;
@@ -366,42 +372,50 @@ public:
 
   /// \cond
 
-  std::size_t peak_node_count()
+  std::size_t peak_node_count() const
   {
     return m_peakNodeCount;
   }
 
-  std::size_t current_node_count()
+  std::size_t current_node_count() const
   {
     return m_currentNodeCount;
   }
 
-  std::size_t peak_queue_size()
+  std::size_t peak_queue_size() const
   {
     return m_peakQueueSize;
   }
 
-  std::size_t current_memory_usage()
+  std::size_t current_memory_usage() const
   {
-    std::size_t baseUsage = m_rootNodes.size() * sizeof(Cone_tree_node*) + m_closestToVertices.size() * sizeof(Node_distance_pair);
+    std::size_t baseUsage = m_rootNodes.size() * sizeof(Cone_tree_node*)
+                            + m_closestToVertices.size() * sizeof(Node_distance_pair);
 
     std::size_t finalUsage = baseUsage + sizeof(Cone_tree_node) * m_currentNodeCount;
 
     for (std::size_t i = 0; i < m_faceOccupiers.size(); ++i)
     {
-      finalUsage += (m_faceOccupiers[i].size() * sizeof(Cone_tree_node*)) + sizeof(std::vector<Cone_tree_node*>);
+      finalUsage += (m_faceOccupiers[i].size() * sizeof(Cone_tree_node*))
+                    + sizeof(std::vector<Cone_tree_node*>);
     }
 
     return finalUsage;
   }
 
-  std::size_t peak_memory_usage()
+  std::size_t peak_memory_usage() const
   {
-    std::size_t baseUsage = m_rootNodes.size() * sizeof(Cone_tree_node*) + m_vertexOccupiers.size() * sizeof(Node_distance_pair) + m_closestToVertices.size() * sizeof(Node_distance_pair);
+    std::size_t baseUsage = m_rootNodes.size() * sizeof(Cone_tree_node*)
+                            + m_vertexOccupiers.size() * sizeof(Node_distance_pair)
+                            + m_closestToVertices.size() * sizeof(Node_distance_pair);
 
-    std::size_t peakNodeUsage = baseUsage + (sizeof(Cone_tree_node) * m_peakNodeCount) + ((sizeof(Cone_expansion_event) + sizeof(Cone_expansion_event*)) * m_queueAtPeakNodes);
+    std::size_t peakNodeUsage = baseUsage + (sizeof(Cone_tree_node) * m_peakNodeCount)
+                                + ((sizeof(Cone_expansion_event)
+                                    + sizeof(Cone_expansion_event*)) * m_queueAtPeakNodes);
 
-    std::size_t peakQueueUsage = baseUsage + (sizeof(Cone_expansion_event) + (sizeof(Cone_expansion_event*)) * m_peakQueueSize) + (sizeof(Cone_tree_node) * m_nodesAtPeakQueue);
+    std::size_t peakQueueUsage = baseUsage
+                                 + (sizeof(Cone_expansion_event) + (sizeof(Cone_expansion_event*)) * m_peakQueueSize)
+                                 + (sizeof(Cone_tree_node) * m_nodesAtPeakQueue);
 
     return std::max(peakNodeUsage, peakQueueUsage);
   }
@@ -452,12 +466,15 @@ private:
 #endif
   }
 
-  Point_2 construct_barycenter_in_triangle_2(const Triangle_2& t, const Barycentric_coordinates& b) const
+  Point_2 construct_barycenter_in_triangle_2(const Triangle_2& t,
+                                             const Barycentric_coordinates& b) const
   {
     return construct_barycenter_in_triangle_2(t, b, m_traits);
   }
 
-  static Point_2 construct_barycenter_in_triangle_2(const Triangle_2& t, const Barycentric_coordinates& b, const Traits& traits)
+  static Point_2 construct_barycenter_in_triangle_2(const Triangle_2& t,
+                                                    const Barycentric_coordinates& b,
+                                                    const Traits& traits)
   {
     typename Traits::Construct_vertex_2 cv2(traits.construct_vertex_2_object());
     typename Traits::Construct_barycentric_coordinates_weight cbcw(traits.construct_barycentric_coordinates_weight_object());
@@ -466,12 +483,15 @@ private:
     return cb2(cv2(t, 0), cbcw(b, 0), cv2(t, 1), cbcw(b, 1), cv2(t, 2), cbcw(b, 2));
   }
 
-  Point_3 construct_barycenter_in_triangle_3(const Triangle_3& t, const Barycentric_coordinates& b) const
+  Point_3 construct_barycenter_in_triangle_3(const Triangle_3& t,
+                                             const Barycentric_coordinates& b) const
   {
     return construct_barycenter_in_triangle_3(t, b, m_traits);
   }
 
-  static Point_3 construct_barycenter_in_triangle_3(const Triangle_3& t, const Barycentric_coordinates& b, const Traits& traits)
+  static Point_3 construct_barycenter_in_triangle_3(const Triangle_3& t,
+                                                    const Barycentric_coordinates& b,
+                                                    const Traits& traits)
   {
     typename Traits::Construct_vertex_3 cv3(traits.construct_vertex_3_object());
     typename Traits::Construct_barycentric_coordinates_weight cbcw(traits.construct_barycentric_coordinates_weight_object());
@@ -480,32 +500,38 @@ private:
     return cb3(cv3(t, 0), cbcw(b, 0), cv3(t, 1), cbcw(b, 1), cv3(t, 2), cbcw(b, 2));
   }
 
-  Triangle_3 triangle_from_halfedge(halfedge_descriptor edge) const
+  Triangle_3 triangle_from_halfedge(const halfedge_descriptor edge) const
   {
     return triangle_from_halfedge(edge, m_graph, m_vertexPointMap);
   }
 
-  static Triangle_3 triangle_from_halfedge(halfedge_descriptor edge, const Triangle_mesh& tm)
+  static Triangle_3 triangle_from_halfedge(const halfedge_descriptor edge,
+                                           const Triangle_mesh& tm)
   {
     return triangle_from_halfedge(edge, tm, get(vertex_point, tm));
   }
 
-  static Triangle_3 triangle_from_halfedge(halfedge_descriptor edge, const Triangle_mesh& tm, Vertex_point_map vertexPointMap)
+  static Triangle_3 triangle_from_halfedge(const halfedge_descriptor edge,
+                                           const Triangle_mesh& tm,
+                                           const Vertex_point_map vertexPointMap)
   {
-    return CGAL::internal::triangle_from_halfedge<Triangle_3, Triangle_mesh, Vertex_point_map>(edge, tm, vertexPointMap);
+    return Surface_mesh_shortest_paths_3::internal::triangle_from_halfedge<Triangle_3, Triangle_mesh, Vertex_point_map>(edge, tm, vertexPointMap);
   }
 
-  Triangle_3 triangle_from_face(face_descriptor f) const
+  Triangle_3 triangle_from_face(const face_descriptor f) const
   {
     return triangle_from_face(f, m_graph, m_vertexPointMap);
   }
 
-  static Triangle_3 triangle_from_face(face_descriptor f, const Triangle_mesh& tm)
+  static Triangle_3 triangle_from_face(const face_descriptor f,
+                                       const Triangle_mesh& tm)
   {
     return triangle_from_halfedge(halfedge(f, tm), tm, get(vertex_point, tm));
   }
 
-  static Triangle_3 triangle_from_face(face_descriptor f, const Triangle_mesh& tm, Vertex_point_map vertexPointMap)
+  static Triangle_3 triangle_from_face(const face_descriptor f,
+                                       const Triangle_mesh& tm,
+                                       const Vertex_point_map vertexPointMap)
   {
     return triangle_from_halfedge(halfedge(f, tm), tm, vertexPointMap);
   }
@@ -514,15 +540,18 @@ private:
     Filtering algorithm described in Xin and Wang (2009) "Improving chen and han's algorithm on the discrete geodesic problem."
     https://dl.acm.org/citation.cfm?doid=1559755.1559761
   */
-  bool window_distance_filter(Cone_tree_node* cone, Segment_2 windowSegment, bool reversed)
+  bool window_distance_filter(Cone_tree_node* cone,
+                              const Segment_2& windowSegment,
+                              const bool reversed)
   {
     typename Traits::Construct_vertex_2 cv2(m_traits.construct_vertex_2_object());
     typename Traits::Compute_squared_distance_2 csd2(m_traits.compute_squared_distance_2_object());
 
-    Segment_2 parentEntrySegment = cone->entry_segment();
-    Point_2 v2 = cone->target_point();
-    Point_2 I = cone->source_image();
-    FT d = cone->distance_from_source_to_root();
+    const Segment_2& parentEntrySegment = cone->entry_segment();
+    const Point_2& v2 = cone->target_point();
+    const Point_2& I = cone->source_image();
+    const FT d = cone->distance_from_source_to_root();
+
     FT d1;
     FT d2;
     FT d3;
@@ -560,51 +589,57 @@ private:
     d2 = v2Distance.second;
     d3 = v3Distance.second;
 
-    bool hasD1 = v1Distance.first != nullptr && v1Distance.first != cone->parent();
-    bool hasD2 = v2Distance.first != nullptr && v2Distance.first != cone->parent();
-    bool hasD3 = v3Distance.first != nullptr && v3Distance.first != cone->parent();
+    const bool hasD1 = v1Distance.first != nullptr && v1Distance.first != cone->parent();
+    const bool hasD2 = v2Distance.first != nullptr && v2Distance.first != cone->parent();
+    const bool hasD3 = v3Distance.first != nullptr && v3Distance.first != cone->parent();
 
-    if (hasD1 && (d + CGAL::internal::select_sqrt(csd2(I, B)) > d1 + CGAL::internal::select_sqrt(csd2(v1, B))))
+    if (hasD1 && (d + CGAL::approximate_sqrt(csd2(I, B)) > d1 + CGAL::approximate_sqrt(csd2(v1, B))))
     {
       if (m_debugOutput)
       {
         std::cout << "Filter: d + |I,B| > d1 + |v1,B|: " << std::endl;
-        std::cout << "v1 = " << v1Index << " , " << d1 << " , v2 = " << v2Index << " , " << d2 << " , v3 = " << v3Index << " , " << d3 << std::endl;
+        std::cout << "v1 = " << v1Index << " , " << d1 << " , v2 = " << v2Index
+                  << " , " << d2 << " , v3 = " << v3Index << " , " << d3 << std::endl;
         std::cout << "d = " << d << std::endl;
-        std::cout << "v1,B = " << CGAL::internal::select_sqrt(csd2(v1, B)) << std::endl;
-        std::cout << "I,B = " << CGAL::internal::select_sqrt(csd2(I, B)) << std::endl;
-        std::cout << "I,A = " << CGAL::internal::select_sqrt(csd2(I, A)) << std::endl;
-        std::cout << (d + CGAL::internal::select_sqrt(csd2(I, B))) << " vs. " << (d1 + CGAL::internal::select_sqrt(csd2(v1, B))) << std::endl;
+        std::cout << "v1,B = " << CGAL::approximate_sqrt(csd2(v1, B)) << std::endl;
+        std::cout << "I,B = " << CGAL::approximate_sqrt(csd2(I, B)) << std::endl;
+        std::cout << "I,A = " << CGAL::approximate_sqrt(csd2(I, A)) << std::endl;
+        std::cout << (d + CGAL::approximate_sqrt(csd2(I, B)))
+                  << " vs. " << (d1 + CGAL::approximate_sqrt(csd2(v1, B))) << std::endl;
       }
 
       return false;
     }
 
-    if (hasD2 && (d + CGAL::internal::select_sqrt(csd2(I, A)) > d2 + CGAL::internal::select_sqrt(csd2(v2, A))))
+    if (hasD2 && (d + CGAL::approximate_sqrt(csd2(I, A)) > d2 + CGAL::approximate_sqrt(csd2(v2, A))))
     {
       if (m_debugOutput)
       {
         std::cout << "Filter: d + |I,A| > d2 + |v2,A|: " << std::endl;
-        std::cout << "v1 = " << v1Index << " , " << d1 << " , v2 = " << v2Index << " , " << d2 << " , v3 = " << v3Index << " , " << d3 << std::endl;
+        std::cout << "v1 = " << v1Index << " , " << d1 << " , v2 = " << v2Index
+                  << " , " << d2 << " , v3 = " << v3Index << " , " << d3 << std::endl;
         std::cout << "d = " << d << std::endl;
-        std::cout << "v2,A = " << CGAL::internal::select_sqrt(csd2(v2, A)) << std::endl;
-        std::cout << "I,A = " << CGAL::internal::select_sqrt(csd2(I, A)) << std::endl;
-        std::cout << (d + CGAL::internal::select_sqrt(csd2(I, A))) << " vs. " << (d2 + CGAL::internal::select_sqrt(csd2(v2, A))) << std::endl;
+        std::cout << "v2,A = " << CGAL::approximate_sqrt(csd2(v2, A)) << std::endl;
+        std::cout << "I,A = " << CGAL::approximate_sqrt(csd2(I, A)) << std::endl;
+        std::cout << (d + CGAL::approximate_sqrt(csd2(I, A)))
+                  << " vs. " << (d2 + CGAL::approximate_sqrt(csd2(v2, A))) << std::endl;
       }
 
       return false;
     }
 
-    if (hasD3 && (d + CGAL::internal::select_sqrt(csd2(I, A)) > d3 + CGAL::internal::select_sqrt(csd2(v3, A))))
+    if (hasD3 && (d + CGAL::approximate_sqrt(csd2(I, A)) > d3 + CGAL::approximate_sqrt(csd2(v3, A))))
     {
       if (m_debugOutput)
       {
         std::cout << "Filter: d + |I,A| > d3 + |v3,A|: " << std::endl;
-        std::cout << "v1 = " << v1Index << " , " << d1 << " , v2 = " << v2Index << " , " << d2 << " , v3 = " << v3Index << " , " << d3 << std::endl;
+        std::cout << "v1 = " << v1Index << " , " << d1 << " , v2 = " << v2Index
+                  << " , " << d2 << " , v3 = " << v3Index << " , " << d3 << std::endl;
         std::cout << "d = " << d << std::endl;
-        std::cout << "v3,A = " << CGAL::internal::select_sqrt(csd2(v3, A)) << std::endl;
-        std::cout << "I,A = " << CGAL::internal::select_sqrt(csd2(I, A)) << std::endl;
-        std::cout << (d + CGAL::internal::select_sqrt(csd2(I, A))) << " vs. " << (d3 + CGAL::internal::select_sqrt(csd2(v3, A))) << std::endl;
+        std::cout << "v3,A = " << CGAL::approximate_sqrt(csd2(v3, A)) << std::endl;
+        std::cout << "I,A = " << CGAL::approximate_sqrt(csd2(I, A)) << std::endl;
+        std::cout << (d + CGAL::approximate_sqrt(csd2(I, A)))
+                  << " vs. " << (d3 + CGAL::approximate_sqrt(csd2(v3, A))) << std::endl;
       }
 
       return false;
@@ -616,11 +651,17 @@ private:
   /*
     Push a new node representing crossing the edge to the left of `cone`'s target vertex
   */
-  void expand_left_child(Cone_tree_node* cone, Segment_2 windowSegment)
+  void expand_left_child(Cone_tree_node* cone,
+                         const Segment_2& windowSegment)
   {
     typename Traits::Construct_vertex_2 cv2(m_traits.construct_vertex_2_object());
     typename Traits::Construct_triangle_3_along_segment_2_flattening ft3as2(m_traits.construct_triangle_3_along_segment_2_flattening_object());
 
+    if (m_debugOutput)
+    {
+      std::cout << std::endl << " >>>>>>>>>>>>>>>>>>> Expanding LEFT CHILD <<<<<<<<<<<<<<<<<<<" <<std::endl;
+    }
+
     CGAL_assertion(cone->m_pendingLeftSubtree != nullptr);
 
     cone->m_pendingLeftSubtree = nullptr;
@@ -629,7 +670,10 @@ private:
     {
       Triangle_3 adjacentFace = triangle_from_halfedge(cone->left_child_edge());
       Triangle_2 layoutFace = ft3as2(adjacentFace, 0, cone->left_child_base_segment());
-      Cone_tree_node* child = new Cone_tree_node(m_traits, m_graph, cone->left_child_edge(), layoutFace, cone->source_image(), cone->distance_from_source_to_root(), cv2(windowSegment, 0), cv2(windowSegment, 1), Cone_tree_node::INTERVAL);
+      Cone_tree_node* child = new Cone_tree_node(m_traits, m_graph, cone->left_child_edge(), layoutFace,
+                                                 cone->source_image(), cone->distance_from_source_to_root(),
+                                                 cv2(windowSegment, 0), cv2(windowSegment, 1),
+                                                 Cone_tree_node::INTERVAL);
       node_created();
       cone->set_left_child(child);
       process_node(child);
@@ -643,20 +687,29 @@ private:
   /*
     Push a new node representing crossing the edge to the right of `cone`'s target vertex
   */
-  void expand_right_child(Cone_tree_node* cone, Segment_2 windowSegment)
+  void expand_right_child(Cone_tree_node* cone,
+                          const Segment_2& windowSegment)
   {
     typename Traits::Construct_vertex_2 cv2(m_traits.construct_vertex_2_object());
     typename Traits::Construct_triangle_3_along_segment_2_flattening ft3as2(m_traits.construct_triangle_3_along_segment_2_flattening_object());
 
     CGAL_assertion(cone->m_pendingRightSubtree != nullptr);
-
     cone->m_pendingRightSubtree = nullptr;
 
+    if (m_debugOutput)
+    {
+      std::cout << std::endl << " >>>>>>>>>>>>>>>>>>> Expanding RIGHT CHILD <<<<<<<<<<<<<<<<<<<" <<std::endl;
+    }
+
+
     if (window_distance_filter(cone, windowSegment, true))
     {
       Triangle_3 adjacentFace = triangle_from_halfedge(cone->right_child_edge());
       Triangle_2 layoutFace = ft3as2(adjacentFace, 0, cone->right_child_base_segment());
-      Cone_tree_node* child = new Cone_tree_node(m_traits, m_graph, cone->right_child_edge(), layoutFace, cone->source_image(), cone->distance_from_source_to_root(), cv2(windowSegment, 0), cv2(windowSegment, 1), Cone_tree_node::INTERVAL);
+      Cone_tree_node* child = new Cone_tree_node(m_traits, m_graph, cone->right_child_edge(), layoutFace,
+                                                 cone->source_image(), cone->distance_from_source_to_root(),
+                                                 cv2(windowSegment, 0), cv2(windowSegment, 1),
+                                                 Cone_tree_node::INTERVAL);
       node_created();
       cone->set_right_child(child);
       process_node(child);
@@ -671,7 +724,9 @@ private:
     Determines whether to expand `location` as a face, edge, or vertex root, depending on
     whether it is near to a given edge or vertex, or is an internal face location
   */
-  void expand_root(face_descriptor f, Barycentric_coordinates location, Source_point_iterator sourcePointIt)
+  void expand_root(const face_descriptor f,
+                   const Barycentric_coordinates& location,
+                   Source_point_iterator sourcePointIt)
   {
     typename Traits::Construct_barycentric_coordinates_weight cbcw(m_traits.construct_barycentric_coordinates_weight_object());
     typename Traits::Classify_barycentric_coordinates classify_barycentric_coordinates(m_traits.classify_barycentric_coordinates_object());
@@ -714,12 +769,14 @@ private:
   /*
     Create source nodes facing each edge of `f`, rooted at the given `faceLocation`
   */
-  void expand_face_root(face_descriptor f, Barycentric_coordinates faceLocation, Source_point_iterator sourcePointIt)
+  void expand_face_root(const face_descriptor f,
+                        const Barycentric_coordinates& faceLocation,
+                        Source_point_iterator sourcePointIt)
   {
     typename Traits::Construct_triangle_3_to_triangle_2_projection pt3t2(m_traits.construct_triangle_3_to_triangle_2_projection_object());
     typename Traits::Construct_vertex_2 cv2(m_traits.construct_vertex_2_object());
 
-    halfedge_descriptor start = halfedge(f, m_graph);
+    const halfedge_descriptor start = halfedge(f, m_graph);
     halfedge_descriptor current = start;
 
     Cone_tree_node* faceRoot = new Cone_tree_node(m_traits, m_graph, m_rootNodes.size());
@@ -729,17 +786,21 @@ private:
     if (m_debugOutput)
     {
       typename Traits::Construct_barycentric_coordinates_weight cbcw(m_traits.construct_barycentric_coordinates_weight_object());
-      std::cout << "\tFace Root Expansion: id = " << get(m_faceIndexMap, f) << " , Location = " << cbcw(faceLocation, 0) << " " << cbcw(faceLocation, 1) << " " << cbcw(faceLocation, 2) << " " << std::endl;
+      std::cout << "\tFace Root Expansion: id = " << get(m_faceIndexMap, f)
+                << " , Location = " << cbcw(faceLocation, 0) << " " << cbcw(faceLocation, 1)
+                << " " << cbcw(faceLocation, 2) << " " << std::endl;
     }
 
     for (std::size_t currentVertex = 0; currentVertex < 3; ++currentVertex)
     {
-      Triangle_3 face3d(triangle_from_halfedge(current));
-      Triangle_2 layoutFace(pt3t2(face3d));
-      Barycentric_coordinates rotatedFaceLocation(shifted_coordinates(faceLocation, currentVertex));
-      Point_2 sourcePoint(construct_barycenter_in_triangle_2(layoutFace, rotatedFaceLocation));
+      const Triangle_3 face3d(triangle_from_halfedge(current));
+      const Triangle_2 layoutFace(pt3t2(face3d));
+      const Barycentric_coordinates rotatedFaceLocation(shifted_coordinates(faceLocation, currentVertex));
+      const Point_2 sourcePoint(construct_barycenter_in_triangle_2(layoutFace, rotatedFaceLocation));
 
-      Cone_tree_node* child = new Cone_tree_node(m_traits, m_graph, current, layoutFace, sourcePoint, FT(0.0), cv2(layoutFace, 0), cv2(layoutFace, 2), Cone_tree_node::FACE_SOURCE);
+      Cone_tree_node* child = new Cone_tree_node(m_traits, m_graph, current, layoutFace, sourcePoint,
+                                                 FT(0), cv2(layoutFace, 0), cv2(layoutFace, 2),
+                                                 Cone_tree_node::FACE_SOURCE);
       node_created();
       faceRoot->push_middle_child(child);
 
@@ -759,7 +820,9 @@ private:
   /*
     Create 'source' nodes to each size of the given edge, rooted at the specified parametric location
   */
-  void expand_edge_root(halfedge_descriptor baseEdge, FT t0, FT t1, Source_point_iterator sourcePointIt)
+  void expand_edge_root(const halfedge_descriptor baseEdge,
+                        const FT t0, const FT t1,
+                        Source_point_iterator sourcePointIt)
   {
     typename Traits::Construct_barycenter_2 cb2(m_traits.construct_barycenter_2_object());
     typename Traits::Construct_vertex_2 cv2(m_traits.construct_vertex_2_object());
@@ -767,7 +830,9 @@ private:
 
     if (m_debugOutput)
     {
-      std::cout << "\tEdge Root Expansion: faceA = " << get(m_faceIndexMap, face(baseEdge, m_graph)) << " , faceB = " << get(m_faceIndexMap, face(opposite(baseEdge, m_graph), m_graph)) << " , t0 = " << t0 << " , t1 = " << t1 << std::endl;
+      std::cout << "\tEdge Root Expansion: faceA = " << get(m_faceIndexMap, face(baseEdge, m_graph))
+                << " , faceB = " << get(m_faceIndexMap, face(opposite(baseEdge, m_graph), m_graph))
+                << " , t0 = " << t0 << " , t1 = " << t1 << std::endl;
     }
 
     halfedge_descriptor baseEdges[2];
@@ -800,7 +865,9 @@ private:
         std::cout << "\t\tLocation = " << sourcePoints[side] << std::endl;
       }
 
-      Cone_tree_node* mainChild = new Cone_tree_node(m_traits, m_graph, baseEdges[side], layoutFaces[side], sourcePoints[side], FT(0.0), cv2(layoutFaces[side], 0), cv2(layoutFaces[side], 1), Cone_tree_node::EDGE_SOURCE);
+      Cone_tree_node* mainChild = new Cone_tree_node(m_traits, m_graph, baseEdges[side], layoutFaces[side],
+                                                     sourcePoints[side], FT(0), cv2(layoutFaces[side], 0),
+                                                     cv2(layoutFaces[side], 1), Cone_tree_node::EDGE_SOURCE);
       node_created();
       edgeRoot->push_middle_child(mainChild);
       process_node(mainChild);
@@ -810,19 +877,21 @@ private:
   /*
     Create a 'source' node for each face surrounding the given vertex.
   */
-  void expand_vertex_root(vertex_descriptor vertex, Source_point_iterator sourcePointIt)
+  void expand_vertex_root(const vertex_descriptor vertex,
+                          Source_point_iterator sourcePointIt)
   {
     if (m_debugOutput)
     {
       std::cout << "\tVertex Root Expansion: Vertex = " << get(m_vertexIndexMap, vertex) << std::endl;
     }
 
-    Cone_tree_node* vertexRoot = new Cone_tree_node(m_traits, m_graph, m_rootNodes.size(), prev(halfedge(vertex, m_graph), m_graph));
+    Cone_tree_node* vertexRoot = new Cone_tree_node(m_traits, m_graph, m_rootNodes.size(),
+                                                    prev(halfedge(vertex, m_graph), m_graph));
 
     node_created();
     m_rootNodes.push_back(std::make_pair(vertexRoot, sourcePointIt));
 
-    m_closestToVertices[get(m_vertexIndexMap, vertex)] = Node_distance_pair(vertexRoot, FT(0.0));
+    m_closestToVertices[get(m_vertexIndexMap, vertex)] = Node_distance_pair(vertexRoot, FT(0));
 
     expand_pseudo_source(vertexRoot);
   }
@@ -846,6 +915,7 @@ private:
 
     if (m_debugOutput)
     {
+      std::cout << "expansionVertex: " << get(m_vertexIndexMap, expansionVertex) << std::endl;
       std::cout << "Distance from target to root: " << distanceFromTargetToRoot << std::endl;
     }
 
@@ -855,24 +925,34 @@ private:
 
     do
     {
-      Triangle_3 face3d(triangle_from_halfedge(currentEdge));
-      Triangle_2 layoutFace(pt3t2(face3d));
+      const Triangle_3 face3d(triangle_from_halfedge(currentEdge));
+      const Triangle_2 layoutFace(pt3t2(face3d));
 
       if (m_debugOutput)
       {
-        std::cout << "Expanding PsuedoSource: id = ";
+        std::cout << std::endl << " >>>>>>>>>>>>>>>>>>> Expanding PseudoSource <<<<<<<<<<<<<<<<<<<" <<std::endl;
+        std::cout << "currentEdge: "
+                  << get(m_vertexIndexMap, source(currentEdge, m_graph)) << " "
+                  << get(m_vertexIndexMap, target(currentEdge, m_graph)) << std::endl;
+        std::cout << "face id = ";
         if (face(currentEdge, m_graph) != Graph_traits::null_face())
         {
-          std::cout << get(m_faceIndexMap, face(currentEdge, m_graph));
+          std::cout << get(m_faceIndexMap, face(currentEdge, m_graph)) << std::endl;
+
+          std::cout << "3D face:" << std::endl << face3d[0] << std::endl << face3d[1] << std::endl << face3d[2] << std::endl;
+          std::cout << "current face: " << get(m_faceIndexMap, face(currentEdge, m_graph)) << " gives 2D layout: " << layoutFace << std::endl;
+          std::cout << "source: " << face3d[1] << std::endl;
         }
         else
         {
-          std::cout << "EXTERNAL";
+          std::cout << "EXTERNAL" << std::endl;
         }
-        std::cout << std::endl;
       }
 
-      Cone_tree_node* child = new Cone_tree_node(m_traits, m_graph, currentEdge, layoutFace, cv2(layoutFace, 1), distanceFromTargetToRoot, cv2(layoutFace, 0), cv2(layoutFace, 2), Cone_tree_node::VERTEX_SOURCE);
+      Cone_tree_node* child = new Cone_tree_node(m_traits, m_graph, currentEdge, layoutFace,
+                                                 cv2(layoutFace, 1), distanceFromTargetToRoot,
+                                                 cv2(layoutFace, 0), cv2(layoutFace, 2),
+                                                 Cone_tree_node::VERTEX_SOURCE);
 
       node_created();
       parent->push_middle_child(child);
@@ -887,7 +967,10 @@ private:
   /*
     Returns the intersection of `segment` and the cone defined by the region to the left `leftBoundary` and right of `rightBoundary`
   */
-  bool clip_to_bounds(const Segment_2& segment, const Ray_2& leftBoundary, const Ray_2& rightBoundary, Segment_2& outSegment)
+  bool clip_to_bounds(const Segment_2& segment,
+                      const Ray_2& leftBoundary,
+                      const Ray_2& rightBoundary,
+                      Segment_2& outSegment) const
   {
     typename Traits::Construct_source_2 cs2(m_traits.construct_source_2_object());
     typename Traits::Construct_segment_2 cseg2(m_traits.construct_segment_2_object());
@@ -905,22 +988,22 @@ private:
 
     if (m_debugOutput)
     {
-      std::cout << "Clipping Segment " << segment << " with left = " << leftBoundary << " and right = " << rightBoundary << std::endl;
+      std::cout << "Clipping Segment " << segment << std::endl
+                << "\t with left = " << leftBoundary << std::endl
+                << "\t with right = " << rightBoundary << std::endl;
     }
 
     FT leftT;
     FT rightT;
 
-    CGAL::Orientation leftOrientation = o2(cs2(leftBoundary), cpo2(leftBoundary, 1), cs2(segment));
-
-    if (leftOrientation == CGAL::RIGHT_TURN || leftOrientation == CGAL::COLLINEAR)
+    if (o2(cs2(leftBoundary), cpo2(leftBoundary, 1), cs2(segment)) != CGAL::LEFT_TURN)
     {
       if (m_debugOutput)
       {
         std::cout << "\tLeft is completely covered." << std::endl;
       }
       leftPoint = cs2(segment);
-      leftT = FT(0.0);
+      leftT = FT(0);
     }
     else
     {
@@ -939,7 +1022,7 @@ private:
         Point_2* result = boost::get<Point_2>(&*cgalIntersection);
         FT t0 = pdas2(cs2(segment), ct2(segment), *result);
 
-        if (t0 >= FT(1.00000))
+        if (t0 >= FT(1))
         {
           if (m_debugOutput)
           {
@@ -948,7 +1031,7 @@ private:
 
           return false;
         }
-        else if (t0 <= FT(0.00000))
+        else if (t0 <= FT(0))
         {
           if (m_debugOutput)
           {
@@ -956,7 +1039,7 @@ private:
           }
 
           leftPoint = cs2(segment);
-          leftT = FT(0.0);
+          leftT = FT(0);
         }
         else
         {
@@ -971,16 +1054,14 @@ private:
       }
     }
 
-    CGAL::Orientation rightOrientation = o2(cs2(rightBoundary), cpo2(rightBoundary, 1), ct2(segment));
-
-    if (rightOrientation == CGAL::LEFT_TURN || rightOrientation == CGAL::COLLINEAR)
+    if (o2(cs2(rightBoundary), cpo2(rightBoundary, 1), ct2(segment)) != CGAL::RIGHT_TURN)
     {
       if (m_debugOutput)
       {
-        std::cout << "Right is completely covered." << std::endl;
+        std::cout << "\tRight is completely covered." << std::endl;
       }
       rightPoint = ct2(segment);
-      rightT = FT(1.0);
+      rightT = FT(1);
     }
     else
     {
@@ -999,7 +1080,7 @@ private:
         Point_2* result = boost::get<Point_2>(&*cgalIntersection);
         FT t0 = pdas2(cs2(segment), ct2(segment), *result);
 
-        if (t0 <= FT(0.00000))
+        if (t0 <= FT(0))
         {
           if (m_debugOutput)
           {
@@ -1007,14 +1088,14 @@ private:
           }
           return false;
         }
-        else if (t0 >= FT(1.00000))
+        else if (t0 >= FT(1))
         {
           if (m_debugOutput)
           {
             std::cout << "\tRight is completely covered (secondary check). " << t0 << std::endl;
           }
           rightPoint = ct2(segment);
-          rightT = FT(1.0);
+          rightT = FT(1);
         }
         else
         {
@@ -1047,7 +1128,31 @@ private:
   */
   void process_node(Cone_tree_node* node)
   {
-    typename Traits::Compare_relative_intersection_along_segment_2 crias2(m_traits.compare_relative_intersection_along_segment_2_object());
+    if (m_debugOutput)
+    {
+      std::cout << std::endl << " ---------------- Processing node ---------------" << std::endl;
+      std::cout << "Node: " << node << std::endl;
+      std::cout << "Node type: " << node->node_type() << std::endl;
+      std::cout << "Tree ID: " << node->tree_id() << " at level = " << node->level() << std::endl;
+      std::cout << "\tParent node: " << node->parent() << std::endl;
+      std::cout << "\tParent node type: " << node->parent()->node_type() << std::endl;
+      std::cout << "\tFace = " << node->layout_face() << std::endl;
+      std::cout << "\tVertices = ";
+      halfedge_descriptor current = node->entry_edge();
+      for (std::size_t i = 0; i<3; ++i)
+      {
+        std::cout << get(m_vertexIndexMap, source(current, m_graph)) << " ";
+        current = next(current, m_graph);
+      }
+      std::cout << std::endl;
+      std::cout << "\tSource Image = " << node->source_image() << std::endl;
+      std::cout << "\tEntry Halfedge = (" << get(m_vertexIndexMap, source(node->entry_edge(), m_graph)) << " "
+                                           << get(m_vertexIndexMap, target(node->entry_edge(), m_graph)) << ")" << std::endl;
+      std::cout << "\tTarget vertex = " << get(m_vertexIndexMap, node->target_vertex()) << std::endl;
+
+      std::cout << "\tWindow Left = " << node->window_left() << std::endl;
+      std::cout << "\tWindow Right = " << node->window_right() << std::endl;
+    }
 
     bool leftSide = false;
     bool rightSide = false;
@@ -1057,39 +1162,30 @@ private:
       leftSide = node->has_left_side();
       rightSide = node->has_right_side();
     }
-    else if (node->node_type() == Cone_tree_node::EDGE_SOURCE) 
+    else // node corresponds to a source
     {
-      leftSide = true;
-      rightSide = true;
+      if (node->node_type() == Cone_tree_node::EDGE_SOURCE)
+      {
+        leftSide = true;
+        rightSide = true;
+      }
+      else
+      {
+        leftSide = true;
+        rightSide = false;
+      }
     }
-    else
+
+    if (m_debugOutput)
     {
-      leftSide = true;
-      rightSide = false;
+      std::cout << "\t Has Left : " << (leftSide ? "yes" : "no")
+                << " , Has Right : " << (rightSide ? "yes" : "no") << std::endl;
     }
 
     bool propagateLeft = false;
     bool propagateRight = false;
     bool propagateMiddle = false;
 
-    if (m_debugOutput)
-    {
-      std::cout << " Processing node " << node << " , level = " << node->level() << std::endl;
-      std::cout << "\tFace = " << node->layout_face() << std::endl;
-      std::cout << "\tVertices = ";
-      halfedge_descriptor current = node->entry_edge();
-      for (std::size_t i = 0; i < 3; ++i)
-      {
-        std::cout << get(m_vertexIndexMap, source(current, m_graph)) << " ";
-        current = next(current, m_graph);
-      }
-      std::cout << std::endl;
-      std::cout << "\tSource Image = " << node->source_image() << std::endl;
-      std::cout << "\tWindow Left = " << node->window_left() << std::endl;
-      std::cout << "\tWindow Right = " << node->window_right() << std::endl;
-      std::cout << "\t Has Left : " << (leftSide ? "yes" : "no") << " , Has Right : " << (rightSide ? "yes" : "no") << std::endl;
-    }
-
     if (node->is_source_node() || (leftSide && rightSide))
     {
       if (m_debugOutput)
@@ -1097,66 +1193,85 @@ private:
         std::cout << "\tContains target vertex" << std::endl;
       }
 
-      std::size_t entryEdgeIndex = get(m_halfedgeIndexMap, node->entry_edge());
+      std::size_t entryHalfEdgeIndex = get(m_halfedgeIndexMap, node->entry_edge());
 
-      Node_distance_pair currentOccupier = m_vertexOccupiers[entryEdgeIndex];
+      Node_distance_pair currentOccupier = m_vertexOccupiers[entryHalfEdgeIndex];
       FT currentNodeDistance = node->distance_from_target_to_root();
 
-      bool isLeftOfCurrent = false;
-
       if (m_debugOutput)
       {
-        std::cout << "\t Entry Edge = " << entryEdgeIndex << std::endl;
-        std::cout << "\t Target vertex = " << get(m_vertexIndexMap, node->target_vertex()) << std::endl;
+        std::cout << "\t Distance to target: " << currentNodeDistance << std::endl;
+        std::cout << "\t Is there a current occupier? " << (currentOccupier.first != nullptr) << std::endl;
       }
 
+      // the relative position of the ray between node.source() and node.target_vertex() and the ray
+      // from occupier.source() (-1 left, 0 collinear, 1 right)
+      CGAL::Comparison_result c = CGAL::EQUAL; // initializing to please weak compilers
+
       if (currentOccupier.first != nullptr)
       {
+        CGAL_assertion(node->entry_edge() == currentOccupier.first->entry_edge());
+        CGAL_assertion(node->target_vertex() == currentOccupier.first->target_vertex());
+
+        // for a vertex source, the ray is along the halfedge pointing towards the target
         if (node->is_vertex_node())
         {
-          isLeftOfCurrent = false;
+          if (currentOccupier.first->is_vertex_node())
+            c = CGAL::EQUAL;
+          else
+            c = CGAL::LARGER;
         }
-        else if (currentOccupier.first->is_vertex_node())
+        else if (currentOccupier.first->is_vertex_node()) // node is not a vertex source
         {
-          isLeftOfCurrent = true;
+          c = CGAL::SMALLER;
         }
-        else
+        else // generic case
         {
-          CGAL::Comparison_result comparison = crias2(
-            node->entry_segment(),
-            node->ray_to_target_vertex().supporting_line(),
-            currentOccupier.first->entry_segment(),
-            currentOccupier.first->ray_to_target_vertex().supporting_line()
-          );
-
-          if (comparison == CGAL::SMALLER)
-          {
-            isLeftOfCurrent = true;
-          }
+          // must compute intersections because although entry edges are identical, their 2D representation is not
+          c = m_traits.compare_relative_intersection_along_segment_2_object()(
+                node->entry_segment(),
+                node->ray_to_target_vertex().supporting_line(),
+                currentOccupier.first->entry_segment(),
+                currentOccupier.first->ray_to_target_vertex().supporting_line());
         }
 
         if (m_debugOutput)
         {
-          std::cout << "\t Current occupier = " << currentOccupier.first << std::endl;
+          std::cout << "\t Current occupier, EH ("
+                    << get(m_vertexIndexMap, source(currentOccupier.first->entry_edge(), m_graph)) << " "
+                    << get(m_vertexIndexMap, target(currentOccupier.first->entry_edge(), m_graph)) << ")" << std::endl;
+          std::cout << "\t Current occupier, Source = " << currentOccupier.first->source_image() << std::endl;
           std::cout << "\t Current Occupier Distance = " << currentOccupier.second << std::endl;
-          std::cout << "\t " << (isLeftOfCurrent ? "Left" : "Right") << " of current" << std::endl;
+          std::cout << "\t smaller (-1)/equal (0)/larger (1) comparison? " << c << std::endl;
         }
       }
 
-      if (m_debugOutput)
+      bool is_node_new_occupier = false;
+      if (currentOccupier.first == nullptr)
       {
-        std::cout << "\t New Distance = " << currentNodeDistance << std::endl;
+        m_vertexOccupiers[entryHalfEdgeIndex] = std::make_pair(node, currentNodeDistance);
+        is_node_new_occupier = true;
+      }
+      else
+      {
+        // Only replace the current occupier if the time is _strictly_ larger
+        // and yield the way to vertex sources (cleaner than manipulating 0-length intervals)
+        if (currentOccupier.second > currentNodeDistance ||
+           (currentOccupier.second == currentNodeDistance && node->node_type() == Cone_tree_node::VERTEX_SOURCE))
+        {
+          m_vertexOccupiers[entryHalfEdgeIndex] = std::make_pair(node, currentNodeDistance);
+          is_node_new_occupier = true;
+        }
       }
 
-      if (currentOccupier.first == nullptr || currentOccupier.second > currentNodeDistance)
+      if (is_node_new_occupier)
       {
         if (m_debugOutput)
         {
-          std::cout << "\t Current node is now the occupier" << std::endl;
+          std::cout << "\t Current node is now the occupier of target vertex "
+                    << get(m_vertexIndexMap, node->target_vertex()) << std::endl;
         }
 
-        m_vertexOccupiers[entryEdgeIndex] = std::make_pair(node, currentNodeDistance);
-
         propagateLeft = true;
         propagateRight = true;
 
@@ -1175,9 +1290,10 @@ private:
           }
         }
 
+        // Some branches from the old occupier that has been superseded can now be pruned
         if (currentOccupier.first != nullptr)
         {
-          if (isLeftOfCurrent)
+          if (c == CGAL::SMALLER) // node's ray is left of occupier's ray
           {
             if (currentOccupier.first->get_left_child())
             {
@@ -1189,7 +1305,7 @@ private:
               currentOccupier.first->m_pendingLeftSubtree = nullptr;
             }
           }
-          else
+          else if (c == CGAL::LARGER) // node's ray is right of occupier's ray
           {
             if (currentOccupier.first->get_right_child())
             {
@@ -1203,9 +1319,8 @@ private:
           }
         }
 
+        // Check if node is now the absolute closest node, and replace the current closest as appropriate
         std::size_t targetVertexIndex = get(m_vertexIndexMap, node->target_vertex());
-
-        // Check if this is now the absolute closest node, and replace the current closest as appropriate
         Node_distance_pair currentClosest = m_closestToVertices[targetVertexIndex];
 
         if (m_debugOutput && currentClosest.first != nullptr)
@@ -1213,16 +1328,25 @@ private:
           std::cout << "\t Current Closest Distance = " << currentClosest.second << std::endl;
         }
 
-        if (currentClosest.first == nullptr || currentClosest.second > currentNodeDistance)
+        // If equal times, give priority to vertex sources since it's cleaner and simpler to handle than interval windows
+        if (currentClosest.first == nullptr ||
+           currentClosest.second > currentNodeDistance ||
+           (currentClosest.second == currentNodeDistance && node->node_type() == Cone_tree_node::VERTEX_SOURCE))
         {
           if (m_debugOutput)
           {
-            std::cout << "\t Current node is now the closest" << std::endl;
+            std::cout << "\t Current node is now the closest at target vertex "
+                      << get(m_vertexIndexMap, node->target_vertex()) << std::endl;
           }
 
           // if this is a saddle vertex, then evict previous closest vertex
           if (m_vertexIsPseudoSource[targetVertexIndex])
           {
+            if (m_debugOutput)
+            {
+              std::cout << "\t Vertex " << targetVertexIndex << " is a pseudo-source" << std::endl;
+            }
+
             if (currentClosest.first != nullptr)
             {
               if (m_debugOutput)
@@ -1253,24 +1377,30 @@ private:
           m_closestToVertices[targetVertexIndex] = Node_distance_pair(node, currentNodeDistance);
         }
       }
-      else
+      else // there is already an occupier, at a strictly smaller distance
       {
-        if (isLeftOfCurrent)
-        {
+        // this is an application of "one angle one split"
+        if (c != CGAL::LARGER) // propage on the left if the node's ray is left of the occupier's
           propagateLeft = true;
-        }
-        else if (!node->is_source_node())
-        {
+        if (c != CGAL::SMALLER && !node->is_source_node()) // by convention a source node only points at the left edge
           propagateRight = true;
-        }
+
+        // No point propagating middle because we know the current occupier has a better time
+        // at the target and if the target is a saddle, middle children have already been spawned
+        // when we were in the current occupier.
       }
     }
-    else
+    else // interval node that does not contain the target vertex
     {
       propagateLeft = leftSide;
       propagateRight = rightSide;
     }
 
+    if (m_debugOutput)
+    {
+      std::cout << "Propagate (L/M/R): " << propagateLeft << " " << propagateMiddle << " " << propagateRight << std::endl;
+    }
+
     if (node->level() <= static_cast<std::size_t>(num_faces(m_graph)))
     {
       if (propagateLeft)
@@ -1309,7 +1439,8 @@ private:
       }
       else
       {
-        bool result = clip_to_bounds(parent->left_child_base_segment(), parent->left_boundary(), parent->right_boundary(), leftWindow);
+        bool result = clip_to_bounds(parent->left_child_base_segment(), parent->left_boundary(),
+                                     parent->right_boundary(), leftWindow);
         if (!result)
         {
           if (m_debugOutput)
@@ -1320,30 +1451,38 @@ private:
         }
       }
 
-      FT distanceEstimate = parent->distance_from_source_to_root() + CGAL::internal::select_sqrt(csd2(parent->source_image(), leftWindow));
+      FT distanceEstimate = parent->distance_from_source_to_root()
+                            + CGAL::approximate_sqrt(csd2(parent->source_image(), leftWindow));
 
       if (m_debugOutput)
       {
-        std::cout << "\tPushing Left Child, Segment = " << parent->left_child_base_segment() << " , clipped = " << leftWindow << " , Estimate = " << distanceEstimate << std::endl;
+        std::cout << ">>> Pushing Left Child, Segment = " << parent->left_child_base_segment()
+                  << " , clipped = " << leftWindow << " , Estimate = " << distanceEstimate << std::endl;
       }
 
-      Cone_expansion_event* event = new Cone_expansion_event(parent, distanceEstimate, Cone_expansion_event::LEFT_CHILD, leftWindow);
+      Cone_expansion_event* event = new Cone_expansion_event(parent, distanceEstimate,
+                                                             Cone_expansion_event::LEFT_CHILD, leftWindow);
       parent->m_pendingLeftSubtree = event;
 
       m_expansionPriqueue.push(event);
-
       queue_pushed();
     }
   }
 
   void push_right_child(Cone_tree_node* parent)
   {
+    if (m_debugOutput)
+    {
+      std::cout << "Tentative push of right child..." << std::endl;
+    }
+
     typename Traits::Compute_squared_distance_2 csd2(m_traits.compute_squared_distance_2_object());
 
     if (face(parent->right_child_edge(), m_graph) != Graph_traits::null_face())
     {
       Segment_2 rightWindow;
-      bool result = clip_to_bounds(parent->right_child_base_segment(), parent->left_boundary(), parent->right_boundary(), rightWindow);
+      bool result = clip_to_bounds(parent->right_child_base_segment(), parent->left_boundary(),
+                                   parent->right_boundary(), rightWindow);
 
       if (!result)
       {
@@ -1354,18 +1493,21 @@ private:
         return;
       }
 
-      FT distanceEstimate = parent->distance_from_source_to_root() + CGAL::internal::select_sqrt(csd2(parent->source_image(), rightWindow));
+      FT distanceEstimate = parent->distance_from_source_to_root()
+                            + CGAL::approximate_sqrt(csd2(parent->source_image(), rightWindow));
 
       if (m_debugOutput)
       {
-        std::cout << "\tPushing Right Child, Segment = " << parent->right_child_base_segment() << " , clipped = " << rightWindow << " , Estimate = " << distanceEstimate << std::endl;
+        std::cout << ">>> Pushing Right Child, Segment = " << parent->right_child_base_segment()
+                  << " , clipped = " << rightWindow << " , Estimate = " << distanceEstimate << std::endl;
       }
 
-      Cone_expansion_event* event = new Cone_expansion_event(parent, distanceEstimate, Cone_expansion_event::RIGHT_CHILD, rightWindow);
-      queue_pushed();
+      Cone_expansion_event* event = new Cone_expansion_event(parent, distanceEstimate,
+                                                             Cone_expansion_event::RIGHT_CHILD, rightWindow);
       parent->m_pendingRightSubtree = event;
 
       m_expansionPriqueue.push(event);
+      queue_pushed();
     }
   }
 
@@ -1373,19 +1515,18 @@ private:
   {
     if (m_debugOutput)
     {
-      std::cout << "\tPushing Middle Child, Estimate = " << parent->distance_from_target_to_root() << std::endl;
+      std::cout << ">>> Pushing Middle Child, Estimate = " << parent->distance_from_target_to_root() << std::endl;
     }
 
     Cone_expansion_event* event = new Cone_expansion_event(parent, parent->distance_from_target_to_root(), Cone_expansion_event::PSEUDO_SOURCE);
-
-    queue_pushed();
-
     parent->m_pendingMiddleSubtree = event;
 
     m_expansionPriqueue.push(event);
+    queue_pushed();
   }
 
-  void delete_node(Cone_tree_node* node, bool destruction = false)
+  void delete_node(Cone_tree_node* node,
+                   const bool destruction = false)
   {
     if (node != nullptr)
     {
@@ -1442,15 +1583,16 @@ private:
         delete_node(node->pop_middle_child(), destruction);
       }
 
-      // At the point of destruction, the `Triangle_mesh` referenced may have gone out of scope, we wish to distinguish between deletion with an assumed reference
+      // At the point of destruction, the `Triangle_mesh` referenced may have gone out of scope,
+      // we wish to distinguish between deletion with an assumed reference
       // to the original `Triangle_mesh`, and deletion without
       if (!node->is_root_node() && !destruction)
       {
-        std::size_t entryEdgeIndex = get(m_halfedgeIndexMap, node->entry_edge());
+        std::size_t entryHalfEdgeIndex = get(m_halfedgeIndexMap, node->entry_edge());
 
-        if (m_vertexOccupiers[entryEdgeIndex].first == node)
+        if (m_vertexOccupiers[entryHalfEdgeIndex].first == node)
         {
-          m_vertexOccupiers[entryEdgeIndex].first = nullptr;
+          m_vertexOccupiers[entryHalfEdgeIndex].first = nullptr;
 
           std::size_t targetVertexIndex = get(m_vertexIndexMap, node->target_vertex());
 
@@ -1474,24 +1616,16 @@ private:
     for (boost::tie(current, end) = vertices(m_graph); current != end; ++current)
     {
       std::size_t vertexIndex = get(m_vertexIndexMap, *current);
-
-      if (is_saddle_vertex(*current) || is_boundary_vertex(*current))
-      {
-        m_vertexIsPseudoSource[vertexIndex] = true;
-      }
-      else
-      {
-        m_vertexIsPseudoSource[vertexIndex] = false;
-      }
+      m_vertexIsPseudoSource[vertexIndex] = (is_saddle_vertex(*current) || is_boundary_vertex(*current));
     }
   }
 
-  bool is_saddle_vertex(vertex_descriptor v)
+  bool is_saddle_vertex(const vertex_descriptor v) const
   {
     return m_traits.is_saddle_vertex_object()(v, m_graph, m_vertexPointMap);
   }
 
-  bool is_boundary_vertex(vertex_descriptor v)
+  bool is_boundary_vertex(const vertex_descriptor v) const
   {
     halfedge_descriptor h = halfedge(v, m_graph);
     halfedge_descriptor first = h;
@@ -1518,13 +1652,13 @@ private:
     }
   }
 
-  void reset_algorithm(bool clearFaceLocations = true)
+  void reset_algorithm(const bool clearFaceLocations = true)
   {
     Cone_tree_node* null_value=nullptr;
     m_closestToVertices.resize(num_vertices(m_graph));
-    std::fill(m_closestToVertices.begin(), m_closestToVertices.end(), Node_distance_pair(null_value, FT(0.0)));
+    std::fill(m_closestToVertices.begin(), m_closestToVertices.end(), Node_distance_pair(null_value, FT(0)));
     m_vertexOccupiers.resize(num_halfedges(m_graph));
-    std::fill(m_vertexOccupiers.begin(), m_vertexOccupiers.end(), Node_distance_pair(null_value, FT(0.0)));
+    std::fill(m_vertexOccupiers.begin(), m_vertexOccupiers.end(), Node_distance_pair(null_value, FT(0)));
 
     while (!m_expansionPriqueue.empty())
     {
@@ -1554,7 +1688,9 @@ private:
   }
 
   template <class Visitor>
-  void visit_shortest_path(Cone_tree_node* startNode, const Point_2& startLocation, Visitor& visitor)
+  void visit_shortest_path(Cone_tree_node* startNode,
+                           const Point_2& startLocation,
+                           Visitor& visitor)
   {
     typename Traits::Compute_parametric_distance_along_segment_2 parametric_distance_along_segment_2(m_traits.compute_parametric_distance_along_segment_2_object());
     typename Traits::Construct_ray_2 construct_ray_2(m_traits.construct_ray_2_object());
@@ -1576,7 +1712,7 @@ private:
         case Cone_tree_node::EDGE_SOURCE:
         {
           Segment_2 entrySegment = current->entry_segment();
-          Point_2 currentSourceImage = current->source_image();
+          const Point_2& currentSourceImage = current->source_image();
           Ray_2 rayToLocation(construct_ray_2(currentSourceImage, currentLocation));
 
           LineLineIntersectResult cgalIntersection = intersect_2(construct_line_2(entrySegment), construct_line_2(rayToLocation));
@@ -1605,10 +1741,13 @@ private:
             std::cout << "Inside cone: " << (current->inside_window(currentLocation) ? "Yes" : "No") << std::endl;
             std::cout << "Current Source: " << current->source_image() << std::endl;
             std::cout << "Current Segment: " << entrySegment << std::endl;
-            std::cout << "Current Left Window: " << current->window_left() << "  ,  " << m_traits.compute_parametric_distance_along_segment_2_object()(entrySegment.start(), entrySegment.end(), current->window_left()) << std::endl;
-            std::cout << "Current Right Window: " << current->window_right() << "  ,  " << m_traits.compute_parametric_distance_along_segment_2_object()(entrySegment.start(), entrySegment.end(), current->window_right()) << std::endl;
+            std::cout << "Current Left Window: " << current->window_left() << "  ,  "
+                      << m_traits.compute_parametric_distance_along_segment_2_object()(entrySegment.start(), entrySegment.end(), current->window_left()) << std::endl;
+            std::cout << "Current Right Window: " << current->window_right() << "  ,  "
+                      << m_traits.compute_parametric_distance_along_segment_2_object()(entrySegment.start(), entrySegment.end(), current->window_right()) << std::endl;
             std::cout << "Current Segment Intersection: " << *result << std::endl;
-            std::cout << "Edge: (" << get(m_vertexIndexMap, source(current->entry_edge(), m_graph)) << "," << get(m_vertexIndexMap, target(current->entry_edge(), m_graph)) << ")  :  " << t0 << std::endl;
+            std::cout << "Edge: (" << get(m_vertexIndexMap, source(current->entry_edge(), m_graph))
+                      << "," << get(m_vertexIndexMap, target(current->entry_edge(), m_graph)) << ")  :  " << t0 << std::endl;
           }
 
           visitor(current->entry_edge(), t0);
@@ -1659,33 +1798,37 @@ private:
     }
   }
 
-  Point_2 face_location_with_normalized_coordinates(Cone_tree_node* node, Barycentric_coordinates location)
+  Point_2 face_location_with_normalized_coordinates(const Cone_tree_node* node,
+                                                    const Barycentric_coordinates& location) const
   {
     return construct_barycenter_in_triangle_2(node->layout_face(), localized_coordiate(node, location));
   }
 
-  Barycentric_coordinates localized_coordiate(Cone_tree_node* node, Barycentric_coordinates location)
+  Barycentric_coordinates localized_coordiate(const Cone_tree_node* node,
+                                              const Barycentric_coordinates& location) const
   {
     return shifted_coordinates(location, node->edge_face_index());
   }
 
-  Barycentric_coordinates shifted_coordinates(Barycentric_coordinates location, std::size_t shift)
+  Barycentric_coordinates shifted_coordinates(const Barycentric_coordinates& location,
+                                              const std::size_t shift) const
   {
     typename Traits::Construct_barycentric_coordinates_weight cbcw(m_traits.construct_barycentric_coordinates_weight_object());
     typename Traits::Construct_barycentric_coordinates cbc(m_traits.construct_barycentric_coordinates_object());
     return cbc(cbcw(location, shift), cbcw(location, (shift + 1) % 3), cbcw(location, (shift + 2) % 3));
   }
 
-  std::pair<Node_distance_pair, Barycentric_coordinates> nearest_on_face(face_descriptor f, Barycentric_coordinates location)
+  std::pair<Node_distance_pair, Barycentric_coordinates> nearest_on_face(const face_descriptor f,
+                                                                         const Barycentric_coordinates& location) const
   {
     typename Traits::Construct_barycentric_coordinates cbc(m_traits.construct_barycentric_coordinates_object());
 
-    std::size_t faceIndex = get(m_faceIndexMap, f);
+    const std::size_t faceIndex = get(m_faceIndexMap, f);
 
     Cone_tree_node* closest = nullptr;
     FT closestDistance = 0;
 
-    std::vector<Cone_tree_node*>& currentFaceList = m_faceOccupiers[faceIndex];
+    const std::vector<Cone_tree_node*>& currentFaceList = m_faceOccupiers[faceIndex];
 
     for (std::size_t i = 0; i < currentFaceList.size(); ++i)
     {
@@ -1696,7 +1839,7 @@ private:
         continue;
       }
 
-      Point_2 locationInContext = face_location_with_normalized_coordinates(current, location);
+      const Point_2 locationInContext = face_location_with_normalized_coordinates(current, location);
 
       if (current->inside_window(locationInContext))
       {
@@ -1716,11 +1859,12 @@ private:
     }
     else
     {
-      return std::make_pair(Node_distance_pair((Cone_tree_node*)nullptr, FT(0.0)), cbc(FT(0.0), FT(0.0), FT(0.0)));
+      return std::make_pair(Node_distance_pair((Cone_tree_node*)nullptr, FT(0)), cbc(FT(0), FT(0), FT(0)));
     }
   }
 
-  std::pair<Node_distance_pair, Barycentric_coordinates> nearest_to_location(face_descriptor f, Barycentric_coordinates location)
+  std::pair<Node_distance_pair, Barycentric_coordinates> nearest_to_location(const face_descriptor f,
+                                                                             const Barycentric_coordinates& location) const
   {
     typename Traits::Construct_barycentric_coordinates_weight cbcw(m_traits.construct_barycentric_coordinates_weight_object());
     typename Traits::Construct_barycentric_coordinates cbc(m_traits.construct_barycentric_coordinates_object());
@@ -1745,11 +1889,11 @@ private:
           std::pair<Node_distance_pair,Barycentric_coordinates> mainFace = nearest_on_face(f, location);
 
           halfedge_descriptor oppositeHalfedge = opposite(he, m_graph);
-          if(!CGAL::is_border(oppositeHalfedge, m_graph))
+          if (!CGAL::is_border(oppositeHalfedge, m_graph))
           {
-              std::size_t oppositeIndex = internal::edge_index(oppositeHalfedge, m_graph);
+              std::size_t oppositeIndex = Surface_mesh_shortest_paths_3::internal::edge_index(oppositeHalfedge, m_graph);
 
-              FT oppositeLocationCoords[3] = { FT(0.0), FT(0.0), FT(0.0) };
+              FT oppositeLocationCoords[3] = { FT(0), FT(0), FT(0) };
               oppositeLocationCoords[oppositeIndex] = cbcw(location, (associatedEdge + 1) % 3);
               oppositeLocationCoords[(oppositeIndex + 1) % 3] = cbcw(location, associatedEdge);
               Barycentric_coordinates oppositeLocation(cbc(oppositeLocationCoords[0], oppositeLocationCoords[1], oppositeLocationCoords[2]));
@@ -1783,7 +1927,7 @@ private:
 
           vertex_descriptor vertex = source(he, m_graph);
 
-          return std::make_pair(m_closestToVertices[get(m_vertexIndexMap, vertex)], cbc(FT(0.0), FT(0.0), FT(1.0)));
+          return std::make_pair(m_closestToVertices[get(m_vertexIndexMap, vertex)], cbc(FT(0), FT(0), FT(1)));
         }
         break;
 
@@ -1799,7 +1943,8 @@ private:
   }
 
   template <class InputIterator>
-  Source_point_iterator add_source_points_internal(InputIterator begin, InputIterator end, vertex_descriptor)
+  Source_point_iterator add_source_points_internal(InputIterator begin, InputIterator end,
+                                                   const vertex_descriptor)
   {
     Source_point_iterator firstAdded;
 
@@ -1817,7 +1962,8 @@ private:
   }
 
   template <class InputIterator>
-  Source_point_iterator add_source_points_internal(InputIterator begin, InputIterator end, Face_location)
+  Source_point_iterator add_source_points_internal(InputIterator begin, InputIterator end,
+                                                   const Face_location&)
   {
     Source_point_iterator firstAdded;
 
@@ -1898,7 +2044,8 @@ private:
     {
       if (m_debugOutput)
       {
-        std::cout << "Root: " << get(m_faceIndexMap, it->first) << " , " << it->second[0] << " " << it->second[1] << " " << it->second[2] << " " << std::endl;
+        std::cout << "Root: " << get(m_faceIndexMap, it->first)
+                  << " , " << it->second[0] << " " << it->second[1] << " " << it->second[2] << " " << std::endl;
       }
 
       expand_root(it->first, it->second, Source_point_iterator(it));
@@ -1915,6 +2062,39 @@ private:
 
     while (m_expansionPriqueue.size() > 0)
     {
+      if (m_debugOutput)
+      {
+        std::cout << " -----------------------------------------------------------------------" << std::endl;
+        std::cout << " -----------------------------------------------------------------------" << std::endl;
+        std::cout << "Num face locations: " << m_faceLocations.size() << std::endl;
+        std::cout << "Num root nodes: " << m_rootNodes.size() << " (Hint: these should be the same size)" << std::endl;
+
+        std::cout << "Prio Queue size = " << m_expansionPriqueue.size() << std::endl;
+        std::cout << "Queue:" << std::endl;
+        auto duplicate_queue = m_expansionPriqueue;
+        while(duplicate_queue.size() > 0)
+        {
+          Cone_expansion_event* event = duplicate_queue.top();
+
+          std::cout << "event type: " << event->m_type << " "
+                    << " time: " << event->m_distanceEstimate << " ";
+          std::cout << "cancelled? " << event->m_cancelled << " " ;
+
+          if (!event->m_cancelled)
+          {
+            std::cout << " ------ Parent (" << event->m_parent << ") INFO: ";
+            std::cout << "EH = (" << get(m_vertexIndexMap, source(event->m_parent->entry_edge(), m_graph)) << " "
+                                              << get(m_vertexIndexMap, target(event->m_parent->entry_edge(), m_graph)) << ") ";
+            std::cout << "S = (" << event->m_parent->source_image() << ") ";
+            std::cout << "T = " << get(m_vertexIndexMap, target(next(event->m_parent->entry_edge(), m_graph), m_graph));
+          }
+
+          std::cout << std::endl;
+
+          duplicate_queue.pop();
+        }
+      }
+
       Cone_expansion_event* event = m_expansionPriqueue.top();
       m_expansionPriqueue.pop();
 
@@ -2018,9 +2198,11 @@ public:
   /*!
   \brief Creates a shortest paths object using `tm` as input.
 
-  Equivalent to `Surface_mesh_shortest_path(tm, get(boost::vertex_index, tm), get(boost::halfedge_index, tm), get(boost::face_index, tm), get(CGAL::vertex_point, tm), traits)`.
+  Equivalent to `Surface_mesh_shortest_path(tm, get(boost::vertex_index, tm), get(boost::halfedge_index, tm),
+                                            get(boost::face_index, tm), get(CGAL::vertex_point, tm), traits)`.
   */
-  Surface_mesh_shortest_path(const Triangle_mesh& tm, const Traits& traits = Traits())
+  Surface_mesh_shortest_path(Triangle_mesh& tm,
+                             const Traits& traits = Traits())
     : m_traits(traits)
     , m_graph(const_cast<Triangle_mesh&>(tm))
     , m_vertexIndexMap(get(boost::vertex_index, tm))
@@ -2078,8 +2260,8 @@ public:
     if (m_debugOutput)
     {
       std::cout << "Final node count: " << m_currentNodeCount << std::endl;
+      std::cout << "Peak node count: " << m_peakNodeCount << std::endl;
     }
-    return;
 #endif
   }
 
@@ -2114,7 +2296,8 @@ public:
   \param location Barycentric coordinates in face `f` specifying the source point.
   \return An iterator to the source point added
   */
-  Source_point_iterator add_source_point(face_descriptor f, Barycentric_coordinates location)
+  Source_point_iterator add_source_point(const face_descriptor f,
+                                         const Barycentric_coordinates& location)
   {
     return add_source_point(std::make_pair(f, location));
   }
@@ -2123,7 +2306,7 @@ public:
   \brief Adds a point inside a face as a source for the shortest path queries,
   equivalent to `Surface_mesh_shortest_path::add_source_point(location.first, location.second);`
   */
-  Source_point_iterator add_source_point(Face_location location)
+  Source_point_iterator add_source_point(const Face_location& location)
   {
     Source_point_underlying_iterator added = m_faceLocations.insert(m_faceLocations.end(), location);
 
@@ -2142,7 +2325,8 @@ public:
   until either `Surface_mesh_shortest_path::build_sequence_tree()` or
   the first shortest path query is done.
 
-  \tparam InputIterator A `ForwardIterator` which dereferences to either `Surface_mesh_shortest_path::Face_location`, or `Surface_mesh_shortest_path::vertex_descriptor`.
+  \tparam InputIterator A `ForwardIterator` which dereferences to either `Surface_mesh_shortest_path::Face_location`,
+                        or `Surface_mesh_shortest_path::vertex_descriptor`.
 
   \param begin iterator to the first in the list of source point locations.
   \param end iterator to one past the end of the list of source point locations.
@@ -2184,7 +2368,8 @@ public:
   */
   void remove_all_source_points()
   {
-    m_deletedSourceLocations.splice(m_deletedSourceLocations.begin(), m_faceLocations, m_faceLocations.begin(), m_faceLocations.end());
+    m_deletedSourceLocations.splice(m_deletedSourceLocations.begin(), m_faceLocations,
+                                    m_faceLocations.begin(), m_faceLocations.end());
     m_firstNewSourcePoint = m_faceLocations.end();
   }
 
@@ -2284,7 +2469,7 @@ public:
     occur when the graph is disconnected), the distance will be a negative
     value and the source point iterator will be equal to `source_points_end()`.
   */
-  Shortest_path_result shortest_distance_to_source_points(vertex_descriptor v)
+  Shortest_path_result shortest_distance_to_source_points(const vertex_descriptor v)
   {
     build_sequence_tree();
 
@@ -2298,7 +2483,7 @@ public:
     }
     else
     {
-      return std::make_pair(FT(-1.0), source_points_end());
+      return std::make_pair(FT(-1), source_points_end());
     }
   }
 
@@ -2312,7 +2497,8 @@ public:
     occur when the graph is disconnected), the distance will be a negative
     value and the source point iterator will be equal to `source_points_end()`.
   */
-  Shortest_path_result shortest_distance_to_source_points(face_descriptor f, Barycentric_coordinates location)
+  Shortest_path_result shortest_distance_to_source_points(const face_descriptor f,
+                                                          const Barycentric_coordinates& location)
   {
     build_sequence_tree();
 
@@ -2326,7 +2512,7 @@ public:
     }
     else
     {
-      return std::make_pair(FT(-1.0), source_points_end());
+      return std::make_pair(FT(-1), source_points_end());
     }
   }
 
@@ -2353,7 +2539,8 @@ public:
   */
   template <class Visitor>
   Shortest_path_result
-  shortest_path_sequence_to_source_points(vertex_descriptor v, Visitor& visitor)
+  shortest_path_sequence_to_source_points(const vertex_descriptor v,
+                                          Visitor& visitor)
   {
     build_sequence_tree();
 
@@ -2368,7 +2555,7 @@ public:
     }
     else
     {
-      return std::make_pair(FT(-1.0), source_points_end());
+      return std::make_pair(FT(-1), source_points_end());
     }
   }
 
@@ -2391,7 +2578,9 @@ public:
   */
   template <class Visitor>
   Shortest_path_result
-  shortest_path_sequence_to_source_points(face_descriptor f, Barycentric_coordinates location, Visitor& visitor)
+  shortest_path_sequence_to_source_points(const face_descriptor f,
+                                          const Barycentric_coordinates& location,
+                                          Visitor& visitor)
   {
     build_sequence_tree();
 
@@ -2407,7 +2596,7 @@ public:
     }
     else
     {
-      return std::make_pair(FT(-1.0), source_points_end());
+      return std::make_pair(FT(-1), source_points_end());
     }
   }
 
@@ -2430,7 +2619,7 @@ public:
   */
   template <class OutputIterator>
   Shortest_path_result
-  shortest_path_points_to_source_points(vertex_descriptor v, OutputIterator output)
+  shortest_path_points_to_source_points(const vertex_descriptor v, OutputIterator output)
   {
     build_sequence_tree();
 
@@ -2453,7 +2642,8 @@ public:
   */
   template <class OutputIterator>
   Shortest_path_result
-  shortest_path_points_to_source_points(face_descriptor f, Barycentric_coordinates location, OutputIterator output)
+  shortest_path_points_to_source_points(const face_descriptor f,
+                                        const Barycentric_coordinates& location, OutputIterator output)
   {
     build_sequence_tree();
 
@@ -2471,27 +2661,33 @@ public:
     of the given face.
 
   \details The following static overloads are also available:
-    - `static Point_3 point(face_descriptor f, Barycentric_coordinates location, const Triangle_mesh& tm, const Traits& traits = Traits())`
-    - `static Point_3 point(face_descriptor f, Barycentric_coordinates location, const Triangle_mesh& tm, Vertex_point_map vertexPointMap, const Traits& traits = Traits())`
+    - `static Point_3 point(face_descriptor f, Barycentric_coordinates location, const Triangle_mesh& tm,
+                            const Traits& traits = Traits())`
+    - `static Point_3 point(face_descriptor f, Barycentric_coordinates location, const Triangle_mesh& tm,
+                            Vertex_point_map vertexPointMap, const Traits& traits = Traits())`
 
   \param f A face of on the input face graph
   \param location The barycentric coordinates of the query point on face `f`
   */
-  Point_3 point(face_descriptor f, Barycentric_coordinates location) const
+  Point_3 point(const face_descriptor f,
+                const Barycentric_coordinates& location) const
   {
     return point(f, location, m_graph, m_vertexPointMap, m_traits);
   }
 
   /// \cond
 
-  static Point_3 point(face_descriptor f, Barycentric_coordinates location, const Triangle_mesh& tm, const Traits& traits = Traits())
+  static Point_3 point(const face_descriptor f,
+                       const Barycentric_coordinates& location,
+                       const Triangle_mesh& tm,
+                       const Traits& traits = Traits())
   {
     using boost::get;
     return point(f, location, tm, get(CGAL::vertex_point, tm), traits);
   }
 
-  static Point_3 point(face_descriptor f,
-                       Barycentric_coordinates location,
+  static Point_3 point(const face_descriptor f,
+                       const Barycentric_coordinates& location,
                        const Triangle_mesh& tm,
                        Vertex_point_map vertexPointMap,
                        const Traits& traits = Traits())
@@ -2512,20 +2708,25 @@ public:
   \param edge An edge of the input face graph
   \param t The parametric distance along edge of the desired point
   */
-  Point_3 point(halfedge_descriptor edge, FT t) const
+  Point_3 point(const halfedge_descriptor edge, const FT t) const
   {
     return point(edge, t, m_graph, m_vertexPointMap, m_traits);
   }
 
   /// \cond
 
-  static Point_3 point(halfedge_descriptor edge, FT t, const Triangle_mesh& tm, const Traits& traits = Traits())
+  static Point_3 point(const halfedge_descriptor edge, const FT t,
+                       const Triangle_mesh& tm,
+                       const Traits& traits = Traits())
   {
     using boost::get;
     return point(edge, t, tm, get(CGAL::vertex_point, tm), traits);
   }
 
-  static Point_3 point(halfedge_descriptor edge, FT t, const Triangle_mesh& tm, Vertex_point_map vertexPointMap, const Traits& traits = Traits())
+  static Point_3 point(const halfedge_descriptor edge, const FT t,
+                       const Triangle_mesh& tm,
+                       Vertex_point_map vertexPointMap,
+                       const Traits& traits = Traits())
   {
     typename Traits::Construct_barycenter_3 construct_barycenter_3(traits.construct_barycenter_3_object());
 
@@ -2540,7 +2741,7 @@ public:
 
   \param vertex A vertex of the input face graph
   */
-  Point_3 point(vertex_descriptor vertex) const
+  Point_3 point(const vertex_descriptor vertex) const
   {
     return get(m_vertexPointMap, vertex);
   }
@@ -2558,23 +2759,29 @@ public:
 
   \param vertex A vertex of the input face graph
   */
-  Face_location face_location(vertex_descriptor vertex) const
+  Face_location face_location(const vertex_descriptor vertex) const
   {
     return face_location(vertex, m_graph, m_traits);
   }
 
   /// \cond
 
-  static Face_location face_location(vertex_descriptor vertex, const Triangle_mesh& tm, const Traits& traits = Traits())
+  static Face_location face_location(const vertex_descriptor vertex,
+                                     const Triangle_mesh& tm,
+                                     const Traits& traits = Traits())
   {
     typename Traits::Construct_barycentric_coordinates construct_barycentric_coordinates(traits.construct_barycentric_coordinates_object());
-    halfedge_descriptor he = next(halfedge(vertex, tm), tm);
+    halfedge_descriptor hinit=halfedge(vertex, tm);
+    while (face(hinit, tm) == Graph_traits::null_face())
+      hinit = opposite(next(hinit, tm), tm);
+
+    halfedge_descriptor he = next(hinit, tm);
     face_descriptor locationFace = face(he, tm);
-    std::size_t edgeIndex = CGAL::internal::edge_index(he, tm);
+    std::size_t edgeIndex = Surface_mesh_shortest_paths_3::internal::edge_index(he, tm);
 
-    FT coords[3] = { FT(0.0), FT(0.0), FT(0.0) };
+    FT coords[3] = { FT(0), FT(0), FT(0) };
 
-    coords[edgeIndex] = FT(1.0);
+    coords[edgeIndex] = FT(1);
 
     return Face_location(locationFace, construct_barycentric_coordinates(coords[0], coords[1], coords[2]));
   }
@@ -2590,22 +2797,24 @@ public:
   \param he A halfedge of the input face graph
   \param t Parametric distance of the desired point along `he`
   */
-  Face_location face_location(halfedge_descriptor he, FT t) const
+  Face_location face_location(const halfedge_descriptor he, const FT t) const
   {
     return face_location(he, t, m_graph, m_traits);
   }
 
   /// \cond
 
-  static Face_location face_location(halfedge_descriptor he, FT t, const Triangle_mesh& tm, const Traits& traits = Traits())
+  static Face_location face_location(const halfedge_descriptor he, FT t,
+                                     const Triangle_mesh& tm,
+                                     const Traits& traits = Traits())
   {
     typename Traits::Construct_barycentric_coordinates cbc(traits.construct_barycentric_coordinates_object());
     face_descriptor locationFace = face(he, tm);
-    std::size_t edgeIndex = CGAL::internal::edge_index(he, tm);
+    std::size_t edgeIndex = Surface_mesh_shortest_paths_3::internal::edge_index(he, tm);
 
-    const FT oneMinusT(FT(1.0) - t);
+    const FT oneMinusT(FT(1) - t);
 
-    FT coords[3] = { FT(0.0), FT(0.0), FT(0.0) };
+    FT coords[3] = { FT(0), FT(0), FT(0) };
 
     coords[edgeIndex] = oneMinusT;
     coords[(edgeIndex + 1) % 3] = t;
@@ -2799,7 +3008,7 @@ public:
     }
     else
     {
-      return Face_location(Graph_traits::null_face(), cbc(FT(0.0), FT(0.0), FT(0.0)));
+      return Face_location(Graph_traits::null_face(), cbc(FT(0), FT(0), FT(0)));
     }
   }
 
@@ -2831,7 +3040,8 @@ public:
   }
 
   template <class AABBTraits>
-  static void build_aabb_tree(const Triangle_mesh& tm, AABB_tree<AABBTraits>& outTree,
+  static void build_aabb_tree(const Triangle_mesh& tm,
+                              AABB_tree<AABBTraits>& outTree,
                               Vertex_point_map vertexPointMap)
   {
     face_iterator facesStart, facesEnd;
diff --git a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/barycentric.h b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/barycentric.h
index 4c9d0f60a45..37b422e5421 100644
--- a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/barycentric.h
+++ b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/barycentric.h
@@ -58,7 +58,9 @@ public:
   {
   }
 
-  Construct_barycentric_coordinates_in_triangle_2(const Construct_barycentric_coordinates& cbc, const Construct_vector_2& cv2, const Compute_scalar_product_2& csp2)
+  Construct_barycentric_coordinates_in_triangle_2(const Construct_barycentric_coordinates& cbc,
+                                                  const Construct_vector_2& cv2,
+                                                  const Compute_scalar_product_2& csp2)
     : m_construct_barycentric_coordinates(cbc)
     , m_construct_vector_2(cv2)
     , m_compute_scalar_product_2(csp2)
@@ -81,7 +83,7 @@ public:
 
     FT v = (d11 * d20 - d01 * d21) / denom;
     FT w = (d00 * d21 - d01 * d20) / denom;
-    return m_construct_barycentric_coordinates(FT(1.0) - v - w, v, w);
+    return m_construct_barycentric_coordinates(FT(1) - v - w, v, w);
   }
 };
 
@@ -109,7 +111,9 @@ public:
   {
   }
 
-  Construct_barycentric_coordinates_in_triangle_3(const Construct_barycentric_coordinates& cbc, const Construct_vector_3& cv3, const Compute_scalar_product_3& csp3)
+  Construct_barycentric_coordinates_in_triangle_3(const Construct_barycentric_coordinates& cbc,
+                                                  const Construct_vector_3& cv3,
+                                                  const Compute_scalar_product_3& csp3)
     : m_construct_barycentric_coordinates(cbc)
     , m_construct_vector_3(cv3)
     , m_compute_scalar_product_3(csp3)
@@ -132,7 +136,7 @@ public:
 
     FT v = (d11 * d20 - d01 * d21) / denom;
     FT w = (d00 * d21 - d01 * d20) / denom;
-    return m_construct_barycentric_coordinates(FT(1.0) - v - w, v, w);
+    return m_construct_barycentric_coordinates(FT(1) - v - w, v, w);
   }
 };
 
@@ -194,7 +198,8 @@ public:
     bool nonZero[3];
     std::size_t numNonZero = 0;
 
-    if (cbcw(baryCoords, 0) + cbcw(baryCoords, 1) + cbcw(baryCoords, 2) > 1.00001 || cbcw(baryCoords, 0) + cbcw(baryCoords, 1) + cbcw(baryCoords, 2) < 0.99999)
+    if (cbcw(baryCoords, 0) + cbcw(baryCoords, 1) + cbcw(baryCoords, 2) > 1.00001 ||
+        cbcw(baryCoords, 0) + cbcw(baryCoords, 1) + cbcw(baryCoords, 2) < 0.99999)
     {
       return std::make_pair(BARYCENTRIC_COORDINATES_ON_UNBOUNDED_SIDE, 0);
     }
diff --git a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/function_objects.h b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/function_objects.h
index e08623c2ce7..de838369d8b 100644
--- a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/function_objects.h
+++ b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/function_objects.h
@@ -18,11 +18,10 @@
 //
 // Author(s)     : Stephen Kiazyk
 
-#include <cstddef>
+#ifndef CGAL_SURFACE_MESH_SHORTEST_PATH_INTERNAL_FUNCTION_OBJECTS_H
+#define CGAL_SURFACE_MESH_SHORTEST_PATH_INTERNAL_FUNCTION_OBJECTS_H
 
-#include <boost/graph/graph_traits.hpp>
-#include <CGAL/result_of.h>
-#include <CGAL/assertions.h>
+#include <CGAL/license/Surface_mesh_shortest_path.h>
 
 #include <CGAL/Surface_mesh_shortest_path/internal/misc_functions.h>
 
@@ -32,13 +31,15 @@
 #endif
 #endif
 
+#include <CGAL/assertions.h>
+#include <CGAL/boost/graph/helpers.h>
+#include <CGAL/number_utils.h>
+#include <CGAL/result_of.h>
 #include <CGAL/Cartesian_converter.h>
 
-#ifndef CGAL_SURFACE_MESH_SHORTEST_PATH_INTERNAL_FUNCTION_OBJECTS_H
-#define CGAL_SURFACE_MESH_SHORTEST_PATH_INTERNAL_FUNCTION_OBJECTS_H
-
-#include <CGAL/license/Surface_mesh_shortest_path.h>
-
+#include <cmath>
+#include <cstddef>
+#include <limits>
 
 namespace CGAL {
 
@@ -246,8 +247,8 @@ public:
 
     Point_3 projectedLocation3d(m_construct_projected_point_3(baseSegment, m_construct_vertex_3(t3, 2)));
     FT scalePoint = m_parametric_distance_along_segment_3(m_construct_vertex_3(t3, 0), m_construct_vertex_3(t3, 1), projectedLocation3d);
-    FT triangleHeight = CGAL::internal::select_sqrt(m_compute_squared_distance_3(projectedLocation3d, t3[2]));
-    FT v01Len = CGAL::internal::select_sqrt(m_compute_squared_distance_3(t3[1], t3[0]));
+    FT triangleHeight = CGAL::approximate_sqrt(m_compute_squared_distance_3(projectedLocation3d, t3[2]));
+    FT v01Len = CGAL::approximate_sqrt(m_compute_squared_distance_3(t3[1], t3[0]));
 
     Point_2 A(m_construct_point_2(0.0, 0.0));
     Point_2 B(m_construct_point_2(v01Len, 0.0));
@@ -367,12 +368,12 @@ public:
   {
     Point_3 projectedLocation3d(m_construct_projected_point_3(m_construct_line_3(m_construct_vertex_3(t3, edgeIndex), m_construct_vertex_3(t3, edgeIndex + 1)), m_construct_vertex_3(t3, edgeIndex + 2)));
     FT scalePoint = m_parametric_distance_along_segment_3(m_construct_segment_3(m_construct_vertex_3(t3, edgeIndex), m_construct_vertex_3(t3, edgeIndex + 1)), projectedLocation3d);
-    FT triangleHeight = CGAL::internal::select_sqrt(m_compute_squared_distance_3(projectedLocation3d, m_construct_vertex_3(t3, edgeIndex + 2)));
+    FT triangleHeight = CGAL::approximate_sqrt(m_compute_squared_distance_3(projectedLocation3d, m_construct_vertex_3(t3, edgeIndex + 2)));
 
     Vector_2 edgeVector(m_construct_vector_2(segment));
 
     Vector_2 perpendicularEdgeVector(m_construct_perpendicular_vector_2(edgeVector, CGAL::COUNTERCLOCKWISE));
-    perpendicularEdgeVector = m_construct_scaled_vector_2(perpendicularEdgeVector, FT(1.0) / CGAL::internal::select_sqrt(m_compute_squared_length_2(perpendicularEdgeVector)));
+    perpendicularEdgeVector = m_construct_scaled_vector_2(perpendicularEdgeVector, FT(1) / CGAL::approximate_sqrt(m_compute_squared_length_2(perpendicularEdgeVector)));
 
     Point_2 points[3];
     points[edgeIndex] = m_construct_source_2(segment);
@@ -433,16 +434,19 @@ public:
 
   typedef typename K::Intersect_2 Intersect_2;
   typedef typename K::Compare_distance_2 Compare_distance_2;
+  typedef typename K::Compute_squared_distance_2 Compute_squared_distance_2;
   typedef typename K::Construct_line_2 Construct_line_2;
   typedef typename K::Construct_source_2 Construct_source_2;
   typedef typename K::Construct_target_2 Construct_target_2;
 
+
   typedef CGAL::Comparison_result result_type;
 
 private:
   Compute_parametric_distance_along_segment_2<K> m_parametric_distance_along_segment_2;
   Intersect_2 m_intersect_2;
   Compare_distance_2 m_compare_distance_2;
+  Compute_squared_distance_2 m_compute_squared_distance_2;
   Construct_line_2 m_construct_line_2;
   Construct_source_2 m_construct_source_2;
   Construct_target_2 m_construct_target_2;
@@ -455,6 +459,7 @@ public:
   Compare_relative_intersection_along_segment_2(const K& kernel)
     : m_intersect_2(kernel.intersect_2_object())
     , m_compare_distance_2(kernel.compare_distance_2_object())
+    , m_compute_squared_distance_2(kernel.compute_squared_distance_2_object())
     , m_construct_line_2(kernel.construct_line_2_object())
     , m_construct_source_2(kernel.construct_source_2_object())
     , m_construct_target_2(kernel.construct_target_2_object())
@@ -466,7 +471,6 @@ public:
     typedef typename CGAL::cpp11::result_of<Intersect_2(Line_2, Line_2)>::type LineLineIntersectResult;
 
     Line_2 s1Line(m_construct_line_2(s1));
-
     Line_2 s2Line(m_construct_line_2(s2));
 
     LineLineIntersectResult intersectResult1(m_intersect_2(s1Line, l1));
@@ -479,7 +483,7 @@ public:
     if (!p1_ptr) return CGAL::SMALLER;
 
     CGAL_assertion_code(FT t1 = m_parametric_distance_along_segment_2(s1, *p1_ptr);)
-    CGAL_assertion(t1 >= FT(-0.00001) && t1 <= FT(1.00001));
+    CGAL_assertion(t1 >= FT(-1)/FT(100000) && t1 <= FT(1)+FT(1)/FT(100000));
 
     LineLineIntersectResult intersectResult2 = m_intersect_2(s2Line, l2);
     CGAL_assertion(bool(intersectResult2));
@@ -491,9 +495,27 @@ public:
     if (!p2_ptr) return CGAL::SMALLER;
 
     CGAL_assertion_code(FT t2 = m_parametric_distance_along_segment_2(s2, *p2_ptr);)
-    CGAL_assertion(t2 >= FT(-0.00001) && t2 <= FT(1.00001));
+    CGAL_assertion(t2 >= FT(-1)/FT(100000) && t2 <= FT(1)+FT(1)/FT(100000));
 
+// #define CGAL_SMSP_DONT_USE_RELAXED_PRUNING
+#ifndef CGAL_SMSP_DONT_USE_RELAXED_PRUNING
+    const FT sqd_1 = m_compute_squared_distance_2(s1.source(), *p1_ptr);
+    const FT sqd_2 = m_compute_squared_distance_2(s2.source(), *p2_ptr);
+
+    // In the case of multiple rays reaching the same target, we want to know their respective position
+    // so that pruning of branches can be done according to the "one angle one split" idiom.
+    // However, the orientation predicate is evaluated in the unfolded 2D plane, which is obtained
+    // via square roots; inconsisnties will exist. We don't want to prune in case it might be wrong,
+    // so we add a little bit of tolerance on the evaluation of the predicate. If it's almost collinear,
+    // return 'collinear' (EQUAL).
+    const FT eps = (FT(100) * std::numeric_limits<FT>::epsilon());
+    if(CGAL::abs(sqd_1 - sqd_2) < eps)
+      return CGAL::EQUAL;
+
+    return CGAL::compare(sqd_1, sqd_2);
+#else
     return m_compare_distance_2(s1.source(), *p1_ptr, s2.source(), *p2_ptr);
+#endif
   }
 };
 
@@ -501,13 +523,14 @@ template <class Kernel, class FaceListGraph>
 class Is_saddle_vertex
 {
 public:
+  typedef typename Kernel::FT FT;
+  typedef typename Kernel::Point_2 Point_2;
   typedef typename Kernel::Point_3 Point_3;
+  typedef typename Kernel::Segment_2 Segment_2;
+  typedef typename Kernel::Vector_2 Vector_2;
   typedef typename Kernel::Vector_3 Vector_3;
   typedef typename Kernel::Triangle_3 Triangle_3;
   typedef typename Kernel::Triangle_2 Triangle_2;
-  typedef typename Kernel::Segment_2 Segment_2;
-  typedef typename Kernel::Vector_2 Vector_2;
-  typedef typename Kernel::Point_2 Point_2;
 
   typedef typename boost::graph_traits<FaceListGraph> Graph_traits;
   typedef typename Graph_traits::vertex_descriptor vertex_descriptor;
@@ -540,7 +563,9 @@ public:
   {
   }
 
-  Is_saddle_vertex(const Kernel& kernel, const Construct_triangle_3_to_triangle_2_projection& pt3tt2, const Construct_triangle_3_along_segment_2_flattening& ft3as2)
+  Is_saddle_vertex(const Kernel& kernel,
+                   const Construct_triangle_3_to_triangle_2_projection& pt3tt2,
+                   const Construct_triangle_3_along_segment_2_flattening& ft3as2)
     : m_orientation_2(kernel.orientation_2_object())
     , m_construct_triangle_3(kernel.construct_triangle_3_object())
     , m_construct_vertex_2(kernel.construct_vertex_2_object())
@@ -557,10 +582,72 @@ public:
     return (*this)(v, g, get(boost::vertex_point, g));
   }
 
+ FT angle(const Vector_3& u, const Vector_3& v) const
+ {
+   typename Kernel::Compute_scalar_product_3 scalar_product;
+
+   double product = CGAL::sqrt(to_double(scalar_product(u,u)) * to_double(scalar_product(v,v)));
+
+   if(product == 0)
+     return 0;
+
+   // cosine
+   double dot = to_double(scalar_product(u,v));
+   double cosine = dot / product;
+
+   if(cosine > 1.)
+     cosine = 1.;
+
+   if(cosine < -1.)
+     cosine = -1.;
+
+   return std::acos(cosine) * 180./CGAL_PI;
+ }
+
+
+ FT angle(const Point_3& p, const Point_3& q, const Point_3& r) const
+ {
+   typename Kernel::Construct_vector_3 cv;
+
+   Vector_3 u = cv(q, p);
+   Vector_3 v = cv(q, r);
+
+   return angle(u, v);
+ }
+
+  template<class VertexPointMap>
+  FT vertex_angle(const vertex_descriptor v,
+                  const FaceListGraph& g,
+                  const VertexPointMap& pointMap) const
+  {
+    FT angle_sum = 0;
+
+    for(halfedge_descriptor h : halfedges_around_target(v, g))
+    {
+      if(is_border(h, g))
+        continue;
+
+      angle_sum += angle(get(pointMap, source(h, g)),
+                         get(pointMap, target(h, g)),
+                         get(pointMap, target(next(h, g), g)));
+    }
+
+    angle_sum *= CGAL_PI / FT(180);
+
+    return angle_sum;
+  }
+
   template<class VertexPointMap>
   result_type operator() (vertex_descriptor v, const FaceListGraph& g, VertexPointMap const& pointMap) const
   {
+#ifndef CGAL_SMSP_DONT_USE_RELAXED_PRUNING
+    const FT ang_sum = vertex_angle(v, g, pointMap);
+    const FT bound = (FT(1) - FT(100) * std::numeric_limits<FT>::epsilon()) * 2 * CGAL_PI;
+    return (ang_sum >= bound);
+#else
     halfedge_descriptor startEdge = halfedge(v, g);
+    while (face(startEdge, g) == Graph_traits::null_face())
+      startEdge=opposite(next(startEdge, g), g);
 
     Point_3 rootPoint(get(pointMap, v));
     Point_3 prevPoint(get(pointMap, source(startEdge, g)));
@@ -607,6 +694,7 @@ public:
     while (currentEdge != startEdge);
 
     return false;
+#endif
   }
 };
 
diff --git a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/Cone_expansion_event.h b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/Cone_expansion_event.h
index fc44dcd9826..43c8394968d 100644
--- a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/Cone_expansion_event.h
+++ b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/Cone_expansion_event.h
@@ -23,9 +23,8 @@
 
 #include <CGAL/license/Surface_mesh_shortest_path.h>
 
-
 namespace CGAL {
-
+namespace Surface_mesh_shortest_paths_3 {
 namespace internal {
 
 template<class Traits>
@@ -54,7 +53,9 @@ public:
   bool m_cancelled;
 
 public:
-  Cone_expansion_event(Cone_tree_node<Traits>* parent, const FT& distanceEstimate, Expansion_type type)
+  Cone_expansion_event(Cone_tree_node<Traits>* parent,
+                       const FT& distanceEstimate,
+                       Expansion_type type)
     : m_parent(parent)
     , m_distanceEstimate(distanceEstimate)
     , m_type(type)
@@ -62,7 +63,10 @@ public:
   {
   }
 
-  Cone_expansion_event(Cone_tree_node<Traits>* parent, const FT& distanceEstimate, Expansion_type type, const Segment_2& windowSegment)
+  Cone_expansion_event(Cone_tree_node<Traits>* parent,
+                       const FT& distanceEstimate,
+                       Expansion_type type,
+                       const Segment_2& windowSegment)
     : m_parent(parent)
     , m_distanceEstimate(distanceEstimate)
     , m_type(type)
@@ -77,14 +81,15 @@ template <class Traits>
 struct Cone_expansion_event_min_priority_queue_comparator
 {
 public:
-  bool operator () (const Cone_expansion_event<Traits>* lhs, const Cone_expansion_event<Traits>* rhs) const
+  bool operator () (const Cone_expansion_event<Traits>* lhs,
+                    const Cone_expansion_event<Traits>* rhs) const
   {
     return rhs->m_distanceEstimate < lhs->m_distanceEstimate;
   }
 };
 
 } // namespace internal
-
+} // namespace Surface_mesh_shortest_paths_3
 } // namespace CGAL
 
 #endif // CGAL_SURFACE_MESH_SHORTEST_PATH_INTERNAL_CONE_EXPANSION_EVENT_H
diff --git a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/Cone_tree.h b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/Cone_tree.h
index 0915c896775..7534017559b 100644
--- a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/Cone_tree.h
+++ b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/Cone_tree.h
@@ -23,7 +23,6 @@
 
 #include <CGAL/license/Surface_mesh_shortest_path.h>
 
-
 #include <vector>
 
 #include <boost/graph/graph_traits.hpp>
@@ -32,11 +31,11 @@
 #include <CGAL/Surface_mesh_shortest_path/internal/Cone_expansion_event.h>
 #include <CGAL/Surface_mesh_shortest_path/internal/misc_functions.h>
 
-namespace CGAL
-{
+#include <CGAL/number_utils.h>
 
-namespace internal
-{
+namespace CGAL {
+namespace Surface_mesh_shortest_paths_3 {
+namespace internal {
 
 template<class Traits>
 class Cone_tree_node
@@ -62,33 +61,38 @@ private:
   typedef typename Graph_traits::face_descriptor face_descriptor;
   typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor;
   typedef typename Graph_traits::vertex_descriptor vertex_descriptor;
-  typedef typename CGAL::internal::Cone_expansion_event<Traits> Cone_expansion_event;
+  typedef typename Surface_mesh_shortest_paths_3::internal::Cone_expansion_event<Traits> Cone_expansion_event;
 
 private:
   // These could be pulled back into a 'context' class to save space
-  Traits& m_traits;
-  Triangle_mesh& m_graph;
+  const Traits& m_traits;
+  const Triangle_mesh& m_graph;
 
-  halfedge_descriptor m_entryEdge;
+  const halfedge_descriptor m_entryEdge;
 
-  Point_2 m_sourceImage;
-  Triangle_2 m_layoutFace;
-  FT m_pseudoSourceDistance;
+  const Point_2 m_sourceImage;
+  const Triangle_2 m_layoutFace;
+  const FT m_pseudoSourceDistance;
 
-  Point_2 m_windowLeft;
-  Point_2 m_windowRight;
+  const Point_2 m_windowLeft;
+  const Point_2 m_windowRight;
 
-  size_t m_level;
-  size_t m_treeId;
+  std::size_t m_level;
+  std::size_t m_treeId;
 
-  Node_type m_nodeType;
+  const Node_type m_nodeType;
 
   Cone_tree_node* m_leftChild;
+  std::vector<Cone_tree_node*> m_middleChildren;
   Cone_tree_node* m_rightChild;
 
-  std::vector<Cone_tree_node*> m_middleChildren;
   Cone_tree_node* m_parent;
 
+public:
+  Cone_expansion_event* m_pendingLeftSubtree;
+  Cone_expansion_event* m_pendingRightSubtree;
+  Cone_expansion_event* m_pendingMiddleSubtree;
+
 private:
   void on_child_link(Cone_tree_node* child)
   {
@@ -98,7 +102,9 @@ private:
   }
 
 public:
-  Cone_tree_node(Traits& traits, Triangle_mesh& g, size_t treeId)
+  Cone_tree_node(const Traits& traits,
+                 const Triangle_mesh& g,
+                 const std::size_t treeId)
     : m_traits(traits)
     , m_graph(g)
     , m_sourceImage(Point_2(CGAL::ORIGIN))
@@ -115,7 +121,10 @@ public:
   {
   }
 
-  Cone_tree_node(Traits& traits, Triangle_mesh& g, size_t treeId, halfedge_descriptor entryEdge)
+  Cone_tree_node(const Traits& traits,
+                 const Triangle_mesh& g,
+                 const std::size_t treeId,
+                 const halfedge_descriptor entryEdge)
     : m_traits(traits)
     , m_graph(g)
     , m_entryEdge(entryEdge)
@@ -133,7 +142,15 @@ public:
   {
   }
 
-  Cone_tree_node(Traits& traits, Triangle_mesh& g, halfedge_descriptor entryEdge, const Triangle_2& layoutFace, const Point_2& sourceImage, const FT& pseudoSourceDistance, const Point_2& windowLeft, const Point_2& windowRight, Node_type nodeType = INTERVAL)
+  Cone_tree_node(const Traits& traits,
+                 const Triangle_mesh& g,
+                 const halfedge_descriptor entryEdge,
+                 const Triangle_2& layoutFace,
+                 const Point_2& sourceImage,
+                 const FT& pseudoSourceDistance,
+                 const Point_2& windowLeft,
+                 const Point_2& windowRight,
+                 const Node_type nodeType = INTERVAL)
     : m_traits(traits)
     , m_graph(g)
     , m_entryEdge(entryEdge)
@@ -151,12 +168,12 @@ public:
   {
   }
 
-  size_t tree_id() const
+  std::size_t tree_id() const
   {
     return m_treeId;
   }
 
-  size_t level() const
+  std::size_t level() const
   {
     return m_level;
   }
@@ -176,7 +193,7 @@ public:
     return m_nodeType == ROOT;
   }
 
-  Triangle_2 layout_face() const
+  const Triangle_2& layout_face() const
   {
     return m_layoutFace;
   }
@@ -191,9 +208,9 @@ public:
     return current_face() == Graph_traits::null_face();
   }
 
-  size_t edge_face_index() const
+  std::size_t edge_face_index() const
   {
-    return CGAL::internal::edge_index(entry_edge(), m_graph);
+    return edge_index(entry_edge(), m_graph);
   }
 
   halfedge_descriptor entry_edge() const
@@ -216,7 +233,7 @@ public:
     return target(next(m_entryEdge, m_graph), m_graph);
   }
 
-  Point_2 source_image() const
+  const Point_2& source_image() const
   {
     return m_sourceImage;
   }
@@ -229,7 +246,7 @@ public:
   FT distance_to_root(const Point_2& point) const
   {
     typename Traits::Compute_squared_distance_2 csd2(m_traits.compute_squared_distance_2_object());
-    return CGAL::internal::select_sqrt(csd2(point, m_sourceImage)) + m_pseudoSourceDistance;
+    return CGAL::approximate_sqrt(csd2(point, m_sourceImage)) + m_pseudoSourceDistance;
   }
 
   FT distance_from_source_to_root() const
@@ -252,12 +269,12 @@ public:
     return Ray_2(source_image(), m_windowRight);
   }
 
-  Point_2 window_left() const
+  const Point_2& window_left() const
   {
     return m_windowLeft;
   }
 
-  Point_2 window_right() const
+  const Point_2& window_right() const
   {
     return m_windowRight;
   }
@@ -270,10 +287,13 @@ public:
   bool inside_window(const Point_2& point) const
   {
     typename Traits::Orientation_2 orientation_2(m_traits.orientation_2_object());
+
     Point_2 sourceImagePoint(source_image());
     CGAL::Orientation leftOrientation = orientation_2(sourceImagePoint, m_windowLeft, point);
     CGAL::Orientation rightOrientation = orientation_2(sourceImagePoint, m_windowRight, point);
-    return (leftOrientation == CGAL::RIGHT_TURN || leftOrientation == CGAL::COLLINEAR) && (rightOrientation == CGAL::LEFT_TURN || rightOrientation == CGAL::COLLINEAR);
+
+    return (leftOrientation == CGAL::RIGHT_TURN || leftOrientation == CGAL::COLLINEAR) &&
+           (rightOrientation == CGAL::LEFT_TURN || rightOrientation == CGAL::COLLINEAR);
   }
 
   Point_2 target_point() const
@@ -289,25 +309,17 @@ public:
 
   bool has_left_side() const
   {
-    typename Traits::Orientation_2 orientation_2(m_traits.orientation_2_object());
-
     if (is_source_node())
     {
       return true;
     }
 
-    CGAL::Orientation orientation = orientation_2(source_image(), m_windowLeft, target_point());
-
-    return (orientation == CGAL::RIGHT_TURN || orientation == CGAL::COLLINEAR);
+    return (m_traits.orientation_2_object()(source_image(), m_windowLeft, target_point()) != CGAL::LEFT_TURN);
   }
 
   bool has_right_side() const
   {
-    typename Traits::Orientation_2 orientation_2(m_traits.orientation_2_object());
-
-    CGAL::Orientation orientation = orientation_2(source_image(), m_windowRight, target_point());
-
-    return (orientation == CGAL::LEFT_TURN || orientation == CGAL::COLLINEAR);
+    return (m_traits.orientation_2_object()(source_image(), m_windowRight, target_point()) != CGAL::RIGHT_TURN);
   }
 
   Segment_2 left_child_base_segment() const
@@ -335,12 +347,12 @@ public:
     return m_middleChildren.size() > 0;
   }
 
-  size_t num_middle_children() const
+  std::size_t num_middle_children() const
   {
     return m_middleChildren.size();
   }
 
-  Cone_tree_node* get_middle_child(size_t i) const
+  Cone_tree_node* get_middle_child(const std::size_t i) const
   {
     return m_middleChildren.at(i);
   }
@@ -426,16 +438,10 @@ public:
   {
     return m_parent != nullptr && m_parent->m_rightChild == this;
   }
-
-public:
-  Cone_expansion_event* m_pendingLeftSubtree;
-  Cone_expansion_event* m_pendingRightSubtree;
-  Cone_expansion_event* m_pendingMiddleSubtree;
-
 };
 
 } // namespace internal
-
+} // namespace Surface_mesh_shortest_paths_3
 } // namespace CGAL
 
 #endif // CGAL_SURFACE_MESH_SHORTEST_PATH_INTERNAL_CONE_TREE_H
diff --git a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/misc_functions.h b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/misc_functions.h
index 8ed7ba9a30c..7123d5fd632 100644
--- a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/misc_functions.h
+++ b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/misc_functions.h
@@ -29,41 +29,45 @@
 #include <CGAL/boost/graph/iterator.h>
 
 namespace CGAL {
-
+namespace Surface_mesh_shortest_paths_3 {
 namespace internal {
 
 template <class Triangle_3, class Triangle_mesh, class VertexPointMap>
-Triangle_3 triangle_from_halfedge(typename boost::graph_traits<Triangle_mesh>::halfedge_descriptor edge, const Triangle_mesh& g, VertexPointMap vertexPointMap)
+Triangle_3 triangle_from_halfedge(typename boost::graph_traits<Triangle_mesh>::halfedge_descriptor edge,
+                                  const Triangle_mesh& g,
+                                  const VertexPointMap vertexPointMap)
 {
   typedef typename boost::graph_traits<Triangle_mesh>::halfedge_descriptor halfedge_descriptor;
 
-  halfedge_descriptor e0 = edge;
-  halfedge_descriptor e1 = next(edge, g);
+  const halfedge_descriptor e0 = edge;
+  const halfedge_descriptor e1 = next(edge, g);
 
-  return Triangle_3(get(vertexPointMap, source(e0, g)), get(vertexPointMap, target(e0, g)), get(vertexPointMap, target(e1, g)));
+  return Triangle_3(get(vertexPointMap, source(e0, g)),
+                    get(vertexPointMap, target(e0, g)),
+                    get(vertexPointMap, target(e1, g)));
 }
 
 template <class Triangle_3, class Triangle_mesh>
-Triangle_3 triangle_from_halfedge(typename boost::graph_traits<Triangle_mesh>::halfedge_descriptor edge, const Triangle_mesh& g)
+Triangle_3 triangle_from_halfedge(typename boost::graph_traits<Triangle_mesh>::halfedge_descriptor edge,
+                                  const Triangle_mesh& g)
 {
-  return triangle_from_halfedge<Triangle_3, Triangle_mesh, typename boost::property_map<Triangle_mesh, boost::vertex_point_t>::type>(edge, g, get(boost::vertex_point, g));
+  return triangle_from_halfedge<Triangle_3>(edge, g, get(boost::vertex_point, g));
 }
 
-
 template <class Triangle_mesh>
-size_t edge_index(typename boost::graph_traits<Triangle_mesh>::halfedge_descriptor he, Triangle_mesh& p)
+std::size_t edge_index(typename boost::graph_traits<Triangle_mesh>::halfedge_descriptor he,
+                       const Triangle_mesh& p)
 {
   typedef typename boost::graph_traits<Triangle_mesh> Graph_traits;
   typedef typename Graph_traits::face_descriptor face_descriptor;
   typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor;
 
-  face_descriptor f = face(he, p);
+  const face_descriptor f = face(he, p);
 
-  halfedge_descriptor start = halfedge(f, p);
+  const halfedge_descriptor start = halfedge(f, p);
   halfedge_descriptor current = start;
 
-  size_t count = 0;
-
+  std::size_t count = 0;
   while (current != he)
   {
     current = next(current, p);
@@ -73,28 +77,8 @@ size_t edge_index(typename boost::graph_traits<Triangle_mesh>::halfedge_descript
   return count;
 }
 
-template <class FT>
-FT internal_sqrt(const FT& x, CGAL::Tag_true)
-{
-  return CGAL::sqrt(x);
-}
-
-template <class FT>
-FT internal_sqrt(const FT& x, CGAL::Tag_false)
-{
-  return FT(std::sqrt(CGAL::to_double(x)));
-}
-
-template <class FT>
-FT select_sqrt(const FT& x)
-{
-  typedef ::CGAL::Algebraic_structure_traits<FT> AST;
-  static const bool has_sqrt = ! ::boost::is_same< ::CGAL::Null_functor, typename AST::Sqrt >::value;
-  return internal_sqrt(x, ::CGAL::Boolean_tag<has_sqrt>());
-}
-
 } // namespace internal
-
+} // namespace Surface_mesh_shortest_paths_3
 } // namespace CGAL
 
 #endif // CGAL_SURFACE_MESH_SHORTEST_PATH_INTERNAL_MISC_H
diff --git a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/CMakeLists.txt b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/CMakeLists.txt
index f018de6f17b..46029b6f917 100644
--- a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/CMakeLists.txt
+++ b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/CMakeLists.txt
@@ -14,6 +14,8 @@ create_single_source_cgal_program( "Surface_mesh_shortest_path_test_1.cpp" )
 create_single_source_cgal_program( "Surface_mesh_shortest_path_test_2.cpp" )
 create_single_source_cgal_program( "Surface_mesh_shortest_path_test_3.cpp" )
 create_single_source_cgal_program( "Surface_mesh_shortest_path_test_4.cpp" )
+  create_single_source_cgal_program( "Surface_mesh_shortest_path_test_5.cpp" )
+  create_single_source_cgal_program( "Surface_mesh_shortest_path_test_6.cpp" )
 create_single_source_cgal_program( "Surface_mesh_shortest_path_traits_test.cpp" )
 
 # Link with Boost.ProgramOptions (optional)
diff --git a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_1.cpp b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_1.cpp
index 457c5f951a0..d952707649d 100644
--- a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_1.cpp
+++ b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_1.cpp
@@ -213,7 +213,7 @@ void test_simple_saddle_vertex_mesh()
   // Now test an internal face location sequence
   halfedge_descriptor firstCrossing = CGAL::halfedge(vertexHandles[4], vertexHandles[7], P).first;
 
-  size_t edgeIndex = CGAL::internal::edge_index(firstCrossing, P);
+  size_t edgeIndex = CGAL::Surface_mesh_shortest_paths_3::internal::edge_index(firstCrossing, P);
 
   Barycentric_coordinates location = construct_barycentric_coordinates(0.25, 0.5, 0.25);
 
@@ -289,6 +289,7 @@ void test_simple_saddle_vertex_mesh()
 
     CHECK_CLOSE(result.first, expectedDistances2[i], Kernel::FT(0.0001));
     assert(result.second == expectedSources2[i]);
+    CGAL_USE(expectedSources2);
     ++currentVertex;
   }
 
diff --git a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_3.cpp b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_3.cpp
index 96af1eed022..31f5198322a 100644
--- a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_3.cpp
+++ b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_3.cpp
@@ -88,7 +88,7 @@ int main(int argc, char* argv[])
     size_t faceIndex = random.get_int(0, static_cast<int>(facesList.size()));
     face_descriptor face = facesList[faceIndex];
 
-    Triangle_3 faceTriangle = CGAL::internal::triangle_from_halfedge<Triangle_3, Polyhedron_3, VPM>(halfedge(face, polyhedron), polyhedron, vertexPointMap);
+    Triangle_3 faceTriangle = CGAL::Surface_mesh_shortest_paths_3::internal::triangle_from_halfedge<Triangle_3, Polyhedron_3, VPM>(halfedge(face, polyhedron), polyhedron, vertexPointMap);
 
     Barycentric_coordinates location = CGAL::test::random_coordinates<Traits>(random);
 
diff --git a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_4.cpp b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_4.cpp
index c6ab30e0fef..3c57487983c 100644
--- a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_4.cpp
+++ b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_4.cpp
@@ -94,7 +94,7 @@ int main(int argc, char* argv[])
     size_t faceIndex = random.get_int(0, static_cast<int>(facesList.size()));
     face_descriptor face = facesList[faceIndex];
 
-    Triangle_3 faceTriangle = CGAL::internal::triangle_from_halfedge<Triangle_3, Polyhedron_3, VPM>(halfedge(face, polyhedron), polyhedron, vertexPointMap);
+    Triangle_3 faceTriangle = CGAL::Surface_mesh_shortest_paths_3::internal::triangle_from_halfedge<Triangle_3, Polyhedron_3, VPM>(halfedge(face, polyhedron), polyhedron, vertexPointMap);
 
     Barycentric_coordinates location = CGAL::test::random_coordinates<Traits>(random);
 
@@ -146,6 +146,7 @@ int main(int argc, char* argv[])
   Ray_3 outsideRay(maxPoint + (awayDir * FT(0.2)), maxPoint + (awayDir * FT(1.0)));
 
   Surface_mesh_shortest_path::Face_location emptyFaceLocation = shortestPaths.locate<AABB_face_graph_traits>(outsideRay);
+  CGAL_USE(emptyFaceLocation);
 
   assert(Graph_traits::null_face() == emptyFaceLocation.first);
 
diff --git a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_5.cpp b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_5.cpp
index aea9ead5bd8..fac75436061 100644
--- a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_5.cpp
+++ b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_5.cpp
@@ -29,11 +29,8 @@ int main(int argc, char* argv[])
   typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
   typedef CGAL::Polyhedron_3<Kernel, CGAL::Polyhedron_items_with_id_3> Polyhedron_3;
   typedef CGAL::Surface_mesh_shortest_path_traits<Kernel, Polyhedron_3> Traits;
-  typedef Traits::Barycentric_coordinates Barycentric_coordinates;
   typedef Traits::FT FT;
   typedef boost::graph_traits<Polyhedron_3> Graph_traits;
-  typedef Graph_traits::vertex_descriptor vertex_descriptor;
-  typedef Graph_traits::vertex_iterator vertex_iterator;
   typedef Graph_traits::face_descriptor face_descriptor;
   typedef Graph_traits::face_iterator face_iterator;
   typedef CGAL::Surface_mesh_shortest_path<Traits> Surface_mesh_shortest_path;
@@ -47,7 +44,6 @@ int main(int argc, char* argv[])
   std::string mesh(argv[1]);
 
   int randSeed = 4983304;
-  const size_t numTests = 15;
 
   if (argc > 2)
   {
@@ -69,6 +65,10 @@ int main(int argc, char* argv[])
   HIM halfedgeIndexMap(get(boost::halfedge_index, polyhedron));
   FIM faceIndexMap(get(boost::face_index, polyhedron));
 
+  CGAL_USE(vertexIndexMap);
+  CGAL_USE(halfedgeIndexMap);
+  CGAL_USE(faceIndexMap);
+
   face_iterator facesStart;
   face_iterator facesEnd;
 
@@ -92,12 +92,12 @@ int main(int argc, char* argv[])
 
   for (size_t i = 0; i < numInitialLocations; ++i)
   {
-    size_t faceId = rand.get_int(0, faces.size());
+    size_t faceId = static_cast<std::size_t>(rand.get_int(0, static_cast<int>(faces.size())));
     sourcePoints.push_back(Face_location(faces[faceId], CGAL::test::random_coordinates<Traits>(rand)));
     shortestPaths.add_source_point(sourcePoints.back().first, sourcePoints.back().second);
   }
 
-  BOOST_CHECK_EQUAL(numInitialLocations, shortestPaths.number_of_source_points());
+  CHECK_EQUAL(numInitialLocations, shortestPaths.number_of_source_points());
 
   size_t checkNumLocations = 0;
 
@@ -107,18 +107,16 @@ int main(int argc, char* argv[])
     ++checkNumLocations;
   }
 
-  BOOST_CHECK_EQUAL(checkNumLocations, shortestPaths.number_of_source_points());
+  CHECK_EQUAL(checkNumLocations, shortestPaths.number_of_source_points());
 
   for (Surface_mesh_shortest_path::Source_point_iterator it = shortestPaths.source_points_begin(); it != shortestPaths.source_points_end(); ++it)
   {
     Surface_mesh_shortest_path::Shortest_path_result result = shortestPaths.shortest_distance_to_source_points(it->first, it->second);
 
-    BOOST_CHECK_CLOSE(FT(0.0), result.first, FT(0.000001));
+    CHECK_CLOSE(FT(0.0), result.first, FT(0.000001));
     assert(result.second == it);
   }
 
-  size_t currentCounter = 0;
-
   // Then, remove half of them
 
   for (size_t i = 0; i < handles.size(); ++i)
@@ -129,7 +127,7 @@ int main(int argc, char* argv[])
     }
   }
 
-  BOOST_CHECK_EQUAL(numInitialLocations / 2, shortestPaths.number_of_source_points());
+  CHECK_EQUAL(numInitialLocations / 2, shortestPaths.number_of_source_points());
 
   // and ensure that they are indeed removed
   for (size_t i = 0; i < sourcePoints.size(); ++i)
@@ -138,12 +136,16 @@ int main(int argc, char* argv[])
 
     if (i % 2 != 0)
     {
-      BOOST_CHECK_CLOSE(FT(0.0), result.first, FT(0.000001));
+      CHECK_CLOSE(FT(0.0), result.first, FT(0.000001));
       assert(handles[i] == result.second);
     }
     else
     {
-      BOOST_CHECK_MESSAGE(result.first < FT(0.0) || result.first > FT(0.00001), "Incorrect resulting distance: " << result.first);
+      if ( !(result.first < FT(0.0) || result.first > FT(0.00001) ) )
+      {
+        std::cerr << "Incorrect resulting distance: " << result.first << "\n";
+        return 1;
+      }
     }
   }
 
@@ -172,12 +174,16 @@ int main(int argc, char* argv[])
 
     if (i % 3 != 0)
     {
-      BOOST_CHECK_CLOSE(FT(0.0), result.first, FT(0.000001));
+      CHECK_CLOSE(FT(0.0), result.first, FT(0.000001));
       assert(handles[i] == result.second);
     }
     else
     {
-      BOOST_CHECK_MESSAGE(result.first < FT(0.0) || result.first > FT(0.00001), "Resulted distance: " << result.first);
+      if( !(result.first < FT(0.0) || result.first > FT(0.00001)) )
+      {
+        std::cerr << "Resulted distance: " << result.first << "\n";
+        return 1;
+      }
     }
   }
    return 0;
diff --git a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_6.cpp b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_6.cpp
new file mode 100644
index 00000000000..71d52570350
--- /dev/null
+++ b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_6.cpp
@@ -0,0 +1,36 @@
+#include <iostream>
+#include <fstream>
+#include <stdlib.h>
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Surface_mesh_shortest_path.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef CGAL::Surface_mesh<Kernel::Point_3> Triangle_mesh;
+typedef CGAL::Surface_mesh_shortest_path_traits<Kernel, Triangle_mesh> Traits;
+typedef CGAL::Surface_mesh_shortest_path<Traits> Surface_mesh_shortest_path;
+typedef boost::graph_traits<Triangle_mesh> Graph_traits;
+typedef Graph_traits::vertex_iterator vertex_iterator;
+typedef Graph_traits::face_iterator face_iterator;
+
+
+int main()
+{
+  CGAL::Surface_mesh<Kernel::Point_3> mesh;
+  std::ifstream input("data/test_mesh_6.off");
+  input >> mesh;
+  input.close();
+
+  for (Triangle_mesh::Vertex_index v1 : vertices(mesh))
+    for (Triangle_mesh::Vertex_index v2 : vertices(mesh))
+    {
+      Surface_mesh_shortest_path shortest_paths(mesh);
+      shortest_paths.add_source_point(v1);
+      double dist = shortest_paths.shortest_distance_to_source_points(v2).first;
+      assert (dist==0 || v1!=v2);
+      CGAL_USE(dist);
+    }
+
+  return 0;
+}
diff --git a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_traits_test.cpp b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_traits_test.cpp
index 04744d1b3ac..5794e3c3e47 100644
--- a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_traits_test.cpp
+++ b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_traits_test.cpp
@@ -269,6 +269,7 @@ void detect_is_saddle_vertex()
 
   Traits traits;
   Traits::Is_saddle_vertex is_saddle_vertex(traits.is_saddle_vertex_object());
+  CGAL_USE(is_saddle_vertex);
 
   std::ifstream inFile("data/saddle_vertex_mesh.off");
 
diff --git a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/TestMesh.cpp b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/TestMesh.cpp
index 89f57b80826..38a1f3b33be 100644
--- a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/TestMesh.cpp
+++ b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/TestMesh.cpp
@@ -70,7 +70,7 @@ struct TestMeshProgramInstance
   TestMeshProgramInstance()
   {
     debugMode = false;
-    randomizer = NULL;
+    randomizer = nullptr;
     numVertices = 0;
     numIterations = 1;
   }
diff --git a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/check.h b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/check.h
index f762f006782..6014156d582 100644
--- a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/check.h
+++ b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/check.h
@@ -2,6 +2,9 @@
 template <typename FT, typename FT2>
 void CHECK_EQUAL(const FT& a, const FT2& b)
 {
+  if (a != b)
+    std::cerr << "ERROR: a (" << a << ") is not equal to b (" << b << ").\n";
+
   assert(a == b);
 }
 
diff --git a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/data/test_mesh_6.off b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/data/test_mesh_6.off
new file mode 100644
index 00000000000..6e64ad682f4
--- /dev/null
+++ b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/data/test_mesh_6.off
@@ -0,0 +1,60 @@
+OFF
+
+25 32 64
+0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00
+2.500000000000000e-01 0.000000000000000e+00 1.000000000000000e+00
+5.000000000000000e-01 0.000000000000000e+00 1.000000000000000e+00
+7.500000000000000e-01 0.000000000000000e+00 1.000000000000000e+00
+1.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00
+0.000000000000000e+00 2.500000000000000e-01 1.000000000000000e+00
+2.500000000000000e-01 2.500000000000000e-01 1.000000000000000e+00
+5.000000000000000e-01 2.500000000000000e-01 1.000000000000000e+00
+7.500000000000000e-01 2.500000000000000e-01 1.000000000000000e+00
+1.000000000000000e+00 2.500000000000000e-01 1.000000000000000e+00
+0.000000000000000e+00 5.000000000000000e-01 1.000000000000000e+00
+2.500000000000000e-01 5.000000000000000e-01 1.000000000000000e+00
+5.000000000000000e-01 5.000000000000000e-01 1.000000000000000e+00
+7.500000000000000e-01 5.000000000000000e-01 1.000000000000000e+00
+1.000000000000000e+00 5.000000000000000e-01 1.000000000000000e+00
+0.000000000000000e+00 7.500000000000000e-01 1.000000000000000e+00
+2.500000000000000e-01 7.500000000000000e-01 1.000000000000000e+00
+5.000000000000000e-01 7.500000000000000e-01 1.000000000000000e+00
+7.500000000000000e-01 7.500000000000000e-01 1.000000000000000e+00
+1.000000000000000e+00 7.500000000000000e-01 1.000000000000000e+00
+0.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00
+2.500000000000000e-01 1.000000000000000e+00 1.000000000000000e+00
+5.000000000000000e-01 1.000000000000000e+00 1.000000000000000e+00
+7.500000000000000e-01 1.000000000000000e+00 1.000000000000000e+00
+1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00
+3 0 1 6
+3 0 6 5
+3 5 6 11
+3 5 11 10
+3 10 11 16
+3 10 16 15
+3 15 16 21
+3 15 21 20
+3 1 2 7
+3 1 7 6
+3 6 7 12
+3 6 12 11
+3 11 12 17
+3 11 17 16
+3 16 17 22
+3 16 22 21
+3 2 3 8
+3 2 8 7
+3 7 8 13
+3 7 13 12
+3 12 13 18
+3 12 18 17
+3 17 18 23
+3 17 23 22
+3 3 4 9
+3 3 9 8
+3 8 9 14
+3 8 14 13
+3 13 14 19
+3 13 19 18
+3 18 19 24
+3 18 24 23
diff --git a/TDS_2/package_info/TDS_2/dependencies b/TDS_2/package_info/TDS_2/dependencies
index 2e02472afbb..acfc6c2a0f2 100644
--- a/TDS_2/package_info/TDS_2/dependencies
+++ b/TDS_2/package_info/TDS_2/dependencies
@@ -2,7 +2,10 @@ Algebraic_foundations
 Circulator
 Hash_map
 Installation
+Interval_support
 Kernel_23
+Modular_arithmetic
+Number_types
 Profiling_tools
 STL_Extension
 Stream_support
diff --git a/Three/package_info/Three/dependencies b/Three/package_info/Three/dependencies
index f49a6e30f24..76331fefe0d 100644
--- a/Three/package_info/Three/dependencies
+++ b/Three/package_info/Three/dependencies
@@ -1,9 +1,2 @@
-Algebraic_foundations
-GraphicsView
 Installation
-Kernel_23
-Number_types
-Profiling_tools
 STL_Extension
-Stream_support
-Three
diff --git a/ccpp.yml b/ccpp.yml
new file mode 100644
index 00000000000..75bbbee4918
--- /dev/null
+++ b/ccpp.yml
@@ -0,0 +1,12 @@
+name: C/C++ CI
+
+on: [push]
+
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+    
+    steps:
+    - name: configure all
+      run: cmake -DWITH_examples=ON -DWITH_tests=ON -DWITH_demos=ON .