Optimal Transport Reconstruction: Fix memory leak

Optimal_transportation_reconstruction_2/examples/Optimal_transportation_reconstruction_2/otr2_simplest_example.cpp

@@ -16,6 +16,7 @@ typedef CGAL::Optimal_transportation_reconstruction_2<K>    Otr;
 
 int main ()
 {
+    CGAL::get_default_random() = CGAL::Random(1671586136);
   // Generate a set of random points on the boundary of a square.
   std::vector<Point> points;
   CGAL::Random_points_on_square_2<Point> point_generator(1.);

Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Cost.h

@@ -50,13 +50,13 @@ public:
 
   const FT total_weight() const { return m_total_weight; }
 
-  template <typename SampleContainer>
+  template <typename Samples, typename SampleContainer>
-  void set_total_weight(const SampleContainer& samples)
+  void set_total_weight(const Samples& m_samples, const SampleContainer& samples)
   {
     m_total_weight = (FT)0;
     for (typename SampleContainer::const_iterator it = samples.begin();
          it != samples.end(); ++ it)
-      m_total_weight += (*it)->mass();
+      m_total_weight += m_samples[*it].mass();
   }
 
   FT finalize(const FT alpha = FT(0.5)) const

Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_face_base_2.h

@@ -16,7 +16,6 @@
 
 
 #include <CGAL/OTR_2/Cost.h>
-#include <CGAL/OTR_2/Sample.h>
 #include <CGAL/Triangulation_face_base_2.h>
 
 #include <vector>
@@ -51,7 +50,7 @@ public:
   typedef typename Traits_::FT    FT;
   typedef OTR_2::Cost<FT>          Cost_;
   typedef OTR_2::Sample<Traits_>   Sample_;
-  typedef std::vector<Sample_*>   Sample_vector;
+  typedef std::vector<int>   Sample_vector;
 
 private:
   Sample_vector m_samples[3];
@@ -176,7 +175,7 @@ public:
   const Sample_vector& samples(int edge) const { return m_samples[edge]; }
   Sample_vector& samples(int edge) { return m_samples[edge]; }
 
-  void add_sample(int edge, Sample_* sample)
+  void add_sample(int edge, int sample)
   {
     m_samples[edge].push_back(sample);
   }

Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_triangulation_2.h

@@ -110,7 +110,7 @@ public:
 
   typedef OTR_2::Cost<FT> Cost_;
   typedef OTR_2::Sample<Traits_> Sample_;
-  typedef std::vector<Sample_*> Sample_vector;
+  typedef std::vector<int> Sample_vector;
   typedef typename Sample_vector::const_iterator Sample_vector_const_iterator;
 
   typedef OTR_2::Sample_with_priority<Sample_> PSample;
@@ -135,12 +135,13 @@ public:
   >
   > MultiIndex;
 
+  std::vector<Sample_>& m_samples;
   FT m_factor; // ghost vs solid
   mutable Random rng;
 
 public:
-  Reconstruction_triangulation_2(Traits_ traits = Traits_())
+  Reconstruction_triangulation_2(std::vector<Sample_>& samples, Traits_ traits = Traits_())
-  : Base(traits), m_factor(1.)
+  : Base(traits), m_samples(samples), m_factor(1.)
   {
   }
 
@@ -360,11 +361,11 @@ public:
       if (cleanup)
         face->clean_all_samples();
     }
-    Sample_* sample = vertex->sample();
+    int sample = vertex->sample();
-    if (sample)
+    if (sample != -1)
       samples.push_back(sample);
     if (cleanup)
-      vertex->set_sample(nullptr);
+      vertex->set_sample(-1);
   }
 
   void collect_all_samples(Sample_vector& samples) const {
@@ -391,7 +392,7 @@ public:
     }
     for (Finite_vertices_iterator vi = Base::finite_vertices_begin();
         vi != Base::finite_vertices_end(); ++vi) {
-      vi->set_sample(nullptr);
+      vi->set_sample(-1);
     }
   }
 
