1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Triangulation_data_structure_2.h>
#include <CGAL/Delaunay_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 <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::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) {}
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;
template <typename V, typename T>
struct Value_function
{
typedef V argument_type;
typedef std::pair<T, bool> result_type;
// read operation
result_type operator()(const argument_type& a) const {
return result_type(a->info().value, true);
}
};
template <typename V, typename G>
struct Gradient_function
{
typedef void value_type;
typedef void Distance;
typedef void Pointer;
typedef void Reference;
typedef std::output_iterator_tag category;
typedef V argument_type;
typedef std::pair<G, bool> result_type;
// read operation
result_type operator()(const argument_type& a) const {
return std::make_pair(a->info().gradient, a->info().gradient != CGAL::NULL_VECTOR);
}
// write operation
const Gradient_function& operator=(const std::pair<V, G>& p) const {
p.first->info().gradient = p.second;
return *this;
}
const Gradient_function& operator++(int) const { return *this; }
const Gradient_function& operator*() const { return *this; }
};
int main()
{
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;
Gradient_function<Vertex_handle, Vector> gradient_function;
// 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) {
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);
vh->info().value = value; // store the value directly in the vertex
}
}
CGAL::sibson_gradient_fitting_nn_2(dt,
gradient_function,
CGAL::Identity<std::pair<Vertex_handle, Vector> >(),
value_function,
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;
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(),
norm,
p,
value_function,
gradient_function,
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())
<< 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;
return EXIT_SUCCESS;
}
|
