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.

The [wikipedia][1] article is a bit daunting but try at least to look at the examples section. [1]: http://en.wikipedia.org/wiki/Spline_interpolation — Mihai Maruseac, Mar 11 '12 at 16:48
For functions of the form f(x,y), it's more common to make an assumption about the form of the underlying data (polynomial of degree K, sum of N Gaussians, etc.) and then determine the coefficients by least squares. Would that work here? Do you know anything about what the data represent? If you really want a spline, then NURBS http://en.wikipedia.org/wiki/NURBS are worth a look. They have nice properties for rendering. Construct a Delaunay triangulation of the (x,y) points to get the basis unless they are on a regular grid. For the plane fitting, you'll want a standard least squares fit. — Gene, May 07 '14 at 06:04
I doesn't necessarily need to be a spline. Linear would suffice for the moment. I cannot say much about the plane fitting but the data points are on a regular grid. Just with some data points missing. — tcurdt, May 07 '14 at 09:19

Teepeemm · Accepted Answer · 2014-05-07T16:16:52.190

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

Determine the Voronoi diagram of your points (for every point in the plane, find the nearest (x_i,y_i))
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).
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)

score 2 · Answer 2 · 2014-05-06T08:31:32.343

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.

edited May 06 '14 at 08:31

answered May 06 '14 at 08:26

Realistically they are not spread that very much. Worst case I could probably even live with a linear interpolation but would prefer a smoother surface. – tcurdt May 06 '14 at 10:14
Approximation sounds like an interesting angle, too. – tcurdt May 06 '14 at 10:15
+1 for radial basis functions. I wrote a paper a few years back working with these objects generalized to functions on manifolds. They work great as long as you don't have the dense clusters and the number of points n doesn't get too big. (If n does get big, you want your RBF to have compact support so that the matrix involved is sparse. But then you'd need to use sparse linear algebra to keep it scaling acceptably.) Nice, smooth interpolation. – Raman Shah May 11 '14 at 14:19

score 1 · Answer 3 · answered Mar 11 '12 at 17:36

1

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).

answered Mar 11 '12 at 17:36

user1071136

15,636
4
42
61

The final implementation will have to be in ruby - so Matlab is not really an option. But the problem is indeed about an rectangle in the XY plane. – tcurdt Mar 11 '12 at 21:50
I never used Ruby, but I'm sure there's a Bezier/Bspline package for it. – user1071136 Mar 11 '12 at 21:58

spline surface interpolation

3 Answers3