songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Update branch from master after trailing whitespaces and tabs removal

This commit is contained in:
Sébastien Loriot 2020-03-26 19:29:05 +01:00
parent 16fc8d1fe2 929e1e4bce
commit a0345b135c
21 changed files with 1323 additions and 1585 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
Optimal_transportation_reconstruction_2/demo/Optimal_transportation_reconstruction_2/scene.h
View File

@ -92,7 +92,7 @@ private:
return K::Point_3 (sample.point().x(), sample.point().y(), 0.);
}
};
// data
std::vector<Sample_> m_samples;
@ -106,10 +106,10 @@ private:
double m_bbox_x;
double m_bbox_y;
double m_bbox_size;
//Random
CGAL::Random random;
template <class Vector>
Vector random_vec(const double scale)
{
@ -407,10 +407,10 @@ public:
Point_3_from_sample()),
boost::make_transform_iterator (m_samples.end(),
Point_3_from_sample())),
3);
3, CGAL::parameters::point_map (CGAL::Identity_property_map_no_lvalue<K::Point_3>()));
std::cerr << "Average spacing = " << spacing << std::endl;
}
void print_vertex(Vertex vertex) {
std::cout << "vertex " << vertex << std::endl;
}
@ -639,7 +639,7 @@ public:
if (view_tolerance)
draw_tolerance(viewer);
if (view_edges)
m_pwsrec->draw_edges(0.5f * line_thickness, 0.9f, 0.9f, 0.9f);

184
Point_set_processing_3/include/CGAL/Point_set_processing_3/internal/Callback_wrapper.h Normal file
View File

@ -0,0 +1,184 @@
// Copyright (c) 2018 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Simon Giraudot
#ifndef CGAL_PSP_INTERNAL_CALLBACK_WRAPPER_H
#define CGAL_PSP_INTERNAL_CALLBACK_WRAPPER_H
#include <CGAL/license/Point_set_processing_3.h>
#include <functional>
#include <CGAL/thread.h>
namespace CGAL {
namespace Point_set_processing_3 {
namespace internal {
template <typename ConcurrencyTag>
class Callback_wrapper
{
std::size_t m_advancement;
bool m_interrupted;
public:
Callback_wrapper (const std::function<bool(double)>&, std::size_t, std::size_t, bool = false) { }
void reset (std::size_t, std::size_t, bool = false) { }
std::size_t& advancement() { return m_advancement; }
bool& interrupted() { return m_interrupted; }
void join() { }
};
template <>
class Callback_wrapper<CGAL::Sequential_tag>
{
const std::function<bool(double)>& m_callback;
std::size_t m_advancement;
bool m_interrupted;
std::size_t m_size;
public:
Callback_wrapper (const std::function<bool(double)>& callback,
std::size_t size,
std::size_t advancement = 0,
bool interrupted = false)
: m_callback (callback)
, m_advancement (advancement)
, m_interrupted (interrupted)
, m_size (size)
{ }
Callback_wrapper (const Callback_wrapper& other)
: m_callback (other.m_callback)
, m_advancement (other.m_advancement)
, m_interrupted (other.m_interrupted)
, m_size (other.m_size)
{ }
~Callback_wrapper () { }
void reset (std::size_t size, std::size_t advancement, bool interrupted = false)
{
m_size = size;
m_advancement = advancement;
m_interrupted = interrupted;
}
std::size_t& advancement()
{
return m_advancement;
}
bool& interrupted()
{
if (m_callback)
m_interrupted = !(m_callback(m_advancement / double(m_size)));
return m_interrupted;
}
void join() { }
};
#ifdef CGAL_LINKED_WITH_TBB
template <>
class Callback_wrapper<CGAL::Parallel_tag>
{
const std::function<bool(double)>& m_callback;
cpp11::atomic<std::size_t>* m_advancement;
cpp11::atomic<bool>* m_interrupted;
std::size_t m_size;
bool m_creator;
cpp11::thread* m_thread;
// assignment operator shouldn't be used (m_callback is const ref)
Callback_wrapper& operator= (const Callback_wrapper&)
{
return *this;
}
public:
Callback_wrapper (const std::function<bool(double)>& callback,
std::size_t size,
std::size_t advancement = 0,
bool interrupted = false)
: m_callback (callback)
, m_advancement (new cpp11::atomic<std::size_t>())
, m_interrupted (new cpp11::atomic<bool>())
, m_size (size)
, m_creator (true)
, m_thread (nullptr)
{
// cpp11::atomic only has default constructor, initialization done in two steps
*m_advancement = advancement;
*m_interrupted = interrupted;
if (m_callback)
m_thread = new cpp11::thread (*this);
}
Callback_wrapper (const Callback_wrapper& other)
: m_callback (other.m_callback)
, m_advancement (other.m_advancement)
, m_interrupted (other.m_interrupted)
, m_size (other.m_size)
, m_creator (false)
, m_thread (nullptr)
{
}
~Callback_wrapper ()
{
if (m_creator)
{
delete m_advancement;
delete m_interrupted;
}
if (m_thread != nullptr)
delete m_thread;
}
void reset (std::size_t size, std::size_t advancement, bool interrupted = false)
{
m_size = size;
*m_advancement = advancement;
*m_interrupted = interrupted;
if (m_callback)
m_thread = new cpp11::thread (*this);
}
cpp11::atomic<std::size_t>& advancement() { return *m_advancement; }
cpp11::atomic<bool>& interrupted() { return *m_interrupted; }
void join()
{
if (m_thread != nullptr)
m_thread->join();
}
void operator()()
{
while (*m_advancement != m_size)
{
if (m_callback && !m_callback (*m_advancement / double(m_size)))
*m_interrupted = true;
if (*m_interrupted)
return;
cpp11::sleep_for (0.00001);
}
if (m_callback)
m_callback (1.);
}
};
#endif
} // namespace internal
} // namespace Point_set_processing_3
} // namespace CGAL
#endif // CGAL_PSP_INTERNAL_CALLBACK_WRAPPER_H

180
Point_set_processing_3/include/CGAL/Point_set_processing_3/internal/Neighbor_query.h Normal file
View File

@ -0,0 +1,180 @@
// Copyright (c) 2019 GeometryFactory Sarl (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Simon Giraudot
#ifndef CGAL_PSP_INTERNAL_NEIGHBOR_QUERY_H
#define CGAL_PSP_INTERNAL_NEIGHBOR_QUERY_H
#include <CGAL/license/Point_set_processing_3.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Fuzzy_sphere.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/point_set_processing_assertions.h>
#include <CGAL/iterator.h>
#include <boost/function_output_iterator.hpp>
namespace CGAL {
namespace Point_set_processing_3 {
namespace internal {
struct Maximum_points_reached_exception : public std::exception { };
template <typename Kernel_, typename PointRangeRef, typename PointMap>
class Neighbor_query
{
public:
typedef Kernel_ Kernel;
typedef PointRangeRef Point_range;
typedef PointMap Point_map;
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point_3;
typedef typename Range_iterator_type<PointRangeRef>::type input_iterator;
typedef typename input_iterator::value_type value_type;
typedef CGAL::Prevent_deref<input_iterator> iterator;
struct Deref_point_map
{
typedef input_iterator key_type;
typedef typename boost::property_traits<PointMap>::reference reference;
typedef typename boost::property_traits<PointMap>::value_type value_type;
typedef typename boost::property_traits<PointMap>::category category;
PointMap point_map;
Deref_point_map () { }
Deref_point_map (PointMap point_map) : point_map(point_map) { }
friend reference get (const Deref_point_map& map, key_type it)
{
return get(map.point_map, *it);
}
};
typedef CGAL::Search_traits_3<Kernel> Tree_traits_base;
typedef CGAL::Search_traits_adapter<input_iterator, Deref_point_map, Tree_traits_base> Tree_traits;
typedef CGAL::Sliding_midpoint<Tree_traits> Splitter;
typedef CGAL::Distance_adapter<input_iterator, Deref_point_map, CGAL::Euclidean_distance<Tree_traits_base> > Distance;
typedef CGAL::Kd_tree<Tree_traits, Splitter, CGAL::Tag_true, CGAL::Tag_true> Tree;
typedef CGAL::Fuzzy_sphere<Tree_traits> Sphere;
typedef CGAL::Orthogonal_k_neighbor_search<Tree_traits, Distance, Splitter, Tree> Neighbor_search;
typedef typename Neighbor_search::iterator Search_iterator;
private:
PointRangeRef m_points;
PointMap m_point_map;
Deref_point_map m_deref_map;
Tree_traits m_traits;
Tree m_tree;
Distance m_distance;
// Forbid copy
Neighbor_query (const Neighbor_query&) { }
public:
Neighbor_query (PointRangeRef points, PointMap point_map)
: m_points (points)
, m_point_map (point_map)
, m_deref_map (point_map)
, m_traits (m_deref_map)
, m_tree (iterator(m_points.begin()), iterator(m_points.end()), Splitter(), m_traits)
, m_distance (m_deref_map)
{
m_tree.build();
}
PointMap point_map() const { return m_point_map; }
template <typename OutputIterator>
void get_iterators (const Point_3& query, unsigned int k, FT neighbor_radius,
OutputIterator output) const
{
if (neighbor_radius != FT(0))
{
Sphere fs (query, neighbor_radius, 0, m_traits);
// if k=0, no limit on the number of neighbors returned
if (k == 0)
k = (std::numeric_limits<unsigned int>::max)();
unsigned int nb = 0;
try
{
std::function<void(const input_iterator&)> output_iterator_with_limit
= [&](const input_iterator& it)
{
*(output ++) = it;
if (++ nb == k)
throw Maximum_points_reached_exception();
};
auto function_output_iterator
= boost::make_function_output_iterator (output_iterator_with_limit);
m_tree.search (function_output_iterator, fs);
}
catch (const Maximum_points_reached_exception&)
{ }
// Fallback, if less than 3 points are return, search for the 3
// first points
if (nb < 3)
k = 3;
// Else, no need to search for K nearest neighbors
else
k = 0;
}
if (k != 0)
{
// Gather set of (k+1) neighboring points.
// Perform k+1 queries (as in point set, the query point is
// output first). Search may be aborted if k is greater
// than number of input points.
Neighbor_search search (m_tree, query, k+1, 0, true, m_distance);
Search_iterator search_iterator = search.begin();
unsigned int i;
for (i = 0; i < (k+1); ++ i)
{
if(search_iterator == search.end())
break; // premature ending
*(output ++) = search_iterator->first;
search_iterator++;
}
}
}
template <typename OutputIterator>
void get_points (const Point_3& query, unsigned int k, FT neighbor_radius,
OutputIterator output) const
{
return get_iterators(query, k, neighbor_radius,
boost::make_function_output_iterator
([&](const input_iterator& it)
{
*(output ++) = get (m_point_map, *it);
}));
}
};
} } } // namespace CGAL::Point_set_processing_3::internal
#endif // CGAL_PSP_INTERNAL_NEIGHBOR_QUERY_H

108
Point_set_processing_3/include/CGAL/Point_set_processing_3/internal/Parallel_callback.h
View File

@ -1,108 +0,0 @@
// Copyright (c) 2018 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Simon Giraudot
#ifndef CGAL_PSP_INTERNAL_PARALLEL_CALLBACK_H
#define CGAL_PSP_INTERNAL_PARALLEL_CALLBACK_H
#include <CGAL/license/Point_set_processing_3.h>
#include <functional>
#include <CGAL/thread.h>
namespace CGAL {
namespace Point_set_processing_3 {
namespace internal {
class Parallel_callback
{
const std::function<bool(double)>& m_callback;
cpp11::atomic<std::size_t>* m_advancement;
cpp11::atomic<bool>* m_interrupted;
std::size_t m_size;
bool m_creator;
cpp11::thread* m_thread;
// assignment operator shouldn't be used (m_callback is const ref)
Parallel_callback& operator= (const Parallel_callback&)
{
return *this;
}
public:
Parallel_callback (const std::function<bool(double)>& callback,
std::size_t size,
std::size_t advancement = 0,
bool interrupted = false)
: m_callback (callback)
, m_advancement (new cpp11::atomic<std::size_t>())
, m_interrupted (new cpp11::atomic<bool>())
, m_size (size)
, m_creator (true)
, m_thread (nullptr)
{
// cpp11::atomic only has default constructor, initialization done in two steps
*m_advancement = advancement;
*m_interrupted = interrupted;
if (m_callback)
m_thread = new cpp11::thread (*this);
}
Parallel_callback (const Parallel_callback& other)
: m_callback (other.m_callback)
, m_advancement (other.m_advancement)
, m_interrupted (other.m_interrupted)
, m_size (other.m_size)
, m_creator (false)
, m_thread (nullptr)
{
}
~Parallel_callback ()
{
if (m_creator)
{
delete m_advancement;
delete m_interrupted;
}
if (m_thread != nullptr)
delete m_thread;
}
cpp11::atomic<std::size_t>& advancement() { return *m_advancement; }
cpp11::atomic<bool>& interrupted() { return *m_interrupted; }
void join()
{
if (m_thread != nullptr)
m_thread->join();
}
void operator()()
{
while (*m_advancement != m_size)
{
if (!m_callback (*m_advancement / double(m_size)))
*m_interrupted = true;
if (*m_interrupted)
return;
cpp11::sleep_for (0.00001);
}
m_callback (1.);
}
};
} // namespace internal
} // namespace Point_set_processing_3
} // namespace CGAL
#endif // CGAL_PSP_INTERNAL_PARALLEL_CALLBACK_H

104
Point_set_processing_3/include/CGAL/Point_set_processing_3/internal/neighbor_query.h
View File

@ -1,104 +0,0 @@
// Copyright (c) 2019 GeometryFactory Sarl (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Simon Giraudot
#ifndef CGAL_PSP_INTERNAL_NEIGHBOR_QUERY_H
#define CGAL_PSP_INTERNAL_NEIGHBOR_QUERY_H
#include <CGAL/license/Point_set_processing_3.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Fuzzy_sphere.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/point_set_processing_assertions.h>
#include <boost/function_output_iterator.hpp>
namespace CGAL {
namespace Point_set_processing_3 {
namespace internal {
struct Maximum_points_reached_exception : public std::exception { };
template <typename Point,
typename TreeTraits,
typename TreeSplitter,
typename TreeUseExtendedNode,
typename FT,
typename PointContainer>
void neighbor_query (const Point& query,
const CGAL::Kd_tree<TreeTraits, TreeSplitter, TreeUseExtendedNode>& tree,
unsigned int k,
FT neighbor_radius,
PointContainer& points)
{
typedef typename CGAL::Orthogonal_k_neighbor_search<TreeTraits> Neighbor_search;
typedef typename Neighbor_search::iterator Search_iterator;
typedef CGAL::Fuzzy_sphere<TreeTraits> Sphere;
if (neighbor_radius != FT(0))
{
Sphere fs (query, neighbor_radius, 0, tree.traits());
// if k=0, no limit on the number of neighbors returned
if (k == 0)
k = (std::numeric_limits<unsigned int>::max)();
try
{
std::function<void(const Point&)> back_insert_with_limit
= [&](const Point& point) -> void
{
points.push_back (point);
if (points.size() == k)
throw Maximum_points_reached_exception();
};
auto function_output_iterator
= boost::make_function_output_iterator (back_insert_with_limit);
tree.search (function_output_iterator, fs);
}
catch (const Maximum_points_reached_exception&)
{ }
// Fallback, if less than 3 points are return, search for the 3
// first points
if (points.size() < 3)
k = 3;
// Else, no need to search for K nearest neighbors
else
k = 0;
}
if (k != 0)
{
// Gather set of (k+1) neighboring points.
// Perform k+1 queries (as in point set, the query point is
// output first). Search may be aborted if k is greater
// than number of input points.
points.reserve(k+1);
Neighbor_search search(tree,query,k+1);
Search_iterator search_iterator = search.begin();
unsigned int i;
for(i=0;i<(k+1);i++)
{
if(search_iterator == search.end())
break; // premature ending
points.push_back(search_iterator->first);
search_iterator++;
}
CGAL_point_set_processing_precondition(points.size() >= 1);
}
}
} } } // namespace CGAL::Point_set_processing_3::internal
#endif // CGAL_PSP_INTERNAL_NEIGHBOR_QUERY_H

473
Point_set_processing_3/include/CGAL/bilateral_smooth_point_set.h
View File

