CGAL 5.1.4 - 2D and Surface Function Interpolation
User Manual

Authors
Julia Flötotto

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}\).

Natural Neighbor Coordinates

Most interpolation methods offered by this package rely on 2D natural and regular neighbor coordinates, which we describe in this section.

Introduction

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.

nn_coords.svg
Figure 95.1 2D example: \( \mathbf{x}\) has five natural neighbors \( \mathbf{p_1},\ldots, \mathbf{p_5}\). The natural neighbor coordinate \( \lambda_3(\mathbf{x})\) is the ratio of the area of the purple polygon, \( \pi_3(\mathbf{x})\), over the area of the total highlighted zone.

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 95.1.

Various papers ([3], [4], [6], [7], [8]) show that the natural neighbor coordinates exhibit the following properties:

  • \( \mathbf{x} = \sum_{i=1}^n \lambda_i(\mathbf{x}) \mathbf{p_i}\) (barycentric coordinate property).
  • For any \( i,j \leq n, \lambda_i(\mathbf{p_j}) = \delta_{ij}\), where \( \delta_{ij}\) is the Kronecker symbol.
  • \( \sum_{i=1}^n \lambda_i(\mathbf{x}) = 1\) (partition of unity property).

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:

  • \( \pi(\mathbf{x}) = \infty\)
  • \( \lambda_i(\mathbf{x}) = 0 \), for all data point \( \mathbf{p_i}\) of \( \mathcal{P}\) except for the two endpoints – say \( \mathbf{p}\) and \( \mathbf{q}\) – of the edge where \( x\) lies.

The natural neighbor coordinate of \( \mathbf{x}\) with respect to these endpoints \( \mathbf{p}\) and \( \mathbf{q}\) will be:

  • \( \lambda_p(\mathbf{x}) = \frac{\|\mathbf{x} - \mathbf{q}\| }{ \|\mathbf{q} - \mathbf{p}\|} \)
  • \( \lambda_q(\mathbf{x}) = \frac{\|\mathbf{x} - \mathbf{p}\| }{ \|\mathbf{q} - \mathbf{p}\|} \)

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}\).

Regular Neighbor Coordinates

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 Triangulation for an in-depth description of power diagrams.

rn_coords.svg
Figure 95.2 Illustration of regular neighbor coordinates. The point set is the same as in Figure 95.1 but weights have been added. These weights are shown using circles, with the radius of each circle being equal to the square root of the weight of the point.

Warning
Contrary to Voronoi diagrams, a weighted point \( p_i\) does not necessarily possess a non-empty cell in the power diagram of \( \mathcal{P}\) (with \( p_i\in\mathcal{P}\)). When this is the case, that point is then said to be hidden and all of its regular neighbor coordinates are null.

Implementation

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.

Computation of the Coordinates

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:

  • points (or vertices) with a non-null coordinate along with these coordinates \( \pi_i(\mathbf{x})\),
  • the normalization factor \( \pi(\mathbf{x})\),
  • a Boolean indicator on whether the coordinate computation was successful (which is equivalent to a Boolean indicating whether the query point lies in the convex hull or not).

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()).

Examples

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.

Example for Natural Neighbor Coordinates

The signature of all coordinate computation functions is about the same.
File Interpolation/nn_coordinates_2.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <iostream>
#include <iterator>
#include <utility>
#include <vector>
typedef K::FT Coord_type;
typedef K::Point_2 Point;
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
// Resulting points-coordinates pairs will be stored in an object of this type
typedef std::vector<std::pair<Point, Coord_type> > Point_coordinate_vector;
int main()
{
Delaunay_triangulation dt;
for(int y=0; y<3; ++y)
for(int x=0; x<3; ++x)
dt.insert(K::Point_2(x, y));
// coordinates computation
K::Point_2 p(1.2, 0.7); // query point
Point_coordinate_vector coords;
CGAL::natural_neighbor_coordinates_2(dt, p, std::back_inserter(coords));
if(!result.third)
{
std::cout << "The coordinate computation was not successful." << std::endl;
std::cout << "The point (" << p << ") lies outside the convex hull." << std::endl;
}
K::FT norm = result.second;
std::cout << "Coordinate computation successful." << std::endl;
std::cout << "Normalization factor: " << norm << std::endl;
std::cout << "Coordinates for point: (" << p << ") are the following: " << std::endl;
for(std::size_t i=0; i<coords.size(); ++i)
std::cout << " Point: (" << coords[i].first << ") coeff: " << coords[i].second << std::endl;
return EXIT_SUCCESS;
}

