diff --git a/Optimal_transportation_reconstruction_2/demo/Optimal_transportation_reconstruction_2/render.cpp b/Optimal_transportation_reconstruction_2/demo/Optimal_transportation_reconstruction_2/render.cpp
index 46837ca5b03..649f280928a 100644
--- a/Optimal_transportation_reconstruction_2/demo/Optimal_transportation_reconstruction_2/render.cpp
+++ b/Optimal_transportation_reconstruction_2/demo/Optimal_transportation_reconstruction_2/render.cpp
@@ -172,9 +172,9 @@ void R_s_k_2::draw_edge_footpoints(const Triangulation& mesh,
   Sample_vector::const_iterator it;
   for (it = samples.begin(); it != samples.end(); ++it)
   {
-    Sample_* sample = *it;
-    Point p = sample->point();
-    FT m = 0.5*(1.0 - sample->mass());
+    const Sample_& sample = this->m_samples[* it];
+    Point p = sample.point();
+    FT m = 0.5*(1.0 - sample.mass());
 
     Point q;
     if (mesh.get_plan(edge) == 0)
@@ -188,7 +188,7 @@ void R_s_k_2::draw_edge_footpoints(const Triangulation& mesh,
     else
     {
       viewer->glColor3f(red + m, green + m, blue + m);
-      FT t = sample->coordinate();
+      FT t = sample.coordinate();
       q = CGAL::ORIGIN + (1.0 - t)*(a - CGAL::ORIGIN) + t*(b - CGAL::ORIGIN);
     }
     draw_segment(p, q);
@@ -396,8 +396,8 @@ void R_s_k_2::draw_bins_plan0(const Edge& edge)
   Sample_vector_const_iterator it;
   for (it = samples.begin(); it != samples.end(); ++it)
   {
-    Sample_* sample = *it;
-    const Point& ps = sample->point();
+    const Sample_& sample = this->m_samples[* it];
+    const Point& ps = sample.point();
 
     Point q = pa;
     FT Da = CGAL::squared_distance(ps, pa);
@@ -423,8 +423,8 @@ void R_s_k_2::draw_bins_plan1(const Edge& edge)
     PSample psample = queue.top();
     queue.pop();
 
-    const FT m = psample.sample()->mass();
-    const Point& ps = psample.sample()->point();
+    const FT m = this->m_samples[psample.sample()].mass();
+    const Point& ps = this->m_samples[psample.sample()].point();
 
     FT bin = m/M;
     FT alpha = start + 0.5*bin;
@@ -511,7 +511,7 @@ void R_s_k_2::draw_one_ring(const float point_size, const float line_width, cons
   bool ok = locate_edge(query, edge);
   if (!ok) return;
 
-  Triangulation copy;
+  Triangulation copy(this->m_samples);
   Edge copy_edge = copy_star(edge, copy);
   draw_mesh_one_ring(point_size, line_width, copy, copy_edge);
 }
@@ -550,7 +550,7 @@ void R_s_k_2::draw_blocking_edges(const float point_size, const float line_width
   bool ok = locate_edge(query, edge);
   if (!ok) return;
 
-  Triangulation copy;
+  Triangulation copy(this->m_samples);
   Edge copy_edge = copy_star(edge, copy);
   draw_mesh_blocking_edges(point_size, line_width, copy, copy_edge);
 }
@@ -590,7 +590,7 @@ void R_s_k_2::draw_collapsible_edge(const float point_size,
   bool ok = locate_edge(query, edge);
   if (!ok) return;
 
-  Triangulation copy;
+  Triangulation copy(this->m_samples);
   Edge copy_edge = copy_star(edge, copy);
   Vertex_handle copy_src = copy.source_vertex(copy_edge);
 
@@ -612,7 +612,7 @@ void R_s_k_2::draw_cost_stencil(const float point_size,
   bool ok = locate_edge(query, edge);
   if (!ok) return;
 
-  Triangulation copy;
+  Triangulation copy(this->m_samples);
   Edge copy_edge = copy_star(edge, copy);
   Vertex_handle copy_src = copy.source_vertex(copy_edge);
 
@@ -709,7 +709,7 @@ void R_s_k_2::draw_push_queue_stencil(const float point_size,
       it++;
   }
 
-  Triangulation copy;
+  Triangulation copy(this->m_samples);
   Edge_vector copy_hull;
   Edge_vector copy_stencil;
   Edge copy_edge = copy_star(edge, copy);
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 e590f0f0489..5ea9b8ba055 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
@@ -74,8 +74,8 @@ public:
   typedef R_s_2::Edge_vector Edge_vector;
 
   typedef R_s_2::Sample_ Sample_;
-  typedef R_s_2::Sample_vector Sample_vector;
-  typedef R_s_2::Sample_vector_const_iterator Sample_vector_const_iterator;
+  typedef std::vector<Sample_> Sample_vector;
+  typedef Sample_vector::const_iterator Sample_vector_const_iterator;
 
   typedef R_s_2::PSample PSample;
   typedef R_s_2::SQueue SQueue;
@@ -454,7 +454,7 @@ public:
     for (std::vector<Sample_>::iterator it = m_samples.begin();
       it != m_samples.end(); ++it) {
       Sample_& s = *it;
-      samples.push_back(&s);
+      samples.push_back(s);
     }
 
     if (filename.contains(".xy", Qt::CaseInsensitive)) {
@@ -471,8 +471,8 @@ public:
     std::ofstream ofs(qPrintable(filename));
     for (Sample_vector_const_iterator it = samples.begin();
       it != samples.end(); ++it) {
-      Sample_* sample = *it;
-      ofs << sample->point() << std::endl;
+      const Sample_& sample = *it;
+      ofs << sample.point() << std::endl;
     }
     ofs.close();
   }
@@ -507,12 +507,12 @@ public:
     Sample_vector_const_iterator it;
     for (it = vertices.begin(); it != vertices.end(); it++) {
       vertices_mass_list.push_back(
-        std::make_pair((*it)->point(), (*it)->mass()));
+        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()));
+        std::make_pair((*it).point(), (*it).mass()));
     }
 
     Point_property_map point_pmap;
@@ -553,10 +553,10 @@ public:
     for (it = m_samples.begin(); it != m_samples.end(); ++it) {
       Sample_& s = *it;
 
-      samples.push_back(&s);
+      samples.push_back(s);
       FT rv = random.get_double(0.0, 1.0);
       if (rv <= percentage)
-        vertices.push_back(&s);
+        vertices.push_back(s);
     }
   }
 
diff --git a/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Sample.h b/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Sample.h
index fd4b6885a56..fee340a5914 100644
--- a/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Sample.h
+++ b/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Sample.h
@@ -39,6 +39,9 @@ private:
   FT m_backup_coord;
 
 public:
+  Sample()
+  {}
+
   Sample(const Point& point,
       const FT mass = FT(1))
   : m_point(point),
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 259b8f556e9..e59ee1de3b8 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
@@ -175,7 +175,6 @@ protected:
   MassPMap  mass_pmap;
 
 public:
-
   /// \name Initialization
   /// @{
 
@@ -319,8 +318,9 @@ public:
 
   /// \cond SKIP_IN_MANUAL
 
+
   Optimal_transportation_reconstruction_2()
-  : m_traits(m_dt.geom_traits())
+    : m_samples(), m_dt(m_samples), m_traits(m_dt.geom_traits())
   {
     initialize_parameters();
   }
@@ -1312,8 +1312,8 @@ public:
       PSample psample = queue.top();
       queue.pop();
 
-      const FT m = psample.sample()->mass();
-      const Point& ps = psample.sample()->point();
+      const FT m = this->m_samples[psample.sample()].mass();
+      const Point& ps = this->m_samples[psample.sample()].point();
 
       // normal + tangnetial
       const FT coord = psample.priority();