|
CGAL 5.6 - 2D and Surface Function Interpolation
|
|
CGAL 5.6 - 2D and Surface Function Interpolation
|
This package implements various neighbor coordinate computation functions as well as different methods for scattered data interpolation. Computation of natural and regular neighbor coordinates in 2D Euclidean space is described in Section Natural Neighbor Coordinates. Coordinate and neighbor computations on surfaces are discussed in Section Surface Natural Neighbor Coordinates and Surface Neighbors. Finally, we describe the different interpolation methods offered by this package in Section Interpolation Methods.
Scattered data interpolation solves the following problem: given measures of a function on a set of discrete data points, how to interpolate this function at an arbitrary query point. More formally, let \( \mathcal{P}=\{\mathbf{p_1},\ldots ,\mathbf{p_n}\}\) be a set of \( n\) points in \( \mathbb{R}^2\) or \( \mathbb{R}^3\) and \( \Phi\) be a scalar function defined on the convex hull of \( \mathcal{P}\). We assume that the function values are known at the points of \( \mathcal{P}\), i.e. to each \( \mathbf{p_i} \in \mathcal{P}\), we associate \( z_i = \Phi(\mathbf{p_i})\). Knowledge of the gradient of \( \Phi\) at \( \mathbf{p_i}\) is also sometimes required. It is then denoted \( \mathbf{g_i}= \nabla \Phi(\mathbf{p_i})\). The interpolation is carried out at an arbitrary query point \( \mathbf{x}\) on the convex hull of \( \mathcal{P}\).
Most interpolation methods offered by this package rely on 2D natural and regular neighbor coordinates, which we describe in this section.
Natural neighbor interpolation was introduced by Sibson [9] to interpolate multivariate scattered data. Given a set of data points \( \mathcal{P}\), the natural neighbor coordinates associated to \( \mathcal{P}\) are defined from the Voronoi diagram of \( \mathcal{P}\). When simulating the insertion of a query point \( \mathbf{x}\) into the Voronoi diagram of \( \mathcal{P}\), the potential Voronoi cell of \( \mathbf{x}\) "steals" some parts from neighboring existing cells.
Let \( \pi(\mathbf{x})\) denote the volume of the potential Voronoi cell of \( \mathbf{x}\) and \( \pi_i(\mathbf{x})\) denote the volume of the sub-cell that would be stolen from the cell of \( \mathbf{p_i}\) by the cell of \( \mathbf{x}\). The natural neighbor coordinate of \( \mathbf{x}\) with respect to the data point \( \mathbf{p_i}\in \mathcal{P}\) is defined by
\[ \lambda_i(\mathbf{x}) = \frac{\pi_i(\mathbf{x})}{\pi(\mathbf{x})}. \]
A two-dimensional example is depicted in Figure 100.1.
Various papers ([3], [4], [6], [7], [8]) show that the natural neighbor coordinates exhibit the following properties:
When the query point x is located on the envelope of the convex hull of \( \mathcal{P}\), the potential Voronoi cell of x becomes infinite and:
The natural neighbor coordinate of \( \mathbf{x}\) with respect to these endpoints \( \mathbf{p}\) and \( \mathbf{q}\) will be:
Furthermore, Piper [7] shows that the coordinate functions are continuous in the convex hull of \( \mathcal{P}\) and continuously differentiable except on the data points \( \mathcal{P}\).
The previous definition naturally extends to weighted Voronoi diagrams. These diagrams, also known as power diagrams, are obtained by considering weighted points (the weight being a scalar) and considering a weighted distance, the power distance, defined between two weighted points \( (p, \omega_p) \) and \( (q, \omega_q) \) by \( \Pi( (p, \omega_p), (q, \omega_q) ) = pq^2 - \omega_p - \omega_q \). See this section of the package 2D Triangulations for an in-depth description of power diagrams.
The interpolation package of CGAL provides functions to compute natural and regular neighbor coordinates for \( 2D\) points. Refer to the reference pages natural_neighbor_coordinates_2() and regular_neighbor_coordinates_2(). In addition, the package provides functions to compute natural neighbor coordinates on well sampled point set surfaces. See Section Surface Natural Neighbor Coordinates and Surface Neighbors and the reference page surface_neighbor_coordinates_3() for further information.
Given a Delaunay triangulation or a regular triangulation, our implementation computes natural and regular neighbor coordinates in two steps. Firstly, the vertices in conflict with the query point (that is, the vertices from which the query point will "steal") are determined. Then, the areas \( \pi_i(\mathbf{x})\) are computed by triangulating the Voronoi sub-cells. The output is threefold:
Note that if the query point has already been located in the triangulation (for example using functions such as locate()) and/or the boundary edges of the conflict zone are already determined, alternative functions allow to avoid the re-computation (see CGAL::natural_neighbor_coordinates_2() and CGAL::regular_neighbor_coordinates_2()).
This section presents some examples of utilization of the natural and regular neighbor coordinates computation methods. Sections Example for Natural Neighbor Coordinates and Example for Regular Neighbor Coordinates show basic examples of the computation of natural and regular neighbor coordinates. Section Example for Regular Neighbor Coordinates shows a more advanced use case of the API, where a functor is passed to change the output.
The signature of all coordinate computation functions is about the same.
File Interpolation/nn_coordinates_2.cpp
For regular neighbor coordinates, it is sufficient to replace the name of the function and the type of triangulation passed as parameter, as shown in the example below:
File Interpolation/rn_coordinates_2.cpp
Internally, natural or regular neighbor coordinates are associated to the vertices of the triangulation via objects of type std::pair<Vertex_handle, Coord_type>, where Coord_type is the number type of the coordinates. By default and for backward compatibility reasons, this output is converted into objects of type std::pair<Point, Coord_type>, where Point is a bare (weightless) point.
It is however possible to collect the output as objects of any desired type, for example the simple Point type, by passing a functor as extra parameter of the coordinates computation function. The argument type of this functor must be std::pair<Vertex_handle, Coord_type> and the result type is chosen by the user, but must be consistent with the output iterator. Usage of this parameter is demonstrated in the example below, where the output is kept as objects of type std::pair<Vertex_handle, Coord_type>, which allows us to then store the coordinate in the vertex, using the class CGAL::Triangulation_vertex_base_with_info_2.
File Interpolation/nn_coordinates_with_info_2.cpp
This section introduces functions to compute natural neighbor coordinates and surface neighbors associated to a set of sample points issued from a surface \( \mathcal{S}\) and given a query point \( \mathbf{x}\) on \( \mathcal{S}\). We assume that \( \mathcal{S}\) is a closed and compact surface of \( \mathbb{R}^3\), and let \( \mathcal{P}= \{\mathbf{p_1}, \ldots,\mathbf{p_n}\}\) be an \( \epsilon\)-sample of \( \mathcal{S}\) (refer to Amenta and Bern [1]). The concepts are based on the definition of Boissonnat and Flötotto [2], [5]. Both references contain a thorough description of the requirements and the mathematical properties.
Two observations lead to the definition of surface neighbors and surface neighbor coordinates: First, it is clear that the tangent plane \( \mathcal{T}_x\) of the surface \( \mathcal{S}\) at the point \( \mathbf{x} \in \mathcal{S}\) approximates \( \mathcal{S}\) in the neighborhood of \( \mathbf{x}\). It has been shown in [2] that, if the surface \( \mathcal{S}\) is well sampled with respect to the curvature and the local thickness of \( \mathcal{S}\), i.e. it is an \( \epsilon\)-sample, then the intersection of the tangent plane \( \mathcal{T}_x\) with the Voronoi cell of \( \mathbf{x}\) in the Voronoi diagram of \( \mathcal{P} \cup \{\mathbf{x}\}\) has a small diameter. Consequently, inside this Voronoi cell, the tangent plane \( \mathcal{T}_x\) is a reasonable approximation of \( \mathcal{S}\). Furthermore, the second observation allows to compute this intersection diagram easily: one can show using Pythagoras' Theorem that the intersection of a three-dimensional Voronoi diagram with a plane \( \mathcal{H}\) is a two-dimensional power diagram. The points defining the power diagram are the projections of the points in \( \mathcal{P}\) onto \( \mathcal{H}\), each point weighted with its negative square distance to \( \mathcal{H}\). Algorithms for the computation of power diagrams via the dual regular triangulation are well known and for example provided by CGAL in the class Regular_triangulation_2<Gt, Tds>.
In CGAL, the regular triangulation dual to the intersection of a \( 3D\) Voronoi diagram with a plane \( \mathcal{H}\) can be computed by instantiating the Regular_triangulation_2<Gt, Tds> class with the traits class Voronoi_intersection_2_traits_3<K>. This traits class contains a point and a vector as class member which define the plane \( \mathcal{H}\). All predicates and constructions used by Regular_triangulation_2<Gt, Tds> are replaced by the corresponding operators on three-dimensional points. For example, the power test predicate (which takes three weighted \( 2D\) points \( p'\), \( q'\), \( r'\) of the regular triangulation and tests the power distance of a fourth point \( t'\) with respect to the power circle orthogonal to \( p\), \( q\), \( r\)) is replaced by a Side_of_plane_centered_sphere_2_3 predicate that tests the position of a \( 3D\) point \( t\) with respect to the sphere centered on the plane \( \mathcal{H}\) passing through the \( 3D\) points \( p\), \( q\), \( r\). This approach allows to avoid the explicit construction of the projected points and the weights, operations which are very prone to rounding errors.
The computation of natural neighbor coordinates on surfaces is based upon the computation of regular neighbor coordinates with respect to the regular triangulation that is dual to \( {\rm Vor}(\mathcal{P}) \cap \mathcal{T}_x\), the intersection of \( \mathcal{T}_x\) and the Voronoi diagram of \( \mathcal{P}\), via the function regular_neighbor_coordinates_2().
Of course, we might introduce all data points \( \mathcal{P}\) into this regular triangulation. However, this is not necessary because we are only interested in the cell of \( \mathbf{x}\). It is sufficient to guarantee that all surface neighbors of the query point \( \mathbf{x}\) are among the input points that are passed as argument to the function. The sample points \( \mathcal{P}\) can be filtered for example by distance, e.g. using range search or \( k\)-nearest neighbor queries, or with the help of the \( 3D\) Delaunay triangulation since the surface neighbors are necessarily a subset of the natural neighbors of the query point in this triangulation. CGAL provides a function that encapsulates the filtering based on the \( 3D\) Delaunay triangulation. For input points filtered by distance, functions are provided that indicate whether or not points that lie outside the input range (i.e.\ points that are further from \( \mathbf{x}\) than the furthest input point) can still influence the result. This allows to iteratively enlarge the set of input points until the range is sufficient to certify the result.
The surface neighbors of the query point are its neighbors in the regular triangulation that is dual to \( {\rm Vor}(\mathcal{P}) \cap \mathcal{T}_x\), the intersection of \( \mathcal{T}_x\) and the Voronoi diagram of \( \mathcal{P}\). As for surface neighbor coordinates, this regular triangulation is computed and the same kind of filtering of the data points as well as the certification described above is provided.
The example below describes the usage of the function CGAL::surface_neighbor_coordinates_3().
File Interpolation/surface_neighbor_coordinates_3.cpp
We describe in this section the different interpolation methods offered by this package. These can be regrouped into two large classes: values and gradient interpolation methods.
The interpolation functions presented below are used to interpolate function values.
Sibson [9] defines a very simple interpolant that re-produces linear functions exactly. The interpolation of \( \Phi(\mathbf{x})\) is given as the linear combination of the neighbors' function values weighted by the coordinates:
\[ Z^0(\mathbf{x}) = \ccSum{i}{}{ \lambda_i(\mathbf{x}) z_i}. \]
Indeed, if \( z_i=a + \mathbf{b}^t \mathbf{p_i}\) for all natural neighbors of \( \mathbf{x}\), we have
\[ Z^0(\mathbf{x}) = \ccSum{i}{}{ \lambda_i(\mathbf{x}) (a + \mathbf{b}^t\mathbf{p_i})} = a+\mathbf{b}^t \mathbf{x} \]
by the barycentric coordinate property.
This interpolation is implemented through the function CGAL::linear_interpolation(). The first example in Subsection Example for Linear Interpolation shows how the function is called.
In [9], Sibson describes a second interpolation method that relies also on the function gradient \( \mathbf{g_i}\) for all \( \mathbf{p_i} \in \mathcal{P}\). It is \( C^1\) continuous with gradient \( \mathbf{g_i}\) at \( \mathbf{p_i}\). Spherical quadrics of the form \( \Phi(\mathbf{x}) =a + \mathbf{b}^t \mathbf{x} +\gamma\ \mathbf{x}^t\mathbf{x}\) are reproduced exactly. The proof relies on the barycentric coordinate property of the natural neighbor coordinates and assumes that the gradient of \( \Phi\) at the data points is known or approximated from the function values as described in [9] (see Section Gradient Fitting).
Sibson's \( Z^1\) interpolant is a combination of the linear interpolant \( Z^0\) and an interpolant \( \xi\) which is the weighted sum of the first degree functions
\[ \xi_i(\mathbf{x}) = z_i +\mathbf{g_i}^t(\mathbf{x}-\mathbf{p_i}),\qquad \xi(\mathbf{x})= \frac{\ccSum{i}{}{ \frac{\lambda_i(\mathbf{x})} {\|\mathbf{x}-\mathbf{p_i}\|}\xi_i(\mathbf{x}) } }{\ccSum{i}{}{ \frac{\lambda_i(\mathbf{x})}{\|\mathbf{x}-\mathbf{p_i}\|}}}. \]
Sibson observed that the combination of \( Z^0\) and \( \xi\) reconstructs exactly a spherical quadric if they are mixed as follows:
\[ Z^1(\mathbf{x}) = \frac{\alpha(\mathbf{x}) Z^0(\mathbf{x}) + \beta(\mathbf{x}) \xi(\mathbf{x})}{\alpha(\mathbf{x}) + \beta(\mathbf{x})} \mbox{ where } \alpha(\mathbf{x}) = \frac{\ccSum{i}{}{ \lambda_i(\mathbf{x}) \frac{\|\mathbf{x} - \mathbf{p_i}\|^2}{f(\|\mathbf{x} - \mathbf{p_i}\|)}}}{\ccSum{i}{}{ \frac{\lambda_i(\mathbf{x})} {f(\|\mathbf{x} - \mathbf{p_i}\|)}}} \mbox{ and } \beta(\mathbf{x})= \ccSum{i}{}{ \lambda_i(\mathbf{x}) \|\mathbf{x} - \mathbf{p_i}\|^2}, \]
where in Sibson's original work, \( f(\|\mathbf{x} - \mathbf{p_i}\|) = \|\mathbf{x} - \mathbf{p_i}\|\).
This interpolation method can be used by calling CGAL::sibson_c1_interpolation().
CGAL contains a second implementation using \( f(\|\mathbf{x} - \mathbf{p_i}\|) = \|\mathbf{x} - \mathbf{p_i}\|^2\), CGAL::sibson_c1_interpolation_square(), which is less demanding on the number type because it avoids the square-root computation needed to compute the distance \( \|\mathbf{x} - \mathbf{p_i}\|\). The theoretical guarantees are the same (see [5]). Simply, the smaller the slope of \( f\) around \( f(0)\), the faster the interpolant approaches \( \xi_i\) as \( \mathbf{x} \) goes to \( \mathbf{p_i}\).
Farin [4] extended Sibson's work and realizes a \( C^1\) continuous interpolant by embedding natural neighbor coordinates in the Bernstein-Bézier representation of a cubic simplex. If the gradient of \( \Phi\) at the data points is known, this interpolant reproduces quadratic functions exactly. The function gradient can be approximated from the function values by Sibson's method [9] (see Section Gradient Fitting) which is exact only for spherical quadrics.
The Farin \( C^1\)-continuous interpolant is implemented in the function CGAL::farin_c1_interpolation().
Knowing the gradient \( \mathbf{g_i}\) for all \( \mathbf{p_i} \in \mathcal{P}\), we formulate a very simple interpolant that reproduces exactly quadratic functions. This interpolant is not \( C^1\) continuous in general. It is defined as follows:
\[ I^1(\mathbf{x}) = \ccSum{i}{}{ \lambda_i(\mathbf{x}) (z_i + \frac{1}{2} \mathbf{g_i}^t (\mathbf{x} - \mathbf{p_i}))} \]
This interpolation can be used with the method CGAL::quadratic_interpolation().
Sibson describes a method to approximate the gradient of the function \( f\) from the function values on the data sites. For the data point \( \mathbf{p_i}\), we determine
\[ \mathbf{g_i} = \min_{\mathbf{g}} \ccSum{j}{}{ \frac{\lambda_j(\mathbf{p_i})}{\|\mathbf{p_i} - \mathbf{p_j}\|^2} \left( z_j - (z_i + \mathbf{g}^t (\mathbf{p_j} -\mathbf{p_i})) \right)}, \]
where \( \lambda_j(\mathbf{p_i})\) is the natural neighbor coordinate of \( \mathbf{p_i}\) with respect to \( \mathbf{p_i}\) associated to \( \mathcal{P} \setminus \{\mathbf{p_i}\}\). For spherical quadrics, the result is exact.
CGAL provides functions to approximate the gradients of all data points that are inside the convex hull. There is one function for each type of neighbor coordinates:
natural_neighbor_coordinates_2(): CGAL::sibson_gradient_fitting_nn_2()regular_neighbor_coordinates_2(): CGAL::sibson_gradient_fitting_rn_2().The example below shows a typical example of the usage of CGAL::linear_interpolation(). Note the use of CGAL::Data_access to transform the std::map in a functor compatible with the function.
File Interpolation/linear_interpolation_2.cpp
The next example shows the interpolation of non-scalar values, namely Vector_3. It may be any type T which provides a multiplication operator with a scalar, and the addition operator for T.
File Interpolation/linear_interpolation_of_vector_3.cpp
File Interpolation/sibson_interpolation_2.cpp
The example interpolation_2.cpp compares numerically the errors of the different interpolation functions with respect to a known function.
In the previous examples, we have stored the values and gradients in a std::map and wrapped them in a CGAL::Data_access object. We can avoid this "external" storage by using the class Triangulation_vertex_base_with_info_2 provided by the triangulation package: values and gradients can be stored directly in the vertices of the triangulation. Functors and output iterators requested by the gradient fitting and interpolation functions are then simply wrappers providing read/write accesses to the info member of the vertices. This approach is shown in the example below.
File Interpolation/sibson_interpolation_vertex_with_info_2.cpp
The original version was written by Julia Flötotto, while working towards her PhD thesis at Inria. The possibility to use values and gradients in a Triangulation_vertex_base_with_info_2 was introduced by Andreas Fabri and Mael Rouxel-Labbé working at GeometryFactory.
1.8.13