Merge remote-tracking branch 'cgal/master' into Iterative_authalic_parameterization

AABB_tree/doc/AABB_tree/aabb_tree.txt

@@ -189,7 +189,7 @@ example illustrates this for two polyhedral surfaces.
 \subsection aabb_tree_examples_7 Trees of Custom Primitives
 
 The AABB tree example folder contains three examples of trees
-constructed with customize primitives. In \ref AABB_tree/AABB_custom_example.cpp "AABB_custom_example.cpp"
+constructed with custom primitives. In \ref AABB_tree/AABB_custom_example.cpp "AABB_custom_example.cpp"
 the primitive contains triangles which are defined by three pointers
 to custom points. In \ref AABB_tree/AABB_custom_triangle_soup_example.cpp "AABB_custom_triangle_soup_example.cpp" all input
 triangles are stored into a single array so as to form a triangle

AABB_tree/include/CGAL/AABB_face_graph_triangle_primitive.h

@@ -112,7 +112,7 @@ public:
   /*!
   Id type:
   - `boost::graph_traits<FaceGraph>::%face_descriptor` if `OneFaceGraphPerTree` is `CGAL::Tag_true`
-  - `std::pair<boost::graph_traits<FaceGraph>::%face_descriptor, FaceGraph>` if `OneFaceGraphPerTree` is `CGAL::Tag_false`
+  - `std::pair<boost::graph_traits<FaceGraph>::%face_descriptor, const FaceGraph*>` if `OneFaceGraphPerTree` is `CGAL::Tag_false`
   */
   unspecified_type Id;
 

AABB_tree/include/CGAL/AABB_halfedge_graph_segment_primitive.h

@@ -125,8 +125,8 @@ public:
   typedef Kernel_traits<Point>::Kernel::Segment_3 Datum;
   /*!
   Id type:
-  - `boost::graph_traits<HalfedgeGraph>::%edge_descriptor if `OneHalfedgeGraphPerTree` is `Tag_true`
+  - `boost::graph_traits<HalfedgeGraph>::%edge_descriptor` if `OneHalfedgeGraphPerTree` is `Tag_true`
-  - `std::pair<boost::graph_traits<HalfedgeGraph>::edge_descriptor, HalfedgeGraph>` if `OneHalfedgeGraphPerTree` is `Tag_false`
+  - `std::pair<boost::graph_traits<HalfedgeGraph>::%edge_descriptor, const HalfedgeGraph*>` if `OneHalfedgeGraphPerTree` is `Tag_false`
   */
   unspecified_type Id;
   /// @}
@@ -202,4 +202,3 @@ public:
 #include <CGAL/enable_warnings.h>
 
 #endif // CGAL_AABB_HALFEDGE_GRAPH_SEGMENT_PRIMITIVE_H
-

AABB_tree/include/CGAL/AABB_traits.h

@@ -171,9 +171,11 @@ class AABB_tree;
 /// \sa `AABBPrimitiveWithSharedData`
 
   template<typename GeomTraits, typename AABBPrimitive, typename BboxMap = Default>
-class AABB_traits:
+class AABB_traits
-  public internal::AABB_tree::AABB_traits_base<AABBPrimitive>,
+#ifndef DOXYGEN_RUNNING
+: public internal::AABB_tree::AABB_traits_base<AABBPrimitive>,
   public internal::AABB_tree::AABB_traits_base_2<GeomTraits>
+#endif
 {
   typedef typename CGAL::Object Object;
 public:

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.58 or later</b>
+<b>Version 1.62 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

Installation/CMakeLists.txt

@@ -382,8 +382,6 @@ include(${CGAL_MODULES_DIR}/CGAL_Macros.cmake)
 include(${CGAL_MODULES_DIR}/CGAL_enable_end_of_configuration_hook.cmake)
 cgal_setup_module_path()
 
-message( STATUS "CGAL_REFERENCE_CACHE_DIR=${CGAL_REFERENCE_CACHE_DIR}" )
-
 if ( RUNNING_CGAL_AUTO_TEST )
   message(STATUS "Operating system:")
   execute_process(COMMAND uname -a
@@ -392,41 +390,6 @@ if ( RUNNING_CGAL_AUTO_TEST )
     ERROR_VARIABLE uname_a)
   message(STATUS "${uname_a}")
   CGAL_display_compiler_version()
-  if ( NOT "${CGAL_REFERENCE_CACHE_DIR}" STREQUAL "" )
-    if ( EXISTS ${CGAL_REFERENCE_CACHE_DIR} )
-      if ( EXISTS ${CGAL_REFERENCE_CACHE_DIR}/CMakeCache.txt )
-        message( STATUS "Loading reference cache from ${CGAL_REFERENCE_CACHE_DIR}" )
-        load_cache( ${CGAL_REFERENCE_CACHE_DIR}
-          EXCLUDE CGAL_Core_LIBRARY
-                  CGAL_CORE_PACKAGE_DIR
-                  WITH_CGAL_Core
-                  CGAL_INSTALLATION_PACKAGE_DIR
-                  CGAL_MAINTENANCE_PACKAGE_DIR
-                  CGAL_PDB_BINARY_DIR
-                  CGAL_PDB_SOURCE_DIR
-                  CGAL_BINARY_DIR
-                  CGAL_SOURCE_DIR)
-#        message("List of cache variables:")
-
-        ## The following lines removes nasty loaded cache values. We do not
-        ## want that the current build tree depends on binaries that were
-        ## build in the reference build tree.
-        get_property(cache_variables DIRECTORY PROPERTY CACHE_VARIABLES)
-        foreach(var ${cache_variables})
-#          get_property(var_value CACHE ${var} PROPERTY VALUE)
-#          get_property(type CACHE ${var} PROPERTY TYPE)
-          string(REGEX MATCH "^CGAL(_.*_(DEPENDS|BINARY_DIR)|_.*LIBRARY)$" var_name_matches ${var})
-          if(var_name_matches)
-            unset(${var} CACHE)
-#          else()
-#            message("${var}:${var_type}=${var_value}")
-          endif()
-        endforeach()
-
-
-      endif()
-    endif()
-  endif()
 endif()
 
 include(CGAL_Common)

Installation/cmake/modules/CGAL_add_test.cmake

@@ -16,6 +16,9 @@ if(NOT POLICY CMP0064)
   # Add a fake function to avoid CMake errors
   function(cgal_add_compilation_test)
   endfunction()
+  # Add a fake function to avoid CMake errors
+  function(cgal_setup_test_properties)
+  endfunction()
 
   # Then return, to exit the file
   return()
@@ -213,7 +216,7 @@ function(cgal_setup_test_properties test_name)
         PROPERTY
         WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/__exec_test_dir)
     endif()
-    
+
     set_property(TEST "${test_name}"
       APPEND PROPERTY FIXTURES_REQUIRED "${PROJECT_NAME}")
     if(exe_name)

Installation/cmake/modules/FindMPFI.cmake

@@ -1,4 +1,5 @@
 find_package( GMP QUIET )
+find_package( MPFR QUIET )
 
 if( GMP_FOUND AND MPFR_FOUND )
 

Installation/include/CGAL/config.h

@@ -36,8 +36,9 @@
 #  define WIN64
 #endif
 
-#ifdef BOOST_MSVC
+#ifdef _MSC_VER
 #define _SILENCE_CXX17_ALLOCATOR_VOID_DEPRECATION_WARNING 1
+#define _SILENCE_CXX17_OLD_ALLOCATOR_MEMBERS_DEPRECATION_WARNING 1
 #endif
 
 #ifdef CGAL_INCLUDE_WINDOWS_DOT_H
