Interpolation / approximation of surfaces based on numerical tensioning and smoothing

Miroslav Dressler

Abstract

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.

Introduction

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.

Description of the ABOS method

Definition of the interpolation function and notations

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)

Matrix of the nearest points and matrix of distances

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

Interpolating algorithm

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.

Flexibility of the ABOS method

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.

Examples

Example 1: Only two points

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.

Example 2:  Smoothness of interpolation and oscillations

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.

Example 3: Smooth and sharp interpolation

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.

Example 4: Digital model of terrain

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.

Example 5: Geological data

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.

Conclusion

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:

• to create and manage projects containing lot of maps
• to filter input data
• to transform coordinates of map objects (points, faults, polylines, boundaries)
• using mouse and keyboard to digitize the map objects
• to model discontinuity in a generated surface
• to double grid using cubic polynomial interpolation between grid nodes
• to calculate and display isolines
• to compute and display gradient lines
• to display raster color map and 3D view of resulting surface including shadowed relief
• to blank grid outside the boundary
• to perform mathematical calculations with surfaces
• to define and display cross-sections through several surfaces
• using special procedures to solve special tasks such as
• wedging out of layers
• preserving extrapolation trend in areas without data
• direct conversion of seismic reflection times into the structural depth of layer boundary
• creation of layer thickness maps and volume calculation

The ABOS method is implemented in the software package SurGe, which is available in the following WWW pages:
http://surgeweb.sweb.cz/surgemain.htm
http://surge.czweb.org/surgemain.htm

Detailed aspects of surface interpolation and approximation is desribed in Art of Surface Interpolation (online eBook in pdf format: http://m.dressler.sweb.cz/AOSIM.pdf)

Author (M.Dressler@centrum.cz) will be grateful for any remarks, opinions or suggestions.