spline surface interpolation

Question

Let's say I have n number of points defining a surface on the z-axis

f(x1,y1) = 10
f(x2,y2) = 12
f(x3,y3) = 5
f(x4,y4) = 2
...
f(xn,yn) = 21

now I want to be able to approximate f(x,y). I am looking for an algorithm for a linear and especially a spline approximation. An example algorithms or at least some pointers would be great.

Answer 1

This is a vague description of an approach to make a linear approximation.

1. Determine the Voronoi diagram of your points (for every point in the plane, find the nearest (x_i,y_i))
2. Take the dual of this to get the Delaunay triangulation: connect (x_i,y_i) and (x_j,y_j) if there is a line segment of points so that (x_i,y_i) and (x_j,y_j) are equidistant (and closer than any other pair).
3. On each triangle, find the plane through the three vertices. This is the linear interpolation you need.

---

The following implements the first two steps in Python. The regularity of your grid may allow you to speed things up (it may also mess up the triangulation).

import itertools

""" Based on http://stackoverflow.com/a/1165943/2336725 """
def is_ccw(tri):
    return ( ( tri[1][0]-tri[0][0] )*( tri[1][1]+tri[0][1] )
            + ( tri[2][0]-tri[1][0] )*( tri[2][1]+tri[1][1] )
            + ( tri[0][0]-tri[2][0] )*( tri[0][1]+tri[2][1] ) ) < 0

def circumcircle_contains_point(triangle,point):
    import numpy as np
    matrix = np.array( [ [p[0],p[1],p[0]**2+p[1]**2,1] for p in triangle+point ] )
    return is_ccw(triangle) == ( np.linalg.det(matrix) > 0 )

triangulation = set()
"""
A triangle is in the Delaunay triangulation if and only if its circumscribing
circle contains no other point.  This implementation is O(n^4).  Faster methods
are at http://en.wikipedia.org/wiki/Delaunay_triangulation
"""
for triangle in itertools.combinations(points,3):
    triangle_empty = True
    for point in points:
        if point in triangle:
            next
        if circumcircle_contains_point(triangle,point):
            triangle_empty = False
            break
    if triangle_empty:
        triangulation.add(triangle)

Answer 2

Interpolation on irregular 2D data is not that easy. I know of no true spline generalization to irregular 2D.

Besides the triangulation-based approaches, you can have a look at Barnes (http://en.wikipedia.org/wiki/Barnes_interpolation) and Inverse Distance Weighting (http://en.wikipedia.org/wiki/Inverse_distance_weighting), or more generally, RBF (http://en.wikipedia.org/wiki/Radial_basis_functions).

If your point are strongly non-uniformly spread (dense clusters), it can be necessary to make the size of the functions adaptive, or resort to approximation rather than interpolation.

Answer 3

You can use your points as the control points of a Bezier (or Bspline) surface, especially if (xi, yi) sample a rectangle in the XY plane. In this respect, there's no fitting involved.

The surface you will get will be in the convex hull of your points, and will intersect (interpolate) the points at the boundary of {xi, yi}.

If you'd like to experiment, This forums posting seems to contain simple code in Matlab, and you can use GuIRIT to do the same if you don't have Matlab (though it requires figuring out the file format of the program).