Miroslav Dressler
The article describes a new method for approximation / interpolation of
irregularly distributed points in 3D Euclidean space by a continuous
function of two independent variables, which can be used for surface map
creation in many fields such as geoscience, geography, geology, meteorology
and so on. The method was named ABOS (Approximation / interpolation
Based On Smoothing) and uses very simple mathematical
tools - numerical tensioning and smoothing. The tensioning and smoothing
are performed so that elements of matrix, which represents surface z-values
at nodes of a regular rectangular grid, are repeatedly replaced by the weighted
average of selected surrounding elements. The selection of elements involved
into weighted average depends on the type of tensioning or smoothing.
Despite the fact that the mathematics of the ABOS method is simple,
the resulting surface can be modified by a few parameters so that it is
comparable with the surface created by sophisticated methods such
as Kriging, Radial Basis Functions or Minimum
Curvature.
There are several methods commonly used for approximation / interpolation
of irregularly distributed points in 3D Euclidean
space R3. Probably the
most advanced are the Kriging and Radial Basis Functions methods,
which provide excellent results and which are suitable for using in almost
any field. However, there may be a problem with using these methods - they
need to solve a system of linear equations having
the number of equations equal to the number of points. If the number
of input points is large (say more than a thousand), the number of equations
is also large and cannot be solved in an acceptable time - in such case,
a search algorithm must be applied (to reduce the number of equations) and
additional parameters influencing the quality of the resulting surface must
be specified - search radius, number of search sectors, maximal number of
points in each sector and so on.
The ABOS method described in this article does not need to solve any
system of equations and does not use any search algorithm even if the number
of input points is very great - say millions.
In contrast to sophisticated interpolation methods, the ABOS method
does not fulfil any mathematical formulation
like linear combination of basis functions, statistical formulation of the
best linear unbiased estimate or requirement of minimal curvature of resulting
surface. Instead, it provides tools for modeling of resulting surface based
on numerical tensioning and smoothing enabling to create broad range of surface
shapes and to accommodate resulting surface according to user's conception.
Let
XYZ={(Xi,Yi,Zi)ε
R3;
i=1,...,n} is a given set of irregularly
distributed points in 3D Euclidean space and
XY={(Xi,Yi)ε
R2;
i=1,...,n} is their horizontal projection.
In this article, the solving an interpolation or approximation problem
will mean the finding such continuous function of two independent variables
f(x,y), for which
f(Xi,Yi)=Zi
or
|f(Xi,Yi)-Zi|<d respectively.
The interpolation (approximation) function used by the ABOS method
is defined by a matrix of real numbers
P, which represents the surface
at the nodes of a square grid covering
points XY. The domain of the interpolation
function is defined by the outer lines of the grid - see the next figure.
Figure 1: Square grid for interpolation.
The value of the surface at any point
(x0,y0)
within the grid can be evaluated from the equation of the bilinear
polynomial
f(x,y)
=
axy
+
bx
+
cy
+
d
which is defined by the corner points of the grid square containing the point
(x0,y0).
The following notation is used in the next text:
| Zmin | minimum of z-coordinates of points XYZ |
| Zmax | maximum of z-coordinates of points XYZ |
| i1,j1 | size of the grid - number of columns and rows of the matrix P |
| Pi,j | elements of the matrix P, i=1,...,i1 j=1,...,j1 |
| DP | auxiliary matrix with the same size as the matrix P |
| Z | vector of z-coordinates of points XYZ |
| DZ | auxiliary vector with the same size as the vector Z |
| K | matrix of integer numbers, which has the same size as the matrix P and determines, for any node in the grid, the Chebyshev distance between the node and the nearest point XYZ measured at steps of grid (so called matrix of the distances) |
| Kmax | maximal element of the matrix K |
| NB | matrix of integer numbers, which has the same size as the matrix P and determines, for any node in the grid, the index of the nearest point XYZ (so called matrix of the nearest points) |
Before any type of tensioning or smoothing is applied, the matrices NB and K are computed using the algorithm based on the "circulation" around the points XYZ as the following figure indicates:
Figure 2: Computation of the matrices
NB and
K
All elements of the matrices NB and K are initially set to zero and the process of circulation continues as long as there are zero values in the matrix NB. The Euclidean distance is compared only if the element Ki,j corresponding to the evaluated node is not zero and IC/SQRT(2)<= Ki,j, where IC is the number of the current circulation. By this way, the number of Euclidean distance computations is significantly reduced.
The computation of the matrix NB defines the natural division of the interpolation function domain into polygons (so called Voronoi or Thiessen polygons, see figure 3). Inside each polygon the values of the matrix NB are equal to the ordinal index of the corresponding point XYZ.
Figure 3: Division of the domain of the interpolation function
The algorithm of the ABOS method can be described by the following scheme:
| step | description |
1 |
Determination of grid size, computation of the matrices NB and K, Z -> DZ |
2 |
Per partes constant interpolation of DZ values into the matrix P |
3 |
Tensioning and smoothing of the matrix P |
4 |
Determination of differences DZ-P -> DZ |
5 |
If the maximal difference (maximal value of elements from the vector DZ) does not exceed defined precision, algorithm is finished |
6 |
P -> DP, per partes constant interpolation of DZ values into the matrix P |
7 |
Tensioning and smoothing of the matrix P |
8 |
P+DP -> P |
9 |
Continue from step 4 again. |
Particular steps of the algorithm are explained in more detail:
Step 1. In the first step, the value of
i1 (grid size in x-direction)
and j1 (grid size in y-direction)
has to be specified. Despite the fact that the grid is referenced as a "square",
in general it is not really possible to set a square grid - the
grid is generally formed by rectangles. The algorithm calculates (according
to the value of i1) the value of
j1 so that the difference between
sizes of rectangle sides is as small as possible.
In this step, the matrices
NB and
K are calculated and the
values of the vector Z are
copied into the vector
DZ.
Step 2. Interpolation of DZ values into the matrix P is performed such that the value Pi,j for every node of grid is equal to the value DZ of the nearest point XYZ (this means, Pi,j = DZNBi,j). Such interpolation (per partes constant, see figure 4) is very fast.
Figure 4: Per partes constant interpolation.
Step 3. Repeated tensioning (see figure 5), linear tensioning (see figure 6) and smoothing (see figure 7) of the matrix P improves the shape of the resulting surface, but smoothing decreases the accuracy - the surface does not pass through the z coordinates of points XYZ exactly.
Figure 5: Tensioning of the matrix
P.
Figure 6: Linear tensioning of the matrix
P.
Figure 7: Smoothing of the matrix
P.
Tensioning replaces elements of the matrix
P by the average value:
Pi,j = (Pi+k,j + Pi,j+k + Pi-k,j + Pi,j-k)/4, k = Ki,j
only of k > 0. The following scheme shows the nodes (marked by grey circles) corresponding to the elements of the matrix P, which are involved to the tensioning.
Tensioning is repeatedly performed in the loop with this pattern:
DO N
=
MAX(4,Kmax/2+2),1,-1
...
ENDDO
If k is greater than the decreasing
loop variable N, then
k
=
N.
Linear tensioning modifies the matrix P according to the formula:
Pi,j = (Q(Pi+u,j+v + Pi-u,j-v)+ Pi-v,j+u+ Pi+v,j-u)/(2Q+2), where
(u,v) is the vector from the node
i,j to the nearest grid node
of point
NBi,j and
Q
=
L(Kmax-Ki,j)2.
The constant L
=
1/((0.107·Kmax-0.714)·Kmax)
is an empirical constant suppressing the influence of
Kmax . This formula is valid
for default degree of linear tensioning 1. Generalized formula for other
degrees is in section Degrees of linear tensioning.
The following scheme shows the nodes (marked by grey circles) corresponding
to the elements of the matrix
P, which are involved to the linear tensioning.
Linear tensioning is repeatedly performed in the loop with this pattern:
DO N
=
MAX(4,Kmax/2+2),1,-1
...
ENDDO
If the
length |(u,v)| of the vector
(u,v) is greater than the decreasing
loop variable N, then the vector
(u,v) is multiplied by the constant
c so that
c
·|(u,v)|
=
N.
Smoothing replaces elements of the matrix P by the weighted average value:
Pi,j
=
(∑
Pk,l+
Pi,j(q
ti,j-1))/(8
+
q
ti,j),
k=i-1,..,i+1,
l=j-1,..,j+1,
where q is the user-defined parameter controlling smoothness of the interpolation (its default value is 0.5) and ti,j are weights, which are zero before the first smoothing and afterwards they are computed according to the formula
ti,j
=
(∑(Pi,j-Pk,l))2,
k=i-2...i+2, l=j-2...j+2
and re-calculated into interval [0,100]. In brief it can be said, the values of ti,j are higher at nodes, where the surface has a local extreme and lower at nodes, where the surface is decreasing / increasing. The following scheme shows the nodes (marked by grey circles) corresponding to the elements of the matrix P, which are involved to the smoothing.
Smoothing is repeatedly performed in the loop with this pattern:
DO N
=
MAX(4,Kmax*Kmax/16)1,-1
···
ENDDO
As an option the ABOS method enables to perform so called LES smoothing - in this case the formula for smoothing is not applied if the decreasing loop variable N is greater than Ki,j+1. This modification of smoothing suppresses oscillations and exceeding of local extremes, if they occur.
Step 4. The tensioned and smoothed surface does not pass through the Z-coordinates of points XYZ exactly, so the differences DZi = Zi-f(Xi,Yi), i=1...n are calculated. The Z-value of the surface can be computed at any point (x0,y0) from the bilinear polynomial equation:
f(x,y) = axy + bx + cy + d,
where a,
b,
c and
d are defined by the corner
points of the grid square containing the point
(x0,y0).
At the beginning of this step, the algorithm offers the option to modify
elements of the matrix P by
the transformation
Pi,j = A·Pi,j + B,
where constants A and
B minimize the term
∑(A·f(Xi,Yi)
+
B
-
DZi)2.
By this way, the number of iteration steps can be reduced by 20-50%.
Step 5. If the maximal difference (maximal value of elements of the
vector DZ ) is less than the specified
accuracy of the interpolation / approximation, the algorithm of the
ABOS method finishes. In the opposite case the algorithm continues
from step 6.
The accuracy is specified as a percentage value from the difference
Zmax
-
Zmin.
Step 6. The matrix P is copied into the matrix DP (element by element) and the vector of differences DZ is interpolated into the matrix P, likewise as in step2.
Step 7. The surface of differences represented by the matrix P is tensioned and smoothed, likewise as in step3.
Step 8. The matrices P and DP are added element by element (sum of two smooth surfaces is a smooth surface) and the result is stored into the matrix P. The algorithm continues from step 4.
Example of whole interpolating process can be viewed here as an animated
GIF - to start the animation click on the next picture or refresh your
browser (shortcut key F5 in Microsoft Internet Explorer).
.
Degrees of linear tensioning
There are four degrees of linear tensioning (0-3) implemented in the ABOS method. The formula for linear tensioning can be expressed in this generalized form:
Pi,j
=
(Q(Pi+u,j+v
+
Pi-u,j-v)+
R(Pi-v,j+u+
Pi+v,j-u))/(2Q+2R)
for all i=1,..,i1, j=1,..,j1;
Ki,j>0
where the weights Q and
R, depending on the degree of linear
tensioning, are calculated as follows:
| Degree | Q |
R | L |
0 |
L(Kmax-Ki,j)2 | 1 |
0.7/((0.107Kmax-0.714)Kmax) |
1 |
L(Kmax-Ki,j)2 | 1 |
1.0/((0.107Kmax-0.714)Kmax) |
2 |
L(Kmax-Ki,j) | 1 |
1.0/(0.0360625Kmax+0.192)) |
3 |
1 |
0 |
- |
Formulas for the computation of the constant L are empirical and their role is to suppress the influence of Kmax .The next figure contains a cross-section plot demonstrating the typical influence of the linear tensioning degree.
Figure 8: Influence of linear tensioning degree on surface shape.
Most of interpolation methods offer some possibility, how to modify constructed
surface. For example, the Radial Basis Functions method enables to
vary smoothness of surface using the smoothing parameter, the Kriging
method may use different variogram models with different parameters and the
Minimum Curvature method uses the tension parameter.
The ABOS method can modify the created surface namely through the
use of smoothness parameter q and
an approximation can be achieved by setting an appropriate accuracy parameter.
Moreover, the number of smoothing cycles suggested by the computer implementation
of the ABOS method can be increased so that suitable trend surface
is obtained. The cross-section through seven surfaces in the next figure
illustrates possible modifications of the surface shape generated by the
ABOS method.
Figure 9: Flexibility of the ABOS method.
In contrast to other interpolation methods, the ABOS method only requires two points with different z-values for creating a surface. The next picture shows the surface created from two points (0,0,0) and (1,1,1):
Picture 8: A surface created from two points.
Depending on the configuration of the points XYZ, a smooth interpolation may cause unwanted oscillations in the generated surface, which is a common problem of most smooth interpolation methods. In this paragraph we will examine oscillations on specially designed example (see figure 9) containing 13 points distributed along both diagonals of the square. All points have z-coordinate equal to zero except the point in the centre of the square, where the z-coordinate is one. We can expect that there will be oscillations between points having zero z-coordinates - that is why we will focus our attention to the cross-section going through points L09, L05 and L01.
Figure 9: Theoretical example of zero based data
Let's now examine details in the cross-section going through points L09, L05 and L01 where the oscillations were expected. As for size of oscillation between points L09 and L05, the Kriging method provides the best result while the Minimum Curvature method provides the worst result. As for oscillation between points L05 and L01, both named methods still produce oscillations but in the opposite order - the oscillation size produced by the Minimum Curvature method is slightly smaller than in the case of Kriging method.
The ABOS method produce worse result than the Kriging and better result than Minimum Curvature method between points L09 and L05, but oscillations between points L05 and L01 are negligible. Moreover, if the LES smoothing is used during the interpolation process (see paragraph Smoothing), the suppression of oscillations and improper extremes is very effective.
Figure 10: Comparison of oscillations.
In any case, smooth interpolation may produce unwanted oscillations and improper extremes, which for example means that there are areas containing negative values in the solution of an interpolation problem, while only positive or zero values are possible. To eliminate this problem, the implementation of the ABOS method contains the possibility to replace the grid node values below the specified constant (zero in our case) by this constant.
The above mentioned smoothness parameter enables to set the shape of a generated
surface in the surroundings of local extremes. The higher the smoothness
parameter, the "sharper" the local extremes will be created. Let us use the
previous example to show different shapes of a generated surface at the point
with the z-coordinate equal to 1.
If the Kriging method is used with the linear variogram and zero nugget
effect, the resulting surface is sharp. The next picture of the cross-section
going through point L13, illustrates that the ABOS method (blue) is
able to produce a sharp surface similar to the surface created by the
Kriging (red) method. If the smoothness parameter is small, for example
0.1, the ABOS method produces a very smooth surface (green) like the
Minimum Curvature method with a tension of 0.2 (grey).
Picture 11: Comparison of the Kriging, Minimum
Curvature and ABOS methods.
If the input data represents measurement from some terrain, it is usually suitable to use only linear tensioning with small number of smoothing cycles (which is close to the Triangulation with Linear Interpolation method) - see the digital model of a stone quarry in figure 12.
Figure 12: Digital model of a stone quarry.
As a rule, characteristic points of terrain are measured - this means that the person performing such a measurement surveys only points where the slope of terrain changes (tops, edges, valleys and so on). For interpretation of such data, the Triangulation with Linear Interpolation method is usually used, so the ABOS method is applicable as well.
The last example shows surfaces created from geological data, which contains
9500 points. The data was interpolated using the Kriging and
ABOS methods to illustrate differences between them. The
Kriging method needed 115 seconds for interpolating these points into
a grid of 200x400 nodes, while the ABOS method needed only 10 seconds
(time was, of course, measured on the same computer).
The result of both interpolations is shown in the next figure. While the
surface in regions with no data is smooth in the ABOS case,
Kriging created, in these regions, strange jumps resulting from
the search algorithm. The Radial Basis Functions method has the same
problem.
Figure 13: Interpolation of geological data using the Kriging and
ABOS methods.
The presented method is a simple, fast and flexible interpolation / approximation method applicable in many fields and can provide results similar to other methods.
Among the mentioned features, the implementation of the method enables:
The ABOS method is implemented in the software package SurGe,
which is available in the following WWW page:
http://surgeweb.mzf.cz/index.htm
Detailed aspects of surface interpolation and approximation is desribed in Art of Surface Interpolation (online eBook in pdf format: http://surgeweb.mzf.cz/AOSIM/AOSIM.pdf)
Author (M.Dressler@centrum.cz) will be grateful for any remarks, opinions or suggestions.