Merge branch 'Interpolation-Use_result_of-GF-old' into Interpolation-Use_result_of-GF

This commit is contained in:
Mael Rouxel-Labbé 2018-06-26 14:27:19 +02:00
commit 68c8b797a6
8 changed files with 524 additions and 168 deletions

View File

@ -86,8 +86,7 @@ function \ref PkgInterpolationNaturalNeighborCoordinates2.
\tparam OutputFunctor must be a functor with argument type `std::pair<Dt::Vertex_handle, Traits::Vector_d>`. \tparam OutputFunctor must be a functor with argument type `std::pair<Dt::Vertex_handle, Traits::Vector_d>`.
Note that the result type of the functor is not specified and is chosen by users to fit their needs. Note that the result type of the functor is not specified and is chosen by users to fit their needs.
\tparam ValueFunctor must be a functor where: \tparam ValueFunctor must be a functor where:
- `ValueFunctor::argument_type` must be either `std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT>` - `ValueFunctor::argument_type` must be either `Dt::Vertex_handle` or `Dt::Point`.
or `std::pair<Dt::Point, Dt::Geom_traits::FT>`.
- `ValueFunctor::result_type` is a pair of the function value type and a Boolean. - `ValueFunctor::result_type` is a pair of the function value type and a Boolean.
The function value type must provide a multiplication and addition operation with the type The function value type must provide a multiplication and addition operation with the type
`Traits::FT` as well as a constructor with argument `0`. `Traits::FT` as well as a constructor with argument `0`.
@ -124,8 +123,7 @@ functions \ref PkgInterpolationRegularNeighborCoordinates2.
\tparam OutputFunctor must be a functor with argument type `std::pair<Rt::Vertex_handle, Traits::Vector_d>`. \tparam OutputFunctor must be a functor with argument type `std::pair<Rt::Vertex_handle, Traits::Vector_d>`.
Note that the result type of the functor is not specified and is chosen by users to fit their needs. Note that the result type of the functor is not specified and is chosen by users to fit their needs.
\tparam ValueFunctor must be a functor where: \tparam ValueFunctor must be a functor where:
- `ValueFunctor::argument_type` must be either `std::pair<Rt::Vertex_handle, Rt::Geom_traits::FT>` - `ValueFunctor::argument_type` must be either `Rt::Vertex_handle` or `Rt::Weighted_point`.
or `std::pair<Rt::Weighted_point, Rt::Geom_traits::FT>`.
- `ValueFunctor::result_type` is a pair of the function value type and a Boolean. - `ValueFunctor::result_type` is a pair of the function value type and a Boolean.
The function value type must provide a multiplication and addition operation with the type The function value type must provide a multiplication and addition operation with the type
`Traits::FT` as well as a constructor with argument `0`. `Traits::FT` as well as a constructor with argument `0`.

View File

@ -94,6 +94,7 @@ int main()
Delaunay_triangulation T; Delaunay_triangulation T;
// Note that a supported alternative to creating the functors below is to use lambdas
Value_function<Vertex_handle, Coord_type> value_function; Value_function<Vertex_handle, Coord_type> value_function;
Gradient_function<Vertex_handle, Vector> gradient_function; Gradient_function<Vertex_handle, Vector> gradient_function;

View File

