Merge remote-tracking branch 'cgal/master' into Demo-Fixes_and_features-maxGimeno

2020-09-24 16:27:54 +02:00 · 2020-09-24 16:27:54 +02:00 · eccde8231b
parent c948b779f9 aaf1e8cd38
commit eccde8231b
106 changed files with 8248 additions and 1068 deletions
--- a/.github/workflows/build_doc.yml
+++ b/.github/workflows/build_doc.yml
@ -0,0 +1,41 @@
+name: Documentation
+
+on: [pull_request_target]
+
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v2.0.0
+    - name: install dependencies
+      run: |
+         set -x
+         sudo apt-get install -y graphviz ssh
+         sudo pip install lxml pyquery
+         wget -O doxygen_exe https://cgal.geometryfactory.com/~mgimeno/doxygen/build_1_8_13/bin/doxygen
+         sudo mv doxygen_exe /usr/bin/doxygen
+         sudo chmod +x /usr/bin/doxygen
+
+         git config --global user.email "maxime.gimeno@geometryfactory.com"
+         git config --global user.name "Maxime Gimeno"
+
+    - name: configure all
+      run: |
+        git clone https://maxGimeno:${{ secrets.PUSH_TO_CGAL_GITHUB_IO_TOKEN }}@github.com/CGAL/cgal.github.io.git --depth=5
+        mkdir -p build_doc && cd build_doc && cmake ../Documentation/doc
+    - name: Upload Doc
+      run: |
+        set -e
+        PR_NUMBER=$(python -c "import json; import os; y = json.load(open(os.environ['GITHUB_EVENT_PATH'])); print(y[\"number\"])")
+        mkdir -p cgal.github.io/${PR_NUMBER}
+        cd build_doc && make -j2 doc && make -j2 doc_with_postprocessing
+        cp -r ./doc_output/* ../cgal.github.io/${PR_NUMBER}/
+        cd ../cgal.github.io
+        if ! egrep -q " ${PR_NUMBER}\."; then
+          echo "<li><a href=https://cgal.github.io/${PR_NUMBER}/Manual/index.html>Manual for PR ${PR_NUMBER}.</a></li>" >> ./index.html
+        fi
+        #rsync --compress -a doc_output/* ${{ secrets.ids }}/cgal.github.io/${PR_NUMBER}/
+        git add ${PR_NUMBER} && git commit -a -m "Add ${PR_NUMBER}" && git push -u origin master
+        #ssh mgimeno@cgal.geometryfactory.com "cd /home/mgimeno/public_html/cgal.github.io && git add ${PR_NUMBER} && git commit -a -m 'Add a new directory' && git push origin master"
--- a/.github/workflows/delete_doc.yml
+++ b/.github/workflows/delete_doc.yml
@ -0,0 +1,25 @@
+name: Documentation Removal
+
+on:
+  pull_request_target:
+    types: [closed, removed]
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v2.0.0
+    - name: delete directory
+      run: |
+         set -x
+         git config --global user.email "maxime.gimeno@geometryfactory.com"
+         git config --global user.name "Maxime Gimeno"
+         git clone https://maxGimeno:${{ secrets.PUSH_TO_CGAL_GITHUB_IO_TOKEN }}@github.com/CGAL/cgal.github.io.git --depth=5
+         PR_NUMBER=$(python -c "import json; import os; y = json.load(open(os.environ['GITHUB_EVENT_PATH'])); print(y[\"number\"])")
+         cd cgal.github.io/
+         egrep -v " ${PR_NUMBER}\." index.html > tmp.html
+         mv tmp.html index.html
+         #git rm -r ${PR_NUMBER} && git commit -a -m "Remove ${PR_NUMBER}" && git push -u origin master
+         git commit -a -m "Remove ${PR_NUMBER}" && git push -u origin master
+
--- a/.travis.yml
+++ b/.travis.yml
@ -1,5 +1,5 @@
 language: cpp
-dist: xenial
+dist: bionic
 sudo: required
 git:
  depth: 3
@ -58,7 +58,7 @@ install:
  - echo "$PWD"
  - if [ -n "$TRAVIS_PULL_REQUEST_BRANCH" ] && [ "$PACKAGE" != CHECK ]; then  DO_IGNORE=FALSE;  for ARG in $(echo "$PACKAGE");do if [ "$ARG" = "Maintenance" ]; then continue; fi; . $PWD/.travis/test_package.sh "$PWD" "$ARG";  echo "DO_IGNORE is $DO_IGNORE";  if [ "$DO_IGNORE" = "FALSE" ]; then  break; fi; done; if [ "$DO_IGNORE" = "TRUE" ]; then travis_terminate 0; fi;fi
  - /usr/bin/time -f 'Spend time of %C -- %E (real)' bash .travis/install.sh
-  - export CXX=clang++ CC=clang;
+  - export CXX=clang++-10 CC=clang-10;
 before_script:
  - wget -O doxygen_exe https://cgal.geometryfactory.com/~mgimeno/doxygen_exe
  - sudo mv doxygen_exe /usr/bin/doxygen
--- a/.travis/build_package.sh
+++ b/.travis/build_package.sh
@ -11,7 +11,7 @@ function mytime {
 function build_examples {
  mkdir -p build-travis
  cd build-travis
-  mytime cmake -DCGAL_DIR="/usr/local/lib/cmake/CGAL" -DCMAKE_CXX_FLAGS="${CXX_FLAGS}" -DCGAL_BUILD_THREE_DOC=TRUE ..
+  mytime cmake -DCGAL_DIR="/usr/local/lib/cmake/CGAL" -DCMAKE_CXX_FLAGS="${CXX_FLAGS}" ..
  mytime make -j2 VERBOSE=1
 }

--- a/.travis/install.sh
+++ b/.travis/install.sh
@ -2,14 +2,15 @@

 [ -n "$CGAL_DEBUG_TRAVIS" ] && set -x
 DONE=0
+sudo add-apt-repository ppa:mikhailnov/pulseeffects -y
 sudo apt-get update
+
 while [ $DONE = 0 ]
 do
-  DONE=1 && sudo -E apt-get -yq --no-install-suggests --no-install-recommends --force-yes install clang zsh \
-flex bison cmake graphviz libgmp-dev libmpfr-dev libmpfi-dev zlib1g-dev libeigen3-dev  libboost-dev \
-libboost-system-dev libboost-program-options-dev libboost-thread-dev libboost-iostreams-dev \
+  DONE=1 && sudo -E apt-get -yq --no-install-suggests --no-install-recommends --force-yes install clang-10 zsh \
+flex bison cmake graphviz libgmp-dev libmpfr-dev libmpfi-dev zlib1g-dev libeigen3-dev  \
 qtbase5-dev libqt5sql5-sqlite libqt5opengl5-dev qtscript5-dev libqt5svg5-dev qttools5-dev qttools5-dev-tools qml-module-qtgraphicaleffects libopencv-dev mesa-common-dev libmetis-dev libglu1-mesa-dev \
-|| DONE=0 && sudo apt-get update
+libboost1.72-dev || DONE=0 && sudo apt-get update
 done
 exit 0

--- a/.travis/template.txt
+++ b/.travis/template.txt
@ -1,5 +1,5 @@
 language: cpp
-dist: xenial
+dist: bionic
 sudo: required
 git:
  depth: 3
@ -12,7 +12,7 @@ install:
  - echo "$PWD"
  - if [ -n "$TRAVIS_PULL_REQUEST_BRANCH" ] && [ "$PACKAGE" != CHECK ]; then  DO_IGNORE=FALSE;  for ARG in $(echo "$PACKAGE");do if [ "$ARG" = "Maintenance" ]; then continue; fi; . $PWD/.travis/test_package.sh "$PWD" "$ARG";  echo "DO_IGNORE is $DO_IGNORE";  if [ "$DO_IGNORE" = "FALSE" ]; then  break; fi; done; if [ "$DO_IGNORE" = "TRUE" ]; then travis_terminate 0; fi;fi
  - /usr/bin/time -f 'Spend time of %C -- %E (real)' bash .travis/install.sh
-  - export CXX=clang++ CC=clang;
+  - export CXX=clang++-10 CC=clang-10;
 before_script:
  - wget -O doxygen_exe https://cgal.geometryfactory.com/~mgimeno/doxygen_exe
  - sudo mv doxygen_exe /usr/bin/doxygen
--- a/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/reconstruction_surface_mesh.cpp
+++ b/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/reconstruction_surface_mesh.cpp
@ -1,10 +1,10 @@
 #include <iostream>
 #include <fstream>
 #include <algorithm>
+#include <array>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Advancing_front_surface_reconstruction.h>
 #include <CGAL/Surface_mesh.h>
-#include <CGAL/array.h>
 #include <CGAL/disable_warnings.h>

 typedef std::array<std::size_t,3> Facet;
--- a/Arrangement_on_surface_2/include/CGAL/Arr_dcel_base.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_dcel_base.h
@ -34,6 +34,7 @@
 #include <CGAL/Arrangement_2/Arrangement_2_iterators.h>
 #include <CGAL/assertions.h>

+#include <boost/pool/pool_alloc.hpp>

 namespace CGAL {

@ -876,7 +877,7 @@ public:
 * The arrangement DCEL class.
 */
 template <class V, class H, class F,
