2D and Surface Function Interpolation

*Julia Flötotto*

This chapter describes CGAL's interpolation package which implements natural neighbor coordinate functions as well as different methods for scattered data interpolation most of which are based on natural neighbor coordinates. The functions for computing natural neighbor coordinates in Euclidean space are described in Section 52.1, the functions concerning the coordinate and neighbor computation on surfaces are discussed in Section 52.2. In Section 52.3, we describe the different interpolation functions.

Scattered data interpolation solves the following problem: given
measures of a function on a set of discrete data points, the task is
to interpolate this function on an arbitrary query point.
More formally, let $$*={ p_{1},...,p_{n}}* be a set of
$$

**Figure 52.1: **$$*2D* example: $$* x* has five natural neighbors
$$

Let $$*( x)* denote the volume of the potential Voronoi cell
of $$

$$*
_{i}(x) =
(_{i}(x))/((x)).*

A two-dimensional example is depicted in Figure 52.1.

Various papers ([Sib80], [Far90], [Pip93], [Bro97], [HS00]) show that the natural neighbor coordinates have the following properties:

**(i)**- $$
(barycentric coordinate property).**x**=_{i=1}^{n}_{i}(**x**)**p**_{i} **(ii)**- For any $$
*i,j n,*, where $$_{i}(**p**)=_{j}_{ij}is the Kronecker symbol._{ij} **(iii)**- $$
_{i=1}^{n}_{i}(**x**) = 1

$$*( x) = *

$$* _{i}(x) = 0 * for all data point $$

The natural neighbor coordinate of $$* x* with respect to these endpoints $$

$$_{p}(**x**) = ( **x** - **q** )/( **q** - **p** )

$$_{q}(**x**) = ( **x** - **p** )/( **q** - **p** )

Furthermore, Piper [Pip93] shows that the coordinate
functions are continuous in the convex hull of $$ and
continuously differentiable except on the data points $$.

The interpolation package of CGAL provides functions to compute
natural neighbor coordinates for $$*2D* and $$*3D* points with respect
to Voronoi diagrams as well as with respect to power diagrams (only
$$*2D*), i.e. for weighted points. Refer to the reference pages
*natural_neighbor_coordinates_2*,
*natural_neighbor_coordinates_3* and
*regular_neighbor_coordinates_2*.

In addition, the package provides functions to compute natural
neighbor coordinates on well sampled point set surfaces. See
Section 52.2 and the reference page
*surface_neighbor_coordinates_3* for further information.

The signature of all coordinate computation functions is about the same.

File:examples/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> struct K : CGAL::Exact_predicates_inexact_constructions_kernel {}; typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation; typedef std::vector< std::pair< K::Point_2, K::FT > > 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)); //coordinate computation K::Point_2 p(1.2, 0.7); Point_coordinate_vector coords; CGAL::Triple< std::back_insert_iterator<Point_coordinate_vector>, K::FT, bool> result = 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 << "done" << std::endl; return 0; }

File:examples/Interpolation/rn_coordinates_2.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Regular_triangulation_2.h> #include <CGAL/Regular_triangulation_euclidean_traits_2.h> #include <CGAL/regular_neighbor_coordinates_2.h> struct K : CGAL::Exact_predicates_inexact_constructions_kernel {}; typedef CGAL::Regular_triangulation_euclidean_traits_2<K> Gt; typedef CGAL::Regular_triangulation_2<Gt> Regular_triangulation; typedef Regular_triangulation::Weighted_point Weighted_point; typedef std::vector< std::pair< Weighted_point, K::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(K::Point_2(x,y), 0)); //coordinate computation Weighted_point wp(K::Point_2(1.2, 0.7),2); Point_coordinate_vector coords; CGAL::Triple< std::back_insert_iterator<Point_coordinate_vector>, K::FT, bool> result = 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 << "done" << std::endl; return 0; }For surface neighbor coordinates, the surface normal at the query point must be provided, see Section 52.2.

This section introduces the functions to compute natural neighbor
coordinates and surface neighbors associated to a set of sample points
issued from a surface $$ and given a query point
$$* x* on $$. We suppose that $$ is a
closed and compact surface of $$

Two observations lead to the definition of surface neighbors and
surface neighbor coordinates: First, it is clear that the tangent
plane $$* _{x}* of the surface $$ at the point
$$

In CGAL, the regular triangulation dual to the intersection of a $$*3D*
Voronoi diagram with a plane $$ 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 $$. 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 $$ passing through the $$*3D* points $$*p*, $$*q*, $$*r*.
This approach allows to avoid the explicit constructions of the
projected points and the weights 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 $$*Vor() _{x}*, the intersection of $$

