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

Fixed trailing whitespace and (some of the) inconsistent indentation 

No real changes.
This commit is contained in:
Mael Rouxel-Labbé 2018-01-10 15:51:05 +01:00
parent 118e5dc9c3
commit 67f99cc53d
17 changed files with 592 additions and 591 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
Interpolation/examples/Interpolation/interpolation_2.cpp
View File

@ -49,7 +49,7 @@ int main()
// Create n+m-4 points within a disc of radius 2
double r_d = 3;
CGAL::Random rng(1513114263);
CGAL::Random_points_in_disc_2<Point> g(r_d,rng );
CGAL::cpp11::copy_n( g, n+m, std::back_inserter(points));
@ -114,7 +114,7 @@ int main()
l_value = CGAL::linear_interpolation(coords.begin(), coords.end(),
norm,
CGAL::Data_access<Point_value_map>(values));
error = CGAL_NTS abs(l_value - exact_value);
l_total += error;
if (error > l_max) l_max = error;
@ -208,6 +208,5 @@ int main()
<< CGAL::to_double(ssquare_total)/n << " max "
<< CGAL::to_double(ssquare_max) << std::endl;
std::cout << "done" << std::endl;
return 0;
return EXIT_SUCCESS;
}

239
Interpolation/examples/Interpolation/interpolation_vertex_with_info_2.cpp
View File

