Details

This Demonstration illustrates distance-based interpolation of five types; to make it visual, we use 2D points, which you can create and drag. You can change the method, some option settings, and the point locations. The controls interact with each other; inactive controls are hidden (e.g. for the method "Delaunay"). One of the control point locators has been selected to have a green dot in its center. If you move the "green point number" slider, or if you move another point locator just a little, that point instead is selected as the green point. There is a twofold use of this green point. For the "function values" setting "Interpolation coefficients", the contour plot shows the weight of the function value of the green control point at different interpolation points. This corresponds to setting the function value to 1 at the green control point and setting it to 0 at all other control points. For the "function values" setting "Random adjustable", the function values are random real numbers connected to the control points. The "green point function value" slider is provided to change the function value at the green point. For the function values setting "sin(x) sin(y)", the function value at each point is . If a control point is moved, the function value is adjusted to the new position. The green point has no special meaning for this last setting.

To keep updates fast, the settings of the contour plot routine are such that sometimes some jaggedness in the curves can be seen.

The "Generate a new point set!" button generates a whole set of new random points and function values.

As with all interpolation functions, interpolated values should be taken with many grains of salt whenever they are used to estimate the values of some underlying function. The basic property of interpolation functions is that they can reproduce the function values in the set of control points for the interpolation.

You might activate the autorun feature of the Manipulate cell. Autorun can be started from the little "+" in a gray circle in the upper-right corner of the Manipulate cell. Autorun here is a bit slow to start, doing nothing for 15 seconds or so, but then 210 different annotated views are shown in sequence.

Distance function choice

The "distance function" choices also require a little explanation. "Automatic" implies here the usual squared Euclidean distance. "Periodic symmetric" implies the distance function

,

and "Periodic unsymmetrical" implies the modulated variant

.

"Cut with asymmetry" means that we have a cut between and . Distances between points are never measured across the cut. They are instead measured along straight lines around the cut (passing one of the endpoints), and they are always measured counterclockwise from the first point to the second. These peculiar distance functions have been chosen to "stress-test" the mathematical routines, but at the same time they illustrate the power of distance-based interpolation, which can be used not only in multidimensional geometry, but also in strange topologies.

Delaunay interpolation

The Delaunay interpolation method is a coordinate-based interpolation method. The method is based on Delaunay triangulation of the control points and triangular interpolation inside the Delaunay triangles. The Mathematica routines were initially written by Tom Wickham–Jones and modified by Daniel Huber (Switzerland) and Frank Iannarilli. I have modified them further, by introducing an extrapolation method to avoid error messages for extrapolated values, and by changing the calls to the routine to agree with the Interpolation command.

The other four implemented methods are distance-based, and they interact with the control points via a distance matrix and a distance function. The distance-based interpolation methods could be applied to a large set of geometries, including multidimensional, periodic, and regions with cuts, as long as the relevant distance function is specified.

Voronoi interpolation

The Voronoi interpolation method to the Interpolation command returns the function value of the nearest control point from the interpolation point. If two control points are at the same minimum distance, the function value of the first is returned. The function is thus constant in regions. The function values could be almost any kind of valid Mathematica statements. This function is useful if you want to construct answers to questions like "What is the name and address of the closest restaurant?" but could also be useful in machine learning.

Shepard interpolation

The Shepard interpolation is included here for comparison and should not be used for serious work without careful consideration. The Shepard interpolation uses the distance to a control point (raised to the inverse ShepardPower parameter) as a weight for the function value in the control point. The weights are then normalized to have 1 as their sum. If the value of ShepardPower is less than the dimension of the point set plus one, the distant points will dominate the interpolated value, and this makes it difficult to use the method if the dimensionality of the control point set is high. However, what is important here is the real dimension of the control point set, which might be lower than the dimension of the embedding coordinate system (if we have such a system). The method gives nice linear interpolation between two points if ShepardPower is one, but only for this value. Thus we might conclude that the Shepard interpolation works best for interpolation in one-dimensional datasets. The Shepard interpolation method does not require any length scale parameter.

RBF interpolation

The RBF interpolation can be seen as an improvement of the Shepard interpolation. We let each control point be the center for one radial base function (a function where the value is determined by the distance to the center), and form a linear combination of all these functions. We choose the weights in this linear combination in such a way that it takes the function values in the control points. Since we have the same number of weights as we have function values, we should in general have a well-specified system of equations to solve, but only if the control points are well separated. Any duplicate points must be removed before applying the method. In this Demonstration the control points are kept apart. If you try to put one on top of another, they are pushed apart by a tiny amount. The radial base functions can be chosen in various ways, but a fundamental property should be that they are finite but nonzero at the origin and smoothly decreasing or increasing outward. In the definition of the radial base function, as a rule there is a length scale parameter, so if the length scale varies strongly, there might be a problem choosing one single appropriate value. For smooth functions the accuracy of the RBF interpolation can be amazingly good.

The RBF interpolation method is translated from routines described by William H. Press in the course Computational Statistics with Application to Bioinformatics at the University of Texas at Austin.

ObtuseAngle interpolation

The ObtuseAngle interpolation is based on the obtuse angle shadowing rule for connecting points in a point set.

The point is NOT connected to point if there is another point in the point set such that is connected to and the angle is obtuse.

