From: watson@madvax (Dave Watson)
Newsgroups: sci.math,sci.math.num-analysis,sci.stat.math
Subject: Re: N dimensional interpolation of ungridded data
Date: 1 Apr 1995 05:28:33 GMT

Yair Lenga (yair@cadence.com) wrote:
: I'm looking for N-Dimensional interpolation algorithm for ungridded
: data. References to public domain source code is my first preference,
: but I will be happy if someone can provide references to articles,
: books, ...

: I am familiar with the family of interpolation functions
: based on weighted-average methods, in which the weight of every
: component is based on the L2 (Or L3, L4, ...) of the distance between
: the interpolated point and the measurement point.

: Unfortunantly, the surfaces generated for my data, based on those
: methods are incorrect - they contain two many "mountains" and "vallies" -
: probably due to the fact that my measurement space is not very dense.

Satisfactory n-dimensional interpolation is quite difficult because
the conventional methods that are used for bivariate interpolation 
mostly have flaws that produce surficial artifacts and these flaws are
accentuated in higher dimensions as you have discovered in the family
of inverse distance weighted methods.

I have published a survey of interpolation and contouring methods -

CONTOURING: A guide to the analysis and display of spatial data,
Pergamon Press, 1992, ISBN 0 08 040286 0

and there are ~ 600 references to the interpolation and contouring
literature.  Direct comparisons between these methods showed that the
most robust, general, and conservative method for interpolating
scattered data is natural neighbor interpolation.  This is a weighted
average method but the weights are area-based (volume-based, etc.).

The interpolation subset (unlike all other methods) is the data which 
are natural neighbors of the interpolation point and this compensates
for the anisotropy and variation in data density which usually is a big
problem for scattered and irregularly spaced data.  With this method, the
weight for a given datum goes to zero exactly as that datum is dropped
from the interpolation subset.  The resulting surface is analogous to a
2D taut rubber sheet stretched to meet all the data.  However gradient
estimates may be used to a user-defined extent to give a more generous
interpretation of the data; the surface is smooth and continuous everywhere.

This method has been successfully applied up to 6D but is not practical
for dimensions greater than that because the number of weights per
interpolation increases factorially with dimension.

This approach can also be used with regularly gridded data and for
trivariate functions it will produce iso-surfaces which are superior to
those produced by Marching Cubes and marching tetrahedra methods.  When
this method is applied to drill core data from mining exploration, the
resulting iso-surfaces provide the shape of the ore body for any
specified cut-off grade value. 

Questions and discussion by email are welcome.
--
Dave Watson                          Internet: watson@maths.uwa.edu.au
Department of Mathematics                          watson@DIALix.oz.au
The University of Western Australia               Tel: (61 9) 380 1359
Nedlands, WA 6907  Australia.                     FAX: (61 9) 380 1028