@@ -465,15 +466,15 @@ public:
     typename Sample_vector::const_iterator it;
     const Sample_vector& samples0 = edge.first->samples(edge.second);
     for (it = samples0.begin(); it != samples0.end(); ++it) {
-      Sample_* sample = *it;
+      const Sample_ & sample = m_samples[* it];
-      mass += sample->mass();
+      mass += sample.mass();
     }
 
     Edge twin = twin_edge(edge);
     const Sample_vector& samples1 = twin.first->samples(twin.second);
     for (it = samples1.begin(); it != samples1.end(); ++it) {
-      Sample_* sample = *it;
+      const Sample_& sample = m_samples[* it];
-      mass += sample->mass();
+      mass += sample.mass();
     }
 
     set_mass(edge, mass);
@@ -511,15 +512,15 @@ public:
     typename Sample_vector::const_iterator it;
     const Sample_vector& samples0 = edge.first->samples(edge.second);
     for (it = samples0.begin(); it != samples0.end(); ++it) {
-      Sample_* sample = *it;
+      const Sample_& sample = m_samples[* it];
-      squeue.push(PSample(sample, sample->coordinate()));
+      squeue.push(PSample(*it, sample.coordinate()));
     }
 
     Edge twin = twin_edge(edge);
     const Sample_vector& samples1 = twin.first->samples(twin.second);
     for (it = samples1.begin(); it != samples1.end(); ++it) {
-      Sample_* sample = *it;
+      const Sample_& sample = m_samples[* it];
-      squeue.push(PSample(sample, 1.0 - sample->coordinate()));
+      squeue.push(PSample(*it, 1.0 - sample.coordinate()));
     }
   }
 
@@ -537,13 +538,13 @@ public:
       PSample psample = squeue.top();
       squeue.pop();
 
-      FT mass = psample.sample()->mass();
+      FT mass = m_samples[psample.sample()].mass();
       FT coord = psample.priority() * L;
       FT bin = mass * coef;
       FT center = start + FT(0.5) * bin;
       FT pos = coord - center;
 
-      FT norm2 = psample.sample()->distance2();
+      FT norm2 = m_samples[psample.sample()].distance2();
       FT tang2 = bin * bin / 12 + pos * pos;
 
       sum.add(Cost_(norm2, tang2), mass);
