1 #include <CGAL/Delaunay_mesher_2.h>
2 #include <CGAL/Interpolation_traits_2.h>
3 #include <CGAL/Projection_traits_xy_3.h>
4 #include <CGAL/interpolation_functions.h>
5 #include <CGAL/Delaunay_mesh_face_base_2.h>
6 #include <CGAL/Delaunay_mesh_size_criteria_2.h>
7 #include <CGAL/Constrained_Delaunay_triangulation_2.h>
8 #include <CGAL/Barycentric_coordinates_2/Mean_value_2.h>
9 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
10 #include <CGAL/Barycentric_coordinates_2/Generalized_barycentric_coordinates_2.h>
11
12 // Some convenient typedefs.
13
14 // General.
15 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
16 typedef CGAL::Projection_traits_xy_3<Kernel> Projection;
17
18 typedef Projection::FT Scalar;
19 typedef Projection::Point_2 Point;
20
21 typedef std::vector<Scalar> Scalar_vector;
22 typedef std::vector<Point> Point_vector;
23
24 // Coordinates related.
25 typedef CGAL::Barycentric_coordinates::Mean_value_2<Projection> Mean_value;
26 typedef CGAL::Barycentric_coordinates::Generalized_barycentric_coordinates_2<Mean_value, Projection> Mean_value_coordinates;
27
28 // Triangulation related.
29 typedef CGAL::Delaunay_mesh_face_base_2<Projection> Face_base;
30 typedef CGAL::Triangulation_vertex_base_2<Projection> Vertex_base;
31 typedef CGAL::Triangulation_data_structure_2<Vertex_base, Face_base> TDS;
32 typedef CGAL::Constrained_Delaunay_triangulation_2<Projection, TDS> CDT;
33 typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Criteria;
34 typedef CGAL::Delaunay_mesher_2<CDT, Criteria> Mesher;
35 typedef CDT::Finite_vertices_iterator Vertex_iterator;
36 typedef CDT::Vertex_handle Vertex_handle;
37
38 // Interpolation related.
39 typedef CGAL::Interpolation_traits_2<Projection> Interpolation_traits;
40 typedef CGAL::Data_access< std::map<Point, Scalar, Projection::Less_xy_2 > > Value_access;
41
42 // STD.
43 using std::cout; using std::endl; using std::string;
44
main()45 int main()
46 {
47 // Construct a polygon bounding a piece of three-dimensional terrain.
48 // Note that z-coordinate of each vertex represents the height function.
49 // Projection in 2D is done automatically by the Projection traits class.
50 const int number_of_vertices = 50;
51 Point_vector vertices(number_of_vertices);
52
53 vertices[0] = Point(0.03, 0.05, 0.000); vertices[1] = Point(0.07, 0.04, 10.00); vertices[2] = Point(0.10, 0.04, 20.00);
54 vertices[3] = Point(0.14, 0.04, 30.00); vertices[4] = Point(0.17, 0.07, 40.00); vertices[5] = Point(0.19, 0.09, 50.00);
55 vertices[6] = Point(0.22, 0.11, 60.00); vertices[7] = Point(0.25, 0.11, 70.00); vertices[8] = Point(0.27, 0.10, 80.00);
56 vertices[9] = Point(0.30, 0.07, 90.00); vertices[10] = Point(0.31, 0.04, 100.0); vertices[11] = Point(0.34, 0.03, 110.0);
57 vertices[12] = Point(0.37, 0.02, 120.0); vertices[13] = Point(0.40, 0.03, 130.0); vertices[14] = Point(0.42, 0.04, 140.0);
58 vertices[15] = Point(0.44, 0.07, 150.0); vertices[16] = Point(0.45, 0.10, 160.0); vertices[17] = Point(0.46, 0.13, 170.0);
59 vertices[18] = Point(0.46, 0.19, 180.0); vertices[19] = Point(0.47, 0.26, 190.0); vertices[20] = Point(0.47, 0.31, 200.0);
60 vertices[21] = Point(0.47, 0.35, 210.0); vertices[22] = Point(0.45, 0.37, 220.0); vertices[23] = Point(0.41, 0.38, 230.0);
61 vertices[24] = Point(0.38, 0.37, 240.0); vertices[25] = Point(0.35, 0.36, 250.0); vertices[26] = Point(0.32, 0.35, 260.0);
62 vertices[27] = Point(0.30, 0.37, 270.0); vertices[28] = Point(0.28, 0.39, 280.0); vertices[29] = Point(0.25, 0.40, 290.0);
63 vertices[30] = Point(0.23, 0.39, 300.0); vertices[31] = Point(0.21, 0.37, 310.0); vertices[32] = Point(0.21, 0.34, 320.0);
64 vertices[33] = Point(0.23, 0.32, 330.0); vertices[34] = Point(0.24, 0.29, 340.0); vertices[35] = Point(0.27, 0.24, 350.0);
65 vertices[36] = Point(0.29, 0.21, 360.0); vertices[37] = Point(0.29, 0.18, 370.0); vertices[38] = Point(0.26, 0.16, 380.0);
66 vertices[39] = Point(0.24, 0.17, 390.0); vertices[40] = Point(0.23, 0.19, 400.0); vertices[41] = Point(0.24, 0.22, 410.0);
67 vertices[42] = Point(0.24, 0.25, 420.0); vertices[43] = Point(0.21, 0.26, 430.0); vertices[44] = Point(0.17, 0.26, 440.0);
68 vertices[45] = Point(0.12, 0.24, 450.0); vertices[46] = Point(0.07, 0.20, 460.0); vertices[47] = Point(0.03, 0.15, 470.0);
69 vertices[48] = Point(0.01, 0.10, 480.0); vertices[49] = Point(0.02, 0.07, 490.0);
70
71 // Mesh this polygon.
72
73 // Create a constrained Delaunay triangulation.
74 CDT cdt;
75
76 std::vector<Vertex_handle> vertex_handle(number_of_vertices);
77
78 // Insert vertices of the polygon as our initial point set.
79 for(int i = 0; i < number_of_vertices; ++i) vertex_handle[i] = cdt.insert(vertices[i]);
80
81 // Insert constraints - edges of the polygon - in order to mesh only the polygon's interior.
82 for(int i = 0; i < number_of_vertices; ++i) cdt.insert_constraint(vertex_handle[i], vertex_handle[(i + 1) % number_of_vertices]);
83
84 Mesher mesher(cdt);
85
86 // Set a criteria on how to mesh.
87 mesher.set_criteria(Criteria(0.01, 0.01));
88
89 // Mesh the polygon.
90 mesher.refine_mesh();
91
92 // Compute mean value coordinates and use them to interpolate data from the polygon's boundary to its interior.
93
94 // Associate each point with the corresponding function value and coordinates.
95 std::map<Point, Scalar, Projection::Less_xy_2> point_function_value;
96 std::vector< std::pair<Point, Scalar> > point_coordinates(number_of_vertices);
97
98 for(int i = 0; i < number_of_vertices; ++i)
99 point_function_value.insert(std::make_pair(vertices[i], vertices[i].z()));
100
101 // Create an instance of the class with mean value coordinates.
102 Mean_value_coordinates mean_value_coordinates(vertices.begin(), vertices.end());
103
104 // Store all new interior points with interpolated data here.
105 std::vector<Point> points(cdt.number_of_vertices());
106
107 cout << endl << "Result of the height interpolation: " << endl << endl;
108
109 // Compute coordinates and interpolate the boundary data to the polygon's interior.
110 int index = 0;
111 for(Vertex_iterator vertex_iterator = cdt.finite_vertices_begin(); vertex_iterator != cdt.finite_vertices_end(); ++vertex_iterator) {
112 Scalar_vector coordinates;
113
114 const Point &point = vertex_iterator->point();
115 mean_value_coordinates(point, std::back_inserter(coordinates));
116
117 for(int j = 0; j < number_of_vertices; ++j)
118 point_coordinates[j] = std::make_pair(vertices[j], coordinates[j]);
119
120 Scalar f = CGAL::linear_interpolation(point_coordinates.begin(), point_coordinates.end(), Scalar(1), Value_access(point_function_value));
121 points[index] = Point(point.x(), point.y(), f);
122 cout << "The interpolated height with index " << index << " is " << f << ";" << endl;
123 ++index;
124 }
125 cout << endl;
126
127 return EXIT_SUCCESS;
128 }
129