@@ -296,6 +297,10 @@
 #if __cplusplus >= 201402L || _MSVC_LANG >= 201402L
 #  define CGAL_CXX14 1
 #endif
+// Same for C++17
+#if __cplusplus >= 201703L || _MSVC_LANG >= 201703L
+#  define CGAL_CXX17 1
+#endif
 
 #if defined(BOOST_NO_CXX11_HDR_FUNCTIONAL) || BOOST_VERSION < 105000
 #define CGAL_CFG_NO_STD_HASH 1
@@ -559,7 +564,7 @@ using std::max;
 #define CGAL_NORETURN  [[noreturn]]
 
 // Macro to specify [[no_unique_address]] if supported
-#if __has_cpp_attribute(no_unique_address)
+#if CGAL_CXX11 && __has_cpp_attribute(no_unique_address)
 #  define CGAL_NO_UNIQUE_ADDRESS [[no_unique_address]]
 #else
 #  define CGAL_NO_UNIQUE_ADDRESS

Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Bounded_3_do_intersect.h

@@ -161,23 +161,40 @@ do_intersect(const typename K::Tetrahedron_3 &tet,
 }
 
 template <class K>
-inline
+inline typename K::Boolean do_intersect(const typename K::Tetrahedron_3 &tet,
-typename K::Boolean
+                                        const CGAL::Bbox_3 &bb, const K &k) {
-do_intersect(const typename K::Tetrahedron_3 &tet,
+  // Swap arguments.
-             const CGAL::Bbox_3 &bb,
+  return do_intersect(bb, tet, k);
-             const K & k)
-{
-  return do_intersect_tetrahedron_bounded(bb, tet, typename K::Point_3(bb.xmin(), bb.ymin(), bb.zmin()), k);
 }
 
-  template <class K>
+// BBox_3 specific code since it is ok for BBox_3 to degenerate.
-inline
+template <class K>
-typename K::Boolean
+inline typename K::Boolean do_intersect(const CGAL::Bbox_3 &aabb,
-do_intersect(const CGAL::Bbox_3 &bb,
+                                        const typename K::Tetrahedron_3 &tet,
-             const typename K::Tetrahedron_3 &tet,
+                                        const K &k) {
-             const K & k)
+  typename K::Construct_triangle_3 tr = k.construct_triangle_3_object();
-{
+  typename K::Boolean result = false;
-  return do_intersect_tetrahedron_bounded(bb, tet, typename K::Point_3(bb.xmin(), bb.ymin(), bb.zmin()), k);
+  typename K::Boolean b = false;
+
+  b = do_intersect(aabb, tr(tet[0], tet[1], tet[2]), k);
+  if (certainly(b)) return b;
+  if (is_indeterminate(b)) result = b;
+  b = do_intersect(aabb, tr(tet[1], tet[2], tet[3]), k);
+  if (certainly(b)) return b;
+  if (is_indeterminate(b)) result = b;
+  b = do_intersect(aabb, tr(tet[2], tet[3], tet[0]), k);
+  if (certainly(b)) return b;
+  if (is_indeterminate(b)) result = b;
+  b = do_intersect(aabb, tr(tet[3], tet[0], tet[1]), k);
+  if (certainly(b)) return b;
+  if (is_indeterminate(b)) result = b;
+
+  b = k.has_on_bounded_side_3_object()(
+      tet, k.construct_point_3_object()(aabb.xmin(), aabb.ymin(), aabb.zmin()));
+  if (certainly(b)) return b;
+  if (is_indeterminate(b)) result = b;
+
+  return result;
 }
 
 } // namespace internal

Maintenance/infrastructure/cgal.geometryfactory.com/crontab

@@ -29,7 +29,7 @@ LC_CTYPE=en_US.UTF-8
 # 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
 # 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 --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 $HOME/CGAL/branches/empty-dir --public --do-it || echo ERROR
 
 ## Older stuff
 # from branch 4.13

Maintenance/test_handling/to_zipped_format

@@ -69,7 +69,7 @@ sub reformat_results($)
 	if (/LEDAWIN_VERSION = '([^']+)'/) {
 	    $LEDA_VERSION="$LEDA_VERSION+win";
 	}
-	if (/^COMPILER_VERSION = '([^']+)'/) {
+	if (/COMPILER_VERSION = '([^']+)'/) {
 	    $COMPILER = $1;
 	}
 	if (/^TESTER_NAME\s+(.*)$/) {

Mesh_3/examples/Mesh_3/CMakeLists.txt

@@ -133,23 +133,26 @@ else()
   message( STATUS "NOTICE: Some examples need the CGAL_ImageIO library, and will not be compiled." )
 endif()
 
-if(CGAL_ACTIVATE_CONCURRENT_MESH_3 AND TARGET CGAL::TBB_support AND TARGET ${target})
+if(CGAL_ACTIVATE_CONCURRENT_MESH_3 AND TARGET CGAL::TBB_support)
   foreach(target
-      mesh_3D_image_with_features
       mesh_3D_image
-      mesh_polyhedral_domain
+      mesh_3D_image_variable_size
       mesh_3D_image_with_custom_initialization
-      mesh_two_implicit_spheres_with_balls
+      mesh_3D_image_with_features
-      mesh_optimization_lloyd_example
-      mesh_optimization_example
       mesh_implicit_sphere
-      mesh_polyhedral_complex_sm
       mesh_implicit_sphere_variable_size
-      mesh_polyhedral_domain_sm
+      mesh_optimization_example
-      mesh_polyhedral_domain_with_lipschitz_sizing
+      mesh_optimization_lloyd_example
       mesh_polyhedral_complex
+      mesh_polyhedral_complex_sm
+      mesh_polyhedral_domain
+      mesh_polyhedral_domain_sm
       mesh_polyhedral_domain_with_features
-      mesh_3D_image_variable_size)
+      mesh_polyhedral_domain_with_features_sm
-    target_link_libraries(${target} PUBLIC CGAL::TBB_support)
+      mesh_polyhedral_domain_with_lipschitz_sizing
+      mesh_two_implicit_spheres_with_balls)
+    if(TARGET ${target})
+      target_link_libraries(${target} PUBLIC CGAL::TBB_support)
+    endif()
   endforeach()
 endif()

Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h

@@ -1572,10 +1572,17 @@ private:
   Cell_vector c3t3_cells(const Cell_vector& cells) const
   {
     Cell_vector c3t3_cells;
+#ifdef CGAL_CXX17
+    std::remove_copy_if(cells.begin(),
+                        cells.end(),
+                        std::back_inserter(c3t3_cells),
+                        std::not_fn(Is_in_c3t3<Cell_handle>(c3t3_)));
+#else
     std::remove_copy_if(cells.begin(),
                         cells.end(),
                         std::back_inserter(c3t3_cells),
                         std::not1(Is_in_c3t3<Cell_handle>(c3t3_)) );
+#endif
     return c3t3_cells;
   }
 

Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h