Example for Regular Neighbor Coordinates

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

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_2.h>
#include <CGAL/regular_neighbor_coordinates_2.h>
#include <iostream>
#include <iterator>
#include <vector>
#include <utility>
typedef K::FT FT;
typedef CGAL::Regular_triangulation_2<K> Regular_triangulation;
typedef Regular_triangulation::Bare_point Bare_point;
typedef Regular_triangulation::Weighted_point Weighted_point;
typedef std::vector<std::pair<Weighted_point, FT> > Point_coordinate_vector;
int main()
{
Regular_triangulation rt;
for (int y=0; y<3; ++y)
for (int x=0; x<3; ++x)
rt.insert(Weighted_point(Bare_point(x, y), 0. /*weight*/));
// coordinate computation
Weighted_point wp(Bare_point(1.2, 0.7), 2.);
Point_coordinate_vector coords;
CGAL::regular_neighbor_coordinates_2(rt, wp, std::back_inserter(coords));
if(!result.third)
{
std::cout << "The coordinate computation was not successful." << std::endl;
std::cout << "The point (" <<wp.point() << ") lies outside the convex hull." << std::endl;
}
K::FT norm = result.second;
std::cout << "Coordinate computation successful." << std::endl;
std::cout << "Normalization factor: " << norm << std::endl;
std::cout << "Coordinates for point: (" << wp << ") are the following: " << std::endl;
for(std::size_t i=0; i<coords.size(); ++i)
std::cout << " Point: (" << coords[i].first << ") coeff: " << coords[i].second << std::endl;
return EXIT_SUCCESS;
}

Formatting the Output Coordinates

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

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <iostream>
#include <iterator>
#include <utility>
#include <vector>
typedef K::FT Coord_type;
typedef CGAL::Triangulation_data_structure_2<Vb> Tds;
typedef CGAL::Delaunay_triangulation_2<K, Tds> Delaunay_triangulation;
// The functor 'Identity' matches anything to itself
typedef Delaunay_triangulation::Vertex_handle Vertex_handle;
// Resulting points-coordinates pairs are here stored in an object of this type:
typedef std::vector<std::pair<Vertex_handle, Coord_type> > Point_coordinate_vector;
int main()
{
Delaunay_triangulation dt;
for(int y=0; y<3; ++y)
for(int x=0; x<3; ++x)
dt.insert(K::Point_2(x, y));
// coordinates computation
K::Point_2 p(1.2, 0.7); // query point
Point_coordinate_vector coords;
// The functor Identity is passed to the method
CGAL::natural_neighbor_coordinates_2(dt, p, std::back_inserter(coords), Identity());
if(!result.third)
{
std::cout << "The coordinate computation was not successful." << std::endl;
std::cout << "The point (" << p << ") lies outside the convex hull." << std::endl;
}
// Assign the coordinates to the vertices
std::cout << "Coordinates for point: (" << p << ") are the following: " << std::endl;
for(std::size_t i=0; i<coords.size(); ++i)
{
Vertex_handle vh = coords[i].first;
vh->info() = coords[i].second;
std::cout << " Vertex: (" << vh->point() << ") coeff: " << vh->info() << std::endl;
}
return EXIT_SUCCESS;
}

Surface Natural Neighbor Coordinates and Surface Neighbors

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.

Introduction

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>.

Implementation

Voronoi Intersection Diagrams

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.

Natural Neighbor Coordinates on Surfaces

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.

Surface Neighbors

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.

Example for Surface Neighbor Coordinates

The example below describes the usage of the function CGAL::surface_neighbor_coordinates_3().


File Interpolation/surface_neighbor_coordinates_3.cpp

// example with random points on a sphere
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/algorithm.h>
#include <CGAL/Origin.h>
#include <CGAL/surface_neighbor_coordinates_3.h>
#include <iostream>
#include <iterator>
#include <vector>
typedef K::FT Coord_type;
typedef K::Point_3 Point_3;
typedef K::Vector_3 Vector_3;
typedef std::vector< std::pair< Point_3, K::FT > > Point_coordinate_vector;
int main()
{
int n = 100;
std::vector< Point_3> points;
points.reserve(n);
std::cout << "Generate " << n << " random points on a sphere." << std::endl;
CGAL::Random_points_on_sphere_3<Point_3> g(1);
std::copy_n(g, n, std::back_inserter(points));
Point_3 p(1, 0, 0);
Vector_3 normal(p - CGAL::ORIGIN);
std::cout << "Compute surface neighbor coordinates for " << p << std::endl;
Point_coordinate_vector coords;
CGAL::surface_neighbor_coordinates_3(points.begin(), points.end(),
p, normal,
std::back_inserter(coords),
K());
if(!result.third)
{
//Undersampling:
std::cout << "The coordinate computation was not successful." << std::endl;
return 0;
}
K::FT norm = result.second;
std::cout << "Testing the barycentric property " << std::endl;
Point_3 b(0, 0, 0);
for(const std::pair< Point_3, Coord_type >& pc : coords)
b = b + (pc.second/norm) * (pc.first - CGAL::ORIGIN);
std::cout << " weighted barycenter: " << b <<std::endl;
std::cout << " squared distance: " << CGAL::squared_distance(p,b) << std::endl;
return EXIT_SUCCESS;
}