@ -16,72 +16,58 @@
#include <cassert>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::FT Coord_type;
typedef K::Vector_2 Vector;
typedef K::Point_2 Point;
typedef K::FT Coord_type;
typedef K::Vector_2 Vector;
typedef K::Point_2 Point;
template <typename V, typename G>
struct Value_and_gradient {
Value_and_gradient()
: value(), gradient(CGAL::NULL_VECTOR)
{}
struct Value_and_gradient
{
Value_and_gradient() : value(), gradient(CGAL::NULL_VECTOR) {}
V value;
G gradient;
};
typedef CGAL::Triangulation_vertex_base_with_info_2<Value_and_gradient<Coord_type,Vector>, K> Vb;
typedef CGAL::Triangulation_data_structure_2<Vb> Tds;
typedef CGAL::Delaunay_triangulation_2<K,Tds> Delaunay_triangulation;
typedef Delaunay_triangulation::Vertex_handle Vertex_handle;
typedef CGAL::Interpolation_traits_2<K> Traits;
typedef CGAL::Triangulation_vertex_base_with_info_2<
Value_and_gradient<Coord_type, Vector>, K> Vb;
typedef CGAL::Triangulation_data_structure_2<Vb> Tds;
typedef CGAL::Delaunay_triangulation_2<K, Tds> Delaunay_triangulation;
typedef Delaunay_triangulation::Vertex_handle Vertex_handle;
typedef CGAL::Interpolation_traits_2<K> Traits;
typedef std::vector< std::pair<Point, Coord_type> > Coordinate_vector;
typedef std::vector< std::pair<Point, Coord_type> > Coordinate_vector;
template <typename V, typename T>
struct Function_value {
typedef V argument_type;
typedef std::pair<T, bool> result_type;
result_type operator()(const argument_type& a)const
{
struct Function_value
{
typedef V argument_type;
typedef std::pair<T, bool> result_type;
result_type operator()(const argument_type& a) const {
return result_type(a->info().value, true);
}
};
template <typename V, typename G>
struct Function_gradient
: public std::iterator<std::output_iterator_tag,void,void,void,void> {
typedef V argument_type;
typedef std::pair<G,bool> result_type;
result_type
operator()(const argument_type& a)const
{
return std::make_pair(a->info().gradient,a->info().gradient != CGAL::NULL_VECTOR) ;
: public std::iterator<std::output_iterator_tag, void, void, void, void>
{
typedef V argument_type;
typedef std::pair<G, bool> result_type;
result_type operator()(const argument_type& a) const {
return std::make_pair(a->info().gradient, a->info().gradient != CGAL::NULL_VECTOR);
}
const Function_gradient& operator=(const std::pair<V, G>& p) const
{
p.first->info().gradient = p.second;
return *this;
}
const Function_gradient& operator++(int) const
{
return *this;
}
const Function_gradient& operator*() const
{
const Function_gradient& operator=(const std::pair<V, G>& p) const {
p.first->info().gradient = p.second;
return *this;
}
const Function_gradient& operator++(int) const { return *this; }
const Function_gradient& operator*() const { return *this; }
};
int main()
@ -103,24 +89,25 @@ int main()
// Create n+m-4 points within a disc of radius 2
double r_d = 3;
CGAL::Random rng(1513114263);
CGAL::Random_points_in_disc_2<Point> g(r_d,rng );
CGAL::Random_points_in_disc_2<Point> g(r_d, rng);
CGAL::cpp11::copy_n( g, n+m, std::back_inserter(points));
Delaunay_triangulation T;
Function_value<Vertex_handle,Coord_type> function_value;
Function_gradient<Vertex_handle,Vector> function_gradient;
Function_value<Vertex_handle, Coord_type> function_value;
Function_gradient<Vertex_handle, Vector> function_gradient;
//parameters for quadratic function:
Coord_type alpha = Coord_type(1.0),
beta1 = Coord_type(2.0),
beta2 = Coord_type(1.0),
gamma1 = Coord_type(0.3),
gamma2 = Coord_type(0.0),
gamma3 = Coord_type(0.0),
gamma4 = Coord_type(0.3);
beta1 = Coord_type(2.0),
beta2 = Coord_type(1.0),
gamma1 = Coord_type(0.3),
gamma2 = Coord_type(0.0),
gamma3 = Coord_type(0.0),
gamma4 = Coord_type(0.3);
for(int j=0; j<n ; j++){
for(int j=0; j<n ; j++)
{
Vertex_handle vh = T.insert(points[j]);
//determine function value/gradient:
@ -129,8 +116,8 @@ int main()
Coord_type value = alpha + beta1*x + beta2*y + gamma1*(x*x) +
gamma4*(y*y) + (gamma2+ gamma3) *(x*y);
Vector gradient(beta1+ (gamma2+ gamma3)*y + Coord_type(2)*(gamma1*x),
beta2+ (gamma2+ gamma3)*x + Coord_type(2)*(gamma4*y));
Vector gradient(beta1+ (gamma2 + gamma3)*y + Coord_type(2)*(gamma1*x),
beta2+ (gamma2 + gamma3)*x + Coord_type(2)*(gamma4*y));
vh->info().value = value;
vh->info().gradient = gradient;
}
@ -138,15 +125,16 @@ int main()
//variables for statistics:
std::pair<Coord_type, bool> res;
Coord_type error, l_total = Coord_type(0),
q_total(l_total), f_total(l_total), s_total(l_total),
ssquare_total(l_total), l_max(l_total),
q_max(l_total), f_max(l_total), s_max(l_total),
ssquare_max(l_total),
total_value(l_total), l_value(l_total);
q_total(l_total), f_total(l_total), s_total(l_total),
ssquare_total(l_total), l_max(l_total),
q_max(l_total), f_max(l_total), s_max(l_total),
ssquare_max(l_total),
total_value(l_total), l_value(l_total);
int failure(0);
//interpolation + error statistics
for(int i=n;i<n+m;i++){
for(int i=n; i<n+m; i++)
{
Coord_type x(points[i].x());
Coord_type y(points[i].y());
@ -155,32 +143,23 @@ int main()
total_value += exact_value;
//Coordinate_vector:
std::vector< std::pair< Vertex_handle, Coord_type > > coords;
typedef CGAL::Identity<std::pair< Vertex_handle, Coord_type > > Identity;
Coord_type norm =
CGAL::natural_neighbor_coordinates_2(T,
points[i],
std::back_inserter(coords),
Identity()).second;
assert(norm>0);
std::vector< std::pair<Vertex_handle, Coord_type> > coords;
typedef CGAL::Identity<std::pair< Vertex_handle, Coord_type> > Identity;
Coord_type norm = CGAL::natural_neighbor_coordinates_2(T,
points[i],
std::back_inserter(coords),
Identity()).second;
assert(norm > 0);
//linear interpolant:
l_value =
CGAL::linear_interpolation(coords.begin(), coords.end(),
norm,
function_value);
l_value = CGAL::linear_interpolation(coords.begin(), coords.end(),
norm, function_value);
error = CGAL_NTS abs(l_value - exact_value);
l_total += error;
if (error > l_max) l_max = error;
//Farin interpolant:
res = CGAL::farin_c1_interpolation(coords.begin(),
coords.end(), norm,points[i],
@ -189,59 +168,79 @@ int main()
Traits());
if(res.second){
if(res.second)
{
error = CGAL_NTS abs(res.first - exact_value);
f_total += error;
if (error > f_max) f_max = error;
} else ++failure;
if (error > f_max)
f_max = error;
}
else
{
++failure;
}
//quadratic interpolant:
res = CGAL::quadratic_interpolation(coords.begin(), coords.end(),
norm,points[i],
function_value,
function_gradient,
Traits());
if(res.second){
res = CGAL::quadratic_interpolation(coords.begin(), coords.end(),
norm, points[i],
function_value,
function_gradient,
Traits());
if(res.second)
{
error = CGAL_NTS abs(res.first - exact_value);
q_total += error;
if (error > q_max) q_max = error;
} else ++failure;
if (error > q_max)
q_max = error;
}
else
{
++failure;
}
//Sibson interpolant: version without sqrt:
res = CGAL::sibson_c1_interpolation_square(coords.begin(),
coords.end(), norm,
points[i],
function_value,
function_gradient,
Traits());
res = CGAL::sibson_c1_interpolation_square(coords.begin(),
coords.end(), norm,
points[i],
function_value,
function_gradient,
Traits());
//error statistics
if(res.second){
if(res.second)
{
error = CGAL_NTS abs(res.first - exact_value);
ssquare_total += error;
if (error > ssquare_max) ssquare_max = error;
} else ++failure;
if (error > ssquare_max)
ssquare_max = error;
}
else
{
++failure;
}
//with sqrt(the traditional):
res = CGAL::sibson_c1_interpolation(coords.begin(),
coords.end(), norm,
points[i],
function_value,
function_gradient,
Traits());
res = CGAL::sibson_c1_interpolation(coords.begin(),
coords.end(), norm,
points[i],
function_value,
function_gradient,
Traits());
//error statistics
if(res.second){
if(res.second)
{
error = CGAL_NTS abs(res.first - exact_value);
s_total += error;
if (error > s_max) s_max = error;
} else ++failure;
if (error > s_max)
s_max = error;
}
else
{
++failure;
}
}
/************** end of Interpolation: dump statistics **************/
std::cout << "Result: -----------------------------------" << std::endl;
std::cout << "Interpolation of '" << alpha <<" + "
@ -270,7 +269,5 @@ int main()
<< CGAL::to_double(ssquare_total)/n << " max "
<< CGAL::to_double(ssquare_max) << std::endl;
std::cout << "done" << std::endl;
return 0;
return EXIT_SUCCESS;
}

33
Interpolation/examples/Interpolation/linear_interpolation_2.cpp
View File

@ -6,42 +6,39 @@
#include <CGAL/interpolation_functions.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
typedef CGAL::Interpolation_traits_2<K> Traits;
typedef K::FT Coord_type;
typedef K::Point_2 Point;
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
typedef CGAL::Interpolation_traits_2<K> Traits;
typedef K::FT Coord_type;
typedef K::Point_2 Point;
int main()
{
Delaunay_triangulation T;
std::map<Point, Coord_type, K::Less_xy_2> function_values;
typedef CGAL::Data_access< std::map<Point, Coord_type, K::Less_xy_2 > >
Value_access;
typedef std::map<Point, Coord_type, K::Less_xy_2> Coord_map;
typedef CGAL::Data_access<Coord_map> Value_access;
Coord_map function_values;
Coord_type a(0.25), bx(1.3), by(-0.7);
for (int y=0 ; y<3 ; y++){
for (int x=0 ; x<3 ; x++){
K::Point_2 p(x,y);
T.insert(p);
function_values.insert(std::make_pair(p,a + bx* x+ by*y));
function_values.insert(std::make_pair(p, a + bx*x + by*y));
}
}
//coordinate computation
K::Point_2 p(1.3,0.34);
std::vector< std::pair< Point, Coord_type > > coords;
Coord_type norm =
CGAL::natural_neighbor_coordinates_2
(T, p,std::back_inserter(coords)).second;
K::Point_2 p(1.3, 0.34);
std::vector<std::pair<Point, Coord_type> > coords;
Coord_type res = CGAL::linear_interpolation(coords.begin(), coords.end(),
norm,
Coord_type norm = CGAL::natural_neighbor_coordinates_2(T, p, std::back_inserter(coords)).second;
Coord_type res = CGAL::linear_interpolation(coords.begin(), coords.end(), norm,
Value_access(function_values));
std::cout << "Tested interpolation on " << p << " interpolation: "
<< res << " exact: " << a + bx* p.x()+ by* p.y()<< std::endl;
<< res << " exact: " << a + bx*p.x() + by*p.y() << std::endl;
std::cout << "done" << std::endl;
return 0;
return EXIT_SUCCESS;
}

36
Interpolation/examples/Interpolation/nn_coordinates_2.cpp
View File

@ -1,36 +1,40 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <iostream>
#include <iterator>
#include <utility>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
typedef std::vector< std::pair< K::Point_2, K::FT > > Point_coordinate_vector;
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
typedef std::vector< std::pair<K::Point_2, K::FT> > Point_coordinate_vector;
int main()
{
Delaunay_triangulation dt;
for (int y=0 ; y<3 ; y++)
for (int x=0 ; x<3 ; x++)
dt.insert(K::Point_2(x,y));
for (int y=0; y<3; y++)
for (int x=0; x<3; x++)
dt.insert(K::Point_2(x, y));
//coordinate computation
K::Point_2 p(1.2, 0.7);
Point_coordinate_vector coords;
CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>,
K::FT, bool> result =
CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>, K::FT, bool> result =
CGAL::natural_neighbor_coordinates_2(dt, p, std::back_inserter(coords));
if(!result.third){
std::cout << "The coordinate computation was not successful."
<< std::endl;
std::cout << "The point (" <<p << ") lies outside the convex hull."
<< std::endl;
if(!result.third)
{
std::cout << "The coordinate computation was not successful." << std::endl;
std::cout << "The point (" << p << ") lies outside the convex hull." << std::endl;
}
K::FT norm = result.second;
K::FT norm = result.second;
std::cout << "Coordinate computation successful." << std::endl;
std::cout << "Normalization factor: " <<norm << std::endl;
std::cout << "done" << std::endl;
return 0;
std::cout << "Normalization factor: " << norm << std::endl;
return EXIT_SUCCESS;
}

46
Interpolation/examples/Interpolation/nn_coordinates_3.cpp
View File

@ -1,27 +1,33 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/natural_neighbor_coordinates_3.h>
#include <fstream>
#include <iostream>
#include <iterator>
#include <utility>
#include <vector>
typedef double NT; //Number Type
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point3;
typedef K::Vector_3 Vector3;
typedef K::Sphere_3 Sphere_3;
typedef K::Point_3 Point3;
typedef K::Vector_3 Vector3;
typedef K::Sphere_3 Sphere_3;
typedef CGAL::Delaunay_triangulation_3<K, CGAL::Fast_location> Dh;
typedef Dh::Facet Facet;
typedef Dh::Vertex_handle Vertex_handle;
typedef Dh::Cell_handle Cell_handle;
typedef Dh::Finite_vertices_iterator Finite_vertices_iterator;
typedef Dh::Vertex_iterator Vertex_iterator;
typedef Dh::Facet_circulator Facet_circulator;
typedef Dh::Cell_iterator Cell_iterator;
typedef Dh::Facet Facet;
typedef Dh::Vertex_handle Vertex_handle;
typedef Dh::Cell_handle Cell_handle;
typedef Dh::Finite_vertices_iterator Finite_vertices_iterator;
typedef Dh::Vertex_iterator Vertex_iterator;
typedef Dh::Facet_circulator Facet_circulator;
typedef Dh::Cell_iterator Cell_iterator;
typedef K::Construct_circumcenter_3 Construct_circumcenter_3;
typedef K::Construct_circumcenter_3 Construct_circumcenter_3;
int main()
{
@ -35,11 +41,11 @@ int main()
Point3 pp[3];
std::cout << "Consider the natural coordinates of P1, P2 and P3 with regard to the triangulation of data/points3 " << std::endl;
pp[0]=Point3(106.55,112.57,110.4); //inside data/points3 convex hull
pp[0] = Point3(106.55,112.57,110.4); //inside data/points3 convex hull
std::cout << "P1 is inside the convex hull" << std::endl;
pp[1]=Point3(250,100,140); //on data/points3 convex hull
pp[1] = Point3(250,100,140); //on data/points3 convex hull
std::cout << "P2 is on a vertex of the triangulation" << std::endl;
pp[2]=Point3(0,0,0); //outside data/points3 convex hull
pp[2] = Point3(0,0,0); //outside data/points3 convex hull
std::cout << "P2 is outside the convex hull" << std::endl;
for(int ii=0; ii<3; ++ii)
@ -48,9 +54,9 @@ int main()
std::vector< std::pair< Vertex_handle,NT> > coor_sibson;
NT norm_coeff_laplace, norm_coeff_sibson;
std::cout << "Point P"<< ii+1 << " : "<<pp[ii].x() << " "
<< pp[ii].y() << " "
<< pp[ii].z() << std::endl;
std::cout << "Point P" << ii+1 << " : " << pp[ii].x() << " "
<< pp[ii].y() << " "
<< pp[ii].z() << std::endl;
laplace_natural_neighbor_coordinates_3(triangulation,pp[ii],
std::back_inserter(coor_laplace),
@ -67,6 +73,7 @@ int main()
sibson_natural_neighbor_coordinates_3(triangulation,pp[ii],
std::back_inserter(coor_sibson),
norm_coeff_sibson);
std::cout << "Linear combination of natural neighbors with Sibson natural coordinates" << std::endl;
std::cout << " + correctness (ensured only with an exact number type)" << std::endl;
std::cout << is_correct_natural_neighborhood(triangulation,pp[ii],
@ -76,6 +83,5 @@ int main()
<< std::endl;
}
std::cout << "done" << std::endl;
return 0;
return EXIT_SUCCESS;
}

36
Interpolation/examples/Interpolation/rn_coordinates_2.cpp
View File

@ -1,40 +1,42 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_2.h>
#include <CGAL/regular_neighbor_coordinates_2.h>
#include <iostream>
#include <iterator>
#include <vector>
#include <utility>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Regular_triangulation_2<K> Regular_triangulation;
typedef Regular_triangulation::Bare_point Bare_point;
typedef Regular_triangulation::Weighted_point Weighted_point;
typedef std::vector< std::pair< Weighted_point, K::FT > > Point_coordinate_vector;
typedef CGAL::Regular_triangulation_2<K> Regular_triangulation;
typedef Regular_triangulation::Bare_point Bare_point;
typedef Regular_triangulation::Weighted_point Weighted_point;
typedef std::vector<std::pair<Weighted_point, K::FT> > Point_coordinate_vector;
int main()
{
Regular_triangulation rt;
for (int y=0 ; y<3 ; y++)
for (int x=0 ; x<3 ; x++)
for (int y=0; y<3; y++)
for (int x=0; x<3; x++)
rt.insert(Weighted_point(Bare_point(x,y), 0));
//coordinate computation
Weighted_point wp(Bare_point(1.2, 0.7),2);
Point_coordinate_vector coords;
CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>,
K::FT, bool> result =
Weighted_point wp(Bare_point(1.2, 0.7), 2);
Point_coordinate_vector coords;
CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>, K::FT, bool> result =
CGAL::regular_neighbor_coordinates_2(rt, wp, std::back_inserter(coords));
if(!result.third){
std::cout << "The coordinate computation was not successful."
<< std::endl;
std::cout << "The point (" <<wp.point() << ") lies outside the convex hull."
<< std::endl;
std::cout << "The coordinate computation was not successful." << std::endl;
std::cout << "The point (" <<wp.point() << ") lies outside the convex hull." << std::endl;
}
K::FT norm = result.second;
std::cout << "Coordinate computation successful." << std::endl;
std::cout << "Normalization factor: " <<norm << std::endl;
std::cout << "Normalization factor: " << norm << std::endl;
std::cout << "done" << std::endl;
return 0;
return EXIT_SUCCESS;
}

50
Interpolation/examples/Interpolation/sibson_interpolation_2.cpp
View File

@ -6,14 +6,20 @@
#include <CGAL/sibson_gradient_fitting.h>
#include <CGAL/interpolation_functions.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
#include <iostream>
#include <iterator>
#include <map>
#include <utility>
#include <vector>
typedef K::FT Coord_type;
typedef K::Point_2 Point;
typedef std::map<Point, Coord_type, K::Less_xy_2> Point_value_map ;
typedef std::map<Point, K::Vector_2 , K::Less_xy_2 > Point_vector_map;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
typedef K::FT Coord_type;
typedef K::Point_2 Point;
typedef std::map<Point, Coord_type, K::Less_xy_2> Point_value_map ;
typedef std::map<Point, K::Vector_2 , K::Less_xy_2> Point_vector_map;
int main()
{
@ -24,16 +30,16 @@ int main()
//parameters for spherical function:
Coord_type a(0.25), bx(1.3), by(-0.7), c(0.2);
for (int y=0 ; y<4 ; y++){
for (int x=0 ; x<4 ; x++){
for (int y=0; y<4; y++) {
for (int x=0; x<4; x++) {
K::Point_2 p(x,y);
T.insert(p);
function_values.insert(std::make_pair(p,a + bx* x+ by*y + c*(x*x+y*y)));
}
}
sibson_gradient_fitting_nn_2(T,std::inserter(function_gradients,
function_gradients.begin()),
sibson_gradient_fitting_nn_2(T, std::inserter(function_gradients,
function_gradients.begin()),
CGAL::Data_access<Point_value_map>(function_values),
Traits());
@ -41,8 +47,8 @@ int main()
{
std::cout << it->first << " " << it->second << std::endl;
}
//coordinate computation
K::Point_2 p(1.6,1.4);
// coordinate computation
K::Point_2 p(1.6, 1.4);
std::vector< std::pair< Point, Coord_type > > coords;
Coord_type norm = CGAL::natural_neighbor_coordinates_2(T, p, std::back_inserter
(coords)).second;
@ -50,23 +56,21 @@ int main()
//Sibson interpolant: version without sqrt:
std::pair<Coord_type, bool> res =
CGAL::sibson_c1_interpolation_square(
coords.begin(),
coords.end(),norm,p,
CGAL::Data_access<Point_value_map>(function_values),
CGAL::Data_access<Point_vector_map>(function_gradients),
Traits());
CGAL::sibson_c1_interpolation_square(coords.begin(),
coords.end(),norm,p,
CGAL::Data_access<Point_value_map>(function_values),
CGAL::Data_access<Point_vector_map>(function_gradients),
Traits());
if(res.second)
std::cout << "Tested interpolation on " << p
<< " interpolation: " << res.first << " exact: "
<< a + bx * p.x()+ by * p.y()+ c*(p.x()*p.x()+p.y()*p.y())
<< a + bx*p.x() + by*p.y() + c*(p.x()*p.x()+p.y()*p.y())
<< std::endl;
else
std::cout << "C^1 Interpolation not successful." << std::endl
<< " not all function_gradients are provided." << std::endl
<< " not all function_gradients are provided." << std::endl
<< " You may resort to linear interpolation." << std::endl;
std::cout << "done" << std::endl;
return 0;
return EXIT_SUCCESS;
}

64
Interpolation/examples/Interpolation/sibson_interpolation_rn_2.cpp
View File

@ -1,28 +1,34 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/Interpolation_gradient_fitting_traits_2.h>
#include <CGAL/sibson_gradient_fitting.h>
#include <CGAL/interpolation_functions.h>
#include <CGAL/Regular_triangulation_2.h>
#include <iostream>
#include <iterator>
#include <map>
#include <utility>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Regular_triangulation_2<K> Regular_triangulation;
typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
typedef CGAL::Regular_triangulation_2<K> Regular_triangulation;
typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
typedef K::FT Coord_type;
typedef K::Weighted_point_2 Point;
typedef K::FT Coord_type;
typedef K::Weighted_point_2 Point;
struct Less {
bool operator()(const Point& p, const Point& q) const
{
struct Less
{
bool operator()(const Point& p, const Point& q) const {
return K::Less_xy_2()(p.point(), q.point());
}
};
typedef std::map<Point, Coord_type, Less> Point_value_map ;
typedef std::map<Point, K::Vector_2 , Less > Point_vector_map;
typedef std::map<Point, Coord_type, Less> Point_value_map ;
typedef std::map<Point, K::Vector_2 , Less> Point_vector_map;
int main()
{
@ -33,11 +39,11 @@ int main()
//parameters for spherical function:
Coord_type a(0.25), bx(1.3), by(-0.7), c(0.2);
for (int y=0 ; y<4 ; y++){
for (int x=0 ; x<4 ; x++){
for (int y=0; y<4; y++) {
for (int x=0; x<4; x++) {
Point p(x,y);
T.insert(p);
function_values.insert(std::make_pair(p,a + bx* x+ by*y + c*(x*x+y*y)));
function_values.insert(std::make_pair(p, a + bx*x + by*y + c*(x*x+y*y)));
}
}
@ -46,36 +52,34 @@ int main()
CGAL::Data_access<Point_value_map>(function_values),
Traits());
for(Point_vector_map::iterator it = function_gradients.begin(); it != function_gradients.end(); ++it)
{
for(Point_vector_map::iterator it = function_gradients.begin();
it != function_gradients.end(); ++it) {
std::cout << it->first << " " << it->second << std::endl;
}
//coordinate computation
Point p(1.6,1.4);
std::vector< std::pair< Point, Coord_type > > coords;
Coord_type norm = CGAL::regular_neighbor_coordinates_2(T, p, std::back_inserter
(coords)).second;
//coordinate computation
Point p(1.6, 1.4);
std::vector<std::pair<Point, Coord_type> > coords;
Coord_type norm = CGAL::regular_neighbor_coordinates_2(T, p, std::back_inserter(coords)).second;
//Sibson interpolant: version without sqrt:
std::pair<Coord_type, bool> res =
CGAL::sibson_c1_interpolation_square(
coords.begin(),
coords.end(),norm,p,
CGAL::Data_access<Point_value_map>(function_values),
CGAL::Data_access<Point_vector_map>(function_gradients),
Traits());
std::pair<Coord_type, bool> res = CGAL::sibson_c1_interpolation_square(coords.begin(),
coords.end(),
norm,
p,
CGAL::Data_access<Point_value_map>(function_values),
CGAL::Data_access<Point_vector_map>(function_gradients),
Traits());
if(res.second)
std::cout << "Tested interpolation on " << p
<< " interpolation: " << res.first << " exact: "
<< a + bx * p.x()+ by * p.y()+ c*(p.x()*p.x()+p.y()*p.y())
<< a + bx*p.x() + by*p.y()+ c*(p.x()*p.x()+p.y()*p.y())
<< std::endl;
else
std::cout << "C^1 Interpolation not successful." << std::endl
<< " not all function_gradients are provided." << std::endl
<< " You may resort to linear interpolation." << std::endl;
std::cout << "done" << std::endl;
return 0;
}

113
Interpolation/examples/Interpolation/sibson_interpolation_rn_vertex_with_info_2.cpp
View File

@ -1,95 +1,87 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Regular_triangulation_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/Interpolation_gradient_fitting_traits_2.h>
#include <CGAL/sibson_gradient_fitting.h>
#include <CGAL/interpolation_functions.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Regular_triangulation_2.h>
typedef K::FT Coord_type;
typedef K::Weighted_point_2 Point;
typedef K::Vector_2 Vector;
#include <iostream>
#include <iterator>
#include <utility>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
typedef K::FT Coord_type;
typedef K::Weighted_point_2 Point;
typedef K::Vector_2 Vector;
template <typename V, typename G>
struct Value_and_gradient {
Value_and_gradient()
: value(), gradient(CGAL::NULL_VECTOR)
{}
struct Value_and_gradient
{
Value_and_gradient() : value(), gradient(CGAL::NULL_VECTOR) {}
V value;
G gradient;
};
typedef CGAL::Triangulation_vertex_base_with_info_2<Value_and_gradient<Coord_type,Vector>, K, CGAL::Regular_triangulation_vertex_base_2<K> > Vb;
typedef CGAL::Triangulation_vertex_base_with_info_2<
Value_and_gradient<Coord_type, Vector>, K,
CGAL::Regular_triangulation_vertex_base_2<K> > Vb;
typedef CGAL::Regular_triangulation_face_base_2<K> Fb;
typedef CGAL::Triangulation_data_structure_2<Vb,Fb> Tds;
typedef CGAL::Regular_triangulation_2<K,Tds> Regular_triangulation;
typedef Regular_triangulation::Vertex_handle Vertex_handle;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
typedef CGAL::Regular_triangulation_2<K, Tds> Regular_triangulation;
typedef Regular_triangulation::Vertex_handle Vertex_handle;
template <typename V, typename T>
struct Function_value {
typedef V argument_type;
typedef std::pair<T, bool> result_type;
result_type operator()(const argument_type& a)const
{
struct Function_value
{
typedef V argument_type;
typedef std::pair<T, bool> result_type;
result_type operator()(const argument_type& a) const {
return result_type(a->info().value, true);
}
};
template <typename V, typename G>
struct Function_gradient
: public std::iterator<std::output_iterator_tag,void,void,void,void> {
typedef V argument_type;
typedef std::pair<G,bool> result_type;
result_type
operator()(const argument_type& a)const
{
return std::make_pair(a->info().gradient,a->info().gradient != CGAL::NULL_VECTOR) ;
: public std::iterator<std::output_iterator_tag, void, void, void, void>
{
typedef V argument_type;
typedef std::pair<G, bool> result_type;
result_type operator()(const argument_type& a) const {
return std::make_pair(a->info().gradient, a->info().gradient != CGAL::NULL_VECTOR);
}
const Function_gradient& operator=(const std::pair<V, G>& p) const
{
p.first->info().gradient = p.second;
return *this;
}
const Function_gradient& operator++(int) const
{
return *this;
}
const Function_gradient& operator*() const
{
const Function_gradient& operator=(const std::pair<V, G>& p) const {
p.first->info().gradient = p.second;
return *this;
}
const Function_gradient& operator++(int) const { return *this; }
const Function_gradient& operator*() const { return *this; }
};
int main()
{
Regular_triangulation rt;
Function_value<Vertex_handle,Coord_type> function_value;
Function_gradient<Vertex_handle,Vector> function_gradient;
Function_value<Vertex_handle, Coord_type> function_value;
Function_gradient<Vertex_handle, Vector> function_gradient;
//parameters for spherical function:
Coord_type a(0.25), bx(1.3), by(-0.7), c(0.2);
for (int y=0 ; y<4 ; y++){
for (int x=0 ; x<4 ; x++){
for (int y=0; y<4; y++) {
for (int x=0; x<4; x++) {
Point p(x,y);
Vertex_handle vh = rt.insert(p);
Coord_type value = a + bx* x+ by*y + c*(x*x+y*y);
Coord_type value = a + bx*x + by*y + c*(x*x+y*y);
vh->info().value = value;
}
}
@ -99,16 +91,16 @@ int main()
function_value,
CGAL::Identity<std::pair<Vertex_handle, Vector> >(),
Traits());
//coordinate computation
Point p(1.6,1.4);
std::vector< std::pair< Vertex_handle, Coord_type > > coords;
typedef CGAL::Identity<std::pair< Vertex_handle, Coord_type > > Identity;
std::vector<std::pair<Vertex_handle, Coord_type> > coords;
typedef CGAL::Identity<std::pair<Vertex_handle, Coord_type> > Identity;
Coord_type norm = CGAL::regular_neighbor_coordinates_2(rt,
p,
std::back_inserter(coords),
Identity()).second;
//Sibson interpolant: version without sqrt:
std::pair<Coord_type, bool> res = CGAL::sibson_c1_interpolation_square(coords.begin(),
coords.end(),
@ -117,7 +109,7 @@ int main()
function_value,
function_gradient,
Traits());
if(res.second)
std::cout << "Tested interpolation on " << p
<< " interpolation: " << res.first << " exact: "
@ -128,6 +120,5 @@ int main()
<< " not all function_gradients are provided." << std::endl
<< " You may resort to linear interpolation." << std::endl;
std::cout << "done" << std::endl;
return 0;
return EXIT_SUCCESS;
}

110
Interpolation/examples/Interpolation/sibson_interpolation_vertex_with_info_2.cpp
View File

@ -7,85 +7,75 @@
#include <CGAL/sibson_gradient_fitting.h>
#include <CGAL/interpolation_functions.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
#include <iostream>
#include <iterator>
#include <utility>
#include <vector>
typedef K::FT Coord_type;
typedef K::Point_2 Point;
typedef K::Vector_2 Vector;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
typedef K::FT Coord_type;
typedef K::Point_2 Point;
typedef K::Vector_2 Vector;
template <typename V, typename G>
struct Value_and_gradient {
Value_and_gradient()
: value(), gradient(CGAL::NULL_VECTOR)
{}
struct Value_and_gradient
{
Value_and_gradient() : value(), gradient(CGAL::NULL_VECTOR) {}
V value;
G gradient;
};
typedef CGAL::Triangulation_vertex_base_with_info_2<Value_and_gradient<Coord_type,Vector>, K> Vb;
typedef CGAL::Triangulation_data_structure_2<Vb> Tds;
typedef CGAL::Delaunay_triangulation_2<K,Tds> Delaunay_triangulation;
typedef Delaunay_triangulation::Vertex_handle Vertex_handle;
typedef CGAL::Triangulation_vertex_base_with_info_2<
Value_and_gradient<Coord_type, Vector>, K> Vb;
typedef CGAL::Triangulation_data_structure_2<Vb> Tds;
typedef CGAL::Delaunay_triangulation_2<K,Tds> Delaunay_triangulation;
typedef Delaunay_triangulation::Vertex_handle Vertex_handle;
template <typename V, typename T>
struct Function_value {
typedef V argument_type;
typedef std::pair<T, bool> result_type;
result_type operator()(const argument_type& a)const
{
struct Function_value
{
typedef V argument_type;
typedef std::pair<T, bool> result_type;
result_type operator()(const argument_type& a) const {
return result_type(a->info().value, true);
}
};
template <typename V, typename G>
struct Function_gradient
: public std::iterator<std::output_iterator_tag,void,void,void,void> {
typedef V argument_type;
typedef std::pair<G,bool> result_type;
result_type
operator()(const argument_type& a)const
{
return std::make_pair(a->info().gradient,a->info().gradient != CGAL::NULL_VECTOR) ;
: public std::iterator<std::output_iterator_tag, void, void, void, void>
{
typedef V argument_type;
typedef std::pair<G, bool> result_type;
result_type operator()(const argument_type& a) const {
return std::make_pair(a->info().gradient, a->info().gradient != CGAL::NULL_VECTOR);
}
const Function_gradient& operator=(const std::pair<V, G>& p) const
{
p.first->info().gradient = p.second;
return *this;
}
const Function_gradient& operator++(int) const
{
return *this;
}
const Function_gradient& operator*() const
{
const Function_gradient& operator=(const std::pair<V, G>& p) const {
p.first->info().gradient = p.second;
return *this;
}
const Function_gradient& operator++(int) const { return *this; }
const Function_gradient& operator*() const { return *this; }
};
int main()
{
Delaunay_triangulation dt;
Function_value<Vertex_handle,Coord_type> function_value;
Function_gradient<Vertex_handle,Vector> function_gradient;
Function_value<Vertex_handle, Coord_type> function_value;
Function_gradient<Vertex_handle, Vector> function_gradient;
//parameters for spherical function:
Coord_type a(0.25), bx(1.3), by(-0.7), c(0.2);
for (int y=0 ; y<4 ; y++){
for (int x=0 ; x<4 ; x++){
for (int y=0 ; y<4 ; y++) {
for (int x=0 ; x<4 ; x++) {
K::Point_2 p(x,y);
Vertex_handle vh = dt.insert(p);
Coord_type value = a + bx* x+ by*y + c*(x*x+y*y);
@ -100,15 +90,14 @@ int main()
Traits());
//coordinate computation
K::Point_2 p(1.6,1.4);
std::vector< std::pair< Vertex_handle, Coord_type > > coords;
typedef CGAL::Identity<std::pair< Vertex_handle, Coord_type > > Identity;
K::Point_2 p(1.6, 1.4);
std::vector<std::pair<Vertex_handle, Coord_type> > coords;
typedef CGAL::Identity<std::pair< Vertex_handle, Coord_type> > Identity;
Coord_type norm = CGAL::natural_neighbor_coordinates_2(dt,
p,
std::back_inserter(coords),
Identity()).second;
//Sibson interpolant: version without sqrt:
std::pair<Coord_type, bool> res = CGAL::sibson_c1_interpolation_square(coords.begin(),
coords.end(),
@ -117,17 +106,16 @@ int main()
function_value,
function_gradient,
Traits());
if(res.second)
std::cout << "Tested interpolation on " << p
<< " interpolation: " << res.first << " exact: "
<< a + bx * p.x()+ by * p.y()+ c*(p.x()*p.x()+p.y()*p.y())
<< a + bx*p.x() + by*p.y() + c*(p.x()*p.x() + p.y()*p.y())
<< std::endl;
else
std::cout << "C^1 Interpolation not successful." << std::endl
<< " not all function_gradients are provided." << std::endl
<< " not all function_gradients are provided." << std::endl
<< " You may resort to linear interpolation." << std::endl;
std::cout << "done" << std::endl;
return 0;
return EXIT_SUCCESS;
}

16
Interpolation/examples/Interpolation/surface_neighbor_coordinates_3.cpp
View File

@ -20,7 +20,7 @@ typedef std::vector< std::pair< Point_3, K::FT > > Point_coordinate_ve
int main()
{
int n=100;
int n = 100;
std::vector< Point_3> points;
points.reserve(n);
@ -28,32 +28,32 @@ int main()
CGAL::Random_points_on_sphere_3<Point_3> g(1);
CGAL::cpp11::copy_n(g, n, std::back_inserter(points));
Point_3 p(1, 0,0);
Point_3 p(1, 0, 0);
Vector_3 normal(p - CGAL::ORIGIN);
std::cout << "Compute surface neighbor coordinates for " << p << std::endl;
Point_coordinate_vector coords;
CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>,
K::FT, bool> result =
CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>, K::FT, bool> result =
CGAL::surface_neighbor_coordinates_3(points.begin(), points.end(),
p, normal,
std::back_inserter(coords),
K());
if(!result.third){
if(!result.third)
{
//Undersampling:
std::cout << "The coordinate computation was not successful." << std::endl;
return 0;
}
K::FT norm = result.second;
std::cout << "Testing the barycentric property " << std::endl;
Point_3 b(0, 0, 0);
for(std::vector< std::pair< Point_3, Coord_type > >::const_iterator
it = coords.begin(); it!=coords.end(); ++it)
b = b + (it->second/norm)* (it->first - CGAL::ORIGIN);
b = b + (it->second/norm) * (it->first - CGAL::ORIGIN);
std::cout << " weighted barycenter: " << b <<std::endl;
std::cout << " squared distance: " << CGAL::squared_distance(p,b) << std::endl;
std::cout << "done" << std::endl;
return 0;
return EXIT_SUCCESS;
}

235
Interpolation/include/CGAL/Interpolation/internal/helpers.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2003,2017 INRIA Sophia-Antipolis (France).
// Copyright (c) 2003, 2017 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
@ -23,134 +23,137 @@
#include <CGAL/license/Interpolation.h>
namespace CGAL {
namespace Interpolation {
namespace internal {
template < class InterpolationTraits>
struct V2P
{
typedef typename InterpolationTraits::Point_d Point;
typedef typename InterpolationTraits::Weighted_point_d Weighted_point;
V2P(const InterpolationTraits& traits)
: traits(traits)
{}
template <typename VH>
const Point& operator()(const VH& vh) const
{
return traits.construct_point_d_object()(vh->point());
}
const Point& operator()(const Point& p) const
{
return p;
}
Point operator()(const Weighted_point& wp) const
{
return traits.construct_point_d_object()(wp);
}
template < class InterpolationTraits >
struct V2P
{
typedef typename InterpolationTraits::Point_d Point;
typedef typename InterpolationTraits::Weighted_point_d Weighted_point;
private:
InterpolationTraits traits;
};
V2P(const InterpolationTraits& traits)
: traits(traits)
{}
template < typename Dt, typename T2>
struct Vertex2Point {
typedef typename Dt::Vertex_handle Vertex_handle;
typedef typename Dt::Point Point;
typedef std::pair<Vertex_handle, T2> argument_type;
typedef std::pair<Point, T2> result_type;
result_type operator()(const argument_type& vp) const
{
return std::make_pair(vp.first->point(), vp.second);
}
};
template < typename Dt, typename T2>
struct Vertex2WPoint {
typedef typename Dt::Vertex_handle Vertex_handle;
typedef typename Dt::Weighted_point Point;
typedef std::pair<Vertex_handle, T2> argument_type;
typedef std::pair<Point, T2> result_type;
result_type operator()(const argument_type& vp) const
{
return std::make_pair(vp.first->point(), vp.second);
}
};
template <typename Dt, typename Map>
struct Vertex2Vertex {
typedef typename Dt::Vertex_handle Vertex_handle;
typedef typename Dt::Geom_traits::FT FT;
typedef std::pair<Vertex_handle,FT> argument_type;
typedef std::pair<Vertex_handle,FT> result_type;
const Map& map;
const Dt& dt;
Vertex2Vertex(const Map& map, const Dt& dt)
: map(map), dt(dt)
{}
result_type operator()(const argument_type& vp) const
{
typename Map::const_iterator it = map.find(vp.first);
CGAL_assertion(it != map.end());
CGAL_assertion(dt.tds().is_vertex(it->second));
return std::make_pair(it->second, vp.second);
}
};
// the struct "Project_vertex_output_iterator"
// is used in the (next two) functions
// as well as in regular_neighbor_coordinates_2 and
// in surface_neighbor_coordinates_3
//
//projection of iterator over std::pair <Vertex_handle, T>
//to iterator over std::pair< Point, T>
template < class OutputIterator, class Fct = void>
struct Project_vertex_output_iterator
template <typename VH>
const Point& operator()(const VH& vh) const
{
// this class wraps the OutputIterator with value type
// std::pair<Vertex_handle,T>
// into an output iterator with value type std::pair<Point, T>
// Conditions: OutputIterator has value type std::pair<Vertex_handle, T>
// and Vertex_handle has a function ->point()
// with return type const Point&
return traits.construct_point_d_object()(vh->point());
}
OutputIterator _base;
Fct fct;
const Point& operator()(const Point& p) const
{
return p;
}
//creation:
Project_vertex_output_iterator(OutputIterator o, Fct fct)
: _base(o), fct(fct)
{}
Point operator()(const Weighted_point& wp) const
{
return traits.construct_point_d_object()(wp);
}
OutputIterator base() {return _base;}
private:
InterpolationTraits traits;
};
Project_vertex_output_iterator& operator++(){_base++; return *this;}
Project_vertex_output_iterator& operator++(int){_base++; return *this;}
Project_vertex_output_iterator& operator*(){return *this;}
template<class Vertex_pair>
Project_vertex_output_iterator&
operator=(const Vertex_pair& vp){
*_base = fct(vp);
return *this;
}
};
template < typename Dt, typename T2 >
struct Vertex2Point
{
typedef typename Dt::Vertex_handle Vertex_handle;
typedef typename Dt::Point Point;
typedef std::pair<Vertex_handle, T2> argument_type;
typedef std::pair<Point, T2> result_type;
result_type operator()(const argument_type& vp) const
{
return std::make_pair(vp.first->point(), vp.second);
}
};
template < typename Dt, typename T2 >
struct Vertex2WPoint
{
typedef typename Dt::Vertex_handle Vertex_handle;
typedef typename Dt::Weighted_point Point;
typedef std::pair<Vertex_handle, T2> argument_type;
typedef std::pair<Point, T2> result_type;
result_type operator()(const argument_type& vp) const
{
return std::make_pair(vp.first->point(), vp.second);
}
};
template < typename Dt, typename Map >
struct Vertex2Vertex
{
typedef typename Dt::Vertex_handle Vertex_handle;
typedef typename Dt::Geom_traits::FT FT;
typedef std::pair<Vertex_handle, FT> argument_type;
typedef std::pair<Vertex_handle, FT> result_type;
const Map& map;
const Dt& dt;
Vertex2Vertex(const Map& map, const Dt& dt)
: map(map), dt(dt)
{}
result_type operator()(const argument_type& vp) const
{
typename Map::const_iterator it = map.find(vp.first);
CGAL_assertion(it != map.end());
CGAL_assertion(dt.tds().is_vertex(it->second));
return std::make_pair(it->second, vp.second);
}
};
// the struct "Project_vertex_output_iterator"
// is used in the (next two) functions
// as well as in regular_neighbor_coordinates_2 and
// in surface_neighbor_coordinates_3
//
//projection of iterator over std::pair <Vertex_handle, T>
//to iterator over std::pair< Point, T>
template < class OutputIterator, class Fct = void >
struct Project_vertex_output_iterator
{
// this class wraps the OutputIterator with value type
// std::pair<Vertex_handle,T>
// into an output iterator with value type std::pair<Point, T>
// Conditions: OutputIterator has value type std::pair<Vertex_handle, T>
// and Vertex_handle has a function ->point()
// with return type const Point&
OutputIterator _base;
Fct fct;
//creation:
Project_vertex_output_iterator(OutputIterator o, Fct fct)
: _base(o), fct(fct)
{}
OutputIterator base() {return _base;}
Project_vertex_output_iterator& operator++(){_base++; return *this;}
Project_vertex_output_iterator& operator++(int){_base++; return *this;}
Project_vertex_output_iterator& operator*(){return *this;}
template<class Vertex_pair>
Project_vertex_output_iterator& operator=(const Vertex_pair& vp)
{
*_base = fct(vp);
return *this;
}
};
} // namespace internal
} // namespace Interpolation
} // namespace CGAL

24
Interpolation/include/CGAL/Voronoi_intersection_2_traits_3.h
View File

@ -163,16 +163,16 @@ public:
Comparison_result operator()(const Point& p, const Point& q) const
{
if(normal.x()!=Coord_type(0))
if(normal.x() != Coord_type(0))
return (Comparison_result) CGAL_NTS
sign(Vector(normal.y(),-normal.x(),Coord_type(0))*(p-q));
if(normal.y()!= Coord_type(0))
sign(Vector(normal.y(), -normal.x(), Coord_type(0))*(p-q));
if(normal.y() != Coord_type(0))
return (Comparison_result) CGAL_NTS
sign(Vector(-normal.y(),normal.x(),Coord_type(0))*(p-q));
sign(Vector(-normal.y(), normal.x(), Coord_type(0))*(p-q));
CGAL_assertion(normal.z()!= Coord_type(0));
CGAL_assertion(normal.z() != Coord_type(0));
return (Comparison_result) CGAL_NTS
sign(Vector(-normal.z(),Coord_type(0),normal.x())*(p-q));
sign(Vector(-normal.z(), Coord_type(0), normal.x())*(p-q));
}
private:
@ -195,16 +195,16 @@ public:
Comparison_result operator()(const Point& p, const Point& q) const
{
if(normal.x()!=Coord_type(0))
if(normal.x() != Coord_type(0))
return (Comparison_result) CGAL_NTS
sign(Vector(normal.z(),Coord_type(0),-normal.x())*(p-q));
if(normal.y()!= Coord_type(0))
sign(Vector(normal.z(), Coord_type(0), -normal.x())*(p-q));
if(normal.y() != Coord_type(0))
return (Comparison_result) CGAL_NTS
sign(Vector(Coord_type(0),normal.z(),-normal.y())*(p-q));
sign(Vector(Coord_type(0), normal.z(), -normal.y())*(p-q));
CGAL_assertion(normal.z()!= Coord_type(0));
CGAL_assertion(normal.z() != Coord_type(0));
return (Comparison_result) CGAL_NTS
sign(Vector(Coord_type(0),-normal.z(),normal.y())*(p-q));
sign(Vector(Coord_type(0), -normal.z(), normal.y())*(p-q));
}
private:

33
Interpolation/include/CGAL/natural_neighbor_coordinates_2.h
View File

@ -22,19 +22,21 @@
#define CGAL_NATURAL_NEIGHBOR_COORDINATES_2_H
#include <CGAL/license/Interpolation.h>
#include <CGAL/Interpolation/internal/helpers.h>
#include <CGAL/Iterator_project.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/number_utils_classes.h>
#include <CGAL/utility.h>
#include <iterator>
#include <list>
#include <utility>
#include <vector>
namespace CGAL {
// The following natural_neighbor_coordinate_2 functions fix the
// traits class to be Dt::Geom_traits. The following signatures could
// be used if one wants to pass a traits class as argument:
@ -60,8 +62,8 @@ template < class Dt, class OutputIterator >
Triple< OutputIterator, typename Dt::Geom_traits::FT, bool >
natural_neighbors_2(const Dt& dt,
const typename Dt::Geom_traits::Point_2& p,
OutputIterator out, typename Dt::Face_handle start
= typename Dt::Face_handle())
OutputIterator out,
typename Dt::Face_handle start = typename Dt::Face_handle())
{
typedef typename Dt::Geom_traits Traits;
typedef typename Traits::FT Coord_type;
@ -96,13 +98,13 @@ natural_neighbors_2(const Dt& dt,
Equal_x_2 equal_x_2;
if(!equal_x_2(p1,p2))
{
coef1 = (p.x() - p2.x()) / (p1.x() - p2.x());
coef1 = (p.x() - p2.x()) / (p1.x() - p2.x());
coef2 = 1 - coef1;
*out++ = std::make_pair(v1,coef1);
*out++ = std::make_pair(v2,coef2);
} else {
coef1 = (p.y() - p2.y()) / (p1.y() - p2.y());
coef2 = 1-coef1;
coef2 = 1 - coef1;
*out++ = std::make_pair(v1,coef1);
*out++ = std::make_pair(v2,coef2);
}
@ -122,7 +124,7 @@ natural_neighbors_2(const Dt& dt,
return natural_neighbors_2(dt, p, out, hole.begin(), hole.end());
}
//function call if the conflict zone is known:
// OutputIterator has value type
// std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT>
@ -148,6 +150,7 @@ natural_neighbors_2(const Dt& dt,
Coord_type area_sum(0);
EdgeIterator hit = hole_end;
--hit;
//in the beginning: prev is the "last" vertex of the hole:
// later: prev is the last vertex processed (previously)
Vertex_handle prev = hit->first->vertex(dt.cw(hit->second));
@ -190,7 +193,7 @@ natural_neighbors_2(const Dt& dt,
return make_triple(out, area_sum, true);
}
template < class Dt, class OutputIterator, class Fct >
Triple< OutputIterator, typename Dt::Geom_traits::FT, bool >
natural_neighbor_coordinates_2(const Dt& dt,
@ -205,8 +208,7 @@ natural_neighbor_coordinates_2(const Dt& dt,
Interpolation::internal::Project_vertex_output_iterator<OutputIterator, Fct> op(out, fct);
Triple<Interpolation::internal::Project_vertex_output_iterator<OutputIterator,Fct>,
typename Dt::Geom_traits::FT, bool > result =
natural_neighbors_2(dt, p, op, start);
typename Dt::Geom_traits::FT, bool > result = natural_neighbors_2(dt, p, op, start);
return make_triple(result.first.base(), result.second, result.third);
}
@ -255,6 +257,7 @@ natural_neighbor_coordinates_2(const Dt& dt,
return make_triple(result.first.base(), result.second, result.third);
}
//OutputIterator has value type
// std::pair< Dt::Geom_traits::Point_2, Dt::Geom_traits::FT>
//function call if the conflict zone is known:
@ -269,7 +272,8 @@ natural_neighbor_coordinates_2(const Dt& dt,
return natural_neighbor_coordinates_2(dt,p, out, hole_begin, hole_end,
Interpolation::internal::Vertex2Point<Dt, typename Dt::Geom_traits::FT>());
}
/**********************************************************/
//compute the coordinates for a vertex of the triangulation
// with respect to the other points in the triangulation
@ -295,11 +299,11 @@ natural_neighbor_coordinates_2(const Dt& dt,
typedef std::map<Vertex_handle /*t2*/, Vertex_handle/*dt*/> V2V;
V2V correspondence_map;
do{
CGAL_assertion(!dt.is_infinite(vc));
Vertex_handle vh2 = t2.insert(vc->point());
correspondence_map[vh2] = vc;
correspondence_map[vh2] = vc;
}
while(++vc!=done);
@ -308,7 +312,7 @@ natural_neighbor_coordinates_2(const Dt& dt,
return natural_neighbor_coordinates_2(t2, vh->point(), out, cfct);
}
template <class Dt, class OutputIterator>
Triple< OutputIterator, typename Dt::Geom_traits::FT, bool >
natural_neighbor_coordinates_2(const Dt& dt,
@ -340,7 +344,6 @@ public:
}
};
} //namespace CGAL
#endif // CGAL_NATURAL_NEIGHBOR_COORDINATES_2_H

18
Interpolation/include/CGAL/natural_neighbor_coordinates_3.h
View File

@ -72,8 +72,8 @@ construct_circumcenter(const typename DT::Facet& f,
template <class Dt, class OutputIterator>
Triple< OutputIterator, // iterator with value type std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT>
typename Dt::Geom_traits::FT, // Should provide 0 and 1
bool >
typename Dt::Geom_traits::FT, // Should provide 0 and 1
bool >
laplace_natural_neighbor_coordinates_3(const Dt& dt,
const typename Dt::Geom_traits::Point_3& Q,
OutputIterator nn_out,
@ -108,23 +108,25 @@ laplace_natural_neighbor_coordinates_3(const Dt& dt,
std::set<Cell_handle> cells;
// To replace the forbidden access to the "in conflict" flag :
// std::find operations on this set
std::vector<Facet> bound_facets; bound_facets.reserve(32);
typename std::vector<Facet>::iterator bound_it;
std::vector<Facet> bound_facets;
bound_facets.reserve(32);
// Find the cells in conflict with Q
dt.find_conflicts(Q, c,
std::back_inserter(bound_facets),
std::inserter(cells,cells.begin()));
std::inserter(cells, cells.begin()));
std::map<Vertex_handle,Coord_type> coordinate;
typename std::map<Vertex_handle,Coord_type>::iterator coor_it;
typename std::vector<Facet>::iterator bound_it;
for (bound_it = bound_facets.begin(); bound_it != bound_facets.end(); ++bound_it)
{
//for each facet on the boundary
Facet f1 = *bound_it;
Cell_handle cc1 = f1.first;
if (dt.is_infinite(cc1))
return make_triple(nn_out,norm_coeff=Coord_type(1), false);//point outside the convex-hull
return make_triple(nn_out, norm_coeff=Coord_type(1), false);//point outside the convex-hull
CGAL_triangulation_assertion_code(Cell_handle cc2 = cc1->neighbor(f1.second);)
CGAL_triangulation_assertion(std::find(cells.begin(),cells.end(),cc1) != cells.end());//TODO : Delete
@ -172,8 +174,8 @@ laplace_natural_neighbor_coordinates_3(const Dt& dt,
template <class Dt, class OutputIterator>
Triple< OutputIterator, // iterator with value type std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT>
typename Dt::Geom_traits::FT, // Should provide 0 and 1
bool >
typename Dt::Geom_traits::FT, // Should provide 0 and 1
bool >
sibson_natural_neighbor_coordinates_3(const Dt& dt,
const typename Dt::Geom_traits::Point_3& Q,
OutputIterator nn_out,

101
Interpolation/include/CGAL/regular_neighbor_coordinates_2.h
View File

@ -120,8 +120,7 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
rt.get_boundary_of_conflicts_and_hidden_vertices(p,
std::back_inserter(hole),
std::back_inserter
(hidden_vertices),
std::back_inserter(hidden_vertices),
fh);
return regular_neighbor_coordinates_vertex_2(rt, p, out, vor_vertices,
hole.begin(),hole.end(),
@ -171,8 +170,8 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
typedef typename Traits::FT Coord_type;
typedef typename Rt::Bare_point Bare_point;
typedef typename Rt::Vertex_handle Vertex_handle;
typedef typename Rt::Face_circulator Face_circulator;
typedef typename Rt::Vertex_handle Vertex_handle;
typedef typename Rt::Face_circulator Face_circulator;
//no hole because only (exactly!) one vertex is hidden:
if(hole_begin==hole_end){
@ -193,54 +192,55 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
Vertex_handle prev = hit->first->vertex(rt.cw(hit->second));
hit = hole_begin;
while(hit != hole_end)
{
Coord_type area(0);
Vertex_handle current = hit->first->vertex(rt.cw(hit->second));
//a first Voronoi vertex of the cell of p:
vor[0] = rt.geom_traits().construct_weighted_circumcenter_2_object()
(current->point(),
hit->first->vertex(rt.ccw(hit->second))->point(), p);
*vor_vertices++= vor[0];
//triangulation of the Voronoi subcell:
//a second vertex as base
Face_circulator fc = rt.incident_faces(current, hit->first);
++fc;
vor[1] = rt.dual(fc);
// iteration over all other "old" Voronoi vertices
while(!fc->has_vertex(prev))
{
Coord_type area(0);
Vertex_handle current = hit->first->vertex(rt.cw(hit->second));
//a first Voronoi vertex of the cell of p:
vor[0] = rt.geom_traits().construct_weighted_circumcenter_2_object()
(current->point(),
hit->first->vertex(rt.ccw(hit->second))->point(), p);
*vor_vertices++= vor[0];
//triangulation of the Voronoi subcell:
//a second vertex as base
Face_circulator fc = rt.incident_faces(current, hit->first);
++fc;
vor[1] = rt.dual(fc);
// iteration over all other "old" Voronoi vertices
while(!fc->has_vertex(prev))
{
++fc;
vor[2] = rt.dual(fc);
area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits());
vor[1] = vor[2];
}
//the second Voronoi vertex of the cell of p:
vor[2] =
rt.geom_traits().construct_weighted_circumcenter_2_object()
(prev->point(),current->point(),p);
*vor_vertices++= vor[2];
vor[2] = rt.dual(fc);
area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits());
*out++= std::make_pair(current,area);
area_sum += area;
//update prev and hit:
prev= current;
++hit;
vor[1] = vor[2];
}
//the second Voronoi vertex of the cell of p:
vor[2] =
rt.geom_traits().construct_weighted_circumcenter_2_object()
(prev->point(),current->point(),p);
*vor_vertices++= vor[2];
area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits());
*out++= std::make_pair(current,area);
area_sum += area;
//update prev and hit:
prev= current;
++hit;
}
//get coordinate for hidden vertices
// <=> the area of their Voronoi cell.
//decomposition of the cell into triangles
// vor1: dual of first triangle
// vor2, vor 3: duals of two consecutive triangles
Face_circulator fc, fc_begin;
for(; hidden_vertices_begin != hidden_vertices_end;
++hidden_vertices_begin){
for(; hidden_vertices_begin != hidden_vertices_end; ++hidden_vertices_begin)
{
Coord_type area(0);
fc_begin = rt.incident_faces(*hidden_vertices_begin);
vor[0] = rt.dual(fc_begin);
@ -248,7 +248,8 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
++fc;
vor[1] = rt.dual(fc);
++fc;
while(fc != fc_begin){
while(fc != fc_begin)
{
vor[2] = rt.dual(fc);
area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits());
@ -263,7 +264,7 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
return make_triple(out, area_sum, true);
}
////////////////////////////////////////////////////////////
//the cast from vertex to point:
// the following functions return an Output_iterator over
@ -289,7 +290,7 @@ regular_neighbor_coordinates_2(const Rt& rt,
Fct fct,
typename boost::disable_if_c<
boost::is_convertible<Fct,
typename Rt::Face_handle>::value
typename Rt::Face_handle>::value
>::type* = 0)
{
return regular_neighbor_coordinates_2(rt, p, out, fct, typename Rt::Face_handle());
@ -346,7 +347,7 @@ regular_neighbor_coordinates_2(const Rt& rt,
return make_triple(result.first.base(), result.second, result.third);
}
template <class Rt, class OutputIterator, class OutputIteratorVorVertices>
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
regular_neighbor_coordinates_2(const Rt& rt,
@ -361,7 +362,7 @@ regular_neighbor_coordinates_2(const Rt& rt,
vor_vertices, start);
}
//OutputIterator has value type
// std::pair< Rt::Geom_traits::Weighted_point_2, Rt::Geom_traits::FT>
template <class Rt, class OutputIterator, class EdgeIterator,
@ -410,7 +411,7 @@ regular_neighbor_coordinates_2(const Rt& rt,
return make_triple(result.first.base(), result.second, result.third);
}
template <class Rt, class OutputIterator, class EdgeIterator,
class VertexIterator , class OutputIteratorVorVertices >
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
@ -454,11 +455,11 @@ regular_neighbor_coordinates_2(const Rt& rt,
typedef std::map<Vertex_handle /*t2*/, Vertex_handle/*dt*/> V2V;
V2V correspondence_map;
do{
CGAL_assertion(!rt.is_infinite(vc));
Vertex_handle vh2 = t2.insert(vc->point());
correspondence_map[vh2] = vc;
correspondence_map[vh2] = vc;
}
while(++vc!=done);
@ -468,7 +469,7 @@ regular_neighbor_coordinates_2(const Rt& rt,
return regular_neighbor_coordinates_2(t2, vh->point(), out, cfct);
}
template <class Rt, class OutputIterator>
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
regular_neighbor_coordinates_2(const Rt& rt,
@ -481,7 +482,7 @@ regular_neighbor_coordinates_2(const Rt& rt,
Interpolation::internal::Vertex2Point<Rt, typename Rt::Geom_traits::FT>());
}
//class providing a function object:
//OutputIterator has value type
// std::pair< Rt::Geom_traits::Weighted_point_2, Rt::Geom_traits::FT>

22
Interpolation/include/CGAL/sibson_gradient_fitting.h
View File

@ -82,7 +82,7 @@ sibson_gradient_fitting(ForwardIterator first,
return Hn.inverse().transform(pn);
}
template < class Triangul, class ForwardIterator, class Functor, class Traits, class VH>
typename Traits::Vector_d
sibson_gradient_fitting(const Triangul& tr,
@ -98,9 +98,9 @@ sibson_gradient_fitting(const Triangul& tr,
{
const typename Traits::Point_d& bare_p = traits.construct_point_d_object()(vh->point());
typename Functor::result_type fn = function_value(bare_p);
CGAL_assertion(fn.second);
return sibson_gradient_fitting(first,
beyond,
norm,
@ -109,7 +109,7 @@ sibson_gradient_fitting(const Triangul& tr,
function_value,
traits);
}
template < class Triangul, class ForwardIterator, class Functor, class Traits, class VH>
typename Traits::Vector_d
sibson_gradient_fitting(const Triangul& tr,
@ -126,7 +126,7 @@ sibson_gradient_fitting(const Triangul& tr,
const typename Traits::Point_d& bare_p = traits.construct_point_d_object()(vh->point());
typename Functor::result_type fn = function_value(vh->point());
CGAL_assertion(fn.second);
return sibson_gradient_fitting(first,
beyond,
norm,
@ -135,7 +135,7 @@ sibson_gradient_fitting(const Triangul& tr,
function_value,
traits);
}
template < class Triangul, class ForwardIterator, class Functor, class Traits, class VH>
typename Traits::Vector_d
sibson_gradient_fitting(const Triangul& tr,
@ -152,7 +152,7 @@ sibson_gradient_fitting(const Triangul& tr,
const typename Traits::Point_d& bare_p = traits.construct_point_d_object()(vh->point());
typename Functor::result_type fn = function_value(vh);
CGAL_assertion(fn.second);
return sibson_gradient_fitting(first,
beyond,
norm,
@ -162,7 +162,7 @@ sibson_gradient_fitting(const Triangul& tr,
traits);
}
template < class Triangul, class OutputIterator, class Functor,
class CoordFunctor, class OIF, class Traits>
OutputIterator
@ -197,14 +197,14 @@ sibson_gradient_fitting_internal(const Triangul& tr,
function_value,
traits,
typename Functor::argument_type())));
coords.clear();
}
}
return out;
}
//the following functions allow to fit the gradients for all points in
// a triangulation except the convex hull points.
// -> _nn2: natural_neighbor_coordinates_2
@ -242,7 +242,7 @@ template < class Dt, class OutputIterator, class Functor, class Traits>
OutputIterator
sibson_gradient_fitting_nn_2(const Dt& dt,
OutputIterator out,
Functor function_value,
Functor function_value,
const Traits& traits)
{
return sibson_gradient_fitting_nn_2(dt, out, function_value,