MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Unstructured Interpolation in V7?

  • To: mathgroup at
  • Subject: [mg109698] Re: Unstructured Interpolation in V7?
  • From: "Ingolf Dahl" <ingolf.dahl at>
  • Date: Wed, 12 May 2010 07:34:38 -0400 (EDT)

Hi Luc and all you other, who have asked for scattered point interpolation.
I am developing an interpolation package for Mathematica, named "Obtuse",
which contains four methods for doing scattered point interpolation. Now I
am looking for "test pilots" for my package. There are still some small
details in the documentation that need to fixed, but I could send a
preliminary version within a few days. Please send an email if you are
interested (to ingolf.dahl at

If I want to interpolate using the testData set by Luc, I write

interpol = 
  Method -> "ObtuseAngle"]; interpol[{-122.000`, 38.000`}]

Running this I obtain

{0.0384003, 0.0490968} 

In the testData set from Luc, only x and y coordinates are specified. It
would be possible to include the z or t coordinates also, or both z and t.

The four interpolation methods are defined as extra options to the
Interpolation function:

The Delaunay interpolation method (option Method->"Delaunay") to the
Interpolation command is a coordinate-based interpolation method, in my
package allowing only 2D-points. The method is based on Delaunay
triangulation (from the Computational Geometry Package) 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 three 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 multidimensional
points or to other types of objects, as long as the relevant distance
function is specified.

The Shepard interpolation (option Method->"Shepard") is included just for
comparison, and should not be used for serious work without careful
consideration. The Shepard interpolation is using 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 dimensionality of the control point set, which might be lower than the
dimensionality 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 is working best for interpolation in
one-dimensional data sets. The Shepard interpolation method does not require
any length scale parameter.

The Radial Basis Function interpolation (option Method->"RBF") 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 weight as we have function values, we should in
general have a well-specified equation system to solve. The radial base
functions could be chosen in various ways, but a fundamental property should
be that they are finite but non-zero at origin and smoothly decreasing or
increasing outwards. In the definition of the radial base function there is
as a rule a length scale parameter included, so if the length scale vary
strongly there might be a problem to choose one single appropriate value.
For smooth functions the accuracy of the RBF interpolation can be amazingly

The RBF interpolation method is translated from another system routines
described by William Press in the course "Computational Statistics with
Application to Bioinformatics".

The ObtuseAngle interpolation (option Method->"ObtuseAngle") is based on the
ObtuseAngle connection list. One 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
two is multiplied by a "shadow" function, varying only in a direction
parallel to the connection. In the case that InterpolationOrder is one, 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
will be zero in 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 weight for all control point are normalized to have 1
as their sum at the interpolation point. For InterpolationOrder equal to
one, 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 InterpolationOrder is two, 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 zero the weight function for more distant points. If
InterpolationOrder is three, the shadow function has second order continuous
derivative, both the first-level, second-level and the third-level
connections are used.

As long as the options CutoffRadius and SmoothenDistance are not set, the
ObtuseAngle interpolation method is scale-independent and does not need any
scale parameter.

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

The function values for the control points could be any kind of mathematical
objects that are such that they can be linearly combined.

The Delaunay and ObtuseAngle interpolation methods are typically local
interpolation methods, where the function value in one control point often
influence only the local neighborhood. The RBF interpolation is more of a
global interpolation method, where the function value in one control point
might influence the interpolated value almost everywhere.

The Obtuse package also provides Mathematica functions to find and
investigate networks created according to the obtuse-angle shadowing rule.

Ingolf Dahl
ingolf.dahl at

