cgal/Interpolation/doc/Interpolation/Interpolation.txt

namespace CGAL {
/*!

\mainpage User Manual 
\anchor Chapter_2D_and_Surface_Function_Interpolation
\anchor chapinterpolation
\cgalAutoToc

\authors 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 \ref seccoordinates, 
the functions concerning the coordinate and neighbor 
computation on surfaces are discussed in Section \ref secsurface. 
In Section \ref secinterpolation, 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 \f$ \mathcal{P}=\{\mathbf{p_1},\ldots ,\mathbf{p_n}\}\f$ be a set of
\f$ n\f$ points in \f$ \mathbb{R}^2\f$ or \f$ \mathbb{R}^3\f$ and \f$ \Phi\f$ be a scalar
function defined on the convex hull of \f$ \mathcal{P}\f$. We assume that
the function values are known at the points of \f$ \mathcal{P}\f$, i.e.\ to
each \f$ \mathbf{p_i} \in \mathcal{P}\f$, we associate \f$ z_i =
\Phi(\mathbf{p_i})\f$. Sometimes, the gradient of \f$ \Phi\f$ is also known
at \f$ \mathbf{p_i}\f$. It is denoted \f$ \mathbf{g_i}= \nabla
\Phi(\mathbf{p_i})\f$. The interpolation is carried out for an arbitrary query point
\f$ \mathbf{x}\f$ on the convex hull of \f$ \mathcal{P}\f$.

\section seccoordinates Natural Neighbor Coordinates 

\subsection InterpolationIntroduction Introduction

Natural neighbor interpolation has been introduced by Sibson
\cgalCite{s-bdnni-81} to interpolate multivariate scattered data. Given
a set of data points \f$ \mathcal{P}\f$, the natural neighbor coordinates
associated to \f$ \mathcal{P}\f$ are defined from the Voronoi diagram of
\f$ \mathcal{P}\f$. When simulating the insertion of a query point
\f$ \mathbf{x}\f$ into the Voronoi diagram of \f$ \mathcal{P}\f$, the potential
Voronoi cell of \f$ \mathbf{x}\f$ "steals" some parts from the existing
cells.

\cgalFigureBegin{fignn_coords,nn_coords.png}
2D example: \f$ \mathbf{x}\f$ has five natural neighbors \f$ \mathbf{p_1},\ldots, \mathbf{p_5}\f$. The natural neighbor coordinate \f$ \lambda_3(\mathbf{x})\f$ is the ratio of the area of the pink polygon, \f$ \pi_3(\mathbf{x})\f$, over the area of the total highlighted zone. 
\cgalFigureEnd

Let \f$ \pi(\mathbf{x})\f$ denote the volume of the potential Voronoi cell
of \f$ \mathbf{x}\f$ and \f$ \pi_i(\mathbf{x})\f$ denote the volume of the
sub-cell that would be stolen from the cell of \f$ \mathbf{p_i}\f$ by the
cell of \f$ \mathbf{x}\f$. The natural neighbor coordinate of \f$ \mathbf{x}\f$
with respect to the data point \f$ \mathbf{p_i}\in \mathcal{P}\f$ is defined by
\f[ 
\lambda_i(\mathbf{x}) =
\frac{\pi_i(\mathbf{x})}{\pi(\mathbf{x})}. \f]
A two-dimensional example
is depicted in \cgalFigureRef{fignn_coords}.

Various papers (\cgalCite{s-vidt-80}, \cgalCite{f-sodt-90},
\cgalCite{cgal:p-plcbd-93}, \cgalCite{b-scaps-97}, \cgalCite{hs-vbihc-00}) show that
the natural neighbor coordinates have the following properties:
<OL>
<LI>\f$ \mathbf{x} = \sum_{i=1}^n \lambda_i(\mathbf{x}) \mathbf{p_i}\f$
(barycentric coordinate property).
<LI>For any \f$ i,j \leq n, \lambda_i(\mathbf{p_j})=
\delta_{ij}\f$, where \f$ \delta_{ij}\f$ is the Kronecker symbol.
<LI>\f$ \sum_{i=1}^n \lambda_i(\mathbf{x}) = 1\f$ (partition of unity
property).
</OL>
For the case where the query point x is located on the envelope of the convex hull of \f$ \mathcal{P}\f$,
the potential Voronoi cell of x becomes infinite and :

\f$ \pi(\mathbf{x}) = \infty\f$

\f$ \lambda_i(\mathbf{x}) = 0 \f$ for all data point \f$ \mathbf{p_i}\f$ of \f$ \mathcal{P}\f$ except for the two endpoints, let's say \f$ \mathbf{p}\f$ and \f$ \mathbf{q}\f$ ,of the edge where x lies.

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

\f$ \lambda_p(\mathbf{x}) = \frac{\|\mathbf{x} - \mathbf{q}\| }{ \|\mathbf{q} - \mathbf{p}\|} \f$

\f$ \lambda_q(\mathbf{x}) = \frac{\|\mathbf{x} - \mathbf{p}\| }{ \|\mathbf{q} - \mathbf{p}\|} \f$

Furthermore, Piper \cgalCite{cgal:p-plcbd-93} shows that the coordinate
functions are continuous in the convex hull of \f$ \mathcal{P}\f$ and
continuously differentiable except on the data points \f$ \mathcal{P}\f$.<BR>

<BR>

The interpolation package of \cgal provides functions to compute
natural neighbor coordinates for \f$ 2D\f$ and \f$ 3D\f$ points with respect
to Voronoi diagrams as well as with respect to power diagrams (only
\f$ 2D\f$), i.e.\ for weighted points. Refer to the reference pages
`natural_neighbor_coordinates_2()`,
`sibson_natural_neighbor_coordinates_3()`
`laplace_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 \ref secsurface and the reference page
`surface_neighbor_coordinates_3()` for further information.

\subsection InterpolationImplementation Implementation

Given a Delaunay triangulation or a Regular triangulation, the
vertices in conflict with the query point are determined. The areas
\f$ \pi_i(\mathbf{x})\f$ are computed by triangulating the Voronoi
sub-cells. The normalization factor \f$ \pi(\mathbf{x})\f$ is also
returned. If the query point is already located and/or the boundary
edges of the conflict zone are already determined, alternative
functions allow to avoid the re-computation.

\subsubsection InterpolationExampleforNaturalNeighborCoordinates Example for Natural Neighbor Coordinates

The signature of all coordinate computation functions is about the
same.
\cgalExample{Interpolation/nn_coordinates_2.cpp}

\subsubsection InterpolationExampleforRegularNeighborCoordinates 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. A
special traits class is needed.
\cgalExample{Interpolation/rn_coordinates_2.cpp}
For surface neighbor coordinates, the surface normal at the query
point must be provided, see Section \ref secsurface.

\section secsurface Surface Natural Neighbor Coordinates and Surface Neighbors 

This section introduces the functions to compute natural neighbor
coordinates and surface neighbors associated to a set of sample points
issued from a surface \f$ \mathcal{S}\f$ and given a query point
\f$ \mathbf{x}\f$ on \f$ \mathcal{S}\f$. We suppose that \f$ \mathcal{S}\f$ is a
closed and compact surface of \f$ \mathbb{R}^3\f$, and let \f$ \mathcal{P}=
\{\mathbf{p_1}, \ldots,\mathbf{p_n}\}\f$ be an \f$ \epsilon\f$-sample of
\f$ \mathcal{S}\f$ (refer to Amenta and Bern \cgalCite{ab-srvf-99}). The
concepts are based on the definition of Boissonnat and Fl&ouml;totto
\cgalCite{bf-lcss-02}, \cgalCite{cgal:f-csapc-03}. Both references
contain a thorough description of the requirements and the
mathematical properties.

\subsection InterpolationIntroduction_1 Introduction

Two observations lead to the definition of surface neighbors and
surface neighbor coordinates: First, it is clear that the tangent
plane \f$ \mathcal{T}_x\f$ of the surface \f$ \mathcal{S}\f$ at the point
\f$ \mathbf{x} \in \mathcal{S}\f$ approximates \f$ \mathcal{S}\f$ in the
neighborhood of \f$ \mathbf{x}\f$. It has been shown in \cgalCite{bf-lcss-02}
that, if the surface \f$ \mathcal{S}\f$ is well sampled with respect to the
curvature and the local thickness of \f$ \mathcal{S}\f$, i.e.\ it is an \f$ \epsilon\f$-sample, the intersection
of the tangent plane \f$ \mathcal{T}_x\f$ with the Voronoi cell of
\f$ \mathbf{x}\f$ in the Voronoi diagram of \f$ \mathcal{P} \cup
\{\mathbf{x}\}\f$ has a small diameter. Consequently, inside this
Voronoi cell, the tangent plane \f$ \mathcal{T}_x\f$ is a reasonable
approximation of \f$ \mathcal{S}\f$. 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 \f$ \mathcal{H}\f$ is a two-dimensional power
diagram. The points defining the power diagram are the projections of
the points in \f$ \mathcal{P}\f$ onto \f$ \mathcal{H}\f$, each point weighted
with its negative square distance to \f$ \mathcal{H}\f$. 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>`.

\subsection InterpolationImplementation_1 Implementation

\subsubsection InterpolationVoronoiIntersectionDiagrams Voronoi Intersection Diagrams

In \cgal, the regular triangulation dual to the intersection of a \f$ 3D\f$
Voronoi diagram with a plane \f$ \mathcal{H}\f$ 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 \f$ \mathcal{H}\f$. 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 \f$ 2D\f$ points
\f$ p'\f$, \f$ q'\f$, \f$ r'\f$ of the regular triangulation and tests the power
distance of a fourth point \f$ t'\f$ with respect to the power circle orthogonal
to \f$ p\f$, \f$ q\f$, \f$ r\f$) is replaced by a
`Side_of_plane_centered_sphere_2_3` predicate that tests the
position of a \f$ 3D\f$ point \f$ t\f$ with respect to the sphere centered on
the plane \f$ \mathcal{H}\f$ passing through the \f$ 3D\f$ points \f$ p\f$, \f$ q\f$, \f$ r\f$.
This approach allows to avoid the explicit constructions of the
projected points and the weights which are very prone to rounding
errors.

\subsubsection InterpolationNaturalNeighborCoordinateson 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 \f$ {\rm Vor}(\mathcal{P}) \cap
\mathcal{T}_x\f$, the intersection of \f$ \mathcal{T}_x\f$ and the Voronoi
diagram of \f$ \mathcal{P}\f$, via the function
`regular_neighbor_coordinates_2()`.

Of course, we might introduce all data points \f$ \mathcal{P}\f$ into this
regular triangulation. However, this is not necessary because we are
only interested in the cell of \f$ \mathbf{x}\f$. It is sufficient to
guarantee that all surface neighbors of the query point \f$ \mathbf{x}\f$
are among the input points that are passed as argument to the
function. The sample points \f$ \mathcal{P}\f$ can be filtered for example
by distance, e.g. using range search or \f$ k\f$-nearest neighbor queries,
or with the help of the \f$ 3D\f$ 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 \f$ 3D\f$ 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 \f$ \mathbf{x}\f$ 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.

\subsubsection InterpolationSurfaceNeighbors Surface Neighbors

The surface neighbors of the query point are its neighbors in the
regular triangulation that is dual to \f$ {\rm Vor}(\mathcal{P}) \cap
\mathcal{T}_x\f$, the intersection of \f$ \mathcal{T}_x\f$ and the Voronoi
diagram of \f$ \mathcal{P}\f$. 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.

\subsection InterpolationExampleforSurfaceNeighborCoordinates Example for Surface Neighbor Coordinates

\cgalExample{Interpolation/surface_neighbor_coordinates_3.cpp}

\section secinterpolation Interpolation Methods 

\subsection InterpolationIntroduction_2 Introduction

\subsubsection InterpolationLinearPrecisionInterpolation Linear Precision Interpolation

Sibson \cgalCite{s-bdnni-81} defines a very simple interpolant that
re-produces linear functions exactly. The interpolation of
\f$ \Phi(\mathbf{x})\f$ is given as the linear combination of the neighbors' function
values weighted by the coordinates:
\f[ 
Z^0(\mathbf{x}) = \ccSum{i}{}{ \lambda_i(\mathbf{x}) z_i}.
\f]
Indeed, if \f$ z_i=a + \mathbf{b}^t \mathbf{p_i}\f$ for all natural
neighbors of \f$ \mathbf{x}\f$, we have
\f[ Z^0(\mathbf{x}) = \ccSum{i}{}{ \lambda_i(\mathbf{x}) (a + \mathbf{b}^t\mathbf{p_i})} = a+\mathbf{b}^t \mathbf{x} \f]
by the barycentric coordinate property. The first example in
Subsection \ref subsecinterpol_examples shows how the function is
called.

\subsubsection InterpolationSibson Sibson's C^1 Continuous Interpolant

In \cgalCite{s-bdnni-81}, Sibson describes a second interpolation method
that relies also on the function gradient \f$ \mathbf{g_i}\f$ for all \f$ \mathbf{p_i} \in \mathcal{P}\f$. It is \f$ C^1\f$ continuous with gradient \f$ \mathbf{g_i}\f$ at
\f$ \mathbf{p_i}\f$. Spherical quadrics of the form \f$ \Phi(\mathbf{x}) =a +
\mathbf{b}^t \mathbf{x} +\gamma\ \mathbf{x}^t\mathbf{x}\f$ are reproduced 
exactly. The
proof relies on the barycentric coordinate property of the natural
neighbor coordinates and assumes that the gradient of \f$ \Phi\f$ at the
data points is known or approximated from the function values as
described in \cgalCite{s-bdnni-81} (see Section \ref sgradient_fitting).

Sibson's \f$ Z^1\f$ interpolant is a combination of the linear interpolant
\f$ Z^0\f$ and an interpolant \f$ \xi\f$ which is the weighted sum of the first
degree functions
\f[ \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}\|}}}. \f]
Sibson observed that the combination of \f$ Z^0\f$ and \f$ \xi\f$ reconstructs exactly
a spherical quadric if they are mixed as follows:
\f[ 
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}, \f]
where in Sibson's original work,
\f$ f(\|\mathbf{x} - \mathbf{p_i}\|) = \|\mathbf{x} - \mathbf{p_i}\|\f$.

\cgal contains a second implementation with \f$ f(\|\mathbf{x} -
\mathbf{p_i}\|) = \|\mathbf{x} - \mathbf{p_i}\|^2\f$ which is less
demanding on the number type because it avoids the square-root
computation needed to compute the distance \f$ \|\mathbf{x} -
\mathbf{p_i}\|\f$. The theoretical guarantees are the same (see
\cgalCite{cgal:f-csapc-03}). Simply, the smaller the slope of \f$ f\f$
around \f$ f(0)\f$, the faster the interpolant approaches \f$ \xi_i\f$ as
\f$ \mathbf{x} \rightarrow \mathbf{p_i}\f$.

\subsubsection InterpolationFarin Farin's C^1 Continuous Interpolant

Farin \cgalCite{f-sodt-90} extended Sibson's work and realizes a \f$ C^1\f$
continuous interpolant by embedding natural neighbor coordinates in
the Bernstein-B&eacute;zier representation of a cubic simplex. If the
gradient of \f$ \Phi\f$ 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
\cgalCite{s-bdnni-81} (see Section \ref sgradient_fitting) which is exact only
for spherical quadrics.

\subsubsection InterpolationQuadraticPrecisionInterpolants Quadratic Precision Interpolants

Knowing the gradient \f$ \mathbf{g_i}\f$ for all \f$ \mathbf{p_i} \in
\mathcal{P}\f$, we formulate a very simple interpolant that reproduces
exactly quadratic functions. This interpolant is not \f$ C^1\f$ continuous
in general. It is defined as follows:
\f[ 
I^1(\mathbf{x}) = \ccSum{i}{}{ \lambda_i(\mathbf{x}) 
(z_i + \frac{1}{2} \mathbf{g_i}^t (\mathbf{x} - \mathbf{p_i}))} 
\f]

\subsection sgradient_fitting Gradient Fitting 

Sibson describes a method to approximate the gradient of the function
\f$ f\f$ from the function values on the data sites. For the data point
\f$ \mathbf{p_i}\f$, we determine
\f[ \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)},
\f]
where \f$ \lambda_j(\mathbf{p_i})\f$ is the natural neighbor coordinate
of \f$ \mathbf{p_i}\f$ with respect to \f$ \mathbf{p_i}\f$ associated to
\f$ \mathcal{P} \setminus \{\mathbf{p_i}\}\f$. 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 natural neighbor coordinate (i.e.\ `natural_neighbor_coordinates_2()`, `regular_neighbor_coordinates_2()`).

\subsection subsecinterpol_examples Example for Linear Interpolation 

\cgalExample{Interpolation/linear_interpolation_2.cpp}

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

\cgalExample{Interpolation/sibson_interpolation_2.cpp}

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.

*/ 
} /* namespace CGAL */