-          class Allocator = CGAL_ALLOCATOR(int) >
+          class Allocator = boost::fast_pool_allocator<int> >
 class Arr_dcel_base {
 public:
  // Define the vertex, halfedge and face types.
--- a/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_cache.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_cache.h
@ -653,7 +653,7 @@ void _Bezier_cache<NtTraits>::_self_intersection_params
  //   II: Y(t) - Y(s) / (t - s) = 0
  //
  Integer                 *coeffs;
-  int                      i, k;
+  int                      i;

  // Consruct the bivariate polynomial that corresponds to Equation I.
  // Note that we represent a bivariate polynomial as a vector of univariate
@ -667,12 +667,11 @@ void _Bezier_cache<NtTraits>::_self_intersection_params
  coeffs = new Integer [degX];

  for (i = 0; i < degX; i++)
-  {
-    for (k = i + 1; k < degX; k++)
-      coeffs[k - i - 1] = nt_traits.get_coefficient (polyX, k);
+    coeffs[i] = nt_traits.get_coefficient(polyX, i + 1);

-    coeffsX_st[i] = nt_traits.construct_polynomial (coeffs, degX - i - 1);
-  }
+  for (i = 0; i < degX; i++)
+    coeffsX_st[degX - i - 1] =
+        nt_traits.construct_polynomial(coeffs + i, degX - i - 1);

  delete[] coeffs;

@ -685,12 +684,11 @@ void _Bezier_cache<NtTraits>::_self_intersection_params
  coeffs = new Integer [degY];

  for (i = 0; i < degY; i++)
-  {
-    for (k = i + 1; k < degY; k++)
-      coeffs[k - i - 1] = nt_traits.get_coefficient (polyY, k);
+    coeffs[i] = nt_traits.get_coefficient(polyY, i + 1);

-    coeffsY_st[i] = nt_traits.construct_polynomial (coeffs, degY - i - 1);
-  }
+  for (i = 0; i < degY; i++)
+    coeffsY_st[degY - i - 1] =
+        nt_traits.construct_polynomial(coeffs + i, degY - i - 1);

  delete[] coeffs;

--- a/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_point_2.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_point_2.h
@ -1641,17 +1641,43 @@ void _Bezier_point_2_rep<RatKer, AlgKer, NtTrt, BndTrt>::_make_exact
  const Algebraic      t_min = nt_traits.convert (orig2.point_bound().t_min);
  const Algebraic      t_max = nt_traits.convert (orig2.point_bound().t_max);

+  bool self_intersecting = (org_it1->curve().id() == org_it2->curve().id());
+
  for (intr_it = intr_list.begin(); intr_it != intr_list.end(); ++intr_it)
  {
-    if (CGAL::compare (intr_it->s, s_min) != SMALLER &&
-        CGAL::compare (intr_it->s, s_max) != LARGER &&
-        CGAL::compare (intr_it->t, t_min) != SMALLER &&
-        CGAL::compare (intr_it->t, t_max) != LARGER)
+    auto in_bounding_interval =
+      [](const auto& s_, const auto& s_min_, const auto& s_max_) -> bool {
+      return CGAL::compare(s_, s_min_) != SMALLER &&
+             CGAL::compare(s_, s_max_) != LARGER;
+    };
+
+    bool st_in_st_range = in_bounding_interval(intr_it->s, s_min, s_max) &&
+                          in_bounding_interval(intr_it->t, t_min, t_max);
+    bool ts_in_st_range = false;
+
+    if (st_in_st_range)
    {
      // Update the originators.
-      orig1.set_parameter (intr_it->s);
-      orig2.set_parameter (intr_it->t);
+      orig1.set_parameter(intr_it->s);
+      orig2.set_parameter(intr_it->t);
+    }
+    else if (self_intersecting)
+    {
+      // check whether s is in t range, and t is in s range
+      // s and t can be interchanged in case of self intersections
+      ts_in_st_range = in_bounding_interval(intr_it->t, s_min, s_max) &&
+                       in_bounding_interval(intr_it->s, t_min, t_max);

+      if (ts_in_st_range)
+      {
+        // Update the originators.
+        orig1.set_parameter(intr_it->t);
+        orig2.set_parameter(intr_it->s);
+      }
+    }
+
+    if (st_in_st_range || ts_in_st_range)
+    {
      // Set the exact point coordinates.
      p_alg_x = new Algebraic (intr_it->x);
      p_alg_y = new Algebraic (intr_it->y);
--- a/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_x_monotone_2.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_x_monotone_2.h
@ -1139,10 +1139,11 @@ _Bezier_x_monotone_2<RatKer, AlgKer, NtTrt, BndTrt>::compare_to_left
    Originator_iterator  org = p.get_originator(_curve, _xid);

    CGAL_assertion(org != p.originators_end());
-    CGAL_assertion(_inc_to_right != cv._inc_to_right);

    if (org->point_bound().type == Bez_point_bound::VERTICAL_TANGENCY_PT)
    {
+      CGAL_assertion(_inc_to_right != cv._inc_to_right);
+
      if (! p.is_exact())
      {
        // Comparison based on the control polygon of the bounded vertical
--- a/Arrangement_on_surface_2/include/CGAL/Arr_segment_traits_2.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_segment_traits_2.h
--- a/Arrangement_on_surface_2/include/CGAL/Arrangement_on_surface_2.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arrangement_on_surface_2.h
@ -46,6 +46,8 @@
 #include <CGAL/Iterator_project.h>
 #include <CGAL/Iterator_transform.h>

+#include <boost/pool/pool_alloc.hpp>
+
 namespace CGAL {

 /*! \class Arrangement_on_surface_2
@ -64,7 +66,7 @@ class Arrangement_on_surface_2 {
 public:
  typedef GeomTraits_                                     Geometry_traits_2;
  typedef TopTraits_                                      Topology_traits;
-  typedef CGAL_ALLOCATOR(int)                             Allocator;
+  typedef boost::fast_pool_allocator<int>                 Allocator;

  // first define adaptor ...
  typedef Arr_traits_basic_adaptor_2<Geometry_traits_2>   Traits_adaptor_2;
--- a/BGL/doc/BGL/BGL.txt
+++ b/BGL/doc/BGL/BGL.txt
@ -634,6 +634,12 @@ Seam meshes are for example used in Chapter \ref PkgSurfaceMeshParameterization
 to parameterize a topological sphere by first virtually cutting it into a topological
 disk.

+\subsection BGLInheriting Inheriting from a Model of a Face Graph Concept
+In order for a type `FG` to be a model of one of the face graph concepts, a specialization of
+`boost::graph_traits<FG>` must exist. When you derive a class from `FG` this graph traits class specialization does not match.
+For such cases, a header file consuming some user defined macros is provided to make the derived class a valid model
+of the same concept. See `CGAL/boost/graph/graph_traits_inheritance_macros.h` for more details.
+
 \section BGLPartitioning Graph Partitioning

 For algorithms that operate locally, partitioning is often an easy way to parallelize
--- a/BGL/doc/BGL/Doxyfile.in
+++ b/BGL/doc/BGL/Doxyfile.in
@ -16,13 +16,16 @@ INPUT += ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/Euler_operations.h \
         ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/partition.h \
         ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/METIS/partition_graph.h \
         ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/METIS/partition_dual_graph.h \
-         ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/alpha_expansion_graphcut.h
+         ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/alpha_expansion_graphcut.h \
+         ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/graph_traits_inheritance_macros.h


 EXAMPLE_PATH               =  ${CGAL_Surface_mesh_skeletonization_EXAMPLE_DIR} \
                              ${CGAL_Surface_mesh_segmentation_EXAMPLE_DIR} \
                              ${CGAL_Polygon_mesh_processing_EXAMPLE_DIR} \
+                              ${CGAL_Surface_mesh_EXAMPLE_DIR} \
                              ${CGAL_Property_map_EXAMPLE_DIR} \
+                              ${CGAL_Polyhedron_EXAMPLE_DIR} \
                              ${CGAL_BGL_EXAMPLE_DIR}

 ALIASES +=  "bgllink{1}=<a href=\"http://www.boost.org/libs/graph/doc/\1.html\"><code>\1</code></a>"
--- a/BGL/doc/BGL/PackageDescription.txt
+++ b/BGL/doc/BGL/PackageDescription.txt
@ -713,6 +713,7 @@ user might encounter.
 - `CGAL::Graph_with_descriptor_with_graph`
 - `CGAL::Graph_with_descriptor_with_graph_property_map`
 - `CGAL::Seam_mesh`
+- `CGAL/boost/graph/graph_traits_inheritance_macros.h`

 \cgalCRPSection{Partitioning Methods}
 - `CGAL::METIS::partition_graph()`
--- a/BGL/doc/BGL/examples.txt
+++ b/BGL/doc/BGL/examples.txt
@ -31,4 +31,6 @@
 \example Surface_mesh_segmentation/extract_segmentation_into_mesh_example.cpp
 \example Polygon_mesh_processing/face_filtered_graph_example.cpp
 \example Property_map/dynamic_properties.cpp
+\example Surface_mesh/sm_derivation.cpp
+\example Polyhedron/poly_derivation.cpp
 */
--- a/BGL/include/CGAL/boost/graph/graph_traits_inheritance_macros.h
+++ b/BGL/include/CGAL/boost/graph/graph_traits_inheritance_macros.h
@ -0,0 +1,171 @@
+// Copyright (c) 2020  GeometryFactory (France).  All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org)
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Sébastien Loriot
+
+// This file is intentionally not protected against re-inclusion.
+// It's aimed at being included from within a user code to
+// make any structure inheriting from a face graph model a face graph
+// model itself
+
+// It is the responsibility of the including file to correctly set the
+// macros CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME, CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME
+// and optionally CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS.
+// They are #undefed at the end of this file.
+
+/**
+* \ingroup PkgBGLRef
+* \file CGAL/boost/graph/graph_traits_inheritance_macros.h
+* Convenience header file defining the necessary specializations and overloads to make a
+* class, inheriting from a model of a face graph concept, a model of that face graph concept itself.
+* Prior to the inclusion of this header, specific macros must be defined and those macros will be
+* undefined automatically when processing to the inclusion of this header.
+* It is possible to include the header several times if the operation must be done for several classes.
+* The macros that must be defined are the following:
+*   - `CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME`: the inheriting class. If it is a template class, it must be instantiated parameters named as in `CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS` or parameters available in the scope including the header;
+*   - `CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME`: the base class. it must be instantiated parameters named as in `CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS` or parameters available in the scope including the header;
+*   - `CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS`: (optional) if the inheriting class, a list of template parameters separated by commas (`,`) including `class/typename/integral type`.
+*
+* Some examples are provided in \ref Surface_mesh/sm_derivation.cpp and \ref Polyhedron/poly_derivation.cpp.
+*
+*/
+
+#include <CGAL/config.h>
+
+#if !defined(CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME) || !defined(CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME)
+
+CGAL_pragma_warning("\nBoth macros CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME and CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME "
+                    "must be defined if you want to use this file\n")
+
+#else
+
+#ifdef CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS
+
+namespace boost {
+
+template <CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS>
+struct graph_traits<CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME> :
+  public graph_traits<CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME>
+{};
+
+template <CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS, typename CGAL_XX_YATP>
+struct property_map<CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME, CGAL_XX_YATP> :
+  public property_map<CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME, CGAL_XX_YATP>
+{};
+
+} // boost namespace
+
+#define CGAL_PM_DT_SPEC(DTAG) \
+namespace boost {\
+template <CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS, typename CGAL_XX_YATP> \
+struct property_map<CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME, DTAG<CGAL_XX_YATP> > \
+  : property_map<CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME, DTAG<CGAL_XX_YATP> > \
+{};\
+} /* boost namespace */\
+\
+namespace CGAL { \
+template <CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS, typename CGAL_XX_YATP>\
+typename boost::property_map<CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME, DTAG<CGAL_XX_YATP> >::type \
+get(DTAG<CGAL_XX_YATP> t, CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME& g) \
+{ \
+  return get(t, static_cast<CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME&>(g)); \
+} \
+\
+template <CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS, typename CGAL_XX_YATP>\
+typename boost::property_map<CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME, DTAG<CGAL_XX_YATP> >::const_type \
+get(DTAG<CGAL_XX_YATP> t, const CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME& g) \
+{ \
+  return get(t, static_cast<const CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME&>(g)); \
+}\
+} //CGAL namespace
+
+CGAL_PM_DT_SPEC(CGAL::dynamic_vertex_property_t)
+CGAL_PM_DT_SPEC(CGAL::dynamic_halfedge_property_t)
+CGAL_PM_DT_SPEC(CGAL::dynamic_face_property_t)
+CGAL_PM_DT_SPEC(CGAL::dynamic_edge_property_t)
+
+#undef CGAL_PM_DT_SPEC
+
+namespace CGAL {
+
+template <CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS, typename CGAL_XX_YATP>
+struct graph_has_property<CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME, CGAL_XX_YATP> :
+  public CGAL::graph_has_property<CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME, CGAL_XX_YATP>
+{};
+
+
+} // CGAL namespace
+
+#undef CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS
+
+#else
+
+namespace boost {
+
+template <>
+struct graph_traits<CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME> :
+  public graph_traits<CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME>
+{};
+
+template <typename CGAL_XX_YATP>
+struct property_map<CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME, CGAL_XX_YATP> :
+  public property_map<CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME, CGAL_XX_YATP>
+{};
+
+} // boost namespace
+
+#define CGAL_PM_DT_SPEC(DTAG) \
+namespace boost {\
+template <typename CGAL_XX_YATP> \
+struct property_map<CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME, DTAG<CGAL_XX_YATP> > \
+  : property_map<CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME, DTAG<CGAL_XX_YATP> > \
+{};\
+} /* boost namespace */\
+\
+namespace CGAL { \
+template <typename CGAL_XX_YATP>\
+typename boost::property_map<CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME, DTAG<CGAL_XX_YATP> >::type \
+get(DTAG<CGAL_XX_YATP> t, CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME& g) \
+{ \
+  return get(t, static_cast<CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME&>(g)); \
+} \
+\
+template <typename CGAL_XX_YATP>\
+typename boost::property_map<CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME, DTAG<CGAL_XX_YATP> >::const_type \
+get(DTAG<CGAL_XX_YATP> t, const CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME& g) \
+{ \
+  return get(t, static_cast<const CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME&>(g)); \
+}\
+} //CGAL namespace
+
+CGAL_PM_DT_SPEC(CGAL::dynamic_vertex_property_t)
+CGAL_PM_DT_SPEC(CGAL::dynamic_halfedge_property_t)
+CGAL_PM_DT_SPEC(CGAL::dynamic_face_property_t)
+CGAL_PM_DT_SPEC(CGAL::dynamic_edge_property_t)
+
+#undef CGAL_PM_DT_SPEC
+
+namespace CGAL {
+
+template <typename CGAL_XX_YATP>
+struct graph_has_property<CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME, CGAL_XX_YATP> :
+  public CGAL::graph_has_property<CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME, CGAL_XX_YATP>
+{};
+
+
+} // CGAL namespace
+
+#endif
+
+
+#undef CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME
+#undef CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME
+
+#endif
+
--- a/BGL/include/CGAL/boost/graph/parameters_interface.h
+++ b/BGL/include/CGAL/boost/graph/parameters_interface.h
@ -101,6 +101,7 @@ CGAL_add_named_parameter(area_threshold_t, area_threshold, area_threshold)
 CGAL_add_named_parameter(halfedges_keeper_t, halfedges_keeper, halfedges_keeper)
 CGAL_add_named_parameter(volume_threshold_t, volume_threshold, volume_threshold)
 CGAL_add_named_parameter(dry_run_t, dry_run, dry_run)
+CGAL_add_named_parameter(do_not_modify_t, do_not_modify, do_not_modify)

 // List of named parameters that we use in the package 'Surface Mesh Simplification'
 CGAL_add_named_parameter(get_cost_policy_t, get_cost_policy, get_cost)
@ -172,6 +173,7 @@ CGAL_add_named_parameter(remesh_boundaries_t, remesh_boundaries, remesh_boundari
 CGAL_add_named_parameter(cell_selector_t, cell_selector, cell_selector)
 CGAL_add_named_parameter(facet_is_constrained_t, facet_is_constrained, facet_is_constrained_map)
 CGAL_add_named_parameter(remeshing_visitor_t, remeshing_visitor, remeshing_visitor)
+CGAL_add_named_parameter(smooth_constrained_edges_t, smooth_constrained_edges, smooth_constrained_edges)

 // output parameters
 CGAL_add_named_parameter(face_proxy_map_t, face_proxy_map, face_proxy_map)
--- a/BGL/test/BGL/CMakeLists.txt
+++ b/BGL/test/BGL/CMakeLists.txt
@ -71,6 +71,8 @@ create_single_source_cgal_program( "graph_concept_Linear_cell_complex.cpp" )

 create_single_source_cgal_program( "graph_concept_Arrangement_2.cpp" )

+create_single_source_cgal_program( "graph_concept_Derived.cpp" )
+
 create_single_source_cgal_program( "test_clear.cpp" )

 create_single_source_cgal_program( "test_helpers.cpp" )
@ -103,6 +105,8 @@ create_single_source_cgal_program( "test_wrl.cpp" )

 create_single_source_cgal_program( "bench_read_from_stream_vs_add_face_and_add_faces.cpp" )

+create_single_source_cgal_program( "graph_traits_inheritance.cpp" )
+
 if(OpenMesh_FOUND)
  target_link_libraries( test_clear PRIVATE ${OPENMESH_LIBRARIES})
  target_link_libraries( test_Euler_operations PRIVATE ${OPENMESH_LIBRARIES})
--- a/BGL/test/BGL/graph_concept_Derived.cpp
+++ b/BGL/test/BGL/graph_concept_Derived.cpp
@ -0,0 +1,112 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Polyhedron_items_with_id_3.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/boost/graph/graph_traits_Linear_cell_complex_for_combinatorial_map.h>
+
+#include <boost/graph/graph_concepts.hpp>
+#include <CGAL/boost/graph/graph_concepts.h>
+
+
+typedef CGAL::Simple_cartesian<double> Kernel;
+typedef Kernel::Point_3 Point_3;
+
+template <typename Traits>
+struct My_mesh_1 : public CGAL::Polyhedron_3<Traits, CGAL::Polyhedron_items_with_id_3> {};
+
+struct My_mesh_2 : public CGAL::Polyhedron_3<Kernel, CGAL::Polyhedron_items_with_id_3> {};
+
+template <typename PT>
+struct My_mesh_3 : public CGAL::Surface_mesh<PT> {};
+
+struct My_mesh_5 : public CGAL::Surface_mesh<Kernel::Point_3> {};
+
+// dim could be hard-coded but for the purpose of the example it is left
+template <int dim, typename K>
+struct My_mesh_4 :
+  CGAL::Linear_cell_complex_for_bgl_combinatorial_map_helper
+         <2, dim, CGAL::Linear_cell_complex_traits<dim, K> >::type
+{};
+
+/// make My_mesh_1 a valid face graph model
+#define CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS typename Traits
+#define CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME My_mesh_1<Traits>
+#define CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME CGAL::Polyhedron_3<Traits, CGAL::Polyhedron_items_with_id_3>
+#include <CGAL/boost/graph/graph_traits_inheritance_macros.h>
+
+/// make My_mesh_2 a valid face graph model
+// no template parameter, CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS is then not defined
+#define CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME My_mesh_2
+#define CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME CGAL::Polyhedron_3<Kernel, CGAL::Polyhedron_items_with_id_3>
+#include <CGAL/boost/graph/graph_traits_inheritance_macros.h>
+
+/// make My_mesh_3 a valid face graph model
+#define CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS typename PT
+#define CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME My_mesh_3<PT>
+#define CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME CGAL::Surface_mesh<PT>
+#include <CGAL/boost/graph/graph_traits_inheritance_macros.h>
+
+/// make My_mesh_4 a valid face graph model
+#define CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS int dim, typename K
+#define CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME My_mesh_4<dim, K>
+#define CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME typename CGAL::Linear_cell_complex_for_bgl_combinatorial_map_helper\
+         <2, dim, CGAL::Linear_cell_complex_traits<dim, K> >::type
+#include <CGAL/boost/graph/graph_traits_inheritance_macros.h>
+
+/// make My_mesh_5 a valid face graph model
+// no template parameter, CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS is then not defined
+#define CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME My_mesh_5
+#define CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME CGAL::Surface_mesh<Kernel::Point_3>
+#include <CGAL/boost/graph/graph_traits_inheritance_macros.h>
+
+
+template <class Graph>
+void concept_check()
+{
+  typedef boost::graph_traits< Graph > Traits;
+  typedef typename Traits::edge_descriptor edge_descriptor;
+  typedef typename Traits::halfedge_descriptor halfedge_descriptor;
+  typedef typename Traits::vertex_descriptor vertex_descriptor;
+  typedef typename Traits::face_descriptor face_descriptor;
+
+  boost::function_requires< boost::GraphConcept<Graph> >();
+  boost::function_requires< boost::VertexListGraphConcept<Graph> >();
+  boost::function_requires< boost::EdgeListGraphConcept<Graph> >();
+  boost::function_requires< boost::IncidenceGraphConcept<Graph> >();
+  boost::function_requires< boost::AdjacencyMatrixConcept<Graph> >();
+  boost::function_requires< boost::BidirectionalGraphConcept<Graph> >();
+  boost::function_requires< CGAL::HalfedgeGraphConcept<Graph> >();
+  boost::function_requires< CGAL::HalfedgeListGraphConcept<Graph> >();
+  boost::function_requires< CGAL::FaceGraphConcept<Graph> >();
+  boost::function_requires< CGAL::FaceListGraphConcept<Graph> >();
+  boost::function_requires< CGAL::MutableHalfedgeGraphConcept<Graph> >();
+  boost::function_requires< CGAL::MutableFaceGraphConcept<Graph> >();
+
+  boost::function_requires< boost::concepts::ReadablePropertyGraph<
+    Graph, halfedge_descriptor, CGAL::halfedge_index_t> >();
+  boost::function_requires< boost::concepts::ReadablePropertyGraph<
+    Graph, edge_descriptor, boost::edge_index_t> >();
+  boost::function_requires< boost::concepts::ReadablePropertyGraph<
+    Graph, edge_descriptor, boost::edge_weight_t> >();
+  boost::function_requires< boost::concepts::PropertyGraph<
+    Graph, vertex_descriptor, CGAL::vertex_point_t> >();
+  boost::function_requires< boost::concepts::ReadablePropertyGraph<
+    Graph, vertex_descriptor, boost::vertex_index_t> >();
+  boost::function_requires< boost::concepts::ReadablePropertyGraph<
+    Graph, face_descriptor, CGAL::face_index_t> >();
+
+  // null
+  boost::graph_traits<Graph>::null_vertex();
+  boost::graph_traits<Graph>::null_halfedge();
+  boost::graph_traits<Graph>::null_face();
+}
+
+int main()
+{
+  concept_check<My_mesh_1<Kernel>>();
+  concept_check<My_mesh_2>();
+  concept_check<My_mesh_3<Point_3> >();
+  concept_check<My_mesh_4<3,Kernel>>();
+  concept_check<My_mesh_5>();
+  return 0;
+}
--- a/BGL/test/BGL/graph_traits_inheritance.cpp
+++ b/BGL/test/BGL/graph_traits_inheritance.cpp
@ -0,0 +1,79 @@
+#include <vector>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/boost/graph/graph_traits_Linear_cell_complex_for_combinatorial_map.h>
+#include <CGAL/convex_hull_3.h>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+
+template <typename Traits>
+struct My_mesh_1 : public CGAL::Polyhedron_3<Traits> {};
+
+struct My_mesh_2 : public CGAL::Polyhedron_3<Kernel> {};
+
+template <typename PT>
+struct My_mesh_3 : public CGAL::Surface_mesh<PT> {};
+
+struct My_mesh_5 : public CGAL::Surface_mesh<Kernel::Point_3> {};
+
+// dim could be hard-coded but for the purpose of the example it is left
+template <int dim, typename K>
+struct My_mesh_4 :
+  CGAL::Linear_cell_complex_for_bgl_combinatorial_map_helper
+         <2, dim, CGAL::Linear_cell_complex_traits<dim, K> >::type
+{};
+
+/// make My_mesh_1 a valid face graph model
+#define CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS typename Traits
+#define CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME My_mesh_1<Traits>
+#define CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME CGAL::Polyhedron_3<Traits>
+#include <CGAL/boost/graph/graph_traits_inheritance_macros.h>
+
+/// make My_mesh_2 a valid face graph model
+// no template parameter, CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS is then not defined
+#define CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME My_mesh_2
+#define CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME CGAL::Polyhedron_3<Kernel>
+#include <CGAL/boost/graph/graph_traits_inheritance_macros.h>
+
+/// make My_mesh_3 a valid face graph model
+#define CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS typename PT
+#define CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME My_mesh_3<PT>
+#define CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME CGAL::Surface_mesh<PT>
+#include <CGAL/boost/graph/graph_traits_inheritance_macros.h>
+
+/// make My_mesh_4 a valid face graph model
+#define CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS int dim, typename K
+#define CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME My_mesh_4<dim, K>
+#define CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME typename CGAL::Linear_cell_complex_for_bgl_combinatorial_map_helper\
+         <2, dim, CGAL::Linear_cell_complex_traits<dim, K> >::type
+#include <CGAL/boost/graph/graph_traits_inheritance_macros.h>
+
+/// make My_mesh_5 a valid face graph model
+// no template parameter, CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS is then not defined
+#define CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME My_mesh_5
+#define CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME CGAL::Surface_mesh<Kernel::Point_3>
+#include <CGAL/boost/graph/graph_traits_inheritance_macros.h>
+
+int main()
+{
+  typedef My_mesh_1<Kernel>       Mesh1;
+
+  std::vector<Kernel::Point_3> points;
+  Mesh1 poly1;
+  CGAL::convex_hull_3(points.begin(), points.end(), poly1);
+
+  My_mesh_2 poly2;
+  CGAL::convex_hull_3(points.begin(), points.end(), poly2);
+
+  My_mesh_3<Kernel::Point_3> poly3;
+  CGAL::convex_hull_3(points.begin(), points.end(), poly3);
+
+  My_mesh_4<3, Kernel> poly4;
+  CGAL::convex_hull_3(points.begin(), points.end(), poly4);
+
+  My_mesh_5 poly5;
+  CGAL::convex_hull_3(points.begin(), points.end(), poly5);
+
+  return 0;
+}
--- a/BGL/test/BGL/test_cgal_bgl_named_params.cpp
+++ b/BGL/test/BGL/test_cgal_bgl_named_params.cpp
@ -97,6 +97,7 @@ void test(const NamedParameters& np)
  assert(get_parameter(np, CGAL::internal_np::do_lock_mesh).v == 61);
  assert(get_parameter(np, CGAL::internal_np::halfedges_keeper).v == 62);
  assert(get_parameter(np, CGAL::internal_np::do_simplify_border).v == 64);
+  assert(get_parameter(np, CGAL::internal_np::do_not_modify).v == 65);
  assert(get_parameter(np, CGAL::internal_np::maximum_number_of_faces).v == 78910);

    // Named parameters that we use in the package 'Surface Mesh Simplification'
@ -184,6 +185,7 @@ void test(const NamedParameters& np)
  check_same_type<54>(get_parameter(np, CGAL::internal_np::use_area_smoothing));
  check_same_type<55>(get_parameter(np, CGAL::internal_np::use_Delaunay_flips));
  check_same_type<56>(get_parameter(np, CGAL::internal_np::use_safety_constraints));
+  check_same_type<65>(get_parameter(np, CGAL::internal_np::do_not_modify));

  check_same_type<12340>(get_parameter(np, CGAL::internal_np::do_self_intersection_tests));
  check_same_type<12341>(get_parameter(np, CGAL::internal_np::do_orientation_tests));
@ -353,6 +355,7 @@ int main()
                         .halfedges_keeper(A<62>(62))
                         .use_convex_hull(A<63>(63))
                         .do_simplify_border(A<64>(64))
+                         .do_not_modify(A<65>(65))
                         .point_map(A<9000>(9000))
                         .query_point_map(A<9001>(9001))
                         .normal_map(A<9002>(9002))
--- a/CGAL_Core/include/CGAL/CORE/poly/Poly.tcc
+++ b/CGAL_Core/include/CGAL/CORE/poly/Poly.tcc
@ -670,19 +670,20 @@ Polynomial<NT> Polynomial<NT>::pseudoRemainder (
  contract();         // Let A = (*this).  Contract A.
  Polynomial<NT> tmpB(B);
  tmpB.contract();    // local copy of B
+  int bTrueDegree = tmpB.degree;
  C = NT(1);  // Initialized to C=1.
-  if (B.degree == -1)  {
+  if (bTrueDegree == -1)  {
    core_error("ERROR in Polynomial<NT>::pseudoRemainder :\n    -- divide by zero polynomial", __FILE__, __LINE__, false);
    return Polynomial(0);  // Unit Polynomial (arbitrary!)
  }
-  if (B.degree > degree) {
+  if (bTrueDegree > degree) {
    return Polynomial(); // Zero Polynomial
    // CHECK: 1*THIS = 0*B + THAT,  deg(THAT) < deg(B)
  }

  Polynomial<NT> Quo;  // accumulate the return polynomial, Quo
  Polynomial<NT> tmpQuo;
-  while (degree >= B.degree) {  // INVARIANT: C*A = B*Quo + (*this)
+  while (degree >= bTrueDegree) {  // INVARIANT: C*A = B*Quo + (*this)
    tmpQuo = reduceStep(tmpB);  // Let (*this) be (*oldthis), which
                                // is transformed into (*newthis). Then,
                                //     c*(*oldthis) = B*m + (*newthis)
--- a/Convex_hull_3/doc/Convex_hull_3/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h
+++ b/Convex_hull_3/doc/Convex_hull_3/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h
@ -4,12 +4,12 @@ namespace CGAL {
 \ingroup PkgConvexHull3Functions

 \brief computes robustly the intersection of the halfspaces defined by the planes contained in the range [`begin`, `end`) without constructing the dual points. The result is stored in the polyhedron `pm`.
-If `origin` is given then it must be a point strictly inside the polyhedron. If an interior point is not given then it is computed using a linear program and thus is slower.
+If `origin` is given then it must be a point strictly inside the polyhedron. If an interior point is not given then it is computed using the function `halfspace_intersection_interior_point_3()` based on solving a linear program and thus is slower.

 This version does not construct the dual points explicitely but uses a special traits class for the function `CGAL::convex_hull_3()` to handle predicates on dual points without constructing them.

-\attention Halfspaces are considered as lower halfspaces that is to say if the plane's equation is \f$ a\, x +b\, y +c\, z + d = 0 \f$ then the corresponding halfspace is defined by \f$ a\, x +b\, y +c\, z + d \le 0 \f$ .
-\attention
+Halfspaces are considered as lower halfspaces, that is if the plane equation is \f$ a\, x +b\, y +c\, z + d = 0 \f$ then the corresponding halfspace is defined by \f$ a\, x +b\, y +c\, z + d \le 0 \f$ .
+

 \pre The point type of `origin` and the point type of the vertices of `PolygonMesh` must come from the same \cgal %Kernel.\pre if provided, `origin` is inside the intersection of halfspaces defined by the range `[begin, end)`.
 \pre The computed intersection must be a bounded convex polyhedron.
--- a/Convex_hull_3/doc/Convex_hull_3/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h
+++ b/Convex_hull_3/doc/Convex_hull_3/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h
@ -4,10 +4,10 @@ namespace CGAL {
 \ingroup PkgConvexHull3Functions

 \brief computes the intersection of the halfspaces defined by the planes contained in the range [`begin`, `end`). The result is stored in the polyhedron `pm`.
-If `origin` is given then it must be a point strictly inside the polyhedron. If an interior point is not given then it is computed using a linear program and thus is slower.
+If `origin` is given then it must be a point strictly inside the polyhedron. If an interior point is not given then it is computed using the function `halfspace_intersection_interior_point_3()` based on solving a linear program and thus is slower.
 This version constructs explicitly the dual points using the convex hull algorithm parametrized with the given traits class.

-\attention Halfspaces are considered as lower halfspaces that is to say if the plane's equation is \f$ a\, x +b\, y +c\, z + d = 0 \f$ then the corresponding halfspace is defined by \f$ a\, x +b\, y +c\, z + d \le 0 \f$ .
+Halfspaces are considered as lower halfspaces, that is if the plane equation is \f$ a\, x +b\, y +c\, z + d = 0 \f$ then the corresponding halfspace is defined by \f$ a\, x +b\, y +c\, z + d \le 0 \f$ .

 \pre The value type of `PlaneIterator` and the point type of `origin` must come from the same \cgal Kernel.
 \pre if provided, `origin` is inside the intersection of halfspaces defined by the range `[begin, end)`.
--- a/Convex_hull_3/doc/Convex_hull_3/Doxyfile.in
+++ b/Convex_hull_3/doc/Convex_hull_3/Doxyfile.in
@ -1,3 +1,4 @@
@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}

 PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - 3D Convex Hulls"
+INPUT += ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Convex_hull_3/dual/halfspace_intersection_interior_point_3.h
--- a/Convex_hull_3/doc/Convex_hull_3/PackageDescription.txt
+++ b/Convex_hull_3/doc/Convex_hull_3/PackageDescription.txt
@ -84,7 +84,8 @@ defining `CGAL_CH_CHECK_EXPENSIVE`.

 \cgalCRPSection{Halfspace Intersection Functions}

- `CGAL::halfspace_intersection_3`
- `CGAL::halfspace_intersection_with_constructions_3`
+- `CGAL::halfspace_intersection_3()`
+- `CGAL::halfspace_intersection_with_constructions_3()`
+- `CGAL::halfspace_intersection_interior_point_3()`

 */
--- a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h
+++ b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h
@ -26,7 +26,7 @@
 #include <CGAL/assertions.h>
 #include <CGAL/boost/graph/Euler_operations.h>
 // For interior_polyhedron_3
-#include <CGAL/Convex_hull_3/dual/interior_polyhedron_3.h>
+#include <CGAL/Convex_hull_3/dual/halfspace_intersection_interior_point_3.h>
 #include <CGAL/internal/Exact_type_selector.h>

 #include <boost/unordered_map.hpp>
@ -241,16 +241,12 @@ namespace CGAL

        // if a point inside is not provided find one using linear programming
        if (!origin) {
-          // choose exact integral type
-          typedef typename internal::Exact_field_selector<void*>::Type ET;
-
          // find a point inside the intersection
-          typedef Interior_polyhedron_3<K, ET> Interior_polyhedron;
-          Interior_polyhedron interior;
-          CGAL_assertion_code(bool interior_point_found = )
-          interior.find(begin, end);
-          CGAL_assertion_msg(interior_point_found, "halfspace_intersection_3: problem when determing a point inside the intersection");
-          origin = boost::make_optional(interior.inside_point());
+          origin = halfspace_intersection_interior_point_3(begin, end);
+
+          CGAL_assertion_msg(origin!=boost::none, "halfspace_intersection_3: problem when determing a point inside the intersection");
+          if (origin==boost::none)
+            return;
        }

        // make sure the origin is on the negative side of all the planes
--- a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_interior_point_3.h
+++ b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_interior_point_3.h
@ -0,0 +1,198 @@
+// Copyright (c) 2014  INRIA Sophia-Antipolis (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Jocelyn Meyron
+//                 Pierre Alliez
+//
+
+#ifndef CGAL_CH3_DUAL_HALFSPACE_INTERIOR_POINT_3_H
+#define CGAL_CH3_DUAL_HALFSPACE_INTERIOR_POINT_3_H
+
+#include <CGAL/license/Convex_hull_3.h>
+
+#include <CGAL/disable_warnings.h>
+
+// LP solver to compute an interior point of a polyhedron
+#include <CGAL/QP_functions.h>
+#include <CGAL/QP_models.h>
+#include <limits>
+#include <CGAL/number_utils.h>
+#include <CGAL/assertions.h>
+#include <CGAL/internal/Exact_type_selector.h>
+
+// Description taken from http://www.qhull.org/html/qhalf.htm
+
+// If you do not know an interior point for the halfspaces, use linear programming
+// to find one. Assume, n halfspaces defined by: aj*x1+bj*x2+cj*x3+dj>=0, j=1..n.
+// Perform the following linear program:
+//                max(x5) aj*x1+bj*x2+cj*x3+dj*x4-x5>=0, j=1..n
+
+// Then, if [x1,x2,x3,x4,x5] is an optimal m_solution with x4,x5>0 we get:
+//                aj*(x1/x4)+bj*(x2/x4)+cj*(x3/x4)+dj>=(x5/x4)>0, j=1..n
+// and conclude that the point [x1/x4,x2/x4,x3/x4] is in the interior of all
+// the halfspaces. Note that x5 is optimal, so this point is "way in" the
+// interior (good for precision errors).
+
+// After finding an interior point, the rest of the intersection algorithm is
+// from Preparata & Shamos ['85, p. 316, "A simple case ..."]. Translate the
+// halfspaces so that the interior point is the origin. Calculate the dual
+// polytope. The dual polytope is the convex hull of the vertices dual to the
+// original faces in regard to the unit sphere (i.e., halfspaces at distance
+// d from the origin are dual to vertices at distance 1/d). Then calculate
+// the resulting polytope, which is the dual of the dual polytope, and
+// translate the origin back to the interior point [S. Spitz and S. Teller].
+
+// NOTE here we change this to max(x4) under constraints aj*x1 + bj*x2 + cj*x3 + dj - x4 >= 0, j=1..n
+//  i.e. aj*x1 + bj*x2 + cj*x3 - x4 >= -dj, j=1..n
+// Then, if [x1,x2,x3,x4] is an optimal m_solution with x3 > 0 we pick
+// the point [x1,x2,x3] as inside point.
+
+namespace CGAL {
+
+namespace internal {
+
+template <class Kernel, class ET>
+class Interior_polyhedron_3 {
+  // 3D
+  typedef typename Kernel::FT FT;
+  typedef typename Kernel::Plane_3 Plane;
+  typedef typename Kernel::Point_3 Point;
+
+  // program and solution types
+  typedef CGAL::Quadratic_program<double> LP;
+  typedef typename CGAL::Quadratic_program_solution<ET> Solution;
+  typedef typename Solution::Variable_value_iterator Variable_value_iterator;
+  typedef CGAL::Real_embeddable_traits<typename Variable_value_iterator::value_type> RE_traits;
+  typename RE_traits::To_double to_double;
+  Solution m_solution;
+  Point m_inside_point;
+  Point m_optimal_point;
+
+public:
+  Point& inside_point() { return m_inside_point; }
+  const Point& inside_point() const { return m_inside_point; }
+  Point& optimal_point() { return m_optimal_point; }
+  const Point& optimal_point() const { return m_optimal_point; }
+
+  // Determines if a value is infinite or not
+  template<typename T>
+  inline bool isinf(T value) {
+      return value == std::numeric_limits<T>::infinity();
+  }
+
+  // Find a point inside the polyhedron defined by a list of planes
+  // InputIterator::value_type = Plane
+  template < class InputIterator >
+  bool find(InputIterator begin, InputIterator end) {
+    // solve linear program
+    LP lp(CGAL::LARGER,false); // with constraints Ax >= b
+
+    // column indices
+    const int index_x1 = 0;
+    const int index_x2 = 1;
+    const int index_x3 = 2;
+    const int index_x4 = 3;
+
+    // assemble linear program
+    int j = 0; // row index
+    // iterate over segments
+    InputIterator it;
+    for(it = begin; it != end; ++it, j++) {
+      const Plane& plane = *it;
+      const double aj = CGAL::to_double(plane.a());
+      const double bj = CGAL::to_double(plane.b());
+      const double cj = CGAL::to_double(plane.c());
+      const double dj = CGAL::to_double(plane.d());
+
+      CGAL_assertion(!isinf(aj));
+      CGAL_assertion(!isinf(bj));
+      CGAL_assertion(!isinf(cj));
+      CGAL_assertion(!isinf(dj));
+
+      // plane defined the halfspace: aj * x1 + bj * x2 + cj * x3 + dj <= 0
+      // <=> - (aj * x1 + bj * x2 + cj * x3 + dj) >= 0
+      // j^th constraint: -(aj * x1 + bj * x2 + cj * x3 + x4) >= dj
+      lp.set_a(index_x1, j,   -aj);
+      lp.set_a(index_x2, j,   -bj);
+      lp.set_a(index_x3, j,   -cj);
+      lp.set_a(index_x4, j, -1.0);
+
+      // right hand side
+      lp.set_b(j, dj);
+    }
+
+    // objective function -> max x4 (negative sign set because
+    // the lp solver always minimizes an objective function)
+    lp.set_c(index_x4,-1.0);
+
+    // solve the linear program
+    m_solution = CGAL::solve_linear_program(lp, ET());
+
+    if(m_solution.is_infeasible())
+        return false;
+
+    if(!m_solution.is_optimal())
+        return false;
+
+    // get variables
+    Variable_value_iterator X = m_solution.variable_values_begin();
+
+    // solution if x4 > 0
+    double x4 = to_double(X[index_x4]);
+    if(x4 <= 0.0)
+        return false;
+
+    // define inside point as (x1;x2;x3)
+    double x1 = to_double(X[index_x1]);
+    double x2 = to_double(X[index_x2]);
+    double x3 = to_double(X[index_x3]);
+    m_inside_point = Point(x1,x2,x3);
+
+    return true;
+  }
+};
+
+} // end of internal namespace
+
+/*!
+\ingroup PkgConvexHull3Functions
+
+computes a point belonging to the intersection of the halfspaces defined by the planes contained in the range `[begin, end)`.
+If the intersection is empty, `boost::none` is returned.
+
+\attention Halfspaces are considered as lower halfspaces that is to say if the plane's equation
+is \f$ a\, x +b\, y +c\, z + d = 0 \f$ then the corresponding halfspace is defined by \f$ a\, x +b\, y +c\, z + d \le 0 \f$ .
+
+\tparam PlaneIterator must be an input iterator with the value type being a `Plane_3` object from \cgal Kernel
+*/
+template <class PlaneIterator>
+boost::optional<typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel::Point_3>
+halfspace_intersection_interior_point_3(PlaneIterator begin, PlaneIterator end)
+{
+  // Types
+  typedef typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel K;
+
+  // choose exact integral type
+  typedef typename internal::Exact_field_selector<void*>::Type ET;
+
+  // find a point inside the intersection
+  internal::Interior_polyhedron_3<K, ET> interior;
+  if (!interior.find(begin, end))
+    return boost::none;
+
+  return boost::make_optional(interior.inside_point());
+}
+
+} // namespace CGAL
+
+#include <CGAL/enable_warnings.h>
+
+#endif // CGAL_CH3_DUAL_HALFSPACE_INTERIOR_POINT_3_H
+
--- a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h
+++ b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h
@ -23,7 +23,7 @@
 #include <CGAL/assertions.h>

 // For interior_polyhedron_3
-#include <CGAL/Convex_hull_3/dual/interior_polyhedron_3.h>
+#include <CGAL/Convex_hull_3/dual/halfspace_intersection_interior_point_3.h>
 #include <CGAL/internal/Exact_type_selector.h>

 #include <boost/unordered_map.hpp>
@ -90,29 +90,24 @@ namespace CGAL
          void halfspace_intersection_with_constructions_3(PlaneIterator pbegin,
                                                           PlaneIterator pend,
                                                           Polyhedron &P,
-                                                           boost::optional<typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel::Point_3> const& origin,
+                                                           boost::optional<typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel::Point_3> origin,
                                                           const Traits & ch_traits) {
          typedef typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel K;
          typedef typename K::Point_3 Point;
          typedef typename K::Plane_3 Plane;

-          Point p_origin;
-
-          if (origin) {
-            p_origin = boost::get(origin);
-          } else {
-            // choose exact integral type
-            typedef typename internal::Exact_field_selector<void*>::Type ET;
-
+          // if a point inside is not provided find one using linear programming
+          if (!origin) {
            // find a point inside the intersection
-            typedef Interior_polyhedron_3<K, ET> Interior_polyhedron;
-            Interior_polyhedron interior;
-            CGAL_assertion_code(bool res = )
-              interior.find(pbegin, pend);
-            CGAL_assertion_msg(res, "halfspace_intersection_with_constructions_3: problem when determing a point inside");
-            p_origin = interior.inside_point();
+            origin = halfspace_intersection_interior_point_3(pbegin, pend);
+
+            CGAL_assertion_msg(origin!=boost::none, "halfspace_intersection_with_constructions_3: problem when determing a point inside the intersection");
+            if (origin==boost::none)
+              return;
          }

+          const Point p_origin = *origin;
+
          // construct dual points to apply the convex hull
          std::vector<Point> dual_points;
          for (PlaneIterator p = pbegin; p != pend; ++p) {
--- a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/interior_polyhedron_3.h
+++ b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/interior_polyhedron_3.h
@ -1,163 +0,0 @@
-// Copyright (c) 2014  INRIA Sophia-Antipolis (France).
-// All rights reserved.
-//
-// This file is part of CGAL (www.cgal.org).
-//
-// $URL$
-// $Id$
-// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
-//
-//
-// Author(s)     : Jocelyn Meyron
-//                 Pierre Alliez
-//
-
-#ifndef CGAL_CH3_INTERIOR_POLYHEDRON_3_H
-#define CGAL_CH3_INTERIOR_POLYHEDRON_3_H
-
-#include <CGAL/license/Convex_hull_3.h>
-
-#include <CGAL/disable_warnings.h>
-
-// LP solver to compute an interior point of a polyhedron
-#include <CGAL/QP_functions.h>
-#include <CGAL/QP_models.h>
-#include <limits>
-#include <CGAL/number_utils.h>
-#include <CGAL/assertions.h>
-
-// Description taken from http://www.qhull.org/html/qhalf.htm
-
-// If you do not know an interior point for the halfspaces, use linear programming
-// to find one. Assume, n halfspaces defined by: aj*x1+bj*x2+cj*x3+dj>=0, j=1..n.
-// Perform the following linear program:
-//                max(x5) aj*x1+bj*x2+cj*x3+dj*x4-x5>=0, j=1..n
-
-// Then, if [x1,x2,x3,x4,x5] is an optimal m_solution with x4,x5>0 we get:
-//                aj*(x1/x4)+bj*(x2/x4)+cj*(x3/x4)+dj>=(x5/x4)>0, j=1..n
-// and conclude that the point [x1/x4,x2/x4,x3/x4] is in the interior of all
-// the halfspaces. Note that x5 is optimal, so this point is "way in" the
-// interior (good for precision errors).
-
-// After finding an interior point, the rest of the intersection algorithm is
-// from Preparata & Shamos ['85, p. 316, "A simple case ..."]. Translate the
-// halfspaces so that the interior point is the origin. Calculate the dual
-// polytope. The dual polytope is the convex hull of the vertices dual to the
-// original faces in regard to the unit sphere (i.e., halfspaces at distance
-// d from the origin are dual to vertices at distance 1/d). Then calculate
-// the resulting polytope, which is the dual of the dual polytope, and
-// translate the origin back to the interior point [S. Spitz and S. Teller].
-
-// NOTE here we change this to max(x4) under constraints aj*x1 + bj*x2 + cj*x3 + dj - x4 >= 0, j=1..n
-//  i.e. aj*x1 + bj*x2 + cj*x3 - x4 >= -dj, j=1..n
-// Then, if [x1,x2,x3,x4] is an optimal m_solution with x3 > 0 we pick
-// the point [x1,x2,x3] as inside point.
-
-namespace CGAL {
-template <class Kernel, class ET>
-class Interior_polyhedron_3 {
-        // 3D
-        typedef typename Kernel::FT FT;
-        typedef typename Kernel::Plane_3 Plane;
-        typedef typename Kernel::Point_3 Point;
-
-        // program and solution types
-        typedef CGAL::Quadratic_program<double> LP;
-        typedef typename CGAL::Quadratic_program_solution<ET> Solution;
-        typedef typename Solution::Variable_value_iterator Variable_value_iterator;
-        typedef CGAL::Real_embeddable_traits<typename Variable_value_iterator::value_type> RE_traits;
-        typename RE_traits::To_double to_double;
-        Solution m_solution;
-        Point m_inside_point;
-        Point m_optimal_point;
-
-    public:
-        Point& inside_point() { return m_inside_point; }
-        const Point& inside_point() const { return m_inside_point; }
-        Point& optimal_point() { return m_optimal_point; }
-        const Point& optimal_point() const { return m_optimal_point; }
-
-        // Determines if a value is infinite or not
-        template<typename T>
-        inline bool isinf(T value) {
-            return value == std::numeric_limits<T>::infinity();
-        }
-
-        // Find a point inside the polyhedron defined by a list of planes
-        // InputIterator::value_type = Plane
-        template < class InputIterator >
-        bool find(InputIterator begin, InputIterator end) {
-            // solve linear program
-            LP lp(CGAL::LARGER,false); // with constraints Ax >= b
-
-            // column indices
-            const int index_x1 = 0;
-            const int index_x2 = 1;
-            const int index_x3 = 2;
-            const int index_x4 = 3;
-
-            // assemble linear program
-            int j = 0; // row index
-            // iterate over segments
-            InputIterator it;
-            for(it = begin; it != end; ++it, j++) {
-                const Plane& plane = *it;
-                const double aj = CGAL::to_double(plane.a());
-                const double bj = CGAL::to_double(plane.b());
-                const double cj = CGAL::to_double(plane.c());
-                const double dj = CGAL::to_double(plane.d());
-
-                CGAL_assertion(!isinf(aj));
-                CGAL_assertion(!isinf(bj));
-                CGAL_assertion(!isinf(cj));
-                CGAL_assertion(!isinf(dj));
-
-                // plane defined the halfspace: aj * x1 + bj * x2 + cj * x3 + dj <= 0
-                // <=> - (aj * x1 + bj * x2 + cj * x3 + dj) >= 0
-                // j^th constraint: -(aj * x1 + bj * x2 + cj * x3 + x4) >= dj
-                lp.set_a(index_x1, j,   -aj);
-                lp.set_a(index_x2, j,   -bj);
-                lp.set_a(index_x3, j,   -cj);
-                lp.set_a(index_x4, j, -1.0);
-
-                // right hand side
-                lp.set_b(j, dj);
-            }
-
-            // objective function -> max x4 (negative sign set because
-            // the lp solver always minimizes an objective function)
-            lp.set_c(index_x4,-1.0);
-
-            // solve the linear program
-            m_solution = CGAL::solve_linear_program(lp, ET());
-
-            if(m_solution.is_infeasible())
-                return false;
-
-            if(!m_solution.is_optimal())
-                return false;
-
-            // get variables
-            Variable_value_iterator X = m_solution.variable_values_begin();
-
-            // solution if x4 > 0
-            double x4 = to_double(X[index_x4]);
-            if(x4 <= 0.0)
-                return false;
-
-            // define inside point as (x1;x2;x3)
-            double x1 = to_double(X[index_x1]);
-            double x2 = to_double(X[index_x2]);
-            double x3 = to_double(X[index_x3]);
-            m_inside_point = Point(x1,x2,x3);
-
-            return true;
-        }
-};
-
-} // namespace CGAL
-
-#include <CGAL/enable_warnings.h>
-
-#endif // CGAL_CH3_INTERIOR_POLYHEDRON_3_H
-
--- a/Documentation/doc/Documentation/Developer_manual/Chapter_portability.txt
+++ b/Documentation/doc/Documentation/Developer_manual/Chapter_portability.txt
@ -101,70 +101,70 @@ the following table.
    <TD ALIGN=LEFT NOWRAP>
 GNU 3.2.1
    <TD ALIGN=LEFT NOWRAP>
-<TT>__GNUC__</TT>
+`__GNUC__`
    <TD ALIGN=LEFT NOWRAP>
 3
  <TR>
    <TD ALIGN=LEFT NOWRAP>
 GNU 3.2.1
    <TD ALIGN=LEFT NOWRAP>
-<TT>__GNUC_MINOR__</TT>
+`__GNUC_MINOR__`
    <TD ALIGN=LEFT NOWRAP>
 2
  <TR>
    <TD ALIGN=LEFT NOWRAP>
 GNU 3.2.1
    <TD ALIGN=LEFT NOWRAP>
-<TT>__GNUC_PATCHLEVEL__</TT>
+`__GNUC_PATCHLEVEL__`
    <TD ALIGN=LEFT NOWRAP>
 1
  <TR>
    <TD ALIGN=LEFT NOWRAP>
 Microsoft VC7.1
    <TD ALIGN=LEFT NOWRAP>
-<TT>_MSC_VER</TT>
+`_MSC_VER`
    <TD ALIGN=LEFT NOWRAP>
 1310
  <TR>
    <TD ALIGN=LEFT NOWRAP>
 Microsoft VC8.0
    <TD ALIGN=LEFT NOWRAP>
-<TT>_MSC_VER</TT>
+`_MSC_VER`
    <TD ALIGN=LEFT NOWRAP>
 1400
  <TR>
    <TD ALIGN=LEFT NOWRAP>
 Intel 11.1
    <TD ALIGN=LEFT NOWRAP>
-<TT>__INTEL_COMPILER</TT>
+`__INTEL_COMPILER`
    <TD ALIGN=LEFT NOWRAP>
 1110
  <TR>
    <TD ALIGN=LEFT NOWRAP>
 Clang 2.9
    <TD ALIGN=LEFT NOWRAP>
-<TT>__clang_major__</TT>
+`__clang_major__`
    <TD ALIGN=LEFT NOWRAP>
 2
  <TR>
    <TD ALIGN=LEFT NOWRAP>
 Clang 2.9
    <TD ALIGN=LEFT NOWRAP>
-<TT>__clang_minor__</TT>
+`__clang_minor__`
    <TD ALIGN=LEFT NOWRAP>
 9
  <TR>
    <TD ALIGN=LEFT NOWRAP>
 SUN 5.3
    <TD ALIGN=LEFT NOWRAP>
-<TT>__SUNPRO_CC</TT>
+`__SUNPRO_CC`
    <TD ALIGN=LEFT NOWRAP>
 0x530
  <TR>
    <TD ALIGN=LEFT NOWRAP>
 SUN 5.10
    <TD ALIGN=LEFT NOWRAP>
-<TT>__SUNPRO_CC</TT>
+`__SUNPRO_CC`
    <TD ALIGN=LEFT NOWRAP>
 0x5100
 </TABLE>
@ -178,17 +178,17 @@ There are also flags to identify the architecture.
    <TD ALIGN=LEFT NOWRAP>
 SGI
    <TD ALIGN=LEFT NOWRAP>
-<TT>__sgi</TT>
+`__sgi`
  <TR>
    <TD ALIGN=LEFT NOWRAP>
 SUN
    <TD ALIGN=LEFT NOWRAP>
-<TT>__sun</TT>
+`__sun`
  <TR>
    <TD ALIGN=LEFT NOWRAP>
 Linux
    <TD ALIGN=LEFT NOWRAP>
-<TT>__linux</TT>
+`__linux`
 </TABLE>
 </TABLE>

--- a/Documentation/doc/Documentation/Third_party.txt
+++ b/Documentation/doc/Documentation/Third_party.txt
@ -56,7 +56,7 @@ or <A HREF="https://msdn.microsoft.com/en-us/library/1fe2x6kt(v=vs.140).aspx">`h
 The \stl comes with the compiler, and as such no installation is required.

 \subsection thirdpartyBoost Boost
-<b>Version 1.62 or later</b>
+<b>Version 1.66 or later</b>

 The \sc{Boost} libraries are a set of portable C++ source libraries.
 Most of \sc{Boost} libraries are header-only, but a few of them need to be compiled or
--- a/Documentation/doc/resources/1.8.13/menu_version.js
+++ b/Documentation/doc/resources/1.8.13/menu_version.js
@ -1,14 +1,14 @@
 (function() {
  'use strict';

-  var url_re =  /(cgal\.geometryfactory\.com\/CGAL\/doc\/|doc\.cgal\.org\/)(master|latest|(\d\.\d+|\d\.\d+\.\d+))\//;
+  var url_re =  /(cgal\.geometryfactory\.com\/CGAL\/doc\/|doc\.cgal\.org\/)(master|latest|(\d\.\d+|\d\.\d+\.\d+)(-beta\d)?)\//;
  var url_local =  /.*\/doc_output\//;
-  var current_version_local = '5.1-beta1'
+  var current_version_local = '5.0.3'
  var all_versions = [
      'master',
      'latest',
-      '5.1-beta1',
-      '5.0.2',
+      '5.1',
+      '5.0.3',
      '4.14.3',
      '4.13.2',
      '4.12.2',
--- a/Documentation/doc/resources/1.8.13/stylesheet.css
+++ b/Documentation/doc/resources/1.8.13/stylesheet.css
@ -1245,10 +1245,12 @@ dl.citelist dt {
        font-weight:bold;
        margin-right:10px;
        padding:5px;
+        text-align:right;
+        width:52px;
 }

 dl.citelist dd {
-        margin:2px 0;
+        margin:2px 0 2px 72px;
        padding:5px 0;
 }

--- a/Documentation/doc/resources/1.8.14/menu_version.js
+++ b/Documentation/doc/resources/1.8.14/menu_version.js
@ -1,14 +1,14 @@
 (function() {
  'use strict';

-  var url_re =  /(cgal\.geometryfactory\.com\/CGAL\/doc\/|doc\.cgal\.org\/)(master|latest|(\d\.\d+|\d\.\d+\.\d+))\//;
+  var url_re =  /(cgal\.geometryfactory\.com\/CGAL\/doc\/|doc\.cgal\.org\/)(master|latest|(\d\.\d+|\d\.\d+\.\d+)(-beta\d)?)\//;
  var url_local =  /.*\/doc_output\//;
-  var current_version_local = '5.1-beta1'
+  var current_version_local = '5.0.3'
  var all_versions = [
      'master',
      'latest',
-      '5.1-beta1',
-      '5.0.2',
+      '5.1',
+      '5.0.3',
      '4.14.3',
      '4.13.2',
      '4.12.2',
--- a/Documentation/doc/resources/1.8.14/stylesheet.css
+++ b/Documentation/doc/resources/1.8.14/stylesheet.css
@ -1245,10 +1245,12 @@ dl.citelist dt {
        font-weight:bold;
        margin-right:10px;
        padding:5px;
+        text-align:right;
+        width:52px;
 }

 dl.citelist dd {
-        margin:2px 0;
+        margin:2px 0 2px 72px;
        padding:5px 0;
 }

--- a/Documentation/doc/resources/1.8.4/menu_version.js
+++ b/Documentation/doc/resources/1.8.4/menu_version.js
@ -1,14 +1,14 @@
 (function() {
  'use strict';

-  var url_re =  /(cgal\.geometryfactory\.com\/CGAL\/doc\/|doc\.cgal\.org\/)(master|latest|(\d\.\d+|\d\.\d+\.\d+))\//;
+  var url_re =  /(cgal\.geometryfactory\.com\/CGAL\/doc\/|doc\.cgal\.org\/)(master|latest|(\d\.\d+|\d\.\d+\.\d+)(-beta\d)?)\//;
  var url_local =  /.*\/doc_output\//;
-  var current_version_local = '5.1-beta1'
+  var current_version_local = '5.0.3'
  var all_versions = [
      'master',
      'latest',
-      '5.1-beta1',
-      '5.0.2',
+      '5.1',
+      '5.0.3',
      '4.14.3',
      '4.13.2',
      '4.12.2',
--- a/Documentation/doc/resources/1.8.4/stylesheet.css
+++ b/Documentation/doc/resources/1.8.4/stylesheet.css
@ -1079,10 +1079,12 @@ dl.citelist dt {
        font-weight:bold;
        margin-right:10px;
        padding:5px;
+        text-align:right;
+        width:52px;
 }

 dl.citelist dd {
-        margin:2px 0;
+        margin:2px 0 2px 72px;
        padding:5px 0;
 }

--- a/Filtered_kernel/include/CGAL/Lazy.h
+++ b/Filtered_kernel/include/CGAL/Lazy.h
@ -25,10 +25,8 @@
 #include <CGAL/Origin.h>
 #include <CGAL/Bbox_2.h>
 #include <CGAL/Bbox_3.h>
-#include <vector>
-#include <type_traits>
 #include <CGAL/Default.h>
-#include<CGAL/tss.h>
+#include <CGAL/tss.h>
 #include <CGAL/is_iterator.h>
 #include <CGAL/transforming_iterator.h>

@ -47,6 +45,14 @@
 #include <boost/preprocessor/repetition/enum_binary_params.hpp>
 #include <boost/preprocessor/repetition/enum.hpp>

+#include <iostream>
+#include <iterator>
+#include <string>
+#include <tuple>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
 namespace CGAL {

 template <class E,
--- a/Installation/CHANGES.md
+++ b/Installation/CHANGES.md
@ -1,10 +1,28 @@
 Release History
 ===============

+[Release 5.2](https://github.com/CGAL/cgal/releases/tag/v5.2)
+-----------
+
+Release date: December 2020
+
+### [CGAL and the Boost Graph Library (BGL)](https://doc.cgal.org/5.2/Manual/packages.html#PkgBGL)
+
+-   Added the convenience header `CGAL/boost/graph/graph_traits_inheritance_macros.h` that allows to easily
+    make any class inheriting from a model of a face graph concept, a model of the same concept.
+
+### [3D Convex Hulls](https://doc.cgal.org/5.2/Manual/packages.html#PkgConvexHull3)
+-   Added the function `CGAL::halfspace_intersection_interior_point_3()` that can be used to retrieve
+    the point that is the most interior a convex closed volume defined by the intersection of a set of halfspaces.
+
+### [Polygon Mesh Processing](https://doc.cgal.org/5.2/Manual/packages.html#PkgPolygonMeshProcessing)
+-   Added an option in `corefine()`, `split()` and `clip()` functions that enables the operation to be done on a mesh with
+    self-intersections present in the intersection area.
+
 [Release 5.1](https://github.com/CGAL/cgal/releases/tag/releases%2FCGAL-5.1)
 -----------

-Release date: July 2020
+Release date: September 2020

 ### [Tetrahedral Remeshing](https://doc.cgal.org/5.1/Manual/packages.html#PkgTetrahedralRemeshing) (new package)

@ -12,6 +30,8 @@ Release date: July 2020
    that improves the quality of tetrahedra in terms of dihedral angles,
    while targeting a given edge length.

+    See also the associated [blog entry](https://www.cgal.org/2020/08/07/Tetrahedral-remeshing/).
+
 ### [Surface Mesh Topology](https://doc.cgal.org/5.1/Manual/packages.html#PkgSurfaceMeshTopologySummary) (new package)

 -   This package enables the computation of some topological invariants of surfaces, such as:
@ -33,6 +53,7 @@ Release date: July 2020
 ### Installation

 -   The CGAL\_Core library no longer requires `Boost.Thread`, even if the g++ compiler is used.
+-   The minimal supported version of Boost is now 1.66.0.

 ### [Tutorials](https://doc.cgal.org/5.1/Manual/tutorials.html)

--- a/Installation/include/CGAL/Mesh_3/Mesh_complex_3_in_triangulation_3_fwd.h
+++ b/Installation/include/CGAL/Mesh_3/Mesh_complex_3_in_triangulation_3_fwd.h
@ -0,0 +1,38 @@
+// Copyright (C) 2020  GeometryFactory Sarl
+//
+// This file is part of CGAL (www.cgal.org)
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
+//
+
+#ifndef CGAL_MESH_3_MESH_COMPLEX_3_IN_TRIANGULATION_3_FWD_H
+#define CGAL_MESH_3_MESH_COMPLEX_3_IN_TRIANGULATION_3_FWD_H
+
+/// \file Mesh_complex_3_in_triangulation_3_fwd.h
+/// Forward declarations of the Mesh_3 package.
+
+#ifndef DOXYGEN_RUNNING
+namespace CGAL {
+
+// fwdS for the public interface
+template <typename Tr,
+          typename CornerIndex = int,
+          typename CurveIndex = int>
+class Mesh_complex_3_in_triangulation_3;
+
+template<class Tr, bool c3t3_loader_failed>
+bool build_triangulation_from_file(std::istream& is,
+                                   Tr& tr,
+                                   bool replace_domain_0);
+
+template<class Tr, bool c3t3_loader_failed>
+bool build_triangulation_from_file(std::istream& is,
+                                   Tr& tr);
+} // CGAL
+#endif
+
+#endif /* CGAL_MESH_3_MESH_COMPLEX_3_IN_TRIANGULATION_3_FWD_H */
+
+
--- a/Installation/include/CGAL/version.h
+++ b/Installation/include/CGAL/version.h
@ -17,12 +17,12 @@
 #define CGAL_VERSION_H

 #ifndef SWIG
-#define CGAL_VERSION 5.1-beta2
+#define CGAL_VERSION 5.2
 #define CGAL_GIT_HASH abcdef
 #endif
-#define CGAL_VERSION_NR 1050100000
+#define CGAL_VERSION_NR 1050200000
 #define CGAL_SVN_REVISION 99999
-#define CGAL_RELEASE_DATE 20200609
+#define CGAL_RELEASE_DATE 20200908

 #include <CGAL/version_macros.h>

--- a/Installation/lib/cmake/CGAL/CGALConfig.cmake
+++ b/Installation/lib/cmake/CGAL/CGALConfig.cmake
@ -150,6 +150,7 @@ foreach(cgal_lib ${CGAL_LIBRARIES})
  endif()
 endforeach()

+include(${CGAL_CONFIG_DIR}/CGALConfigVersion.cmake)

 #
 #
--- a/Installation/lib/cmake/CGAL/CGALConfigVersion.cmake
+++ b/Installation/lib/cmake/CGAL/CGALConfigVersion.cmake
@ -1,14 +1,8 @@
 set(CGAL_MAJOR_VERSION 5)
-set(CGAL_MINOR_VERSION 1)
+set(CGAL_MINOR_VERSION 2)
 set(CGAL_BUGFIX_VERSION 0)
 include(${CMAKE_CURRENT_LIST_DIR}/CGALConfigBuildVersion.cmake)
-if(NOT CMAKE_SCRIPT_MODE_FILE)
-  set(CGAL_MAJOR_VERSION  ${CGAL_MAJOR_VERSION} PARENT_SCOPE)
-  set(CGAL_MINOR_VERSION  ${CGAL_MINOR_VERSION} PARENT_SCOPE)
-  set(CGAL_BUGFIX_VERSION ${CGAL_BUGFIX_VERSION} PARENT_SCOPE)
-  set(CGAL_BUILD_VERSION  ${CGAL_BUILD_VERSION} PARENT_SCOPE)
-endif()
-set(CGAL_VERSION_PUBLIC_RELEASE_VERSION "5.1-beta2")
+set(CGAL_VERSION_PUBLIC_RELEASE_VERSION "5.2-dev")
 set(CGAL_VERSION_PUBLIC_RELEASE_NAME "CGAL-${CGAL_VERSION_PUBLIC_RELEASE_VERSION}")

 if (CGAL_BUGFIX_VERSION AND CGAL_BUGFIX_VERSION GREATER 0)
--- a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Bounded_3_do_intersect.h
+++ b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Bounded_3_do_intersect.h
@ -153,18 +153,11 @@ do_intersect(const typename K::Sphere_3 &sp,
 template <class K>
 inline
 typename K::Boolean
-do_intersect(const typename K::Tetrahedron_3 &tet,
-             const typename K::Tetrahedron_3 &sp,
+do_intersect(const typename K::Tetrahedron_3 &lh_tet,
+             const typename K::Tetrahedron_3 &rh_tet,
             const K & k)
 {
-  return do_intersect_tetrahedron_bounded(sp, tet, tet[0], k);
-}
-
-template <class K>
-inline typename K::Boolean do_intersect(const typename K::Tetrahedron_3 &tet,
-                                        const CGAL::Bbox_3 &bb, const K &k) {
-  // Swap arguments.
-  return do_intersect(bb, tet, k);
+  return do_intersect_tetrahedron_bounded(lh_tet, rh_tet, lh_tet[0], k);
 }

 // BBox_3 specific code since it is ok for BBox_3 to degenerate.
@ -197,6 +190,14 @@ inline typename K::Boolean do_intersect(const CGAL::Bbox_3 &aabb,
  return result;
 }

+
+template <class K>
+inline typename K::Boolean do_intersect(const typename K::Tetrahedron_3 &tet,
+                                        const CGAL::Bbox_3 &bb, const K &k) {
+  // Swap arguments.
+  return do_intersect(bb, tet, k);
+}
+
 } // namespace internal
 } // namespace Intersections
 } //namespace CGAL
--- a/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp
+++ b/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp
@ -0,0 +1,36 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Kernel/global_functions_3.h>
+#include <iostream>
+#include <utility>
+using K = CGAL::Simple_cartesian<double>;
+using Point_3 = K::Point_3;
+
+struct query {
+  Point_3 p;
+  double expected_angle;
+};
+
+int main() {
+  Point_3 a = {0, 0, 0};
+  Point_3 b = {0, 1, 0};
+  Point_3 c = {1, 0, 0};
+
+  const query queries[] = {
+    { {  1, 0,  0},    0.},
+    { {  1, 0,  1},   45.},
+    { {  0, 0,  1},   90.},
+    { { -1, 0,  1},  135.},
+    { { -1, 0,  0},  180.},
+    { { -1, 0, -1}, -135.},
+    { {  0, 0, -1},  -90.},
+    { {  1, 0, -1},  -45.},
+  };
+
+  for(auto query: queries) {
+    const auto& expected = query.expected_angle;
+    const auto& p = query.p;
+    auto approx = CGAL::approximate_dihedral_angle(a, b, c, p);
+    std::cout << approx << "  -- " << expected << '\n';
+    assert( std::abs(approx - expected) < 0.1 );
+  }
+};
--- a/Maintenance/infrastructure/cgal.geometryfactory.com/crontab
+++ b/Maintenance/infrastructure/cgal.geometryfactory.com/crontab
@ -23,13 +23,15 @@ 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 scl enable rh-git29 -- $HOME/bin/create_release $HOME/CGAL/branches/master.git --do-it --beta 2 --public || echo ERROR
+0 21  *  *  Sun                      cd $HOME/CGAL/create_internal_release && /usr/bin/time scl enable rh-git29 -- $HOME/bin/create_release $HOME/CGAL/branches/master.git --do-it || echo ERROR
 # "integration"
-0 21  *  *  Mon,Tue,Wed,Thu          cd $HOME/CGAL/create_internal_release && /usr/bin/time scl enable rh-git29 -- $HOME/bin/create_release $HOME/CGAL/branches/integration.git $HOME/CGAL/branches/empty-dir --do-it --beta 2 --public || echo ERROR
+0 21  *  *  Mon,Tue,Wed,Thu          cd $HOME/CGAL/create_internal_release && /usr/bin/time scl enable rh-git29 -- $HOME/bin/create_release $HOME/CGAL/branches/integration.git $HOME/CGAL/branches/empty-dir --do-it || echo ERROR
+# from branch 5.1
+0 21  *  *  Sat                      cd $HOME/CGAL/create_internal_release-5.1-branch && /usr/bin/time scl enable rh-git29 -- $HOME/bin/create_release $HOME/CGAL/branches/CGAL-5.1-branch.git --public --do-it || echo ERROR
 # from branch 5.0
-0 21  *  *  Sat                      cd $HOME/CGAL/create_internal_release-5.0-branch && /usr/bin/time scl enable rh-git29 -- $HOME/bin/create_release $HOME/CGAL/branches/CGAL-5.0-branch.git --public --do-it || echo ERROR
+0 21  *  *  Fri                      cd $HOME/CGAL/create_internal_release-5.0-branch && /usr/bin/time scl enable rh-git29 -- $HOME/bin/create_release $HOME/CGAL/branches/CGAL-5.0-branch.git --public --do-it || echo ERROR
 # from branch 4.14
-0 21  *  *  Fri                      cd $HOME/CGAL/create_internal_release-4.14-branch && /usr/bin/time scl enable rh-git29 -- $HOME/bin/create_release $HOME/CGAL/branches/CGAL-4.14-branch.git $HOME/CGAL/branches/empty-dir --public --do-it || echo ERROR
+#0 21  *  *  Fri                      cd $HOME/CGAL/create_internal_release-4.14-branch && /usr/bin/time scl enable rh-git29 -- $HOME/bin/create_release $HOME/CGAL/branches/CGAL-4.14-branch.git --public --do-it || echo ERROR

 ## Older stuff
 # from branch 4.13
--- a/Maintenance/public_release/announcement/mailing-beta.eml
+++ b/Maintenance/public_release/announcement/mailing-beta.eml
@ -1,14 +1,14 @@
-Subject: CGAL 5.1 Beta 1 Released, Computational Geometry Algorithms Library
+Subject: CGAL 5.1 Beta 2 Released, Computational Geometry Algorithms Library
 Content-Type: text/plain; charset="utf-8"
 Body:

-The CGAL Open Source Project is pleased to announce the release 5.1 Beta 1
+The CGAL Open Source Project is pleased to announce the release 5.1 Beta 2
 of CGAL, the Computational Geometry Algorithms Library.


-CGAL version 5.1 Beta 1 is a public testing release. It should provide a
+CGAL version 5.1 Beta 2 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 5.1 in July.
+of the final version of CGAL 5.1 in September.


 Besides fixes and general enhancement to existing packages, the following
@ -45,6 +45,11 @@ Optimal Bounding Box (new package)
    See also the associated blog entry:
      https://www.cgal.org/2020/04/20/Optimal_bounding_box/

+Installation
+
+-   The CGAL_Core library no longer requires Boost.Thread, even if the
+    g++ compiler is used.
+
 Tutorials

 -   Two new, detailed tutorials have been added:
@ -118,7 +123,7 @@ CGAL and the Boost Graph Library (BGL)
    uses this alpha expansion graphcut to regularize the borders of a
    selected faces on a triangle mesh.

-See https://www.cgal.org/2020/06/09/cgal51-beta1/ for a complete list of
+See https://www.cgal.org/2020/07/28/cgal51-beta2/ for a complete list of
 changes.


--- a/Maintenance/public_release/announcement/mailing.eml
+++ b/Maintenance/public_release/announcement/mailing.eml
@ -1,51 +1,129 @@
-Subject: CGAL 4.14 Released, Computational Geometry Algorithms Library
+Subject: CGAL 5.1 Released, Computational Geometry Algorithms Library
 Content-Type: text/plain; charset="utf-8"
-Body: 
+Body:

-The CGAL Open Source Project is pleased to announce the release 4.14
+The CGAL Open Source Project is pleased to announce the release 5.1
 of CGAL, the Computational Geometry Algorithms Library.

-This release will be the last one supporting C++98/C++03 compilers,
-and the next release, CGAL version 5.0, will require support for
-C++14.

-Besides fixes and general enhancement to existing packages, the
-following has changed since CGAL 4.13:
+Besides fixes and general enhancement to existing packages, the following
+has changed since CGAL 5.0:


-### 2D Periodic Hyperbolic Triangulations (new package)
+Tetrahedral Remeshing (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.
+-   This package implements a tetrahedral isotropic remeshing algorithm,
+    that improves the quality of tetrahedra in terms of dihedral angles,
+    while targeting a given edge length.

-### 2D Hyperbolic Triangulations (new package)
+    See also the associated blog entry:
+      https://www.cgal.org/2020/08/07/Tetrahedral-remeshing/

-   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.
+Surface Mesh Topology (new package)

-### The Heat Method (new package)
+-   This package enables the computation of some topological invariants
+    of surfaces, such as:
+    -   test if two (closed) curves on a combinatorial surface are
+        homotopic. Users can choose between free homotopy and homotopy
+        with fixed endpoints;
+    -   test is a curve is contractible;
+    -   compute shortest non-contractible cycles on a surface, with or
+        without weights on edges.

-   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.
+    See also the associated blog entry:
+      https://www.cgal.org/2020/05/08/Surface_mesh_topology/

-### Triangulated Surface Mesh Approximation (new package)
+Optimal Bounding Box (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 implements an optimization algorithm that aims to
+    construct a close approximation of the _optimal bounding box_ of a
+    mesh or a point set, which is defined as the smallest (in terms of
+    volume) bounding box that contains a given mesh or point set.
+
+    See also the associated blog entry:
+      https://www.cgal.org/2020/04/20/Optimal_bounding_box/
+
+Installation
+
+-   The CGAL_Core library no longer requires Boost.Thread, even if the
+    g++ compiler is used.
+
+-   The minimal supported version of Boost is now 1.66.0.
+
+Tutorials
+
+-   Two new, detailed tutorials have been added:
+    -   Surface Reconstruction from Point Clouds, which goes over a
+        typical full processing pipeline in a CGAL environment.
+    -   Geographic Information Systems (GIS), which demonstrates usage
+        of CGAL data structures and algorithms in the context of a
+        typical GIS application.
+
+    Both tutorials provide complete code.
+
+    See https://doc.cgal.org/5.1/Manual/tutorials.html
+
+Point Set Processing
+
+-   Added wrapper functions for registration, using the Super4PCS and
+    ICP algorithms implemented in the third party libraries OpenGR and
+    libpointmatcher.


-See https://www.cgal.org/2019/03/29/cgal414/ for a complete list of
+Surface Mesh Simplification
+
+-   Added a new simplification method based on the quadric error defined
+    by Garland and Heckbert.
+
+dD Spatial Searching
+
+-   The kd-tree can now be built in parallel: CGAL::Kd_tree::build() is
+    given an optional template parameter ConcurrencyTag (default value
+    remains CGAL::Sequential_tag for backward compatibility).
+
+Intersecting Sequences of dD Iso-oriented Boxes
+
+-   Added parallel versions of the functions CGAL::box_intersection_d()
+    and CGAL::box_self_intersection_d().
+
+Polygon Mesh Processing
+
+-   Added the function CGAL::Polygon_mesh_processing::split(), which can
+    be used to split meshes along a mesh or a plane.
+-   Added the function
+    CGAL::Polygon_mesh_processing::split_connected_components() to split
+    a single mesh containing several connected components into several
+    meshes containing one connected component.
+-   Added parallel versions of the functions
+    CGAL::Polygon_mesh_processing::does_self_intersect() and
+    CGAL::Polygon_mesh_processing::self_intersections().
+-   Added several mesh repair functions (see the complete changelog for
+    more information).
+
+3D Fast Intersection and Distance Computation
+
+-   The behavior of the internal search tree used to accelerate distance
+    queries has changed: usage of the internal search tree will now be
+    enabled by default, and its construction will be triggered by the
+    first distance query. Automatic construction and usage can be
+    disabled by calling
+    CGAL::AABB_tree::do_not_accelerate_distance_queries() before the
+    first distance query, and the tree can be built at any moment by
+    calling CGAL::AABB_tree::accelerate_distance_queries().
+-   BREAKING CHANGE: CGAL::AABB_tree::accelerate_distance_queries() and
+    CGAL::AABB_tree::do_not_accelerate_distance_queries() are no longer
+    const functions.
+
+
+CGAL and the Boost Graph Library (BGL)
+
+-   Added the function CGAL::alpha_expansion_graphcut(), which
+    regularizes a multi-label partition over a user-defined graph.
+-   Added the function CGAL::regularize_face_selection_borders(), which
+    uses this alpha expansion graphcut to regularize the borders of a
+    selected faces on a triangle mesh.
+
+See https://www.cgal.org/2020/09/08/cgal51/ for a complete list of
 changes.


--- a/Maintenance/public_release/scripts/github-release-upload
+++ b/Maintenance/public_release/scripts/github-release-upload
@ -9,5 +9,5 @@
 source $HOME/private/github-token.txt

 for f in *.(zip|xz|txt|exe); do
-    github-release upload -u CGAL -r cgal -t releases/${PWD:t} --name $f -f $f;
+    github-release upload -u CGAL -r cgal -t v${${PWD:t}/CGAL-/} --name $f -f $f;
 done
--- a/Maintenance/public_release/scripts/prepare_release
+++ b/Maintenance/public_release/scripts/prepare_release
@ -3,12 +3,14 @@
 set -e
 setopt -o glob
 setopt -o extended_glob
+setopt -o nullglob

 if [ -z "$1" ]; then
  printf "Usage:\n  prepare_release <internal_release_public_dir>\n"
  exit 1
 fi

+error_code=0

 RELEASE_CANDIDATES_DIR=/srv/CGAL/www/Members/release-candidates
 MANUAL_TESTS_DIR=/srv/CGAL/www/Members/Manual_doxygen_test
@ -25,7 +27,7 @@ PUBLIC_RELEASE_NAME=${PUBLIC_RELEASE_NAME/-library.zip/}
 PUBLIC_RELEASE_NAME=${PUBLIC_RELEASE_NAME/.zip/}
 DEST_DIR="${RELEASE_CANDIDATES_DIR}/$PUBLIC_RELEASE_NAME"

-if [ -z "$PUBLIC_RELEASE_NAME" ]; then 
+if [ -z "$PUBLIC_RELEASE_NAME" ]; then
    echo The first argument must be a public release directory, like CGAL-3.7-I-167-public.
    exit 1
 fi
@ -40,30 +42,36 @@ printf "Copy documentation to doc_html/ and doc_html_online/...\n"
 [ -d "/srv/CGAL/www/${PUBLIC_RELEASE_NAME#CGAL-}/Manual" ] || mkdir -p "/srv/CGAL/www/${PUBLIC_RELEASE_NAME#CGAL-}/Manual"

 cp "$PUBLIC_RELEASE_DIR"/*(.) "${RELEASE_CANDIDATES_DIR}/$PUBLIC_RELEASE_NAME"
-rsync -a "$MANUAL_TESTS_DIR/$INTERNAL_RELEASE"/output2/* "$DEST_DIR/doc_html/"
-pushd "$DEST_DIR/doc_html/Manual/search"
-  for i in g n c s i; do sed -i "s/..\/BGL$i/..\/BGL\/$i/g" *; done
-popd
-rsync -a "$DEST_DIR/doc_html/" "$DEST_DIR/doc_html_online/"
+files=("$MANUAL_TESTS_DIR/$INTERNAL_RELEASE"/output2/*)
+if ((${#files[@]} == 0)); then
+    printf "ERROR: documentation files are missing\n"
+    error_code=1
+else
+    rsync -a "$MANUAL_TESTS_DIR/$INTERNAL_RELEASE"/output2/* "$DEST_DIR/doc_html/"
+    pushd "$DEST_DIR/doc_html/Manual/search"
+    for i in g n c s i; do sed -i "s/..\/BGL$i/..\/BGL\/$i/g" *; done
+    popd
+    rsync -a "$DEST_DIR/doc_html/" "$DEST_DIR/doc_html_online/"

-printf "Modify \"../../MathJax\" to MathJax CDN in doc_html...\n"
-pushd "$DEST_DIR/doc_html"
-find -name '*.html' -exec sed -i -e 's|../../MathJax|https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.2|g' {} +
-popd
+    printf "Modify \"../../MathJax\" to MathJax CDN in doc_html...\n"
+    pushd "$DEST_DIR/doc_html"
+    find -name '*.html' -exec sed -i -e 's|../../MathJax|https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.2|g' {} +
+    popd

-printf "Create doc_html tarball...\n"
-tar cf "$DEST_DIR/${PUBLIC_RELEASE_NAME}-doc_html.tar" -C "$DEST_DIR" doc_html
+    printf "Create doc_html tarball...\n"
+    tar cf "$DEST_DIR/${PUBLIC_RELEASE_NAME}-doc_html.tar" -C "$DEST_DIR" doc_html

-printf "Create doc_html zip...\n"
-pushd "$DEST_DIR"
-zip -q -r "$DEST_DIR/${PUBLIC_RELEASE_NAME}-doc_html.zip" doc_html
-popd
+    printf "Create doc_html zip...\n"
+    pushd "$DEST_DIR"
+    zip -q -r "$DEST_DIR/${PUBLIC_RELEASE_NAME}-doc_html.zip" doc_html
+    popd

-printf "xz doc_html tarball...\n"
-xz --best < "$DEST_DIR/${PUBLIC_RELEASE_NAME}-doc_html.tar" > "$DEST_DIR/${PUBLIC_RELEASE_NAME}-doc_html.tar.xz"
+    printf "xz doc_html tarball...\n"
+    xz --best < "$DEST_DIR/${PUBLIC_RELEASE_NAME}-doc_html.tar" > "$DEST_DIR/${PUBLIC_RELEASE_NAME}-doc_html.tar.xz"

-printf "gzip doc_html tarball...\n"
-gzip -f --best "$DEST_DIR/${PUBLIC_RELEASE_NAME}-doc_html.tar"
+    printf "gzip doc_html tarball...\n"
+    gzip -f --best "$DEST_DIR/${PUBLIC_RELEASE_NAME}-doc_html.tar"
+fi

 printf "Create the internal_release file\n"
 printf "%s\n" "${INTERNAL_RELEASE}" > "$DEST_DIR/internal_release"
@ -92,3 +100,4 @@ printf "DONE!\n"
 pushd "$DEST_DIR"
 ls -l
 popd
+exit $error_code
--- a/Mesh_3/include/CGAL/IO/File_medit.h
+++ b/Mesh_3/include/CGAL/IO/File_medit.h
@ -17,6 +17,7 @@
 #include <CGAL/license/Triangulation_3.h>

 #include <CGAL/Mesh_3/config.h>
+#include <CGAL/Mesh_3/Mesh_complex_3_in_triangulation_3_fwd.h>

 #include <CGAL/utility.h>
 #include <CGAL/basic.h>
@ -893,6 +894,23 @@ output_to_medit(std::ostream& os,
  }
 }

+
+template<typename T3>
+void write_MEDIT(std::ostream& os, const T3& t3)
+{
+  CGAL::Mesh_complex_3_in_triangulation_3<T3, int, int> c3t3;
+  c3t3.triangulation() = t3;
+  c3t3.rescan_after_load_of_triangulation();
+  output_to_medit(os, c3t3);
+}
+
+template<typename T3>
+bool read_MEDIT(std::istream& in, T3& t3)
+{
+  CGAL_assertion(!(!in));
+  return CGAL::build_triangulation_from_file<T3, true>(in, t3);
+}
+
 } // end namespace CGAL

 #endif // CGAL_IO_FILE_MEDIT_H
--- a/Mesh_3/include/CGAL/Mesh_3/tet_soup_to_c3t3.h
+++ b/Mesh_3/include/CGAL/Mesh_3/tet_soup_to_c3t3.h
@ -106,7 +106,8 @@ bool build_finite_cells(Tr& tr,
    boost::unordered_map<std::array<typename Tr::Vertex_handle, 3>,
                        std::vector<std::pair<typename Tr::Cell_handle, int> > >& incident_cells_map,
    const std::map<std::array<int,3>, typename Tr::Cell::Surface_patch_index>& border_facets,
-    const bool verbose)
+    const bool verbose,
+                        bool replace_domain_0 = false)
 {
  typedef std::array<int, 5>              Tet_with_ref; // 4 ids + 1 reference

@ -118,7 +119,16 @@ bool build_finite_cells(Tr& tr,
    typename Tr::Geom_traits::Construct_point_3 cp =
      tr.geom_traits().construct_point_3_object();
  )
-
+  int max_domain = 0;
+  if(replace_domain_0)
+  {
+    for(std::size_t i=0; i<finite_cells.size(); ++i)
+    {
+      const Tet_with_ref& tet = finite_cells[i];
+      if(tet[4] > max_domain)
+        max_domain=tet[4];
+    }
+  }
  // build the finite cells
  for(std::size_t i=0; i<finite_cells.size(); ++i)
  {
@ -141,7 +151,10 @@ bool build_finite_cells(Tr& tr,

    Cell_handle c = tr.tds().create_cell(vs[0], vs[1], vs[2], vs[3]);
    c->set_subdomain_index(tet[4]); // the cell's info keeps the reference of the tetrahedron
-
+    if(replace_domain_0 && tet[4] == 0)
+    {
+      c->set_subdomain_index(max_domain+1); // the cell's info keeps the reference of the tetrahedron
+    }
    // assign cells to vertices
    for(int j=0; j<4; ++j)
    {
@ -328,7 +341,8 @@ bool build_triangulation(Tr& tr,
                         const std::vector<std::array<int,5> >& finite_cells,
                         const std::map<std::array<int,3>, typename Tr::Cell::Surface_patch_index>& border_facets,
                         std::vector<typename Tr::Vertex_handle>& vertex_handle_vector,
-                         const bool verbose = false)
+                         const bool verbose = false,
+                         bool replace_domain_0 = false)
 {
  typedef typename Tr::Vertex_handle            Vertex_handle;
  typedef typename Tr::Cell_handle              Cell_handle;
@ -360,7 +374,7 @@ bool build_triangulation(Tr& tr,
  if (!finite_cells.empty())
  {
    if(!build_finite_cells<Tr>(tr, finite_cells, vertex_handle_vector, incident_cells_map,
-                           border_facets, verbose))
+                           border_facets, verbose, replace_domain_0))
      return false;
    if(!build_infinite_cells<Tr>(tr, incident_cells_map, verbose))
      return false;
@ -384,7 +398,8 @@ bool build_triangulation(Tr& tr,

 template<class Tr, bool c3t3_loader_failed>
 bool build_triangulation_from_file(std::istream& is,
-                                   Tr& tr)
+                                   Tr& tr,
+                                   bool replace_domain_0)
 {
  typedef typename Tr::Point                                  Point_3;

@ -473,10 +488,16 @@ bool build_triangulation_from_file(std::istream& is,

  std::vector<typename Tr::Vertex_handle> vertices(points.size() + 1);
  bool is_well_built = build_triangulation<Tr, c3t3_loader_failed>(tr,
-    points, finite_cells, border_facets, vertices);
+    points, finite_cells, border_facets, vertices, false, replace_domain_0);
  return is_well_built;
 }

+template<class Tr, bool c3t3_loader_failed>
+bool build_triangulation_from_file(std::istream& is,
+                                   Tr& tr)
+{
+  return build_triangulation_from_file<Tr, c3t3_loader_failed>(is, tr, false);
+}
 }  // namespace CGAL

 #include <CGAL/enable_warnings.h>
--- a/Mesh_3/include/CGAL/Mesh_complex_3_in_triangulation_3.h
+++ b/Mesh_3/include/CGAL/Mesh_complex_3_in_triangulation_3.h
@ -20,6 +20,7 @@

 #include <CGAL/license/Triangulation_3.h>

+#include <CGAL/Mesh_3/Mesh_complex_3_in_triangulation_3_fwd.h>
 #include <CGAL/disable_warnings.h>
 #include <CGAL/iterator.h>
 #include <CGAL/Mesh_3/utilities.h>
@ -38,8 +39,8 @@ namespace CGAL {


 template <typename Tr,
-          typename CornerIndex = int,
-          typename CurveIndex = int>
+          typename CornerIndex,
+          typename CurveIndex>
 class Mesh_complex_3_in_triangulation_3 :
  public Mesh_3::Mesh_complex_3_in_triangulation_3_base<
    Tr, typename Tr::Concurrency_tag>
--- a/Nef_3/include/CGAL/Nef_3/Bounding_box_3.h
+++ b/Nef_3/include/CGAL/Nef_3/Bounding_box_3.h
@ -73,17 +73,13 @@ public:
  }

  void extend(const Point_3& p) {
-    FT q[3];
-    q[0] = p.x();
-    q[1] = p.y();
-    q[2] = p.z();
+    FT q[3] = { p.x(), p.y(), p.z() };

    if(initialized)
      Base::extend(q);
    else {
      initialized = true;
-      std::copy( q, q + 3, Base::lo );
-      std::copy( q, q + 3, Base::hi );
+      *this = Bounding_box_3(q);
    }
  }
 };
--- a/Nef_3/include/CGAL/Nef_3/SNC_constructor.h
+++ b/Nef_3/include/CGAL/Nef_3/SNC_constructor.h
@ -229,9 +229,10 @@ public:
    Vertex_handle v = this->sncp()->new_vertex(p , boundary);
    CGAL_NEF_TRACEN( v->point());
    SM_decorator SD(&*v);
-    Sphere_point sp[] = { Sphere_point(NT(-x), 0, 0),
-                          Sphere_point(0, NT(-y), 0),
-                          Sphere_point(0, 0, NT(-z)) };
+    const NT zero(0);
+    Sphere_point sp[] = { Sphere_point(NT(-x), zero, zero),
+                          Sphere_point(zero, NT(-y), zero),
+                          Sphere_point(zero, zero, NT(-z)) };

  /* create box vertices */
    SVertex_handle sv[3];
--- a/Nef_3/include/CGAL/Nef_3/SNC_io_parser.h
+++ b/Nef_3/include/CGAL/Nef_3/SNC_io_parser.h
@ -56,7 +56,8 @@ namespace Nef_3_internal{

 BOOST_MPL_HAS_XXX_TRAIT_NAMED_DEF(Has_nested_Exact_kernel,Exact_kernel,false)

-template <class R,bool has_exact_kernel=Has_nested_Exact_kernel<R>::value, class FT = typename R::RT, class Kernel_tag=typename R::Kernel_tag>
+template <class R,bool has_exact_kernel=Has_nested_Exact_kernel<R>::value &&
+                                        !std::is_floating_point<typename R::RT>::value, class FT = typename R::RT, class Kernel_tag=typename R::Kernel_tag>
 struct Type_converter{
  typedef const CGAL::Point_3<R>& Point_3;
  typedef const CGAL::Vector_3<R>& Vector_3;
@ -670,8 +671,9 @@ template<typename Tag, typename Kernel> class Geometry_io;
 template<typename Kernel>
 class Geometry_io<Cartesian_tag, Kernel> {
 public:
-  template <typename EK, typename K> static
-  typename EK::Point_3 read_point(std::istream& in) {
+
+  template <typename EK, typename K, typename Compose_> static
+  typename EK::Point_3 read_point_impl(std::istream& in, Compose_) {
    typedef Fraction_traits<typename K::FT> FracTraits;
    typename FracTraits::Type hx, hy, hz, hw;
    typename FracTraits::Numerator_type num;
@ -688,8 +690,8 @@ class Geometry_io<Cartesian_tag, Kernel> {
    return typename EK::Point_3(hx,hy,hz,hw);
  }

-  template <typename EK, typename K> static
-  typename EK::Plane_3 read_plane(std::istream& in) {
+  template <typename EK, typename K, typename Compose_> static
+  typename EK::Plane_3 read_plane_impl(std::istream& in, Compose_) {
    typedef Fraction_traits<typename K::FT> FracTraits;
    typename FracTraits::Type a, b, c, d;
    typename FracTraits::Numerator_type num;
@ -706,8 +708,32 @@ class Geometry_io<Cartesian_tag, Kernel> {
    return typename EK::Plane_3(a,b,c,d);
  }

-  template <typename R> static
-    void print_point_impl(std::ostream& out, const CGAL::Point_3<R> p) {
+  template <typename EK, typename K> static
+  typename EK::Point_3 read_point_impl(std::istream& in, Null_functor) {
+    typename K::FT hx, hy, hz, hw;
+    in >> hx >> hy >> hz >> hw;
+    return typename EK::Point_3(hx,hy,hz,hw);
+  }
+
+  template <typename EK, typename K> static
+  typename EK::Plane_3 read_plane_impl(std::istream& in, Null_functor) {
+    typename K::FT a, b, c, d;
+    in >> a >> b >> c >> d;
+    return typename EK::Plane_3(a,b,c,d);
+  }
+
+  template <typename EK, typename K> static
+  typename EK::Point_3 read_point(std::istream& in) {
+    return read_point_impl<EK,K>(in, typename Fraction_traits<typename K::FT>::Compose());
+  }
+
+  template <typename EK, typename K> static
+  typename EK::Plane_3 read_plane(std::istream& in) {
+    return read_plane_impl<EK,K>(in, typename Fraction_traits<typename K::FT>::Compose());
+  }
+
+  template <typename R, typename Decompose_> static
+  void print_point_impl(std::ostream& out, const CGAL::Point_3<R> p, Decompose_) {
    typedef Fraction_traits<typename R::FT> FracTraits;
    typedef std::vector<typename FracTraits::Numerator_type> NV;

@ -737,8 +763,8 @@ class Geometry_io<Cartesian_tag, Kernel> {
        << vec[2] << " " << vec[3];
  }

-  template <typename R> static
-    void print_vector_impl(std::ostream& out, const CGAL::Vector_3<R> p) {
+  template <typename R, typename Decompose_> static
+  void print_vector_impl(std::ostream& out, const CGAL::Vector_3<R> p, Decompose_) {
    typedef Fraction_traits<typename R::FT> FracTraits;
    typedef typename FracTraits::Numerator_type NumType;
    typedef std::vector<NumType> NV;
@ -766,8 +792,8 @@ class Geometry_io<Cartesian_tag, Kernel> {
        << vec[2] << " " << NumType(1);
  }

-  template <typename R> static
-  void print_plane_impl(std::ostream& out, const CGAL::Plane_3<R> p) {
+  template <typename R, typename Decompose_> static
+  void print_plane_impl(std::ostream& out, const CGAL::Plane_3<R> p, Decompose_) {
    typedef Fraction_traits<typename R::FT> FracTraits;
    typedef std::vector<typename FracTraits::Numerator_type> NV;

@ -803,19 +829,37 @@ class Geometry_io<Cartesian_tag, Kernel> {
        << vec[2] << " " << vec[3];
  }

+  template <typename R> static
+  void print_point_impl(std::ostream& out, const CGAL::Point_3<R> p, Null_functor)
+  {
+    out << p.x() << " " << p.y() << " " << p.z() << " " << 1;
+  }
+
+  template <typename R> static
+  void print_vector_impl(std::ostream& out, const CGAL::Vector_3<R> v, Null_functor)
+  {
+    out << v.x() << " " << v.y() << " " << v.z() << " " << 1;
+  }
+
+  template <typename R> static
+  void print_plane_impl(std::ostream& out, const CGAL::Plane_3<R> p, Null_functor)
+  {
+    out << p.a() << " " << p.b() << " " << p.c() << " " << p.d();
+  }
+
  template <class R> static
  void print_point(std::ostream& out, const CGAL::Point_3<R>& p) {
-    print_point_impl(out, Nef_3_internal::get_point(p) );
+    print_point_impl(out, Nef_3_internal::get_point(p), typename Fraction_traits<typename R::FT>::Decompose() );
  }

  template <class R> static
  void print_vector(std::ostream& out, const CGAL::Vector_3<R>& v) {
-    print_vector_impl(out, Nef_3_internal::get_vector(v) );
+    print_vector_impl(out, Nef_3_internal::get_vector(v), typename Fraction_traits<typename R::FT>::Decompose() );
  }

  template <class R> static
  void print_plane(std::ostream& out, const CGAL::Plane_3<R>& p) {
-    print_plane_impl(out, Nef_3_internal::get_plane(p) );
+    print_plane_impl(out, Nef_3_internal::get_plane(p), typename Fraction_traits<typename R::FT>::Decompose() );
  }

 };
--- a/Nef_3/include/CGAL/Nef_3/SNC_point_locator.h
+++ b/Nef_3/include/CGAL/Nef_3/SNC_point_locator.h
@ -563,7 +563,7 @@ public:
    Partial_facet pf;
 #endif
    bool hit = false;
-    Point_3 eor; // 'end of ray', the latest ray's hit point
+    Point_3 eor = CGAL::ORIGIN; // 'end of ray', the latest ray's hit point
    Objects_along_ray objects = candidate_provider->objects_along_ray(ray);
    Objects_along_ray_iterator objects_iterator = objects.begin();
    while( !hit && objects_iterator != objects.end()) {
--- a/Nef_3/test/Nef_3/CMakeLists.txt
+++ b/Nef_3/test/Nef_3/CMakeLists.txt
@ -24,6 +24,7 @@ if ( CGAL_FOUND )
        execution___of__test_nef_3_io_Cartesian
        execution___of__test_nef_3_io_Cartesian_Lazy
        execution___of__test_nef_3_io_EPEC
+        execution___of__test_nef_3_io_EPIC
        execution___of__test_nef_3_io_Homogeneous
        execution___of__test_nef_3_io_Homogenoeus_Lazy
        execution___of__test_with_extended_homogeneous
--- a/Nef_3/test/Nef_3/test_nef_3_io_EPIC.cpp
+++ b/Nef_3/test/Nef_3/test_nef_3_io_EPIC.cpp
@ -0,0 +1,57 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Nef_polyhedron_3.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/IO/Nef_polyhedron_iostream_3.h>
+
+#include <fstream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel EPEC;
+
+template <class Kernel>
+void test_write_read()
+{
+  typedef CGAL::Nef_polyhedron_3< Kernel > Nef_polyhedron;
+  typedef CGAL::Polyhedron_3< Kernel >     Polyhedron;
+  typedef typename Kernel::Point_3         Point;
+
+  typename Kernel::RT n, d;
+  std::istringstream str_n("5");
+  str_n >> n;
+  std::istringstream str_d("1");
+  str_d >> d;
+
+  Point p(n, 0, 0, d);
+  Point q(0, n, 0, d);
+  Point r(0, 0, n, d);
+  Point s(0, 0, 0, 1);
+
+  std::cout << "    build...\n";
+  Polyhedron P;
+  P.make_tetrahedron( p, q, r, s);
+  Nef_polyhedron nef_1( P );
+
+  std::cout << "    write...\n";
+  std::ofstream out ("temp.nef");
+  out << std::setprecision(17) << nef_1;
+  out.close();
+
+  std::cout << "    read...\n";
+  std::ifstream in ("temp.nef");
+  Nef_polyhedron nef_2;
+  in >> nef_2;
+  in.close();
+
+  std::cout << "    check...\n";
+  assert( nef_1 == nef_2);
+}
+
+int main()
+{
+  std::cout << "Testing Exact_predicates_inexact_constructions_kernel\n";
+  test_write_read<EPEC>();
+  std::cout << "Testing Simple_cartesian<double>\n";
+  test_write_read<CGAL::Simple_cartesian<double>>();
+
+  return 0;
+}
--- a/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_functors.h
+++ b/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_functors.h
@ -224,6 +224,7 @@ template<class R_> struct Scalar_product {
 };

 template<class R_> struct Squared_distance_to_origin_stored {
+        // What about weighted points, should they store sdo-w?
        CGAL_FUNCTOR_INIT_IGNORE(Squared_distance_to_origin_stored)
        typedef R_ R;
        typedef typename R::LA_vector LA;
--- a/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h
+++ b/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h
@ -149,6 +149,7 @@ template <class R_> struct Power_center : Store_kernel<R_> {
  typedef typename Get_type<R_, Weighted_point_tag>::type        WPoint;
  typedef WPoint result_type;
  typedef typename Get_type<R_, Point_tag>::type        Point;
+  typedef typename Get_type<R_, Vector_tag>::type       Vector;
  typedef typename Get_type<R_, FT_tag>::type FT;
  template <class Iter>
  result_type operator()(Iter f, Iter e)const{
@ -168,11 +169,26 @@ template <class R_> struct Power_center : Store_kernel<R_> {

    WPoint const& wp0 = *f;
    Point const& p0 = pdw(wp0);
+    FT const& w0 = pw(wp0);
    int d = pd(p0);
    int k = static_cast<int>(std::distance(f,e));
-    if (d+1 == k)
-    {
-      FT const& n0 = sdo(p0) - pw(wp0);
+    if (k == 1) return cwp(p0, -w0);
+    // TODO: check for degenerate cases?
+    if (k == 2) {
+      typename Get_functor<R_, Difference_of_points_tag>::type dp(this->kernel());
+      typename Get_functor<R_, Squared_length_tag>::type sl(this->kernel());
+      typename Get_functor<R_, Translated_point_tag>::type tp(this->kernel());
+      typename Get_functor<R_, Scaled_vector_tag>::type sv(this->kernel());
+      WPoint const& wp1 = *++f;
+      Point const& p1 = pdw(wp1);
+      FT const& w1 = pw(wp1);
+      Vector v01 = dp(p1, p0);
+      FT l01 = sl(v01);
+      FT coef = ((w0 - w1) / l01 + 1) / 2;
+      return cwp(tp(p0, sv(v01, coef)), CGAL::square(coef) * l01 - w0);
+    }
+    if (d+1 == k) {
+      FT const& n0 = sdo(p0) - w0;
      Matrix m(d,d);
      Vec b = typename CVec::Dimension()(d);
      // Write the point coordinates in lines.
@ -194,64 +210,54 @@ template <class R_> struct Power_center : Store_kernel<R_> {
      FT const& r2 = pdp (wp0, center);
      return cwp(std::move(center), r2);
    }
-    else
+
    {
-      /*
-       * Matrix P=(p1, p2, ...) (each point as a column)
-       * Matrix Q=2*t(p2-p1,p3-p1, ...) (each vector as a line)
-       * Matrix M: QP, adding a line of 1 at the top
-       * Vector B: (1, p2^2-p1^2, p3^2-p1^2, ...) plus weights
-       * Solve ML=B, the center of the sphere is PL
-       *
-       * It would likely be faster to write P then transpose, multiply,
-       * etc instead of doing it by hand.
-       */
-      // TODO: check for degenerate cases?
-
-      typedef typename R_::Max_ambient_dimension D2;
-      typedef typename R_::LA::template Rebind_dimension<Dynamic_dimension_tag,D2>::Other LAd;
-      typedef typename LAd::Square_matrix Matrix;
-      typedef typename LAd::Vector Vec;
+      // The general case. ui=p(i+1)-p0, center-p0=c=sum ai*ui, c.2ui=ui²+w-w, M*a=b with M symmetric
+      typedef typename Increment_dimension<typename R_::Max_ambient_dimension>::type D2;
+      typedef typename R_::LA::template Rebind_dimension<Dynamic_dimension_tag,D2>::Other LA;
+      typedef typename LA::Square_matrix Matrix;
+      typedef typename LA::Vector Vec;
+      typedef typename LA::Construct_vector CVec;
+      typename Get_functor<R_, Translated_point_tag>::type tp(this->kernel());
+      typename Get_functor<R_, Scaled_vector_tag>::type sv(this->kernel());
+      typename Get_functor<R_, Difference_of_points_tag>::type dp(this->kernel());
      typename Get_functor<R_, Scalar_product_tag>::type sp(this->kernel());
-      Matrix m(k,k);
-      Vec b(k);
-      Vec l(k);
-      int j,i=0;
-      for(Iter f2=f; f2!=e; ++f2,++i){
-        WPoint const& wp = *f2;
+      typename Get_functor<R_, Sum_of_vectors_tag>::type pv(this->kernel());
+      typename Get_functor<R_, Squared_length_tag>::type sl(this->kernel());
+
+      Matrix m(k-1,k-1);
+      Vec b = typename CVec::Dimension()(k-1);
+      std::vector<Vector> vecs; vecs.reserve(k-1);
+      for(int i=0; ++f!=e; ++i) {
+        WPoint const& wp = *f;
        Point const& p = pdw(wp);
-        b(i) = m(i,i) = sdo(p) - pw(wp);
-        j=0;
-        for(Iter f3=f; f3!=e; ++f3,++j){
-          // FIXME: scalar product of points ???
-          m(j,i) = m(i,j) = sp(p,pdw(*f3));
-        }
+        vecs.emplace_back(dp(p, p0));
+        b[i] = w0 - pw(wp);
      }
-      for(i=1;i<k;++i){
-        b(i)-=b(0);
-        for(j=0;j<k;++j){
-          m(i,j)=2*(m(i,j)-m(0,j));
+      // Only need to fill the lower half
+      for(int i = 0; i < k-1; ++i){
+        for(int j = i; j < k-1; ++j){
+          m(j, i) = sp(vecs[i], vecs[j]);
+#if ! EIGEN_VERSION_AT_LEAST(3, 3, 5)
+          m(i, j) = m(j, i);
+#endif
        }
+        b[i] += m(i, i);
+        b[i] /= 2;
      }
-      for(j=0;j<k;++j) m(0,j)=1;
-      b(0)=1;
-
-      LAd::solve(l,std::move(m),std::move(b));
-
-      typename LA::Vector center=typename LA::Construct_vector::Dimension()(d);
-      for(i=0;i<d;++i) center(i)=0;
-      j=0;
-      for(Iter f2=f;f2!=e;++f2,++j){
-        WPoint const& wp = *f2;
-        Point const& p = pdw(wp);
-        for(i=0;i<d;++i){
-          center(i)+=l(j)*c(p,i);
-        }
-      }
-
-      Point c = cp(d, LA::vector_begin(center), LA::vector_end(center));
-      FT r2 = pdp (wp0, c);
-      return cwp(std::move(c), std::move(r2));
+      // Assumes Eigen...
+#if EIGEN_VERSION_AT_LEAST(3, 3, 5)
+      Vec res = m.ldlt().solve(b);
+#else
+      // Older versions of Eigen use 1/highest as tolerance,
+      // which we have no way to set to 0 for exact types.
+      // Use something slow but that should work.
+      Vec res = m.fullPivLu().solve(b);
+#endif
+      Vector to_center = sv(vecs[0], res[0]);
+      for(int i=1;i<k-1;++i)
+        to_center = pv(to_center, sv(vecs[i],res[i]));
+      return cwp(tp(p0, to_center), sl(to_center) - w0);
    }
  }
 };
@ -279,7 +285,7 @@ CGAL_KD_DEFAULT_FUNCTOR(Power_side_of_power_sphere_tag,(CartesianDKernelFunctors
 CGAL_KD_DEFAULT_FUNCTOR(In_flat_power_side_of_power_sphere_tag,(CartesianDKernelFunctors::In_flat_power_side_of_power_sphere<K>),(Weighted_point_tag),(In_flat_power_side_of_power_sphere_raw_tag,Point_drop_weight_tag,Point_weight_tag));
 CGAL_KD_DEFAULT_FUNCTOR(Power_distance_tag,(CartesianDKernelFunctors::Power_distance<K>),(Weighted_point_tag,Point_tag),(Squared_distance_tag,Point_drop_weight_tag,Point_weight_tag));
 CGAL_KD_DEFAULT_FUNCTOR(Power_distance_to_point_tag,(CartesianDKernelFunctors::Power_distance_to_point<K>),(Weighted_point_tag,Point_tag),(Squared_distance_tag,Point_drop_weight_tag,Point_weight_tag));
-CGAL_KD_DEFAULT_FUNCTOR(Power_center_tag,(CartesianDKernelFunctors::Power_center<K>),(Weighted_point_tag,Point_tag),(Compute_point_cartesian_coordinate_tag,Construct_ttag<Point_tag>,Construct_ttag<Weighted_point_tag>,Point_dimension_tag,Squared_distance_to_origin_tag,Point_drop_weight_tag,Point_weight_tag,Power_distance_to_point_tag));
+CGAL_KD_DEFAULT_FUNCTOR(Power_center_tag,(CartesianDKernelFunctors::Power_center<K>),(Weighted_point_tag,Point_tag,Vector_tag),(Compute_point_cartesian_coordinate_tag,Construct_ttag<Point_tag>,Construct_ttag<Weighted_point_tag>,Point_dimension_tag,Squared_distance_to_origin_tag,Point_drop_weight_tag,Point_weight_tag,Power_distance_to_point_tag,Translated_point_tag,Scaled_vector_tag,Difference_of_points_tag,Scalar_product_tag,Sum_of_vectors_tag,Squared_length_tag));
 CGAL_KD_DEFAULT_FUNCTOR(Power_side_of_bounded_power_circumsphere_tag,(CartesianDKernelFunctors::Power_side_of_bounded_power_circumsphere<K>),(Weighted_point_tag),(Power_distance_tag,Power_center_tag));
 } // namespace CGAL
 #endif
--- a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h
+++ b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h
@ -27,6 +27,7 @@
 #include <CGAL/predicates/sign_of_determinant.h>
 #include <functional>
 #include <initializer_list>
+#include <vector>

 namespace CGAL {
 namespace CartesianDKernelFunctors {
@ -603,25 +604,57 @@ namespace CartesianDKernelFunctors {
 template <class R_> struct Construct_circumcenter : Store_kernel<R_> {
  CGAL_FUNCTOR_INIT_STORE(Construct_circumcenter)
  typedef typename Get_type<R_, Point_tag>::type Point;
+  typedef typename Get_type<R_, Vector_tag>::type Vector;
  typedef Point result_type;
  typedef typename Get_type<R_, FT_tag>::type FT;
  template <class Iter>
  result_type operator()(Iter f, Iter e)const{
-    typedef typename Get_type<R_, Point_tag>::type      Point;
-    typedef typename R_::LA LA;
-    typename Get_functor<R_, Compute_point_cartesian_coordinate_tag>::type c(this->kernel());
-    typename Get_functor<R_, Construct_ttag<Point_tag> >::type cp(this->kernel());
    typename Get_functor<R_, Point_dimension_tag>::type pd(this->kernel());
-    typename Get_functor<R_, Squared_distance_to_origin_tag>::type sdo(this->kernel());

    Point const& p0=*f;
    int d = pd(p0);
-    if (d+1 == std::distance(f,e))
+    int k = static_cast<int>(std::distance(f,e));
+    // Sorted from fastest to slowest, whether the dimension is 2, 3 or 4. It may be worth checking at some point.
+    CGAL_assume(k>=1);
+    if(k==1) return p0;
+    if(k==2){
+      typename Get_functor<R_, Midpoint_tag>::type mid(this->kernel());
+      return mid(p0, *++f);
+    }
+    // TODO: check for degenerate cases in all the following cases?
+    if(k==3){
+      // Same equations as in the general case, but solved by hand (Cramer)
+      // (c-r).(p-r)=(p-r)²/2
+      // (c-r).(q-r)=(q-r)²/2
+      typename Get_functor<R_, Squared_length_tag>::type sl(this->kernel());
+      typename Get_functor<R_, Scalar_product_tag>::type sp(this->kernel());
+      typename Get_functor<R_, Scaled_vector_tag>::type sv(this->kernel());
+      typename Get_functor<R_, Difference_of_points_tag>::type dp(this->kernel());
+      typename Get_functor<R_, Translated_point_tag>::type tp(this->kernel());
+      Iter f2=f;
+      Point const& q=*++f2;
+      Point const& r=*++f2;
+      Vector u = dp(p0, r);
+      Vector v = dp(q, r);
+      FT uv = sp(u, v);
+      FT u2 = sl(u);
+      FT v2 = sl(v);
+      FT den = 2 * (u2 * v2 - CGAL::square(uv));
+      FT a = (u2 - uv) * v2 / den;
+      FT b = (v2 - uv) * u2 / den;
+      // Wasteful if we only want the radius
+      return tp(tp(r, sv(u, a)), sv(v, b));
+    }
+    if (k == d+1)
    {
      // 2*(x-y).c == x^2-y^2
+      typedef typename R_::LA LA;
      typedef typename LA::Square_matrix Matrix;
      typedef typename LA::Vector Vec;
      typedef typename LA::Construct_vector CVec;
+      typename Get_functor<R_, Compute_point_cartesian_coordinate_tag>::type c(this->kernel());
+      typename Get_functor<R_, Construct_ttag<Point_tag> >::type cp(this->kernel());
+      typename Get_functor<R_, Squared_distance_to_origin_tag>::type sdo(this->kernel());
      FT const& n0 = sdo(p0);
      Matrix m(d,d);
      Vec b = typename CVec::Dimension()(d);
@ -635,72 +668,59 @@ template <class R_> struct Construct_circumcenter : Store_kernel<R_> {
        b[i] = sdo(p) - n0;
      }
      CGAL_assertion (i == d);
-      Vec res = typename CVec::Dimension()(d);;
-      //std::cout << "Mat: " << m << "\n Vec: " << one << std::endl;
-      LA::solve(res, std::move(m), std::move(b));
-      //std::cout << "Sol: " << res << std::endl;
+      //Vec res = typename CVec::Dimension()(d);;
+      //LA::solve(res, std::move(m), std::move(b));
+      // We already assume Eigen below...
+      Vec res=m.partialPivLu().solve(b);
      return cp(d,LA::vector_begin(res),LA::vector_end(res));
    }
-    else
+
    {
-      /*
-       * Matrix P=(p1, p2, ...) (each point as a column)
-       * Matrix Q=2*t(p2-p1,p3-p1, ...) (each vector as a line)
-       * Matrix M: QP, adding a line of 1 at the top
-       * Vector B: (1, p2^2-p1^2, p3^2-p1^2, ...)
-       * Solve ML=B, the center of the sphere is PL
-       *
-       * It would likely be faster to write P then transpose, multiply,
-       * etc instead of doing it by hand.
-       */
-      // TODO: check for degenerate cases?
-
-      typedef typename R_::Max_ambient_dimension D2;
-      typedef typename R_::LA::template Rebind_dimension<Dynamic_dimension_tag,D2>::Other LAd;
-      typedef typename LAd::Square_matrix Matrix;
-      typedef typename LAd::Vector Vec;
+      // The general case. ui=p(i+1)-p0, center-p0=c=sum ai*ui, c.ui=ui²/2, M*a=b with M symmetric
+      typedef typename Increment_dimension<typename R_::Max_ambient_dimension>::type D2;
+      typedef typename R_::LA::template Rebind_dimension<Dynamic_dimension_tag,D2>::Other LA;
+      typedef typename LA::Square_matrix Matrix;
+      typedef typename LA::Vector Vec;
+      typedef typename LA::Construct_vector CVec;
+      typename Get_functor<R_, Translated_point_tag>::type tp(this->kernel());
+      typename Get_functor<R_, Scaled_vector_tag>::type sv(this->kernel());
+      typename Get_functor<R_, Difference_of_points_tag>::type dp(this->kernel());
      typename Get_functor<R_, Scalar_product_tag>::type sp(this->kernel());
-      int k=static_cast<int>(std::distance(f,e));
-      Matrix m(k,k);
-      Vec b(k);
-      Vec l(k);
-      int j,i=0;
-      // We are doing a quadratic number of *f, which can be costly with transforming_iterator.
-      for(Iter f2=f;f2!=e;++f2,++i){
-        Point const& p2=*f2;
-        b(i)=m(i,i)=sdo(p2);
-        j=0;
-        for(Iter f3=f;f3!=e;++f3,++j){
-          m(j,i)=m(i,j)=sp(p2,*f3);
+      Matrix m(k-1,k-1);
+      Vec b = typename CVec::Dimension()(k-1);
+      std::vector<Vector> vecs; vecs.reserve(k-1);
+      while(++f!=e)
+        vecs.emplace_back(dp(*f,p0));
+      // Only need to fill the lower half
+      for(int i=0;i<k-1;++i){
+        for(int j=i;j<k-1;++j) {
+          m(j,i)=sp(vecs[i],vecs[j]);
+#if ! EIGEN_VERSION_AT_LEAST(3, 3, 5)
+          m(i,j)=m(j,i);
+#endif
        }
+        b[i]=m(i,i)/2;
      }
-      for(i=1;i<k;++i){
-        b(i)-=b(0);
-        for(j=0;j<k;++j){
-          m(i,j)=2*(m(i,j)-m(0,j));
-        }
-      }
-      for(j=0;j<k;++j) m(0,j)=1;
-      b(0)=1;
-
-      LAd::solve(l,std::move(m),std::move(b));
-
-      typename LA::Vector center=typename LA::Construct_vector::Dimension()(d);
-      for(i=0;i<d;++i) center(i)=0;
-      j=0;
-      for(Iter f2=f;f2!=e;++f2,++j){
-        for(i=0;i<d;++i){
-          center(i)+=l(j)*c(*f2,i);
-        }
-      }
-
-      return cp(LA::vector_begin(center),LA::vector_end(center));
+      // Assumes Eigen...
+#if EIGEN_VERSION_AT_LEAST(3, 3, 5)
+      Vec res=m.ldlt().solve(b);
+#else
+      // Older versions of Eigen use 1/highest as tolerance,
+      // which we have no way to set to 0 for exact types.
+      // Use something slow but that should work.
+      Vec res=m.fullPivLu().solve(b);
+#endif
+      Point center=p0;
+      // Wasteful if we only want the radius
+      for(int i=0;i<k-1;++i)
+        center=tp(center,sv(vecs[i],res[i]));
+      return center;
    }
  }
 };
 }

-CGAL_KD_DEFAULT_FUNCTOR(Construct_circumcenter_tag,(CartesianDKernelFunctors::Construct_circumcenter<K>),(Point_tag),(Construct_ttag<Point_tag>,Compute_point_cartesian_coordinate_tag,Scalar_product_tag,Squared_distance_to_origin_tag,Point_dimension_tag));
+CGAL_KD_DEFAULT_FUNCTOR(Construct_circumcenter_tag,(CartesianDKernelFunctors::Construct_circumcenter<K>),(Point_tag,Vector_tag),(Construct_ttag<Point_tag>,Compute_point_cartesian_coordinate_tag,Scalar_product_tag,Squared_distance_to_origin_tag,Point_dimension_tag,Translated_point_tag,Scaled_vector_tag,Difference_of_points_tag,Squared_length_tag));

 namespace CartesianDKernelFunctors {
 template <class R_> struct Squared_circumradius : Store_kernel<R_> {
--- a/NewKernel_d/test/NewKernel_d/Epick_d.cpp
+++ b/NewKernel_d/test/NewKernel_d/Epick_d.cpp
@ -712,6 +712,38 @@ void test3(){
  std::ostringstream sv1; sv1 << v1; assert(sv1.str()=="3 3 2 1");
  std::istringstream sv2("3 4 5 6"); sv2 >> v1; assert(v1[0]==4&&v1[1]==5);
 }
+template<class Ker>
+void test4(){
+  typedef typename Ker::Point_d P;
+  typedef typename Ker::Weighted_point_d WP;
+  typedef typename Ker::Construct_circumcenter_d CCc;
+  typedef typename Ker::Equal_d E;
+  typedef typename Ker::Power_center_d PC;
+  typedef typename Ker::Power_distance_d PoD;
+  typedef typename Ker::Affine_rank_d AR;
+  Ker k(4);
+  CCc ccc Kinit(construct_circumcenter_d_object);
+  E ed Kinit(equal_d_object);
+  PC pc Kinit(power_center_d_object);
+  PoD pod Kinit(power_distance_d_object);
+  AR ar Kinit(affine_rank_d_object);
+  auto mkpt=[](auto...x){double l[]{(double)x...};return P(std::begin(l), std::end(l));};
+  P tab1[]={mkpt(15,20,40,80),mkpt(10,23,36,80),mkpt(10,20,40,85),mkpt(10,15,40,80),mkpt(13,20,40,76)};
+  assert(ed(ccc(tab1+0, tab1+5),mkpt(10,20,40,80)));
+  P tab2[]={mkpt(15,20,40,80),mkpt(13,24,40,80),mkpt(10,25,40,80),mkpt(10,20,43,84)};
+  assert(ed(ccc(tab2+0, tab2+4),mkpt(10,20,40,80)));
+  P tab3[]={mkpt(15,20,35,80),mkpt(10,25,40,75),mkpt(13,24,37,76)};
+  assert(ed(ccc(tab3+0, tab3+3),mkpt(10,20,40,80)));
+  auto mkwpt=[](auto...x){double l[]{(double)x...};auto last=std::prev(std::end(l));return WP(P(std::begin(l), last),*last);};
+  WP tab4[]={mkwpt(89,17,29,97,14),mkwpt(86,99,64,26,44),mkwpt(40,9,13,91,20),mkwpt(41,30,93,13,10),mkwpt(45,6,98,9,0),mkwpt(0,0,0,0,0)};
+  for(int i=5;i>=1;--i){
+    tab4[i]=pc(tab4+0, tab4+i);
+    for(int j=0;j<i;++j)
+      assert(pod(tab4[i],tab4[j])==0);
+    auto drop=[](WP const&x){return x.point();};
+    assert(ar(CGAL::make_transforming_iterator(tab4+0,drop), CGAL::make_transforming_iterator(tab4+i+1,drop))==i-1);
+  }
+}
 template struct CGAL::Epick_d<CGAL::Dimension_tag<2> >;
 template struct CGAL::Epick_d<CGAL::Dimension_tag<3> >;
 template struct CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
@ -733,6 +765,7 @@ int main(){
  test2<CGAL::Epeck_d<CGAL::Dimension_tag<2>>>();
  test3<CGAL::Epeck_d<CGAL::Dimension_tag<3>>>();
  test3<CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>>();
+  test4<CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>>();
 #endif
 }

--- a/Number_types/include/CGAL/Interval_nt.h
+++ b/Number_types/include/CGAL/Interval_nt.h
@ -1573,6 +1573,8 @@ namespace Eigen {

    static inline Real epsilon() { return 0; }
    static inline Real dummy_precision() { return 0; }
+    static inline Real highest() { return Real((std::numeric_limits<double>::max)(), std::numeric_limits<double>::infinity()); }
+    static inline Real lowest() { return Real(-std::numeric_limits<double>::infinity(), std::numeric_limits<double>::lowest()); }

    // Costs could depend on b.
    enum {
--- a/Point_set_processing_3/include/CGAL/edge_aware_upsample_point_set.h
+++ b/Point_set_processing_3/include/CGAL/edge_aware_upsample_point_set.h
@ -393,7 +393,6 @@ edge_aware_upsample_point_set(
                                       &&sharpness_angle <= 90);
  CGAL_point_set_processing_precondition(edge_sensitivity >= 0
                                       &&edge_sensitivity <= 1);
-  CGAL_point_set_processing_precondition(neighbor_radius > 0);

  edge_sensitivity *= 10;  // just project [0, 1] to [0, 10].

--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/repair_polygon_soup_example.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/repair_polygon_soup_example.cpp
@ -5,30 +5,55 @@
 #include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
 #include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>

+#include <algorithm>
+#include <array>
 #include <iostream>
 #include <vector>

 typedef CGAL::Exact_predicates_inexact_constructions_kernel     K;
+typedef K::FT                                                   FT;
 typedef K::Point_3                                              Point_3;

-typedef std::vector<std::size_t>                                CGAL_Polygon;
 typedef CGAL::Surface_mesh<Point_3>                             Mesh;

+typedef std::array<FT, 3>                                       Custom_point;
+typedef std::vector<std::size_t>                                CGAL_Polygon;
+
 namespace PMP = CGAL::Polygon_mesh_processing;

+struct Array_traits
+{
+  struct Equal_3
+  {
+    bool operator()(const Custom_point& p, const Custom_point& q) const {
+      return (p == q);
+    }
+  };
+
+  struct Less_xyz_3
+  {
+    bool operator()(const Custom_point& p, const Custom_point& q) const {
+      return std::lexicographical_compare(p.begin(), p.end(), q.begin(), q.end());
+    }
+  };
+
+  Equal_3 equal_3_object() const { return Equal_3(); }
+  Less_xyz_3 less_xyz_3_object() const { return Less_xyz_3(); }
+};
+
 int main(int, char**)
 {
  // First, construct a polygon soup with some problems
-  std::vector<Point_3> points;
+  std::vector<std::array<FT, 3> > points;
  std::vector<CGAL_Polygon> polygons;

-  points.push_back(Point_3(0,0,0));
-  points.push_back(Point_3(1,0,0));
-  points.push_back(Point_3(0,1,0));
-  points.push_back(Point_3(-1,0,0));
-  points.push_back(Point_3(0,-1,0));
-  points.push_back(Point_3(0,1,0)); // duplicate point
-  points.push_back(Point_3(0,-2,0)); // unused point
+  points.push_back(CGAL::make_array<FT>(0,0,0));
+  points.push_back(CGAL::make_array<FT>(1,0,0));
+  points.push_back(CGAL::make_array<FT>(0,1,0));
+  points.push_back(CGAL::make_array<FT>(-1,0,0));
+  points.push_back(CGAL::make_array<FT>(0,-1,0));
+  points.push_back(CGAL::make_array<FT>(0,1,0)); // duplicate point
+  points.push_back(CGAL::make_array<FT>(0,-2,0)); // unused point

  CGAL_Polygon p;
  p.push_back(0); p.push_back(1); p.push_back(2);
@ -67,10 +92,12 @@ int main(int, char**)
  p.push_back(4); p.push_back(3); p.push_back(2); p.push_back(1);
  polygons.push_back(p);

-  PMP::repair_polygon_soup(points, polygons);
-  PMP::orient_polygon_soup(points, polygons);
+  std::cout << "Before reparation, the soup has " << points.size() << " vertices and " << polygons.size() << " faces" << std::endl;
+  PMP::repair_polygon_soup(points, polygons, CGAL::parameters::geom_traits(Array_traits()));
+  std::cout << "After reparation, the soup has " << points.size() << " vertices and " << polygons.size() << " faces" << std::endl;

  Mesh mesh;
+  PMP::orient_polygon_soup(points, polygons);
  PMP::polygon_soup_to_polygon_mesh(points, polygons, mesh);

  std::cout << "Mesh has " << num_vertices(mesh) << " vertices and " << num_faces(mesh) << " faces" << std::endl;
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h
@ -23,6 +23,7 @@
 #include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
 #include <CGAL/Polygon_mesh_processing/border.h>
 #include <CGAL/Polygon_mesh_processing/repair.h>
+#include <CGAL/Polygon_mesh_processing/internal/Corefinement/Generic_clip_output_builder.h>
 #include <CGAL/iterator.h>

 #include <CGAL/AABB_triangle_primitive.h>
@ -400,6 +401,97 @@ void split_along_edges(TriangleMesh& tm,
  CGAL_assertion(is_valid_polygon_mesh(tm));
 }

+template <class TriangleMesh,
+          class NamedParameters1,
+          class NamedParameters2>
+void
+generic_clip_impl(
+        TriangleMesh& tm1,
+        TriangleMesh& tm2,
+  const NamedParameters1& np1,
+  const NamedParameters2& np2)
+{
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+// Vertex point maps
+  //for input meshes
+  typedef typename GetVertexPointMap<TriangleMesh,
+                                     NamedParameters1>::type Vpm;
+  typedef typename GetVertexPointMap<TriangleMesh,
+                                     NamedParameters2>::type Vpm2;
+  CGAL_USE_TYPE(Vpm2);
+  CGAL_assertion_code(
+    static const bool same_vpm = (boost::is_same<Vpm,Vpm2>::value); )
+  CGAL_static_assertion(same_vpm);
+
+  Vpm vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point),
+                              get_property_map(boost::vertex_point, tm1));
+
+  Vpm vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point),
+                              get_property_map(boost::vertex_point, tm2));
+
+  if (&tm1==&tm2)
+  {
+    // TODO mark all edges
+    return;
+  }
+
+  // handle case of empty meshes (isolated vertices are ignored)
+  if (faces(tm1).empty())
+    return;
+
+// Edge is-constrained maps
+  //for input meshes
+  typedef typename internal_np::Lookup_named_param_def <
+    internal_np::edge_is_constrained_t,
+    NamedParameters1,
+    Corefinement::No_mark<TriangleMesh>//default
+  > ::type User_ecm1;
+
+  // User and internal edge is-constrained map
+  typedef typename boost::template property_map<TriangleMesh, CGAL::dynamic_edge_property_t<bool> >::type Algo_ecm1;
+  typedef Corefinement::No_mark<TriangleMesh> Ecm2;
+  typedef OR_property_map<Algo_ecm1, User_ecm1> Ecm1;
+  typedef Corefinement::Ecm_bind<TriangleMesh, Ecm1, Ecm2> Ecm_in;
+
+  Algo_ecm1 algo_ecm1  = get(CGAL::dynamic_edge_property_t<bool>(), tm1);
+  Ecm1 ecm1 = Ecm1(algo_ecm1, choose_parameter<User_ecm1>(get_parameter(np1, internal_np::edge_is_constrained)));
+  Ecm2 ecm2;
+
+  // Face index point maps
+  typedef typename CGAL::GetInitializedFaceIndexMap<TriangleMesh, NamedParameters1>::type FaceIndexMap1;
+  FaceIndexMap1 fid_map1 = get_initialized_face_index_map(tm1, np1);
+
+  const bool use_compact_clipper =
+    choose_parameter(get_parameter(np1, internal_np::use_compact_clipper), true);
+
+
+  // User visitor
+  typedef typename internal_np::Lookup_named_param_def <
+    internal_np::graph_visitor_t,
+    NamedParameters1,
+    Corefinement::Default_visitor<TriangleMesh>//default
+  > ::type User_visitor;
+  User_visitor uv(choose_parameter<User_visitor>(get_parameter(np1, internal_np::graph_visitor)));
+
+  // surface intersection algorithm call
+  typedef Corefinement::Generic_clip_output_builder<TriangleMesh,
+                                                    Vpm,
+                                                    Algo_ecm1,
+                                                    FaceIndexMap1,
+                                                    Default> Ob;
+
+  typedef Corefinement::Surface_intersection_visitor_for_corefinement<
+    TriangleMesh, Vpm, Ob, Ecm_in, User_visitor> Algo_visitor;
+  Ecm_in ecm_in(tm1,tm2,ecm1,ecm2);
+  Ob ob(tm1, tm2, vpm1, vpm2, algo_ecm1, fid_map1, use_compact_clipper);
+
+  Corefinement::Intersection_of_triangle_meshes<TriangleMesh, Vpm, Algo_visitor >
+    functor(tm1, tm2, vpm1, vpm2, Algo_visitor(uv,ob,ecm_in,&tm2));
+  functor(CGAL::Emptyset_iterator(), false, true);
+}
+
 } // end of internal namespace

 /**
@ -473,6 +565,14 @@ void split_along_edges(TriangleMesh& tm,
  *     \cgalParamExtra{This option has an effect only if a surface and not a volume is clipped,
  *                     (i.e., if `clip_volume` is `false` or if `tm` is open).}
  *   \cgalParamNEnd
+  *   \cgalParamNBegin{do_not_modify}
+  *     \cgalParamDescription{(`np_c` only) if `true`, `clipper` will not be modified.}
+  *     \cgalParamType{Boolean}
+  *     \cgalParamDefault{`false`}
+  *     \cgalParamExtra{If this option is set to `true`, `tm` is no longer required to be without self-intersection.
+  *                     Setting this option to `true` will automatically set `throw_on_self_intersection` to `false`
+  *                     and `clip_volume` to `false`.}
+  *   \cgalParamNEnd
  * \cgalNamedParamsEnd
  *
  * @return `true` if the output surface mesh is manifold.
@ -487,6 +587,14 @@ clip(TriangleMesh& tm,
     const NamedParameters1& np_tm,
     const NamedParameters2& np_c)
 {
+  if (parameters::choose_parameter(parameters::get_parameter(np_c, internal_np::do_not_modify), false))
+  {
+    CGAL_assertion(is_closed(clipper));
+
+    internal::generic_clip_impl(tm, clipper, np_tm, np_c);
+    return true;
+  }
+
  const bool clip_volume =
    parameters::choose_parameter(parameters::get_parameter(np_tm, internal_np::clip_volume), false);

@ -710,6 +818,13 @@ bool clip(TriangleMesh& tm,
  *                           will be thrown if at least one self-intersection is found.}
  *     \cgalParamType{Boolean}
  *     \cgalParamDefault{`false`}
+  *   \cgalParamNBegin{do_not_modify}
+  *     \cgalParamDescription{(`np_s` only) if `true`, `splitter` will not be modified.}
+  *     \cgalParamType{Boolean}
+  *     \cgalParamDefault{`false`}
+  *     \cgalParamExtra{If this option is set to `true`, `tm` is no longer required to be without self-intersection.
+  *                     Setting this option to `true` will automatically set `throw_on_self_intersection` to `false`
+  *                     and `clip_volume` to `false`.}
  *   \cgalParamNEnd
  *
  * \cgalNamedParamsEnd
@ -742,9 +857,11 @@ void split(TriangleMesh& tm,
  // create a constrained edge map and corefine input mesh with the splitter,
  // and mark edges

+  const bool do_not_modify_splitter = choose_parameter(get_parameter(np_s, internal_np::do_not_modify), false);
+
  PMP::corefine(tm, splitter,
                CGAL::parameters::vertex_point_map(vpm_tm).edge_is_constrained_map(ecm),
-                CGAL::parameters::vertex_point_map(vpm_s));
+                CGAL::parameters::vertex_point_map(vpm_s).do_not_modify(do_not_modify_splitter));

  //split mesh along marked edges
  internal::split_along_edges(tm, ecm, vpm_tm);
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/corefinement.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/corefinement.h
@ -695,6 +695,14 @@ corefine_and_compute_difference(      TriangleMesh& tm1,
 *     \cgalParamDefault{`false`}
 *     \cgalParamExtra{`np1` only}
 *   \cgalParamNEnd
+ *   \cgalParamNBegin{do_not_modify}
+ *     \cgalParamDescription{if `true`, the corresponding mesh will not be updated.}
+ *     \cgalParamType{Boolean}
+ *     \cgalParamDefault{`false`}
+ *     \cgalParamExtra{If this parameter is set to `true` for both meshes nothing will be done.
+ *                      If this option is set to `true` for one mesh,
+ *                      the other mesh is no longer required to be without self-intersection.}
+ *   \cgalParamNEnd
 * \cgalNamedParamsEnd
 *
 */
@ -710,6 +718,19 @@ corefine(      TriangleMesh& tm1,
  using parameters::choose_parameter;
  using parameters::get_parameter;

+  TriangleMesh* const_mesh_ptr=nullptr;
+  if (choose_parameter(get_parameter(np1, internal_np::do_not_modify), false))
+  {
+    if (choose_parameter(get_parameter(np2, internal_np::do_not_modify), false))
+      return;
+    const_mesh_ptr=&tm1;
+  }
+  else
+  {
+    if (choose_parameter(get_parameter(np2, internal_np::do_not_modify), false))
+      const_mesh_ptr=&tm2;
+  }
+
  const bool throw_on_self_intersection =
    choose_parameter(get_parameter(np1, internal_np::throw_on_self_intersection), false);

@ -769,7 +790,7 @@ corefine(      TriangleMesh& tm1,
  Ob ob;
  Ecm ecm(tm1,tm2,ecm1,ecm2);
  Corefinement::Intersection_of_triangle_meshes<TriangleMesh, Vpm, Algo_visitor>
-    functor(tm1, tm2, vpm1, vpm2, Algo_visitor(uv,ob,ecm));
+    functor(tm1, tm2, vpm1, vpm2, Algo_visitor(uv,ob,ecm,const_mesh_ptr));
  functor(CGAL::Emptyset_iterator(), throw_on_self_intersection, true);
 }

--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Generic_clip_output_builder.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Generic_clip_output_builder.h
@ -0,0 +1,243 @@
+// Copyright (c) 2020 GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Sebastien Loriot
+
+#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_GENERIC_CLIP_OUTPUT_BUILDER_H
+#define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_GENERIC_CLIP_OUTPUT_BUILDER_H
+
+#include <CGAL/license/Polygon_mesh_processing/corefinement.h>
+
+
+#include <CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h>
+
+#include <CGAL/Polygon_mesh_processing/orientation.h>
+#include <CGAL/Polygon_mesh_processing/connected_components.h>
+#include <CGAL/Side_of_triangle_mesh.h>
+
+#include <CGAL/property_map.h>
+#include <CGAL/Default.h>
+#include <CGAL/use.h>
+
+#include <boost/dynamic_bitset.hpp>
+
+
+namespace CGAL {
+namespace Polygon_mesh_processing {
+namespace Corefinement {
+
+
+namespace PMP=Polygon_mesh_processing;
+namespace params=PMP::parameters;
+
+template <class TriangleMesh,
+          class VertexPointMap,
+          class Ecm1,
+          class FaceIdMap1,
+          class Kernel_=Default>
+class Generic_clip_output_builder
+{
+//Default typedefs
+  typedef typename Default::Get<
+    Kernel_,
+    typename Kernel_traits<
+      typename boost::property_traits<VertexPointMap>::value_type
+    >::Kernel >::type                                           Kernel;
+
+// graph_traits typedefs
+  typedef TriangleMesh                                              TM;
+  typedef boost::graph_traits<TM>                                   GT;
+  typedef typename GT::edge_descriptor                 edge_descriptor;
+  typedef typename GT::face_descriptor                 face_descriptor;
+  typedef typename GT::halfedge_descriptor         halfedge_descriptor;
+  typedef typename GT::vertex_descriptor             vertex_descriptor;
+// Internal typedefs
+  typedef std::size_t                                          Node_id;
+  typedef std::pair<Node_id,Node_id>                      Node_id_pair;
+  // to maintain a halfedge on each polyline per TriangleMesh + pair<bool,size_t>
+  // with first = "is the key (pair<Node_id,Node_id>) was reversed?" and
+  // second is the number of edges -1 in the polyline
+  typedef std::map< Node_id_pair,
+                    std::pair< std::map<TriangleMesh*,
+                                        halfedge_descriptor>,
+                               std::pair<bool,std::size_t> > >
+                                              An_edge_per_polyline_map;
+
+  typedef boost::unordered_map<vertex_descriptor, Node_id> Node_id_map;
+  typedef boost::unordered_map<edge_descriptor,
+                               edge_descriptor>               Edge_map;
+//Data members
+  TriangleMesh &tm1, &tm2;
+  // property maps of input meshes
+  const VertexPointMap vpm1;
+  const VertexPointMap vpm2;
+  Ecm1 ecm1;
+  FaceIdMap1 fids1;
+  bool use_compact_clipper;
+
+  // mapping vertex to node id
+  Node_id_map vertex_to_node_id1;
+
+  // orientation of input surface meshes
+  bool is_tm2_inside_out;
+  // constants
+  const Node_id NID;
+
+  typename An_edge_per_polyline_map::iterator last_polyline;
+
+  Node_id get_node_id(vertex_descriptor v,
+                      const Node_id_map& node_ids)
+  {
+    typename Node_id_map::const_iterator it = node_ids.find(v);
+    if (it == node_ids.end())
+      return NID;
+    return it->second;
+  }
+
+public:
+
+  Generic_clip_output_builder(TriangleMesh& tm1,
+                              TriangleMesh& tm2,
+                              const VertexPointMap vpm1,
+                              const VertexPointMap vpm2,
+                              const Ecm1& ecm1,
+                              FaceIdMap1 fids1,
+                              bool use_compact_clipper)
+    : tm1(tm1), tm2(tm2)
+    , vpm1(vpm1), vpm2(vpm2)
+    , ecm1(ecm1)
+    , fids1(fids1)
+    , use_compact_clipper(use_compact_clipper)
+    , is_tm2_inside_out( !PMP::is_outward_oriented(tm2, parameters::vertex_point_map(vpm2)) )
+    , NID((std::numeric_limits<Node_id>::max)())
+  {}
+
+// functions called by the intersection visitor
+  void start_new_polyline(Node_id, Node_id) {}
+  void add_node_to_polyline(Node_id) {}
+  void set_edge_per_polyline(TriangleMesh&, Node_id_pair, halfedge_descriptor){}
+
+  void set_vertex_id(vertex_descriptor v, Node_id node_id, const TriangleMesh& tm)
+  {
+    CGAL_USE(tm);
+    CGAL_assertion(&tm == &tm1);
+    vertex_to_node_id1.insert( std::make_pair(v, node_id) );
+  }
+
+  template <class Nodes_vector, class Mesh_to_map_node>
+  void operator()(
+    const Nodes_vector& nodes,
+    bool /* input_have_coplanar_faces */,
+    const boost::dynamic_bitset<>& /* is_node_of_degree_one */,
+    const Mesh_to_map_node&)
+  {
+    // The property map must be either writable or well-initialized
+    if( CGAL::internal::Is_writable_property_map<FaceIdMap1>::value &&
+        !BGL::internal::is_index_map_valid(fids1, num_faces(tm1), faces(tm1)) )
+    {
+      BGL::internal::initialize_face_index_map(fids1, tm1);
+    }
+    CGAL_assertion(BGL::internal::is_index_map_valid(fids1, num_faces(tm1), faces(tm1)));
+
+    // (1) Assign a patch id to each face indicating in which connected
+    // component limited by intersection edges of the surface they are.
+    std::vector<std::size_t> tm1_patch_ids( num_faces(tm1),NID );
+
+    std::size_t nb_patches_tm1 =
+      PMP::connected_components(tm1,
+                                bind_property_maps(fids1,make_property_map(&tm1_patch_ids[0])),
+                                params::edge_is_constrained_map(ecm1)
+                                       .face_index_map(fids1));
+
+    std::vector <std::size_t> tm1_patch_sizes(nb_patches_tm1, 0);
+    for(std::size_t i : tm1_patch_ids)
+      if(i!=NID)
+        ++tm1_patch_sizes[i];
+
+    // Use the class Side_of_triangle_mesh to classify each patch
+    boost::dynamic_bitset<> patch_status_not_set_tm1(nb_patches_tm1);
+    patch_status_not_set_tm1.set();
+
+    typedef Side_of_triangle_mesh<TriangleMesh,
+                                  Kernel,
+                                  VertexPointMap> Inside_poly_test;
+
+    CGAL::Bounded_side in_tm2 = is_tm2_inside_out
+                              ? ON_UNBOUNDED_SIDE : ON_BOUNDED_SIDE;
+
+    Inside_poly_test inside_tm2(tm2, vpm2);
+    std::size_t nb_classified=0;
+    std::vector<std::size_t> cc_to_keep;
+    for(face_descriptor f : faces(tm1))
+    {
+      const std::size_t f_id = get(fids1, f);
+      const std::size_t patch_id = tm1_patch_ids[ f_id ];
+      if ( patch_status_not_set_tm1.test( patch_id ) )
+      {
+        patch_status_not_set_tm1.reset( patch_id );
+        halfedge_descriptor h = halfedge(f, tm1);
+        Node_id index_p1 = get_node_id(target(h, tm1), vertex_to_node_id1);
+        std::array<Node_id, 3> fnids = { index_p1, index_p1, index_p1 };
+        if (index_p1 != NID)
+        {
+          h=next(h, tm1);
+          index_p1 = get_node_id(target(h, tm1), vertex_to_node_id1);
+          fnids[1]=index_p1;
+          if (index_p1 != NID)
+          {
+            h=next(h, tm1);
+            index_p1 = get_node_id(target(h, tm1), vertex_to_node_id1);
+            fnids[2]=index_p1;
+          }
+        }
+        if (index_p1 != NID)
+        {
+          typename Nodes_vector::Exact_kernel ek;
+          typedef typename Nodes_vector::Exact_kernel::Point_3 Exact_point_3;
+          Exact_point_3 e_centroid = centroid(nodes.exact_node(fnids[0]),
+                                              nodes.exact_node(fnids[1]),
+                                              nodes.exact_node(fnids[2]));
+
+          Bounded_side position = inside_tm2(e_centroid, ek);
+          if ( position==ON_BOUNDARY )
+          {
+            if (use_compact_clipper)
+            {
+              cc_to_keep.push_back(patch_id);
+            }
+          }
+          else
+            if ( position == in_tm2 )
+            {
+              cc_to_keep.push_back(patch_id);
+            }
+        }
+        else
+        {
+          Bounded_side position = inside_tm2( get(vpm1, target(h, tm1)));
+          CGAL_assertion( position != ON_BOUNDARY);
+          if ( position == in_tm2 )
+          {
+            cc_to_keep.push_back(patch_id);
+          }
+        }
+        if ( ++nb_classified==nb_patches_tm1) break;
+      }
+    }
+
+    PMP::keep_connected_components(tm1, cc_to_keep, bind_property_maps(fids1,make_property_map(&tm1_patch_ids[0])));
+  }
+};
+
+
+} } } // CGAL::Polygon_mesh_processing::Corefinement
+
+#undef CGAL_COREF_FUNCTION_CALL
+#endif // CGAL_POLYGON_MESH_PROCESSING_INTERNAL_GENERIC_CLIP_OUTPUT_BUILDER_H
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h
@ -170,6 +170,7 @@ private:
  OutputBuilder& output_builder;
  EdgeMarkMapBind marks_on_edges;
  bool input_with_coplanar_faces;
+  TriangleMesh* const_mesh_ptr;

  template <class Ecm1, class Ecm2>
  void call_put(Ecm_bind<TriangleMesh, Ecm1, Ecm2>& ecm,
@ -200,12 +201,13 @@ private:
 // visitor public functions
 public:
  Surface_intersection_visitor_for_corefinement(
-    UserVisitor& uv, OutputBuilder& o, const EdgeMarkMapBind& emm)
+    UserVisitor& uv, OutputBuilder& o, const EdgeMarkMapBind& emm, TriangleMesh* const_mesh_ptr=nullptr)
    : number_coplanar_vertices(0)
    , user_visitor(uv)
    , output_builder(o)
    , marks_on_edges(emm)
    , input_with_coplanar_faces(false)
+    , const_mesh_ptr(const_mesh_ptr)
  {}

  template<class Graph_node>
@ -332,34 +334,39 @@ public:

    //forward to the visitor
 //    user_visitor.new_node_added(node_id, type, h_1, h_2, is_target_coplanar, is_source_coplanar); // NODE_VISITOR_TAG
-    switch(type)
+    if (tm2_ptr!=const_mesh_ptr)
    {
-      case ON_FACE: //Face intersected by an edge
-        on_face[tm2_ptr][face(h_2,tm2)].push_back(node_id);
-      break;
-      case ON_EDGE: //Edge intersected by an edge
+      switch(type)
      {
-        on_edge[tm2_ptr][edge(h_2,tm2)].push_back(node_id);
-      //   check_node_on_non_manifold_edge(node_id,h_2,tm2);
+        case ON_FACE: //Face intersected by an edge
+          on_face[tm2_ptr][face(h_2,tm2)].push_back(node_id);
+        break;
+        case ON_EDGE: //Edge intersected by an edge
+        {
+          on_edge[tm2_ptr][edge(h_2,tm2)].push_back(node_id);
+        //   check_node_on_non_manifold_edge(node_id,h_2,tm2);
+        }
+        break;
+        case ON_VERTEX:
+        {
+          //grab original vertex that is on commom intersection
+          mesh_to_vertices_on_inter[tm2_ptr].insert(std::make_pair(node_id,h_2));
+          Node_id_to_vertex& node_id_to_vertex=mesh_to_node_id_to_vertex[tm2_ptr];
+          if (node_id_to_vertex.size()<=node_id)
+            node_id_to_vertex.resize(node_id+1,Graph_traits::null_vertex());
+          node_id_to_vertex[node_id]=target(h_2,tm2);
+          all_incident_faces_got_a_node_as_vertex(h_2,node_id,*tm2_ptr);
+        //   check_node_on_non_manifold_vertex(node_id,h_2,tm2);
+          output_builder.set_vertex_id(target(h_2, tm2), node_id, tm2);
+        }
+        break;
+        default:
+        return;
      }
-      break;
-      case ON_VERTEX:
-      {
-        //grab original vertex that is on commom intersection
-        mesh_to_vertices_on_inter[tm2_ptr].insert(std::make_pair(node_id,h_2));
-        Node_id_to_vertex& node_id_to_vertex=mesh_to_node_id_to_vertex[tm2_ptr];
-        if (node_id_to_vertex.size()<=node_id)
-          node_id_to_vertex.resize(node_id+1,Graph_traits::null_vertex());
-        node_id_to_vertex[node_id]=target(h_2,tm2);
-        all_incident_faces_got_a_node_as_vertex(h_2,node_id,*tm2_ptr);
-      //   check_node_on_non_manifold_vertex(node_id,h_2,tm2);
-        output_builder.set_vertex_id(target(h_2, tm2), node_id, tm2);
-      }
-      break;
-      default:
-      return;
    }

+    if (tm1_ptr==const_mesh_ptr) return;
+
    CGAL_assertion(!is_target_coplanar || !is_source_coplanar); //coplanar edge are not forwarded

    if ( is_target_coplanar )
@ -624,6 +631,7 @@ public:
          ++it)
    {
      TriangleMesh& tm=*it->first;
+      CGAL_assertion(&tm!=const_mesh_ptr);
    //   Face_boundaries& face_boundaries=mesh_to_face_boundaries[&tm];

      Node_to_target_of_hedge_map& nodes_to_hedge=it->second;
@ -693,6 +701,7 @@ public:
      it=on_edge.begin(); it!=on_edge.end(); ++it)
    {
      TriangleMesh& tm=*it->first;
+      CGAL_assertion(&tm!=const_mesh_ptr);
      const VertexPointMap& vpm=vpms[&tm];
      On_edge_map& on_edge_map=it->second;
      On_face_map& on_face_map=on_face[&tm];
@ -788,6 +797,7 @@ public:
      it=on_face.begin(); it!=on_face.end(); ++it)
    {
      TriangleMesh& tm=*it->first;
+      CGAL_assertion(&tm!=const_mesh_ptr);
      const VertexPointMap& vpm=vpms[&tm];
      On_face_map& on_face_map=it->second;
      Face_boundaries& face_boundaries=mesh_to_face_boundaries[&tm];
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Side_of_triangle_mesh/Point_inside_vertical_ray_cast.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Side_of_triangle_mesh/Point_inside_vertical_ray_cast.h
@ -75,13 +75,6 @@ public:
  {
    typename Traits::Bounding_box bbox = m_helper.get_tree_bbox(tree);

-    if(   point.x() < bbox.xmin() || point.x() > bbox.xmax()
-      || point.y() < bbox.ymin() || point.y() > bbox.ymax()
-      || point.z() < bbox.zmin() || point.z() > bbox.zmax() )
-    {
-      return ON_UNBOUNDED_SIDE;
-    }
-
    //the direction of the vertical ray depends on the position of the point in the bbox
    //in order to limit the expected number of nodes visited.
    Ray query = ray_functor(point, vector_functor(0,0,(2*point.z() <  bbox.zmax()+bbox.zmin()?-1:1)));
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Side_of_triangle_mesh/Ray_3_Triangle_3_traversal_traits.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Side_of_triangle_mesh/Ray_3_Triangle_3_traversal_traits.h
@ -216,6 +216,60 @@ public:
  }
 };

+//special case when ray query is from another Kernel K1 is the kernel compatible with the AABB-tree
+template<typename AABBTraits, class K1, class K2, class Helper>
+class K2_Ray_3_K1_Triangle_3_traversal_traits
+{
+  //the status indicates whether the query point is strictly inside the polyhedron, and the number of intersected triangles if yes
+  std::pair<boost::logic::tribool,std::size_t>& m_status;
+  bool m_stop;
+  const AABBTraits& m_aabb_traits;
+  typedef typename AABBTraits::Primitive Primitive;
+  typedef CGAL::AABB_node<AABBTraits> Node;
+  Helper m_helper;
+  CGAL::Cartesian_converter<K1,K2> to_K2;
+
+public:
+  K2_Ray_3_K1_Triangle_3_traversal_traits(std::pair<boost::logic::tribool,std::size_t>& status,
+                                          const AABBTraits& aabb_traits,
+                                          const Helper& h)
+    :m_status(status), m_stop(false), m_aabb_traits(aabb_traits), m_helper(h)
+  {m_status.first=true;}
+
+  bool go_further() const { return !m_stop; }
+
+  template<class Query>
+  void intersection(const Query& query, const Primitive& primitive)
+  {
+    Intersections::internal::r3t3_do_intersect_endpoint_position_visitor visitor;
+    std::pair<bool,Intersections::internal::R3T3_intersection::type> res=
+      Intersections::internal::do_intersect(to_K2(m_helper.get_primitive_datum(primitive, m_aabb_traits)),
+                                            query, K2(), visitor);
+
+    if (res.first){
+      switch (res.second){
+        case Intersections::internal::R3T3_intersection::CROSS_FACET:
+          ++m_status.second;
+        break;
+        case Intersections::internal::R3T3_intersection::ENDPOINT_IN_TRIANGLE:
+          m_status.first=false;
+          m_stop=true;
+        break;
+        default:
+          m_status.first=boost::logic::indeterminate;
+          m_stop=true;
+      }
+    }
+  }
+
+  template<class Query>
+  bool do_intersect(const Query& query, const Node& node) const
+  {
+    return CGAL::do_intersect(query, m_helper.get_node_bbox(node));
+  }
+};
+
+
 }// namespace internal
 }// namespace CGAL

--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h
@ -29,6 +29,7 @@
 #include <boost/range/value_type.hpp>
 #include <boost/range/reference.hpp>

+#include <array>
 #include <set>
 #include <type_traits>
 #include <vector>
@ -37,13 +38,26 @@ namespace CGAL {
 namespace Polygon_mesh_processing {
 namespace internal {

+template <typename PM_Point, typename PS_Point>
+PM_Point convert_to_pm_point(const PS_Point& p)
+{
+  CGAL_static_assertion((std::is_convertible<PS_Point, PM_Point>::value));
+  return PM_Point(p);
+}
+
+// just for backward compatibility reasons
+template <typename PM_Point, typename PS_FT>
+PM_Point convert_to_pm_point(const std::array<PS_FT, 3>& p)
+{
+  return PM_Point(p[0], p[1], p[2]);
+}
+
 template <typename PointRange,
          typename PolygonRange,
          typename PointMap = typename CGAL::GetPointMap<PointRange>::const_type>
 class PS_to_PM_converter
 {
  typedef typename boost::range_value<PolygonRange>::type                 Polygon;
-  typedef typename boost::property_traits<PointMap>::value_type           Point;

 public:
  /**
@ -68,8 +82,6 @@ public:

    typedef typename boost::property_traits<VertexPointMap>::value_type     PM_Point;

-    CGAL_static_assertion((std::is_convertible<Point, PM_Point>::value));
-
    reserve(pmesh, static_cast<typename boost::graph_traits<PolygonMesh>::vertices_size_type>(m_points.size()),
            static_cast<typename boost::graph_traits<PolygonMesh>::edges_size_type>(2*m_polygons.size()),
            static_cast<typename boost::graph_traits<PolygonMesh>::faces_size_type>(m_polygons.size()));
@ -94,7 +106,7 @@ public:
        continue;

      vertices[i] = add_vertex(pmesh);
-      PM_Point pi(get(m_pm, m_points[i]));
+      PM_Point pi = convert_to_pm_point<PM_Point>(get(m_pm, m_points[i]));
      put(vpm, vertices[i], pi);
    }

--- 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
@ -118,43 +118,21 @@ bool simplify_polygon(PointRange& points,
 {
  const std::size_t ini_polygon_size = polygon.size();

-  // Start at the last since if two points are identical, the second one gets removed.
-  // By starting at 'last', we ensure that 'to_remove' is ordered from closest to .begin()
-  // to closest to .end()
-  std::size_t last = ini_polygon_size - 1, i = last;
-  bool stop = false;
-  std::vector<std::size_t> to_remove;
-
-  do
+  for(std::size_t i=0; i<polygon.size(); ++i)
  {
-    std::size_t next_i = (i == last) ? 0 : i+1;
-    stop = (next_i == last);
+    const std::size_t s = polygon.size();
+    if(s == 1)
+      break;

-    while(polygon[i] == polygon[next_i] || // combinatorial equality
-          traits.equal_3_object()(points[polygon[i]], points[polygon[next_i]])) // geometric equality
+    const std::size_t ni = (i + 1) % s;
+    if(polygon[i] == polygon[ni] ||
+       traits.equal_3_object()(points[polygon[i]], points[polygon[ni]]))
    {
-      to_remove.push_back(next_i);
 #ifdef CGAL_PMP_REPAIR_POLYGON_SOUP_VERBOSE_PP
-      std::cout << "Duplicate point: polygon[" << next_i << "] = " << polygon[next_i] << std::endl;
+      std::cout << "Duplicate point: polygon[" << ni << "] = " << polygon[ni] << std::endl;
 #endif
-      next_i = (next_i == last) ? 0 : next_i+1;
-
-      // Every duplicate in front of 'last' (circularly-speaking) has already been cleared
-      if(next_i == last)
-      {
-        stop = true;
-        break;
-      }
+      polygon.erase(polygon.begin() + i--);
    }
-
-    i = next_i;
-  }
-  while(!stop);
-
-  while(!to_remove.empty())
-  {
-    polygon.erase(polygon.begin() + to_remove.back());
-    to_remove.pop_back();
  }

  const std::size_t removed_points_n = ini_polygon_size - polygon.size();
--- a/Polygon_mesh_processing/include/CGAL/Side_of_triangle_mesh.h
+++ b/Polygon_mesh_processing/include/CGAL/Side_of_triangle_mesh.h
@ -242,6 +242,80 @@ public:
    }
  }

+#ifndef DOXYGEN_RUNNING
+  template <class K2>
+  Bounded_side operator()(const typename K2::Point_3& point, const K2& k2) const
+  {
+    if(point.x() < box.xmin()
+       || point.x() > box.xmax()
+       || point.y() < box.ymin()
+       || point.y() > box.ymax()
+       || point.z() < box.zmin()
+       || point.z() > box.zmax())
+    {
+      return CGAL::ON_UNBOUNDED_SIDE;
+    }
+
+#ifdef CGAL_HAS_THREADS
+    AABB_tree_* tree_ptr =
+      const_cast<AABB_tree_*>(atomic_tree_ptr.load(std::memory_order_acquire));
+#endif
+    // Lazily build the tree only when needed
+    if (tree_ptr==nullptr)
+    {
+#ifdef CGAL_HAS_THREADS
+      CGAL_SCOPED_LOCK(tree_mutex);
+      tree_ptr = const_cast<AABB_tree_*>(atomic_tree_ptr.load(std::memory_order_relaxed));
+#endif
+      CGAL_assertion(tm_ptr != nullptr && opt_vpm!=boost::none);
+      if (tree_ptr==nullptr)
+      {
+        tree_ptr = new AABB_tree(faces(*tm_ptr).first,
+                                 faces(*tm_ptr).second,
+                                 *tm_ptr, *opt_vpm);
+        const_cast<AABB_tree_*>(tree_ptr)->build();
+#ifdef CGAL_HAS_THREADS
+      atomic_tree_ptr.store(tree_ptr, std::memory_order_release);
+#endif
+      }
+    }
+
+    typedef typename Kernel_traits<Point>::Kernel K1;
+    typedef typename AABB_tree::AABB_traits AABB_traits;
+    typedef internal::Default_tree_helper<AABB_tree> Helper;
+    Helper helper;
+
+    static const unsigned int seed = 1340818006;
+    CGAL::Random rg(seed); // seed some value for make it easy to debug
+    Random_points_on_sphere_3<typename K2::Point_3> random_point(1.,rg);
+
+    typename K2::Construct_ray_3 ray = k2.construct_ray_3_object();
+    typename K2::Construct_vector_3 vector = k2.construct_vector_3_object();
+
+    do { //retry with a random ray
+      typename K2::Ray_3 query = ray(point, vector(CGAL::ORIGIN,*random_point++));
+
+       std::pair<boost::logic::tribool,std::size_t>
+          status( boost::logic::tribool(boost::logic::indeterminate), 0);
+
+      internal::K2_Ray_3_K1_Triangle_3_traversal_traits<AABB_traits, K1, K2, Helper>
+        traversal_traits(status, tree_ptr->traits(), helper);
+
+      tree_ptr->traversal(query, traversal_traits);
+
+      if ( !boost::logic::indeterminate(status.first) )
+      {
+        if (status.first)
+          return (status.second&1) == 1 ? ON_BOUNDED_SIDE : ON_UNBOUNDED_SIDE;
+        //otherwise the point is on the facet
+        return ON_BOUNDARY;
+      }
+    } while (true);
+    return ON_BOUNDARY; // should never be reached
+  }
+
+#endif
+
 };

 } // namespace CGAL
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_corefine.cpp
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_corefine.cpp
@ -71,11 +71,71 @@ void test(const char* f1, const char* f2)
  assert(P.is_valid());
  assert(Q.is_valid());
 }
+
+void test_no_modifications(const char* f1, const char* f2)
+{
+  std::cout << "Corefining " << f1
+            << " and " << f2 << "\n";
+
+  std::cout << "  with Surface_mesh\n";
+  Surface_mesh sm1, sm2;
+  std::ifstream input(f1);
+  assert(input);
+  input >> sm1;
+  input.close();
+  input.open(f2);
+  assert(input);
+  input >> sm2;
+  input.close();
+  My_visitor<Surface_mesh> sm_v;
+
+  std::size_t nb_v_before1 = vertices(sm1).size();
+  std::size_t nb_v_before2 = vertices(sm2).size();
+
+  CGAL::Polygon_mesh_processing::corefine(sm1, sm2,
+    CGAL::parameters::visitor(sm_v),
+    CGAL::parameters::do_not_modify(true));
+
+  std::size_t nb_v_after1 = vertices(sm1).size();
+  std::size_t nb_v_after2 = vertices(sm2).size();
+
+  assert(sm1.is_valid());
+  assert(sm2.is_valid());
+  assert(nb_v_after2==nb_v_before2);
+
+  assert((*(sm_v.i) != 0) == (nb_v_before1!=nb_v_after1));
+
+  std::cout << "  with Polyhedron_3\n";
+  Polyhedron_3 P, Q;
+  input.open(f1);
+  assert(input);
+  input >> P;
+  input.close();
+  input.open(f2);
+  assert(input);
+  input >> Q;
+  My_visitor<Polyhedron_3> sm_p;
+
+  nb_v_before1 = vertices(P).size();
+  nb_v_before2 = vertices(Q).size();
+  CGAL::Polygon_mesh_processing::corefine(P, Q,
+    CGAL::parameters::visitor(sm_p).do_not_modify(true));
+  nb_v_after1 = vertices(P).size();
+  nb_v_after2 = vertices(Q).size();
+
+  assert(nb_v_after1==nb_v_before1);
+  assert((*(sm_p.i) != 0)  == (nb_v_before2!=nb_v_after2));
+
+  assert(P.is_valid());
+  assert(Q.is_valid());
+}
+
 int main(int argc, char** argv)
 {
  for(int i=0; i< (argc-1)/2;++i)
  {
    test(argv[2*i+1], argv[2*(i+1)]);
    test(argv[2*(i+1)], argv[2*i+1]);
+    test_no_modifications(argv[2*(i+1)], argv[2*i+1]);
  }
 }
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_repair_polygon_soup.cpp
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_repair_polygon_soup.cpp
@ -289,14 +289,24 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
  points.push_back(Point_3(1,3,0)); // #3
  points.push_back(Point_3(0,1,0)); // #4
  points.push_back(Point_3(1,1,0)); // #5
-  points.push_back(Point_3(0,0,0)); // #6
+  points.push_back(Point_3(0,0,0)); // #6 == #0

  // ------
  CGAL_polygon polygon;
-  polygon.push_back(0); polygon.push_back(2); polygon.push_back(4);
+  polygon.push_back(0); polygon.push_back(2); polygon.push_back(4); polygon.push_back(0); polygon.push_back(0);
  polygons.push_back(polygon);

  std::size_t res = PMP::internal::simplify_polygons_in_polygon_soup<K>(points, polygons);
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
+  assert(res == 1 && polygons.back().size() == 3);
+
+  // ------
+  polygon.clear();
+  polygon.push_back(0); polygon.push_back(2); polygon.push_back(4);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::simplify_polygons_in_polygon_soup<K>(points, polygons);
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
  assert(res == 0 && polygons.back().size() == 3);

  // ------
@ -305,6 +315,7 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
  polygons.push_back(polygon);

  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
  assert(res == 1 && polygons.back().size() == 1);

  // ------
@ -313,6 +324,7 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
  polygons.push_back(polygon);

  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
  assert(res == 1 && polygons.back().size() == 1);

  // ------
@ -321,8 +333,18 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
  polygons.push_back(polygon);

  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
  assert(res == 1 && polygons.back().size() == 2);

+  // ------
+  polygon.clear();
+  polygon.push_back(0); polygon.push_back(2); polygon.push_back(0); polygon.push_back(4);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
+  assert(res == 0 && polygons.back().size() == 4);
+
  // ------
  // Now with the same geometric positions, but different combinatorial information
  polygon.clear();
@ -330,6 +352,7 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
  polygons.push_back(polygon);

  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
  assert(res == 1 && polygons.back().size() == 3);

  // ------
@ -338,6 +361,7 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
  polygons.push_back(polygon);

  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
  assert(res == 1 && polygons.back().size() == 2);

  // ------
@ -346,6 +370,7 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
  polygons.push_back(polygon);

  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
  assert(res == 1 && polygons.back().size() == 3);
 }

--- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/C3t3_io_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/C3t3_io_plugin.cpp
@ -127,7 +127,7 @@ Polyhedron_demo_c3t3_binary_io_plugin::load(
      item->setName(fileinfo.baseName());
      item->set_valid(false);

-      if(CGAL::build_triangulation_from_file<C3t3::Triangulation, true>(in, item->c3t3().triangulation()))
+      if(CGAL::build_triangulation_from_file<C3t3::Triangulation, true>(in, item->c3t3().triangulation(), true))
      {
        item->c3t3().rescan_after_load_of_triangulation();
        for( C3t3::Triangulation::Finite_cells_iterator
--- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp
@ -87,6 +87,14 @@ Meshing_thread* cgal_code_mesh_3(QList<const SMesh*> pMeshes,
  std::cerr << " done (" << timer.time() * 1000 << " ms)" << std::endl;

  Scene_c3t3_item* p_new_item = new Scene_c3t3_item(surface_only);
+  if(polylines.empty()) {
+    if(protect_features) {
+      p_new_item->set_sharp_edges_angle(sharp_edges_angle);
+    }
+    else if (protect_borders) {
+      p_new_item->set_detect_borders(true);
+    }
+  }

  QString tooltip = QString("<div>From \"") + filename +
    QString("\" with the following mesh parameters"
--- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Optimization_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Optimization_plugin.cpp
@ -524,10 +524,9 @@ treat_result(Scene_c3t3_item& source_item,
                            static_cast<float>(bbox.ymin() + bbox.ymax())/2.f,
                            static_cast<float>(bbox.zmin() + bbox.zmax())/2.f);

-    result_item.setColor(QColor(59,74,226));
+    result_item.setColor(source_item.color());
    result_item.setRenderingMode(source_item.renderingMode());
    result_item.set_data_item(source_item.data_item());
-
    source_item.setVisible(false);

    const Scene_interface::Item_id index = scene->mainSelectionIndex();
@ -539,6 +538,7 @@ treat_result(Scene_c3t3_item& source_item,
  else
  {
    result_item.update_histogram();
+    result_item.invalidateOpenGLBuffers();

    const Scene_interface::Item_id index = scene->mainSelectionIndex();
    scene->itemChanged(index);
--- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Optimization_plugin_cgal_code.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Optimization_plugin_cgal_code.cpp
@ -141,6 +141,10 @@ Optimizer_thread* cgal_code_optimization(Scene_c3t3_item& c3t3_item,
      return NULL;
    }
    Polyhedral_mesh_domain* sm_domain = new Polyhedral_mesh_domain(*smesh);
+    if(c3t3_item.get_sharp_edges_angle() != -1 )
+      sm_domain->detect_features(c3t3_item.get_sharp_edges_angle());
+    else if(c3t3_item.get_detect_borders())
+      sm_domain->detect_borders();

    // Create thread
    typedef Optimization_function<Polyhedral_mesh_domain,Parameters> Opt_function;
--- a/Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_dialog.ui
+++ b/Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_dialog.ui
@ -43,17 +43,7 @@
     <property name="title">
      <string>Tetrahedral remeshing</string>
     </property>
-     <layout class="QGridLayout" name="gridLayout_2" rowstretch="0,0,0,0" columnstretch="0,0">
-      <item row="3" column="0">
-       <widget class="QLabel" name="protect_label">
-        <property name="text">
-         <string>Protect boundaries</string>
-        </property>
-        <property name="alignment">
-         <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
-        </property>
-       </widget>
-      </item>
+     <layout class="QGridLayout" name="gridLayout_2" rowstretch="0,0,0,0,0" columnstretch="0,0">
      <item row="1" column="1">
       <widget class="QSpinBox" name="nbIterations_spinbox">
        <property name="minimumSize">
@ -70,13 +60,26 @@
      <item row="3" column="1">
       <widget class="QCheckBox" name="protect_checkbox">
        <property name="text">
-         <string/>
+         <string>Protect boundaries</string>
        </property>
        <property name="checked">
         <bool>false</bool>
        </property>
       </widget>
      </item>
+      <item row="1" column="0">
+       <widget class="QLabel" name="nbIterations_label">
+        <property name="text">
+         <string>Number of Main iterations</string>
+        </property>
+        <property name="alignment">
+         <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
+        </property>
+        <property name="buddy">
+         <cstring>nbIterations_spinbox</cstring>
+        </property>
+       </widget>
+      </item>
      <item row="2" column="0">
       <spacer name="verticalSpacer_3">
        <property name="orientation">
@ -93,16 +96,19 @@
        </property>
       </spacer>
      </item>
-      <item row="1" column="0">
-       <widget class="QLabel" name="nbIterations_label">
-        <property name="text">
-         <string>Number of Main iterations</string>
+      <item row="0" column="1">
+       <widget class="QDoubleSpinBox" name="edgeLength_dspinbox">
+        <property name="minimumSize">
+         <size>
+          <width>110</width>
+          <height>0</height>
+         </size>
        </property>
-        <property name="alignment">
-         <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
+        <property name="maximum">
+         <double>1000.000000000000000</double>
        </property>
-        <property name="buddy">
-         <cstring>nbIterations_spinbox</cstring>
+        <property name="singleStep">
+         <double>0.100000000000000</double>
        </property>
       </widget>
      </item>
@ -119,19 +125,23 @@
        </property>
       </widget>
      </item>
-      <item row="0" column="1">
-       <widget class="QDoubleSpinBox" name="edgeLength_dspinbox">
-        <property name="minimumSize">
-         <size>
-          <width>110</width>
-          <height>0</height>
-         </size>
+      <item row="4" column="1">
+       <widget class="QCheckBox" name="smoothEdges_checkBox">
+        <property name="toolTip">
+         <string extracomment="Smoothing of constrained polylines.\nWarning : this may reduce quality of dihedral angles."/>
        </property>
-        <property name="maximum">
-         <double>1000.000000000000000</double>
+        <property name="text">
+         <string>Smooth constrained edges</string>
        </property>
-        <property name="singleStep">
-         <double>0.100000000000000</double>
+       </widget>
+      </item>
+      <item row="3" column="0">
+       <widget class="QLabel" name="options_label">
+        <property name="text">
+         <string>Options :</string>
+        </property>
+        <property name="alignment">
+         <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
        </property>
       </widget>
      </item>
--- a/Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_plugin.cpp
@ -3,7 +3,6 @@
 //#define CGAL_TETRAHEDRAL_REMESHING_DEBUG
 //#define CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS
 //#define CGAL_TETRAHEDRAL_REMESHING_PROFILE
-//#define CGAL_TETRAHEDRAL_REMESHING_SMOOTH_SHARP_EDGES

 #include <QtCore/qglobal.h>

@ -99,6 +98,7 @@ public Q_SLOTS:
      double target_length = ui.edgeLength_dspinbox->value();
      unsigned int nb_iter = ui.nbIterations_spinbox->value();
      bool protect = ui.protect_checkbox->isChecked();
+      bool smooth_edges = ui.smoothEdges_checkBox->isChecked();

      // wait cursor
      QApplication::setOverrideCursor(Qt::WaitCursor);
@ -110,7 +110,8 @@ public Q_SLOTS:
          c3t3_item->c3t3(),
          target_length,
          CGAL::parameters::remesh_boundaries(!protect)
-          .number_of_iterations(nb_iter));
+          .number_of_iterations(nb_iter)
+          .smooth_constrained_edges(smooth_edges));

      std::cout << "Remeshing done (" << time.elapsed() << " ms)" << std::endl;

@ -204,6 +205,10 @@ private:
    ui.nbIterations_spinbox->setValue(1);

    ui.protect_checkbox->setChecked(false);
+    ui.smoothEdges_checkBox->setChecked(false);
+
+    connect(ui.protect_checkbox, SIGNAL(toggled(bool)),
+            ui.smoothEdges_checkBox, SLOT(setDisabled(bool)));

    return ui;
  }
--- a/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp
+++ b/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp
@ -347,6 +347,8 @@ struct Scene_c3t3_item_priv {
    cnc_are_shown = false;
    is_aabb_tree_built = false;
    alphaSlider = NULL;
+    sharp_edges_angle = -1;
+    detect_borders = false;
  }
  void computeIntersection(const Primitive& facet);
  void fill_aabb_tree() {
@ -501,6 +503,9 @@ struct Scene_c3t3_item_priv {
  bool is_valid;
  bool is_surface;
  bool last_intersection;
+  //only for optimizers
+  double sharp_edges_angle;
+  bool detect_borders;

  void push_normal(std::vector<float>& normals, const EPICK::Vector_3& n) const
  {
@ -2106,5 +2111,12 @@ Scene_c3t3_item* Scene_c3t3_item::clone() const
  return new Scene_c3t3_item(d->c3t3, d->is_surface);
 }

+void Scene_c3t3_item::set_sharp_edges_angle(double a) { d->sharp_edges_angle = a; }
+double Scene_c3t3_item::get_sharp_edges_angle() { return d->sharp_edges_angle; }
+
+void Scene_c3t3_item::set_detect_borders(bool b) { d->detect_borders = b;}
+bool Scene_c3t3_item::get_detect_borders() { return d->detect_borders; }
+
+
 #include "Scene_c3t3_item.moc"

--- a/Polyhedron/demo/Polyhedron/Scene_c3t3_item.h
+++ b/Polyhedron/demo/Polyhedron/Scene_c3t3_item.h
@ -155,6 +155,12 @@ public:

  QColor get_histogram_color(const double v) const;

+  void set_sharp_edges_angle(double d);
+  double get_sharp_edges_angle();
+
+  void set_detect_borders(bool b);
+  bool get_detect_borders();
+
  void itemAboutToBeDestroyed(Scene_item *) Q_DECL_OVERRIDE;

  void initializeBuffers(Viewer_interface *) const Q_DECL_OVERRIDE;
--- a/Polyhedron/demo/Polyhedron/Viewer.cpp
+++ b/Polyhedron/demo/Polyhedron/Viewer.cpp
@ -1596,6 +1596,7 @@ void Viewer_impl::showDistance(QPoint pixel)
                                        .arg(BPoint.x/scaler.x()-viewer->offset().x, 0, 'g', 10)
                                        .arg(BPoint.y/scaler.y()-viewer->offset().y, 0, 'g', 10)
                                        .arg(BPoint.z/scaler.z()-viewer->offset().z, 0, 'g', 10), true, font, Qt::red, true);
+
        distance_text.append(BCoord);
        CGAL::qglviewer::Vec centerPoint = 0.5*(BPoint+APoint);
        TextItem *centerCoord = new TextItem(float(centerPoint.x),
--- a/Polyhedron/examples/Polyhedron/poly_derivation.cpp
+++ b/Polyhedron/examples/Polyhedron/poly_derivation.cpp
@ -0,0 +1,43 @@
+#include <iostream>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+#include <CGAL/boost/graph/generators.h>
+
+
+template <class Traits>
+struct Mesh: public CGAL::Polyhedron_3<Traits> {
+  std::string name;
+};
+
+#define CGAL_GRAPH_TRAITS_INHERITANCE_TEMPLATE_PARAMS typename Traits
+#define CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME Mesh<Traits>
+#define CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME CGAL::Polyhedron_3<Traits>
+#include <CGAL/boost/graph/graph_traits_inheritance_macros.h>
+
+typedef CGAL::Simple_cartesian<double> K;
+typedef K::Point_3 Point_3;
+
+int main()
+{
+  Mesh<K> mesh;
+  CGAL::make_triangle(Point_3(0,0,0), Point_3(1,0,0), Point_3(1,1,1), mesh);
+  typedef boost::graph_traits<Mesh<K>>::vertex_descriptor vertex_descriptor;
+
+  typedef boost::property_map<Mesh<K>,CGAL::vertex_point_t>::type Point_property_map;
+  Point_property_map ppm = get(CGAL::vertex_point, mesh);
+
+  for(vertex_descriptor vd : vertices(mesh)){
+    if (vd != boost::graph_traits<Mesh<K>>::null_vertex()){
+      std::cout << get(ppm, vd) << std::endl;
+    }
+  }
+  std::cout << CGAL::Polygon_mesh_processing::bbox(mesh) << std::endl;
+
+  return 0;
+}
+
+
+
+
+
--- a/Scripts/developer_scripts/cgal_create_release_with_cmake.cmake
+++ b/Scripts/developer_scripts/cgal_create_release_with_cmake.cmake
@ -86,7 +86,7 @@ if (CGAL_VERSION_FOUND)
  if (NOT CGAL_VERSION)
    set(CGAL_VERSION "${CGAL_VERSION_INPUT}")
  endif()
-  set (GITHUB_PREFIX "https://github.com/CGAL/cgal/blob/releases/CGAL-${CGAL_VERSION}")
+  set (GITHUB_PREFIX "https://github.com/CGAL/cgal/blob/v${CGAL_VERSION}")
 else()
  message(FATAL_ERROR "Cannot extract CGAL version number.")
 endif()
--- a/Scripts/developer_scripts/tag_pr_per_release.sh
+++ b/Scripts/developer_scripts/tag_pr_per_release.sh
@ -32,7 +32,7 @@ REMOTE=`git config branch.releases/CGAL-${PREVIOUS_MAJOR_RELEASE}-branch.remote`
 # refs/pull/*/head as well.
 git fetch --tags "${REMOTE}" `git config "remote.${REMOTE}.fetch"` 'refs/pull/*/head:refs/pull/*/head'

-PR_LIST=`git log --pretty='%D' releases/CGAL-${PREVIOUS_MAJOR_RELEASE}..releases/CGAL-${CURRENT_RELEASE} | awk 'match($0, /refs\/pull\/([0-9]+)\/head/, a) {print a[1]}' | sort -u`
+PR_LIST=`git log --pretty='%D' v${PREVIOUS_MAJOR_RELEASE}..v${CURRENT_RELEASE} | awk 'match($0, /refs\/pull\/([0-9]+)\/head/, a) {print a[1]}' | sort -u`

 for i in ${PR_LIST}; do
  echo ghi label $i -a Merged_in_${CURRENT_RELEASE} -- CGAL/cgal
--- a/Stream_support/include/CGAL/IO/Color.h
+++ b/Stream_support/include/CGAL/IO/Color.h
@ -221,8 +221,6 @@ public:
  /*!
    replaces the rgb values of the colors by the conversion to rgb of
    the hsv values given as parameters.
-
-    Double values given as parameters should take range between 0 and 1.
  */
  void set_hsv (double hue,
                double saturation,
@ -232,8 +230,8 @@ public:
    saturation /= 100.;
    value /= 100.;
    double C = value*saturation;
-    int hh = (int)(hue/60.);
-    double X = C * (1-std::abs (hh % 2 - 1));
+    double hh = (hue/60.);
+    double X = C * (1-std::abs(std::fmod(hh, 2) - 1));
    double r = 0, g = 0, b = 0;

    if( hh>=0 && hh<1 )
--- a/Stream_support/include/CGAL/IO/write_vtk.h
+++ b/Stream_support/include/CGAL/IO/write_vtk.h
@ -13,8 +13,6 @@
 #ifndef CGAL_WRITE_VTK_IO_H
 #define CGAL_WRITE_VTK_IO_H

-#include <CGAL/license/Polyhedron.h>
-
 #include <fstream>
 #include <vector>
 template <class FT>
--- a/Surface_mesh/examples/Surface_mesh/sm_derivation.cpp
+++ b/Surface_mesh/examples/Surface_mesh/sm_derivation.cpp
@ -16,26 +16,9 @@ namespace My {
 } // namespace My


-namespace boost {
-
-  template <>
-  struct graph_traits<My::Mesh>
-    : public boost::graph_traits<My::Mesh::Base>
-  {};
-
-  template <typename T>
-  struct property_map<My::Mesh, T>
-    : public boost::property_map<My::Mesh::Base, T>
-  {};
-
-}
-
-namespace CGAL{
-  template <typename T>
-  struct graph_has_property<My::Mesh, T>
-      : public CGAL::graph_has_property<My::Mesh::Base, T>
-  {};
-}
+#define CGAL_GRAPH_TRAITS_INHERITANCE_CLASS_NAME My::Mesh
+#define CGAL_GRAPH_TRAITS_INHERITANCE_BASE_CLASS_NAME CGAL::Surface_mesh<::Point_3>
+#include <CGAL/boost/graph/graph_traits_inheritance_macros.h>

 int main()
 {
--- a/Tetrahedral_remeshing/doc/Tetrahedral_remeshing/Tetrahedral_remeshing.txt
+++ b/Tetrahedral_remeshing/doc/Tetrahedral_remeshing/Tetrahedral_remeshing.txt
@ -100,12 +100,34 @@ setting the named parameter `remesh_boundaries` to `false`.

 The tetrahedral remeshing algorithm is designed as a post-processing for
 mesh generation algorithms. The API allows to generate a tetrahedral mesh
-with `Mesh_3`, the tetrahedral mesh generation package, and
+with the \cgal \ref PkgMesh3 package, and
 further improve it with the tetrahedral remeshing algorithm.
 This example shows how to use tetrahedral mesh generation and remeshing in sequence,
 from a polyhedral domain with features.

 \cgalExample{Tetrahedral_remeshing/mesh_and_remesh_polyhedral_domain_with_features.cpp}

+
+\subsection ssecEx5 Tetrahedral Remeshing from Any Tetrahedral Mesh
+
+The following example shows how to read a mesh from a triangulation stored in a
+Medit file, perform tetrahedral remeshing, and save the output triangulation.
+The input triangulation should follow the validity requirements of a
+`CGAL::Triangulation_3` (valid connectivity, positive orientation of the cells, and
+coverage of the convex hull of the vertices).
+
+\cgalExample{Tetrahedral_remeshing/tetrahedral_remeshing_from_mesh.cpp}
+
+
+\section secTetRemeshingHistory Implementation History
+
+This package implements the uniform version of the "Multi-Material Adaptive Volume Remesher"
+algorithm for quality tetrahedral remeshing, described by Noura Faraj et al. in \cgalCite{faraj2016mvr}.
+
+A first version of the code was written by Noura Faraj, Jean-Marc Thiery, and Tamy Boubekeur.
+Jane Tournois worked on the finalization of the code, the API, and documentation.
+
+It was initially published in CGAL-5.1.
+
 */
 } /* namespace CGAL */
--- a/Tetrahedral_remeshing/doc/Tetrahedral_remeshing/examples.txt
+++ b/Tetrahedral_remeshing/doc/Tetrahedral_remeshing/examples.txt
@ -4,4 +4,5 @@
 \example Tetrahedral_remeshing/tetrahedral_remeshing_of_one_subdomain.cpp
 \example Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp
 \example Tetrahedral_remeshing/mesh_and_remesh_polyhedral_domain_with_features.cpp
+\example Tetrahedral_remeshing/tetrahedral_remeshing_from_mesh.cpp
 */
--- a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt
+++ b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt
@ -28,6 +28,7 @@ include(CGAL_Eigen_support)
 create_single_source_cgal_program( "tetrahedral_remeshing_example.cpp" )
 create_single_source_cgal_program( "tetrahedral_remeshing_with_features.cpp")
 create_single_source_cgal_program( "tetrahedral_remeshing_of_one_subdomain.cpp")
+create_single_source_cgal_program( "tetrahedral_remeshing_from_mesh.cpp")

 create_single_source_cgal_program( "mesh_and_remesh_polyhedral_domain_with_features.cpp" )
 target_link_libraries(mesh_and_remesh_polyhedral_domain_with_features PUBLIC CGAL::Eigen_support)
--- a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/data/sphere.mesh
+++ b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/data/sphere.mesh
--- a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/tetrahedral_remeshing_example.cpp
+++ b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/tetrahedral_remeshing_example.cpp
@ -3,7 +3,6 @@
 #include <CGAL/Tetrahedral_remeshing/Remeshing_triangulation_3.h>
 #include <CGAL/tetrahedral_remeshing.h>

-#include <CGAL/Tetrahedral_remeshing/tetrahedral_remeshing_io.h>
 #include "tetrahedral_remeshing_generate_input.h"

 #include <iostream>
--- a/Show More
+++ b/Show More