Interpolation Methods

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.

Interpolation of Function Values

The interpolation functions presented below are used to interpolate function values.

Linear Precision Interpolation

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.

Sibson's C^1 Continuous Interpolant

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's C^1 Continuous Interpolant

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().

Quadratic Precision Interpolants

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().

Gradient Fitting

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:

Example for Linear Interpolation

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

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Interpolation_traits_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/interpolation_functions.h>
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
typedef K::FT Coord_type;
typedef K::Point_2 Point;
int main()
{
Delaunay_triangulation T;
typedef std::map<Point, Coord_type, K::Less_xy_2> Coord_map;
typedef CGAL::Data_access<Coord_map> Value_access;
Coord_map value_function;
Coord_type a(0.25), bx(1.3), by(-0.7);
for (int y=0 ; y<3 ; y++){
for (int x=0 ; x<3 ; x++){
K::Point_2 p(x,y);
T.insert(p);
value_function.insert(std::make_pair(p, a + bx*x + by*y));
}
}
//coordinate computation
K::Point_2 p(1.3, 0.34);
std::vector<std::pair<Point, Coord_type> > coords;
Coord_type norm = CGAL::natural_neighbor_coordinates_2(T, p, std::back_inserter(coords)).second;
Coord_type res = CGAL::linear_interpolation(coords.begin(), coords.end(), norm,
Value_access(value_function));
std::cout << "Tested interpolation on " << p << " interpolation: "
<< res << " exact: " << a + bx*p.x() + by*p.y() << std::endl;
return EXIT_SUCCESS;
}

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

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Interpolation_traits_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/interpolation_functions.h>
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
typedef K::Vector_3 Vector_3;
typedef K::Point_2 Point_2;
int main()
{
Delaunay_triangulation T;
typedef std::map<Point_2, Vector_3, K::Less_xy_2> Coord_map;
typedef CGAL::Data_access<Coord_map> Value_access;
Coord_map value_function;
for (int y=0 ; y<255 ; y++){
for (int x=0 ; x<255 ; x++){
K::Point_2 p(x,y);
T.insert(p);
value_function.insert(std::make_pair(p, Vector_3(x,y,1)));
}
}
//coordinate computation
K::Point_2 p(1.3, 0.34);
std::vector<std::pair<Point_2, double> > coords;
double norm = CGAL::natural_neighbor_coordinates_2(T, p, std::back_inserter(coords)).second;
Vector_3 res = CGAL::linear_interpolation(coords.begin(), coords.end(), norm,
Value_access(value_function));
std::cout << "Tested interpolation on " << p << " interpolation: "
<< res << std::endl;
return EXIT_SUCCESS;
}

Example for Sibson's C^1 Interpolation Scheme with Gradient Estimation


File Interpolation/sibson_interpolation_2.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.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 <map>
#include <utility>
#include <vector>
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
typedef K::FT Coord_type;
typedef K::Point_2 Point;
typedef std::map<Point, Coord_type, K::Less_xy_2> Point_value_map ;
typedef std::map<Point, K::Vector_2 , K::Less_xy_2> Point_vector_map;
int main()
{
Delaunay_triangulation T;
Point_value_map value_function;
Point_vector_map 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);
T.insert(p);
value_function.insert(std::make_pair(p,a + bx* x+ by*y + c*(x*x+y*y)));
}
}
sibson_gradient_fitting_nn_2(T, std::inserter(gradient_function,
gradient_function.begin()),
Traits());
for(const Point_vector_map::value_type& pv : gradient_function)
{
std::cout << pv.first << " " << pv.second << std::endl;
}
// coordinate computation
K::Point_2 p(1.6, 1.4);
std::vector< std::pair< Point, Coord_type > > coords;
Coord_type norm = CGAL::natural_neighbor_coordinates_2(T, p, std::back_inserter(coords)).second;
//Sibson interpolant: version without sqrt:
std::pair<Coord_type, bool> res =
coords.end(), norm, p,
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 gradients are provided." << std::endl
<< " You may resort to linear interpolation." << std::endl;
return EXIT_SUCCESS;
}

The example interpolation_2.cpp compares numerically the errors of the different interpolation functions with respect to a known function.

Example for Storing Values and Gradients in Vertices

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

#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 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;
};
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
: public std::iterator<std::output_iterator_tag, void, void, void, void>
{
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
}
}
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;
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;
}

Design and Implementation History

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.