Merge pull request #3869 from sgiraudot/PSP-Range_neighborhood_NP-GF

[Small Feature] PSP: Spherical Neighborhood Queries
This commit is contained in:
Laurent Rineau 2019-05-06 15:50:10 +02:00
commit 554672f363
11 changed files with 312 additions and 222 deletions

View File

@ -93,7 +93,7 @@ along sharp features. The range of possible values is `[0, 1]`. \n
\cgalNPBegin{neighbor_radius} \anchor PSP_neighbor_radius
is the spherical neighborhood radius.\n
\b Type: floating scalar value\n
<b>Default value</b>: `-1`, in which case the functions that use it will compute an adapted value automatically.
<b>Default value</b>: `-1`, in which case it is not used.
\cgalNPEnd
\cgalNPBegin{number_of_output_points} \anchor PSP_number_of_output_points

View File

@ -281,11 +281,10 @@ points has the appearance of a curve in 2D or a surface in 3D
- `estimate_global_k_neighbor_scale()`
- `estimate_global_range_scale()`
Functions such as `jet_estimate_normals()` or
`remove_outliers()` require a K neighbor scale while others such as
`grid_simplify_point_set()` require a range
scale. `vcm_estimate_normals()` is an example of a function that
accepts both.
Functions such as `grid_simplify_point_set()` require a range scale
while `jet_estimate_normals()`, `remove_outliers()` or
`vcm_estimate_normals()` are examples of functions that accepts both a
K neighbor scale or a range scale.
In some specific cases, the scale of a point set might not be
homogeneous (for example if the point set contains variable
@ -325,7 +324,9 @@ points in the domain.
Function `remove_outliers()` deletes a user-specified fraction
of outliers from an input point set. More specifically, it sorts the
input points in increasing order of average squared distances to their
`k` nearest neighbors and deletes the points with largest value.
nearest neighbors and deletes the points with largest value. The user
can either specify a fixed number of nearest neighbors or a fixed
spherical neighborhood radius.
\subsection Point_set_processing_3Example_outlier_removal Example
@ -446,13 +447,16 @@ Two smoothing functions are devised to smooth an input point set.
Function `jet_smooth_point_set()` smooths the input point set by
projecting each point onto a smooth parametric surface patch
(so-called jet surface) fitted over its `k` nearest neighbors.
(so-called jet surface) fitted over its nearest neighbors.
Function `bilateral_smooth_point_set()` smooths the input point set by
iteratively projecting each point onto the implicit surface patch fitted over its `k` nearest neighbors.
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.
Normals are thus required as input. For more details, see section 4 of \cgalCite{ear-2013}.
For both functions, the user can either specify a fixed number of
nearest neighbors or a fixed spherical neighborhood radius.
\subsection Point_set_processing_3Example_jet_smoothing Jet Smoothing Example
The following example generates a set of 9 points close to the `xy`
@ -487,13 +491,13 @@ functions provide an estimate of the normal to \b S at each
point. The result is an unoriented normal vector for each input point.
Function `jet_estimate_normals()` estimates the normal direction
at each point from the input set by fitting a jet surface over its `k`
at each point from the input set by fitting a jet surface over its
nearest neighbors. The default jet is a quadric surface. This
algorithm is well suited to point sets scattered over curved surfaces.
Function `pca_estimate_normals()` estimates the normal direction
at each point from the set by linear least squares fitting of a plane
over its `k` nearest neighbors. This algorithm is simpler and
over its nearest neighbors. This algorithm is simpler and
faster than `jet_estimate_normals()`.
Function `vcm_estimate_normals()` estimates the normal direction
@ -502,13 +506,16 @@ of the point set. This algorithm is more complex and slower than
the previous algorithms. It is based on the
article \cgalCite{cgal:mog-vbcfe-11}.
For these three functions, the user can either specify a fixed number
of nearest neighbors or a fixed spherical neighborhood radius.
\section Point_set_processing_3NormalOrientation Normal Orientation
Function `mst_orient_normals()` orients the normals of a set of
points with unoriented normals using the method described by Hoppe et
al. in <I>Surface reconstruction from unorganized points</I> \cgalCite{cgal:hddms-srup-92}.
More specifically, this method constructs a
Riemannian graph over the input points (the graph of the `k`
Riemannian graph over the input points (the graph of the
nearest neighbor points) and propagates a seed normal orientation
within a minimum spanning tree computed over this graph. The result is
an oriented normal vector for each input unoriented normal, except for
@ -521,8 +528,9 @@ Normal orientation of a sampled cube surface. Left: unoriented normals. Right: o
\subsection Point_set_processing_3Example_normals Example
The following example reads a point set from a file, estimates the
normals through PCA over the 6 nearest neighbors and orients the
normals:
normals through PCA (either over the 18 nearest neighbors or using a
spherical neighborhood radius of twice the average spacing) and
orients the normals:
\cgalExample{Point_set_processing_3/normals_example.cpp}

View File

@ -1,4 +1,5 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/compute_average_spacing.h>
#include <CGAL/pca_estimate_normals.h>
#include <CGAL/mst_orient_normals.h>
#include <CGAL/property_map.h>
@ -42,11 +43,31 @@ int main(int argc, char*argv[])
// Note: pca_estimate_normals() requiresa range of points
// as well as property maps to access each point's position and normal.
const int nb_neighbors = 18; // K-nearest neighbors = 3 rings
CGAL::pca_estimate_normals<Concurrency_tag>
(points, nb_neighbors,
CGAL::parameters::point_map(CGAL::First_of_pair_property_map<PointVectorPair>()).
normal_map(CGAL::Second_of_pair_property_map<PointVectorPair>()));
if (argc > 2 && std::strcmp(argv[2], "-r") == 0) // Use a fixed neighborhood radius
{
// First compute a spacing using the K parameter
double spacing
= CGAL::compute_average_spacing<Concurrency_tag>
(points, nb_neighbors,
CGAL::parameters::point_map(CGAL::First_of_pair_property_map<PointVectorPair>()));
// Then, estimate normals with a fixed radius
CGAL::pca_estimate_normals<Concurrency_tag>
(points,
0, // when using a neighborhood radius, K=0 means no limit on the number of neighbors returns
CGAL::parameters::point_map(CGAL::First_of_pair_property_map<PointVectorPair>()).
normal_map(CGAL::Second_of_pair_property_map<PointVectorPair>()).
neighbor_radius(2. * spacing)); // use 2*spacing as neighborhood radius
}
else // Use a fixed number of neighbors
{
CGAL::pca_estimate_normals<Concurrency_tag>
(points, nb_neighbors,
CGAL::parameters::point_map(CGAL::First_of_pair_property_map<PointVectorPair>()).
normal_map(CGAL::Second_of_pair_property_map<PointVectorPair>()));
}
// Orients normals.
// Note: mst_orient_normals() requires a range of points
// as well as property maps to access each point's position and normal.

View File

@ -0,0 +1,113 @@
// Copyright (c) 2019 GeometryFactory Sarl (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0+
//
// 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

View File

@ -27,6 +27,7 @@
#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/property_map.h>
#include <CGAL/point_set_processing_assertions.h>
#include <CGAL/Point_with_normal_3.h>
@ -185,58 +186,6 @@ compute_denoise_projection(
return Pwn(update_point, update_normal);
}
/// Computes neighbors from kdtree.
///
/// \pre `k >= 2`.
///
/// @tparam Kernel Geometric traits class.
/// @tparam Tree KD-tree.
///
/// @return neighbors pwn of query point.
template < typename Kernel,
typename Tree>
std::vector<CGAL::Point_with_normal_3<Kernel>,
CGAL_PSP3_DEFAULT_ALLOCATOR<CGAL::Point_with_normal_3<Kernel> > >
compute_kdtree_neighbors(
const CGAL::Point_with_normal_3<Kernel>& query, ///< 3D point
const Tree& tree, ///< KD-tree
unsigned int k ///< number of neighbors
)
{
// basic geometric types
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;
unsigned int i;
std::vector<CGAL::Point_with_normal_3<Kernel>
, CGAL_PSP3_DEFAULT_ALLOCATOR<CGAL::Point_with_normal_3<Kernel> >
> neighbor_pwns;
for(i = 0; i < (k+1); ++i)
{
if(search_iterator == search.end())
break; // premature ending
Pwn pwn = search_iterator->first;
neighbor_pwns.push_back(pwn);
++search_iterator;
}
// output
return neighbor_pwns;
}
/// Computes max-spacing of one query point from K nearest neighbors.
///
/// \pre `k >= 2`.
@ -302,6 +251,7 @@ class Compute_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;
@ -309,11 +259,12 @@ class Compute_pwns_neighbors
cpp11::atomic<bool>& interrupted;
public:
Compute_pwns_neighbors(unsigned int k, const Tree &tree,
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_tree(tree), m_pwns(pwns), m_pwns_neighbors(neighbors)
: 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
@ -322,9 +273,10 @@ public:
{
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]);
m_pwns_neighbors[i] = bilateral_smooth_point_set_internal::
compute_kdtree_neighbors<Kernel, Tree>(m_pwns[i], m_tree, m_k);
++ advancement;
}
}
@ -392,7 +344,7 @@ public:
\ingroup PkgPointSetProcessing3Algorithms
This function smooths an input point set by iteratively projecting each
point onto the implicit surface patch fitted over its k nearest neighbors.
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}.
@ -422,6 +374,12 @@ 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{neighbor_radius} spherical neighborhood radius. If
provided, the neighborhood of a query point is computed with a fixed spherical
radius instead of a fixed number of neighbors. In that case, the parameter
`k` is used as a limit on the number of points returned by each spherical
query (to avoid overly large number of points in high density areas). If no
limit is wanted, use `k=0`.\cgalParamEnd
\cgalParamBegin{sharpness_angle} controls the sharpness of the result.\cgalParamEnd
\cgalParamBegin{callback} an instance of
`std::function<bool(double)>`. It is called regularly when the
@ -477,7 +435,8 @@ bilateral_smooth_point_set(
PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap());
NormalMap normal_map = choose_param(get_param(np, internal_np::normal_map), NormalMap());
FT neighbor_radius = choose_param(get_param(np, internal_np::neighbor_radius), FT(0));
// copy points and normals
Pwns pwns;
for(typename PointRange::iterator it = points.begin(); it != points.end(); ++it)
@ -546,7 +505,7 @@ bilateral_smooth_point_set(
Point_set_processing_3::internal::Parallel_callback
parallel_callback (callback, 2 * nb_points);
Compute_pwns_neighbors<Kernel, Tree> f(k, tree, pwns, pwns_neighbors,
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);
@ -570,8 +529,9 @@ bilateral_smooth_point_set(
std::size_t nb = 0;
for(pwn_iter = pwns.begin(); pwn_iter != pwns.end(); ++pwn_iter, ++pwns_iter, ++ nb)
{
*pwns_iter = bilateral_smooth_point_set_internal::
compute_kdtree_neighbors<Kernel, Tree>(*pwn_iter, tree, k);
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();
}

View File

@ -28,6 +28,7 @@
#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/Monge_via_jet_fitting.h>
#include <CGAL/property_map.h>
#include <CGAL/point_set_processing_assertions.h>
@ -74,39 +75,22 @@ typename Kernel::Vector_3
jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute the normal at
Tree& tree, ///< KD-tree
unsigned int k, ///< number of neighbors
typename Kernel::FT neighbor_radius,
unsigned int degree_fitting)
{
// basic geometric types
typedef typename Kernel::Point_3 Point;
// types for K nearest neighbors search
typedef typename CGAL::Search_traits_3<Kernel> Tree_traits;
typedef typename CGAL::Orthogonal_k_neighbor_search<Tree_traits> Neighbor_search;
typedef typename Neighbor_search::iterator Search_iterator;
// types for jet fitting
typedef Monge_via_jet_fitting< Kernel,
Simple_cartesian<double>,
SvdTraits> Monge_jet_fitting;
typedef typename Monge_jet_fitting::Monge_form Monge_form;
// 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.
std::vector<Point> 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);
std::vector<Point> points;
CGAL::Point_set_processing_3::internal::neighbor_query
(query, tree, k, neighbor_radius, points);
// performs jet fitting
Monge_jet_fitting monge_fit;
const unsigned int degree_monge = 1; // we seek for normal and not more.
@ -120,10 +104,12 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
#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;
@ -131,11 +117,13 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
cpp11::atomic<bool>& interrupted;
public:
Jet_estimate_normals(Tree& tree, unsigned int k, std::vector<Point>& points,
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), degree_fitting (degree_fitting), input (points), output (output)
: tree(tree), k (k), neighbor_radius (neighbor_radius)
, degree_fitting (degree_fitting), input (points), output (output)
, advancement (advancement)
, interrupted (interrupted)
{ }
@ -146,7 +134,7 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
{
if (interrupted)
break;
output[i] = CGAL::internal::jet_estimate_normal<Kernel,SvdTraits>(input[i], tree, k, degree_fitting);
output[i] = CGAL::internal::jet_estimate_normal<Kernel,SvdTraits>(input[i], tree, k, neighbor_radius, degree_fitting);
++ advancement;
}
}
@ -168,7 +156,7 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
/**
\ingroup PkgPointSetProcessing3Algorithms
Estimates normal directions of the range of `points`
using jet fitting on the k nearest neighbors.
using jet fitting on the nearest neighbors.
The output normals are randomly oriented.
\pre `k >= 2`
@ -188,6 +176,12 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
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{neighbor_radius} spherical neighborhood radius. If
provided, the neighborhood of a query point is computed with a fixed spherical
radius instead of a fixed number of neighbors. In that case, the parameter
`k` is used as a limit on the number of points returned by each spherical
query (to avoid overly large number of points in high density areas). If no
limit is wanted, use `k=0`.\cgalParamEnd
\cgalParamBegin{degree_fitting} degree of jet fitting.\cgalParamEnd
\cgalParamBegin{svd_traits} template parameter for the class `Monge_via_jet_fitting`. If
\ref thirdpartyEigen "Eigen" 3.2 (or greater) is available and `CGAL_EIGEN3_ENABLED` is defined,
@ -220,6 +214,7 @@ jet_estimate_normals(
typedef typename Point_set_processing_3::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::FT FT;
typedef typename GetSvdTraits<NamedParameters>::type SvdTraits;
CGAL_static_assertion_msg(!(boost::is_same<NormalMap,
@ -232,6 +227,8 @@ jet_estimate_normals(
PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap());
NormalMap normal_map = choose_param(get_param(np, internal_np::normal_map), NormalMap());
unsigned int degree_fitting = choose_param(get_param(np, internal_np::degree_fitting), 2);
FT neighbor_radius = choose_param(get_param(np, internal_np::neighbor_radius), FT(0));
const std::function<bool(double)>& callback = choose_param(get_param(np, internal_np::callback),
std::function<bool(double)>());
@ -251,7 +248,7 @@ jet_estimate_normals(
CGAL_point_set_processing_precondition(points.begin() != points.end());
// precondition: at least 2 nearest neighbors
CGAL_point_set_processing_precondition(k >= 2);
CGAL_point_set_processing_precondition(k >= 2 || neighbor_radius > FT(0));
std::size_t memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20);
CGAL_TRACE(" Creates KD-tree\n");
@ -282,7 +279,8 @@ jet_estimate_normals(
std::vector<Vector> normals (kd_tree_points.size (),
CGAL::NULL_VECTOR);
CGAL::internal::Jet_estimate_normals<Kernel, SvdTraits, Tree>
f (tree, k, kd_tree_points, degree_fitting, normals,
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);
@ -301,7 +299,7 @@ jet_estimate_normals(
{
Vector normal = internal::jet_estimate_normal<Kernel,SvdTraits,Tree>(
get(point_map,*it),
tree, k, degree_fitting);
tree, k, neighbor_radius, degree_fitting);
put(normal_map, *it, normal); // normal_map[it] = normal
if (callback && !callback ((nb+1) / double(kd_tree_points.size())))

View File

@ -28,6 +28,7 @@
#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/Monge_via_jet_fitting.h>
#include <CGAL/property_map.h>
#include <CGAL/point_set_processing_assertions.h>
@ -74,39 +75,22 @@ jet_smooth_point(
const typename Kernel::Point_3& query, ///< 3D point to project
Tree& tree, ///< KD-tree
const unsigned int k, ///< number of neighbors.
typename Kernel::FT neighbor_radius,
const unsigned int degree_fitting,
const unsigned int degree_monge)
{
// basic geometric types
typedef typename Kernel::Point_3 Point;
// types for K nearest neighbors search
typedef typename CGAL::Search_traits_3<Kernel> Tree_traits;
typedef typename CGAL::Orthogonal_k_neighbor_search<Tree_traits> Neighbor_search;
typedef typename Neighbor_search::iterator Search_iterator;
// types for jet fitting
typedef Monge_via_jet_fitting< Kernel,
Simple_cartesian<double>,
SvdTraits> Monge_jet_fitting;
typedef typename Monge_jet_fitting::Monge_form Monge_form;
// Gather set of (k+1) neighboring points.
// Performs k + 1 queries (if unique the query point is
// output first). Search may be aborted if k is greater
// than number of input points.
std::vector<Point> 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);
std::vector<Point> points;
CGAL::Point_set_processing_3::internal::neighbor_query
(query, tree, k, neighbor_radius, points);
// performs jet fitting
Monge_jet_fitting monge_fit;
@ -123,6 +107,7 @@ jet_smooth_point(
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;
@ -131,12 +116,14 @@ jet_smooth_point(
cpp11::atomic<bool>& interrupted;
public:
Jet_smooth_pwns (Tree& tree, unsigned int k, std::vector<Point>& points,
unsigned int degree_fitting, unsigned int degree_monge, std::vector<Point>& output,
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), degree_fitting (degree_fitting),
degree_monge (degree_monge), input (points), output (output)
: tree(tree), k (k), neighbor_radius(neighbor_radius)
, degree_fitting (degree_fitting)
, degree_monge (degree_monge), input (points), output (output)
, advancement (advancement)
, interrupted (interrupted)
{ }
@ -148,6 +135,7 @@ jet_smooth_point(
if (interrupted)
break;
output[i] = CGAL::internal::jet_smooth_point<Kernel, SvdTraits>(input[i], tree, k,
neighbor_radius,
degree_fitting,
degree_monge);
++ advancement;
@ -170,7 +158,7 @@ jet_smooth_point(
/**
\ingroup PkgPointSetProcessing3Algorithms
Smoothes the range of `points` using jet fitting on the k
Smoothes the range of `points` using jet fitting on the
nearest neighbors and reprojection onto the jet.
As this method relocates the points, it
should not be called on containers sorted w.r.t. point locations.
@ -190,6 +178,12 @@ jet_smooth_point(
\cgalNamedParamsBegin
\cgalParamBegin{point_map} a model of `ReadablePropertyMap` with value type `geom_traits::Point_3`.
If this parameter is omitted, `CGAL::Identity_property_map<geom_traits::Point_3>` is used.\cgalParamEnd
\cgalParamBegin{neighbor_radius} spherical neighborhood radius. If
provided, the neighborhood of a query point is computed with a fixed spherical
radius instead of a fixed number of neighbors. In that case, the parameter
`k` is used as a limit on the number of points returned by each spherical
query (to avoid overly large number of points in high density areas). If no
limit is wanted, use `k=0`.\cgalParamEnd
\cgalParamBegin{degree_fitting} degree of jet fitting.\cgalParamEnd
\cgalParamBegin{degree_monge} Monge degree.\cgalParamEnd
\cgalParamBegin{svd_traits} template parameter for the class `Monge_via_jet_fitting`. If
@ -228,6 +222,8 @@ jet_smooth_point_set(
"Error: no SVD traits");
PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap());
typename Kernel::FT neighbor_radius = choose_param(get_param(np, internal_np::neighbor_radius),
typename Kernel::FT(0));
unsigned int degree_fitting = choose_param(get_param(np, internal_np::degree_fitting), 2);
unsigned int degree_monge = choose_param(get_param(np, internal_np::degree_monge), 2);
const std::function<bool(double)>& callback = choose_param(get_param(np, internal_np::callback),
@ -271,8 +267,8 @@ jet_smooth_point_set(
std::vector<Point> mutated_points (kd_tree_points.size (), CGAL::ORIGIN);
CGAL::internal::Jet_smooth_pwns<Kernel, SvdTraits, Tree>
f (tree, k, kd_tree_points, degree_fitting, degree_monge,
mutated_points,
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);
@ -293,7 +289,7 @@ jet_smooth_point_set(
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,degree_fitting,degree_monge) );
p,tree,k,neighbor_radius,degree_fitting,degree_monge) );
if (callback && !callback ((nb+1) / double(kd_tree_points.size())))
break;
}

View File

@ -28,6 +28,7 @@
#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/property_map.h>
#include <CGAL/Index_property_map.h>
@ -319,6 +320,7 @@ create_riemannian_graph(
IndexMap index_map, ///< property map ForwardIterator -> index
ConstrainedMap constrained_map, ///< property map ForwardIterator -> bool
unsigned int k, ///< number of neighbors
typename Kernel::FT neighbor_radius,
const Kernel& /*kernel*/) ///< geometric traits.
{
// Input points types
@ -331,7 +333,6 @@ create_riemannian_graph(
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 Neighbor_search::iterator Search_iterator;
// Riemannian_graph types
typedef internal::Riemannian_graph<ForwardIterator> Riemannian_graph;
@ -395,21 +396,15 @@ create_riemannian_graph(
std::size_t it_index = get(index_map,it);
Vector_ref it_normal_vector = get(normal_map,*it);
// 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.
Point_ref point = get(point_map, *it);
Point_vertex_handle_3 point_wrapper(point.x(), point.y(), point.z(), it);
Neighbor_search search(*tree, point_wrapper, k+1);
Search_iterator search_iterator = search.begin();
for(std::size_t i=0;i<(k+1);i++)
{
if(search_iterator == search.end())
break; // premature ending
std::vector<Point_vertex_handle_3> neighbor_points;
CGAL::Point_set_processing_3::internal::neighbor_query
(point_wrapper, *tree, k, neighbor_radius, neighbor_points);
ForwardIterator neighbor = search_iterator->first;
for (std::size_t i = 0; i < neighbor_points.size(); ++ i)
{
ForwardIterator neighbor = neighbor_points[i];
std::size_t neighbor_index = get(index_map,neighbor);
if (neighbor_index > it_index) // undirected graph
{
@ -431,8 +426,6 @@ create_riemannian_graph(
weight = 0; // safety check
riemannian_graph_weight_map[e] = (float)weight;
}
search_iterator++;
}
// Check if point is source
@ -595,6 +588,12 @@ create_mst_graph(
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{neighbor_radius} spherical neighborhood radius. If
provided, the neighborhood of a query point is computed with a fixed spherical
radius instead of a fixed number of neighbors. In that case, the parameter
`k` is used as a limit on the number of points returned by each spherical
query (to avoid overly large number of points in high density areas). If no
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,
@ -628,6 +627,8 @@ mst_orient_normals(
PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap());
NormalMap normal_map = choose_param(get_param(np, internal_np::normal_map), NormalMap());
typename Kernel::FT neighbor_radius = choose_param(get_param(np, internal_np::neighbor_radius),
typename Kernel::FT(0));
ConstrainedMap constrained_map = choose_param(get_param(np, internal_np::point_is_constrained), ConstrainedMap());
Kernel kernel;
@ -677,12 +678,14 @@ mst_orient_normals(
point_map, normal_map,
kernel)),
k,
neighbor_radius,
kernel);
else
riemannian_graph = create_riemannian_graph(points.begin(), points.end(),
point_map, normal_map, index_map,
constrained_map,
k,
neighbor_radius,
kernel);
// Creates a Minimum Spanning Tree starting at source_point

View File

@ -29,6 +29,7 @@
#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/linear_least_squares_fitting_3.h>
#include <CGAL/property_map.h>
#include <CGAL/point_set_processing_assertions.h>
@ -57,7 +58,6 @@ namespace CGAL {
/// \cond SKIP_IN_MANUAL
namespace internal {
/// Estimates normal direction using linear least
/// squares fitting of a plane on the K nearest neighbors.
///
@ -73,34 +73,17 @@ template < typename Kernel,
typename Kernel::Vector_3
pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute the normal at
const Tree& tree, ///< KD-tree
unsigned int k) ///< number of neighbors
unsigned int k, ///< number of neighbors
typename Kernel::FT neighbor_radius)
{
// basic geometric types
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::Plane_3 Plane;
// types for K nearest neighbors search
typedef typename CGAL::Search_traits_3<Kernel> Tree_traits;
typedef typename CGAL::Orthogonal_k_neighbor_search<Tree_traits> Neighbor_search;
typedef typename Neighbor_search::iterator Search_iterator;
// 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.
std::vector<Point> 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);
std::vector<Point> points;
CGAL::Point_set_processing_3::internal::neighbor_query
(query, tree, k, neighbor_radius, points);
// performs plane fitting by point-based PCA
Plane plane;
linear_least_squares_fitting_3(points.begin(),points.end(),plane,Dimension_tag<0>());
@ -113,21 +96,25 @@ pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
#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, std::vector<Point>& points,
std::vector<Vector>& output,
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), input (points), output (output)
: tree(tree), k (k), neighbor_radius (neighbor_radius)
, input (points), output (output)
, advancement (advancement)
, interrupted (interrupted)
{ }
@ -138,7 +125,7 @@ pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
{
if (interrupted)
break;
output[i] = CGAL::internal::pca_estimate_normal<Kernel,Tree>(input[i], tree, k);
output[i] = CGAL::internal::pca_estimate_normal<Kernel,Tree>(input[i], tree, k, neighbor_radius);
++ advancement;
}
}
@ -159,7 +146,7 @@ pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
/**
\ingroup PkgPointSetProcessing3Algorithms
Estimates normal directions of the range of `points`
by linear least squares fitting of a plane over the k nearest neighbors.
by linear least squares fitting of a plane over the nearest neighbors.
The output normals are randomly oriented.
\pre `k >= 2`
@ -174,20 +161,29 @@ pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
\param k number of neighbors
\param np optional sequence of \ref psp_namedparameters "Named Parameters" among the ones listed below.
\cgalNamedParamsBegin
\cgalParamBegin{point_map} a model of `ReadablePropertyMap` with value type `geom_traits::Point_3`.
If this parameter is omitted, `CGAL::Identity_property_map<geom_traits::Point_3>` is used.\cgalParamEnd
\cgalParamBegin{normal_map} a model of `WritablePropertyMap` with value type
\cgalNamedParamsBegin \cgalParamBegin{point_map} a model of
`ReadablePropertyMap` with value type `geom_traits::Point_3`. If
this parameter is omitted,
`CGAL::Identity_property_map<geom_traits::Point_3>` is
used.\cgalParamEnd \cgalParamBegin{normal_map} a model of
`WritablePropertyMap` with value type
`geom_traits::Vector_3`.\cgalParamEnd
\cgalParamBegin{neighbor_radius} spherical neighborhood
radius. If provided, the neighborhood of a query point is
computed with a fixed spherical radius instead of a fixed number
of neighbors. In that case, the parameter `k` is used as a limit
on the number of points returned by each spherical query (to
avoid overly large number of points in high density areas). If no
limit is wanted, use `k=0`.\cgalParamEnd
\cgalParamBegin{callback} an instance of
`std::function<bool(double)>`. It is called regularly when the
algorithm is running: the current advancement (between 0. and
1.) is passed as parameter. If it returns `true`, then the
algorithm continues its execution normally; if it returns
`false`, the algorithm is stopped and the remaining normals are
left unchanged.\cgalParamEnd
\cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
\cgalNamedParamsEnd
`std::function<bool(double)>`. It is called regularly when the
algorithm is running: the current advancement (between 0. and 1.)
is passed as parameter. If it returns `true`, then the algorithm
continues its execution normally; if it returns `false`, the
algorithm is stopped and the remaining normals are left
unchanged.\cgalParamEnd \cgalParamBegin{geom_traits} an instance
of a geometric traits class, model of `Kernel`\cgalParamEnd
\cgalNamedParamsEnd
*/
template <typename ConcurrencyTag,
typename PointRange,
@ -206,6 +202,7 @@ pca_estimate_normals(
typedef typename Point_set_processing_3::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::FT FT;
CGAL_static_assertion_msg(!(boost::is_same<NormalMap,
typename Point_set_processing_3::GetNormalMap<PointRange, NamedParameters>::NoMap>::value),
@ -213,6 +210,7 @@ pca_estimate_normals(
PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap());
NormalMap normal_map = choose_param(get_param(np, internal_np::normal_map), NormalMap());
FT neighbor_radius = choose_param(get_param(np, internal_np::neighbor_radius), FT(0));
const std::function<bool(double)>& callback = choose_param(get_param(np, internal_np::callback),
std::function<bool(double)>());
@ -263,7 +261,7 @@ pca_estimate_normals(
std::vector<Vector> normals (kd_tree_points.size (),
CGAL::NULL_VECTOR);
CGAL::internal::PCA_estimate_normals<Kernel, Tree>
f (tree, k, kd_tree_points, normals,
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);
@ -283,7 +281,7 @@ pca_estimate_normals(
Vector normal = internal::pca_estimate_normal<Kernel,Tree>(
get(point_map,*it),
tree,
k);
k, neighbor_radius);
put(normal_map, *it, normal); // normal_map[it] = normal
if (callback && !callback ((nb+1) / double(kd_tree_points.size())))

View File

@ -27,6 +27,7 @@
#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/property_map.h>
#include <CGAL/point_set_processing_assertions.h>
#include <functional>
@ -63,33 +64,16 @@ typename Kernel::FT
compute_avg_knn_sq_distance_3(
const typename Kernel::Point_3& query, ///< 3D point to project
Tree& tree, ///< KD-tree
unsigned int k) ///< number of neighbors
unsigned int k, ///< number of neighbors
typename Kernel::FT neighbor_radius)
{
// geometric types
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point;
// types for K nearest neighbors search
typedef typename CGAL::Search_traits_3<Kernel> Tree_traits;
typedef typename CGAL::Orthogonal_k_neighbor_search<Tree_traits> Neighbor_search;
typedef typename Neighbor_search::iterator Search_iterator;
// Gather set of (k+1) neighboring points.
// Perform k+1 queries (if in point set, the query point is
// output first). Search may be aborted if k is greater
// than number of input points.
std::vector<Point> 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);
std::vector<Point> points;
CGAL::Point_set_processing_3::internal::neighbor_query
(query, tree, k, neighbor_radius, points);
// compute average squared distance
typename Kernel::Compute_squared_distance_3 sqd;
@ -111,7 +95,7 @@ compute_avg_knn_sq_distance_3(
/**
\ingroup PkgPointSetProcessing3Algorithms
Removes outliers:
- computes average squared distance to the K nearest neighbors,
- computes average squared distance to the nearest neighbors,
- and sorts the points in increasing order of average distance.
This method modifies the order of input points so as to pack all remaining points first,
@ -130,6 +114,13 @@ compute_avg_knn_sq_distance_3(
\cgalNamedParamsBegin
\cgalParamBegin{point_map} a model of `ReadablePropertyMap` with value type `geom_traits::Point_3`.
If this parameter is omitted, `CGAL::Identity_property_map<geom_traits::Point_3>` is used.\cgalParamEnd
\cgalParamBegin{neighbor_radius} spherical neighborhood
radius. If provided, the neighborhood of a query point is
computed with a fixed spherical radius instead of a fixed number
of neighbors. In that case, the parameter `k` is used as a limit
on the number of points returned by each spherical query (to
avoid overly large number of points in high density areas). If no
limit is wanted, use `k=0`.\cgalParamEnd
\cgalParamBegin{threshold_percent} maximum percentage of points to remove.\cgalParamEnd
\cgalParamBegin{threshold_distance} minimum distance for a point to be considered as outlier
(distance here is the square root of the average squared distance to K nearest neighbors).\cgalParamEnd
@ -169,6 +160,8 @@ remove_outliers(
typedef typename Point_set_processing_3::GetK<PointRange, NamedParameters>::Kernel Kernel;
PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap());
typename Kernel::FT neighbor_radius = choose_param(get_param(np, internal_np::neighbor_radius),
typename Kernel::FT(0));
double threshold_percent = choose_param(get_param(np, internal_np::threshold_percent), 10.);
double threshold_distance = choose_param(get_param(np, internal_np::threshold_distance), 0.);
const std::function<bool(double)>& callback = choose_param(get_param(np, internal_np::callback),
@ -213,7 +206,7 @@ remove_outliers(
{
FT sq_distance = internal::compute_avg_knn_sq_distance_3<Kernel>(
get(point_map,*it),
tree, k);
tree, k, neighbor_radius);
sorted_points.insert( std::make_pair(sq_distance, *it) );
if (callback && !callback ((nb+1) / double(kd_tree_points.size())))
return points.end();

View File

@ -442,8 +442,8 @@ public:
\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
of neighboring sample points gives satisfactory result. The default value is set to 8 times the average
spacing of the point set.\cgalParamEnd
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
value is 35. More iterations give a more regular result but increase the runtime.\cgalParamEnd
\cgalParamBegin{require_uniform_sampling} an optional preprocessing, which will give better result if the