@@ -1155,6 +1155,7 @@ insert_balls(const Vertex_handle& vp,
   //   n = 2(d-sq) / (sp+sq)
   // =======================
 
+  const FT d_signF = static_cast<FT>(d_sign);
   int n = static_cast<int>(std::floor(FT(2)*(d-sq) / (sp+sq))+.5);
   // if( minimal_weight_ != 0 && n == 0 ) return;
 
@@ -1180,7 +1181,7 @@ insert_balls(const Vertex_handle& vp,
       const Bare_point new_point =
         domain_.construct_point_on_curve(cp(vp_wp),
                                          curve_index,
-                                         d_sign * d / 2);
+                                         d_signF * d / 2);
       const int dim = 1; // new_point is on edge
       const Index index = domain_.index_from_curve_index(curve_index);
       const FT point_weight = CGAL::square(size_(new_point, dim, index));
@@ -1230,7 +1231,6 @@ insert_balls(const Vertex_handle& vp,
   FT norm_step_size = dleft_frac * step_size;
 
   // Initial distance
-  FT d_signF = static_cast<FT>(d_sign);
   FT pt_dist = d_signF * norm_step_size;
   Vertex_handle prev = vp;
   const Bare_point& p = cp(c3t3_.triangulation().point(vp));

Mesh_3/include/CGAL/Mesh_3/tet_soup_to_c3t3.h

@@ -105,13 +105,14 @@ bool build_finite_cells(Tr& tr,
     const std::vector<typename Tr::Vertex_handle>& vertex_handle_vector,
     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>, int>& border_facets,
+    const std::map<std::array<int,3>, typename Tr::Cell::Surface_patch_index>& border_facets,
     const bool verbose)
 {
   typedef std::array<int, 5>              Tet_with_ref; // 4 ids + 1 reference
 
   typedef typename Tr::Vertex_handle                            Vertex_handle;
   typedef typename Tr::Cell_handle                              Cell_handle;
+  typedef typename Tr::Cell::Surface_patch_index                Surface_patch_index;
 
   CGAL_assertion_code(
     typename Tr::Geom_traits::Construct_point_3 cp =
@@ -171,7 +172,8 @@ bool build_finite_cells(Tr& tr,
           ++k;
         } while(f[0] != n0);
 
-        typename std::map<std::array<int,3>, int>::const_iterator it = border_facets.find(f);
+        typename std::map<std::array<int,3>, Surface_patch_index>::const_iterator
+          it = border_facets.find(f);
         if(it != border_facets.end())
         {
           c->set_surface_patch_index(j, it->second);
@@ -186,7 +188,7 @@ bool build_finite_cells(Tr& tr,
           if(it != border_facets.end())
             c->set_surface_patch_index(j, it->second);
           else
-            c->set_surface_patch_index(j, 0);
+            c->set_surface_patch_index(j, Surface_patch_index());
         }
       }
     }
@@ -259,6 +261,8 @@ bool build_infinite_cells(Tr& tr,
     // the only finite facet
     it->second.push_back(std::make_pair(opp_c, 0));
     CGAL_assertion(it->second.size() == 2);
+
+    opp_c->set_surface_patch_index(0, c->surface_patch_index(i));
   }
 
 #ifdef CGAL_TET_SOUP_TO_C3T3_DEBUG
@@ -322,7 +326,7 @@ template<class Tr, bool c3t3_loader_failed>
 bool build_triangulation(Tr& tr,
                          const std::vector<typename Tr::Point>& points,
                          const std::vector<std::array<int,5> >& finite_cells,
-                         const std::map<std::array<int,3>, int>& border_facets,
+                         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)
 {
@@ -389,7 +393,7 @@ bool build_triangulation_from_file(std::istream& is,
 
   std::vector<Tet_with_ref> finite_cells;
   std::vector<Point_3> points;
-  std::map<Facet, int> border_facets;
+  std::map<Facet, typename Tr::Cell::Surface_patch_index> border_facets;
 
   // grab the vertices
   int dim;
@@ -420,7 +424,8 @@ bool build_triangulation_from_file(std::istream& is,
       is >> nf;
       for(int i=0; i<nf; ++i)
       {
-        int n1, n2, n3, surface_patch_id;
+        int n1, n2, n3;
+        typename Tr::Cell::Surface_patch_index surface_patch_id;
         is >> n1 >> n2 >> n3 >> surface_patch_id;
         Facet facet;
         facet[0] = n1 - 1;

Partition_2/include/CGAL/Partition_2/Triangulation_indirect_traits_2.h

@@ -113,6 +113,7 @@ template <class Circulator>
 class Construct_circulator_2
 {
 public:
+   typedef Circulator result_type;
    Circulator operator()(Circulator p1) const { return p1; }
 };
 

Periodic_3_triangulation_3/include/CGAL/internal/Periodic_3_Delaunay_triangulation_remove_traits_3.h

@@ -40,6 +40,14 @@ public:
 
   using Base::operator();
 
+  template<typename F>
+  struct result : Base::template result<F> {};
+
+  template<typename F>
+  struct result<F(Point_3)> {
+    typedef const Point_3& type;
+  };
+
   const Point_3& operator()(const Point_3& p) const { return p; }
 };
 

Periodic_3_triangulation_3/include/CGAL/internal/Periodic_3_regular_triangulation_remove_traits_3.h

@@ -49,6 +49,14 @@ public:
 
   using Base::operator(); // for K::Weighted_point_3 to Point_3
 
+  template<typename F>
+  struct result : Base::template result<F> {};
+
+  template<typename F>
+  struct result<F(Weighted_point_3)> {
+    typedef Point_3 type;
+  };
+
   Point_3 operator()(const Weighted_point_3& wp) const {
     return std::make_pair(operator()(wp.first), wp.second /* offset */);
   }

Point_set_processing_3/include/CGAL/Point_set_processing_3/internal/Neighbor_query.h

@@ -112,7 +112,7 @@ public:
 
   template <typename OutputIterator>
   void get_iterators (const Point& query, unsigned int k, FT neighbor_radius,
-                      OutputIterator output, bool fallback_k_if_sphere_empty = true) const
+                      OutputIterator output, unsigned int fallback_k_if_sphere_empty = 3) const
   {
     if (neighbor_radius != FT(0))
     {
@@ -142,11 +142,10 @@ public:
       catch (const Maximum_points_reached_exception&)
       { }
 
-      // Fallback, if less than 3 points are return, search for the 3
+      // Fallback, if not enough  points are return, search for the knn
       // first points
-      if (fallback_k_if_sphere_empty && nb < 3)
+      if (nb < fallback_k_if_sphere_empty)
-        k = 3;
+        k = fallback_k_if_sphere_empty;
-      // Else, no need to search for K nearest neighbors
       else
         k = 0;
     }
@@ -173,14 +172,14 @@ public:
 
   template <typename OutputIterator>
   void get_points (const Point& query, unsigned int k, FT neighbor_radius,
-                   OutputIterator output) const
+                   OutputIterator output, unsigned int fallback_k_if_sphere_empty = 3) const
   {
     return get_iterators(query, k, neighbor_radius,
                          boost::make_function_output_iterator
                          ([&](const input_iterator& it)
                           {
                             *(output ++) = get (m_point_map, *it);
-                          }));
+                          }), fallback_k_if_sphere_empty);
   }
 };
 

Point_set_processing_3/include/CGAL/cluster_point_set.h

@@ -209,7 +209,7 @@ std::size_t cluster_point_set (PointRange& points,
 
       neighbors.clear();
       neighbor_query.get_iterators (get (point_map, p), 0, neighbor_radius,
-                                    std::back_inserter (neighbors), false);
+                                    std::back_inserter (neighbors), 0);
 
       for (const iterator& it : neighbors)
       {

Point_set_processing_3/include/CGAL/jet_estimate_normals.h

@@ -70,7 +70,10 @@ jet_estimate_normal(const typename NeighborQuery::Point_3& query, ///< point to
   typedef typename Monge_jet_fitting::Monge_form Monge_form;
 
   std::vector<Point> points;
-  neighbor_query.get_points (query, k, neighbor_radius, std::back_inserter(points));
+
+  // query using as fallback minimum requires nb points for jet fitting (d+1)*(d+2)/2
+  neighbor_query.get_points (query, k, neighbor_radius, std::back_inserter(points),
+                             (degree_fitting + 1) * (degree_fitting + 2) / 2);
 
   // performs jet fitting
   Monge_jet_fitting monge_fit;

Point_set_processing_3/include/CGAL/jet_smooth_point_set.h

@@ -75,7 +75,10 @@ jet_smooth_point(
   typedef typename Monge_jet_fitting::Monge_form Monge_form;
 
   std::vector<Point> points;
-  neighbor_query.get_points (query, k, neighbor_radius, std::back_inserter(points));
+
+  // query using as fallback minimum requires nb points for jet fitting (d+1)*(d+2)/2
+  neighbor_query.get_points (query, k, neighbor_radius, std::back_inserter(points),
+                             (degree_fitting + 1) * (degree_fitting + 2) / 2);
 
   // performs jet fitting
   Monge_jet_fitting monge_fit;

Scripts/developer_scripts/autotest_cgal

@@ -490,13 +490,7 @@ build_cgal_on_host()
 
       CGAL_BINARY_DIR="${CGAL_BINARY_DIR_BASE}/${PLATFORM}"
 
-      if [ -d "${REFERENCE_PLATFORMS_DIR}/${PLATFORM}" ] ; then
+      log "${ACTUAL_LOGFILE}" "Building cgal libs on host ${HOST} and platform ${PLATFORM}\nUnder ${CGAL_BINARY_DIR}\n"
-        CGAL_REFERENCE_CACHE_DIR="${REFERENCE_PLATFORMS_DIR}/${PLATFORM}"
-      else
-        CGAL_REFERENCE_CACHE_DIR=""
-      fi
-
-      log "${ACTUAL_LOGFILE}" "Building cgal libs on host ${HOST} and platform ${PLATFORM}\nUnder ${CGAL_BINARY_DIR}\nUsing reference cache directory ${CGAL_REFERENCE_CACHE_DIR}"
 
       if [ -f "${CGAL_BINARY_DIR}/localbuildscript" ] ; then
         log "${ACTUAL_LOGFILE}" "WARNING! Already built on platform ${PLATFORM}."
@@ -514,19 +508,12 @@ if [ -z "\$CMAKE_GENERATOR" ]; then
   CMAKE_GENERATOR='${CMAKE_GENERATOR}';
 fi
 MAKE_CMD='${MAKE_CMD}';
-if [ -n "${IS_CYGWIN}" ]; then
-  CGAL_REFERENCE_CACHE_DIR=\$( cygpath -w '${CGAL_REFERENCE_CACHE_DIR}' );
-else
-  CGAL_REFERENCE_CACHE_DIR='${CGAL_REFERENCE_CACHE_DIR}';\
-fi
 
 export CMAKE_GENERATOR;
 export MAKE_CMD;
 export CGAL_BINARY_DIR;
-export CGAL_REFERENCE_CACHE_DIR;
 cd '${CGAL_BINARY_DIR}';
 cmake \${INIT_FILE:+"-C\${INIT_FILE}"} '${CMAKE_GENERATOR}' -DRUNNING_CGAL_AUTO_TEST=TRUE \\
-                           -DCGAL_REFERENCE_CACHE_DIR="\$CGAL_REFERENCE_CACHE_DIR" \\
                            VERBOSE=1 \\
                            ../../..;
 ${MAKE_CMD} VERBOSE=ON -k -fMakefile  ;

Tetrahedral_remeshing/doc/Tetrahedral_remeshing/PackageDescription.txt

@@ -56,5 +56,6 @@ while preserving the input geometric curve and surface features.
 \cgalCRPSection{Function Templates}
 
 - `CGAL::tetrahedral_isotropic_remeshing()`
+- `CGAL::convert_to_triangulation_3()`
 
 */

Tetrahedral_remeshing/doc/Tetrahedral_remeshing/Tetrahedral_remeshing.txt

@@ -39,17 +39,18 @@ The tetrahedral remeshing algorithm improves the quality of dihedral angles,
 while targeting the user-defined uniform sizing field and preserving the
 topology of the feature complex, as highlighted by Figure \cgalFigureRef{Remesh_liver}.
 
-\cgalFigureBegin{Remesh_liver, tetrahedral_remeshing_before_after.png}
-Tetrahedral mesh, improved by the uniform tetrahedral remeshing algorithm.
-(Left) Before remeshing, dihedral angles were in the interval [0.7; 178.9].
-(Right) After remeshing and keeping the same density,
-dihedral angles are the interval [12,7; 157.7].
-\cgalFigureEnd
-
 Experimental evidence shows that a higher number of remeshing iterations
 leads to a mesh with a improved fidelity to the sizing criterion,
 and higher quality dihedral angles.
 
+\cgalFigureBegin{Remesh_liver, tetrahedral_remeshing_before_after.png}
+Tetrahedral mesh, modified by our uniform tetrahedral remeshing method.
+(Left) Before remeshing, dihedral angles were in the interval [1.3; 177.8].
+(Right) After remeshing and keeping the same density,
+dihedral angles are in the interval [12,7; 157.7].
+\cgalFigureEnd
+
+
 \section secTetRemeshingAPI API
 
 The tetrahedral remeshing algorithm is implemented as a single free function
@@ -97,5 +98,17 @@ setting the named parameter `remesh_boundaries` to `false`.
 
 \cgalExample{Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp }
 
+
+\subsection ssecEx4 Tetrahedral Remeshing After Mesh Generation
+
+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
+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}
+
 */
 } /* namespace CGAL */

Tetrahedral_remeshing/doc/Tetrahedral_remeshing/examples.txt

@@ -3,5 +3,5 @@
 \example Tetrahedral_remeshing/tetrahedral_remeshing_example.cpp
 \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
 */

Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt

@@ -19,9 +19,15 @@ if ( NOT Boost_FOUND )
   return()
 endif()
 
+# Use Eigen for Mesh_3
+find_package(Eigen3 3.1.0 REQUIRED) #(3.1.0 or greater)
+include(CGAL_Eigen_support)
 
 # Creating entries for all C++ files with "main" routine
 # ##########################################################
 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( "mesh_and_remesh_polyhedral_domain_with_features.cpp" )
+target_link_libraries(mesh_and_remesh_polyhedral_domain_with_features PUBLIC CGAL::Eigen_support)

Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_polyhedral_domain_with_features.cpp

@@ -0,0 +1,75 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+
+#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
+#include <CGAL/make_mesh_3.h>
+
+#include <CGAL/tetrahedral_remeshing.h>
+
+// Domain
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;
+typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;
+
+// Triangulation for Meshing
+typedef CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, CGAL::Default>::type Tr;
+typedef CGAL::Mesh_complex_3_in_triangulation_3<
+  Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index> C3t3;
+
+// Criteria
+typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+
+// Triangulation for Remeshing
+typedef CGAL::Triangulation_3<typename Tr::Geom_traits,
+  typename Tr::Triangulation_data_structure> Triangulation_3;
+
+// To avoid verbose function and named parameters call
+using namespace CGAL::parameters;
+
+int main(int argc, char* argv[])
+{
+  const char* fname = (argc > 1) ? argv[1] : "data/fandisk.off";
+  std::ifstream input(fname);
+  Polyhedron polyhedron;
+  input >> polyhedron;
+  if (input.fail()) {
+    std::cerr << "Error: Cannot read file " << fname << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  if (!CGAL::is_triangle_mesh(polyhedron)) {
+    std::cerr << "Input geometry is not triangulated." << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  // Create domain
+  Mesh_domain domain(polyhedron);
+
+  // Get sharp features
+  domain.detect_features();
+
+  // Mesh criteria
+  Mesh_criteria criteria(edge_size = 0.025,
+    facet_angle = 25, facet_size = 0.05, facet_distance = 0.005,
+    cell_radius_edge_ratio = 3, cell_size = 0.05);
+
+  // Mesh generation
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
+
+  Triangulation_3 tr = CGAL::convert_to_triangulation_3(std::move(c3t3));
+  //note we use the move semantic, with std::move(c3t3),
+  //  to avoid a copy of the triangulation by the function
+  //  `CGAL::convert_to_triangulation_3()`
+  //  After the call to this function, c3t3 is an empty and valid C3t3.
+  //It is possible to use :  CGAL::convert_to_triangulation_3(c3t3),
+  //  Then the triangulation is copied and duplicated, and c3t3 remains as is.
+
+  const double target_edge_length = 0.1;//coarsen the mesh
+  CGAL::tetrahedral_isotropic_remeshing(tr, target_edge_length,
+    CGAL::parameters::number_of_iterations(3));
+
+  return EXIT_SUCCESS;
+}

Tetrahedral_remeshing/examples/Tetrahedral_remeshing/tetrahedral_remeshing_example.cpp

@@ -14,7 +14,6 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
 
 typedef CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3<K> Remeshing_triangulation;
 
-
 int main(int argc, char* argv[])
 {
   const double target_edge_length = (argc > 1) ? atof(argv[1]) : 0.1;

Tetrahedral_remeshing/examples/Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp

@@ -79,9 +79,7 @@ int main(int argc, char* argv[])
   boost::unordered_set<std::pair<Vertex_handle, Vertex_handle> > constraints;
 
   CGAL::Tetrahedral_remeshing::generate_input_cube(nbv, t3, constraints);
-
   make_constraints_from_cube_edges(t3, constraints);
-
   CGAL_assertion(t3.is_valid());
 
   CGAL::tetrahedral_isotropic_remeshing(t3, target_edge_length,

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h

@@ -54,12 +54,13 @@ class CollapseTriangulation
   typedef typename C3t3::Cell_handle                          Cell_handle;
   typedef typename C3t3::Vertex_handle                        Vertex_handle;
   typedef typename C3t3::Subdomain_index                      Subdomain_index;
+  typedef typename C3t3::Surface_patch_index                  Surface_patch_index;
   typedef typename C3t3::Triangulation::Point                 Point_3;
   typedef typename C3t3::Triangulation::Geom_traits::Vector_3 Vector_3;
 
 public:
-  CollapseTriangulation(C3t3& c3t3,
+  CollapseTriangulation(const Edge& e,
-                        const Edge& e,
+                        const std::unordered_set<Cell_handle>& cells_to_insert,
                         Collapse_type _collapse_type)
     : collapse_type(_collapse_type)
     , v0_init(e.first->vertex(e.second))
@@ -68,14 +69,6 @@ public:
     typedef std::array<int, 3> Facet; // 3 = id
     typedef std::array<int, 5> Tet_with_ref; // first 4 = id, fifth = reference
 
-    std::unordered_set<Cell_handle> cells_to_insert;
-    c3t3.triangulation().finite_incident_cells(v0_init,
-      std::inserter(cells_to_insert, cells_to_insert.end()));
-    c3t3.triangulation().finite_incident_cells(v1_init,
-      std::inserter(cells_to_insert, cells_to_insert.end()));
-
-    make_cells_set_manifold(c3t3, cells_to_insert);
-
     std::unordered_set<Vertex_handle> vertices_to_insert;
     for (Cell_handle ch : cells_to_insert)
     {
@@ -115,7 +108,7 @@ public:
 
     // finished
     std::vector<Vertex_handle> new_vertices;
-    std::map<Facet, int> border_facets;
+    std::map<Facet, typename C3t3::Surface_patch_index> border_facets;
      if (CGAL::build_triangulation<Tr, false>(triangulation,
                                               points, finite_cells, border_facets,
                                               new_vertices, false/*verbose*/))
@@ -144,58 +137,6 @@ public:
 #endif
   }
 
-  void make_cells_set_manifold(const C3t3& c3t3,
-                               std::unordered_set<Cell_handle>& cells)
-  {
-    typedef Vertex_handle Vh;
-    typedef std::array<Vh, 3> FV;
-    typedef std::pair<Vh, Vh> EV;
-
-    boost::unordered_map<FV, int> facets;
-    for (Cell_handle c : cells)
-    {
-      for (int i = 0; i < 4; ++i)
-      {
-        const FV fvi = make_vertex_array(c->vertex((i + 1) % 4),
-                                         c->vertex((i + 2) % 4),
-                                         c->vertex((i + 3) % 4));
-        typename boost::unordered_map<FV, int>::iterator fit = facets.find(fvi);
-        if(fit == facets.end())
-          facets.insert(std::make_pair(fvi, 1));
-        else
-          fit->second++;
-      }
-    }
-
-    boost::unordered_map<EV, int> edges;
-    for (const std::pair<FV, int>& fvv : facets)
-    {
-      if(fvv.second != 1)
-        continue;
-
-      for (int i = 0; i < 3; ++i)
-      {
-        const EV evi = make_vertex_pair(fvv.first[i], fvv.first[(i + 1) % 3]);
-        typename boost::unordered_map<EV, int>::iterator eit = edges.find(evi);
-        if (eit == edges.end())
-          edges.insert(std::make_pair(evi, 1));
-        else
-          eit->second++;
-      }
-    }
-
-    for (const std::pair<EV, int>& evv : edges)
-    {
-      if (evv.second != 2)
-      {
-        c3t3.triangulation().finite_incident_cells(evv.first.first,
-          std::inserter(cells, cells.begin()));
-        c3t3.triangulation().finite_incident_cells(evv.first.second,
-          std::inserter(cells, cells.begin()));
-      }
-    }
-  }
-
   Result_type collapse()
   {
     if (not_an_edge)
@@ -1038,6 +979,55 @@ typename C3t3::Vertex_handle collapse(typename C3t3::Edge& edge,
   return vh;
 }
 
+template<typename C3t3>
+bool is_cells_set_manifold(const C3t3&,
+    std::unordered_set<typename C3t3::Cell_handle>& cells)
+{
+  typedef typename C3t3::Cell_handle Cell_handle;
+  typedef typename C3t3::Vertex_handle Vh;
+  typedef std::array<Vh, 3> FV;
+  typedef std::pair<Vh, Vh> EV;
+
+  boost::unordered_map<FV, int> facets;
+  for (Cell_handle c : cells)
+  {
+    for (int i = 0; i < 4; ++i)
+    {
+      const FV fvi = make_vertex_array(c->vertex((i + 1) % 4),
+        c->vertex((i + 2) % 4),
+        c->vertex((i + 3) % 4));
+      typename boost::unordered_map<FV, int>::iterator fit = facets.find(fvi);
+      if (fit == facets.end())
+        facets.insert(std::make_pair(fvi, 1));
+      else
+        fit->second++;
+    }
+  }
+
+  boost::unordered_map<EV, int> edges;
+  for (const std::pair<FV, int>& fvv : facets)
+  {
+    if (fvv.second != 1)
+      continue;
+
+    for (int i = 0; i < 3; ++i)
+    {
+      const EV evi = make_vertex_pair(fvv.first[i], fvv.first[(i + 1) % 3]);
+      typename boost::unordered_map<EV, int>::iterator eit = edges.find(evi);
+      if (eit == edges.end())
+        edges.insert(std::make_pair(evi, 1));
+      else
+        eit->second++;
+    }
+  }
+
+  for (const std::pair<EV, int>& evv : edges)
+    if (evv.second != 2)
+      return false;
+
+  return true;
+}
+
 template<typename C3t3, typename CellSelector, typename Visitor>
 typename C3t3::Vertex_handle collapse_edge(typename C3t3::Edge& edge,
     C3t3& c3t3,
@@ -1049,6 +1039,7 @@ typename C3t3::Vertex_handle collapse_edge(typename C3t3::Edge& edge,
   typedef typename C3t3::Triangulation   Tr;
   typedef typename Tr::Point             Point;
   typedef typename Tr::Vertex_handle     Vertex_handle;
+  typedef typename Tr::Cell_handle       Cell_handle;
 
   const Vertex_handle v0 = edge.first->vertex(edge.second);
   const Vertex_handle v1 = edge.first->vertex(edge.third);
@@ -1121,7 +1112,19 @@ typename C3t3::Vertex_handle collapse_edge(typename C3t3::Edge& edge,
     CGAL_assertion(c3t3.triangulation().is_edge(edge.first->vertex(edge.second),
                                                 edge.first->vertex(edge.third), dc, di, dj));
 
-    CollapseTriangulation<C3t3> local_tri(c3t3, edge, collapse_type);
+    Vertex_handle v0_init = edge.first->vertex(edge.second);
+    Vertex_handle v1_init = edge.first->vertex(edge.third);
+
+    std::unordered_set<Cell_handle> cells_to_insert;
+    c3t3.triangulation().finite_incident_cells(v0_init,
+      std::inserter(cells_to_insert, cells_to_insert.end()));
+    c3t3.triangulation().finite_incident_cells(v1_init,
+      std::inserter(cells_to_insert, cells_to_insert.end()));
+
+    if(!is_cells_set_manifold(c3t3, cells_to_insert))
+      return Vertex_handle();
+
+    CollapseTriangulation<C3t3> local_tri(edge, cells_to_insert, collapse_type);
 
     Result_type res = local_tri.collapse();
     if (res == VALID)

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h

@@ -53,7 +53,7 @@ typename C3t3::Vertex_handle split_edge(const typename C3t3::Edge& e,
     (point(v1->point()), point(v2->point()));
 
   //backup subdomain info of incident cells before making changes
-  short dimension = (c3t3.is_in_complex(e)) ? 1 : 3;
+  short dimension = (std::max)((std::max)(1, c3t3.in_dimension(v1)), c3t3.in_dimension(v2));
   boost::unordered_map<Facet, Subdomain_index> cells_info;
   boost::unordered_map<Facet, std::pair<Vertex_handle, Surface_patch_index> > facets_info;
 
@@ -110,11 +110,6 @@ typename C3t3::Vertex_handle split_edge(const typename C3t3::Edge& e,
 
     // surface data for facets of the cells to be split
     const int findex = CGAL::Triangulation_utils_3::next_around_edge(index_v1, index_v2);
-    if (c3t3.is_in_complex(c, findex))
-    {
-      if (dimension == 3)
-        dimension = 2;
-    }
     Surface_patch_index patch = c3t3.surface_patch_index(c, findex);
     Vertex_handle opp_vertex = c->vertex(findex);
     facets_info.insert(std::make_pair(opp_facet1,
@@ -179,6 +174,8 @@ typename C3t3::Vertex_handle split_edge(const typename C3t3::Edge& e,
     // will have its patch tagged from the other side, if needed
   }
 
+  set_index(new_v, c3t3);
+
   return new_v;
 }
 

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h

@@ -30,6 +30,7 @@
 #include <CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h>
 #include <CGAL/Tetrahedral_remeshing/internal/compute_c3t3_statistics.h>
 
+#include <boost/optional.hpp>
 
 namespace CGAL
 {
@@ -103,6 +104,8 @@ class Adaptive_remesher
   typedef typename C3t3::Vertex_handle       Vertex_handle;
   typedef typename C3t3::Subdomain_index     Subdomain_index;
   typedef typename C3t3::Surface_patch_index Surface_patch_index;
+  typedef typename C3t3::Curve_index         Curve_index;
+  typedef typename C3t3::Corner_index        Corner_index;
 
   typedef Tetrahedral_remeshing_smoother<C3t3> Smoother;
 
@@ -141,6 +144,8 @@ public:
 
 #ifdef CGAL_DUMP_REMESHING_STEPS
     CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "00-init");
+    CGAL::Tetrahedral_remeshing::debug::dump_facets_in_complex(m_c3t3,
+      "00-facets_in_complex_after_init.off");
 #endif
   }
 
@@ -167,6 +172,8 @@ public:
 
 #ifdef CGAL_DUMP_REMESHING_STEPS
     CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "00-init");
+    CGAL::Tetrahedral_remeshing::debug::dump_facets_in_complex(m_c3t3,
+      "00-facets_in_complex_after_init.off");
 #endif
   }
 
@@ -184,8 +191,14 @@ public:
     split_long_edges(m_c3t3, emax, m_protect_boundaries,
                      m_cell_selector, m_visitor);
 
+#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
     CGAL_assertion(tr().tds().is_valid(true));
     CGAL_assertion(debug::are_cell_orientations_valid(tr()));
+    CGAL::Tetrahedral_remeshing::debug::dump_facets_in_complex(m_c3t3,
+      "1-facets_in_complex_after_split.off");
+    CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension(
+      m_c3t3.triangulation(), "1-c3t3_vertices_after_split");
+#endif
 #ifdef CGAL_DUMP_REMESHING_STEPS
     CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "1-split");
 #endif
@@ -201,8 +214,12 @@ public:
     collapse_short_edges(m_c3t3, emin, emax, m_protect_boundaries,
                          m_cell_selector, m_visitor);
 
+#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
     CGAL_assertion(tr().tds().is_valid(true));
     CGAL_assertion(debug::are_cell_orientations_valid(tr()));
+    CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension(
+      m_c3t3.triangulation(), "2-c3t3_vertices_after_collapse");
+#endif
 #ifdef CGAL_DUMP_REMESHING_STEPS
     CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "2-collapse");
 #endif
@@ -213,8 +230,12 @@ public:
     flip_edges(m_c3t3, m_protect_boundaries,
                m_cell_selector, m_visitor);
 
+#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
     CGAL_assertion(tr().tds().is_valid(true));
     CGAL_assertion(debug::are_cell_orientations_valid(tr()));
+    CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension(
+      m_c3t3.triangulation(), "3-c3t3_vertices_after_flip");
+#endif
 #ifdef CGAL_DUMP_REMESHING_STEPS
     CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "3-flip");
 #endif
@@ -224,8 +245,12 @@ public:
   {
     m_vertex_smoother.smooth_vertices(m_c3t3, m_protect_boundaries, m_cell_selector);
 
+#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
     CGAL_assertion(tr().tds().is_valid(true));
     CGAL_assertion(debug::are_cell_orientations_valid(tr()));
+    CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension(
+      m_c3t3.triangulation(), "4-c3t3_vertices_after_smooth");
+#endif
 #ifdef CGAL_DUMP_REMESHING_STEPS
     CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "4-smooth");
 #endif
@@ -294,31 +319,33 @@ public:
       }
     }
 
+#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
+    std::cout << "Peelable cells : " << peelable_cells.size() << std::endl;
+#endif
+
     for (auto c_i : peelable_cells)
     {
       Cell_handle c = c_i.first;
       std::array<bool, 4> f_on_surface = c_i.second;
 
-      bool found = false;
+      boost::optional<Surface_patch_index> patch;
-      Surface_patch_index patch;
       for (int i = 0; i < 4; ++i)
       {
         if (f_on_surface[i])
         {
           Surface_patch_index spi = m_c3t3.surface_patch_index(c, i);
-          if (found && patch != spi)
+          if (patch != boost::none && patch != spi)
           {
-            found = false;
+            patch = boost::none;
             break;
           }
           else
           {
-            found = true;
             patch = spi;
           }
         }
       }
-      if(!found)
+      if(patch == boost::none)
         continue;
 
       for (int i = 0; i < 4; ++i)
@@ -326,16 +353,17 @@ public:
         if(f_on_surface[i])
           m_c3t3.remove_from_complex(c, i);
         else
-          m_c3t3.add_to_complex(c, i, patch);
+          m_c3t3.add_to_complex(c, i, patch.get());
       }
 
       m_c3t3.remove_from_complex(c);
       ++nb_slivers_peel;
     }
 
+#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
     CGAL_assertion(tr().tds().is_valid(true));
     CGAL_assertion(debug::are_cell_orientations_valid(tr()));
-
+#endif
 #ifdef CGAL_DUMP_REMESHING_STEPS
     CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "99-postprocess");
 #endif
@@ -359,6 +387,7 @@ private:
                  const FacetIsConstrainedMap& fcmap)
   {
 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
+    debug_c3t3();
     std::size_t nbc = 0;
     std::size_t nbf = 0;
     std::size_t nbe = 0;
@@ -466,7 +495,7 @@ private:
 #endif
 
     //tag vertices
-    unsigned int corner_id = 0;
+    Corner_index corner_id = 0;
     for (Vertex_handle vit : tr().finite_vertex_handles())
     {
       if ( vit->in_dimension() == 0
@@ -478,12 +507,17 @@ private:
         if (vit->in_dimension() == -1 || vit->in_dimension() > 0)
           vit->set_dimension(0);
 
+        vit->set_index(corner_id);
+
 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
         ++nbv;
 #endif
       }
     }
 
+    for (Vertex_handle v : tr().finite_vertex_handles())
+      set_index(v, m_c3t3);
+
 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
     std::cout << "C3t3 ready :" << std::endl;
     std::cout << "\t cells    = " << nbc << std::endl;
@@ -508,6 +542,14 @@ private:
     }
     return true;
   }
+  void debug_c3t3()
+  {
+    for (typename Tr::Facet f : tr().finite_facets())
+    {
+      typename Tr::Facet mf = tr().mirror_facet(f);
+      CGAL_assertion(m_c3t3.is_in_complex(f) == m_c3t3.is_in_complex(mf));
+    }
+  }
 
   template<typename PatchIndex>
   void set_surface_patch_index_to_default(const Subdomain_index&,

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h

@@ -307,6 +307,28 @@ typename C3t3::Surface_patch_index surface_patch_index(const typename C3t3::Vert
   return Surface_patch_index();
 }
 
+template<typename C3t3>
+void set_index(typename C3t3::Vertex_handle v, const C3t3& c3t3)
+{
+  switch (v->in_dimension())
+  {
+  case 3:
+    v->set_index(v->cell()->subdomain_index());
+    break;
+  case 2:
+    v->set_index(surface_patch_index(v, c3t3));
+    break;
+  case 1:
+    v->set_index(typename C3t3::Curve_index(1));
+    break;
+  case 0:
+    v->set_index(boost::get<typename C3t3::Corner_index>(v->index()));
+    break;
+  default:
+    CGAL_assertion(false);
+  }
+}
+
 template<typename C3t3>
 bool is_edge_in_complex(const typename C3t3::Vertex_handle& v0,
                         const typename C3t3::Vertex_handle& v1,

Tetrahedral_remeshing/include/CGAL/tetrahedral_remeshing.h

@@ -252,10 +252,54 @@ void tetrahedral_isotropic_remeshing(
                                  CGAL::parameters::all_default());
 }
 
+/*!
+* \ingroup PkgTetrahedralRemeshingRef
+* converts the triangulation contained in the input to a `Triangulation_3`.
+*
+* This function should be used to generate a valid triangulation
+* for tetrahedral remeshing, when the input triangulation is generated with the
+* tetrahedral mesh generation package.
+*
+* @tparam Tr is the underlying triangulation for `Mesh_complex_3_in_triangulation_3`.
+*            It can be instantiated with any 3D regular triangulation of CGAL provided
+*            that its vertex and cell base classes are models of the concepts
+*            `MeshVertexBase_3` (refined by `RemeshingCellBase_3`)
+*            and `MeshCellBase_3` (refined by `RemeshingVertexBase_3`), respectively.
+* @tparam CornerIndex is the type of the indices for feature corners.
+*            If `c3t3` has been generated using `CGAL::make_mesh_3()`, it must match
+*            the `Corner_index` type of the model of the `MeshDomainWithFeatures_3` concept used for mesh generation.
+* @tparam CurveIndex is the type of the indices for feature curves.
+*            If `c3t3` has been generated using `CGAL::make_mesh_3()`, it must match
+*            the `Curve_index` type of the model of the `MeshDomainWithFeatures_3` concept used for mesh generation.
+*
+* @param c3t3 the complex containing the triangulation to be remeshed.
+*/
+
+template<typename Tr,
+         typename CornerIndex,
+         typename CurveIndex>
+CGAL::Triangulation_3<typename Tr::Geom_traits,
+                      typename Tr::Triangulation_data_structure>
+convert_to_triangulation_3(
+  CGAL::Mesh_complex_3_in_triangulation_3<Tr, CornerIndex, CurveIndex> c3t3)
+{
+  using GT   = typename Tr::Geom_traits;
+  using TDS  = typename Tr::Triangulation_data_structure;
+
+  CGAL::Triangulation_3<GT, TDS> tr;
+  tr.swap(c3t3.triangulation());
+  return tr;
+}
+
 ///////////////////////////////////////////////////
 /////// MESH_COMPLEX_3_IN_TRIANGULATION_3 /////////
 ///////////////////////////////////////////////////
 
+///////
+////// Warning with using this version :
+////// the triangulation after remeshing is not regular anymore
+////// the empty-sphere property is not respected
+///////
 template<typename Tr,
          typename CornerIndex, typename CurveIndex,
          typename NamedParameters>

Tetrahedral_remeshing/test/Tetrahedral_remeshing/CMakeLists.txt

@@ -6,7 +6,7 @@ cmake_minimum_required(VERSION 3.1...3.14)
 project( Tetrahedral_remeshing_Tests )
 
 # CGAL and its components
-find_package( CGAL REQUIRED )
+find_package(CGAL COMPONENTS ImageIO)
 if ( NOT CGAL_FOUND )
   message(STATUS "This project requires the CGAL library, and will not be compiled.")
   return()
@@ -19,10 +19,26 @@ if ( NOT Boost_FOUND )
   return()
 endif()
 
-
 # Creating entries for all C++ files with "main" routine
 # ##########################################################
 create_single_source_cgal_program( "test_tetrahedral_remeshing.cpp" )
 create_single_source_cgal_program( "test_tetrahedral_remeshing_with_features.cpp")
 create_single_source_cgal_program( "test_tetrahedral_remeshing_of_one_subdomain.cpp")
 create_single_source_cgal_program( "test_tetrahedral_remeshing_io.cpp")
+create_single_source_cgal_program( "test_tetrahedral_remeshing_from_mesh_file.cpp")
+
+# Tests using Mesh_3 require Eigen
+find_package(Eigen3 3.1.0 REQUIRED) #(3.1.0 or greater)
+include(CGAL_Eigen_support)
+if (NOT TARGET CGAL::Eigen_support)
+  message(STATUS "This project requires the Eigen library, and will not be compiled.")
+  return()
+endif()
+
+create_single_source_cgal_program("test_mesh_and_remesh_polyhedral_domain_with_features.cpp")
+target_link_libraries(test_mesh_and_remesh_polyhedral_domain_with_features PUBLIC CGAL::Eigen_support)
+
+if( CGAL_ImageIO_USE_ZLIB )
+  create_single_source_cgal_program( "test_mesh_and_remesh_image.cpp")
+  target_link_libraries(test_mesh_and_remesh_image PUBLIC CGAL::Eigen_support)
+endif()

Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_mesh_and_remesh_image.cpp

@@ -0,0 +1,72 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+
+#include <CGAL/Labeled_mesh_domain_3.h>
+#include <CGAL/make_mesh_3.h>
+#include <CGAL/Image_3.h>
+
+#include <CGAL/tetrahedral_remeshing.h>
+
+// Domain
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+typedef CGAL::Image_3 Image;
+typedef CGAL::Labeled_mesh_domain_3<K> Mesh_domain;
+
+// Triangulation
+typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
+typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
+
+// Mesh Criteria
+typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+typedef Mesh_criteria::Facet_criteria    Facet_criteria;
+typedef Mesh_criteria::Cell_criteria     Cell_criteria;
+
+typedef CGAL::Triangulation_3<typename Tr::Geom_traits,
+                              typename Tr::Triangulation_data_structure> T3;
+
+using namespace CGAL::parameters;
+
+int main()
+{
+  const char* filename = "data/liver.inr.gz";
+
+  CGAL::Image_3 image;
+  if (!image.read(filename)) {
+    std::cerr << "Error: Cannot read file " << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+  Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image, 1e-9);
+
+  // Mesh criteria
+  Facet_criteria facet_criteria(25, 20, 2); // angle, size, approximation
+  Cell_criteria cell_criteria(3, 20); // radius-edge ratio, size
+  Mesh_criteria criteria(facet_criteria, cell_criteria);
+
+  // Meshing
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_exude(), no_perturb());
+
+  std::cout << "Meshing done." << std::endl;
+
+  //Remeshing : extract triangulation
+  T3 t3 = CGAL::convert_to_triangulation_3(c3t3);
+
+  //Remeshing : coarsen
+  double target_edge_length = 15.;
+  CGAL::tetrahedral_isotropic_remeshing(t3, target_edge_length,
+    CGAL::parameters::number_of_iterations(2));
+
+  std::cout << "Remeshing 1 done." << std::endl;
+
+  //Remeshing : refine
+  target_edge_length = 20.;
+  CGAL::tetrahedral_isotropic_remeshing(t3, target_edge_length,
+    CGAL::parameters::number_of_iterations(2).remesh_boundaries(false));
+
+  std::cout << "Remeshing 2 done." << std::endl;
+
+  return 0;
+}

Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_mesh_and_remesh_polyhedral_domain_with_features.cpp

@@ -0,0 +1,81 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
+#include <CGAL/make_mesh_3.h>
+
+#include <CGAL/tetrahedral_remeshing.h>
+
+// Domain
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Surface_mesh<K::Point_3> Polyhedron;
+typedef CGAL::Polyhedral_mesh_domain_with_features_3<K, Polyhedron> Mesh_domain;
+
+#ifdef CGAL_CONCURRENT_MESH_3
+typedef CGAL::Parallel_tag Concurrency_tag;
+#else
+typedef CGAL::Sequential_tag Concurrency_tag;
+#endif
+
+// Triangulation
+typedef CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type Tr;
+
+typedef CGAL::Mesh_complex_3_in_triangulation_3<
+  Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index> C3t3;
+
+// Criteria
+typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+
+// To avoid verbose function and named parameters call
+using namespace CGAL::parameters;
+
+int main(int argc, char* argv[])
+{
+  const char* fname = (argc > 1) ? argv[1] : "data/fandisk.off";
+  std::ifstream input(fname);
+  Polyhedron polyhedron;
+  input >> polyhedron;
+  if (input.fail()) {
+    std::cerr << "Error: Cannot read file " << fname << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  if (!CGAL::is_triangle_mesh(polyhedron)) {
+    std::cerr << "Input geometry is not triangulated." << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  // Create domain
+  Mesh_domain domain(polyhedron);
+
+  // Get sharp features
+  domain.detect_features();
+
+  // Mesh criteria
+  Mesh_criteria criteria(edge_size = 0.025,
+    facet_angle = 25, facet_size = 0.05, facet_distance = 0.005,
+    cell_radius_edge_ratio = 3, cell_size = 0.05);
+
+  // Mesh generation
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
+
+  auto t3 = CGAL::convert_to_triangulation_3(std::move(c3t3));
+
+  double target_edge_length = 0.1;//coarsen the mesh
+  CGAL::tetrahedral_isotropic_remeshing(t3, target_edge_length,
+    CGAL::parameters::number_of_iterations(3).remesh_boundaries(false));
+
+  assert(t3.is_valid());
+
+  target_edge_length = 0.04;//re-refine the mesh
+  CGAL::tetrahedral_isotropic_remeshing(t3, target_edge_length,
+    CGAL::parameters::number_of_iterations(1));//remesh_boundaries(true)
+
+  assert(t3.is_valid());
+
+  return EXIT_SUCCESS;
+}

Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_tetrahedral_remeshing.cpp

@@ -1,6 +1,6 @@
 #define CGAL_TETRAHEDRAL_REMESHING_VERBOSE
-#define CGAL_DUMP_REMESHING_STEPS
+//#define CGAL_DUMP_REMESHING_STEPS
-#define CGAL_TETRAHEDRAL_REMESHING_DEBUG
+//#define CGAL_TETRAHEDRAL_REMESHING_DEBUG
 //#define CGAL_TETRAHEDRAL_REMESHING_GENERATE_INPUT_FILES
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_tetrahedral_remeshing_from_mesh_file.cpp

@@ -0,0 +1,32 @@
+//#define CGAL_TETRAHEDRAL_REMESHING_DEBUG
+//#define CGAL_TETRAHEDRAL_REMESHING_VERBOSE
+//#define CGAL_DUMP_REMESHING_STEPS
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Tetrahedral_remeshing/Remeshing_triangulation_3.h>
+#include <CGAL/tetrahedral_remeshing.h>
+
+#include <iostream>
+#include <fstream>
+
+#include <CGAL/Mesh_3/tet_soup_to_c3t3.h>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+typedef CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3<K> Remeshing_triangulation;
+
+int main(int argc, char* argv[])
+{
+  const double target_edge_length = (argc > 1) ? atof(argv[1]) : 0.1;
+
+  Remeshing_triangulation tr;
+
+  std::ifstream in("data/sphere.mesh");
+  if (CGAL::build_triangulation_from_file<Remeshing_triangulation, true>(in, tr))
+    std::cout << "build triangulation ok" << std::endl;
+
+  CGAL::tetrahedral_isotropic_remeshing(tr, target_edge_length);
+
+  return EXIT_SUCCESS;
+}

Triangulation_3/include/CGAL/Triangulation_3.h

@@ -567,7 +567,7 @@ protected:
 
 public:
   template<typename P> // Point or Point_3
-  typename boost::result_of<const Construct_point_3(const P&)>::type
+  typename boost::result_of<Construct_point_3(P)>::type
   construct_point(const P& p) const
   {
     return geom_traits().construct_point_3_object()(p);