@ -73,6 +73,7 @@ int main()
{ {
Regular_triangulation rt; Regular_triangulation rt;
// Note that a supported alternative to creating the functors below is to use lambdas
Value_function<Vertex_handle, Coord_type> value_function; Value_function<Vertex_handle, Coord_type> value_function;
Gradient_function<Vertex_handle, Vector> gradient_function; Gradient_function<Vertex_handle, Vector> gradient_function;

View File

@ -74,6 +74,7 @@ int main()
{ {
Delaunay_triangulation dt; Delaunay_triangulation dt;
// Note that a supported alternative to creating the functors below is to use lambdas
Value_function<Vertex_handle, Coord_type> value_function; Value_function<Vertex_handle, Coord_type> value_function;
Gradient_function<Vertex_handle, Vector> gradient_function; Gradient_function<Vertex_handle, Vector> gradient_function;

View File

@ -27,6 +27,8 @@
#include <CGAL/double.h> #include <CGAL/double.h>
#include <CGAL/use.h> #include <CGAL/use.h>
#include <boost/utility/result_of.hpp>
#include <iterator> #include <iterator>
#include <utility> #include <utility>
#include <vector> #include <vector>
@ -58,17 +60,21 @@ struct Data_access
//the interpolation functions: //the interpolation functions:
template < class ForwardIterator, class ValueFunctor > template < class ForwardIterator, class ValueFunctor >
typename ValueFunctor::result_type::first_type typename boost::result_of<
ValueFunctor(typename std::iterator_traits<ForwardIterator>::value_type::first_type)>
::type::first_type
linear_interpolation(ForwardIterator first, ForwardIterator beyond, linear_interpolation(ForwardIterator first, ForwardIterator beyond,
const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm, const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm,
ValueFunctor value_function) ValueFunctor value_function)
{ {
CGAL_precondition(norm > 0); CGAL_precondition(norm > 0);
typedef typename ValueFunctor::result_type::first_type Value_type; typedef typename std::iterator_traits<ForwardIterator>::value_type::first_type arg_type;
typedef typename boost::result_of<ValueFunctor(arg_type)>::type result_type;
typedef typename result_type::first_type Value_type;
Value_type result(0); Value_type result(0);
typename ValueFunctor::result_type val; result_type val;
for(; first!=beyond; ++first) for(; first!=beyond; ++first)
{ {
val = value_function(first->first); val = value_function(first->first);
@ -78,9 +84,12 @@ linear_interpolation(ForwardIterator first, ForwardIterator beyond,
return result; return result;
} }
template < class ForwardIterator, class ValueFunctor, class GradFunctor, class Traits, class Point > template < class ForwardIterator, class ValueFunctor, class GradFunctor, class Traits, class Point >
std::pair< typename ValueFunctor::result_type::first_type, bool> std::pair<
typename boost::result_of<
ValueFunctor(typename std::iterator_traits<ForwardIterator>::value_type::first_type)>
::type::first_type,
bool>
quadratic_interpolation(ForwardIterator first, ForwardIterator beyond, quadratic_interpolation(ForwardIterator first, ForwardIterator beyond,
const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm, const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm,
const Point& p, const Point& p,
@ -89,15 +98,20 @@ quadratic_interpolation(ForwardIterator first, ForwardIterator beyond,
const Traits& traits) const Traits& traits)
{ {
CGAL_precondition(norm > 0); CGAL_precondition(norm > 0);
typedef typename ValueFunctor::result_type::first_type Value_type;
typedef typename Traits::Point_d Bare_point; typedef typename std::iterator_traits<ForwardIterator>::value_type::first_type arg_type;
typedef typename boost::result_of<ValueFunctor(arg_type)>::type value_functor_result_type;
typedef typename boost::result_of<GradFunctor(arg_type)>::type gradient_functor_result_type;
typedef typename value_functor_result_type::first_type Value_type;
typedef typename Traits::Point_d Bare_point;
Interpolation::internal::Extract_bare_point<Traits> cp(traits); Interpolation::internal::Extract_bare_point<Traits> cp(traits);
const Bare_point& bp = cp(p); const Bare_point& bp = cp(p);
Value_type result(0); Value_type result(0);
typename ValueFunctor::result_type f; value_functor_result_type f;
typename GradFunctor::result_type grad; gradient_functor_result_type grad;
for(; first!=beyond; ++first) for(; first!=beyond; ++first)
{ {
@ -119,7 +133,11 @@ quadratic_interpolation(ForwardIterator first, ForwardIterator beyond,
template < class ForwardIterator, class ValueFunctor, class GradFunctor, class Traits, class Point > template < class ForwardIterator, class ValueFunctor, class GradFunctor, class Traits, class Point >
std::pair< typename ValueFunctor::result_type::first_type, bool > std::pair<
typename boost::result_of<
ValueFunctor(typename std::iterator_traits<ForwardIterator>::value_type::first_type)>
::type::first_type,
bool>
sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond, sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond,
const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm, const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm,
const Point& p, const Point& p,
@ -129,17 +147,21 @@ sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond,
{ {
CGAL_precondition(norm >0); CGAL_precondition(norm >0);
typedef typename ValueFunctor::result_type::first_type Value_type; typedef typename std::iterator_traits<ForwardIterator>::value_type::first_type arg_type;
typedef typename Traits::FT Coord_type; typedef typename boost::result_of<ValueFunctor(arg_type)>::type value_functor_result_type;
typedef typename Traits::Point_d Bare_point; typedef typename boost::result_of<GradFunctor(arg_type)>::type gradient_functor_result_type;
typedef typename value_functor_result_type::first_type Value_type;
typedef typename Traits::FT Coord_type;
typedef typename Traits::Point_d Bare_point;
Interpolation::internal::Extract_bare_point<Traits> cp(traits); Interpolation::internal::Extract_bare_point<Traits> cp(traits);
const Bare_point& bp = cp(p); const Bare_point& bp = cp(p);
Coord_type term1(0), term2(term1), term3(term1), term4(term1); Coord_type term1(0), term2(term1), term3(term1), term4(term1);
Value_type linear_int(0), gradient_int(0); Value_type linear_int(0), gradient_int(0);
typename ValueFunctor::result_type f; value_functor_result_type f;
typename GradFunctor::result_type grad; gradient_functor_result_type grad;
for(; first!=beyond; ++first) for(; first!=beyond; ++first)
{ {
@ -193,7 +215,11 @@ sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond,
// gradient_int += (coeff/inv_weight) * (vh->get_value()+ vh->get_gradient() * (p - vh->point())); // gradient_int += (coeff/inv_weight) * (vh->get_value()+ vh->get_gradient() * (p - vh->point()));
template < class ForwardIterator, class ValueFunctor, class GradFunctor, class Traits, class Point > template < class ForwardIterator, class ValueFunctor, class GradFunctor, class Traits, class Point >
std::pair< typename ValueFunctor::result_type::first_type, bool > std::pair<
typename boost::result_of<
ValueFunctor(typename std::iterator_traits<ForwardIterator>::value_type::first_type)>
::type::first_type,
bool>
sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator beyond, sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator beyond,
const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm, const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm,
const Point& p, const Point& p,
@ -203,17 +229,21 @@ sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator beyond,
{ {
CGAL_precondition(norm > 0); CGAL_precondition(norm > 0);
typedef typename ValueFunctor::result_type::first_type Value_type; typedef typename std::iterator_traits<ForwardIterator>::value_type::first_type arg_type;
typedef typename Traits::FT Coord_type; typedef typename boost::result_of<ValueFunctor(arg_type)>::type value_functor_result_type;
typedef typename Traits::Point_d Bare_point; typedef typename boost::result_of<GradFunctor(arg_type)>::type gradient_functor_result_type;
typedef typename value_functor_result_type::first_type Value_type;
typedef typename Traits::FT Coord_type;
typedef typename Traits::Point_d Bare_point;
Interpolation::internal::Extract_bare_point<Traits> cp(traits); Interpolation::internal::Extract_bare_point<Traits> cp(traits);
const Bare_point& bp = cp(p); const Bare_point& bp = cp(p);
Coord_type term1(0), term2(term1), term3(term1), term4(term1); Coord_type term1(0), term2(term1), term3(term1), term4(term1);
Value_type linear_int(0), gradient_int(0); Value_type linear_int(0), gradient_int(0);
typename ValueFunctor::result_type f; value_functor_result_type f;
typename GradFunctor::result_type grad; gradient_functor_result_type grad;
for(; first!=beyond; ++first) for(; first!=beyond; ++first)
{ {
@ -257,7 +287,11 @@ sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator beyond,
template < class RandomAccessIterator, class ValueFunctor, class GradFunctor, template < class RandomAccessIterator, class ValueFunctor, class GradFunctor,
class Traits, class Point_> class Traits, class Point_>
std::pair< typename ValueFunctor::result_type::first_type, bool> std::pair<
typename boost::result_of<
ValueFunctor(typename std::iterator_traits<RandomAccessIterator>::value_type::first_type)>
::type::first_type,
bool>
farin_c1_interpolation(RandomAccessIterator first, farin_c1_interpolation(RandomAccessIterator first,
RandomAccessIterator beyond, RandomAccessIterator beyond,
const typename std::iterator_traits<RandomAccessIterator>::value_type::second_type& norm, const typename std::iterator_traits<RandomAccessIterator>::value_type::second_type& norm,
@ -268,14 +302,16 @@ farin_c1_interpolation(RandomAccessIterator first,
{ {
CGAL_precondition(norm >0); CGAL_precondition(norm >0);
// the function value is available for all points typedef typename std::iterator_traits<RandomAccessIterator>::value_type::first_type arg_type;
// if a gradient value is not availble: function returns false typedef typename boost::result_of<ValueFunctor(arg_type)>::type value_functor_result_type;
typedef typename ValueFunctor::result_type::first_type Value_type; typedef typename boost::result_of<GradFunctor(arg_type)>::type gradient_functor_result_type;
typedef typename Traits::FT Coord_type; typedef typename value_functor_result_type::first_type Value_type;
typedef typename Traits::FT Coord_type;
Interpolation::internal::Extract_bare_point<Traits> cp(traits); Interpolation::internal::Extract_bare_point<Traits> cp(traits);
typename ValueFunctor::result_type f; value_functor_result_type f;
typename GradFunctor::result_type grad; gradient_functor_result_type grad;
int n = static_cast<int>(beyond - first); int n = static_cast<int>(beyond - first);
if(n == 1) if(n == 1)

View File

@ -29,14 +29,22 @@
#include <CGAL/regular_neighbor_coordinates_2.h> #include <CGAL/regular_neighbor_coordinates_2.h>
#include <CGAL/Origin.h> #include <CGAL/Origin.h>
#include <CGAL/function.h>
#include <boost/any.hpp>
#include <boost/mpl/if.hpp> #include <boost/mpl/if.hpp>
#include <boost/type_traits/is_same.hpp> #include <boost/utility/enable_if.hpp>
#include <boost/utility/result_of.hpp>
#include <iterator> #include <iterator>
#include <utility> #include <utility>
#include <vector> #include <vector>
#ifdef CGAL_CXX11
#include <type_traits>
#include <functional>
#endif
namespace CGAL { namespace CGAL {
template < class ForwardIterator, class ValueFunctor, class Traits, class Point > template < class ForwardIterator, class ValueFunctor, class Traits, class Point >
@ -45,12 +53,17 @@ sibson_gradient_fitting(ForwardIterator first,
ForwardIterator beyond, ForwardIterator beyond,
const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm, const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm,
const Point& p, const Point& p,
const typename ValueFunctor::result_type::first_type fn, const typename boost::result_of<
ValueFunctor(typename std::iterator_traits<ForwardIterator>::value_type::first_type)>
::type::first_type fn,
ValueFunctor value_function, ValueFunctor value_function,
const Traits& traits) const Traits& traits)
{ {
CGAL_precondition( first != beyond && norm != 0); CGAL_precondition( first != beyond && norm != 0);
typedef typename std::iterator_traits<ForwardIterator>::value_type::first_type arg_type;
typedef typename boost::result_of<ValueFunctor(arg_type)>::type value_functor_result_type;
typedef typename Traits::Aff_transformation_d Aff_transformation; typedef typename Traits::Aff_transformation_d Aff_transformation;
typedef typename Traits::FT Coord_type; typedef typename Traits::FT Coord_type;
typedef typename Traits::Point_d Bare_point; typedef typename Traits::Point_d Bare_point;
@ -69,7 +82,7 @@ sibson_gradient_fitting(ForwardIterator first,
typename Traits::Vector_d d = traits.construct_vector_d_object()(bp, bare_f); typename Traits::Vector_d d = traits.construct_vector_d_object()(bp, bare_f);
// compute the vector pn: // compute the vector pn:
typename ValueFunctor::result_type f = value_function(first->first); value_functor_result_type f = value_function(first->first);
CGAL_assertion(f.second); // function value of first->first is valid CGAL_assertion(f.second); // function value of first->first is valid
pn = pn + traits.construct_scaled_vector_d_object()(d, scale * (f.first - fn)); pn = pn + traits.construct_scaled_vector_d_object()(d, scale * (f.first - fn));
@ -91,7 +104,10 @@ sibson_gradient_fitting(ForwardIterator first,
ValueFunctor value_function, ValueFunctor value_function,
const Traits& traits) const Traits& traits)
{ {
typename ValueFunctor::result_type fn = value_function(p); typedef typename std::iterator_traits<ForwardIterator>::value_type::first_type arg_type;
typedef typename boost::result_of<ValueFunctor(arg_type)>::type value_functor_result_type;
value_functor_result_type fn = value_function(p);
CGAL_assertion(fn.second); CGAL_assertion(fn.second);
return sibson_gradient_fitting(first, beyond, norm, p, fn.first, value_function, traits); return sibson_gradient_fitting(first, beyond, norm, p, fn.first, value_function, traits);
@ -110,8 +126,11 @@ sibson_gradient_fitting_internal(ForwardIterator first,
const Traits& traits, const Traits& traits,
const typename Traits::Point_d& /*dummy*/) const typename Traits::Point_d& /*dummy*/)
{ {
typedef typename std::iterator_traits<ForwardIterator>::value_type::first_type arg_type;
typedef typename boost::result_of<ValueFunctor(arg_type)>::type value_functor_result_type;
const typename Traits::Point_d& bare_p = traits.construct_point_d_object()(vh->point()); const typename Traits::Point_d& bare_p = traits.construct_point_d_object()(vh->point());
typename ValueFunctor::result_type fn = value_function(bare_p); value_functor_result_type fn = value_function(bare_p);
CGAL_assertion(fn.second); CGAL_assertion(fn.second);
return sibson_gradient_fitting(first, beyond, norm, bare_p, fn.first, value_function, traits); return sibson_gradient_fitting(first, beyond, norm, bare_p, fn.first, value_function, traits);
@ -129,7 +148,10 @@ sibson_gradient_fitting_internal(ForwardIterator first,
const Traits& traits, const Traits& traits,
const typename Traits::Weighted_point_d& /*dummy*/) const typename Traits::Weighted_point_d& /*dummy*/)
{ {
typename ValueFunctor::result_type fn = value_function(vh->point()); typedef typename std::iterator_traits<ForwardIterator>::value_type::first_type arg_type;
typedef typename boost::result_of<ValueFunctor(arg_type)>::type value_functor_result_type;
value_functor_result_type fn = value_function(vh->point());
CGAL_assertion(fn.second); CGAL_assertion(fn.second);
return sibson_gradient_fitting(first, beyond, norm, vh->point(), fn.first, value_function, traits); return sibson_gradient_fitting(first, beyond, norm, vh->point(), fn.first, value_function, traits);
@ -147,15 +169,19 @@ sibson_gradient_fitting_internal(ForwardIterator first,
const Traits& traits, const Traits& traits,
VH /*dummy*/) VH /*dummy*/)
{ {
typedef typename std::iterator_traits<ForwardIterator>::value_type::first_type arg_type;
typedef typename boost::result_of<ValueFunctor(arg_type)>::type value_functor_result_type;
const typename Traits::Point_d& bare_p = traits.construct_point_d_object()(vh->point()); const typename Traits::Point_d& bare_p = traits.construct_point_d_object()(vh->point());
typename ValueFunctor::result_type fn = value_function(vh); value_functor_result_type fn = value_function(vh);
CGAL_assertion(fn.second); CGAL_assertion(fn.second);
return sibson_gradient_fitting(first, beyond, norm, bare_p, fn.first, value_function, traits); return sibson_gradient_fitting(first, beyond, norm, bare_p, fn.first, value_function, traits);
} }
template < class Tr, class OutputIterator, class OutputFunctor, template < class ValueFunctorArgType,
class Tr, class OutputIterator, class OutputFunctor,
class ValueFunctor, class CoordFunctor, class Traits > class ValueFunctor, class CoordFunctor, class Traits >
OutputIterator OutputIterator
sibson_gradient_fitting_internal(const Tr& tr, sibson_gradient_fitting_internal(const Tr& tr,
@ -170,7 +196,7 @@ sibson_gradient_fitting_internal(const Tr& tr,
typedef typename Tr::Vertex_handle Vertex_handle; typedef typename Tr::Vertex_handle Vertex_handle;
Coord_type norm; Coord_type norm;
std::vector<std::pair<typename ValueFunctor::argument_type, Coord_type> > coords; std::vector<std::pair<ValueFunctorArgType, Coord_type> > coords;
typename Tr::Finite_vertices_iterator vit = tr.finite_vertices_begin(); typename Tr::Finite_vertices_iterator vit = tr.finite_vertices_begin();
for(; vit != tr.finite_vertices_end(); ++vit) for(; vit != tr.finite_vertices_end(); ++vit)
@ -187,7 +213,7 @@ sibson_gradient_fitting_internal(const Tr& tr,
Vertex_handle(vit), Vertex_handle(vit),
value_function, value_function,
traits, traits,
typename ValueFunctor::argument_type()))); ValueFunctorArgType())));
coords.clear(); coords.clear();
} }
@ -196,12 +222,73 @@ sibson_gradient_fitting_internal(const Tr& tr,
return out; return out;
} }
// The following functions allow to fit the gradients for all points in // The following functions allow to fit the gradients for all points in
// a triangulation except the convex hull points. // a triangulation except the convex hull points.
// -> _nn2: natural_neighbor_coordinates_2 // -> _nn2: natural_neighbor_coordinates_2
// -> _rn2: regular_neighbor_coordinates_2 // -> _rn2: regular_neighbor_coordinates_2
// -> _sn2_3: surface_neighbor_coordinates_2_3 // -> _sn2_3: surface_neighbor_coordinates_2_3
// The ugly distinction below is needed to make it work with lambdas for C++11 because std::is_constructible
// is used, which is C++11 (there is a boost equivalent, but it is said (by boost) to be relying on C++11 features
// to properly work...)
#ifdef CGAL_CXX11
template < class Dt, class OutputIterator, class OutputFunctor, class ValueFunctor, class Traits >
OutputIterator
sibson_gradient_fitting_nn_2(const Dt& dt,
OutputIterator out,
OutputFunctor fct,
ValueFunctor value_function,
const Traits& traits,
// Some SFINAE to distinguish whether the argument type
// of the value functor is 'DT::Point' or 'DT::Vertex_handle'
typename boost::enable_if_c<
std::is_constructible<
std::function<boost::any(typename Dt::Point)>,
ValueFunctor
>::value>::type* = NULL)
{
typedef typename Traits::FT FT;
typedef typename Dt::Point VF_arg_type;
typedef typename std::back_insert_iterator<std::vector<
std::pair<VF_arg_type, FT> > > CoordInserter;
typedef Interpolation::internal::Extract_point_in_pair<Dt, FT> Coord_OutputFunctor;
return sibson_gradient_fitting_internal<VF_arg_type>(dt, out, fct, value_function,
natural_neighbor_coordinates_2_object<Dt,
CoordInserter,
Coord_OutputFunctor>(),
traits);
}
template < class Dt, class OutputIterator, class OutputFunctor, class ValueFunctor, class Traits >
OutputIterator
sibson_gradient_fitting_nn_2(const Dt& dt,
OutputIterator out,
OutputFunctor fct,
ValueFunctor value_function,
const Traits& traits,
typename boost::enable_if_c<
std::is_constructible<
std::function<boost::any(typename Dt::Vertex_handle)>,
ValueFunctor
>::value>::type* = NULL)
{
typedef typename Traits::FT FT;
typedef typename Dt::Vertex_handle VF_arg_type;
typedef typename std::back_insert_iterator<std::vector<
std::pair<VF_arg_type, FT> > > CoordInserter;
typedef CGAL::Identity<std::pair<VF_arg_type, FT> > Coord_OutputFunctor;
return sibson_gradient_fitting_internal<VF_arg_type>(dt, out, fct, value_function,
natural_neighbor_coordinates_2_object<Dt,
CoordInserter,
Coord_OutputFunctor>(),
traits);
}
#else // not CGAL_CXX11
template < class Dt, class OutputIterator, class OutputFunctor, class ValueFunctor, class Traits > template < class Dt, class OutputIterator, class OutputFunctor, class ValueFunctor, class Traits >
OutputIterator OutputIterator
sibson_gradient_fitting_nn_2(const Dt& dt, sibson_gradient_fitting_nn_2(const Dt& dt,
@ -224,15 +311,16 @@ sibson_gradient_fitting_nn_2(const Dt& dt,
CGAL::Identity<std::pair<VF_arg_type, FT> > CGAL::Identity<std::pair<VF_arg_type, FT> >
>::type Coord_OutputFunctor; >::type Coord_OutputFunctor;
return sibson_gradient_fitting_internal(dt, out, fct, value_function, return sibson_gradient_fitting_internal<VF_arg_type>(dt, out, fct, value_function,
natural_neighbor_coordinates_2_object<Dt, natural_neighbor_coordinates_2_object<Dt,
CoordInserter, CoordInserter,
Coord_OutputFunctor>(), Coord_OutputFunctor>(),
traits); traits);
} }
#endif // CGAL_CXX11
// Same as above but without OutputFunctor.
// Same as above but without OutputFunctor. Default to extracting the point, for backward compatibility. // Defaults to extracting the point, for backward compatibility.
template < class Dt, class OutputIterator, class ValueFunctor, class Traits > template < class Dt, class OutputIterator, class ValueFunctor, class Traits >
OutputIterator OutputIterator
sibson_gradient_fitting_nn_2(const Dt& dt, sibson_gradient_fitting_nn_2(const Dt& dt,
@ -246,6 +334,64 @@ sibson_gradient_fitting_nn_2(const Dt& dt,
return sibson_gradient_fitting_nn_2(dt, out, OutputFunctor(), value_function, traits); return sibson_gradient_fitting_nn_2(dt, out, OutputFunctor(), value_function, traits);
} }
// See above for the explanation.
#ifdef CGAL_CXX11
template < class Rt, class OutputIterator, class OutputFunctor, class ValueFunctor, class Traits >
OutputIterator
sibson_gradient_fitting_rn_2(const Rt& rt,
OutputIterator out,
OutputFunctor fct,
ValueFunctor value_function,
const Traits& traits,
// Some SFINAE to distinguish whether the argument type
// of the value functor is 'Rt::Point' (weighted point) or 'Rt::Vertex_handle'
typename boost::enable_if_c<
std::is_constructible<
std::function<boost::any(typename Rt::Point)>,
ValueFunctor
>::value>::type* = NULL)
{
typedef typename Traits::FT FT;
typedef typename Rt::Point VF_arg_type;
typedef typename std::back_insert_iterator<std::vector<
std::pair<VF_arg_type, FT> > > CoordInserter;
typedef Interpolation::internal::Extract_point_in_pair<Rt, FT> Coord_OutputFunctor;
return sibson_gradient_fitting_internal<VF_arg_type>(rt, out, fct, value_function,
regular_neighbor_coordinates_2_object<Rt,
CoordInserter,
Coord_OutputFunctor>(),
traits);
}
template < class Rt, class OutputIterator, class OutputFunctor, class ValueFunctor, class Traits >
OutputIterator
sibson_gradient_fitting_rn_2(const Rt& rt,
OutputIterator out,
OutputFunctor fct,
ValueFunctor value_function,
const Traits& traits,
typename boost::enable_if_c<
std::is_constructible<
std::function<boost::any(typename Rt::Vertex_handle)>,
ValueFunctor
>::value>::type* = NULL)
{
typedef typename Traits::FT FT;
typedef typename Rt::Vertex_handle VF_arg_type;
typedef typename std::back_insert_iterator<std::vector<
std::pair<VF_arg_type, FT> > > CoordInserter;
typedef CGAL::Identity<std::pair<VF_arg_type, FT> > Coord_OutputFunctor;
return sibson_gradient_fitting_internal<VF_arg_type>(rt, out, fct, value_function,
regular_neighbor_coordinates_2_object<Rt,
CoordInserter,
Coord_OutputFunctor>(),
traits);
}
#else // CGAL_CXX11
template < class Rt, class OutputIterator, class OutputFunctor, class ValueFunctor, class Traits > template < class Rt, class OutputIterator, class OutputFunctor, class ValueFunctor, class Traits >
OutputIterator OutputIterator
@ -269,13 +415,14 @@ sibson_gradient_fitting_rn_2(const Rt& rt,
CGAL::Identity<std::pair<VF_arg_type, FT> > CGAL::Identity<std::pair<VF_arg_type, FT> >
>::type Coord_OutputFunctor; >::type Coord_OutputFunctor;
return sibson_gradient_fitting_internal(rt, out, fct, value_function, return sibson_gradient_fitting_internal<VF_arg_type>(rt, out, fct, value_function,
regular_neighbor_coordinates_2_object<Rt, regular_neighbor_coordinates_2_object<Rt,
CoordInserter, CoordInserter,
Coord_OutputFunctor>(), Coord_OutputFunctor>(),
traits); traits);
} }
#endif
// Same as above but without OutputFunctor. Default to extracting the point, for backward compatibility. // Same as above but without OutputFunctor. Default to extracting the point, for backward compatibility.
template < class Rt, class OutputIterator, class ValueFunctor, class Traits > template < class Rt, class OutputIterator, class ValueFunctor, class Traits >

View File

@ -31,11 +31,14 @@
#include <CGAL/algorithm.h> #include <CGAL/algorithm.h>
#include <CGAL/double.h> #include <CGAL/double.h>
#include <CGAL/function_objects.h> #include <CGAL/function_objects.h>
#include <CGAL/function.h>
#include <CGAL/Origin.h> #include <CGAL/Origin.h>
#include <CGAL/point_generators_2.h> #include <CGAL/point_generators_2.h>
#include <CGAL/Random.h> #include <CGAL/Random.h>
#include <CGAL/squared_distance_2.h> #include <CGAL/squared_distance_2.h>
#include <boost/utility/result_of.hpp>
#include <iostream> #include <iostream>
#include <cassert> #include <cassert>
#include <utility> #include <utility>
@ -47,6 +50,8 @@ struct Extract_point
typedef typename Traits::Point_2 Point_2; typedef typename Traits::Point_2 Point_2;
typedef typename Traits::Weighted_point_2 Weighted_point_2; typedef typename Traits::Weighted_point_2 Weighted_point_2;
typedef typename Traits::Construct_point_2 Construct_point_2;
Extract_point(const Traits& traits = Traits()) : traits(traits) {} Extract_point(const Traits& traits = Traits()) : traits(traits) {}
const Point_2& operator()(const Point_2& p) const { return p; } const Point_2& operator()(const Point_2& p) const { return p; }
@ -56,7 +61,8 @@ struct Extract_point
} }
template <typename VH> template <typename VH>
const Point_2& operator()(const VH& vh) const { typename boost::result_of<const Construct_point_2(const Point_2&)>::type
operator()(const VH& vh) const {
return traits.construct_point_2_object()(vh->point()); return traits.construct_point_2_object()(vh->point());
} }
@ -64,6 +70,42 @@ private:
Traits traits; Traits traits;
}; };
template <typename V, typename T>
struct Value_function
{
typedef V argument_type;
typedef std::pair<T, bool> result_type;
Value_function(std::size_t i) : index(i) { }
result_type operator()(const argument_type& a) const {
return result_type(a->info()[index].value, true);
}
private:
std::size_t index;
};
template <typename V, typename G>
struct Gradient_function
: public CGAL::iterator<std::output_iterator_tag, void, void, void, void>
{
typedef V argument_type;
typedef std::pair<G, bool> result_type;
Gradient_function(std::size_t i) : index(i) { }
result_type operator()(const argument_type& a) const {
return std::make_pair(a->info()[index].gradient,
a->info()[index].gradient != CGAL::NULL_VECTOR);
}
private:
std::size_t index;
};
template < class ForwardIterator > template < class ForwardIterator >
bool test_norm(ForwardIterator first, ForwardIterator beyond, bool test_norm(ForwardIterator first, ForwardIterator beyond,
typename std::iterator_traits<ForwardIterator>::value_type::second_type norm) typename std::iterator_traits<ForwardIterator>::value_type::second_type norm)
@ -83,10 +125,10 @@ bool test_norm(ForwardIterator first, ForwardIterator beyond,
} }
} }
template < class Tr, class ForwardIterator > template < class Tr, class ForwardIterator, class Point >
bool test_barycenter(ForwardIterator first, ForwardIterator beyond, bool test_barycenter(ForwardIterator first, ForwardIterator beyond,
typename std::iterator_traits<ForwardIterator>::value_type::second_type norm, typename std::iterator_traits<ForwardIterator>::value_type::second_type norm,
const typename std::iterator_traits<ForwardIterator>::value_type::first_type& p, const Point& p,
const typename std::iterator_traits<ForwardIterator>::value_type::second_type& tolerance) const typename std::iterator_traits<ForwardIterator>::value_type::second_type& tolerance)
{ {
typedef typename Tr::Geom_traits Gt; typedef typename Tr::Geom_traits Gt;
@ -133,9 +175,13 @@ bool _test_sibson_c1_interpolation_sqrt(ForwardIterator first, ForwardIterator b
const typename std::iterator_traits<ForwardIterator>::value_type::second_type& exact_value, const typename std::iterator_traits<ForwardIterator>::value_type::second_type& exact_value,
CGAL::Field_with_sqrt_tag) CGAL::Field_with_sqrt_tag)
{ {
typename ValueFunctor::result_type res = CGAL::sibson_c1_interpolation(first, beyond, typedef typename std::iterator_traits<ForwardIterator>::value_type::first_type arg_type;
norm, p, f, typedef typename boost::result_of<ValueFunctor(arg_type)>::type value_functor_result_type;
grad_f, geom_traits);
value_functor_result_type res = CGAL::sibson_c1_interpolation(first, beyond,
norm, p, f,
grad_f, geom_traits);
return res.second && (CGAL_NTS abs(res.first-exact_value) <= tolerance); return res.second && (CGAL_NTS abs(res.first-exact_value) <= tolerance);
} }
@ -143,14 +189,18 @@ template < class ForwardIterator, class ValueFunctor, class GradFunctor, class G
bool test_interpolation_with_value(ForwardIterator first, ForwardIterator beyond, bool test_interpolation_with_value(ForwardIterator first, ForwardIterator beyond,
const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm, const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm,
const Point& p, const Point& p,
const typename ValueFunctor::result_type::first_type exact_value, const typename boost::result_of<
ValueFunctor(typename std::iterator_traits<ForwardIterator>::value_type::first_type)>
::type::first_type exact_value,
ValueFunctor f, ValueFunctor f,
GradFunctor grad_f, GradFunctor grad_f,
const Gt& geom_traits, const Gt& geom_traits,
const int& i, const int& i,
const typename std::iterator_traits<ForwardIterator>::value_type::second_type& tolerance) const typename std::iterator_traits<ForwardIterator>::value_type::second_type& tolerance)
{ {
typedef typename ValueFunctor::result_type::first_type Value_type; typedef typename std::iterator_traits<ForwardIterator>::value_type::first_type arg_type;
typedef typename boost::result_of<ValueFunctor(arg_type)>::type value_functor_result_type;
typedef typename value_functor_result_type::first_type Value_type;
if(i == 0) if(i == 0)
{ {
@ -158,8 +208,8 @@ bool test_interpolation_with_value(ForwardIterator first, ForwardIterator beyond
assert(CGAL_NTS abs(val - exact_value) <= tolerance); assert(CGAL_NTS abs(val - exact_value) <= tolerance);
} }
typename ValueFunctor::result_type res = CGAL::quadratic_interpolation(first, beyond, norm, p, f, value_functor_result_type res = CGAL::quadratic_interpolation(first, beyond, norm, p, f,
grad_f, geom_traits); grad_f, geom_traits);
assert(res.second && (CGAL_NTS abs(res.first - exact_value) <= tolerance)); assert(res.second && (CGAL_NTS abs(res.first - exact_value) <= tolerance));
if(i<2) if(i<2)
@ -185,20 +235,17 @@ bool test_interpolation_with_value(ForwardIterator first, ForwardIterator beyond
return true; return true;
} }
template < class ForwardIterator, class ValueFunctor, class GradFunctor, class Gt, class Point> template < class ForwardIterator, class ValueFunctor, class GradFunctor, class Gt, class Point, class Value_type>
bool test_interpolation(ForwardIterator first, ForwardIterator beyond, bool test_interpolation(ForwardIterator first, ForwardIterator beyond,
const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm, const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm,
const Point& p, const Point& p,
const Value_type exact_value,
ValueFunctor f, ValueFunctor f,
GradFunctor grad_f, GradFunctor grad_f,
const Gt& geom_traits, const Gt& geom_traits,
const int& i, const int& i,
const typename std::iterator_traits<ForwardIterator>::value_type::second_type& tolerance) const typename std::iterator_traits<ForwardIterator>::value_type::second_type& tolerance)
{ {
typedef typename ValueFunctor::result_type::first_type Value_type;
assert(f(p).second);
Value_type exact_value = f(p).first;
return test_interpolation_with_value(first, beyond, norm, p, exact_value, f, grad_f, geom_traits, i, tolerance); return test_interpolation_with_value(first, beyond, norm, p, exact_value, f, grad_f, geom_traits, i, tolerance);
} }
@ -222,7 +269,7 @@ void _test_interpolation_functions_2_Delaunay_without_OutputFunctor(const Dt&, c
typedef typename Gt::FT Coord_type; typedef typename Gt::FT Coord_type;
typedef typename Gt::Vector_2 Vector; typedef typename Gt::Vector_2 Vector;
typedef std::map<Point, Coord_type, typename Gt::Less_xy_2> Point_value_map ; typedef std::map<Point, Coord_type, typename Gt::Less_xy_2> Point_value_map;
typedef std::map<Point, Vector, typename Gt::Less_xy_2> Point_vector_map; typedef std::map<Point, Vector, typename Gt::Less_xy_2> Point_vector_map;
typedef std::vector<std::pair<Point, Coord_type> > Point_coordinate_vector; typedef std::vector<std::pair<Point, Coord_type> > Point_coordinate_vector;
@ -308,7 +355,8 @@ void _test_interpolation_functions_2_Delaunay_without_OutputFunctor(const Dt&, c
for(int i=0; i<3; ++i) for(int i=0; i<3; ++i)
{ {
assert(test_interpolation(coords.begin(), coords.end(), norm, points[j], assert(test_interpolation(coords.begin(), coords.end(), norm,
points[j], values[i][points[j]],
CGAL::Data_access< Point_value_map >(values[i]), CGAL::Data_access< Point_value_map >(values[i]),
CGAL::Data_access< Point_vector_map >(gradients[i]), CGAL::Data_access< Point_vector_map >(gradients[i]),
Traits(), i, tolerance)); Traits(), i, tolerance));
@ -368,7 +416,8 @@ void _test_interpolation_functions_2_Delaunay_without_OutputFunctor(const Dt&, c
for(int j=0; j<3; ++j) for(int j=0; j<3; ++j)
{ {
assert(test_interpolation(coords.begin(), coords.end(), norm, points[n/2], assert(test_interpolation(coords.begin(), coords.end(), norm,
points[n/2], values[j][points[n/2]],
CGAL::Data_access<Point_value_map>(values[j]), CGAL::Data_access<Point_value_map>(values[j]),
CGAL::Data_access<Point_vector_map>(gradients[j]), CGAL::Data_access<Point_vector_map>(gradients[j]),
Traits(), j, tolerance)); Traits(), j, tolerance));
@ -391,16 +440,18 @@ void _test_interpolation_functions_2_Delaunay_with_OutputFunctor(const Dt&, cons
typedef typename Dt::Geom_traits Gt; typedef typename Dt::Geom_traits Gt;
typedef CGAL::Interpolation_traits_2<Gt> Traits; typedef CGAL::Interpolation_traits_2<Gt> Traits;
typedef typename Dt::Vertex_handle Vertex_handle;
typedef typename Gt::FT Coord_type; typedef typename Gt::FT Coord_type;
typedef typename Dt::Point Point; typedef typename Dt::Point Point;
typedef typename Gt::Vector_2 Vector; typedef typename Gt::Vector_2 Vector;
typedef std::map<Point, Coord_type, typename Gt::Less_xy_2> Point_value_map ; typedef std::vector<std::pair<Vertex_handle, Coord_type> > Coordinate_vector;
typedef std::map<Point, Vector, typename Gt::Less_xy_2> Point_vector_map; typedef typename Coordinate_vector::const_iterator CV_cit;
typedef CGAL::Identity<std::pair<Vertex_handle, Coord_type> > Output_functor;
typedef std::vector<std::pair<Point, Coord_type> > Point_coordinate_vector; typedef std::map<Point, Coord_type> Point_value_map;
typedef typename Point_coordinate_vector::const_iterator PCV_cit; typedef std::map<Point, Vector> Point_vector_map;
typedef CGAL::Interpolation::internal::Extract_point_in_pair<Dt, Coord_type> Point_output_functor;
std::cout << "NN2: Testing random points." << std::endl; std::cout << "NN2: Testing random points." << std::endl;
@ -420,9 +471,6 @@ void _test_interpolation_functions_2_Delaunay_with_OutputFunctor(const Dt&, cons
CGAL::Random random; CGAL::Random random;
Point_value_map values[3];
Point_vector_map gradients[3];
Coord_type alpha = Coord_type(random.get_double(-max_value, max_value)), Coord_type alpha = Coord_type(random.get_double(-max_value, max_value)),
beta1 = Coord_type(random.get_double(-max_value, max_value)), beta1 = Coord_type(random.get_double(-max_value, max_value)),
beta2 = Coord_type(random.get_double(-max_value, max_value)), beta2 = Coord_type(random.get_double(-max_value, max_value)),
@ -431,67 +479,110 @@ void _test_interpolation_functions_2_Delaunay_with_OutputFunctor(const Dt&, cons
gamma3 = Coord_type(random.get_double(-max_value, max_value)); gamma3 = Coord_type(random.get_double(-max_value, max_value));
//INSERTION + DET. of GRADIENT for n DATA POINTS : //INSERTION + DET. of GRADIENT for n DATA POINTS :
for(int j=0; j<n; ++j)
{
T.insert(points[j]);
gradients[0].insert(std::make_pair(points[j], Vector(beta1, beta2)));
gradients[1].insert(std::make_pair(points[j],
Vector(beta1 + Coord_type(2)*gamma1*points[j].x(),
beta2 + Coord_type(2)*gamma1*points[j].y())));
gradients[2].insert(std::make_pair(points[j],
Vector(beta1 + Coord_type(2)*gamma1*points[j].x() + gamma3*points[j].y(),
beta2 + Coord_type(2)*gamma2*points[j].y() + gamma3*points[j].x())));
}
//DETERMINE VALUES FOR n DATA POINTS AND m RANDOM TEST POINTS: //DETERMINE VALUES FOR n DATA POINTS AND m RANDOM TEST POINTS:
for(int j=0; j<n+m; j++) Point_value_map exact_values[3];
std::map<Point, Vertex_handle> p_to_vh;
for(int j=0; j<n+m; ++j)
{ {
// linear function Vector gradient0(beta1, beta2);
values[0].insert(std::make_pair(points[j], alpha + beta1*points[j].x() + beta2*points[j].y())); Vector gradient1(beta1 + Coord_type(2)*gamma1*points[j].x(),
beta2 + Coord_type(2)*gamma1*points[j].y());
Vector gradient2(beta1 + Coord_type(2)*gamma1*points[j].x() + gamma3*points[j].y(),
beta2 + Coord_type(2)*gamma2*points[j].y() + gamma3*points[j].x());
// spherical function: Coord_type value0 = alpha + beta1*points[j].x() + beta2*points[j].y();
values[1].insert(std::make_pair(points[j], alpha + beta1*points[j].x() + Coord_type value1 = alpha + beta1*points[j].x()
beta2*points[j].y() + + beta2*points[j].y()
gamma1*points[j].x()*points[j].x()+ + gamma1*points[j].x()*points[j].x()
gamma1*points[j].y()*points[j].y())); + gamma1*points[j].y()*points[j].y();
Coord_type value2 = alpha + beta1*points[j].x()
+ beta2*points[j].y()
+ gamma1*points[j].x()*points[j].x()
+ gamma2*points[j].y()*points[j].y()
+ gamma3*points[j].x()*points[j].y();
// quadratic function if(j<n) // only insert n points
values[2].insert(std::make_pair(points[j], alpha + beta1*points[j].x() + {
beta2*points[j].y() + Vertex_handle vh = T.insert(points[j]);
gamma1*points[j].x()*points[j].x() + p_to_vh[points[j]] = vh;
gamma2*points[j].y()*points[j].y() +
gamma3*points[j].x()*points[j].y())); vh->info()[0].gradient = gradient0;
vh->info()[1].gradient = gradient1;
vh->info()[2].gradient = gradient2;
vh->info()[0].value = value0;
vh->info()[1].value = value1;
vh->info()[2].value = value2;
}
else
{
exact_values[0][points[j]] = value0;
exact_values[1][points[j]] = value1;
exact_values[2][points[j]] = value2;
}
} }
//INTERPOLATION OF RANDOM POINTS: //INTERPOLATION OF RANDOM POINTS:
Coord_type norm; Coord_type norm;
Point_coordinate_vector pt_coords; Coordinate_vector coords;
Point_output_functor pt_fct; Output_functor out_fct;
for(int j=n; j<n+m; ++j) for(int j=n; j<n+m; ++j)
{ {
CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>, Coord_type, bool> coordinate_result = CGAL::Triple<std::back_insert_iterator<Coordinate_vector>, Coord_type, bool> coordinate_result =
CGAL::natural_neighbor_coordinates_2(T, points[j], std::back_inserter(pt_coords), pt_fct); CGAL::natural_neighbor_coordinates_2(T, points[j], std::back_inserter(coords), out_fct);
assert(coordinate_result.third); assert(coordinate_result.third);
norm = coordinate_result.second; norm = coordinate_result.second;
bool is_equal = test_norm(pt_coords.begin(), pt_coords.end(), norm); bool is_equal = test_norm(coords.begin(), coords.end(), norm);
assert(norm > 0); assert(norm > 0);
assert(is_equal); assert(is_equal);
is_equal = test_barycenter<Dt>(pt_coords.begin(), pt_coords.end(), norm, points[j], tolerance); is_equal = test_barycenter<Dt>(coords.begin(), coords.end(), norm, points[j], tolerance);
assert(is_equal); assert(is_equal);
#ifndef CGAL_CFG_NO_CPP0X_LAMBDAS
assert(test_interpolation(coords.begin(), coords.end(), norm,
points[j], exact_values[0][points[j]],
[](const Vertex_handle vh) -> std::pair<Coord_type, bool> { return std::make_pair(vh->info()[0].value, true); },
[](const Vertex_handle vh) -> std::pair<Vector, bool> { return std::make_pair(vh->info()[0].gradient, true); },
Traits(), 0, tolerance));
// wrapping the lambda in a std function
CGAL::cpp11::function<std::pair<Coord_type, bool>(const Vertex_handle)> value_function_1 =
[](const Vertex_handle vh) -> std::pair<Coord_type, bool> { return std::make_pair(vh->info()[1].value, true); };
std::function<std::pair<Vector, bool>(const Vertex_handle)> gradient_function_1 =
[](const Vertex_handle vh) -> std::pair<Vector, bool> { return std::make_pair(vh->info()[1].gradient, true); };
assert(test_interpolation(coords.begin(), coords.end(), norm,
points[j], exact_values[1][points[j]],
value_function_1, gradient_function_1,
Traits(), 1, tolerance));
assert(test_interpolation(coords.begin(), coords.end(), norm,
points[j], exact_values[2][points[j]],
[](const Vertex_handle vh) -> std::pair<Coord_type, bool> { return std::make_pair(vh->info()[2].value, true); },
[](const Vertex_handle vh) -> std::pair<Vector, bool> { return std::make_pair(vh->info()[2].gradient, true); },
Traits(), 2, tolerance));
#else
for(int i=0; i<3; ++i) for(int i=0; i<3; ++i)
{ {
assert(test_interpolation(pt_coords.begin(), pt_coords.end(), norm, points[j], Value_function<Vertex_handle, Coord_type> value_function(i);
CGAL::Data_access< Point_value_map >(values[i]), Gradient_function<Vertex_handle, Vector> gradient_function(i);
CGAL::Data_access< Point_vector_map >(gradients[i]),
assert(test_interpolation(coords.begin(), coords.end(), norm,
points[j], exact_values[i][points[j]],
value_function, gradient_function,
Traits(), i, tolerance)); Traits(), i, tolerance));
} }
pt_coords.clear(); #endif
coords.clear();
} }
//TESTING THE GRADIENT APPRXIMATION METHOD: //TESTING THE GRADIENT APPRXIMATION METHOD:
@ -499,18 +590,44 @@ void _test_interpolation_functions_2_Delaunay_with_OutputFunctor(const Dt&, cons
std::cout << "Testing gradient estimation method on random points." << std::endl; std::cout << "Testing gradient estimation method on random points." << std::endl;
typedef CGAL::Interpolation_gradient_fitting_traits_2<Gt> GradTraits; typedef CGAL::Interpolation_gradient_fitting_traits_2<Gt> GradTraits;
Point_vector_map approx_gradients[2]; Point_vector_map approx_gradients[2];
#ifndef CGAL_CFG_NO_CPP0X_LAMBDAS
{
CGAL::sibson_gradient_fitting_nn_2(T,
std::inserter(approx_gradients[0], approx_gradients[0].begin()), // OutputIterator
CGAL::Interpolation::internal::Extract_point_in_pair<Dt, Vector>(), // OutputFunctor
[](const Vertex_handle vh)
-> std::pair<Coord_type, bool>
{ return std::make_pair(vh->info()[0].value, true); },
GradTraits());
std::function<std::pair<Coord_type, bool>(const Vertex_handle)> value_function_1 =
[](const Vertex_handle vh) -> std::pair<Coord_type, bool> { return std::make_pair(vh->info()[1].value, true); };
CGAL::sibson_gradient_fitting_nn_2(T,
std::inserter(approx_gradients[1], approx_gradients[1].begin()),
CGAL::Interpolation::internal::Extract_point_in_pair<Dt, Vector>(),
value_function_1,
GradTraits());
}
#else
Value_function<Vertex_handle, Coord_type> value_function_0(0);
Value_function<Vertex_handle, Coord_type> value_function_1(1);
CGAL::sibson_gradient_fitting_nn_2(T, CGAL::sibson_gradient_fitting_nn_2(T,
std::inserter(approx_gradients[0], approx_gradients[0].begin()), // OutputIterator std::inserter(approx_gradients[0], approx_gradients[0].begin()), // OutputIterator
CGAL::Interpolation::internal::Extract_point_in_pair<Dt, Vector>(), // OutputFunctor CGAL::Interpolation::internal::Extract_point_in_pair<Dt, Vector>(), // OutputFunctor
CGAL::Data_access<Point_value_map>(values[0]), // ValueFunctor value_function_0,
GradTraits()); GradTraits());
CGAL::sibson_gradient_fitting_nn_2(T, CGAL::sibson_gradient_fitting_nn_2(T,
std::inserter(approx_gradients[1], approx_gradients[1].begin()), std::inserter(approx_gradients[1], approx_gradients[1].begin()),
CGAL::Interpolation::internal::Extract_point_in_pair<Dt, Vector>(), CGAL::Interpolation::internal::Extract_point_in_pair<Dt, Vector>(),
CGAL::Data_access<Point_value_map>(values[1]), value_function_1,
GradTraits()); GradTraits());
#endif
for(int j=0; j<n; ++j) for(int j=0; j<n; ++j)
{ {
@ -518,16 +635,19 @@ void _test_interpolation_functions_2_Delaunay_with_OutputFunctor(const Dt&, cons
if(res.second) if(res.second)
{ {
Gradient_function<Vertex_handle, Vector> gradient_function_0(0);
Gradient_function<Vertex_handle, Vector> gradient_function_1(1);
// if it is the exact computation kernel: test the equality: // if it is the exact computation kernel: test the equality:
assert(tolerance > Coord_type(0) || assert(tolerance > Coord_type(0) ||
res.first == CGAL::Data_access<Point_vector_map>(gradients[0])(points[j]).first); res.first == (gradient_function_0(p_to_vh[points[j]])).first);
res = CGAL::Data_access<Point_vector_map>(approx_gradients[1])(points[j]); res = CGAL::Data_access<Point_vector_map>(approx_gradients[1])(points[j]);
// if one exists->the other must also exist // if one exists->the other must also exist
assert(res.second); assert(res.second);
assert(tolerance > Coord_type(0) || assert(tolerance > Coord_type(0) ||
res.first == CGAL::Data_access<Point_vector_map>(gradients[1])(points[j]).first); res.first == gradient_function_1(p_to_vh[points[j]]).first);
} }
else else
{ {
@ -536,31 +656,66 @@ void _test_interpolation_functions_2_Delaunay_with_OutputFunctor(const Dt&, cons
} }
//TESTING A POINT == A DATA POINT: //TESTING A POINT == A DATA POINT:
CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>, Coord_type, bool> coordinate_result = CGAL::Triple<std::back_insert_iterator<Coordinate_vector>, Coord_type, bool> coordinate_result =
CGAL::natural_neighbor_coordinates_2(T, points[n/2], std::back_inserter(pt_coords), pt_fct); CGAL::natural_neighbor_coordinates_2(T, points[n/2], std::back_inserter(coords), out_fct);
assert(coordinate_result.third); assert(coordinate_result.third);
norm = coordinate_result.second; norm = coordinate_result.second;
assert(norm == Coord_type(1)); assert(norm == Coord_type(1));
PCV_cit ci = pt_coords.begin(); CV_cit ci = coords.begin();
assert(ci->first == points[n/2]); assert(ci->first == p_to_vh[points[n/2]]);
assert(ci->second == Coord_type(1)); assert(ci->second == Coord_type(1));
ci++; ci++;
assert(ci == pt_coords.end()); assert(ci == coords.end());
#ifndef CGAL_CFG_NO_CPP0X_LAMBDAS
Value_function<Vertex_handle, Coord_type> value_function_0(0);
Value_function<Vertex_handle, Coord_type> value_function_2(2);
assert(test_interpolation(coords.begin(), coords.end(), norm,
points[n/2], value_function_0(p_to_vh[points[n/2]]).first,
[](const Vertex_handle vh) -> std::pair<Coord_type, bool> { return std::make_pair(vh->info()[0].value, true); },
[](const Vertex_handle vh) -> std::pair<Vector, bool> { return std::make_pair(vh->info()[0].gradient, true); },
Traits(), 0, tolerance));
// wrapping the lambda in a std function
CGAL::cpp11::function<std::pair<Coord_type, bool>(const Vertex_handle)> value_function_1 =
[](const Vertex_handle vh) -> std::pair<Coord_type, bool> { return std::make_pair(vh->info()[1].value, true); };
std::function<std::pair<Vector, bool>(const Vertex_handle)> gradient_function_1 =
[](const Vertex_handle vh) -> std::pair<Vector, bool> { return std::make_pair(vh->info()[1].gradient, true); };
assert(test_interpolation(coords.begin(), coords.end(), norm,
points[n/2], value_function_1(p_to_vh[points[n/2]]).first,
value_function_1, gradient_function_1,
Traits(), 1, tolerance));
assert(test_interpolation(coords.begin(), coords.end(), norm,
points[n/2], value_function_2(p_to_vh[points[n/2]]).first,
[](const Vertex_handle vh) -> std::pair<Coord_type, bool> { return std::make_pair(vh->info()[2].value, true); },
[](const Vertex_handle vh) -> std::pair<Vector, bool> { return std::make_pair(vh->info()[2].gradient, true); },
Traits(), 2, tolerance));
#else
for(int j=0; j<3; ++j) for(int j=0; j<3; ++j)
{ {
assert(test_interpolation(pt_coords.begin(), pt_coords.end(), norm, points[n/2], Value_function<Vertex_handle, Coord_type> value_function(j);
CGAL::Data_access<Point_value_map>(values[j]), Gradient_function<Vertex_handle, Vector> gradient_function(j);
CGAL::Data_access<Point_vector_map>(gradients[j]),
assert(test_interpolation(coords.begin(), coords.end(), norm,
points[n/2], value_function(p_to_vh[points[n/2]]).first,
value_function, gradient_function,
Traits(), j, tolerance)); Traits(), j, tolerance));
} }
pt_coords.clear(); #endif
coords.clear();
} }
template <class Rt> template <class Rt>
void _test_interpolation_functions_2_regular_without_OutputFunctor(const Rt&, const typename Rt::Geom_traits::FT& tolerance) void _test_interpolation_functions_2_regular_without_OutputFunctor(const Rt&, const typename Rt::Geom_traits::FT& tolerance)
{ {
std::cout << "Testing backward compatibility..." << std::endl;
CGAL::Set_ieee_double_precision pfr; CGAL::Set_ieee_double_precision pfr;
Rt T; Rt T;
@ -577,12 +732,12 @@ void _test_interpolation_functions_2_regular_without_OutputFunctor(const Rt&, co
typedef typename Gt::FT Coord_type; typedef typename Gt::FT Coord_type;
typedef typename Gt::Vector_2 Vector; typedef typename Gt::Vector_2 Vector;
typedef std::map<Weighted_point, Coord_type> Point_value_map ; typedef std::map<Weighted_point, Coord_type> Point_value_map;
typedef std::map<Weighted_point, Vector> Point_vector_map; typedef std::map<Weighted_point, Vector> Point_vector_map;
typedef std::vector<std::pair<Weighted_point, Coord_type> > Point_coordinate_vector; typedef std::vector<std::pair<Weighted_point, Coord_type> > Point_coordinate_vector;
std::cout << "NN2: Testing random points." << std::endl; std::cout << "RN2: Testing random points." << std::endl;
// test random points in a square of length r: // test random points in a square of length r:
std::vector<Weighted_point> points; std::vector<Weighted_point> points;
@ -674,7 +829,8 @@ void _test_interpolation_functions_2_regular_without_OutputFunctor(const Rt&, co
for(int i=0; i<3; ++i) for(int i=0; i<3; ++i)
{ {
assert(test_interpolation(coords.begin(), coords.end(), norm, points[j], assert(test_interpolation(coords.begin(), coords.end(), norm,
points[j], values[i][points[j]],
CGAL::Data_access< Point_value_map >(values[i]), CGAL::Data_access< Point_value_map >(values[i]),
CGAL::Data_access< Point_vector_map >(gradients[i]), CGAL::Data_access< Point_vector_map >(gradients[i]),
Traits(), i, tolerance)); Traits(), i, tolerance));
@ -735,7 +891,8 @@ void _test_interpolation_functions_2_regular_without_OutputFunctor(const Rt&, co
for(int j=0; j<3; ++j) for(int j=0; j<3; ++j)
{ {
assert(test_interpolation(coords.begin(), coords.end(), norm, points[n/2], assert(test_interpolation(coords.begin(), coords.end(), norm,
points[n/2], values[j][points[n/2]],
CGAL::Data_access<Point_value_map>(values[j]), CGAL::Data_access<Point_value_map>(values[j]),
CGAL::Data_access<Point_vector_map>(gradients[j]), CGAL::Data_access<Point_vector_map>(gradients[j]),
Traits(), j, tolerance)); Traits(), j, tolerance));
@ -771,7 +928,7 @@ void _test_interpolation_functions_2_regular_with_OutputFunctor(const Rt&, const
typedef typename Rt::Vertex_handle Vertex_handle; typedef typename Rt::Vertex_handle Vertex_handle;
// These are the values at points which won't be inserted in the triangulation // These are the values at points which won't be inserted in the triangulation
typedef std::map<Weighted_point, Coord_type> Point_value_map ; typedef std::map<Weighted_point, Coord_type> Point_value_map;
typedef std::map<Vertex_handle, Coord_type> Vertex_value_map; typedef std::map<Vertex_handle, Coord_type> Vertex_value_map;
typedef std::map<Vertex_handle, Vector> Vertex_vector_map; typedef std::map<Vertex_handle, Vector> Vertex_vector_map;
@ -782,7 +939,7 @@ void _test_interpolation_functions_2_regular_with_OutputFunctor(const Rt&, const
Identity_output_functor vh_fct; Identity_output_functor vh_fct;
std::cout << "NN2: Testing random points." << std::endl; std::cout << "RN2: Testing random points." << std::endl;
// test random points in a square of length r: // test random points in a square of length r:
std::vector<Weighted_point> points; std::vector<Weighted_point> points;
@ -979,8 +1136,8 @@ void _test_interpolation_functions_2_regular_with_OutputFunctor(const Rt&, const
std::pair<FT, bool> ev = CGAL::Data_access<Vertex_value_map>(values[j])(vh); std::pair<FT, bool> ev = CGAL::Data_access<Vertex_value_map>(values[j])(vh);
assert(ev.second); assert(ev.second);
assert(test_interpolation_with_value(vh_coords.begin(), vh_coords.end(), norm, vh->point(), assert(test_interpolation_with_value(vh_coords.begin(), vh_coords.end(), norm,
ev.first /*exact value*/, vh->point(), ev.first /*exact value*/,
CGAL::Data_access<Vertex_value_map>(values[j]), CGAL::Data_access<Vertex_value_map>(values[j]),
CGAL::Data_access<Vertex_vector_map>(gradients[j]), CGAL::Data_access<Vertex_vector_map>(gradients[j]),
Traits(), j, tolerance)); Traits(), j, tolerance));

View File

@ -21,44 +21,59 @@
#include <CGAL/Exact_predicates_exact_constructions_kernel.h> #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/array.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Delaunay_triangulation_2.h> #include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Regular_triangulation_2.h> #include <CGAL/Regular_triangulation_2.h>
#include <CGAL/Origin.h>
#include <iostream> #include <iostream>
typedef CGAL::Exact_predicates_exact_constructions_kernel K; typedef CGAL::Exact_predicates_exact_constructions_kernel EPECK;
typedef CGAL::Delaunay_triangulation_2<K> Dt; typedef CGAL::Exact_predicates_inexact_constructions_kernel EPICK;
typedef CGAL::Regular_triangulation_2<K> Rt;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K2; template <typename V, typename G>
typedef CGAL::Delaunay_triangulation_2<K2> Dt2; struct Value_and_gradient
typedef CGAL::Regular_triangulation_2<K2> Rt2; {
Value_and_gradient() : value(), gradient(CGAL::NULL_VECTOR) {}
V value;
G gradient;
};
template<typename Kernel>
void test_interpolation_functions()
{
// For the Delaunay triangulation, values and gradients (three different data sets)
// are stored directly in the vertices
typedef typename Kernel::FT Coord_type;
typedef typename Kernel::Vector_2 Vector;
typedef CGAL::Triangulation_vertex_base_with_info_2<
CGAL::cpp11::array<
Value_and_gradient<Coord_type, Vector>, 3>,
Kernel> Vb;
typedef CGAL::Triangulation_data_structure_2<Vb> Tds;
typedef CGAL::Delaunay_triangulation_2<Kernel, Tds> Delaunay_triangulation;
typedef CGAL::Regular_triangulation_2<Kernel> Regular_triangulation;
std::cout << "Testing interpolation functions with 2D NN neighbors " << std::endl;
_test_interpolation_functions_2_Delaunay(Delaunay_triangulation(), Coord_type(1e-10));
std::cout << "Testing interpolation functions with 2D RN neighbors " << std::endl;
_test_interpolation_functions_2_regular(Regular_triangulation(), Coord_type(1e-10));
}
int main() int main()
{ {
std::cout << "Testing interpolation functions with 2D NN neighbors " std::cout << "--------------------------------------------" << std::endl;
<< std::endl; std::cout << "Testing with EPECK" << std::endl;
std::cout << " using Exact_predicates_exact_constructions_kernel: " test_interpolation_functions<EPECK>();
<< std::endl ;
_test_interpolation_functions_2_Delaunay(Dt(), K::FT(1e-10));
std::cout << "Testing interpolation functions with 2D NN neighbors " std::cout << "--------------------------------------------" << std::endl;
<< std::endl; std::cout << "Testing with EPICK" << std::endl;
std::cout << " using Exact_predicates_inexact_constructions_kernel: " test_interpolation_functions<EPICK>();
<< std::endl ;
_test_interpolation_functions_2_Delaunay(Dt2(), K2::FT(1e-10));
std::cout << "Testing interpolation functions with 2D RN neighbors "
<< std::endl;
std::cout << " using Exact_predicates_exact_constructions_kernel: "
<< std::endl ;
_test_interpolation_functions_2_regular(Rt(), K::FT(1e-10));
std::cout << "Testing interpolation functions with 2D RN neighbors "
<< std::endl;
std::cout << " using Exact_predicates_inexact_constructions_kernel: "
<< std::endl ;
_test_interpolation_functions_2_regular(Rt2(), K2::FT(1e-10));
std::cout << "test_interpolation_functions_2 is finished" << std::endl; std::cout << "test_interpolation_functions_2 is finished" << std::endl;