diff --git a/Optimal_transportation_reconstruction_2/examples/Optimal_transportation_reconstruction_2/otr2_simplest_example.cpp b/Optimal_transportation_reconstruction_2/examples/Optimal_transportation_reconstruction_2/otr2_simplest_example.cpp index cef2ecc34b0..a9f49e12d7b 100644 --- a/Optimal_transportation_reconstruction_2/examples/Optimal_transportation_reconstruction_2/otr2_simplest_example.cpp +++ b/Optimal_transportation_reconstruction_2/examples/Optimal_transportation_reconstruction_2/otr2_simplest_example.cpp @@ -16,6 +16,7 @@ typedef CGAL::Optimal_transportation_reconstruction_2 Otr; int main () { + CGAL::get_default_random() = CGAL::Random(1671586136); // Generate a set of random points on the boundary of a square. std::vector points; CGAL::Random_points_on_square_2 point_generator(1.); diff --git a/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Cost.h b/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Cost.h index 937de05c73a..0c60d82e830 100644 --- a/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Cost.h +++ b/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 - void set_total_weight(const SampleContainer& samples) + template + 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 diff --git a/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_face_base_2.h b/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_face_base_2.h index 6903891e6a4..39c994d671d 100644 --- a/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_face_base_2.h +++ b/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_face_base_2.h @@ -16,7 +16,6 @@ #include -#include #include #include @@ -51,7 +50,7 @@ public: typedef typename Traits_::FT FT; typedef OTR_2::Cost Cost_; typedef OTR_2::Sample Sample_; - typedef std::vector Sample_vector; + typedef std::vector 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); } diff --git a/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_triangulation_2.h b/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_triangulation_2.h index f952c9232e8..16fe600bb34 100644 --- a/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_triangulation_2.h +++ b/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_triangulation_2.h @@ -110,7 +110,7 @@ public: typedef OTR_2::Cost Cost_; typedef OTR_2::Sample Sample_; - typedef std::vector Sample_vector; + typedef std::vector Sample_vector; typedef typename Sample_vector::const_iterator Sample_vector_const_iterator; typedef OTR_2::Sample_with_priority PSample; @@ -135,12 +135,13 @@ public: > > MultiIndex; + std::vector& m_samples; FT m_factor; // ghost vs solid mutable Random rng; public: - Reconstruction_triangulation_2(Traits_ traits = Traits_()) - : Base(traits), m_factor(1.) + Reconstruction_triangulation_2(std::vector& samples, Traits_ traits = Traits_()) + : Base(traits), m_samples(samples), m_factor(1.) { } @@ -360,11 +361,11 @@ public: if (cleanup) face->clean_all_samples(); } - Sample_* sample = vertex->sample(); - if (sample) + int sample = vertex->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; - mass += sample->mass(); + const Sample_ & sample = m_samples[* it]; + 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; - mass += sample->mass(); + const Sample_& sample = m_samples[* it]; + 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; - squeue.push(PSample(sample, sample->coordinate())); + const Sample_& sample = m_samples[* it]; + 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; - squeue.push(PSample(sample, 1.0 - sample->coordinate())); + const Sample_& sample = m_samples[* it]; + 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; - FT mass = sample->mass(); - const Point& query = sample->point(); + const Sample_& sample = m_samples[* it]; + FT mass = sample.mass(); + 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 // 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 // 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) { - const Point& point = sample->point(); + bool assign_sample(int sample) { + 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) { - const Point& point = sample->point(); + bool assign_sample_brute_force(int sample_index) { + 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->distance2() = FT(0); - sample->coordinate() = FT(0); - vertex->set_sample(sample); + Sample_& sample = m_samples[sample_index]; + sample.distance2() = FT(0); + sample.coordinate() = FT(0); + 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(); - sample->distance2() = compute_distance2(query, segment); - sample->coordinate() = compute_coordinate(query, segment); - edge.first->add_sample(edge.second, sample); + const Point& query = sample.point(); + sample.distance2() = compute_distance2(query, segment); + sample.coordinate() = compute_coordinate(query, segment); + edge.first->add_sample(edge.second, sample_index); } FT compute_distance2(const Point& query, const Segment& segment) const { diff --git a/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_vertex_base_2.h b/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_vertex_base_2.h index 6155f65e5fe..7df7b10b047 100644 --- a/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Reconstruction_vertex_base_2.h +++ b/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 Sample_; + typedef OTR_2::Sample 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; } - void set_sample(Sample_* sample) { m_sample = sample; } + int sample() const { return m_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 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 7972b69f02c..fd4b6885a56 100644 --- a/Optimal_transportation_reconstruction_2/include/CGAL/OTR_2/Sample.h +++ b/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)) - { - m_sample = sample; - m_priority = priority; - } + Sample_with_priority(int sample, const FT priority = FT(0)) + : 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; } }; 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 9617930d52f..42445c8708d 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 @@ -117,7 +117,6 @@ public: The Output simplex. */ typedef OTR_2::Reconstruction_triangulation_2 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 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 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 @@ -398,7 +403,7 @@ public: insert_loose_bbox(bbox); init(vertices_start, vertices_beyond); - std::vector 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() { - 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; - } - } + void clear() + {} // INIT // @@ -494,7 +495,7 @@ public: m_dt.cleanup_assignments(); } - template // value_type = Sample_* + template // 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 // 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->backup(); + Sample_& sample = m_samples[* it]; + sample.backup(); } } template // 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->restore(); + Sample_& sample = m_samples[* it]; + 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 (sample) { - const FT m = sample->mass(); - const Point& ps = sample->point(); + if (vertex->sample() != -1) { + const Sample_& sample = m_samples[vertex->sample()]; + const FT m = sample.mass(); + 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 FT m = sample->mass(); - const Point& ps = sample->point(); + const Sample_& sample = m_samples[* it]; + const FT m = sample.mass(); + 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 FT m = sample->mass(); - const Point& ps = sample->point(); + const Sample_& sample = m_samples[* it]; + const FT m = sample.mass(); + 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 Point& ps = psample.sample()->point(); + const FT m = m_samples[psample.sample()].mass(); + const Point& ps = m_samples[psample.sample()].point(); const FT coord = psample.priority(); const FT one_minus_coord = 1.0 - coord;