> -----Original Message-----
> From: Luc Barthelet [mailto:luc at]
> Sent: den 11 maj 2010 12:25
> To: mathgroup at
> Subject: [mg109658] Unstructured Interpolation in V7?
> I have a small data set  (wind in the bay area, automatically extracted
> by Mathematica)
> testData={{{-121.082`,38.955`},{-
> 0.017250000000004206`,0.02987787643056=
> 0946`}},{{-
> 122.05`,37.98`},{0.03175000000000239`,0.054992613140314006`}},=
> {{-122.12`,37.66`},{0.06795173495784468`,-0.011981724259015891`}},{{-
> 121.=
> 82`,37.7`},{0.06795173495784468`,0.011981724259015891`}},{{-
> 121.3`,38.58`=
> },{0.02599999999999625`,0.04503332099679369`}},{{-
> 122.05`,37.43`},{0.0380=
> 57551141832846`,0.01385181580469208`}},{{-
> 122.28`,38.21`},{0.070476946558=
> 94077`,0.02565151074942662`}},{{-
> 122.556`,38.1436`},{0.05121000315664048`=
> ,0.009029705238681629`}},{{-
> 121.957`,38.378`},{0.03241939541710792`,0.011=
> 79969494473454`}},{{-
> 122.22`,37.7`},{0.09060231327711676`,0.0159756323453=
> 5926`}},{{-122.12`,37.48`},{0.05662644579820153`,-
> 0.009984770215851313`}}=
> ,{{-121.6`,38.7`},{0.018640840680916426`,0.022215288850453874`}},{{-
> 121.5=
> `,38.52`},{0.02603289819229815`,0.031024799946315795`}},{{-
> 122.25`,37.53`=
> },{0.05121000315664048`,0.009029705238681629`}},{{-
> 122.36`,37.62`},{0.051=
> 21000315664048`,0.009029705238681629`}},{{-
> 121.93`,37.37`},{0.06495190528=
> 383432`,0.03750000000000142`}},{{-
> 121.82`,37.34`},{0.06795173495784468`,0=
> =2E011981724259015891`}},{{-
> 122.81`,38.51`},{0.03983716857408126`,0.02300=
> 000000000324`}},{{-
> 121.25`,37.91`},{0.05121000315664048`,0.00902970523868=
> 1629`}},{{-
> 121.93`,38.28`},{0.09725818625133797`,0.03539908483420362`}},{=
> {-122.465`,37.807`},{0.03464101615138304`,0.020000000000003126`}},{{-
> 122.=
> 21`,37.507`},{0.07386058147591257`,0.013023613325017891`}},{{-
> 122.298`,37=
> =2E772`},{0.05967048141990006`,-0.021718279101179405`}},{{-
> 122.33`,37.8`}=
> ,{0.1039999999999992`,0.`}},{{-
> 122.4`,37.928`},{0.029877876430560946`,0.0=
> 172499999999971`}},{{-122.975`,37.997`},{0.04530115663855838`,-
> 0.00798781=
> 617267963`}},{{-
> 122.038`,38.057`},{0.048643822138060955`,0.04081701321509=
> 712`}},{{-122.881`,37.361`},{0.06598211945181731`,-
> 0.011634427903686628`}=
> },{{-122.833`,37.759`},{0.055999999999997385`,0.`}},{{-
> 122.833`,37.759`},=
> {0.055999999999997385`,0.`}},{{-
> 122.026`,38.223`},{0.07924132444627219`,0=
> =2E04574999999999818`}},{{-
> 122.447`,37.891`},{0.005035797152345367`,0.028=
> 559424837354186`}}};
> That interpolate, extrapolate nicely with
> ListVectorPlot[testData]
> Now, I would like to build the interpolated data outside of
> ListVectorPlot,
> but Interpolation[testData] returns an error:
> Interpolation::indim: The coordinates do not lie on a structured tensor
> product grid
> Which is an old issue (for MathGroup at least).
> The solution are to use the IMS package imsUnstructuredInterpolation or
> TWJ's TriangularInterpolation. (or Cases on ListVectorPlot.., which is
> my current solution)
> But is there a new better way to do that in version 7? the
> documentation
> avoids the issue and nicely just uses Table to create dataset from a
> function! (Which I guess is a hint that this is not resolved)
> Thanks,
> Luc

  • Prev by Date: Re: Announce: O'Reilly Mathematica Cookbook Published
  • Next by Date: Re: Part specification... is neither an integer nor a list of integers
  • Previous by thread: Unstructured Interpolation in V7?
  • Next by thread: Scroll content inside InputField?