1 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
2
3 #include <CGAL/natural_neighbor_coordinates_2.h>
4 #include <CGAL/Interpolation_gradient_fitting_traits_2.h>
5 #include <CGAL/sibson_gradient_fitting.h>
6 #include <CGAL/interpolation_functions.h>
7
8 #include <CGAL/Triangulation_vertex_base_with_info_2.h>
9 #include <CGAL/Regular_triangulation_2.h>
10
11 #include <boost/iterator/function_output_iterator.hpp>
12
13 #include <iostream>
14 #include <iterator>
15 #include <utility>
16 #include <vector>
17
18 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
19 typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
20
21 typedef K::FT Coord_type;
22 typedef K::Point_2 Bare_point;
23 typedef K::Weighted_point_2 Weighted_point;
24 typedef K::Vector_2 Vector;
25
26 template <typename V, typename G>
27 struct Value_and_gradient
28 {
Value_and_gradientValue_and_gradient29 Value_and_gradient() : value(), gradient(CGAL::NULL_VECTOR) {}
30
31 V value;
32 G gradient;
33 };
34
35 typedef CGAL::Triangulation_vertex_base_with_info_2<
36 Value_and_gradient<Coord_type, Vector>, K,
37 CGAL::Regular_triangulation_vertex_base_2<K> > Vb;
38 typedef CGAL::Regular_triangulation_face_base_2<K> Fb;
39 typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
40 typedef CGAL::Regular_triangulation_2<K, Tds> Regular_triangulation;
41 typedef Regular_triangulation::Vertex_handle Vertex_handle;
42
main()43 int main()
44 {
45 Regular_triangulation rt;
46
47 auto value_function = [](const Vertex_handle& a) -> std::pair<Coord_type, bool>
48 {
49 return std::make_pair(a->info().value, true);
50 };
51
52 auto gradient_function = [](const Vertex_handle& a) -> std::pair<Vector, bool>
53 {
54 return std::make_pair(a->info().gradient, a->info().gradient != CGAL::NULL_VECTOR);
55 };
56
57 auto gradient_output_iterator
58 = boost::make_function_output_iterator
59 ([](const std::pair<Vertex_handle, Vector>& p)
60 {
61 p.first->info().gradient = p.second;
62 });
63
64 // parameters for spherical function:
65 Coord_type a(0.25), bx(1.3), by(-0.7), c(0.2);
66 for (int y=0; y<4; y++) {
67 for (int x=0; x<4; x++) {
68 Weighted_point p(Bare_point(x,y), (x-y)/3. /*weight*/);
69 Vertex_handle vh = rt.insert(p);
70 Coord_type value = a + bx*x + by*y + c*(x*x+y*y);
71 vh->info().value = value;
72 }
73 }
74
75 CGAL::sibson_gradient_fitting_rn_2(rt,
76 gradient_output_iterator,
77 CGAL::Identity<std::pair<Vertex_handle, Vector> >(),
78 value_function,
79 Traits());
80
81 // coordinate computation
82 Weighted_point p(Bare_point(1.6, 1.4), -0.3 /*weight*/);
83 std::vector<std::pair<Vertex_handle, Coord_type> > coords;
84 typedef CGAL::Identity<std::pair<Vertex_handle, Coord_type> > Identity;
85 Coord_type norm = CGAL::regular_neighbor_coordinates_2(rt,
86 p,
87 std::back_inserter(coords),
88 Identity()).second;
89
90 // Sibson interpolant: version without sqrt:
91 std::pair<Coord_type, bool> res = CGAL::sibson_c1_interpolation_square(coords.begin(),
92 coords.end(),
93 norm,
94 p,
95 value_function,
96 gradient_function,
97 Traits());
98
99 if(res.second)
100 std::cout << "Tested interpolation on " << p
101 << " interpolation: " << res.first << " exact: "
102 << a + bx * p.x()+ by * p.y()+ c*(p.x()*p.x()+p.y()*p.y())
103 << std::endl;
104 else
105 std::cout << "C^1 Interpolation not successful." << std::endl
106 << " not all gradients are provided." << std::endl
107 << " You may resort to linear interpolation." << std::endl;
108
109 return EXIT_SUCCESS;
110 }
111