Of course, we might introduce all data points $$ into this
regular triangulation. However, this is not necessary because we are
only interested in the cell of $$* x*. It is sufficient to
guarantee that all surface neighbors of the query point $$

The surface neighbors of the query point are its neighbors in the
regular triangulation that is dual to $$*Vor() _{x}*, the intersection of $$

File:examples/Interpolation/surface_neighbor_coordinates_3.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/point_generators_3.h> #include <CGAL/copy_n.h> #include <CGAL/Origin.h> #include <CGAL/surface_neighbor_coordinates_3.h> struct K : CGAL::Exact_predicates_inexact_constructions_kernel {}; 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); CGAL::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::Triple< std::back_insert_iterator<Point_coordinate_vector>, K::FT, bool> result = 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(std::vector< std::pair< Point_3, Coord_type > >::const_iterator it = coords.begin(); it!=coords.end(); ++it) b = b + (it->second/norm)* (it->first - CGAL::ORIGIN); std::cout <<" weighted barycenter: " << b <<std::endl; std::cout << " squared distance: " << CGAL::squared_distance(p,b) <<std::endl; std::cout << "done" << std::endl; return 0; }

Sibson [Sib81] defines a very simple interpolant that
re-produces linear functions exactly. The interpolation of
$$*( x)* is given as the linear combination of the neighbors' function
values weighted by the coordinates:

$$*
Z ^{0}(x) = _{i} _{i}(x) z_{i}.
*

Indeed, if $$*z _{i}=a + b^{t} p_{i}* for all natural
neighbors of $$

$$* Z ^{0}(x) = _{i} _{i}(x) (a + b^{t}p_{i}) = a+b^{t} x*

by the barycentric coordinate property. The first example in Subsection 52.3.3 shows how the function is called.

Sibson's $$*Z ^{1}* interpolant is a combination of the linear interpolant
$$

$$_{i}(**x**) = z_{i}
+**g _{i}**

Sibson observed that the combination of $$*Z ^{0}* and $$ reconstructs exactly
a spherical quadric if they are mixed as follows:

$$*
Z ^{1}(x) = ((x) Z^{0}(x) +
(x) (x))/((x) +
(x)) where (x) =
( _{i} _{i}(x) ( x -
p_{i} ^{2})/(f( x - p_{i} )))/( _{i}
(_{i}(x))/(f( x - p_{i} )))
and (x)= _{i} _{i}(x)
x - p_{i} ^{2},*

where in Sibson's original work,
$$*f( x - p_{i} ) = x - p_{i} *.

CGAL contains a second implementation with $$*f( x -
p_{i} ) = x - p_{i} ^{2}* which is less
demanding on the number type because it avoids the square-root
computation needed to compute the distance $$

Farin [Far90] 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 $$ 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
[Sib81] (see Section 52.3.2) which is exact only
for spherical quadrics.

Knowing the gradient $$* g_{i}* for all $$

$$*
I ^{1}(x) = _{i} _{i}(x)
(z_{i} + (1)/(2) g_{i}^{t} (x - p_{i}))
*

$$* g_{i}
= *min$$

where $$* _{j}(p_{i})* is the natural neighbor coordinate
of $$

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 natural neighbor coordinate (i.e. *natural_neighbor_coordinates_2*, *regular_neighbor_coordinates_2*).

File:examples/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> struct K : CGAL::Exact_predicates_inexact_constructions_kernel {}; typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation; typedef CGAL::Interpolation_traits_2<K> Traits; typedef K::FT Coord_type; typedef K::Point_2 Point; int main() { Delaunay_triangulation T; std::map<Point, Coord_type, K::Less_xy_2> function_values; typedef CGAL::Data_access< std::map<Point, Coord_type, K::Less_xy_2 > > Value_access; 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); function_values.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(function_values)); std::cout << " Tested interpolation on " << p << " interpolation: " << res << " exact: " << a + bx* p.x()+ by* p.y()<< std::endl; std::cout << "done" << std::endl; return 0; }

File:examples/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> struct K : CGAL::Exact_predicates_inexact_constructions_kernel {}; typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation; typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits; 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 function_values; Point_vector_map function_gradients; //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); function_values.insert(std::make_pair(p,a + bx* x+ by*y + c*(x*x+y*y))); } sibson_gradient_fitting_nn_2(T,std::inserter(function_gradients, function_gradients.begin()), CGAL::Data_access<Point_value_map> (function_values), Traits()); //coordiante 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 = CGAL::sibson_c1_interpolation_square (coords.begin(), coords.end(),norm,p, CGAL::Data_access<Point_value_map>(function_values), CGAL::Data_access<Point_vector_map>(function_gradients), 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; std::cout << "done" << std::endl; return 0; }

An additional example in the distribution compares numerically the errors of the different interpolation functions with respect to a known function. It is distributed in the examples directory.