@@ -566,15 +567,15 @@ public:
     Cost_ sum;
     for (Sample_vector_const_iterator it = samples.begin();
         it != samples.end(); ++it) {
-      Sample_* sample = *it;
+      const Sample_& sample =  m_samples[* it];
-      FT mass = sample->mass();
+      FT mass = sample.mass();
-      const Point& query = sample->point();
+      const Point& query = sample.point();
 
       FT Ds = geom_traits().compute_squared_distance_2_object()(query, ps);
       FT Dt = geom_traits().compute_squared_distance_2_object()(query, pt);
       FT dist2 = ((std::min))(Ds, Dt);
 
-      FT norm2 = sample->distance2();
+      FT norm2 = sample.distance2();
       FT tang2 = dist2 - norm2;
 
       sum.add(Cost_(norm2, tang2), mass);
@@ -589,7 +590,7 @@ public:
   template<class Iterator> // value_type = Sample_*
   void assign_samples(Iterator begin, Iterator end) {
     for (Iterator it = begin; it != end; ++it) {
-      Sample_* sample = *it;
+      int sample = *it;
       assign_sample(sample);
     }
   }
@@ -597,13 +598,13 @@ public:
   template<class Iterator> // value_type = Sample_*
   void assign_samples_brute_force(Iterator begin, Iterator end) {
     for (Iterator it = begin; it != end; ++it) {
-      Sample_* sample = *it;
+      int sample = *it;
       assign_sample_brute_force(sample);
     }
   }
 
-  bool assign_sample(Sample_* sample) {
+  bool assign_sample(int sample) {
-    const Point& point = sample->point();
+    const Point& point = m_samples[sample].point();
     Face_handle face = Base::locate(point);
 
     if (face == Face_handle() || Base::is_infinite(face)) {
@@ -622,8 +623,9 @@ public:
     return true;
   }
 
-  bool assign_sample_brute_force(Sample_* sample) {
+  bool assign_sample_brute_force(int sample_index) {
-    const Point& point = sample->point();
+    const Sample_& sample = m_samples[sample_index];
+    const Point& point = sample.point();
     Face_handle nearest_face = Face_handle();
     for (Finite_faces_iterator fi = Base::finite_faces_begin();
         fi != Base::finite_faces_end(); ++fi) {
@@ -641,12 +643,12 @@ public:
 
     Vertex_handle vertex = find_nearest_vertex(point, nearest_face);
     if (vertex != Vertex_handle()) {
-      assign_sample_to_vertex(sample, vertex);
+      assign_sample_to_vertex(sample_index, vertex);
       return true;
     }
 
     Edge edge = find_nearest_edge(point, nearest_face);
-    assign_sample_to_edge(sample, edge);
+    assign_sample_to_edge(sample_index, edge);
     return true;
   }
 
@@ -688,23 +690,24 @@ public:
     return nearest;
   }
 
-  void assign_sample_to_vertex(Sample_* sample, Vertex_handle vertex) const {
+  void assign_sample_to_vertex(int sample_index, Vertex_handle vertex) const {
     /*if (vertex->sample()) {
       std::cout << "assign to vertex: vertex already has sample"
           << std::endl;
     }*/
-
+    Sample_& sample = m_samples[sample_index];
-    sample->distance2() = FT(0);
+    sample.distance2() = FT(0);
-    sample->coordinate() = FT(0);
+    sample.coordinate() = FT(0);
-    vertex->set_sample(sample);
+    vertex->set_sample(sample_index);
   }
 
-  void assign_sample_to_edge(Sample_* sample, const Edge& edge) const {
+  void assign_sample_to_edge(int sample_index, const Edge& edge) const {
+    Sample_& sample = m_samples[sample_index];
     Segment segment = get_segment(edge);
-    const Point& query = sample->point();
+    const Point& query = sample.point();
-    sample->distance2() = compute_distance2(query, segment);
+    sample.distance2() = compute_distance2(query, segment);
-    sample->coordinate() = compute_coordinate(query, segment);
+    sample.coordinate() = compute_coordinate(query, segment);
-    edge.first->add_sample(edge.second, sample);
+    edge.first->add_sample(edge.second, sample_index);
   }
 
   FT compute_distance2(const Point& query, const Segment& segment) const {

Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_vertex_base_2.h

@@ -37,7 +37,7 @@ class Reconstruction_vertex_base_2 : public Vb
 public:
   typedef Vb Base;
   typedef typename Traits_::FT        FT;
-  typedef OTR_2::Sample<Traits_>       Sample_;
+  typedef OTR_2::Sample<Traits_>      Sample_;
   typedef typename Traits_::Point_2   Point;
   typedef typename Base::Face_handle  Face_handle;
 
@@ -50,7 +50,7 @@ public:
 private:
   int       m_id;
   bool      m_pinned;
-  Sample_*  m_sample;
+  int       m_sample;
   Point     m_relocated;
   FT        m_relevance;
 
@@ -60,7 +60,7 @@ public:
   : Base(),
     m_id(-1),
     m_pinned(false),
-    m_sample(nullptr),
+    m_sample(-1),
     m_relevance(0)
 {
 }
@@ -69,7 +69,7 @@ public:
   : Base(p),
     m_id(-1),
     m_pinned(false),
-    m_sample(nullptr),
+    m_sample(-1),
     m_relevance(0)
   {
   }
@@ -78,7 +78,7 @@ public:
   : Base(f),
     m_id(-1),
     m_pinned(false),
-    m_sample(nullptr),
+    m_sample(-1),
     m_relevance(0)
   {
   }
@@ -87,7 +87,7 @@ public:
   : Base(p, f),
     m_id(-1),
     m_pinned(false),
-    m_sample(nullptr),
+    m_sample(-1),
     m_relevance(0)
   {
   }
@@ -103,13 +103,13 @@ public:
   FT relevance() const { return m_relevance; }
   void set_relevance(FT relevance) { m_relevance = relevance; }
 
-  Sample_* sample() const { return m_sample; }
+  int sample() const { return m_sample; }
-  void set_sample(Sample_* sample) { m_sample = sample; }
+  void set_sample(int sample) { m_sample = sample; }
 
   const Point& relocated() const { return m_relocated; }
   Point& relocated() { return m_relocated; }
 
-  bool  has_sample_assigned() const { return sample() != nullptr; }
+  bool  has_sample_assigned() const { return sample() != -1; }
 };
 //---------------STRUCT LESS VERTEX_HANDLE---------------------
 template <class T>

Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Sample.h

@@ -94,23 +94,20 @@ public:
   typedef typename Sample_::FT FT;
 
 private:
-  Sample_* m_sample;
+  int m_sample;
   FT m_priority;
 
 public:
-  Sample_with_priority(Sample_* sample, const FT priority = FT(0))
+  Sample_with_priority(int sample, const FT priority = FT(0))
-  {
+    : m_sample(sample), m_priority(priority)
-    m_sample   = sample;
+  {}
-    m_priority = priority;
-  }
 
   Sample_with_priority(const Sample_with_priority& psample)
-  {
+    : m_sample(psample.sample()), m_priority(psample.priority())
-    m_sample   = psample.sample();
+  {}
-    m_priority = psample.priority();
-  }
 
-  ~Sample_with_priority() { }
+  ~Sample_with_priority()
+  {}
 
   Sample_with_priority& operator = (const Sample_with_priority& psample)
   {
@@ -119,7 +116,7 @@ public:
     return *this;
   }
 
-  Sample_* sample() const { return m_sample; }
+  int sample() const { return m_sample; }
 
   const FT priority() const { return m_priority; }
 };

Optimal_transportation_reconstruction_2/include/CGAL/Optimal_transportation_reconstruction_2.h

@@ -117,7 +117,6 @@ public:
     The Output simplex.
    */
   typedef OTR_2::Reconstruction_triangulation_2<Traits>  Triangulation;
-
   typedef typename Triangulation::Vertex                Vertex;
   typedef typename Triangulation::Vertex_handle         Vertex_handle;
   typedef typename Triangulation::Vertex_iterator       Vertex_iterator;
@@ -159,6 +158,7 @@ public:
   /// @}
 
 protected:
+  std::vector<Sample_> m_samples;
   Triangulation m_dt;
   Traits const& m_traits;
   MultiIndex m_mindex;
@@ -212,7 +212,8 @@ public:
     unsigned int relocation = 2,
     int verbose = 0,
     Traits traits = Traits())
-  : m_dt(traits),
+    : m_samples(),
+    m_dt(m_samples, traits),
     m_traits(m_dt.geom_traits()),
     m_ignore(0),
     m_verbose(verbose),
@@ -369,14 +370,18 @@ public:
     insert_loose_bbox(bbox);
     init(start, beyond);
 
-    std::vector<Sample_*> m_samples;
+    m_samples.reserve(std::distance(start,beyond));
     for (InputIterator it = start; it != beyond; it++) {
       Point point = get(point_pmap, *it);
       FT    mass  = get( mass_pmap, *it);
-      Sample_* s = new Sample_(point, mass);
+      Sample_ s(point, mass);
       m_samples.push_back(s);
     }
-    assign_samples(m_samples.begin(), m_samples.end());
+    Sample_vector sv(m_samples.size());
+    for(int i = 0; i < sv.size(); ++i){
+      sv[i] = i;
+    }
+    assign_samples(sv.begin(), sv.end());
   }
 
   template <class InputIterator>
@@ -398,7 +403,7 @@ public:
     insert_loose_bbox(bbox);
     init(vertices_start, vertices_beyond);
 
-    std::vector<Sample_*> m_samples;
+    m_samples.reserve(std::distance(start,beyond));
     for (InputIterator it = samples_start; it != samples_beyond; it++) {
 #ifdef CGAL_USE_PROPERTY_MAPS_API_V1
       Point point = get(point_pmap, it);
@@ -407,10 +412,14 @@ public:
       Point point = get(point_pmap, *it);
       FT    mass  = get( mass_pmap, *it);
 #endif
-      Sample_* s = new Sample_(point, mass);
+      Sample_s(point, mass);
       m_samples.push_back(s);
     }
-    assign_samples(m_samples.begin(), m_samples.end());
+    Sample_vector sv(m_samples.size());
+    for(int i = 0; i < sv.size(); ++i){
+      sv[i] = i;
+    }
+    assign_samples(sv.begin(), sv.end());
   }
 
 
@@ -422,16 +431,8 @@ public:
     return m_traits.construct_vector_2_object()(dx, dy);
   }
 
-  void clear() {
+  void clear()
-    Sample_vector samples;
+  {}
-    m_dt.collect_all_samples(samples);
-    // Deallocate samples
-    for (Sample_vector_const_iterator s_it = samples.begin();
-        s_it != samples.end(); ++s_it)
-    {
-      delete *s_it;
-    }
-  }
 
 
   // INIT //
@@ -494,7 +495,7 @@ public:
     m_dt.cleanup_assignments();
   }
 
