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