diff --git a/Optimal_transportation_reconstruction_2/demo/Optimal_transportation_reconstruction_2/scene.h b/Optimal_transportation_reconstruction_2/demo/Optimal_transportation_reconstruction_2/scene.h
index d3e73c63800..1f749f0a310 100644
--- a/Optimal_transportation_reconstruction_2/demo/Optimal_transportation_reconstruction_2/scene.h
+++ b/Optimal_transportation_reconstruction_2/demo/Optimal_transportation_reconstruction_2/scene.h
@@ -415,10 +415,15 @@ public:
     Sample_vector vertices, samples;
     select_samples(percentage, vertices, samples);
 
-    PointMassList point_mass_list;
+    PointMassList vertices_mass_list;
     Sample_vector_const_iterator it;
     for (it = vertices.begin(); it != vertices.end(); it++) {
-      point_mass_list.push_back(
+      vertices_mass_list.push_back(
+        std::make_pair((*it)->point(), (*it)->mass()));
+    }
+    PointMassList samples_mass_list;
+    for (it = samples.begin(); it != samples.end(); it++) {
+      samples_mass_list.push_back(
         std::make_pair((*it)->point(), (*it)->mass()));
     }
 
@@ -426,8 +431,10 @@ public:
     Mass_property_map mass_pmap;
     MassPoint mp;
 
-    m_pwsrec->initialize(point_mass_list.begin(), point_mass_list.end(),
-      point_pmap, mass_pmap);
+    m_pwsrec->initialize_with_custom_vertices
+      (samples_mass_list.begin(), samples_mass_list.end(),
+       vertices_mass_list.begin(), vertices_mass_list.end(),
+       point_pmap, mass_pmap);
 
     m_init_done = true;
 
diff --git a/Optimal_transportation_reconstruction_2/include/CGAL/Optimal_transportation_reconstruction_2.h b/Optimal_transportation_reconstruction_2/include/CGAL/Optimal_transportation_reconstruction_2.h
index e60d18e4492..7a2a1014d7c 100644
--- a/Optimal_transportation_reconstruction_2/include/CGAL/Optimal_transportation_reconstruction_2.h
+++ b/Optimal_transportation_reconstruction_2/include/CGAL/Optimal_transportation_reconstruction_2.h
@@ -380,6 +380,34 @@ public:
     assign_samples(m_samples.begin(), m_samples.end());
   }
 
+  template <class InputIterator>
+  void initialize_with_custom_vertices(InputIterator samples_start,
+                                       InputIterator samples_beyond,
+                                       InputIterator vertices_start,
+                                       InputIterator vertices_beyond,
+                                       PointPMap point_map,
+                                       MassPMap  mass_map) {
+    point_pmap = point_map;
+    mass_pmap  = mass_map;
+    clear();
+    insert_loose_bbox(m_bbox_x, m_bbox_y, 2 * m_bbox_size);
+    init(vertices_start, vertices_beyond);
+
+    std::vector<Sample_*> m_samples;
+    for (InputIterator it = samples_start; it != samples_beyond; it++) {
+#ifdef CGAL_USE_PROPERTY_MAPS_API_V1
+      Point point = get(point_pmap, it);
+      FT    mass  = get( mass_pmap, it);
+#else
+      Point point = get(point_pmap, *it);
+      FT    mass  = get( mass_pmap, *it);
+#endif
+      Sample_* s = new Sample_(point, mass);
+      m_samples.push_back(s);
+    }
+    assign_samples(m_samples.begin(), m_samples.end());
+  }
+
 
   template <class Vector>
   Vector random_vec(const FT scale) const