-  template<class Iterator>  // value_type = Sample_*
+  template<class Iterator>  // value_type = int
   void assign_samples(Iterator begin, Iterator end) {
     CGAL::Real_timer timer;
     if (m_verbose > 0)
@@ -587,7 +588,7 @@ public:
         << s->id() << "->" << t->id() << ") ... " << std::endl;
     }
 
-    Triangulation copy;
+    Triangulation copy(m_samples);
     Edge copy_edge = copy_star(edge, copy);
     Vertex_handle copy_source = copy.source_vertex(copy_edge);
 
@@ -632,7 +633,7 @@ public:
     copy.assign_samples_brute_force(samples.begin(), samples.end());
     copy.reset_all_costs();
     cost = copy.compute_total_cost();
-    cost.set_total_weight (samples);
+    cost.set_total_weight (m_samples, samples);
     restore_samples(samples.begin(), samples.end());
 
     if (m_verbose > 1) {
@@ -643,18 +644,18 @@ public:
   }
 
   template<class Iterator> // value_type = Sample_*
-  void backup_samples(Iterator begin, Iterator end) const {
+  void backup_samples(Iterator begin, Iterator end) {
     for (Iterator it = begin; it != end; ++it) {
-      Sample_* sample = *it;
+      Sample_& sample = m_samples[* it];
-      sample->backup();
+      sample.backup();
     }
   }
 
   template<class Iterator> // value_type = Sample_*
-  void restore_samples(Iterator begin, Iterator end) const {
+  void restore_samples(Iterator begin, Iterator end) {
     for (Iterator it = begin; it != end; ++it) {
-      Sample_* sample = *it;
+      Sample_& sample = m_samples[* it];
-      sample->restore();
+      sample.restore();
     }
   }
 
@@ -1089,7 +1090,7 @@ public:
       m_dt.collect_samples_from_edge(twin, samples);
       copy_twin.first->samples(copy_twin.second) = samples;
     }
-    copy_vertex->set_sample(nullptr);
+    copy_vertex->set_sample(-1);
   }
 
   Edge get_copy_edge(
@@ -1230,10 +1231,10 @@ public:
   void compute_relocation_for_vertex(
     Vertex_handle vertex, FT& coef, Vector& rhs) const
   {
-    Sample_* sample = vertex->sample();
+    if (vertex->sample() != -1) {
-    if (sample) {
+      const Sample_& sample = m_samples[vertex->sample()];
-      const FT m = sample->mass();
+      const FT m = sample.mass();
-      const Point& ps = sample->point();
+      const Point& ps = sample.point();
       rhs = m_traits.construct_sum_of_vectors_2_object()(rhs,
         m_traits.construct_scaled_vector_2_object()(
           m_traits.construct_vector_2_object()(CGAL::ORIGIN, ps), m));
@@ -1253,9 +1254,9 @@ public:
     Vector grad = m_traits.construct_vector_2_object()(FT(0), FT(0));
     Sample_vector_const_iterator it;
     for (it = samples.begin(); it != samples.end(); ++it) {
-      Sample_* sample = *it;
+      const Sample_& sample = m_samples[* it];
-      const FT m = sample->mass();
+      const FT m = sample.mass();
-      const Point& ps = sample->point();
+      const Point& ps = sample.point();
 
       FT Da = m_traits.compute_squared_distance_2_object()(ps, pa);
       FT Db = m_traits.compute_squared_distance_2_object()(ps, pb);
@@ -1281,9 +1282,9 @@ public:
 
     Sample_vector_const_iterator it;
     for (it = samples.begin(); it != samples.end(); ++it) {
-      Sample_* sample = *it;
+      const Sample_& sample = m_samples[* it];
-      const FT m = sample->mass();
+      const FT m = sample.mass();
-      const Point& ps = sample->point();
+      const Point& ps = sample.point();
 
       FT Da = m_traits.compute_squared_distance_2_object()(ps, pa);
       FT Db = m_traits.compute_squared_distance_2_object()(ps, pb);
@@ -1357,8 +1358,8 @@ public:
       PSample psample = queue.top();
       queue.pop();
 
-      const FT m = psample.sample()->mass();
+      const FT m = m_samples[psample.sample()].mass();
-      const Point& ps = psample.sample()->point();
+      const Point& ps = m_samples[psample.sample()].point();
 
       const FT coord = psample.priority();
       const FT one_minus_coord = 1.0 - coord;