diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h
index f5f24d096b6..64cecfb5ea0 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h
@@ -29,6 +29,7 @@
 #include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
 #include <CGAL/algorithm.h>
 #include <set>
+#include <boost/dynamic_bitset.hpp>
 
 #include <boost/range/size.hpp>
 #include <boost/range/value_type.hpp>
@@ -67,13 +68,29 @@ public:
       _polygons(polygons)
   { }
 
-  void operator()(PM& pmesh)
+  void operator()(PM& pmesh, const bool insert_isolated_vertices = true)
   {
     Vpmap vpmap = get(CGAL::vertex_point, pmesh);
 
+    boost::dynamic_bitset<> not_isolated;
+    if (!insert_isolated_vertices)
+    {
+      not_isolated.resize(_points.size());
+      for (std::size_t i = 0, end = _polygons.size(); i < end; ++i)
+      {
+        const Polygon& polygon = _polygons[i];
+        const std::size_t size = polygon.size();
+        for (std::size_t j = 0; j < size; ++j)
+          not_isolated.set(polygon[j], true);
+      }
+    }
+
     std::vector<vertex_descriptor> vertices(_points.size());
     for (std::size_t i = 0, end = _points.size(); i < end; ++i)
     {
+      if (!insert_isolated_vertices && !not_isolated.test(i))
+        continue;
+
       Point_3 pi(_points[i][0], _points[i][1], _points[i][2]);
       vertices[i] = add_vertex(pmesh);
       put(vpmap, vertices[i], pi);
@@ -162,6 +179,8 @@ public:
     typename Orienter::Marked_edges marked_edges;
     Orienter::fill_edge_map(edges, marked_edges, polygons);
     //returns false if duplication is necessary
+    if (!marked_edges.empty())
+      return false;
     return Orienter::has_singular_vertices(static_cast<std::size_t>(max_id+1),polygons,edges,marked_edges);
   }
 
diff --git a/Polyhedron/demo/Polyhedron/Plugins/IO/Surf_io_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/IO/Surf_io_plugin.cpp
index 7574bdb077f..0f274e8f610 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/IO/Surf_io_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/IO/Surf_io_plugin.cpp
@@ -15,6 +15,7 @@
 
 #include "Color_map.h"
 #include <fstream>
+#include <boost/container/flat_set.hpp>
 
 using namespace CGAL::Three;
 class Surf_io_plugin:
@@ -62,13 +63,19 @@ CGAL::Three::Scene_item* Surf_io_plugin::load(QFileInfo fileinfo)
   std::vector<MaterialData> material_data;
   CGAL::Bbox_3 grid_box;
   CGAL::cpp11::array<unsigned int, 3> grid_size = {{1, 1, 1}};
-  read_surf(in, patches, material_data, grid_box, grid_size);
+  boost::container::flat_set<Polyhedron::Point_3> duplicated_points;
+  read_surf(in, patches, material_data, grid_box, grid_size
+          , std::inserter(duplicated_points, duplicated_points.end()));
+
   for(std::size_t i=0; i<material_data.size(); ++i)
   {
    std::cout<<"The patch #"<<i<<":\n  -inner region : material's id = "<<material_data[i].innerRegion.first<<" material's name = "
            <<material_data[i].innerRegion.second<<"\n  -outer region: material's id = "<<material_data[i].outerRegion.first<<" material's name = "
              <<material_data[i].outerRegion.second<<std::endl;
   }
+  if (!duplicated_points.empty())
+    std::cout << duplicated_points.size() << " points have been duplicated." << std::endl;
+  
   std::vector<QColor> colors_;
   compute_color_map(QColor(100, 100, 255), static_cast<unsigned>(patches.size()),
                     std::back_inserter(colors_));
diff --git a/Polyhedron/demo/Polyhedron/Plugins/IO/Surf_to_sm_io_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/IO/Surf_to_sm_io_plugin.cpp
index bca3516e5b1..3dcc124a932 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/IO/Surf_to_sm_io_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/IO/Surf_to_sm_io_plugin.cpp
@@ -14,6 +14,7 @@
 
 #include "Color_map.h"
 #include <fstream>
+#include <boost/container/flat_set.hpp>
 
 using namespace CGAL::Three;
 class Surf_io_plugin:
@@ -51,6 +52,7 @@ public:
 CGAL::Three::Scene_item* Surf_io_plugin::load(QFileInfo fileinfo)
 {
   typedef Scene_surface_mesh_item::SMesh SMesh;
+  typedef Scene_surface_mesh_item::Point Point;
   // Open file
   std::ifstream in(fileinfo.filePath().toUtf8());
   if(!in) {
@@ -62,13 +64,19 @@ CGAL::Three::Scene_item* Surf_io_plugin::load(QFileInfo fileinfo)
   std::vector<MaterialData> material_data;
   CGAL::Bbox_3 grid_box;
   CGAL::cpp11::array<unsigned int, 3> grid_size = {{1, 1, 1}};
-  read_surf(in, patches, material_data, grid_box, grid_size);
-  for(std::size_t i=0; i<material_data.size(); ++i)
+  boost::container::flat_set<Point> duplicated_points;
+  read_surf(in, patches, material_data, grid_box, grid_size
+    , std::inserter(duplicated_points, duplicated_points.end()));
+
+  for (std::size_t i = 0; i<material_data.size(); ++i)
   {
    std::cout<<"The patch #"<<i<<":\n  -inner region : material's id = "<<material_data[i].innerRegion.first<<" material's name = "
            <<material_data[i].innerRegion.second<<"\n  -outer region: material's id = "<<material_data[i].outerRegion.first<<" material's name = "
              <<material_data[i].outerRegion.second<<std::endl;
   }
+  if (!duplicated_points.empty())
+    std::cout << duplicated_points.size() << " points have been duplicated." << std::endl;
+
   std::vector<QColor> colors_;
   compute_color_map(QColor(100, 100, 255), static_cast<unsigned>(patches.size()),
                     std::back_inserter(colors_));
diff --git a/Polyhedron/demo/Polyhedron/include/CGAL/IO/read_surf_trianglemesh.h b/Polyhedron/demo/Polyhedron/include/CGAL/IO/read_surf_trianglemesh.h
index b4c784f2fd0..48895652521 100644
--- a/Polyhedron/demo/Polyhedron/include/CGAL/IO/read_surf_trianglemesh.h
+++ b/Polyhedron/demo/Polyhedron/include/CGAL/IO/read_surf_trianglemesh.h
@@ -77,16 +77,18 @@ bool line_starts_with(const std::string& line, const char* cstr)
  * read_surf reads a file which extension is .surf and fills `output` with one `Mesh` per patch.
  * Mesh is a model of FaceListGraph.
  */
-template<class Mesh, class NamedParameters>
+template<class Mesh, typename DuplicatedPointsOutIterator, class NamedParameters>
 bool read_surf(std::istream& input, std::vector<Mesh>& output,
     std::vector<MaterialData>& metadata,
     CGAL::Bbox_3& grid_box,
     CGAL::cpp11::array<unsigned int, 3>& grid_size,
+    DuplicatedPointsOutIterator out,
     const NamedParameters&)
 {
   typedef typename CGAL::GetGeomTraits<Mesh,
       NamedParameters>::type Kernel;
   typedef typename Kernel::Point_3 Point_3;
+  typedef typename boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
   std::vector<Point_3> points;
   std::string line;
   std::istringstream iss;
@@ -264,33 +266,45 @@ bool read_surf(std::istream& input, std::vector<Mesh>& output,
     {
       std::cout << "Orientation of patch #" << (i + 1) << "...";
       std::cout.flush();
+
+      const std::size_t nbp_init = points.size();
       bool no_duplicates =
       PMP::orient_polygon_soup(points, polygons);//returns false if some points
                                                  //were duplicated
+
       std::cout << "\rOrientation of patch #" << (i + 1) << " done";
-      if (!no_duplicates)
-        std::cout << " (non manifold -> duplicated vertices)";
+
+      if(!no_duplicates) //collect duplicates
+      { 
+        for (std::size_t i = nbp_init; i < points.size(); ++i)
+          *out++ = points[i];
+        std::cout << " (non manifold -> "
+                  << (points.size() - nbp_init) << " duplicated vertices)";
+      }
       std::cout << "." << std::endl;
     }
+
     Mesh& mesh = output[i];
 
-    CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(
-      points, polygons, mesh);
-    CGAL::Polygon_mesh_processing::remove_isolated_vertices(mesh);
+    PMP::internal::Polygon_soup_to_polygon_mesh<Mesh, Point_3, Triangle_ind>
+      converter(points, polygons);
+    converter(mesh, false/*insert_isolated_vertices*/);
 
+    CGAL_assertion(PMP::remove_isolated_vertices(mesh) == 0);
     CGAL_assertion(is_valid(mesh));
   } // end loop on patches
 
   return true;
 }
 
-template<class Mesh>
+template<class Mesh, typename DuplicatedPointsOutIterator>
 bool read_surf(std::istream& input, std::vector<Mesh>& output,
   std::vector<MaterialData>& metadata,
   CGAL::Bbox_3& grid_box,
-  CGAL::cpp11::array<unsigned int, 3>& grid_size)
+  CGAL::cpp11::array<unsigned int, 3>& grid_size,
+  DuplicatedPointsOutIterator out)
 {
-  return read_surf(input, output, metadata, grid_box, grid_size,
+  return read_surf(input, output, metadata, grid_box, grid_size, out,
     CGAL::Polygon_mesh_processing::parameters::all_default());
 }