This first variant of the rule is of the type "Directed", and the connections are shown as arrows in the Demonstration. The type "Undirected" is a second variant of the obtuse angle shadowing rule and reads:

"Two points and in the network are NOT connected if there are two other points and such that is connected to and is connected to and both the angles and are obtuse."

Undirected connections are shown as lines without arrowheads.

Further, two points in the network are not connected if the distance between them is zero, smaller than the smoothing distance, or larger than twice the cutoff radius. The cutoff radius is an optional parameter that you can vary with the "cutoff radius" slider. For small systems the cutoff radius can be large or infinite, but for large systems we might get better speed and scaling properties by an appropriate choice of the cutoff radius. With a hierarchical organization of the point, we could then avoid calculating and storing the complete distance matrix. (That is not included in this Demonstration.)

For normal Euclidean distances the two variants "Directed" and "Undirected" often result in exactly the same graphs, but especially for asymmetric distance functions, the difference can be large. For asymmetric distance functions, think about geometries with one timelike coordinate, or think about biking in Norway.

Note the difference between this kind of network and Delaunay triangulation: here some of the points are only connected to one of the other points. Occasionally two connections might cross. We can estimate the dimensionality of the local environment to a given point by looking at the number of connected points. If you look at a random infinite network, you can estimate the number of connected neighbors:

To experiment with these networks without interpolation, choose "method" to be "ObtuseAngle" and "interpolation order" to be "No interpolation". In this mode the response to changes in the point set is rapid.

The interpolation algorithm starts with a weight function for each control point that is one everywhere. Then, if there is a connection from control point 1 to control point 2, the weight function for control point 2 is multiplied by a "shadow" function, varying only in a direction parallel to the connection. In the case that "interpolation order" is 1, this function is zero from point 1 and beyond, one from point 2 and beyond, and increases linearly from point 1 to point 2. In this way, the weight function is zero at all points shadowed by point 1 from point 2, according to the obtuse angle shadowing rule. The same operation is applied to all connections. Then the weights for all control points are normalized to have 1 as their sum at the interpolation point. For "interpolation order" set to 1, we always know that the interpolation coefficients are in the closed interval 0 to 1. For some applications this is essential.

In the case that "interpolation order" is set to 2, the shadow function is chosen in such a way that it has a continuous derivative. At the same time it loses its ability to shadow points behind the from-point in the connections, and we have to introduce a second-level shadowing function for the second-level connections to make the weight function zero for more distant points. If "interpolation order" is set to 3, the shadow function has a second-order continuous derivative; the first-level, second-level, and third-level connections are used.

If there are two control points in the same position, and the function values are different, it is difficult to find interpolation functions that pass through the function values of both the points. If two control points instead are very close to each other, we might obtain very strong gradients. For the "ObtuseAngle" method we might specify that such points are not connected, and then the algorithm will do an averaging action in the neighborhood of the two close points. The "smoothing distance" option specifies a maximum distance for when we should consider two points to be close, and when the algorithm thus should consider them as unconnected.

As long as the setting "cutoff radius" is large and "smoothing distance" is zero, the ObtuseAngle interpolation method is scale-independent and does not need any scale parameter.

The zeroth-order interpolated function is piecewise constant; it is essentially the average of surrounding neighbor points.

The first-order interpolated function is continuous, but does not have a continuous derivative. The interpolation coefficients are always positive or strictly zero, so the interpolation always gives very "conservative" results, never beyond the control point values.

The second-order interpolated function has a continuous derivative with interpolation coefficients that might be negative. This function is able to suggest positions for minima and maxima. To calculate the second-order interpolation, both the first level and the second level of connected points are used. The second level of connected points is calculated in the same way as the first level, just ignoring all first-level connections.

The third-order interpolated function has a continuous second-order derivative.

The ObtuseAngle interpolation method is my invention, and the method has some qualities that should make it an interesting alternative to RBF interpolation in some applications.

Applications

The function values for the control points could be any kind of mathematical objects that are such that they can be linearly combined. For the Voronoi method, the function values could be any kind of mathematical objects. As long as we can define an appropriate distance function, we might find interpolation functions for periodic systems and systems with strange topology, where it might be difficult to define good coordinate systems.

These methods have probably been invented several times before, but I have not found any information about that in the easily available surface layer. Please inform me if you know! The programming is to be seen as "proof-of-principle", and I have tried to keep the programming simple. The algorithms could surely be improved in several aspects, especially the speed for a large number of points.

I find the simplicity and the generality of the methods attractive. One could think of different applications:

• Comparing DNA sequences

• Finding spell-checking suggestions

• Finding contextually close words from word correlation analysis

• Organizing sets of color points, for example, to be able to do dithering with a finite number of colors or reduction of the number of colors to a given set

• Finding topological features of experimental datasets: interior, exterior, and connecting points, clusters, voids, and so on

• Solving optimization problems

• Using the application in machine learning as a complement to neural networks: just add relevant points to the network

• Organizing datasets for Wolfram|Alpha and finding connections in them

You might download the Obtuse Mathematica package (requires Mathematica) and also browse the online documentation, which has many examples, at http://www.familydahl.se/mathematica, where the concepts are discussed in more detail.