@ -7,7 +7,7 @@
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Shihao Wu, Clement Jamin, Pierre Alliez
// Author(s) : Shihao Wu, Clement Jamin, Pierre Alliez
#ifndef CGAL_BILATERAL_SMOOTH_POINT_SET_H
#define CGAL_BILATERAL_SMOOTH_POINT_SET_H
@ -17,18 +17,19 @@
#include <CGAL/disable_warnings.h>
#include <CGAL/number_type_config.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Point_set_processing_3/internal/neighbor_query.h>
#include <CGAL/Point_set_processing_3/internal/Neighbor_query.h>
#include <CGAL/Point_set_processing_3/internal/Callback_wrapper.h>
#include <CGAL/for_each.h>
#include <CGAL/property_map.h>
#include <CGAL/point_set_processing_assertions.h>
#include <CGAL/Point_with_normal_3.h>
#include <CGAL/squared_distance_3.h>
#include <functional>
#include <CGAL/boost/graph/Named_function_parameters.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <boost/iterator/zip_iterator.hpp>
#include <iterator>
#include <set>
#include <algorithm>
@ -38,24 +39,7 @@
#include <CGAL/Memory_sizer.h>
#include <CGAL/property_map.h>
#ifdef CGAL_LINKED_WITH_TBB
#include <CGAL/Point_set_processing_3/internal/Parallel_callback.h>
#include <tbb/parallel_for.h>
#include <tbb/blocked_range.h>
#include <tbb/scalable_allocator.h>
#include <atomic>
#endif // CGAL_LINKED_WITH_TBB
// Default allocator: use TBB allocators if available
#ifdef CGAL_LINKED_WITH_TBB
# define CGAL_PSP3_DEFAULT_ALLOCATOR tbb::scalable_allocator
#else // CGAL_LINKED_WITH_TBB
# define CGAL_PSP3_DEFAULT_ALLOCATOR std::allocator
#endif // CGAL_LINKED_WITH_TBB
//#define CGAL_PSP3_VERBOSE
//#define CGAL_PSP3_VERBOSE
namespace CGAL {
@ -66,63 +50,24 @@ namespace CGAL {
namespace bilateral_smooth_point_set_internal{
// Item in the Kd-tree: position (Point_3) + index
template <typename Kernel>
class Kd_tree_element : public Point_with_normal_3<Kernel>
{
public:
unsigned int index;
// basic geometric types
typedef typename CGAL::Origin Origin;
typedef CGAL::Point_with_normal_3<Kernel> Base;
Kd_tree_element(const Origin& o = ORIGIN, unsigned int id=0)
: Base(o), index(id)
{}
Kd_tree_element(const Base& p, unsigned int id=0)
: Base(p), index(id)
{}
Kd_tree_element(const Kd_tree_element& other)
: Base(other), index(other.index)
{}
Kd_tree_element& operator=(const Kd_tree_element&)=default;
};
// Helper class for the Kd-tree
template <typename Kernel>
class Kd_tree_gt : public Kernel
{
public:
typedef Kd_tree_element<Kernel> Point_3;
};
template <typename Kernel>
class Kd_tree_traits : public CGAL::Search_traits_3<Kd_tree_gt<Kernel> >
{
public:
typedef typename Kernel::Point_3 PointType;
};
/// Compute bilateral projection for each point
/// according to their KNN neighborhood points
///
///
/// \pre `k >= 2`, radius > 0 , sharpness_angle > 0 && sharpness_angle < 90
///
/// @tparam Kernel Geometric traits class.
/// @tparam Tree KD-tree.
///
/// @return
/// @return
template <typename Kernel>
CGAL::Point_with_normal_3<Kernel>
template <typename Kernel, typename PointRange,
typename PointMap, typename VectorMap>
std::pair<typename Kernel::Point_3, typename Kernel::Vector_3>
compute_denoise_projection(
const CGAL::Point_with_normal_3<Kernel>& query, ///< 3D point to project
const std::vector<CGAL::Point_with_normal_3<Kernel>,
CGAL_PSP3_DEFAULT_ALLOCATOR<CGAL::Point_with_normal_3<Kernel> > >& neighbor_pwns, //
const typename PointRange::iterator::value_type& vt,
PointMap point_map,
VectorMap normal_map,
const std::vector<typename PointRange::iterator>& neighbor_pwns,
typename Kernel::FT radius, ///< accept neighborhood radius
typename Kernel::FT sharpness_angle ///< control sharpness(0-90)
)
@ -133,7 +78,6 @@ compute_denoise_projection(
// basic geometric types
typedef typename Kernel::FT FT;
typedef CGAL::Point_with_normal_3<Kernel> Pwn;
typedef typename Kernel::Vector_3 Vector;
typedef typename Kernel::Point_3 Point;
@ -143,28 +87,26 @@ compute_denoise_projection(
FT iradius16 = -(FT)4.0/radius2;
FT project_dist_sum = FT(0.0);
FT project_weight_sum = FT(0.0);
Vector normal_sum = CGAL::NULL_VECTOR;
Vector normal_sum = CGAL::NULL_VECTOR;
FT cos_sigma = cos(sharpness_angle * CGAL_PI / 180.0);
FT sharpness_bandwidth = std::pow((CGAL::max)(1e-8, 1 - cos_sigma), 2);
typename std::vector<Pwn,CGAL_PSP3_DEFAULT_ALLOCATOR<Pwn> >::const_iterator
pwn_iter = neighbor_pwns.begin();
for (; pwn_iter != neighbor_pwns.end(); ++pwn_iter)
for (typename PointRange::iterator it : neighbor_pwns)
{
const Point& np = pwn_iter->position();
const Vector& nn = pwn_iter->normal();
const Point& np = get(point_map, *it);
const Vector& nn = get(normal_map, *it);
FT dist2 = CGAL::squared_distance(query.position(), np);
FT dist2 = CGAL::squared_distance(get(point_map, vt), np);
if (dist2 < radius2)
{
FT theta = std::exp(dist2 * iradius16);
FT psi = std::exp(-std::pow(1 - query.normal() * nn, 2)
FT psi = std::exp(-std::pow(1 - get(normal_map, vt) * nn, 2)
/ sharpness_bandwidth);
weight = theta * psi;
project_dist_sum += ((query.position() - np) * nn) * weight;
project_dist_sum += ((get(point_map, vt) - np) * nn) * weight;
project_weight_sum += weight;
normal_sum = normal_sum + nn * weight;
}
@ -173,10 +115,10 @@ compute_denoise_projection(
Vector update_normal = normal_sum / project_weight_sum;
update_normal = update_normal / sqrt(update_normal.squared_length());
Point update_point = query.position() - update_normal *
(project_dist_sum / project_weight_sum);
Point update_point = get(point_map, vt) - update_normal *
(project_dist_sum / project_weight_sum);
return Pwn(update_point, update_normal);
return std::make_pair (update_point, update_normal);
}
/// Computes max-spacing of one query point from K nearest neighbors.
@ -187,41 +129,30 @@ compute_denoise_projection(
/// @tparam Tree KD-tree.
///
/// @return max spacing.
template < typename Kernel,
typename Tree >
typename Kernel::FT
template <typename NeighborQuery>
typename NeighborQuery::Kernel::FT
compute_max_spacing(
const CGAL::Point_with_normal_3<Kernel>& query, ///< 3D point
Tree& tree, ///< KD-tree
const typename NeighborQuery::value_type& vt,
typename NeighborQuery::Point_map point_map,
NeighborQuery& neighbor_query, ///< KD-tree
unsigned int k) ///< number of neighbors
{
// basic geometric types
typedef typename NeighborQuery::Kernel Kernel;
typedef typename Kernel::FT FT;
typedef CGAL::Point_with_normal_3<Kernel> Pwn;
// types for K nearest neighbors search
typedef bilateral_smooth_point_set_internal::Kd_tree_traits<Kernel> Tree_traits;
typedef CGAL::Orthogonal_k_neighbor_search<Tree_traits> Neighbor_search;
typedef typename Neighbor_search::iterator Search_iterator;
// performs k + 1 queries (if unique the query point is
// output first). search may be aborted when k is greater
// than number of input points
Neighbor_search search(tree,query,k+1);
Search_iterator search_iterator = search.begin();
++search_iterator;
FT max_distance = (FT)0.0;
unsigned int i;
for(i = 0; i < (k+1) ; ++i)
{
if(search_iterator == search.end())
break; // premature ending
Pwn pwn = search_iterator->first;
double dist2 = CGAL::squared_distance(query.position(), pwn.position());
max_distance = (CGAL::max)(dist2, max_distance);
++search_iterator;
}
neighbor_query.get_iterators
(get(point_map, vt), k, (FT)(0.0),
boost::make_function_output_iterator
([&](const typename NeighborQuery::input_iterator& it)
{
double dist2 = CGAL::squared_distance (get(point_map, vt), get(point_map, *it));
max_distance = (CGAL::max)(dist2, max_distance);
}));
// output max spacing
return std::sqrt(max_distance);
@ -231,102 +162,6 @@ compute_max_spacing(
/// \endcond
#ifdef CGAL_LINKED_WITH_TBB
/// \cond SKIP_IN_MANUAL
/// This is for parallelization of function: bilateral_smooth_point_set()
template <typename Kernel, typename Tree>
class Compute_pwns_neighbors
{
typedef typename CGAL::Point_with_normal_3<Kernel> Pwn;
typedef typename std::vector<Pwn,CGAL_PSP3_DEFAULT_ALLOCATOR<Pwn> > Pwns;
typedef typename std::vector<Pwns,CGAL_PSP3_DEFAULT_ALLOCATOR<Pwns> >
Pwns_neighbors;
typedef typename Kernel::FT FT;
unsigned int m_k;
FT m_neighbor_radius;
const Tree & m_tree;
const Pwns & m_pwns;
Pwns_neighbors & m_pwns_neighbors;
cpp11::atomic<std::size_t>& advancement;
cpp11::atomic<bool>& interrupted;
public:
Compute_pwns_neighbors(unsigned int k, FT neighbor_radius, const Tree &tree,
const Pwns &pwns, Pwns_neighbors &neighbors,
cpp11::atomic<std::size_t>& advancement,
cpp11::atomic<bool>& interrupted)
: m_k(k), m_neighbor_radius (neighbor_radius), m_tree(tree)
, m_pwns(pwns), m_pwns_neighbors(neighbors)
, advancement (advancement), interrupted (interrupted) {}
void operator() ( const tbb::blocked_range<size_t>& r ) const
{
for (size_t i = r.begin(); i!=r.end(); i++)
{
if (interrupted)
break;
CGAL::Point_set_processing_3::internal::neighbor_query
(m_pwns[i], m_tree, m_k, m_neighbor_radius, m_pwns_neighbors[i]);
++ advancement;
}
}
};
/// \endcond
/// \cond SKIP_IN_MANUAL
/// This is for parallelization of function: compute_denoise_projection()
template <typename Kernel>
class Pwn_updater
{
typedef typename CGAL::Point_with_normal_3<Kernel> Pwn;
typedef typename std::vector<Pwn,CGAL_PSP3_DEFAULT_ALLOCATOR<Pwn> > Pwns;
typedef typename Kernel::FT FT;
FT sharpness_angle;
FT radius;
Pwns* pwns;
Pwns* update_pwns;
std::vector<Pwns,CGAL_PSP3_DEFAULT_ALLOCATOR<Pwns> >* pwns_neighbors;
cpp11::atomic<std::size_t>& advancement;
cpp11::atomic<bool>& interrupted;
public:
Pwn_updater(FT sharpness,
FT r,
Pwns *in,
Pwns *out,
std::vector<Pwns,CGAL_PSP3_DEFAULT_ALLOCATOR<Pwns> >* neighbors,
cpp11::atomic<std::size_t>& advancement,
cpp11::atomic<bool>& interrupted):
sharpness_angle(sharpness),
radius(r),
pwns(in),
update_pwns(out),
pwns_neighbors(neighbors),
advancement (advancement),
interrupted (interrupted) {}
void operator() ( const tbb::blocked_range<size_t>& r ) const
{
for (size_t i = r.begin(); i != r.end(); ++i)
{
if (interrupted)
break;
(*update_pwns)[i] = bilateral_smooth_point_set_internal::
compute_denoise_projection<Kernel>((*pwns)[i],
(*pwns_neighbors)[i],
radius,
sharpness_angle);
++ advancement;
}
}
};
/// \endcond
#endif // CGAL_LINKED_WITH_TBB
// ----------------------------------------------------------------------------
// Public section
@ -335,17 +170,17 @@ public:
/**
\ingroup PkgPointSetProcessing3Algorithms
This function smooths an input point set by iteratively projecting each
This function smooths an input point set by iteratively projecting each
point onto the implicit surface patch fitted over its nearest neighbors.
Bilateral projection preserves sharp features according to the normal
(gradient) information. Both point positions and normals will be modified.
For more details, please see section 4 in \cgalCite{ear-2013}.
(gradient) information. Both point positions and normals will be modified.
For more details, please see section 4 in \cgalCite{ear-2013}.
A parallel version of this function is provided and requires the executable to be
A parallel version of this function is provided and requires the executable to be
linked against the <a href="https://www.threadingbuildingblocks.org">Intel TBB library</a>.
To control the number of threads used, the user may use the tbb::task_scheduler_init class.
See the <a href="https://www.threadingbuildingblocks.org/documentation">TBB documentation</a>
See the <a href="https://www.threadingbuildingblocks.org/documentation">TBB documentation</a>
for more details.
\pre Normals must be unit vectors
@ -398,79 +233,58 @@ bilateral_smooth_point_set(
{
using parameters::choose_parameter;
using parameters::get_parameter;
// basic geometric types
typedef typename PointRange::iterator iterator;
typedef typename iterator::value_type value_type;
typedef typename CGAL::GetPointMap<PointRange, NamedParameters>::type PointMap;
typedef typename Point_set_processing_3::GetNormalMap<PointRange, NamedParameters>::type NormalMap;
typedef typename Point_set_processing_3::GetK<PointRange, NamedParameters>::Kernel Kernel;
typedef typename Kernel::Point_3 Point_3;
typedef typename Kernel::Vector_3 Vector_3;
CGAL_static_assertion_msg(!(boost::is_same<NormalMap,
typename Point_set_processing_3::GetNormalMap<PointRange, NamedParameters>::NoMap>::value),
"Error: no normal map");
typedef typename CGAL::Point_with_normal_3<Kernel> Pwn;
typedef typename std::vector<Pwn,CGAL_PSP3_DEFAULT_ALLOCATOR<Pwn> > Pwns;
typedef typename Kernel::FT FT;
double sharpness_angle = choose_parameter(get_parameter(np, internal_np::sharpness_angle), 30.);
const std::function<bool(double)>& callback = choose_parameter(get_parameter(np, internal_np::callback),
std::function<bool(double)>());
CGAL_point_set_processing_precondition(points.begin() != points.end());
CGAL_point_set_processing_precondition(k > 1);
// types for K nearest neighbors search structure
typedef bilateral_smooth_point_set_internal::
Kd_tree_element<Kernel> Kd_tree_element;
typedef bilateral_smooth_point_set_internal::Kd_tree_traits<Kernel> Tree_traits;
typedef CGAL::Orthogonal_k_neighbor_search<Tree_traits> Neighbor_search;
typedef typename Neighbor_search::Tree Tree;
typedef Point_set_processing_3::internal::Neighbor_query<Kernel, PointRange&, PointMap> Neighbor_query;
PointMap point_map = choose_parameter<PointMap>(get_parameter(np, internal_np::point_map));
NormalMap normal_map = choose_parameter<NormalMap>(get_parameter(np, internal_np::normal_map));
FT neighbor_radius = choose_parameter(get_parameter(np, internal_np::neighbor_radius), FT(0));
// copy points and normals
Pwns pwns;
for(typename PointRange::iterator it = points.begin(); it != points.end(); ++it)
{
typename boost::property_traits<PointMap>::reference p = get(point_map, *it);
typename boost::property_traits<NormalMap>::reference n = get(normal_map, *it);
CGAL_point_set_processing_precondition(n.squared_length() > 1e-10);
pwns.push_back(Pwn(p, n));
}
std::size_t nb_points = pwns.size();
std::size_t nb_points = points.size();
#ifdef CGAL_PSP3_VERBOSE
std::cout << "Initialization and compute max spacing: " << std::endl;
#endif
// initiate a KD-tree search for points
std::vector<Kd_tree_element,
CGAL_PSP3_DEFAULT_ALLOCATOR<Kd_tree_element> > treeElements;
treeElements.reserve(pwns.size());
typename std::vector<Pwn,CGAL_PSP3_DEFAULT_ALLOCATOR<Pwn> >::iterator
pwn_iter = pwns.begin();
for (unsigned int i = 0; pwn_iter != pwns.end(); ++pwn_iter)
{
treeElements.push_back(Kd_tree_element(*pwn_iter, i));
}
Tree tree(treeElements.begin(), treeElements.end());
Neighbor_query neighbor_query (points, point_map);
// Guess spacing
#ifdef CGAL_PSP3_VERBOSE
CGAL::Real_timer task_timer;
task_timer.start();
#endif
FT guess_neighbor_radius = 0.0;
FT guess_neighbor_radius = 0.0;
for(pwn_iter = pwns.begin(); pwn_iter != pwns.end(); ++pwn_iter)
for (const value_type& vt : points)
{
FT max_spacing = bilateral_smooth_point_set_internal::
compute_max_spacing<Kernel,Tree>(*pwn_iter, tree, k);
guess_neighbor_radius = (CGAL::max)(max_spacing, guess_neighbor_radius);
compute_max_spacing (vt, point_map, neighbor_query, k);
guess_neighbor_radius = (CGAL::max)(max_spacing, guess_neighbor_radius);
}
#ifdef CGAL_PSP3_VERBOSE
task_timer.stop();
#endif
@ -486,50 +300,41 @@ bilateral_smooth_point_set(
task_timer.start();
#endif
// compute all neighbors
std::vector<Pwns,CGAL_PSP3_DEFAULT_ALLOCATOR<Pwns> > pwns_neighbors;
typedef std::vector<iterator> iterators;
std::vector<iterators> pwns_neighbors;
pwns_neighbors.resize(nb_points);
Point_set_processing_3::internal::Callback_wrapper<ConcurrencyTag>
callback_wrapper (callback, 2 * nb_points);
#ifndef CGAL_LINKED_WITH_TBB
CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
"Parallel_tag is enabled but TBB is unavailable.");
#else
if (boost::is_convertible<ConcurrencyTag,Parallel_tag>::value)
{
Point_set_processing_3::internal::Parallel_callback
parallel_callback (callback, 2 * nb_points);
typedef boost::zip_iterator<boost::tuple<iterator, typename std::vector<iterators>::iterator> > Zip_iterator;
Compute_pwns_neighbors<Kernel, Tree> f(k, neighbor_radius, tree, pwns, pwns_neighbors,
parallel_callback.advancement(),
parallel_callback.interrupted());
tbb::parallel_for(tbb::blocked_range<size_t>(0, nb_points), f);
CGAL::for_each<ConcurrencyTag>
(CGAL::make_range (boost::make_zip_iterator (boost::make_tuple (points.begin(), pwns_neighbors.begin())),
boost::make_zip_iterator (boost::make_tuple (points.end(), pwns_neighbors.end()))),
[&](const typename Zip_iterator::reference& t)
{
if (callback_wrapper.interrupted())
return false;
bool interrupted = parallel_callback.interrupted();
neighbor_query.get_iterators (get(point_map, get<0>(t)), k, neighbor_radius,
std::back_inserter (get<1>(t)));
++ callback_wrapper.advancement();
// We interrupt by hand as counter only goes halfway and won't terminate by itself
parallel_callback.interrupted() = true;
parallel_callback.join();
return true;
});
// If interrupted during this step, nothing is computed, we return NaN
if (interrupted)
return std::numeric_limits<double>::quiet_NaN();
}
else
#endif
{
typename std::vector<Pwns,CGAL_PSP3_DEFAULT_ALLOCATOR<Pwns> >::iterator
pwns_iter = pwns_neighbors.begin();
std::size_t nb = 0;
for(pwn_iter = pwns.begin(); pwn_iter != pwns.end(); ++pwn_iter, ++pwns_iter, ++ nb)
{
CGAL::Point_set_processing_3::internal::neighbor_query
(*pwn_iter, tree, k, neighbor_radius, *pwns_iter);
if (callback && !callback ((nb+1) / double(2. * nb_points)))
return std::numeric_limits<double>::quiet_NaN();
}
}
bool interrupted = callback_wrapper.interrupted();
// We interrupt by hand as counter only goes halfway and won't terminate by itself
callback_wrapper.interrupted() = true;
callback_wrapper.join();
// If interrupted during this step, nothing is computed, we return NaN
if (interrupted)
return std::numeric_limits<double>::quiet_NaN();
#ifdef CGAL_PSP3_VERBOSE
task_timer.stop();
memory = CGAL::Memory_sizer().virtual_size();
@ -541,75 +346,67 @@ bilateral_smooth_point_set(
task_timer.start();
#endif
// update points and normals
Pwns update_pwns(nb_points);
std::vector<std::pair<Point_3, Vector_3> > update_pwns(nb_points);
#ifdef CGAL_LINKED_WITH_TBB
if(boost::is_convertible<ConcurrencyTag, CGAL::Parallel_tag>::value)
{
Point_set_processing_3::internal::Parallel_callback
parallel_callback (callback, 2 * nb_points, nb_points);
callback_wrapper.reset (2 * nb_points, nb_points);
//tbb::task_scheduler_init init(4);
tbb::blocked_range<size_t> block(0, nb_points);
Pwn_updater<Kernel> pwn_updater(sharpness_angle,
guess_neighbor_radius,
&pwns,
&update_pwns,
&pwns_neighbors,
parallel_callback.advancement(),
parallel_callback.interrupted());
tbb::parallel_for(block, pwn_updater);
typedef boost::zip_iterator
<boost::tuple<iterator,
typename std::vector<iterators>::iterator,
typename std::vector<std::pair<Point_3, Vector_3> >::iterator> > Zip_iterator_2;
parallel_callback.join();
CGAL::for_each<ConcurrencyTag>
(CGAL::make_range (boost::make_zip_iterator (boost::make_tuple
(points.begin(), pwns_neighbors.begin(), update_pwns.begin())),
boost::make_zip_iterator (boost::make_tuple
(points.end(), pwns_neighbors.end(), update_pwns.end()))),
[&](const typename Zip_iterator_2::reference& t)
{
if (callback_wrapper.interrupted())
return false;
// If interrupted during this step, nothing is computed, we return NaN
if (parallel_callback.interrupted())
return std::numeric_limits<double>::quiet_NaN();
}
else
#endif // CGAL_LINKED_WITH_TBB
{
std::size_t nb = nb_points;
get<2>(t) = bilateral_smooth_point_set_internal::
compute_denoise_projection<Kernel, PointRange>
(get<0>(t),
point_map, normal_map,
get<1>(t),
guess_neighbor_radius,
sharpness_angle);
++ callback_wrapper.advancement();
return true;
});
callback_wrapper.join();
typename std::vector<Pwn,CGAL_PSP3_DEFAULT_ALLOCATOR<Pwn> >::iterator
update_iter = update_pwns.begin();
typename std::vector<Pwns,CGAL_PSP3_DEFAULT_ALLOCATOR<Pwns> >::iterator
neighbor_iter = pwns_neighbors.begin();
for(pwn_iter = pwns.begin(); pwn_iter != pwns.end();
++pwn_iter, ++update_iter, ++neighbor_iter, ++ nb)
{
*update_iter = bilateral_smooth_point_set_internal::
compute_denoise_projection<Kernel>
(*pwn_iter,
*neighbor_iter,
guess_neighbor_radius,
sharpness_angle);
if (callback && !callback ((nb+1) / double(2. * nb_points)))
return std::numeric_limits<double>::quiet_NaN();
}
}
// If interrupted during this step, nothing is computed, we return NaN
if (callback_wrapper.interrupted())
return std::numeric_limits<double>::quiet_NaN();
#ifdef CGAL_PSP3_VERBOSE
task_timer.stop();
task_timer.stop();
memory = CGAL::Memory_sizer().virtual_size();
std::cout << "done: " << task_timer.time() << " seconds, "
<< (memory>>20) << " Mb allocated" << std::endl;
#endif
// save results
FT sum_move_error = 0;
typename PointRange::iterator it = points.begin();
for(unsigned int i = 0 ; it != points.end(); ++it, ++i)
std::size_t nb = 0;
for (value_type& vt : points)
{
typename boost::property_traits<PointMap>::reference p = get(point_map, *it);
sum_move_error += CGAL::squared_distance(p, update_pwns[i].position());
put (point_map, *it, update_pwns[i].position());
put (normal_map, *it, update_pwns[i].normal());
sum_move_error += CGAL::squared_distance(get(point_map, vt), update_pwns[nb].first);
put (point_map, vt, update_pwns[nb].first);
put (normal_map, vt, update_pwns[nb].second);
++ nb;
}
return sum_move_error / nb_points;
}
/// \cond SKIP_IN_MANUAL
// variant with default NP
// variant with default NP
template <typename ConcurrencyTag,
typename PointRange>
double

170
Point_set_processing_3/include/CGAL/compute_average_spacing.h
View File

@ -18,7 +18,9 @@
#include <CGAL/Search_traits_3.h>
#include <CGAL/squared_distance_3.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Point_set_processing_3/internal/Neighbor_query.h>
#include <CGAL/Point_set_processing_3/internal/Callback_wrapper.h>
#include <CGAL/for_each.h>
#include <CGAL/property_map.h>
#include <CGAL/point_set_processing_assertions.h>
#include <CGAL/assertions.h>
@ -27,15 +29,12 @@
#include <CGAL/boost/graph/Named_function_parameters.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <boost/iterator/zip_iterator.hpp>
#include <iterator>
#include <list>
#ifdef CGAL_LINKED_WITH_TBB
#include <CGAL/Point_set_processing_3/internal/Parallel_callback.h>
#include <tbb/parallel_for.h>
#include <tbb/blocked_range.h>
#include <tbb/scalable_allocator.h>
#endif // CGAL_LINKED_WITH_TBB
#ifdef DOXYGEN_RUNNING
#define CGAL_BGL_NP_TEMPLATE_PARAMETERS NamedParameters
@ -60,81 +59,37 @@ namespace internal {
/// @tparam Tree KD-tree.
///
/// @return average spacing (scalar).
template < typename Kernel,
typename Tree >
typename Kernel::FT
compute_average_spacing(const typename Kernel::Point_3& query, ///< 3D point whose spacing we want to compute
const Tree& tree, ///< KD-tree
template <typename NeighborQuery>
typename NeighborQuery::Kernel::FT
compute_average_spacing(const typename NeighborQuery::Kernel::Point_3& query, ///< 3D point whose spacing we want to compute
const NeighborQuery& neighbor_query, ///< KD-tree
unsigned int k) ///< number of neighbors
{
// basic geometric types
typedef typename NeighborQuery::Kernel Kernel;
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point;
// types for K nearest neighbors search
typedef Search_traits_3<Kernel> Tree_traits;
typedef Orthogonal_k_neighbor_search<Tree_traits> Neighbor_search;
typedef typename Neighbor_search::iterator Search_iterator;
// performs k + 1 queries (if unique the query point is
// output first). search may be aborted when k is greater
// than number of input points
Neighbor_search search(tree,query,k+1);
Search_iterator search_iterator = search.begin();
FT sum_distances = (FT)0.0;
unsigned int i;
for(i=0;i<(k+1);i++)
{
if(search_iterator == search.end())
break; // premature ending
Point p = search_iterator->first;
sum_distances += std::sqrt(CGAL::squared_distance(query,p));
search_iterator++;
}
unsigned int i = 0;
neighbor_query.get_points
(query, k, 0,
boost::make_function_output_iterator
([&](const Point& p)
{
sum_distances += std::sqrt(CGAL::squared_distance (query,p));
++ i;
}));
// output average spacing
return sum_distances / (FT)i;
}
#ifdef CGAL_LINKED_WITH_TBB
template <typename Kernel, typename Tree>
class Compute_average_spacings {
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::FT FT;
const Tree& tree;
const unsigned int k;
const std::vector<Point>& input;
std::vector<FT>& output;
cpp11::atomic<std::size_t>& advancement;
cpp11::atomic<bool>& interrupted;
public:
Compute_average_spacings(Tree& tree, unsigned int k, std::vector<Point>& points,
std::vector<FT>& output,
cpp11::atomic<std::size_t>& advancement,
cpp11::atomic<bool>& interrupted)
: tree(tree), k (k), input (points), output (output)
, advancement (advancement)
, interrupted (interrupted)
{ }
void operator()(const tbb::blocked_range<std::size_t>& r) const
{
for( std::size_t i = r.begin(); i != r.end(); ++i)
{
if (interrupted)
break;
output[i] = CGAL::internal::compute_average_spacing<Kernel,Tree>(input[i], tree, k);
++ advancement;
}
}
};
#endif // CGAL_LINKED_WITH_TBB
} /* namespace internal */
/// \endcond
@ -178,7 +133,7 @@ compute_average_spacing(const typename Kernel::Point_3& query, ///< 3D point who
of `points`.
*/
template <typename ConcurrencyTag,
typename PointRange,
typename PointRange,
typename CGAL_BGL_NP_TEMPLATE_PARAMETERS
>
#ifdef DOXYGEN_RUNNING
@ -195,20 +150,17 @@ compute_average_spacing(
using parameters::get_parameter;
// basic geometric types
typedef typename PointRange::const_iterator iterator;
typedef typename CGAL::GetPointMap<PointRange, CGAL_BGL_NP_CLASS>::const_type PointMap;
typedef typename Point_set_processing_3::GetK<PointRange, CGAL_BGL_NP_CLASS>::Kernel Kernel;
typedef typename Kernel::Point_3 Point;
PointMap point_map = choose_parameter<PointMap>(get_parameter(np, internal_np::point_map));
PointMap point_map = choose_parameter(get_parameter(np, internal_np::point_map), PointMap());
const std::function<bool(double)>& callback = choose_parameter(get_parameter(np, internal_np::callback),
std::function<bool(double)>());
// types for K nearest neighbors search structure
typedef typename Kernel::FT FT;
typedef Search_traits_3<Kernel> Tree_traits;
typedef Orthogonal_k_neighbor_search<Tree_traits> Neighbor_search;
typedef typename Neighbor_search::Tree Tree;
typedef Point_set_processing_3::internal::Neighbor_query<Kernel, const PointRange&, PointMap> Neighbor_query;
// precondition: at least one element in the container.
// to fix: should have at least three distinct points
@ -219,60 +171,46 @@ compute_average_spacing(
CGAL_point_set_processing_precondition(k >= 2);
// Instanciate a KD-tree search.
// Note: We have to convert each input iterator to Point_3.
std::vector<Point> kd_tree_points;
for(typename PointRange::const_iterator it = points.begin(); it != points.end(); it++)
kd_tree_points.push_back(get(point_map, *it));
Tree tree(kd_tree_points.begin(), kd_tree_points.end());
Neighbor_query neighbor_query (points, point_map);
// iterate over input points, compute and output normal
// vectors (already normalized)
FT sum_spacings = (FT)0.0;
std::size_t nb = 0;
std::size_t nb_points = std::distance(points.begin(), points.end());
#ifndef CGAL_LINKED_WITH_TBB
CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
"Parallel_tag is enabled but TBB is unavailable.");
#else
if (boost::is_convertible<ConcurrencyTag,Parallel_tag>::value)
{
Point_set_processing_3::internal::Parallel_callback
parallel_callback (callback, kd_tree_points.size());
Point_set_processing_3::internal::Callback_wrapper<ConcurrencyTag>
callback_wrapper (callback, nb_points);
std::vector<FT> spacings (nb_points, -1);
std::vector<FT> spacings (kd_tree_points.size (), -1);
CGAL::internal::Compute_average_spacings<Kernel, Tree>
f (tree, k, kd_tree_points, spacings,
parallel_callback.advancement(),
parallel_callback.interrupted());
tbb::parallel_for(tbb::blocked_range<size_t>(0, kd_tree_points.size ()), f);
for (unsigned int i = 0; i < spacings.size (); ++ i)
if (spacings[i] >= 0.)
{
sum_spacings += spacings[i];
++ nb;
}
parallel_callback.join();
}
else
#endif
typedef boost::zip_iterator<boost::tuple<iterator, typename std::vector<FT>::iterator> > Zip_iterator;
CGAL::for_each<ConcurrencyTag>
(CGAL::make_range (boost::make_zip_iterator (boost::make_tuple (points.begin(), spacings.begin())),
boost::make_zip_iterator (boost::make_tuple (points.end(), spacings.end()))),
[&](const typename Zip_iterator::reference& t)
{
for(typename PointRange::const_iterator it = points.begin(); it != points.end(); it++, nb++)
{
sum_spacings += internal::compute_average_spacing<Kernel,Tree>(
get(point_map,*it),
tree,k);
if (callback && !callback ((nb+1) / double(kd_tree_points.size())))
{
++ nb;
break;
}
}
}
if (callback_wrapper.interrupted())
return false;
get<1>(t) = CGAL::internal::compute_average_spacing<Neighbor_query>
(get(point_map, get<0>(t)), neighbor_query, k);
++ callback_wrapper.advancement();
return true;
});
for (unsigned int i = 0; i < spacings.size (); ++ i)
if (spacings[i] >= 0.)
{
sum_spacings += spacings[i];
++ nb;
}
callback_wrapper.join();
// return average spacing
return sum_spacings / (FT)(nb);
return sum_spacings / (FT)(nb);
}
/// \cond SKIP_IN_MANUAL

156
Point_set_processing_3/include/CGAL/jet_estimate_normals.h
View File

@ -17,9 +17,9 @@
#include <CGAL/disable_warnings.h>
#include <CGAL/IO/trace.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Point_set_processing_3/internal/neighbor_query.h>
#include <CGAL/Point_set_processing_3/internal/Neighbor_query.h>
#include <CGAL/Point_set_processing_3/internal/Callback_wrapper.h>
#include <CGAL/for_each.h>
#include <CGAL/Monge_via_jet_fitting.h>
#include <CGAL/property_map.h>
#include <CGAL/point_set_processing_assertions.h>
@ -32,13 +32,6 @@
#include <iterator>
#include <list>
#ifdef CGAL_LINKED_WITH_TBB
#include <CGAL/Point_set_processing_3/internal/Parallel_callback.h>
#include <tbb/parallel_for.h>
#include <tbb/blocked_range.h>
#include <tbb/scalable_allocator.h>
#endif // CGAL_LINKED_WITH_TBB
namespace CGAL {
@ -58,18 +51,16 @@ namespace internal {
/// @tparam Tree KD-tree.
///
/// @return Computed normal. Orientation is random.
template < typename Kernel,
typename SvdTraits,
typename Tree
>
typename Kernel::Vector_3
jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute the normal at
Tree& tree, ///< KD-tree
template <typename SvdTraits, typename NeighborQuery>
typename NeighborQuery::Kernel::Vector_3
jet_estimate_normal(const typename NeighborQuery::Point_3& query, ///< point to compute the normal at
const NeighborQuery& neighbor_query, ///< KD-tree
unsigned int k, ///< number of neighbors
typename Kernel::FT neighbor_radius,
typename NeighborQuery::FT neighbor_radius,
unsigned int degree_fitting)
{
// basic geometric types
typedef typename NeighborQuery::Kernel Kernel;
typedef typename Kernel::Point_3 Point;
// types for jet fitting
@ -79,9 +70,8 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
typedef typename Monge_jet_fitting::Monge_form Monge_form;
std::vector<Point> points;
CGAL::Point_set_processing_3::internal::neighbor_query
(query, tree, k, neighbor_radius, points);
neighbor_query.get_points (query, k, neighbor_radius, std::back_inserter(points));
// performs jet fitting
Monge_jet_fitting monge_fit;
const unsigned int degree_monge = 1; // we seek for normal and not more.
@ -92,49 +82,6 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
return monge_form.normal_direction();
}
#ifdef CGAL_LINKED_WITH_TBB
template <typename Kernel, typename SvdTraits, typename Tree>
class Jet_estimate_normals {
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::Vector_3 Vector;
const Tree& tree;
const unsigned int k;
const FT neighbor_radius;
const unsigned int degree_fitting;
const std::vector<Point>& input;
std::vector<Vector>& output;
cpp11::atomic<std::size_t>& advancement;
cpp11::atomic<bool>& interrupted;
public:
Jet_estimate_normals(Tree& tree, unsigned int k, FT neighbor_radius,
std::vector<Point>& points,
unsigned int degree_fitting, std::vector<Vector>& output,
cpp11::atomic<std::size_t>& advancement,
cpp11::atomic<bool>& interrupted)
: tree(tree), k (k), neighbor_radius (neighbor_radius)
, degree_fitting (degree_fitting), input (points), output (output)
, advancement (advancement)
, interrupted (interrupted)
{ }
void operator()(const tbb::blocked_range<std::size_t>& r) const
{
for( std::size_t i = r.begin(); i != r.end(); ++i)
{
if (interrupted)
break;
output[i] = CGAL::internal::jet_estimate_normal<Kernel,SvdTraits>(input[i], tree, k, neighbor_radius, degree_fitting);
++ advancement;
}
}
};
#endif // CGAL_LINKED_WITH_TBB
} /* namespace internal */
/// \endcond
@ -187,7 +134,7 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
\cgalNamedParamsEnd
*/
template <typename ConcurrencyTag,
typename PointRange,
typename PointRange,
typename NamedParameters
>
void
@ -198,10 +145,12 @@ jet_estimate_normals(
{
using parameters::choose_parameter;
using parameters::get_parameter;
CGAL_TRACE("Calls jet_estimate_normals()\n");
// basic geometric types
typedef typename PointRange::iterator iterator;
typedef typename iterator::value_type value_type;
typedef typename CGAL::GetPointMap<PointRange, NamedParameters>::type PointMap;
typedef typename Point_set_processing_3::GetNormalMap<PointRange, NamedParameters>::type NormalMap;
typedef typename Point_set_processing_3::GetK<PointRange, NamedParameters>::Kernel Kernel;
@ -219,19 +168,12 @@ jet_estimate_normals(
NormalMap normal_map = choose_parameter<NormalMap>(get_parameter(np, internal_np::normal_map));
unsigned int degree_fitting = choose_parameter(get_parameter(np, internal_np::degree_fitting), 2);
FT neighbor_radius = choose_parameter(get_parameter(np, internal_np::neighbor_radius), FT(0));
const std::function<bool(double)>& callback = choose_parameter(get_parameter(np, internal_np::callback),
std::function<bool(double)>());
typedef typename Kernel::Point_3 Point;
// Input points types
typedef typename boost::property_traits<NormalMap>::value_type Vector;
// types for K nearest neighbors search structure
typedef typename CGAL::Search_traits_3<Kernel> Tree_traits;
typedef typename CGAL::Orthogonal_k_neighbor_search<Tree_traits> Neighbor_search;
typedef typename Neighbor_search::Tree Tree;
typedef Point_set_processing_3::internal::Neighbor_query<Kernel, PointRange&, PointMap> Neighbor_query;
// precondition: at least one element in the container.
// to fix: should have at least three distinct points
@ -244,60 +186,32 @@ jet_estimate_normals(
std::size_t memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20);
CGAL_TRACE(" Creates KD-tree\n");
typename PointRange::iterator it;
// Instanciate a KD-tree search.
// Note: We have to convert each input iterator to Point_3.
std::vector<Point> kd_tree_points;
for(it = points.begin(); it != points.end(); it++)
kd_tree_points.push_back(get(point_map, *it));
Tree tree(kd_tree_points.begin(), kd_tree_points.end());
Neighbor_query neighbor_query (points, point_map);
memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20);
CGAL_TRACE(" Computes normals\n");
// iterate over input points, compute and output normal
// vectors (already normalized)
#ifndef CGAL_LINKED_WITH_TBB
CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
"Parallel_tag is enabled but TBB is unavailable.");
#else
if (boost::is_convertible<ConcurrencyTag,Parallel_tag>::value)
{
Point_set_processing_3::internal::Parallel_callback
parallel_callback (callback, kd_tree_points.size());
std::size_t nb_points = points.size();
std::vector<Vector> normals (kd_tree_points.size (),
CGAL::NULL_VECTOR);
CGAL::internal::Jet_estimate_normals<Kernel, SvdTraits, Tree>
f (tree, k, neighbor_radius,
kd_tree_points, degree_fitting, normals,
parallel_callback.advancement(),
parallel_callback.interrupted());
tbb::parallel_for(tbb::blocked_range<size_t>(0, kd_tree_points.size ()), f);
std::size_t i = 0;
for(it = points.begin(); it != points.end(); ++ it, ++ i)
if (normals[i] != CGAL::NULL_VECTOR)
put (normal_map, *it, normals[i]);
Point_set_processing_3::internal::Callback_wrapper<ConcurrencyTag>
callback_wrapper (callback, nb_points);
parallel_callback.join();
}
else
#endif
CGAL::for_each<ConcurrencyTag>
(points,
[&](value_type& vt)
{
std::size_t nb = 0;
for(it = points.begin(); it != points.end(); it++, ++ nb)
{
Vector normal = internal::jet_estimate_normal<Kernel,SvdTraits,Tree>(
get(point_map,*it),
tree, k, neighbor_radius, degree_fitting);
if (callback_wrapper.interrupted())
return false;
put (normal_map, vt,
CGAL::internal::jet_estimate_normal<SvdTraits>
(get(point_map, vt), neighbor_query, k, neighbor_radius, degree_fitting));
++ callback_wrapper.advancement();
put(normal_map, *it, normal); // normal_map[it] = normal
if (callback && !callback ((nb+1) / double(kd_tree_points.size())))
break;
}
}
return true;
});
callback_wrapper.join();
memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20);
CGAL_TRACE("End of jet_estimate_normals()\n");
@ -307,7 +221,7 @@ jet_estimate_normals(
/// \cond SKIP_IN_MANUAL
// variant with default NP
template <typename ConcurrencyTag,
typename PointRange>
typename PointRange>
void
jet_estimate_normals(
PointRange& points,

157
Point_set_processing_3/include/CGAL/jet_smooth_point_set.h
View File

@ -17,9 +17,9 @@
#include <CGAL/disable_warnings.h>
#include <CGAL/IO/trace.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Point_set_processing_3/internal/neighbor_query.h>
#include <CGAL/Point_set_processing_3/internal/Neighbor_query.h>
#include <CGAL/Point_set_processing_3/internal/Callback_wrapper.h>
#include <CGAL/for_each.h>
#include <CGAL/Monge_via_jet_fitting.h>
#include <CGAL/property_map.h>
#include <CGAL/point_set_processing_assertions.h>
@ -31,13 +31,6 @@
#include <iterator>
#include <list>
#ifdef CGAL_LINKED_WITH_TBB
#include <CGAL/Point_set_processing_3/internal/Parallel_callback.h>
#include <tbb/parallel_for.h>
#include <tbb/blocked_range.h>
#include <tbb/scalable_allocator.h>
#endif // CGAL_LINKED_WITH_TBB
namespace CGAL {
@ -57,20 +50,20 @@ namespace internal {
/// @tparam Tree KD-tree.
///
/// @return computed point
template <typename Kernel,
typename SvdTraits,
typename Tree
template <typename SvdTraits,
typename NeighborQuery
>
typename Kernel::Point_3
typename NeighborQuery::Kernel::Point_3
jet_smooth_point(
const typename Kernel::Point_3& query, ///< 3D point to project
Tree& tree, ///< KD-tree
const typename NeighborQuery::Kernel::Point_3& query, ///< 3D point to project
NeighborQuery& neighbor_query, ///< KD-tree
const unsigned int k, ///< number of neighbors.
typename Kernel::FT neighbor_radius,
typename NeighborQuery::Kernel::FT neighbor_radius,
const unsigned int degree_fitting,
const unsigned int degree_monge)
{
// basic geometric types
typedef typename NeighborQuery::Kernel Kernel;
typedef typename Kernel::Point_3 Point;
// types for jet fitting
@ -78,10 +71,9 @@ jet_smooth_point(
Simple_cartesian<double>,
SvdTraits> Monge_jet_fitting;
typedef typename Monge_jet_fitting::Monge_form Monge_form;
std::vector<Point> points;
CGAL::Point_set_processing_3::internal::neighbor_query
(query, tree, k, neighbor_radius, points);
neighbor_query.get_points (query, k, neighbor_radius, std::back_inserter(points));
// performs jet fitting
Monge_jet_fitting monge_fit;
@ -92,50 +84,6 @@ jet_smooth_point(
return monge_form.origin();
}
#ifdef CGAL_LINKED_WITH_TBB
template <typename Kernel, typename SvdTraits, typename Tree>
class Jet_smooth_pwns {
typedef typename Kernel::Point_3 Point;
const Tree& tree;
const unsigned int k;
const typename Kernel::FT neighbor_radius;
unsigned int degree_fitting;
unsigned int degree_monge;
const std::vector<Point>& input;
std::vector<Point>& output;
cpp11::atomic<std::size_t>& advancement;
cpp11::atomic<bool>& interrupted;
public:
Jet_smooth_pwns (Tree& tree, unsigned int k, typename Kernel::FT neighbor_radius,
std::vector<Point>& points,
unsigned int degree_fitting, unsigned int degree_monge, std::vector<Point>& output,
cpp11::atomic<std::size_t>& advancement,
cpp11::atomic<bool>& interrupted)
: tree(tree), k (k), neighbor_radius(neighbor_radius)
, degree_fitting (degree_fitting)
, degree_monge (degree_monge), input (points), output (output)
, advancement (advancement)
, interrupted (interrupted)
{ }
void operator()(const tbb::blocked_range<std::size_t>& r) const
{
for( std::size_t i = r.begin(); i != r.end(); ++i)
{
if (interrupted)
break;
output[i] = CGAL::internal::jet_smooth_point<Kernel, SvdTraits>(input[i], tree, k,
neighbor_radius,
degree_fitting,
degree_monge);
++ advancement;
}
}
};
#endif // CGAL_LINKED_WITH_TBB
} /* namespace internal */
@ -191,7 +139,7 @@ jet_smooth_point(
*/
template <typename ConcurrencyTag,
typename PointRange,
typename PointRange,
typename NamedParameters
>
void
@ -202,8 +150,10 @@ jet_smooth_point_set(
{
using parameters::choose_parameter;
using parameters::get_parameter;
// basic geometric types
typedef typename PointRange::iterator iterator;
typedef typename iterator::value_type value_type;
typedef typename CGAL::GetPointMap<PointRange, NamedParameters>::type PointMap;
typedef typename Point_set_processing_3::GetK<PointRange, NamedParameters>::Kernel Kernel;
typedef typename GetSvdTraits<NamedParameters>::type SvdTraits;
@ -220,12 +170,8 @@ jet_smooth_point_set(
const std::function<bool(double)>& callback = choose_parameter(get_parameter(np, internal_np::callback),
std::function<bool(double)>());
typedef typename Kernel::Point_3 Point;
// types for K nearest neighbors search structure
typedef typename CGAL::Search_traits_3<Kernel> Tree_traits;
typedef typename CGAL::Orthogonal_k_neighbor_search<Tree_traits> Neighbor_search;
typedef typename Neighbor_search::Tree Tree;
typedef Point_set_processing_3::internal::Neighbor_query<Kernel, PointRange&, PointMap> Neighbor_query;
// precondition: at least one element in the container.
// to fix: should have at least three distinct points
@ -235,63 +181,44 @@ jet_smooth_point_set(
// precondition: at least 2 nearest neighbors
CGAL_point_set_processing_precondition(k >= 2);
typename PointRange::iterator it;
// Instanciate a KD-tree search.
// Note: We have to convert each input iterator to Point_3.
std::vector<Point> kd_tree_points;
for(it = points.begin(); it != points.end(); it++)
kd_tree_points.push_back(get(point_map, *it));
Tree tree(kd_tree_points.begin(), kd_tree_points.end());
Neighbor_query neighbor_query (points, point_map);
// Iterates over input points and mutates them.
// Implementation note: the cast to Point& allows to modify only the point's position.
#ifndef CGAL_LINKED_WITH_TBB
CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
"Parallel_tag is enabled but TBB is unavailable.");
#else
if (boost::is_convertible<ConcurrencyTag,Parallel_tag>::value)
{
Point_set_processing_3::internal::Parallel_callback
parallel_callback (callback, kd_tree_points.size());
std::size_t nb_points = points.size();
std::vector<Point> mutated_points (kd_tree_points.size (), CGAL::ORIGIN);
CGAL::internal::Jet_smooth_pwns<Kernel, SvdTraits, Tree>
f (tree, k, neighbor_radius, kd_tree_points, degree_fitting, degree_monge,
mutated_points,
parallel_callback.advancement(),
parallel_callback.interrupted());
tbb::parallel_for(tbb::blocked_range<size_t>(0, kd_tree_points.size ()), f);
unsigned int i = 0;
for(it = points.begin(); it != points.end(); ++ it, ++ i)
if (mutated_points[i] != CGAL::ORIGIN)
put(point_map, *it, mutated_points[i]);
Point_set_processing_3::internal::Callback_wrapper<ConcurrencyTag>
callback_wrapper (callback, nb_points);
parallel_callback.join();
}
else
#endif
CGAL::for_each<ConcurrencyTag>
(points,
[&](value_type vt)
{
std::size_t nb = 0;
for(it = points.begin(); it != points.end(); it++, ++ nb)
{
const typename boost::property_traits<PointMap>::reference p = get(point_map, *it);
put(point_map, *it ,
internal::jet_smooth_point<Kernel, SvdTraits>(
p,tree,k,neighbor_radius,degree_fitting,degree_monge) );
if (callback && !callback ((nb+1) / double(kd_tree_points.size())))
break;
}
}
if (callback_wrapper.interrupted())
return false;
put (point_map, vt,
CGAL::internal::jet_smooth_point<SvdTraits>
(get (point_map, vt), neighbor_query,
k,
neighbor_radius,
degree_fitting,
degree_monge));
++ callback_wrapper.advancement();
return true;
});
callback_wrapper.join();
}
/// \cond SKIP_IN_MANUAL
// variant with default NP
template <typename ConcurrencyTag,
typename PointRange>
typename PointRange>
void
jet_smooth_point_set(
PointRange& points,
@ -301,7 +228,7 @@ jet_smooth_point_set(
(points, k, CGAL::Point_set_processing_3::parameters::all_default(points));
}
/// \endcond
} //namespace CGAL

82
Point_set_processing_3/include/CGAL/mst_orient_normals.h
View File

@ -17,10 +17,7 @@
#include <CGAL/disable_warnings.h>
#include <CGAL/IO/trace.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Point_set_processing_3/internal/neighbor_query.h>
#include <CGAL/Point_set_processing_3/internal/Search_traits_vertex_handle_3.h>
#include <CGAL/Point_set_processing_3/internal/Neighbor_query.h>
#include <CGAL/property_map.h>
#include <CGAL/Index_property_map.h>
#include <CGAL/Memory_sizer.h>
@ -125,7 +122,7 @@ public:
typedef value_type reference;
private:
ForwardIterator m_source_point;
public:
@ -188,7 +185,7 @@ struct Propagate_normal_orientation
// Gets target
vertex_descriptor target_vertex = target(edge, mst_graph);
bool& target_normal_is_oriented = ((MST_graph&)mst_graph)[target_vertex].is_oriented;
// special case if vertex is source vertex (and thus has no related point/normal)
if (source_vertex == m_source)
{
@ -202,7 +199,7 @@ struct Propagate_normal_orientation
// Gets target
Vector_ref target_normal = get( mst_graph.m_normal_map, *(mst_graph[target_vertex].input_point) );
if ( ! target_normal_is_oriented )
{
// -> ->
@ -261,10 +258,10 @@ mst_find_source(
ForwardIterator top_point = first;
for (ForwardIterator v = ++first; v != beyond; v++)
{
double top_z = get(point_map,*top_point).z(); // top_point's Z coordinate
double z = get(point_map,*v).z();
if (top_z < z)
top_point = v;
}
@ -296,17 +293,16 @@ mst_find_source(
/// @tparam Kernel Geometric traits class.
///
/// @return the Riemannian graph
template <typename ForwardIterator,
template <typename PointRange,
typename PointMap,
typename NormalMap,
typename IndexMap,
typename ConstrainedMap,
typename Kernel
>
Riemannian_graph<ForwardIterator>
Riemannian_graph<typename PointRange::iterator>
create_riemannian_graph(
ForwardIterator first, ///< iterator over the first input point.
ForwardIterator beyond, ///< past-the-end iterator over the input points.
PointRange& points, ///< input points
PointMap point_map, ///< property map: value_type of ForwardIterator -> Point_3
NormalMap normal_map, ///< property map: value_type of ForwardIterator -> Vector_3
IndexMap index_map, ///< property map ForwardIterator -> index
@ -320,43 +316,23 @@ create_riemannian_graph(
typedef typename boost::property_traits<NormalMap>::reference Vector_ref;
// Types for K nearest neighbors search structure
typedef Point_vertex_handle_3<ForwardIterator> Point_vertex_handle_3;
typedef internal::Search_traits_vertex_handle_3<ForwardIterator> Traits;
typedef Euclidean_distance_vertex_handle_3<ForwardIterator> KDistance;
typedef Orthogonal_k_neighbor_search<Traits,KDistance> Neighbor_search;
typedef typename Neighbor_search::Tree Tree;
typedef typename PointRange::iterator ForwardIterator;
typedef Point_set_processing_3::internal::Neighbor_query<Kernel, PointRange&, PointMap> Neighbor_query;
// Riemannian_graph types
typedef internal::Riemannian_graph<ForwardIterator> Riemannian_graph;
typedef typename boost::property_map<Riemannian_graph, boost::edge_weight_t>::type Riemannian_graph_weight_map;
// Precondition: at least one element in the container.
CGAL_point_set_processing_precondition(first != beyond);
// Precondition: at least 2 nearest neighbors
CGAL_point_set_processing_precondition(k >= 2);
// Number of input points
const std::size_t num_input_points = distance(first, beyond);
const std::size_t num_input_points = points.size();
std::size_t memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20);
CGAL_TRACE(" Creates KD-tree\n");
// Instanciate a KD-tree search.
// Notes: We have to wrap each input point by a Point_vertex_handle_3.
// The KD-tree is allocated dynamically to recover RAM as soon as possible.
std::vector<Point_vertex_handle_3> kd_tree_points; kd_tree_points.reserve(num_input_points);
for (ForwardIterator it = first; it != beyond; it++)
{
Point_ref point = get(point_map, *it);
Point_vertex_handle_3 point_wrapper(point.x(), point.y(), point.z(), it);
kd_tree_points.push_back(point_wrapper);
}
boost::shared_ptr<Tree> tree( new Tree(kd_tree_points.begin(), kd_tree_points.end()) );
// Recover RAM
kd_tree_points.clear();
Neighbor_query neighbor_query (points, point_map);
memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20);
CGAL_TRACE(" Creates Riemannian Graph\n");
@ -369,7 +345,7 @@ create_riemannian_graph(
Riemannian_graph riemannian_graph;
//
// add vertices
for (ForwardIterator it = first; it != beyond; it++)
for (ForwardIterator it = points.begin(); it != points.end(); it++)
{
typename Riemannian_graph::vertex_descriptor v = add_vertex(riemannian_graph);
CGAL_point_set_processing_assertion(v == get(index_map,it));
@ -383,16 +359,14 @@ create_riemannian_graph(
//
// add edges
Riemannian_graph_weight_map riemannian_graph_weight_map = get(boost::edge_weight, riemannian_graph);
for (ForwardIterator it = first; it != beyond; it++)
for (ForwardIterator it = points.begin(); it != points.end(); it++)
{
std::size_t it_index = get(index_map,it);
Vector_ref it_normal_vector = get(normal_map,*it);
Point_ref point = get(point_map, *it);
Point_vertex_handle_3 point_wrapper(point.x(), point.y(), point.z(), it);
std::vector<Point_vertex_handle_3> neighbor_points;
CGAL::Point_set_processing_3::internal::neighbor_query
(point_wrapper, *tree, k, neighbor_radius, neighbor_points);
std::vector<ForwardIterator> neighbor_points;
neighbor_query.get_iterators (point, k, neighbor_radius, std::back_inserter(neighbor_points));
for (std::size_t i = 0; i < neighbor_points.size(); ++ i)
{
@ -411,7 +385,7 @@ create_riemannian_graph(
// -> ->
// Computes edge weight = 1 - | normal1 * normal2 |
// where normal1 and normal2 are the normal at the edge extremities.
Vector_ref neighbor_normal_vector = get(normal_map,*neighbor);
double weight = 1.0 - std::abs(it_normal_vector * neighbor_normal_vector);
if (weight < 0)
@ -432,7 +406,7 @@ create_riemannian_graph(
riemannian_graph_weight_map[e] = 0.;
}
}
return riemannian_graph;
@ -496,7 +470,7 @@ create_mst_graph(
// Computes Minimum Spanning Tree.
std::size_t source_point_index = num_input_points;
Riemannian_graph_weight_map riemannian_graph_weight_map = get(boost::edge_weight, riemannian_graph);
typedef std::vector<typename Riemannian_graph::vertex_descriptor> PredecessorMap;
PredecessorMap predecessor(num_input_points + 1);
@ -524,11 +498,11 @@ create_mst_graph(
mst_graph[v].input_point = it;
mst_graph[v].is_oriented = false;
}
typename MST_graph::vertex_descriptor v = add_vertex(mst_graph);
CGAL_point_set_processing_assertion(v == source_point_index);
mst_graph[v].is_oriented = true;
// add edges
for (std::size_t i=0; i < predecessor.size(); i++) // add edges
{
@ -553,7 +527,7 @@ create_mst_graph(
// Public section
// ----------------------------------------------------------------------------
/**
/**
\ingroup PkgPointSetProcessing3Algorithms
Orients the normals of the range of `points` using the propagation
of a seed orientation through a minimum spanning tree of the Riemannian graph.
@ -588,7 +562,7 @@ create_mst_graph(
limit is wanted, use `k=0`.\cgalParamEnd
\cgalParamBegin{point_is_constrained_map} a model of `ReadablePropertyMap` with value type
`bool`. Points with a `true` value will be used as seed points: their normal will be considered as already
oriented, it won't be altered and it will be propagated to its neighbors. If this parameter is omitted,
oriented, it won't be altered and it will be propagated to its neighbors. If this parameter is omitted,
the highest point (highest Z coordinate) will be used as the unique seed with an upward oriented
normal\cgalParamEnd
\cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
@ -665,7 +639,7 @@ mst_orient_normals(
if (boost::is_same<ConstrainedMap,
typename CGAL::Point_set_processing_3::GetIsConstrainedMap<PointRange, NamedParameters>::NoMap>::value)
riemannian_graph = create_riemannian_graph(points.begin(), points.end(),
riemannian_graph = create_riemannian_graph(points,
point_map, normal_map, index_map,
Default_constrained_map<typename PointRange::iterator>
(mst_find_source(points.begin(), points.end(),
@ -675,7 +649,7 @@ mst_orient_normals(
neighbor_radius,
kernel);
else
riemannian_graph = create_riemannian_graph(points.begin(), points.end(),
riemannian_graph = create_riemannian_graph(points,
point_map, normal_map, index_map,
constrained_map,
k,
@ -694,7 +668,7 @@ mst_orient_normals(
const std::size_t num_input_points = distance(points.begin(), points.end());
std::size_t source_point_index = num_input_points;
// Traverse the point set along the MST to propagate source_point's orientation
Propagate_normal_orientation<typename PointRange::iterator, NormalMap, Kernel> orienter(source_point_index);

152
Point_set_processing_3/include/CGAL/pca_estimate_normals.h
View File

@ -18,9 +18,9 @@
#include <CGAL/IO/trace.h>
#include <CGAL/Dimension.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Point_set_processing_3/internal/neighbor_query.h>
#include <CGAL/Point_set_processing_3/internal/Neighbor_query.h>
#include <CGAL/Point_set_processing_3/internal/Callback_wrapper.h>
#include <CGAL/for_each.h>
#include <CGAL/linear_least_squares_fitting_3.h>
#include <CGAL/property_map.h>
#include <CGAL/point_set_processing_assertions.h>
@ -33,13 +33,6 @@
#include <iterator>
#include <list>
#ifdef CGAL_LINKED_WITH_TBB
#include <CGAL/Point_set_processing_3/internal/Parallel_callback.h>
#include <tbb/parallel_for.h>
#include <tbb/blocked_range.h>
#include <tbb/scalable_allocator.h>
#endif // CGAL_LINKED_WITH_TBB
namespace CGAL {
@ -58,23 +51,21 @@ namespace internal {
/// @tparam Tree KD-tree.
///
/// @return Computed normal. Orientation is random.
template < typename Kernel,
typename Tree
>
typename Kernel::Vector_3
pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute the normal at
const Tree& tree, ///< KD-tree
template <typename NeighborQuery>
typename NeighborQuery::Kernel::Vector_3
pca_estimate_normal(const typename NeighborQuery::Kernel::Point_3& query, ///< point to compute the normal at
const NeighborQuery& neighbor_query, ///< KD-tree
unsigned int k, ///< number of neighbors
typename Kernel::FT neighbor_radius)
typename NeighborQuery::Kernel::FT neighbor_radius)
{
// basic geometric types
typedef typename NeighborQuery::Kernel Kernel;
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::Plane_3 Plane;
std::vector<Point> points;
CGAL::Point_set_processing_3::internal::neighbor_query
(query, tree, k, neighbor_radius, points);
neighbor_query.get_points (query, k, neighbor_radius, std::back_inserter(points));
// performs plane fitting by point-based PCA
Plane plane;
linear_least_squares_fitting_3(points.begin(),points.end(),plane,Dimension_tag<0>());
@ -83,48 +74,6 @@ pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
return plane.orthogonal_vector();
}
#ifdef CGAL_LINKED_WITH_TBB
template <typename Kernel, typename Tree>
class PCA_estimate_normals {
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::Vector_3 Vector;
const Tree& tree;
const unsigned int k;
const FT neighbor_radius;
const std::vector<Point>& input;
std::vector<Vector>& output;
cpp11::atomic<std::size_t>& advancement;
cpp11::atomic<bool>& interrupted;
public:
PCA_estimate_normals(Tree& tree, unsigned int k, FT neighbor_radius,
std::vector<Point>& points,
std::vector<Vector>& output,
cpp11::atomic<std::size_t>& advancement,
cpp11::atomic<bool>& interrupted)
: tree(tree), k (k), neighbor_radius (neighbor_radius)
, input (points), output (output)
, advancement (advancement)
, interrupted (interrupted)
{ }
void operator()(const tbb::blocked_range<std::size_t>& r) const
{
for( std::size_t i = r.begin(); i != r.end(); ++i)
{
if (interrupted)
break;
output[i] = CGAL::internal::pca_estimate_normal<Kernel,Tree>(input[i], tree, k, neighbor_radius);
++ advancement;
}
}
};
#endif // CGAL_LINKED_WITH_TBB
} /* namespace internal */
/// \endcond
@ -176,7 +125,7 @@ pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
\cgalNamedParamsEnd
*/
template <typename ConcurrencyTag,
typename PointRange,
typename PointRange,
typename NamedParameters
>
void
@ -206,15 +155,12 @@ pca_estimate_normals(
const std::function<bool(double)>& callback = choose_parameter(get_parameter(np, internal_np::callback),
std::function<bool(double)>());
typedef typename Kernel::Point_3 Point;
// Input points types
typedef typename boost::property_traits<NormalMap>::value_type Vector;
typedef typename PointRange::iterator iterator;
typedef typename iterator::value_type value_type;
// types for K nearest neighbors search structure
typedef typename CGAL::Search_traits_3<Kernel> Tree_traits;
typedef typename CGAL::Orthogonal_k_neighbor_search<Tree_traits> Neighbor_search;
typedef typename Neighbor_search::Tree Tree;
typedef Point_set_processing_3::internal::Neighbor_query<Kernel, PointRange&, PointMap> Neighbor_query;
// precondition: at least one element in the container.
// to fix: should have at least three distinct points
@ -227,60 +173,34 @@ pca_estimate_normals(
std::size_t memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20);
CGAL_TRACE(" Creates KD-tree\n");
typename PointRange::iterator it;
// Instanciate a KD-tree search.
// Note: We have to convert each input iterator to Point_3.
std::vector<Point> kd_tree_points;
for(it = points.begin(); it != points.end(); it++)
kd_tree_points.push_back(get(point_map, *it));
Tree tree(kd_tree_points.begin(), kd_tree_points.end());
Neighbor_query neighbor_query (points, point_map);
memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20);
CGAL_TRACE(" Computes normals\n");
// iterate over input points, compute and output normal
// vectors (already normalized)
#ifndef CGAL_LINKED_WITH_TBB
CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
"Parallel_tag is enabled but TBB is unavailable.");
#else
if (boost::is_convertible<ConcurrencyTag,Parallel_tag>::value)
{
Point_set_processing_3::internal::Parallel_callback
parallel_callback (callback, kd_tree_points.size());
std::size_t nb_points = points.size();
std::vector<Vector> normals (kd_tree_points.size (),
CGAL::NULL_VECTOR);
CGAL::internal::PCA_estimate_normals<Kernel, Tree>
f (tree, k, neighbor_radius, kd_tree_points, normals,
parallel_callback.advancement(),
parallel_callback.interrupted());
tbb::parallel_for(tbb::blocked_range<size_t>(0, kd_tree_points.size ()), f);
unsigned int i = 0;
for(it = points.begin(); it != points.end(); ++ it, ++ i)
if (normals[i] != CGAL::NULL_VECTOR)
put (normal_map, *it, normals[i]);
Point_set_processing_3::internal::Callback_wrapper<ConcurrencyTag>
callback_wrapper (callback, nb_points);
parallel_callback.join();
}
else
#endif
{
std::size_t nb = 0;
for(it = points.begin(); it != points.end(); it++, ++ nb)
{
Vector normal = internal::pca_estimate_normal<Kernel,Tree>(
get(point_map,*it),
tree,
k, neighbor_radius);
CGAL::for_each<ConcurrencyTag>
(points,
[&](value_type& vt)
{
if (callback_wrapper.interrupted())
return false;
put (normal_map, vt,
CGAL::internal::pca_estimate_normal
(get(point_map, vt), neighbor_query, k, neighbor_radius));
++ callback_wrapper.advancement();
put(normal_map, *it, normal); // normal_map[it] = normal
if (callback && !callback ((nb+1) / double(kd_tree_points.size())))
break;
}
}
return true;
});
callback_wrapper.join();
memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20);
CGAL_TRACE("End of pca_estimate_normals()\n");
}
@ -288,7 +208,7 @@ pca_estimate_normals(
/// \cond SKIP_IN_MANUAL
// variant with default NP
template <typename ConcurrencyTag,
typename PointRange
typename PointRange
>
void
pca_estimate_normals(

54
Point_set_processing_3/include/CGAL/remove_outliers.h
View File

@ -16,9 +16,7 @@
#include <CGAL/disable_warnings.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Point_set_processing_3/internal/neighbor_query.h>
#include <CGAL/Point_set_processing_3/internal/Neighbor_query.h>
#include <CGAL/property_map.h>
#include <CGAL/point_set_processing_assertions.h>
#include <functional>
@ -49,22 +47,21 @@ namespace internal {
/// @tparam Tree KD-tree.
///
/// @return computed distance.
template < typename Kernel,
typename Tree >
typename Kernel::FT
template <typename NeighborQuery>
typename NeighborQuery::Kernel::FT
compute_avg_knn_sq_distance_3(
const typename Kernel::Point_3& query, ///< 3D point to project
Tree& tree, ///< KD-tree
const typename NeighborQuery::Kernel::Point_3& query, ///< 3D point to project
NeighborQuery& neighbor_query, ///< KD-tree
unsigned int k, ///< number of neighbors
typename Kernel::FT neighbor_radius)
typename NeighborQuery::Kernel::FT neighbor_radius)
{
// geometric types
typedef typename NeighborQuery::Kernel Kernel;
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point;
std::vector<Point> points;
CGAL::Point_set_processing_3::internal::neighbor_query
(query, tree, k, neighbor_radius, points);
neighbor_query.get_points (query, k, neighbor_radius, std::back_inserter(points));
// compute average squared distance
typename Kernel::Compute_squared_distance_3 sqd;
@ -146,7 +143,7 @@ remove_outliers(
{
using parameters::choose_parameter;
using parameters::get_parameter;
// geometric types
typedef typename CGAL::GetPointMap<PointRange, NamedParameters>::type PointMap;
typedef typename Point_set_processing_3::GetK<PointRange, NamedParameters>::Kernel Kernel;
@ -158,19 +155,18 @@ remove_outliers(
double threshold_distance = choose_parameter(get_parameter(np, internal_np::threshold_distance), 0.);
const std::function<bool(double)>& callback = choose_parameter(get_parameter(np, internal_np::callback),
std::function<bool(double)>());
typedef typename Kernel::FT FT;
// basic geometric types
typedef typename Kernel::Point_3 Point;
typedef typename PointRange::iterator iterator;
typedef typename iterator::value_type value_type;
// actual type of input points
typedef typename std::iterator_traits<typename PointRange::iterator>::value_type Enriched_point;
// types for K nearest neighbors search structure
typedef typename CGAL::Search_traits_3<Kernel> Tree_traits;
typedef typename CGAL::Orthogonal_k_neighbor_search<Tree_traits> Neighbor_search;
typedef typename Neighbor_search::Tree Tree;
typedef Point_set_processing_3::internal::Neighbor_query<Kernel, PointRange&, PointMap> Neighbor_query;
// precondition: at least one element in the container.
// to fix: should have at least three distinct points
@ -182,26 +178,22 @@ remove_outliers(
CGAL_point_set_processing_precondition(threshold_percent >= 0 && threshold_percent <= 100);
typename PointRange::iterator it;
Neighbor_query neighbor_query (points, point_map);
// Instanciate a KD-tree search.
// Note: We have to convert each input iterator to Point_3.
std::vector<Point> kd_tree_points;
for(it = points.begin(); it != points.end(); it++)
kd_tree_points.push_back( get(point_map, *it) );
Tree tree(kd_tree_points.begin(), kd_tree_points.end());
std::size_t nb_points = points.size();
// iterate over input points and add them to multimap sorted by distance to k
std::multimap<FT,Enriched_point> sorted_points;
std::size_t nb = 0;
for(it = points.begin(); it != points.end(); it++, ++ nb)
for(const value_type& vt : points)
{
FT sq_distance = internal::compute_avg_knn_sq_distance_3<Kernel>(
get(point_map,*it),
tree, k, neighbor_radius);
sorted_points.insert( std::make_pair(sq_distance, *it) );
if (callback && !callback ((nb+1) / double(kd_tree_points.size())))
FT sq_distance = internal::compute_avg_knn_sq_distance_3(
get(point_map, vt),
neighbor_query, k, neighbor_radius);
sorted_points.insert( std::make_pair(sq_distance, vt) );
if (callback && !callback ((nb+1) / double(nb_points)))
return points.end();
++ nb;
}
// Replaces [points.begin(), points.end()) range by the multimap content.

243
Point_set_processing_3/include/CGAL/wlop_simplify_and_regularize_point_set.h
View File

@ -7,7 +7,7 @@
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Shihao Wu, Clement Jamin, Pierre Alliez
// Author(s) : Shihao Wu, Clement Jamin, Pierre Alliez
#ifndef CGAL_wlop_simplify_and_regularize_point_set_H
#define CGAL_wlop_simplify_and_regularize_point_set_H
@ -32,12 +32,8 @@
#include <cmath>
#include <ctime>
#ifdef CGAL_LINKED_WITH_TBB
#include <CGAL/Point_set_processing_3/internal/Parallel_callback.h>
#include <tbb/parallel_for.h>
#include <tbb/blocked_range.h>
#include <tbb/scalable_allocator.h>
#endif // CGAL_LINKED_WITH_TBB
#include <CGAL/Point_set_processing_3/internal/Callback_wrapper.h>
#include <CGAL/for_each.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Fuzzy_sphere.h>
@ -87,9 +83,9 @@ public:
typedef typename Kernel::Point_3 PointType;
};
/// Compute average and repulsion term, then
/// Compute average and repulsion term, then
/// compute and update sample point locations
///
///
/// \pre `radius > 0`
///
/// @tparam Kernel Geometric traits class.
@ -105,8 +101,8 @@ compute_update_sample_point(
const Tree& original_kd_tree, ///< original Kd-tree
const Tree& sample_kd_tree, ///< sample Kd-tree
const typename Kernel::FT radius, ///< neighborhood radius square
const std::vector<typename Kernel::FT>& original_densities, ///<
const std::vector<typename Kernel::FT>& sample_densities ///<
const std::vector<typename Kernel::FT>& original_densities, ///<
const std::vector<typename Kernel::FT>& sample_densities ///<
)
{
CGAL_point_set_processing_precondition(radius > 0);
@ -130,7 +126,7 @@ compute_update_sample_point(
//Compute average term
FT radius2 = radius * radius;
Vector average = CGAL::NULL_VECTOR;
Vector average = CGAL::NULL_VECTOR;
FT weight = (FT)0.0, average_weight_sum = (FT)0.0;
FT iradius16 = -(FT)4.0 / radius2;
@ -162,10 +158,10 @@ compute_update_sample_point(
}
else
{
average = average / average_weight_sum;
average = average / average_weight_sum;
}
neighbor_original_points.clear();
//Compute repulsion term
@ -175,7 +171,7 @@ compute_update_sample_point(
weight = (FT)0.0;
FT repulsion_weight_sum = (FT)0.0;
Vector repulsion = CGAL::NULL_VECTOR;
Vector repulsion = CGAL::NULL_VECTOR;
iter = neighbor_sample_points.begin();
for(; iter != neighbor_sample_points.end(); ++iter)
@ -188,9 +184,9 @@ compute_update_sample_point(
FT dist2 = CGAL::squared_distance(query, np);
if (dist2 < 1e-10) continue;
FT dist = std::sqrt(dist2);
weight = std::exp(dist2 * iradius16) * std::pow(FT(1.0) / dist, 2); // L1
if (!is_sample_densities_empty)
{
weight *= sample_densities[sample_index];
@ -208,7 +204,7 @@ compute_update_sample_point(
}
else
{
repulsion = repulsion / repulsion_weight_sum;
repulsion = repulsion / repulsion_weight_sum;
}
neighbor_sample_points.clear();
@ -220,7 +216,7 @@ compute_update_sample_point(
/// Compute density weight for each original points,
/// according to their neighbor original points
///
///
/// \pre `k >= 2`, radius > 0
///
/// @tparam Kernel Geometric traits class.
@ -240,7 +236,7 @@ compute_density_weight_for_original_point(
// basic geometric types
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::FT FT;
//types for range search
typedef simplify_and_regularize_internal::Kd_tree_element<Kernel> Kd_tree_point;
typedef simplify_and_regularize_internal::Kd_tree_traits<Kernel> Traits;
@ -266,7 +262,7 @@ compute_density_weight_for_original_point(
FT dist2 = CGAL::squared_distance(query, np);
if (dist2 < 1e-8) continue;
density_weight += std::exp(dist2 * iradius16);
}
@ -277,7 +273,7 @@ compute_density_weight_for_original_point(
/// Compute density weight for sample point,
/// according to their neighbor sample points
///
///
/// \pre `k >= 2`, radius > 0
///
/// @tparam Kernel Geometric traits class.
@ -321,7 +317,7 @@ compute_density_weight_for_sample_point(
FT dist2 = CGAL::squared_distance(query, np);
density_weight += std::exp(dist2 * iradius16);
}
return density_weight;
}
@ -329,67 +325,6 @@ compute_density_weight_for_sample_point(
/// \endcond
#ifdef CGAL_LINKED_WITH_TBB
/// \cond SKIP_IN_MANUAL
/// This is for parallelization of function: compute_denoise_projection()
template <typename Kernel, typename Tree, typename RandomAccessIterator>
class Sample_point_updater
{
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::FT FT;
std::vector<Point> &update_sample_points;
std::vector<Point> &sample_points;
const Tree &original_kd_tree;
const Tree &sample_kd_tree;
const typename Kernel::FT radius;
const std::vector<typename Kernel::FT> &original_densities;
const std::vector<typename Kernel::FT> &sample_densities;
cpp11::atomic<std::size_t>& advancement;
cpp11::atomic<bool>& interrupted;
public:
Sample_point_updater(
std::vector<Point> &out,
std::vector<Point> &in,
const Tree &_original_kd_tree,
const Tree &_sample_kd_tree,
const typename Kernel::FT _radius,
const std::vector<typename Kernel::FT> &_original_densities,
const std::vector<typename Kernel::FT> &_sample_densities,
cpp11::atomic<std::size_t>& advancement,
cpp11::atomic<bool>& interrupted):
update_sample_points(out),
sample_points(in),
original_kd_tree(_original_kd_tree),
sample_kd_tree(_sample_kd_tree),
radius(_radius),
original_densities(_original_densities),
sample_densities(_sample_densities),
advancement (advancement),
interrupted (interrupted) {}
void operator() ( const tbb::blocked_range<size_t>& r ) const
{
for (size_t i = r.begin(); i != r.end(); ++i)
{
if (interrupted)
break;
update_sample_points[i] = simplify_and_regularize_internal::
compute_update_sample_point<Kernel, Tree, RandomAccessIterator>(
sample_points[i],
original_kd_tree,
sample_kd_tree,
radius,
original_densities,
sample_densities);
++ advancement;
}
}
};
/// \endcond
#endif // CGAL_LINKED_WITH_TBB
// ----------------------------------------------------------------------------
// Public section
@ -398,24 +333,24 @@ public:
/**
\ingroup PkgPointSetProcessing3Algorithms
This is an implementation of the Weighted Locally Optimal Projection (WLOP) simplification algorithm.
The WLOP simplification algorithm can produce a set of
denoised, outlier-free and evenly distributed particles over the original
dense point cloud.
The WLOP simplification algorithm can produce a set of
denoised, outlier-free and evenly distributed particles over the original
dense point cloud.
The core of the algorithm is a Weighted Locally Optimal Projection operator
with a density uniformization term.
with a density uniformization term.
For more details, please refer to \cgalCite{wlop-2009}.
A parallel version of WLOP is provided and requires the executable to be
A parallel version of WLOP is provided and requires the executable to be
linked against the <a href="https://www.threadingbuildingblocks.org">Intel TBB library</a>.
To control the number of threads used, the user may use the tbb::task_scheduler_init class.
See the <a href="https://www.threadingbuildingblocks.org/documentation">TBB documentation</a>
See the <a href="https://www.threadingbuildingblocks.org/documentation">TBB documentation</a>
for more details.
\tparam ConcurrencyTag enables sequential versus parallel algorithm. Possible values are `Sequential_tag`,
`Parallel_tag`, and `Parallel_if_available_tag`.
\tparam PointRange is a model of `Range`. The value type of
its iterator is the key type of the named parameter `point_map`.
\tparam OutputIterator Type of the output iterator.
\tparam OutputIterator Type of the output iterator.
It must accept objects of type `geom_traits::Point_3`.
\param points input point range.
@ -427,11 +362,11 @@ public:
If this parameter is omitted, `CGAL::Identity_property_map<geom_traits::Point_3>` is used.\cgalParamEnd
\cgalParamBegin{normal_map} a model of `ReadWritePropertyMap` with value type
`geom_traits::Vector_3`.\cgalParamEnd
\cgalParamBegin{select_percentage} percentage of points to retain. The default value is set to
\cgalParamBegin{select_percentage} percentage of points to retain. The default value is set to
5 (\%).\cgalParamEnd
\cgalParamBegin{neighbor_radius} spherical neighborhood radius. This is a key parameter that needs to be
finely tuned. The result will be irregular if too small, but a larger value will impact the runtime. In
practice, choosing a radius such that the neighborhood of each sample point includes at least two rings
finely tuned. The result will be irregular if too small, but a larger value will impact the runtime. In
practice, choosing a radius such that the neighborhood of each sample point includes at least two rings
of neighboring sample points gives satisfactory result. If this parameter is not provided, it is
automatically set to 8 times the average spacing of the point set.\cgalParamEnd
\cgalParamBegin{number_of_iterations} number of iterations to solve the optimsation problem. The default
@ -462,7 +397,7 @@ wlop_simplify_and_regularize_point_set(
{
using parameters::choose_parameter;
using parameters::get_parameter;
// basic geometric types
typedef typename CGAL::GetPointMap<PointRange, NamedParameters>::type PointMap;
typedef typename Point_set_processing_3::GetK<PointRange, NamedParameters>::Kernel Kernel;
@ -488,15 +423,15 @@ wlop_simplify_and_regularize_point_set(
// to fix: should have at least three distinct points
// but this is costly to check
CGAL_point_set_processing_precondition(points.begin() != points.end());
CGAL_point_set_processing_precondition(select_percentage >= 0
CGAL_point_set_processing_precondition(select_percentage >= 0
&& select_percentage <= 100);
// Random shuffle
CGAL::cpp98::random_shuffle (points.begin(), points.end());
// Computes original(input) and sample points size
// Computes original(input) and sample points size
std::size_t number_of_original = std::distance(points.begin(), points.end());
std::size_t number_of_sample = (std::size_t)(FT(number_of_original) *
std::size_t number_of_sample = (std::size_t)(FT(number_of_original) *
(select_percentage / FT(100.0)));
std::size_t first_index_to_sample = number_of_original - number_of_sample;
@ -509,13 +444,13 @@ wlop_simplify_and_regularize_point_set(
//Copy sample points
std::vector<Point> sample_points;
sample_points.reserve(number_of_sample);
unsigned int i;
unsigned int i;
for(it = first_sample_iter; it != points.end(); ++it)
{
sample_points.push_back(get(point_map, *it));
}
//compute default neighbor_radius, if no radius in
if (radius < 0)
{
@ -529,20 +464,19 @@ wlop_simplify_and_regularize_point_set(
#endif
}
FT radius2 = radius * radius;
CGAL_point_set_processing_precondition(radius > 0);
// Initiate a KD-tree search for original points
std::vector<Kd_tree_element> original_treeElements;
for (it = first_original_iter, i=0 ; it != points.end() ; ++it, ++i)
original_treeElements.push_back( Kd_tree_element(get(point_map, *it), i) );
Kd_Tree original_kd_tree(original_treeElements.begin(),
Kd_Tree original_kd_tree(original_treeElements.begin(),
original_treeElements.end());
std::vector<Point> update_sample_points(number_of_sample);
typename std::vector<Point>::iterator sample_iter;
// Compute original density weight for original points if user needed
std::vector<FT> original_density_weights;
@ -555,7 +489,7 @@ wlop_simplify_and_regularize_point_set(
compute_density_weight_for_original_point<Kernel, Kd_Tree>
(
get(point_map, *it),
original_kd_tree,
original_kd_tree,
radius);
original_density_weights.push_back(density);
@ -581,86 +515,57 @@ wlop_simplify_and_regularize_point_set(
{
FT density = simplify_and_regularize_internal::
compute_density_weight_for_sample_point<Kernel, Kd_Tree>
(*sample_iter,
sample_kd_tree,
(*sample_iter,
sample_kd_tree,
radius);
sample_density_weights.push_back(density);
}
typedef boost::zip_iterator<boost::tuple<typename std::vector<Point>::iterator,
typename std::vector<Point>::iterator> > Zip_iterator;
typename std::vector<Point>::iterator update_iter = update_sample_points.begin();
#ifndef CGAL_LINKED_WITH_TBB
CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
"Parallel_tag is enabled but TBB is unavailable.");
#else
//parallel
if (boost::is_convertible<ConcurrencyTag, Parallel_tag>::value)
{
Point_set_processing_3::internal::Parallel_callback
parallel_callback (callback, iter_number * number_of_sample, iter_n * number_of_sample);
Point_set_processing_3::internal::Callback_wrapper<ConcurrencyTag>
callback_wrapper (callback, iter_number * number_of_sample, iter_n * number_of_sample);
tbb::blocked_range<size_t> block(0, number_of_sample);
Sample_point_updater<Kernel, Kd_Tree, typename PointRange::iterator> sample_updater(
update_sample_points,
sample_points,
original_kd_tree,
sample_kd_tree,
radius2,
original_density_weights,
sample_density_weights,
parallel_callback.advancement(),
parallel_callback.interrupted());
CGAL::for_each<ConcurrencyTag>
(CGAL::make_range (boost::make_zip_iterator (boost::make_tuple (sample_points.begin(), update_sample_points.begin())),
boost::make_zip_iterator (boost::make_tuple (sample_points.end(), update_sample_points.end()))),
[&](const typename Zip_iterator::reference& t)
{
if (callback_wrapper.interrupted())
return false;
tbb::parallel_for(block, sample_updater);
get<1>(t) = simplify_and_regularize_internal::
compute_update_sample_point<Kernel, Kd_Tree, typename PointRange::iterator>(
get<0>(t),
original_kd_tree,
sample_kd_tree,
radius,
original_density_weights,
sample_density_weights);
++ callback_wrapper.advancement();
bool interrupted = parallel_callback.interrupted();
return true;
});
// We interrupt by hand as counter only goes halfway and won't terminate by itself
parallel_callback.interrupted() = true;
parallel_callback.join();
// If interrupted during this step, nothing is computed, we return NaN
if (interrupted)
return output;
}else
#endif
{
//sequential
std::size_t nb = iter_n * number_of_sample;
for (sample_iter = sample_points.begin();
sample_iter != sample_points.end(); ++sample_iter, ++update_iter, ++ nb)
{
*update_iter = simplify_and_regularize_internal::
compute_update_sample_point<Kernel,
Kd_Tree,
typename PointRange::iterator>
(*sample_iter,
original_kd_tree,
sample_kd_tree,
radius2,
original_density_weights,
sample_density_weights);
if (callback && !callback ((nb+1) / double(iter_number * number_of_sample)))
return output;
}
}
bool interrupted = callback_wrapper.interrupted();
// We interrupt by hand as counter only goes halfway and won't terminate by itself
callback_wrapper.interrupted() = true;
callback_wrapper.join();
// If interrupted during this step, nothing is computed, we return NaN
if (interrupted)
return output;
sample_iter = sample_points.begin();
for (update_iter = update_sample_points.begin();
update_iter != update_sample_points.end();
++update_iter, ++sample_iter)
{
*sample_iter = *update_iter;
}
for (std::size_t i = 0; i < sample_points.size(); ++ i)
sample_points[i] = update_sample_points[i];
}
// final output
for(sample_iter = sample_points.begin();
sample_iter != sample_points.end(); ++sample_iter)
{
*output++ = *sample_iter;
}
std::copy (sample_points.begin(), sample_points.end(), output);
return output;
}

124
Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h
View File

@ -48,7 +48,7 @@
#include <boost/type_traits/is_convertible.hpp>
#include <boost/utility/enable_if.hpp>
/*!
/*!
\file Poisson_reconstruction_function.h
*/
@ -104,15 +104,15 @@ struct Poisson_visitor {
{}
};
struct Poisson_skip_vertices {
struct Poisson_skip_vertices {
double ratio;
Random& m_random;
Poisson_skip_vertices(const double ratio, Random& random)
: ratio(ratio), m_random(random) {}
Poisson_skip_vertices(const double ratio = 0)
: ratio(ratio) {}
template <typename Iterator>
bool operator()(Iterator) const {
return m_random.get_double() < ratio;
bool operator()(Iterator it) const {
// make the result deterministic for each iterator
return Random((std::size_t)(&*it)).get_double() < ratio;
}
};
@ -123,7 +123,7 @@ template <typename F1, typename F2>
struct Special_wrapper_of_two_functions_keep_pointers {
F1 *f1;
F2 *f2;
Special_wrapper_of_two_functions_keep_pointers(F1* f1, F2* f2)
Special_wrapper_of_two_functions_keep_pointers(F1* f1, F2* f2)
: f1(f1), f2(f2) {}
template <typename X>
@ -136,16 +136,16 @@ struct Special_wrapper_of_two_functions_keep_pointers {
return (std::max)((*f1)(x), CGAL::square((*f2)(x)));
}
}; // end struct Special_wrapper_of_two_functions_keep_pointers<F1, F2>
/// \endcond
/// \endcond
/*!
\ingroup PkgPoissonSurfaceReconstruction3Ref
\brief Implementation of the Poisson Surface Reconstruction method.
Given a set of 3D points with oriented normals sampled on the boundary
of a 3D solid, the Poisson Surface Reconstruction method \cgalCite{Kazhdan06}
of a 3D solid, the Poisson Surface Reconstruction method \cgalCite{Kazhdan06}
solves for an approximate indicator function of the inferred
solid, whose gradient best matches the input normals. The output
scalar function, represented in an adaptive octree, is then
@ -155,7 +155,7 @@ iso-contoured using an adaptive marching cubes.
algorithm which solves for a piecewise linear function on a 3D
Delaunay triangulation instead of an adaptive octree.
\tparam Gt Geometric traits class.
\tparam Gt Geometric traits class.
\cgalModels `ImplicitFunction`
@ -166,7 +166,7 @@ class Poisson_reconstruction_function
// Public types
public:
/// \name Types
/// \name Types
/// @{
typedef Gt Geom_traits; ///< Geometric traits class
@ -180,7 +180,7 @@ public:
typedef typename Geom_traits::FT FT; ///< number type.
typedef typename Geom_traits::Point_3 Point; ///< point type.
typedef typename Geom_traits::Vector_3 Vector; ///< vector type.
typedef typename Geom_traits::Sphere_3 Sphere;
typedef typename Geom_traits::Sphere_3 Sphere;
/// @}
@ -269,22 +269,22 @@ private:
// Public methods
public:
/// \name Creation
/// \name Creation
/// @{
/*!
Creates a Poisson implicit function from the range of points `[first, beyond)`.
/*!
Creates a Poisson implicit function from the range of points `[first, beyond)`.
\tparam InputIterator iterator over input points.
\tparam InputIterator iterator over input points.
\tparam PointPMap is a model of `ReadablePropertyMap` with
a `value_type = Point`. It can be omitted if `InputIterator`
`value_type` is convertible to `Point`.
`value_type` is convertible to `Point`.
\tparam NormalPMap is a model of `ReadablePropertyMap`
with a `value_type = Vector`.
*/
*/
template <typename InputIterator,
typename PointPMap,
typename NormalPMap
@ -337,7 +337,7 @@ public:
: m_tr(new Triangulation), m_Bary(new std::vector<boost::array<double,9> > )
, average_spacing(CGAL::compute_average_spacing<CGAL::Sequential_tag>(CGAL::make_range(first, beyond), 6))
{
forward_constructor(first, beyond,
forward_constructor(first, beyond,
make_identity_property_map(
typename std::iterator_traits<InputIterator>::value_type()),
normal_pmap, Poisson_visitor());
@ -355,12 +355,12 @@ public:
{
return m_tr->bounding_sphere();
}
/// \cond SKIP_IN_MANUAL
const Triangulation& tr() const {
return *m_tr;
}
// This variant requires all parameters.
template <class SparseLinearAlgebraTraits_d,
class Visitor>
@ -368,7 +368,7 @@ public:
SparseLinearAlgebraTraits_d solver,// = SparseLinearAlgebraTraits_d(),
Visitor visitor,
double approximation_ratio = 0,
double average_spacing_ratio = 5)
double average_spacing_ratio = 5)
{
CGAL::Timer task_timer; task_timer.start();
CGAL_TRACE_STREAM << "Delaunay refinement...\n";
@ -382,7 +382,7 @@ public:
internal::Poisson::Constant_sizing_field<Triangulation> sizing_field(CGAL::square(cell_radius_bound));
std::vector<int> NB;
std::vector<int> NB;
NB.push_back( delaunay_refinement(radius_edge_ratio_bound,sizing_field,max_vertices,enlarge_ratio));
@ -391,7 +391,7 @@ public:
NB.push_back( delaunay_refinement(radius_edge_ratio_bound,sizing_field,max_vertices,enlarge_ratio));
}
if(approximation_ratio > 0. &&
if(approximation_ratio > 0. &&
approximation_ratio * std::distance(m_tr->input_points_begin(),
m_tr->input_points_end()) > 20) {
@ -417,16 +417,14 @@ public:
typedef Filter_iterator<typename Triangulation::Input_point_iterator,
Poisson_skip_vertices> Some_points_iterator;
//make it deterministic
Random random(0);
Poisson_skip_vertices skip(1.-approximation_ratio,random);
Poisson_skip_vertices skip(1.-approximation_ratio);
CGAL_TRACE_STREAM << "SPECIAL PASS that uses an approximation of the result (approximation ratio: "
<< approximation_ratio << ")" << std::endl;
CGAL::Timer approximation_timer; approximation_timer.start();
CGAL::Timer sizing_field_timer; sizing_field_timer.start();
Poisson_reconstruction_function<Geom_traits>
Poisson_reconstruction_function<Geom_traits>
coarse_poisson_function(Some_points_iterator(m_tr->input_points_end(),
skip,
m_tr->input_points_begin()),
@ -435,18 +433,18 @@ public:
Normal_of_point_with_normal_map<Geom_traits>() );
coarse_poisson_function.compute_implicit_function(solver, Poisson_visitor(),
0.);
internal::Poisson::Constant_sizing_field<Triangulation>
internal::Poisson::Constant_sizing_field<Triangulation>
min_sizing_field(CGAL::square(average_spacing));
internal::Poisson::Constant_sizing_field<Triangulation>
internal::Poisson::Constant_sizing_field<Triangulation>
sizing_field_ok(CGAL::square(average_spacing*average_spacing_ratio));
Special_wrapper_of_two_functions_keep_pointers<
internal::Poisson::Constant_sizing_field<Triangulation>,
Poisson_reconstruction_function<Geom_traits> > sizing_field2(&min_sizing_field,
&coarse_poisson_function);
sizing_field_timer.stop();
std::cerr << "Construction time of the sizing field: " << sizing_field_timer.time()
std::cerr << "Construction time of the sizing field: " << sizing_field_timer.time()
<< " seconds" << std::endl;
NB.push_back( delaunay_refinement(radius_edge_ratio_bound,
@ -458,12 +456,12 @@ public:
CGAL_TRACE_STREAM << "SPECIAL PASS END (" << approximation_timer.time() << " seconds)" << std::endl;
}
// Prints status
CGAL_TRACE_STREAM << "Delaunay refinement: " << "added ";
for(std::size_t i = 0; i < NB.size()-1; i++){
CGAL_TRACE_STREAM << NB[i] << " + ";
}
CGAL_TRACE_STREAM << NB[i] << " + ";
}
CGAL_TRACE_STREAM << NB.back() << " Steiner points, "
<< task_timer.time() << " seconds, "
<< std::endl;
@ -510,12 +508,12 @@ public:
If \ref thirdpartyEigen "Eigen" 3.1 (or greater) is available and `CGAL_EIGEN3_ENABLED`
is defined, an overload with \link Eigen_solver_traits <tt>Eigen_solver_traits<Eigen::ConjugateGradient<Eigen_sparse_symmetric_matrix<double>::EigenType> ></tt> \endlink
as default solver is provided.
\param solver sparse linear solver.
\param smoother_hole_filling controls if the Delaunay refinement is done for the input points, or for an approximation of the surface obtained from a first pass of the algorithm on a sample of the points.
\return `false` if the linear solver fails.
*/
\return `false` if the linear solver fails.
*/
template <class SparseLinearAlgebraTraits_d>
bool compute_implicit_function(SparseLinearAlgebraTraits_d solver, bool smoother_hole_filling = false)
{
@ -555,14 +553,14 @@ public:
}
/// \endcond
/*!
`ImplicitFunction` interface: evaluates the implicit function at a
given 3D query point. The function `compute_implicit_function()` must be
called before the first call to `operator()`.
*/
/*!
`ImplicitFunction` interface: evaluates the implicit function at a
given 3D query point. The function `compute_implicit_function()` must be
called before the first call to `operator()`.
*/
FT operator()(const Point& p) const
{
m_hint = m_tr->locate(p ,m_hint);
m_hint = m_tr->locate(p ,m_hint);
if(m_tr->is_infinite(m_hint)) {
int i = m_hint->index(m_tr->infinite_vertex());
@ -576,7 +574,7 @@ public:
c * m_hint->vertex(2)->f() +
d * m_hint->vertex(3)->f();
}
/// \cond SKIP_IN_MANUAL
void initialize_cell_indices()
{
@ -616,7 +614,7 @@ public:
void initialize_duals() const
{
Dual.resize(m_tr->number_of_cells());
Dual.resize(m_tr->number_of_cells());
int i = 0;
for(Finite_cells_iterator fcit = m_tr->finite_cells_begin();
fcit != m_tr->finite_cells_end();
@ -642,18 +640,18 @@ public:
const Point& pb = ch->vertex(1)->point();
const Point& pc = ch->vertex(2)->point();
const Point& pd = ch->vertex(3)->point();
Vector va = pa - pd;
Vector vb = pb - pd;
Vector vc = pc - pd;
internal::invert(va.x(), va.y(), va.z(),
vb.x(), vb.y(), vb.z(),
vc.x(), vc.y(), vc.z(),
entry[0],entry[1],entry[2],entry[3],entry[4],entry[5],entry[6],entry[7],entry[8]);
}
/// \endcond
/// Returns a point located inside the inferred surface.
Point get_inner_point() const
{
@ -684,7 +682,7 @@ private:
internal::Poisson::Constant_sizing_field<Triangulation>());
}
template <typename Sizing_field,
template <typename Sizing_field,
typename Second_sizing_field>
unsigned int delaunay_refinement(FT radius_edge_ratio_bound, ///< radius edge ratio bound (ignored if zero)
Sizing_field sizing_field, ///< cell radius bound (ignored if zero)
@ -794,12 +792,12 @@ private:
// it is closed using the smallest part of the sphere).
std::vector<Vertex_handle> convex_hull;
m_tr->adjacent_vertices (m_tr->infinite_vertex (),
std::back_inserter (convex_hull));
std::back_inserter (convex_hull));
unsigned int nb_negative = 0;
for (std::size_t i = 0; i < convex_hull.size (); ++ i)
if (convex_hull[i]->f() < 0.0)
++ nb_negative;
if(nb_negative > convex_hull.size () / 2)
flip_f();
@ -816,7 +814,7 @@ private:
Finite_vertices_iterator v, e;
for(v = m_tr->finite_vertices_begin(),
e= m_tr->finite_vertices_end();
v != e;
v != e;
v++)
if(v->type() == Triangulation::INPUT)
values.push_back(v->f());
@ -935,14 +933,14 @@ private:
// TODO: Some entities are computed too often
// - nn and area should not be computed for the face and its opposite face
//
//
// divergent
FT div_normalized(Vertex_handle v)
{
std::vector<Cell_handle> cells;
cells.reserve(32);
m_tr->incident_cells(v,std::back_inserter(cells));
FT div = 0;
typename std::vector<Cell_handle>::iterator it;
for(it = cells.begin(); it != cells.end(); it++)
@ -988,7 +986,7 @@ private:
std::vector<Cell_handle> cells;
cells.reserve(32);
m_tr->incident_cells(v,std::back_inserter(cells));
FT div = 0.0;
typename std::vector<Cell_handle>::iterator it;
for(it = cells.begin(); it != cells.end(); it++)
@ -996,14 +994,14 @@ private:
Cell_handle cell = *it;
if(m_tr->is_infinite(cell))
continue;
const int index = cell->index(v);
const Point& a = cell->vertex(m_tr->vertex_triple_index(index, 0))->point();
const Point& b = cell->vertex(m_tr->vertex_triple_index(index, 1))->point();
const Point& c = cell->vertex(m_tr->vertex_triple_index(index, 2))->point();
const Vector nn = CGAL::cross_product(b-a,c-a);
div+= nn * (//v->normal() +
div+= nn * (//v->normal() +
cell->vertex((index+1)%4)->normal() +
cell->vertex((index+2)%4)->normal() +
cell->vertex((index+3)%4)->normal());
@ -1200,7 +1198,7 @@ private:
A.set_coef(vi->index(),vi->index(), diagonal, true /*new*/);
}
}
/// Computes enlarged geometric bounding sphere of the embedded triangulation.
Sphere enlarged_bounding_sphere(FT ratio) const

48
Property_map/include/CGAL/property_map.h
View File

@ -184,7 +184,7 @@ struct Dereference_property_map
/// Free function to create a `Dereference_property_map` property map.
///
/// \relates Dereference_property_map
/// \relates Dereference_property_map
template <class Iter> // Type convertible to `key_type`
Dereference_property_map<typename CGAL::value_type_traits<Iter>::type>
make_dereference_property_map(Iter)
@ -216,9 +216,25 @@ struct Identity_property_map
/// @}
};
/// \cond SKIP_IN_MANUAL
template <typename T>
struct Identity_property_map_no_lvalue
{
typedef T key_type; ///< typedef to `T`
typedef T value_type; ///< typedef to `T`
typedef T reference; ///< typedef to `T`
typedef boost::readable_property_map_tag category; ///< `boost::readable_property_map_tag`
typedef Identity_property_map_no_lvalue<T> Self;
friend reference get(const Self&, const key_type& k) {return k;}
};
/// \endcond
/// Free function to create a `Identity_property_map` property map.
///
/// \relates Identity_property_map
/// \relates Identity_property_map
template <class T> // Key and value type
Identity_property_map<T>
make_identity_property_map(T)
@ -228,8 +244,8 @@ Identity_property_map<T>
/// \ingroup PkgPropertyMapRef
/// Property map that accesses the first item of a `std::pair`.
/// \tparam Pair Instance of `std::pair`.
/// Property map that accesses the first item of a `std::pair`.
/// \tparam Pair Instance of `std::pair`.
/// \cgalModels `LvaluePropertyMap`
///
/// \sa `CGAL::Second_of_pair_property_map<Pair>`
@ -253,9 +269,9 @@ struct First_of_pair_property_map
/// @}
};
/// Free function to create a `First_of_pair_property_map` property map.
/// Free function to create a `First_of_pair_property_map` property map.
///
/// \relates First_of_pair_property_map
/// \relates First_of_pair_property_map
template <class Pair> // Pair type
First_of_pair_property_map<Pair>
make_first_of_pair_property_map(Pair)
@ -264,13 +280,13 @@ First_of_pair_property_map<Pair>
}
/// \ingroup PkgPropertyMapRef
///
/// Property map that accesses the second item of a `std::pair`.
///
/// \tparam Pair Instance of `std::pair`.
///
///
/// Property map that accesses the second item of a `std::pair`.
///
/// \tparam Pair Instance of `std::pair`.
///
/// \cgalModels `LvaluePropertyMap`
///
///
/// \sa `CGAL::First_of_pair_property_map<Pair>`
template <typename Pair>
struct Second_of_pair_property_map
@ -294,7 +310,7 @@ struct Second_of_pair_property_map
/// Free function to create a Second_of_pair_property_map property map.
///
/// \relates Second_of_pair_property_map
/// \relates Second_of_pair_property_map
template <class Pair> // Pair type
Second_of_pair_property_map<Pair>
make_second_of_pair_property_map(Pair)
@ -303,12 +319,12 @@ Second_of_pair_property_map<Pair>
}
/// \ingroup PkgPropertyMapRef
///
///
/// Property map that accesses the Nth item of a `boost::tuple` or a `std::tuple`.
///
///
/// \tparam N %Index of the item to access.
/// \tparam Tuple Instance of `boost::tuple` or `std::tuple`.
///
///
/// \cgalModels `LvaluePropertyMap`
template <int N, typename Tuple>
struct Nth_of_tuple_property_map

136
STL_Extension/include/CGAL/for_each.h Normal file
View File

@ -0,0 +1,136 @@
// Copyright (c) 2020 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Simon Giraudot
#ifndef CGAL_FOR_EACH_H
#define CGAL_FOR_EACH_H
#include <CGAL/iterator.h>
#ifdef CGAL_LINKED_WITH_TBB
#include <tbb/parallel_for.h>
#include <tbb/blocked_range.h>
#include <tbb/scalable_allocator.h>
#endif // CGAL_LINKED_WITH_TBB
/*
CGAL::for_each<ConcurrencyTag>(Range, Function) does the following:
- if Sequential_tag is used, apply Function to all elements of Range
- if Parallel_tag is used:
* static assert that TBB is available
* if Range has random access iterators, use a TBB parallel_for
loop to apply Function to all elements of Range
* if Range doesn't have random access iterators, first copy the
iterators in a vector and then use a TBB parallel_for loop to
apply Function to all elements of Range
The loop is interrupted if `functor` returns false (it carries on
until the end otherwise).
*/
namespace CGAL {
namespace internal {
template <typename RangeRef, typename IteratorCategory>
void for_each (RangeRef range,
const std::function<bool(typename std::iterator_traits
<typename Range_iterator_type<RangeRef>::type>::reference)>& functor,
const Sequential_tag&,
IteratorCategory)
{
for (typename std::iterator_traits
<typename Range_iterator_type<RangeRef>::type>::reference r : range)
if (!functor(r))
break;
}
#ifdef CGAL_LINKED_WITH_TBB
template <typename RangeRef, typename IteratorCategory>
void for_each (RangeRef range,
const std::function<bool(typename std::iterator_traits
<typename Range_iterator_type<RangeRef>::type>::reference)>& functor,
const Parallel_tag&,
IteratorCategory)
{
std::size_t range_size = std::distance (range.begin(), range.end());
std::vector<typename Range_iterator_type<RangeRef>::type> iterators;
iterators.reserve (range_size);
for (typename Range_iterator_type<RangeRef>::type it = range.begin(); it != range.end(); ++ it)
iterators.push_back (it);
tbb::parallel_for (tbb::blocked_range<std::size_t>(0, range_size),
[&](const tbb::blocked_range<std::size_t>& r)
{
for (std::size_t i = r.begin(); i != r.end(); ++ i)
if (!functor (*(iterators[i])))
break;
});
}
template <typename RangeRef>
void for_each (RangeRef range,
const std::function<bool(typename std::iterator_traits
<typename Range_iterator_type<RangeRef>::type>::reference)>& functor,
const Parallel_tag&,
std::random_access_iterator_tag)
{
std::size_t range_size = std::distance (range.begin(), range.end());
tbb::parallel_for (tbb::blocked_range<std::size_t>(0, range_size),
[&](const tbb::blocked_range<std::size_t>& r)
{
for (std::size_t i = r.begin(); i != r.end(); ++ i)
if (!functor (*(range.begin() + i)))
break;
});
}
#endif
} // namespace internal
template <typename ConcurrencyTag, typename Range>
void for_each (const Range& range,
const std::function<bool(typename std::iterator_traits
<typename Range::const_iterator>::reference)>& functor)
{
#ifndef CGAL_LINKED_WITH_TBB
CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
"Parallel_tag is enabled but TBB is unavailable.");
#endif
internal::for_each<const Range&>
(range, functor,
ConcurrencyTag(),
typename std::iterator_traits<typename Range::const_iterator>::iterator_category());
}
template <typename ConcurrencyTag, typename Range>
void for_each (Range& range,
const std::function<bool(typename std::iterator_traits
<typename Range::iterator>::reference)>& functor)
{
#ifndef CGAL_LINKED_WITH_TBB
CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
"Parallel_tag is enabled but TBB is unavailable.");
#endif
internal::for_each<Range&>
(range, functor,
ConcurrencyTag(),
typename std::iterator_traits<typename Range::iterator>::iterator_category());
}
} // namespace CGAL
#endif // CGAL_FOR_EACH_H

139
STL_Extension/include/CGAL/iterator.h
View File

@ -1,16 +1,16 @@
// Copyright (c) 2003
// Copyright (c) 2003
// Utrecht University (The Netherlands),
// ETH Zurich (Switzerland),
// INRIA Sophia-Antipolis (France),
// Max-Planck-Institute Saarbruecken (Germany),
// and Tel-Aviv University (Israel). All rights reserved.
// and Tel-Aviv University (Israel). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Michael Hoffmann <hoffmann@inf.ethz.ch>
// Lutz Kettner <kettner@mpi-sb.mpg.de>
@ -172,11 +172,11 @@ class Oneset_iterator
void, void, void, void >
{
T* t;
public:
// types
typedef Oneset_iterator<T> Self;
public:
Oneset_iterator(T& t) : t(&t) {}
@ -202,65 +202,65 @@ public:
template < typename T >
class Const_oneset_iterator {
public:
// types
typedef std::random_access_iterator_tag iterator_category;
typedef std::ptrdiff_t difference_type;
typedef T value_type;
typedef value_type* pointer;
typedef value_type& reference;
typedef Const_oneset_iterator<T> Self;
typedef difference_type Diff;
typedef value_type Val;
typedef pointer Ptr;
typedef reference Ref;
// construction
Const_oneset_iterator( const T& t = T(), Diff n = 0)
: value( t), index( n)
{ }
// access
Ref operator * ( ) { return value; }
const value_type& operator * ( ) const { return value; }
Ptr operator -> ( ) { return &value; }
const value_type* operator -> ( ) const { return &value; }
// equality operator
bool operator == ( const Self& x) const { return ( index==x.index); }
bool operator != ( const Self& x) const { return ( index!=x.index); }
// forward operations
// ------------------
Self& operator ++ ( ) { ++index; return *this; }
Self operator ++ ( int) { Self tmp = *this; ++index; return tmp; }
// bidirectional operations
// ------------------------
Self& operator -- ( ) { --index; return *this; }
Self operator -- ( int) { Self tmp = *this; --index; return tmp; }
// random access operations
// ------------------------
// access
Ref operator [] ( Diff ) { return value;}
const value_type& operator [] ( Diff ) const { return value;}
// less operator
bool operator < ( const Self& x) const { return ( index < x.index);}
// arithmetic operations
Self& operator += ( Diff n) { index += n; return *this; }
Self& operator -= ( Diff n) { index -= n; return *this; }
Self operator + ( Diff n) const { Self tmp = *this; return tmp+=n; }
Self operator - ( Diff n) const { Self tmp = *this; return tmp-=n; }
Diff operator - ( const Self& x) const { return index - x.index; }
private:
// data members
Val value;
Diff index;
@ -443,12 +443,12 @@ public:
};
// Microsoft 1300 cannot handle the default template parameters. Hence, ...
template < class I, int N, class Ref, class Ptr,
class Val, class Dist, class Ctg >
template < class I, int N, class Ref, class Ptr,
class Val, class Dist, class Ctg >
inline
N_step_adaptor<I,N,Ref,Ptr,Val,Dist,Ctg>
operator+(typename N_step_adaptor<I,N,Ref,Ptr,Val,Dist,Ctg>::difference_type n,
N_step_adaptor<I,N,Ref,Ptr,Val,Dist,Ctg> i)
N_step_adaptor<I,N,Ref,Ptr,Val,Dist,Ctg> i)
{ return i += n; }
template < class I, int N>
@ -605,7 +605,7 @@ public:
--(*this);
return tmp;
}
reference operator*() const { return *c_; }
pointer operator->() const { return &*c_; }
const Predicate& predicate() const { return p_; }
@ -678,9 +678,9 @@ public:
: i1(it.i1), op(it.op) {}
Join_input_iterator_1(I1 i,const Op& o=Op())
: i1(i), op(o) {}
I1 current_iterator1() const { return i1; }
bool operator==(const Self& i) const {
return i1 == i.i1;
}
@ -695,12 +695,12 @@ public:
op = it.op;
return *this;
}
const value_type& operator*() const {
const value_type& operator*() const {
val = op(*i1);
return val;
}
Self& operator++( ) {
++i1;
return *this;
@ -711,12 +711,12 @@ public:
return *this;
}
Self operator--(int) { Self tmp = *this; --(*this); return tmp; }
const value_type& operator[](difference_type i) const {
val = op(i1[i]);
return val;
}
Self& operator+=(difference_type n) {
i1 += n;
return *this;
@ -758,17 +758,17 @@ protected:
mutable value_type val; // Note: mutable is needed because we want to
// return a reference in operator*() and
// operator[](int) below.
public:
Join_input_iterator_2() {}
Join_input_iterator_2(const Join_input_iterator_2& it)
: i1(it.i1), i2(it.i2), op(it.op) {}
Join_input_iterator_2(I1 i1,I2 i2,const Op& op=Op())
: i1(i1), i2(i2), op(op) {}
I1 current_iterator1() const { return i1; }
I2 current_iterator2() const { return i2; }
bool operator==(const Self& i) const {
return i1 == i.i1 && i2 == i.i2;
}
@ -776,7 +776,7 @@ public:
bool operator< (const Self& i) const {
return i1 < i.i1 && i2 < i.i2;
}
Join_input_iterator_2& operator=(const Join_input_iterator_2& it)
{
i1 = it.i1;
@ -785,11 +785,11 @@ public:
return *this;
}
const value_type& operator*() const {
const value_type& operator*() const {
val = op(*i1,*i2);
return val;
}
Self& operator++( ) {
++i1;
++i2;
@ -802,12 +802,12 @@ public:
return *this;
}
Self operator--(int) { Self tmp = *this; --(*this); return tmp; }
const value_type& operator[](difference_type i) const {
val = op(i1[i],i2[i]);
return val;
}
Self& operator+=(difference_type n) {
i1 += n;
i2 += n;
@ -845,7 +845,7 @@ public:
typedef typename std::iterator_traits<I1>::difference_type difference_type;
typedef value_type* pointer;
typedef value_type& reference;
protected:
I1 i1;
I2 i2;
@ -854,18 +854,18 @@ protected:
mutable value_type val; // Note: mutable is needed because we want to
// return a reference in operator*() and
// operator[](int) below.
public:
Join_input_iterator_3() {}
Join_input_iterator_3(const Join_input_iterator_3& it)
: i1(it.i1), i2(it.i2), i3(it.i3), op(it.op) {}
Join_input_iterator_3(I1 i1,I2 i2,I3 i3,const Op& op=Op())
: i1(i1), i2(i2), i3(i3), op(op) {}
I1 current_iterator1() const { return i1; }
I2 current_iterator2() const { return i2; }
I2 current_iterator3() const { return i3; }
bool operator==(const Self& i) const {
return i1 == i.i1 && i2 == i.i2 && i3 == i.i3;
}
@ -873,7 +873,7 @@ public:
bool operator< (const Self& i) const {
return i1 < i.i1 && i2 < i.i2 && i3 < i.i3;
}
Join_input_iterator_3& operator=(const Join_input_iterator_3& it)
{
i1 = it.i1;
@ -883,11 +883,11 @@ public:
return *this;
}
const value_type& operator*() const {
const value_type& operator*() const {
val = op(*i1,*i2,*i3);
return val;
}
Self& operator++( ) {
++i1;
++i2;
@ -902,12 +902,12 @@ public:
return *this;
}
Self operator--(int) { Self tmp = *this; --(*this); return tmp; }
const value_type& operator[](difference_type i) const {
val = op(i1[i],i2[i],i3[i]);
return val;
}
Self& operator+=(difference_type n) {
i1 += n;
i2 += n;
@ -1220,17 +1220,17 @@ template<typename _Iterator, typename Predicate>
public:
typedef _Iterator iterator_type;
explicit Filter_output_iterator(_Iterator& __x, const Predicate& pred)
: iterator(__x), predicate(pred)
explicit Filter_output_iterator(_Iterator& __x, const Predicate& pred)
: iterator(__x), predicate(pred)
{}
template <typename T>
Filter_output_iterator&
operator=(const T& t)
{
if(! predicate(t))
*iterator = t;
return *this;
if(! predicate(t))
*iterator = t;
return *this;
}
Filter_output_iterator&
@ -1239,9 +1239,9 @@ template<typename _Iterator, typename Predicate>
Filter_output_iterator&
operator++()
{
{
++iterator;
return *this;
return *this;
}
Filter_output_iterator
@ -1249,7 +1249,7 @@ template<typename _Iterator, typename Predicate>
{
Filter_output_iterator res(*this);
++iterator;
return res;
return res;
}
};
@ -1264,7 +1264,7 @@ template<typename OutputIterator>
struct Output_visitor : boost::static_visitor<OutputIterator&> {
Output_visitor(OutputIterator* it) : out(it) {}
OutputIterator* out;
template<typename T>
OutputIterator& operator()(const T& t) {
*(*out)++ = t;
@ -1298,7 +1298,7 @@ struct Derivator<D, std::tuple<V1, V...>, std::tuple<O1, O...> >
Self& operator=(const Self&) = delete;
using Base::operator=;
D& operator=(const V1& v)
{
* std::get< D::size - sizeof...(V) - 1 >(static_cast<typename D::Iterator_tuple&>(static_cast<D&>(*this))) ++ = v;
@ -1308,7 +1308,7 @@ struct Derivator<D, std::tuple<V1, V...>, std::tuple<O1, O...> >
template <class Tuple>
void tuple_dispatch(const Tuple& t)
{
* std::get< D::size - sizeof...(V) - 1 >(static_cast<typename D::Iterator_tuple&>(static_cast<D&>(*this))) ++ =
* std::get< D::size - sizeof...(V) - 1 >(static_cast<typename D::Iterator_tuple&>(static_cast<D&>(*this))) ++ =
std::get< D::size - sizeof...(V) - 1 >(t);
static_cast<Base&>(*this).tuple_dispatch(t);
}
@ -1366,7 +1366,7 @@ public:
Dispatch_output_iterator(O... o) : std::tuple<O...>(o...) {}
Dispatch_output_iterator(const Dispatch_output_iterator&)=default;
Self& operator=(const Self& s)
@ -1374,7 +1374,7 @@ public:
static_cast<Iterator_tuple&>(*this) = static_cast<const Iterator_tuple&>(s);
return *this;
}
template<BOOST_VARIANT_ENUM_PARAMS(typename T)>
Self& operator=(const boost::variant<BOOST_VARIANT_ENUM_PARAMS(T) >& t) {
internal::Output_visitor<Self> visitor(this);
@ -1402,18 +1402,18 @@ public:
Self& operator*() { return *this; }
const Iterator_tuple& get_iterator_tuple() const { return *this; }
Self& operator=(const std::tuple<V...>& t)
{
tuple_dispatch(t);
return *this;
}
operator std::tuple<O&...>()
{
return tuple_internal::to_tuple(*this, std::index_sequence_for<O...>{});
}
operator std::tuple<const O&...>()const
{
return tuple_internal::to_tuple(*this, std::index_sequence_for<O...>{});
@ -1447,7 +1447,7 @@ class Dispatch_or_drop_output_iterator < std::tuple<V...>, std::tuple<O...> >
public:
Dispatch_or_drop_output_iterator(O... o) : Base(o...) {}
Dispatch_or_drop_output_iterator(const Dispatch_or_drop_output_iterator&)=default;
Dispatch_or_drop_output_iterator& operator=(const Dispatch_or_drop_output_iterator&)=default;
@ -1471,6 +1471,17 @@ dispatch_or_drop_output(O... o)
return Dispatch_or_drop_output_iterator<std::tuple<V...>, std::tuple<O...> >(o...);
}
// Trick to select iterator or const_iterator depending on the range constness
template <typename RangeRef>
struct Range_iterator_type;
template <typename RangeRef>
struct Range_iterator_type<RangeRef&> { typedef typename RangeRef::iterator type; };
template <typename RangeRef>
struct Range_iterator_type<const RangeRef&> { typedef typename RangeRef::const_iterator type; };
} //namespace CGAL
#include <CGAL/enable_warnings.h>

8
STL_Extension/test/STL_Extension/CMakeLists.txt
View File

@ -42,9 +42,13 @@ if ( CGAL_FOUND )
create_single_source_cgal_program( "test_Uncertain.cpp" )
create_single_source_cgal_program( "test_vector.cpp" )
create_single_source_cgal_program( "test_join_iterators.cpp" )
create_single_source_cgal_program( "test_for_each.cpp" )
if(TBB_FOUND)
CGAL_target_use_TBB(test_for_each)
endif()
else()
message(STATUS "This program requires the CGAL library, and will not be compiled.")
endif()

48
STL_Extension/test/STL_Extension/test_for_each.cpp Normal file
View File

@ -0,0 +1,48 @@
#include <CGAL/for_each.h>
#include <vector>
#include <list>
int main()
{
std::vector<int> vec { 0, 1, 2, 3, 4, 5 };
std::list<int> list { 0, 1, 2, 3, 4, 5 };
std::cerr << "Testing sequential random access" << std::endl;
CGAL::for_each<CGAL::Sequential_tag>
(vec, [](int& i) -> bool { i *= 2; return true; });
for (const int& i : vec)
std::cerr << i << " ";
std::cerr << std::endl;
#ifdef CGAL_LINKED_WITH_TBB
std::cerr << "Testing parallel random access" << std::endl;
CGAL::for_each<CGAL::Parallel_tag>
(vec, [](int& i) -> bool { i *= 2; return true; });
for (const int& i : vec)
std::cerr << i << " ";
std::cerr << std::endl;
#endif
std::cerr << "Testing sequential non-random access" << std::endl;
CGAL::for_each<CGAL::Sequential_tag>
(list, [](int& i) -> bool { i *= 2; return true; });
for (const int& i : list)
std::cerr << i << " ";
std::cerr << std::endl;
#ifdef CGAL_LINKED_WITH_TBB
std::cerr << "Testing parallel non-random access" << std::endl;
CGAL::for_each<CGAL::Parallel_tag>
(list, [](int& i) -> bool { i *= 2; return true; });
for (const int& i : list)
std::cerr << i << " ";
std::cerr << std::endl;
#endif
return 0;
}

258
Spatial_searching/include/CGAL/Kd_tree_rectangle.h
View File

@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Hans Tangelder (<hanst@cs.uu.nl>)
@ -31,42 +31,42 @@ namespace CGAL {
T *lower;
T *upper;
Construct_cartesian_const_iterator_d construct_it;
set_bounds_from_pointer(int d, T *l, T *u,Construct_cartesian_const_iterator_d construct_it_)
set_bounds_from_pointer(int d, T *l, T *u,Construct_cartesian_const_iterator_d construct_it_)
: dim(d), lower(l), upper(u), construct_it(construct_it_)
{}
void
operator()(P p)
void
operator()(P p)
{
T h;
typename Construct_cartesian_const_iterator_d::result_type pit = construct_it(*p);
for (int i = 0; i < dim; ++i, ++pit) {
h=(*pit);
if (h < lower[i]) lower[i] = h;
if (h > upper[i]) upper[i] = h;
h=(*pit);
if (h < lower[i]) lower[i] = h;
if (h > upper[i]) upper[i] = h;
}
}
};
template <class FT_, typename D = Dynamic_dimension_tag>
template <class FT_, typename D = Dynamic_dimension_tag>
class Kd_tree_rectangle {
public:
typedef FT_ FT;
typedef FT T;
private:
//int dim;
std::array<T,D::value> lower_;
std::array<T,D::value> upper_;
int max_span_coord_;
public:
inline void
set_upper_bound(int i, const FT& x)
inline void
set_upper_bound(int i, const FT& x)
{
CGAL_assertion(i >= 0 && i < D::value);
CGAL_assertion(x >= lower_[i]);
@ -74,8 +74,8 @@ namespace CGAL {
set_max_span();
}
inline void
set_lower_bound(int i, const FT& x)
inline void
set_lower_bound(int i, const FT& x)
{
CGAL_assertion(i >= 0 && i < D::value);
CGAL_assertion(x <= upper_[i]);
@ -83,20 +83,20 @@ namespace CGAL {
set_max_span();
}
inline void
set_max_span()
inline void
set_max_span()
{
FT span = upper_[0]-lower_[0];
max_span_coord_ = 0;
for (int i = 1; i < D::value; ++i) {
FT tmp = upper_[i] - lower_[i];
if (span < tmp) {
span = tmp;
max_span_coord_ = i;
}
FT tmp = upper_[i] - lower_[i];
if (span < tmp) {
span = tmp;
max_span_coord_ = i;
}
}
}
Kd_tree_rectangle(int)
: max_span_coord_(0)
{
@ -104,29 +104,30 @@ namespace CGAL {
upper_.fill(FT(0));
}
Kd_tree_rectangle()
Kd_tree_rectangle()
: max_span_coord_(-1)
{}
explicit
explicit
Kd_tree_rectangle(const Kd_tree_rectangle& r)
: max_span_coord_(r.max_span_coord_)
: max_span_coord_(r.max_span_coord_)
{
lower_ = r.lower_;
upper_ = r.upper_;
}
template <class Construct_cartesian_const_iterator_d,class PointPointerIter>
void update_from_point_pointers(PointPointerIter begin,
void update_from_point_pointers(PointPointerIter begin,
PointPointerIter end,
const Construct_cartesian_const_iterator_d& construct_it
)
)
{
if (begin ==end)
return;
// initialize with values of first point
typename Construct_cartesian_const_iterator_d::result_type bit = construct_it(**begin);
for (int i=0; i < D::value; ++i, ++bit) {
lower_[i]= *bit; upper_[i]=lower_[i];
}
@ -135,77 +136,78 @@ namespace CGAL {
std::for_each(begin, end,set_bounds_from_pointer<Construct_cartesian_const_iterator_d,P,T>(D::value, &(lower_[0]), &(upper_[0]), construct_it));
set_max_span();
}
template <class Construct_cartesian_const_iterator_d,class PointPointerIter> // was PointIter
Kd_tree_rectangle(int, PointPointerIter begin, PointPointerIter end,const Construct_cartesian_const_iterator_d& construct_it)
: max_span_coord_(-1)
{
update_from_point_pointers<Construct_cartesian_const_iterator_d>(begin,end,construct_it);
}
inline int
max_span_coord() const
{
return max_span_coord_;
inline int
max_span_coord() const
{
return max_span_coord_;
}
inline FT
max_span() const
inline FT
max_span() const
{
return upper_[max_span_coord_] - lower_[max_span_coord_];
}
inline FT
min_coord(int i) const
inline FT
min_coord(int i) const
{
CGAL_assume(i<D::value);
CGAL_assertion(lower_.size() != 0);
return lower_[i];
}
inline FT
max_coord(int i) const
inline FT
max_coord(int i) const
{
CGAL_assume(i<D::value);
return upper_[i];
}
std::ostream&
print(std::ostream& s) const
std::ostream&
print(std::ostream& s) const
{
s << "Rectangle dimension = " << D::value;
s << "\n lower: ";
for (int i=0; i < D::value; ++i)
s << lower_[i] << " ";
s << lower_[i] << " ";
// std::copy(lower_, lower_ + D,
// std::ostream_iterator<FT>(s," "));
// std::ostream_iterator<FT>(s," "));
s << "\n upper: ";
for (int j=0; j < D::value; ++j)
s << upper_[j] << " ";
s << upper_[j] << " ";
// std::copy(upper_, upper_ + D,
// std::ostream_iterator<FT>(s," "));
// std::ostream_iterator<FT>(s," "));
s << "\n maximum span " << max_span() <<
" at coordinate " << max_span_coord() << std::endl;
" at coordinate " << max_span_coord() << std::endl;
return s;
}
// Splits rectangle by modifying itself to lower half
// Splits rectangle by modifying itself to lower half
// and returns upper half
// Kd_tree_rectangle*
// Kd_tree_rectangle*
void
split(Kd_tree_rectangle& r, int d, FT value)
{
CGAL_assertion(d >= 0 && d < D::value);
CGAL_assertion(lower_[d] <= value && value <= upper_[d]);
//Kd_tree_rectangle* r = new Kd_tree_rectangle(*this);
upper_[d]=value;
r.lower_[d]=value;
//return r;
}
int
dimension() const
int
dimension() const
{
return D::value;
}
@ -214,7 +216,7 @@ namespace CGAL {
T* upper() {return upper_.data();}
const T* lower() const {return lower_.data();}
const T* upper() const {return upper_.data();}
Kd_tree_rectangle&
operator=(const Kd_tree_rectangle& r)
{
@ -222,35 +224,35 @@ namespace CGAL {
if (this != &r) {
lower_ = r.lower_;
upper_ = r.upper_;
set_max_span();
set_max_span();
}
return *this;
}
}; // of class Kd_tree_rectangle
// Partial specialization for dynamic dimension, which means dimension at runtime
template <class FT_>
template <class FT_>
class Kd_tree_rectangle<FT_, Dynamic_dimension_tag> {
public:
typedef FT_ FT;
typedef FT T;
private:
T* coords_;
int dim;
int max_span_coord_;
public:
inline void
set_upper_bound(int i, const FT& x)
inline void
set_upper_bound(int i, const FT& x)
{
CGAL_assertion(i >= 0 && i < dim);
CGAL_assertion(x >= lower()[i]);
@ -258,8 +260,8 @@ namespace CGAL {
set_max_span();
}
inline void
set_lower_bound(int i, const FT& x)
inline void
set_lower_bound(int i, const FT& x)
{
CGAL_assertion(i >= 0 && i < dim);
CGAL_assertion(x <= upper()[i]);
@ -267,51 +269,51 @@ namespace CGAL {
set_max_span();
}
inline void
set_max_span()
inline void
set_max_span()
{
FT span = upper()[0]-lower()[0];
max_span_coord_ = 0;
for (int i = 1; i < dim; ++i) {
FT tmp = upper()[i] - lower()[i];
if (span < tmp) {
span = tmp;
max_span_coord_ = i;
}
FT tmp = upper()[i] - lower()[i];
if (span < tmp) {
span = tmp;
max_span_coord_ = i;
}
}
}
Kd_tree_rectangle(int d)
Kd_tree_rectangle(int d)
: coords_(new FT[2*d]), dim(d), max_span_coord_(0)
{
std::fill(coords_, coords_ + 2*dim, FT(0));
}
Kd_tree_rectangle()
: coords_(0), dim(0)
Kd_tree_rectangle()
: coords_(0), dim(0), max_span_coord_(-1)
{
}
explicit
explicit
Kd_tree_rectangle(const Kd_tree_rectangle& r)
: coords_(new FT[2*r.dim]), dim(r.dim),
max_span_coord_(r.max_span_coord_)
max_span_coord_(r.max_span_coord_)
{
std::copy(r.coords_, r.coords_+2*dim, lower());
}
template <class Construct_cartesian_const_iterator_d,class PointPointerIter>
void update_from_point_pointers(PointPointerIter begin,
void update_from_point_pointers(PointPointerIter begin,
PointPointerIter end,
const Construct_cartesian_const_iterator_d& construct_it
)
)
{
if (begin ==end)
return;
// initialize with values of first point
typename Construct_cartesian_const_iterator_d::result_type bit = construct_it(**begin);
for (int i=0; i < dim; ++i, ++bit) {
lower()[i]= *bit; upper()[i]=lower()[i];
}
@ -320,83 +322,83 @@ namespace CGAL {
std::for_each(begin, end,set_bounds_from_pointer<Construct_cartesian_const_iterator_d,P,T>(dim, lower(), upper(),construct_it));
set_max_span();
}
template <class Construct_cartesian_const_iterator_d,class PointPointerIter> // was PointIter
Kd_tree_rectangle(int d, PointPointerIter begin, PointPointerIter end,const Construct_cartesian_const_iterator_d& construct_it)
: coords_(new FT[2*d]), dim(d)
: coords_(new FT[2*d]), dim(d), max_span_coord_(-1)
{
update_from_point_pointers<Construct_cartesian_const_iterator_d>(begin,end,construct_it);
}
inline int
max_span_coord() const
{
return max_span_coord_;
inline int
max_span_coord() const
{
return max_span_coord_;
}
inline FT
max_span() const
inline FT
max_span() const
{
return upper()[max_span_coord_] - lower()[max_span_coord_];
}
inline FT
min_coord(int i) const
inline FT
min_coord(int i) const
{
CGAL_assertion(coords_ != nullptr);
return lower()[i];
}
inline FT
max_coord(int i) const
inline FT
max_coord(int i) const
{
return upper()[i];
}
std::ostream&
print(std::ostream& s) const
std::ostream&
print(std::ostream& s) const
{
s << "Rectangle dimension = " << dim;
s << "\n lower: ";
for (int i=0; i < dim; ++i)
s << lower()[i] << " ";
s << lower()[i] << " ";
// std::copy(lower(), lower() + dim,
// std::ostream_iterator<FT>(s," "));
// std::ostream_iterator<FT>(s," "));
s << "\n upper: ";
for (int j=0; j < dim; ++j)
s << upper()[j] << " ";
s << upper()[j] << " ";
// std::copy(upper(), upper() + dim,
// std::ostream_iterator<FT>(s," "));
// std::ostream_iterator<FT>(s," "));
s << "\n maximum span " << max_span() <<
" at coordinate " << max_span_coord() << std::endl;
" at coordinate " << max_span_coord() << std::endl;
return s;
}
// Splits rectangle by modifying itself to lower half
// Splits rectangle by modifying itself to lower half
// and returns upper half
// Kd_tree_rectangle*
// Kd_tree_rectangle*
void
split(Kd_tree_rectangle& r, int d, FT value)
{
CGAL_assertion(d >= 0 && d < dim);
CGAL_assertion(lower()[d] <= value && value <= upper()[d]);
//Kd_tree_rectangle* r = new Kd_tree_rectangle(*this);
upper()[d]=value;
r.lower()[d]=value;
//return r;
}
~Kd_tree_rectangle()
~Kd_tree_rectangle()
{
if (dim) {
if (coords_) delete [] coords_;
if (coords_) delete [] coords_;
}
}
int
dimension() const
int
dimension() const
{
return dim;
}
@ -405,30 +407,30 @@ namespace CGAL {
T* upper() {return coords_ + dim;}
const T* lower() const {return coords_;}
const T* upper() const {return coords_ + dim;}
Kd_tree_rectangle&
operator=(const Kd_tree_rectangle& r)
{
CGAL_assertion(dimension() == r.dimension());
if (this != &r) {
std::copy(r.coords_, r.coords_+2*dim, coords_);
set_max_span();
set_max_span();
}
return *this;
}
}; // of partial specialization of class Kd_tree_rectangle<FT,0>
template <class FT, typename D>
std::ostream&
operator<<(std::ostream& s, const Kd_tree_rectangle<FT,D>& r)
std::ostream&
operator<<(std::ostream& s, const Kd_tree_rectangle<FT,D>& r)
{
return r.print(s);
}
} // namespace CGAL
#endif // CGAL_KD_TREE_RECTANGLE_H

72
Spatial_searching/include/CGAL/Search_traits_adapter.h
View File

@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Sebastien Loriot
@ -30,13 +30,13 @@
namespace CGAL{
using ::get;
template <class Point_with_info,class PointPropertyMap,class Base_traits>
class Search_traits_adapter;
template <class Point_with_info,class PointPropertyMap,class Base_distance>
class Distance_adapter;
namespace internal{
BOOST_MPL_HAS_XXX_TRAIT_NAMED_DEF(Has_typedef_Iso_box_d,Iso_box_d,false)
@ -55,7 +55,7 @@ struct Get_iso_box_d<T,true>
{
typedef typename T::Iso_box_d type;
};
template <class Point_with_info,class PointPropertyMap,class Base_traits>
struct Spatial_searching_default_distance< ::CGAL::Search_traits_adapter<Point_with_info,PointPropertyMap,Base_traits> >{
typedef ::CGAL::Distance_adapter<Point_with_info,
@ -64,8 +64,8 @@ struct Get_iso_box_d<T,true>
};
} //namespace internal
template <class Point_with_info,class PointPropertyMap,class Base_traits>
class Search_traits_adapter : public Base_traits{
PointPropertyMap ppmap;
@ -78,21 +78,21 @@ public:
const Base_traits& base=Base_traits()
):Base_traits(base),ppmap(ppmap_){}
typedef Point_with_info Point_d;
typedef typename Base_traits::FT FT;
typedef typename Base_traits::Dimension Dimension;
// Default if point map is lvalue: use Construct_cartesian_const_iterator_d
struct Construct_cartesian_const_iterator_d_lvalue: public Base_traits::Construct_cartesian_const_iterator_d{
PointPropertyMap ppmap;
typedef typename Base_traits::Construct_cartesian_const_iterator_d Base;
Construct_cartesian_const_iterator_d_lvalue
(const typename Base_traits::Construct_cartesian_const_iterator_d& base, const PointPropertyMap& ppmap_)
:Base_traits::Construct_cartesian_const_iterator_d(base), ppmap(ppmap_){}
typename Base_traits::Cartesian_const_iterator_d operator()(const Point_with_info& p) const
{ return Base::operator() (get(ppmap,p)); }
@ -100,9 +100,9 @@ public:
{ return Base::operator() (get(ppmap,p),0); }
// These 2 additional operators forward the call to Base_traits.
// This is needed because of an undocumented requirement of
// Orthogonal_k_neighbor_search and Orthogonal_incremental_neighbor_search:
// Traits::Construct_cartesian_const_iterator should be callable
// This is needed because of an undocumented requirement of
// Orthogonal_k_neighbor_search and Orthogonal_incremental_neighbor_search:
// Traits::Construct_cartesian_const_iterator should be callable
// on the query point type. If the query point type is the same as
// Point_with_info, we disable it.
@ -122,7 +122,7 @@ public:
) const
{ return Base::operator() (p,0); }
};
// If point map is not lvalue, use this work-around that stores a
// Point object in a shared pointer to avoid iterating on a temp
// object
@ -136,7 +136,7 @@ public:
typename std::iterator_traits<typename Base::Cartesian_const_iterator_d>::value_type,
std::random_access_iterator_tag
> Facade;
typedef typename std::iterator_traits<typename Base::Cartesian_const_iterator_d>::value_type
Dereference_type;
typedef typename boost::property_traits<PointPropertyMap>::value_type
@ -146,13 +146,13 @@ public:
std::size_t idx;
public:
No_lvalue_iterator() : point(NULL), idx(0) { }
No_lvalue_iterator(const Point& point) : point(new Point(point)), idx(0) { }
No_lvalue_iterator(const Point& point, int) : point(new Point(point)), idx(Base::Dimension::value) { }
private:
friend class boost::iterator_core_access;
void increment()
{
@ -164,7 +164,7 @@ public:
--idx;
CGAL_assertion(point != boost::shared_ptr<Point>());
}
void advance(std::ptrdiff_t n)
{
idx += n;
@ -174,7 +174,7 @@ public:
std::ptrdiff_t distance_to(const No_lvalue_iterator& other) const
{
return other.idx - this->idx;
}
bool equal(const No_lvalue_iterator& other) const
{
@ -182,7 +182,11 @@ public:
}
Dereference_type&
dereference() const { return const_cast<Dereference_type&>((*point)[idx]); }
dereference() const
{
// Point::operator[] takes an int as parameter...
return const_cast<Dereference_type&>((*point)[static_cast<int>(idx)]);
}
};
@ -191,11 +195,11 @@ public:
struct Construct_cartesian_const_iterator_d_no_lvalue {
typedef No_lvalue_iterator result_type;
PointPropertyMap ppmap;
Construct_cartesian_const_iterator_d_no_lvalue
(const typename Base_traits::Construct_cartesian_const_iterator_d&, const PointPropertyMap& ppmap_)
: ppmap(ppmap_) { }
No_lvalue_iterator operator()(const Point_with_info& p) const
{ return No_lvalue_iterator(get(ppmap, p)); }
@ -203,9 +207,9 @@ public:
{ return No_lvalue_iterator(get(ppmap, p),0); }
// These 2 additional operators forward the call to Base_traits.
// This is needed because of an undocumented requirement of
// Orthogonal_k_neighbor_search and Orthogonal_incremental_neighbor_search:
// Traits::Construct_cartesian_const_iterator should be callable
// This is needed because of an undocumented requirement of
// Orthogonal_k_neighbor_search and Orthogonal_incremental_neighbor_search:
// Traits::Construct_cartesian_const_iterator should be callable
// on the query point type. If the query point type is the same as
// Point_with_info, we disable it.
@ -242,7 +246,7 @@ public:
Construct_cartesian_const_iterator_d_lvalue,
Construct_cartesian_const_iterator_d_no_lvalue>::type
Construct_cartesian_const_iterator_d;
struct Construct_iso_box_d: public Base::Construct_iso_box_d{
PointPropertyMap ppmap;
typedef typename Base_traits::FT FT; // needed for VC++, because otherwise it is taken from the private typedef of the base class
@ -256,9 +260,9 @@ public:
return Base_functor::operator() (get(ppmap,p),get(ppmap,q));
}
};
const PointPropertyMap& point_property_map() const {return ppmap;}
Construct_cartesian_const_iterator_d construct_cartesian_const_iterator_d_object() const {
return Construct_cartesian_const_iterator_d(
Base::construct_cartesian_const_iterator_d_object(),
@ -271,7 +275,7 @@ class Distance_adapter : public Base_distance {
PointPropertyMap ppmap;
public:
Distance_adapter( const PointPropertyMap& ppmap_=PointPropertyMap(),
const Base_distance& distance=Base_distance()
):Base_distance(distance),ppmap(ppmap_){}
@ -282,7 +286,7 @@ public:
typedef Point_with_info Point_d;
typedef typename Base_distance::Query_item Query_item;
const PointPropertyMap& point_property_map() const {return ppmap;}
const PointPropertyMap& point_property_map() const {return ppmap;}
FT transformed_distance(const Query_item& p1, const Point_with_info& p2) const
{
@ -305,12 +309,12 @@ public:
FT max_distance_to_rectangle(const Query_item& p,const CGAL::Kd_tree_rectangle<FT,Dimension>& b) const
{
return Base_distance::max_distance_to_rectangle(p,b);
}
}
template <class FT,class Dimension>
FT max_distance_to_rectangle(const Query_item& p,const CGAL::Kd_tree_rectangle<FT,Dimension>& b,std::vector<FT>& dists) const
{
return Base_distance::max_distance_to_rectangle(p,b,dists);
}
}
};
}//namespace CGAL