diff --git a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_approximation_test.cpp b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_approximation_test.cpp
index 8a01d8b254e..c657eaaf00d 100644
--- a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_approximation_test.cpp
+++ b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_approximation_test.cpp
@@ -2,29 +2,46 @@
 #include <fstream>
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Polyhedron_3.h>
-#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
-#include <CGAL/IO/Polyhedron_iostream.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/remesh.h>
+
 #include <CGAL/Surface_mesh_approximation/approximate_triangle_mesh.h>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
+typedef CGAL::Surface_mesh<Kernel::Point_3> Mesh;
+typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
 
 /**
  * This file tests the free function CGAL::Surface_mesh_approximation::approximate_triangle_mesh.
  */
 int main()
 {
-  Polyhedron mesh;
-  std::ifstream input("./data/cube_meshed.off");
-  if (!input || !(input >> mesh) || mesh.empty()) {
-    std::cerr << "Invalid off file." << std::endl;
+  Mesh mesh;
+  std::ifstream input("./data/cube.off");
+  if (!input || !(input >> mesh) || !CGAL::is_triangle_mesh(mesh)) {
+    std::cerr << "Invalid input file." << std::endl;
     return EXIT_FAILURE;
   }
 
-  Polyhedron out_mesh;
-  std::map<Polyhedron::Facet_handle, std::size_t> fidxmap;
-  boost::associative_property_map<std::map<Polyhedron::Facet_handle, std::size_t> > fpxmap(fidxmap);
+  const double target_edge_length = 0.05;
+  const unsigned int nb_iter = 3;
+
+  std::cout << "Start remeshing. "
+    << " (" << num_faces(mesh) << " faces)..." << std::endl;
+  PMP::isotropic_remeshing(
+    faces(mesh),
+    target_edge_length,
+    mesh,
+    PMP::parameters::number_of_iterations(nb_iter));
+  std::cout << "Remeshing done. "
+    << " (" << num_faces(mesh) << " faces)..." << std::endl;
+
+  Mesh out_mesh;
+  std::map<face_descriptor, std::size_t> fidxmap;
+  boost::associative_property_map<std::map<face_descriptor, std::size_t> > fpxmap(fidxmap);
   std::vector<Kernel::Vector_3> proxies;
   std::vector<Kernel::Point_3> points;
   std::vector<CGAL::cpp11::array<std::size_t, 3> > triangles;
diff --git a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_class_interface_test.cpp b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_class_interface_test.cpp
index 96bfedc24e1..8a158ac1d1b 100644
--- a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_class_interface_test.cpp
+++ b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_class_interface_test.cpp
@@ -2,52 +2,65 @@
 #include <fstream>
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Polyhedron_3.h>
-#include <CGAL/IO/Polyhedron_iostream.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/remesh.h>
 
-#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
 #include <CGAL/Variational_shape_approximation.h>
 #include <CGAL/Surface_mesh_approximation/L2_metric_plane_proxy.h>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-typedef Kernel::FT FT;
-typedef Kernel::Point_3 Point;
+typedef CGAL::Surface_mesh<Kernel::Point_3> Mesh;
+typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
+typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
 
-typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
-typedef Polyhedron::Facet_handle Facet_handle;
-typedef boost::associative_property_map<std::map<Facet_handle, std::size_t> > Face_proxy_map;
-typedef boost::property_map<Polyhedron, boost::vertex_point_t>::type Vertex_point_map;
+typedef boost::associative_property_map<std::map<face_descriptor, std::size_t> > Face_proxy_map;
+typedef boost::property_map<Mesh, boost::vertex_point_t>::type Vertex_point_map;
 
-typedef CGAL::Surface_mesh_approximation::L2_metric_plane_proxy<Polyhedron> L2_metric_plane_proxy;
-typedef CGAL::Variational_shape_approximation<Polyhedron, Vertex_point_map, L2_metric_plane_proxy> L2_approx;
+typedef CGAL::Surface_mesh_approximation::L2_metric_plane_proxy<Mesh> L2_metric_plane_proxy;
+typedef CGAL::Variational_shape_approximation<Mesh, Vertex_point_map, L2_metric_plane_proxy> L2_approx;
 typedef L2_approx::Proxy Plane_proxy;
 
+namespace PMP = CGAL::Polygon_mesh_processing;
+
 /**
  * This file tests the main class API and the L2 metric.
  * It should cover all the APIs.
  */
 int main()
 {
-  Polyhedron mesh;
-  std::ifstream input("./data/sphere_iso.off");
-  if (!input || !(input >> mesh) || mesh.empty()) {
-    std::cerr << "Invalid off file." << std::endl;
+  Mesh mesh;
+  std::ifstream input("./data/sphere.off");
+  if (!input || !(input >> mesh) || !CGAL::is_triangle_mesh(mesh)) {
+    std::cerr << "Invalid input file." << std::endl;
     return EXIT_FAILURE;
   }
 
+  const double target_edge_length = 0.05;
+  const unsigned int nb_iter = 3;
+
+  std::cout << "Start remeshing. "
+    << " (" << num_faces(mesh) << " faces)..." << std::endl;
+  PMP::isotropic_remeshing(
+    faces(mesh),
+    target_edge_length,
+    mesh,
+    PMP::parameters::number_of_iterations(nb_iter));
+  std::cout << "Remeshing done. "
+    << " (" << num_faces(mesh) << " faces)..." << std::endl;
+
   // face area map
-  std::map<Facet_handle, std::size_t> face_index;
-  for (Polyhedron::Facet_iterator fitr = mesh.facets_begin();
-    fitr != mesh.facets_end(); ++fitr)
-    face_index.insert(std::pair<Facet_handle, std::size_t>(fitr, 0));
+  std::map<face_descriptor, std::size_t> face_index;
+  BOOST_FOREACH(face_descriptor f, faces(mesh))
+    face_index.insert(std::pair<face_descriptor, std::size_t>(f, 0));
   Face_proxy_map proxy_pmap(face_index);
 
   // create L2_approx L2 metric approximation algorithm instance
   std::cout << "setup algorithm instance" << std::endl;
   L2_metric_plane_proxy error_metric(mesh,
-    get(boost::vertex_point, const_cast<Polyhedron &>(mesh)));
+    get(boost::vertex_point, const_cast<Mesh &>(mesh)));
   L2_approx approx(mesh,
-    get(boost::vertex_point, const_cast<Polyhedron &>(mesh)),
+    get(boost::vertex_point, const_cast<Mesh &>(mesh)),
     error_metric);
 
   // random seeding and run
@@ -87,7 +100,7 @@ int main()
   if (approx.number_of_proxies() != 17)
     return EXIT_FAILURE;
 
-  // extract the approximation polyhedron
+  // extract the approximation Mesh
   std::cout << "meshing" << std::endl;
   if (approx.extract_mesh(CGAL::parameters::subdivision_ratio(1.0)))
     std::cout << "manifold." << std::endl;
@@ -99,17 +112,17 @@ int main()
   approx.proxy_map(proxy_pmap);
 
   for (std::size_t i = 0; i < approx.number_of_proxies(); ++i) {
-    std::list<Facet_handle> patch;
+    std::list<face_descriptor> patch;
     approx.proxy_region(i, std::back_inserter(patch));
   }
 
   std::vector<Plane_proxy> proxies;
   approx.proxies(std::back_inserter(proxies));
 
-  std::vector<Point> anchor_pos;
+  std::vector<Kernel::Point_3> anchor_pos;
   approx.anchor_points(std::back_inserter(anchor_pos));
 
-  std::vector<Polyhedron::Vertex_handle> anchor_vtx;
+  std::vector<vertex_descriptor> anchor_vtx;
   approx.anchor_vertices(std::back_inserter(anchor_vtx));
 
   std::vector<CGAL::cpp11::array<std::size_t, 3> > tris;
@@ -118,7 +131,7 @@ int main()
   std::vector<std::vector<std::size_t> > boundary;
   approx.indexed_boundary_polygons(std::back_inserter(boundary));
 
-  const FT drop(0.001);
+  const Kernel::FT drop(0.001);
   const std::size_t iterations = 5;
   std::cout << "re-initialize and hierarchical seeding" << std::endl;
   approx.initialize_seeds(CGAL::parameters::seeding_method(CGAL::Surface_mesh_approximation::HIERARCHICAL)
@@ -134,7 +147,7 @@ int main()
   approx.run(10);
   std::cout << "#proxies " << approx.number_of_proxies() << std::endl;
 
-  // extract the approximation polyhedron
+  // extract the approximation Mesh
   std::cout << "meshing" << std::endl;
   if (approx.extract_mesh(CGAL::parameters::subdivision_ratio(1.0)))
     std::cout << "manifold." << std::endl;
diff --git a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_correctness_test.cpp b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_correctness_test.cpp
index f42dfa7fcd1..a8a961994cd 100644
--- a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_correctness_test.cpp
+++ b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_correctness_test.cpp
@@ -2,40 +2,40 @@
 #include <fstream>
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Polyhedron_3.h>
-#include <CGAL/IO/Polyhedron_iostream.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/remesh.h>
 
-#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
 #include <CGAL/Variational_shape_approximation.h>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-typedef Kernel::FT FT;
-typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
-typedef boost::property_map<Polyhedron, boost::vertex_point_t>::type Vertex_point_map;
+typedef CGAL::Surface_mesh<Kernel::Point_3> Mesh;
+typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
 
-typedef CGAL::Variational_shape_approximation<Polyhedron, Vertex_point_map> L21_approx;
+typedef boost::property_map<Mesh, boost::vertex_point_t>::type Vertex_point_map;
+typedef CGAL::Variational_shape_approximation<Mesh, Vertex_point_map> L21_approx;
 typedef L21_approx::Error_metric L21_metric;
 
 bool test_shape(const char *file_name, const std::size_t target_num_proxies)
 {
-  Polyhedron mesh;
+  Mesh mesh;
   std::ifstream input(file_name);
-  if (!input || !(input >> mesh) || mesh.empty()) {
-    std::cout << "Invalid off file." << std::endl;
+  if (!input || !(input >> mesh) || !CGAL::is_triangle_mesh(mesh)) {
+    std::cout << "Invalid input file." << std::endl;
     return false;
   }
 
   std::cout << "Testing \"" << file_name << '\"' << std::endl;
   // algorithm instance
   L21_metric error_metric(mesh,
-    get(boost::vertex_point, const_cast<Polyhedron &>(mesh)));
+    get(boost::vertex_point, const_cast<Mesh &>(mesh)));
   L21_approx approx(mesh,
-    get(boost::vertex_point, const_cast<Polyhedron &>(mesh)),
+    get(boost::vertex_point, const_cast<Mesh &>(mesh)),
     error_metric);
 
   // approximation, seeding from error, drop to the target error incrementally
   // should reach targeted number of proxies gradually
-  const FT drop(1e-8);
+  const Kernel::FT drop(1e-8);
   const std::size_t num_iterations = 20;
   const std::size_t inner_iterations = 10;
   approx.initialize_seeds(CGAL::parameters::seeding_method(CGAL::Surface_mesh_approximation::INCREMENTAL)
@@ -45,8 +45,7 @@ bool test_shape(const char *file_name, const std::size_t target_num_proxies)
 
   // eliminate redundant area (local minima) by merging
   boost::optional<std::pair<std::size_t, std::size_t> > best_pair = boost::none;
-  while ( ( best_pair = approx.find_best_merge(true) ) != boost::none )
-  {
+  while ((best_pair = approx.find_best_merge(true)) != boost::none) {
     approx.merge(best_pair->first, best_pair->second);
     approx.run(num_iterations);
   }
@@ -73,14 +72,14 @@ bool test_shape(const char *file_name, const std::size_t target_num_proxies)
 int main()
 {
   std::cout << "Correctness test." << std::endl;
-  if (!test_shape("./data/cube_meshed.off", 6))
+  if (!test_shape("./data/cube.off", 6))
     return EXIT_FAILURE;
 
-  if (!test_shape("./data/cube_meshed_open.off", 5))
+  if (!test_shape("./data/cube-ouvert.off", 5))
     return EXIT_FAILURE;
 
   std::cout << "Surface with disconnected components test." << std::endl;
-  if (!test_shape("./data/cubes_disconnected.off", 11))
+  if (!test_shape("./data/cubes-merged.off", 11))
     return EXIT_FAILURE;
 
   return EXIT_SUCCESS;
diff --git a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_error_decrease_test.cpp b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_error_decrease_test.cpp
index 6717487aea3..8f62d3815de 100644
--- a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_error_decrease_test.cpp
+++ b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_error_decrease_test.cpp
@@ -2,20 +2,23 @@
 #include <fstream>
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Polyhedron_3.h>
-#include <CGAL/IO/Polyhedron_iostream.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/remesh.h>
 
-#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
 #include <CGAL/Variational_shape_approximation.h>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
 typedef Kernel::FT FT;
-typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
-typedef boost::property_map<Polyhedron, boost::vertex_point_t>::type Vertex_point_map;
 
-typedef CGAL::Variational_shape_approximation<Polyhedron, Vertex_point_map> L21_approx;
+typedef CGAL::Surface_mesh<Kernel::Point_3> Mesh;
+typedef boost::property_map<Mesh, boost::vertex_point_t>::type Vertex_point_map;
+
+typedef CGAL::Variational_shape_approximation<Mesh, Vertex_point_map> L21_approx;
 typedef L21_approx::Error_metric L21_metric;
 
+namespace PMP = CGAL::Polygon_mesh_processing;
+
 bool check_strict_ordering(const std::vector<FT> &error)
 {
   if (error.empty()) {
@@ -35,18 +38,31 @@ bool check_strict_ordering(const std::vector<FT> &error)
  */
 int main()
 {
-  Polyhedron mesh;
-  std::ifstream input("./data/sphere_iso.off");
-  if (!input || !(input >> mesh) || mesh.empty()) {
-    std::cerr << "Invalid off file." << std::endl;
+  Mesh mesh;
+  std::ifstream input("./data/sphere.off");
+  if (!input || !(input >> mesh) || !CGAL::is_triangle_mesh(mesh)) {
+    std::cerr << "Invalid input file." << std::endl;
     return EXIT_FAILURE;
   }
 
+  const double target_edge_length = 0.05;
+  const unsigned int nb_iter = 3;
+
+  std::cout << "Start remeshing. "
+    << " (" << num_faces(mesh) << " faces)..." << std::endl;
+  PMP::isotropic_remeshing(
+    faces(mesh),
+    target_edge_length,
+    mesh,
+    PMP::parameters::number_of_iterations(nb_iter));
+  std::cout << "Remeshing done. "
+    << " (" << num_faces(mesh) << " faces)..." << std::endl;
+
   // algorithm instance
   L21_metric error_metric(mesh,
-    get(boost::vertex_point, const_cast<Polyhedron &>(mesh)));
+    get(boost::vertex_point, const_cast<Mesh &>(mesh)));
   L21_approx approx(mesh,
-    get(boost::vertex_point, const_cast<Polyhedron &>(mesh)),
+    get(boost::vertex_point, const_cast<Mesh &>(mesh)),
     error_metric);
 
   approx.initialize_seeds(CGAL::parameters::seeding_method(CGAL::Surface_mesh_approximation::RANDOM)
diff --git a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_kernel_test.cpp b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_kernel_test.cpp
index f0adc0a3b07..55efc29c749 100644
--- a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_kernel_test.cpp
+++ b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_kernel_test.cpp
@@ -10,33 +10,93 @@
 
 #include <CGAL/Surface_mesh.h>
 
+#include <CGAL/Polygon_mesh_processing/remesh.h>
+
 #include <CGAL/Surface_mesh_approximation/approximate_triangle_mesh.h>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic;
 typedef CGAL::Simple_cartesian<double> Sckernel;
 
-template <typename K, typename TM>
-int test() {
-  TM tm;
-  std::ifstream input("./data/cube_meshed.off");
-  if (!input || !(input >> tm) || num_vertices(tm) == 0) {
-    std::cerr << "Invalid off file." << std::endl;
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+template <typename TM>
+int load_and_remesh_sm(TM &mesh) {
+  std::ifstream input("./data/cube.off");
+  if (!input || !(input >> mesh) || !CGAL::is_triangle_mesh(mesh)) {
+    std::cerr << "Invalid input file." << std::endl;
     return EXIT_FAILURE;
   }
 
-  typedef CGAL::Polyhedron_3<K> Polyhedron;
-  Polyhedron out_mesh;
+  const double target_edge_length = 0.05;
+  const unsigned int nb_iter = 3;
+
+  std::cout << "Start remeshing. "
+    << " (" << num_faces(mesh) << " faces)..." << std::endl;
+  PMP::isotropic_remeshing(
+    faces(mesh),
+    target_edge_length,
+    mesh,
+    PMP::parameters::number_of_iterations(nb_iter));
+  std::cout << "Remeshing done. "
+    << " (" << num_faces(mesh) << " faces)..." << std::endl;
+
+  return EXIT_SUCCESS;
+}
+
+template <typename TM>
+int load_and_remesh_poly(TM &mesh) {
+  std::ifstream input("./data/cube.off");
+  if (!input || !(input >> mesh) || !CGAL::is_triangle_mesh(mesh)) {
+    std::cerr << "Invalid input file." << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  const double target_edge_length = 0.05;
+  const unsigned int nb_iter = 3;
+
+  std::cout << "Start remeshing. "
+    << " (" << num_faces(mesh) << " faces)..." << std::endl;
+  PMP::isotropic_remeshing(
+    faces(mesh),
+    target_edge_length,
+    mesh,
+    PMP::parameters::number_of_iterations(nb_iter).
+    face_index_map(get(boost::face_external_index, mesh)));
+  std::cout << "Remeshing done. "
+    << " (" << num_faces(mesh) << " faces)..." << std::endl;
+
+  return EXIT_SUCCESS;
+}
+
+template <typename K, typename TM>
+void run_approximation(const TM &mesh) {
   std::vector<typename K::Point_3> points;
   std::vector<CGAL::cpp11::array<std::size_t, 3> > triangles;
 
-  CGAL::Surface_mesh_approximation::approximate_triangle_mesh(tm,
+  CGAL::Surface_mesh_approximation::approximate_triangle_mesh(mesh,
     CGAL::parameters::max_number_of_proxies(6).
       number_of_iterations(30).
       number_of_relaxations(5).
       subdivision_ratio(0.5).
       anchors(std::back_inserter(points)).
       triangles(std::back_inserter(triangles)));
+}
 
+template <typename K, typename TM>
+int test_sm() {
+  TM mesh;
+  if (load_and_remesh_sm(mesh) == EXIT_FAILURE)
+    return EXIT_FAILURE;
+  run_approximation<K, TM>(mesh);
+  return EXIT_SUCCESS;
+}
+
+template <typename K, typename TM>
+int test_poly() {
+  TM mesh;
+  if (load_and_remesh_poly(mesh) == EXIT_FAILURE)
+    return EXIT_FAILURE;
+  run_approximation<K, TM>(mesh);
   return EXIT_SUCCESS;
 }
 
@@ -45,16 +105,16 @@ int test() {
  */
 int main()
 {
-  if (test<Epic, CGAL::Polyhedron_3<Epic> >() == EXIT_FAILURE)
+  if (test_poly<Epic, CGAL::Polyhedron_3<Epic> >() == EXIT_FAILURE)
     return EXIT_FAILURE;
 
-  if (test<Epic, CGAL::Surface_mesh<Epic::Point_3> >() == EXIT_FAILURE)
+  if (test_sm<Epic, CGAL::Surface_mesh<Epic::Point_3> >() == EXIT_FAILURE)
     return EXIT_FAILURE;
 
-  if (test<Sckernel, CGAL::Polyhedron_3<Sckernel> >() == EXIT_FAILURE)
+  if (test_poly<Sckernel, CGAL::Polyhedron_3<Sckernel> >() == EXIT_FAILURE)
     return EXIT_FAILURE;
 
-  if (test<Sckernel, CGAL::Surface_mesh<Sckernel::Point_3> >() == EXIT_FAILURE)
+  if (test_sm<Sckernel, CGAL::Surface_mesh<Sckernel::Point_3> >() == EXIT_FAILURE)
     return EXIT_FAILURE;
 
   return EXIT_SUCCESS;
diff --git a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_meshing_manifold_test.cpp b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_meshing_manifold_test.cpp
index 7ee22e164ca..e79ef5dd449 100644
--- a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_meshing_manifold_test.cpp
+++ b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_meshing_manifold_test.cpp
@@ -2,35 +2,50 @@
 #include <fstream>
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Polyhedron_3.h>
-#include <CGAL/IO/Polyhedron_iostream.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/remesh.h>
 
-#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
 #include <CGAL/Variational_shape_approximation.h>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
 typedef Kernel::FT FT;
-typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
-typedef boost::property_map<Polyhedron, boost::vertex_point_t>::type Vertex_point_map;
+typedef CGAL::Surface_mesh<Kernel::Point_3> Mesh;
+typedef boost::property_map<Mesh, boost::vertex_point_t>::type Vertex_point_map;
 
-typedef CGAL::Variational_shape_approximation<Polyhedron, Vertex_point_map> L21_approx;
+typedef CGAL::Variational_shape_approximation<Mesh, Vertex_point_map> L21_approx;
 typedef L21_approx::Error_metric L21_metric;
 
+namespace PMP = CGAL::Polygon_mesh_processing;
+
 bool test_manifold(const char *file_name, const FT drop = FT(1e-8))
 {
-  Polyhedron mesh;
+  Mesh mesh;
   std::ifstream input(file_name);
-  if (!input || !(input >> mesh) || mesh.empty()) {
-    std::cout << "Invalid off file." << std::endl;
+  if (!input || !(input >> mesh) || !CGAL::is_triangle_mesh(mesh)) {
+    std::cout << "Invalid input file." << std::endl;
     return false;
   }
 
+  const double target_edge_length = 0.05;
+  const unsigned int nb_iter = 3;
+
+  std::cout << "Start remeshing. "
+    << " (" << num_faces(mesh) << " faces)..." << std::endl;
+  PMP::isotropic_remeshing(
+    faces(mesh),
+    target_edge_length,
+    mesh,
+    PMP::parameters::number_of_iterations(nb_iter));
+  std::cout << "Remeshing done. "
+    << " (" << num_faces(mesh) << " faces)..." << std::endl;
+
   std::cout << "Testing \"" << file_name << '\"' << std::endl;
   // algorithm instance
   L21_metric error_metric(mesh,
-    get(boost::vertex_point, const_cast<Polyhedron &>(mesh)));
+    get(boost::vertex_point, const_cast<Mesh &>(mesh)));
   L21_approx approx(mesh,
-    get(boost::vertex_point, const_cast<Polyhedron &>(mesh)),
+    get(boost::vertex_point, const_cast<Mesh &>(mesh)),
     error_metric);
 
   // approximation, seeding from error, drop to the target error incrementally
@@ -59,13 +74,13 @@ bool test_manifold(const char *file_name, const FT drop = FT(1e-8))
 int main()
 {
   std::cout << "Meshing manifold test." << std::endl;
-  if (!test_manifold("./data/cube_meshed.off"))
+  if (!test_manifold("./data/cube.off"))
     return EXIT_FAILURE;
 
-  if (!test_manifold("./data/cube_meshed_open.off"))
+  if (!test_manifold("./data/cube-ouvert.off"))
     return EXIT_FAILURE;
 
-  if (!test_manifold("./data/sphere_iso.off", FT(1e-2)))
+  if (!test_manifold("./data/sphere.off", FT(1e-2)))
     return EXIT_FAILURE;
 
   return EXIT_SUCCESS;
diff --git a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_metric_test.cpp b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_metric_test.cpp
index b70701f4aa0..aa6ada16a03 100644
--- a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_metric_test.cpp
+++ b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_metric_test.cpp
@@ -2,30 +2,31 @@
 #include <fstream>
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Polyhedron_3.h>
-#include <CGAL/IO/Polyhedron_iostream.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/remesh.h>
 
-#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
-#include <CGAL/property_map.h>
 #include <CGAL/Variational_shape_approximation.h>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
 typedef Kernel::FT FT;
-typedef Kernel::Vector_3 Vector;
-typedef Kernel::Point_3 Point;
+typedef Kernel::Vector_3 Vector_3;
+typedef Kernel::Point_3 Point_3;
 
-typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
-typedef Polyhedron::Facet_handle Facet_handle;
-typedef Polyhedron::Halfedge_handle Halfedge_handle;
-typedef Polyhedron::Facet_iterator Facet_iterator;
-typedef boost::associative_property_map<std::map<Facet_handle, FT> > Face_area_map;
-typedef boost::associative_property_map<std::map<Facet_handle, Point> > Face_center_map;
-typedef boost::property_map<Polyhedron, boost::vertex_point_t>::type Vertex_point_map;
+typedef CGAL::Surface_mesh<Point_3> Mesh;
+typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
+typedef boost::graph_traits<Mesh>::halfedge_descriptor halfedge_descriptor;
+
+typedef Mesh::Property_map<face_descriptor, FT> Face_area_map;
+typedef Mesh::Property_map<face_descriptor, Point_3> Face_center_map;
+typedef boost::property_map<Mesh, boost::vertex_point_t>::type Vertex_point_map;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
 
 // user defined point-wise compact metric
 struct Compact_metric_point_proxy {
 // use point as proxy
-  typedef Point Proxy;
+  typedef Point_3 Proxy;
 
   // we keep a precomputed property map to speed up computations
   Compact_metric_point_proxy(const Face_center_map &_center_pmap, const Face_area_map &_area_pmap)
@@ -34,8 +35,7 @@ struct Compact_metric_point_proxy {
   // compute and return error from a face to a proxy,
   // defined as the Euclidean distance between
   // the face center of mass and proxy point.
-  FT compute_error(const Facet_handle &f, const Polyhedron &tm, const Proxy &px) const {
-    (void)(tm);
+  FT compute_error(const face_descriptor &f, const Mesh &, const Proxy &px) const {
     return FT(std::sqrt(CGAL::to_double(
       CGAL::squared_distance(center_pmap[f], px))));
   }
@@ -43,12 +43,11 @@ struct Compact_metric_point_proxy {
   // template functor to compute a best-fit
   // proxy from a range of faces
   template <typename FaceRange>
-  Proxy fit_proxy(const FaceRange &faces, const Polyhedron &tm) const {
-    (void)(tm);
+  Proxy fit_proxy(const FaceRange &faces, const Mesh &) const {
     // fitting center
-    Vector center = CGAL::NULL_VECTOR;
+    Vector_3 center = CGAL::NULL_VECTOR;
     FT sum_areas = FT(0.0);
-    BOOST_FOREACH(const Facet_handle &f, faces) {
+    BOOST_FOREACH(const face_descriptor &f, faces) {
       center = center + (center_pmap[f] - CGAL::ORIGIN) * area_pmap[f];
       sum_areas += area_pmap[f];
     }
@@ -61,42 +60,40 @@ struct Compact_metric_point_proxy {
 };
 
 typedef CGAL::Variational_shape_approximation<
-  Polyhedron, Vertex_point_map, Compact_metric_point_proxy> Compact_approx;
+  Mesh, Vertex_point_map, Compact_metric_point_proxy> Compact_approx;
 
 /**
  * This file tests the user defined metric.
  */
 int main()
 {
-  Polyhedron mesh;
-  std::ifstream input("./data/cube_meshed_open.off");
-  if (!input || !(input >> mesh) || mesh.empty()) {
-    std::cerr << "Invalid off file." << std::endl;
+  Mesh mesh;
+  std::ifstream input("./data/cube-ouvert.off");
+  if (!input || !(input >> mesh) || !CGAL::is_triangle_mesh(mesh)) {
+    std::cerr << "Invalid input file." << std::endl;
     return EXIT_FAILURE;
   }
 
-  // construct face normal & area map
-  std::map<Facet_handle, FT> face_areas;
-  std::map<Facet_handle, Point> face_centers;
-  for(Facet_iterator fitr = mesh.facets_begin(); fitr != mesh.facets_end(); ++fitr) {
-    const Halfedge_handle he = fitr->halfedge();
-    const Point &p0 = he->opposite()->vertex()->point();
-    const Point &p1 = he->vertex()->point();
-    const Point &p2 = he->next()->vertex()->point();
-    FT area(std::sqrt(CGAL::to_double(CGAL::squared_area(p0, p1, p2))));
-    face_areas.insert(std::pair<Facet_handle, FT>(fitr, area));
-    face_centers.insert(std::pair<Facet_handle, Point>(fitr, CGAL::centroid(p0, p1, p2)));
+  // construct face normal and area map
+  Vertex_point_map vpmap = get(boost::vertex_point, const_cast<Mesh &>(mesh));
+  Face_area_map area_pmap =
+    mesh.add_property_map<face_descriptor, FT>("f:area", FT(0.0)).first;
+  Face_center_map center_pmap =
+    mesh.add_property_map<face_descriptor, Point_3>("f:center", CGAL::ORIGIN).first;
+  BOOST_FOREACH (face_descriptor f, faces(mesh)) {
+    const halfedge_descriptor he = halfedge(f, mesh);
+    const Point_3 &p0 = vpmap[source(he, mesh)];
+    const Point_3 &p1 = vpmap[target(he, mesh)];
+    const Point_3 &p2 = vpmap[target(next(he, mesh), mesh)];
+    put(area_pmap, f, FT(std::sqrt(CGAL::to_double(CGAL::squared_area(p0, p1, p2)))));
+    put(center_pmap, f, CGAL::centroid(p0, p1, p2));
   }
-  Face_area_map area_pmap(face_areas);
-  Face_center_map center_pmap(face_centers);
 
   // create compact metric approximation algorithm instance
   std::cout << "create compact vas instance" << std::endl;
   Compact_metric_point_proxy error_metric(center_pmap, area_pmap);
 
-  Compact_approx approx(mesh,
-    get(boost::vertex_point, const_cast<Polyhedron &>(mesh)),
-    error_metric);
+  Compact_approx approx(mesh, vpmap, error_metric);
 
   std::cout << "random seeding and run" << std::endl;
   approx.initialize_seeds(CGAL::parameters::seeding_method(CGAL::Surface_mesh_approximation::RANDOM)
diff --git a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_segmentation_test.cpp b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_segmentation_test.cpp
index e3970a2766a..ff1a9ac3474 100644
--- a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_segmentation_test.cpp
+++ b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_segmentation_test.cpp
@@ -2,36 +2,35 @@
 #include <fstream>
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/remesh.h>
+
 #include <CGAL/Surface_mesh_approximation/approximate_triangle_mesh.h>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
-typedef typename boost::graph_traits<Polyhedron>::face_descriptor face_descriptor;
-typedef boost::unordered_map<face_descriptor, std::size_t> Face_index_map;
-typedef boost::associative_property_map<Face_index_map> Face_proxy_pmap;
+typedef CGAL::Surface_mesh<Kernel::Point_3> Mesh;
+typedef typename boost::graph_traits<Mesh>::face_descriptor face_descriptor;
+typedef Mesh::Property_map<face_descriptor, std::size_t> Face_proxy_pmap;
 
 /**
  * This file tests the free function CGAL::Surface_mesh_approximation::approximate_triangle_mesh.
  */
 int main()
 {
-  Polyhedron input;
-  std::ifstream file("data/sphere_iso.off");
-  if (!file || !(file >> input) || input.empty()) {
-    std::cerr << "Invalid off file." << std::endl;
+  Mesh mesh;
+  std::ifstream file("data/sphere.off");
+  if (!file || !(file >> mesh) || !CGAL::is_triangle_mesh(mesh)) {
+    std::cerr << "Invalid input file." << std::endl;
     return EXIT_FAILURE;
   }
 
-  Face_index_map fidx_map;
-  BOOST_FOREACH(face_descriptor f, faces(input))
-    fidx_map[f] = 0;
-  Face_proxy_pmap fpxmap(fidx_map);
-
+  Face_proxy_pmap fpxmap =
+    mesh.add_property_map<face_descriptor, std::size_t>("f:proxy_id", 0).first;
   std::vector<Kernel::Vector_3> proxies;
 
   // free function interface with named parameters
-  CGAL::Surface_mesh_approximation::approximate_triangle_mesh(input,
+  CGAL::Surface_mesh_approximation::approximate_triangle_mesh(mesh,
     CGAL::parameters::seeding_method(CGAL::Surface_mesh_approximation::HIERARCHICAL). // hierarchical seeding
     max_number_of_proxies(200). // both maximum number of proxies stop criterion,
     min_error_drop(0.05). // and minimum error drop stop criterion are specified
diff --git a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_teleportation_test.cpp b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_teleportation_test.cpp
index a77382bb9b1..3ba28a79a72 100644
--- a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_teleportation_test.cpp
+++ b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_teleportation_test.cpp
@@ -3,27 +3,27 @@
 #include <ctime>
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Polyhedron_3.h>
-#include <CGAL/IO/Polyhedron_iostream.h>
+#include <CGAL/Surface_mesh.h>
 #include <CGAL/Bbox_3.h>
 
-#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
+#include <CGAL/Polygon_mesh_processing/remesh.h>
+
 #include <CGAL/Variational_shape_approximation.h>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
 typedef Kernel::FT FT;
-typedef Kernel::Point_3 Point;
-typedef Kernel::Vector_3 Vector;
+typedef Kernel::Point_3 Point_3;
+typedef Kernel::Vector_3 Vector_3;
 
-typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
-typedef Polyhedron::Facet_handle Facet_handle;
-typedef Polyhedron::Facet_iterator Facet_iterator;
-typedef Polyhedron::Halfedge_handle Halfedge_handle;
+typedef CGAL::Surface_mesh<Kernel::Point_3> Mesh;
+typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
+typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
+typedef boost::graph_traits<Mesh>::halfedge_descriptor halfedge_descriptor;
 
-typedef boost::property_map<Polyhedron, boost::vertex_point_t>::type Vertex_point_map;
-typedef boost::associative_property_map<std::map<Facet_handle, std::size_t> > Face_proxy_map;
+typedef boost::property_map<Mesh, boost::vertex_point_t>::type Vertex_point_map;
+typedef Mesh::Property_map<face_descriptor, std::size_t> Face_proxy_map;
 
-typedef CGAL::Variational_shape_approximation<Polyhedron, Vertex_point_map> L21_approx;
+typedef CGAL::Variational_shape_approximation<Mesh, Vertex_point_map> L21_approx;
 typedef L21_approx::Error_metric L21_metric;
 typedef L21_approx::Proxy Plane_proxies;
 
@@ -35,10 +35,12 @@ bool check_strict_ordering(const std::vector<FT> &error)
     std::cout << "Empty error sequence." << std::endl;
     return false;
   }
+
   FT pre = error.front();
-  for (std::vector<FT>::const_iterator itr = error.begin(); itr != error.end(); ++itr)
-    if (pre < *itr)
+  BOOST_FOREACH(const FT &e, error) {
+    if (pre < e)
       return false;
+  }
 
   return true;
 }
@@ -50,20 +52,18 @@ bool check_strict_ordering(const std::vector<FT> &error)
  */
 int main()
 {
-  Polyhedron mesh;
+  Mesh mesh;
   std::ifstream input("./data/plane-sphere-high.off");
-  if (!input || !(input >> mesh) || mesh.empty()) {
-    std::cerr << "Invalid off file." << std::endl;
+  if (!input || !(input >> mesh) || !CGAL::is_triangle_mesh(mesh)) {
+    std::cerr << "Invalid input file." << std::endl;
     return EXIT_FAILURE;
   }
   std::cout << "Teleportation test." << std::endl;
 
   // algorithm instance
-  L21_metric error_metric(mesh,
-    get(boost::vertex_point, const_cast<Polyhedron &>(mesh)));
-  L21_approx approx(mesh,
-    get(boost::vertex_point, const_cast<Polyhedron &>(mesh)),
-    error_metric);
+  Vertex_point_map vpmap = get(boost::vertex_point, const_cast<Mesh &>(mesh));
+  L21_metric error_metric(mesh, vpmap);
+  L21_approx approx(mesh, vpmap, error_metric);
 
   std::cout << "Random seeding by number." << std::endl;
   std::srand(static_cast<unsigned int>(std::time(0)));
@@ -92,15 +92,16 @@ int main()
 
   // test partition placement
   std::cout << "Test partition placement." << std::endl;
-  std::map<Facet_handle, std::size_t> internal_fidxmap;
-  for (Facet_iterator fitr = mesh.facets_begin(); fitr != mesh.facets_end(); ++fitr)
-    internal_fidxmap[fitr] = 0;
-  Face_proxy_map fproxymap(internal_fidxmap);
+  Face_proxy_map fproxymap =
+    mesh.add_property_map<face_descriptor, std::size_t>("f:porxy_id", 0).first;
   approx.proxy_map(fproxymap);
   std::vector<Plane_proxies> proxies;
   approx.proxies(std::back_inserter(proxies));
 
-  CGAL::Bbox_3 bbox = CGAL::bbox_3(mesh.points_begin(), mesh.points_end());
+
+  CGAL::Bbox_3 bbox;
+  BOOST_FOREACH(const vertex_descriptor v, vertices(mesh))
+    bbox += vpmap[v].bbox();
   const FT ymin = bbox.ymin(), ymax = bbox.ymax(), yrange = ymax - ymin;
   std::cout << "Range along y axis: [" << ymin << ", " << ymax << "]" << std::endl;
 
@@ -108,21 +109,20 @@ int main()
   std::size_t planar_pxidx = static_cast<std::size_t>(-1);
   std::size_t num_planar_faces = 0;
   bool first = true;
-  for (Facet_iterator fitr = mesh.facets_begin(); fitr != mesh.facets_end(); ++fitr) {
-    Halfedge_handle he = fitr->halfedge();
-    const Point &p0 = he->opposite()->vertex()->point();
-    const Point &p1 = he->vertex()->point();
-    const Point &p2 = he->next()->vertex()->point();
-    const Point fcenter = CGAL::centroid(p0, p1, p2);
-    Vector fnormal = CGAL::normal(p0, p1, p2);
-    fnormal = fnormal / FT(std::sqrt(CGAL::to_double(fnormal.squared_length())));
+  BOOST_FOREACH(const face_descriptor f, faces(mesh)) {
+    const halfedge_descriptor he = halfedge(f, mesh);
+    const Point_3 &p0 = vpmap[source(he, mesh)];
+    const Point_3 &p1 = vpmap[target(he, mesh)];
+    const Point_3 &p2 = vpmap[target(next(he, mesh), mesh)];
+    const Point_3 fcenter = CGAL::centroid(p0, p1, p2);
+    const Vector_3 fnormal = CGAL::unit_normal(p0, p1, p2);
 
     // check the face center and normal to see if it is on the planar part of the geometry
     double dis_var = std::abs(CGAL::to_double((fcenter.y() - ymin) / yrange));
     double dir_var = std::abs(CGAL::to_double(fnormal.y()) - 1.0);
     if (dis_var < CGAL_VSA_TEST_TOLERANCE && dir_var < CGAL_VSA_TEST_TOLERANCE) {
       ++num_planar_faces;
-      const std::size_t pxidx = fproxymap[fitr];
+      const std::size_t pxidx = fproxymap[f];
       if (first) {
         first = false;
         planar_pxidx = pxidx;