Spline surface interpolation

Question

Say that I have n the 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 linear and especially spline approximation. Sample algorithms, or at least some pointers, would be great.

Answer 1

This is an indefinite description of the linear approximation approach.

• Define the Voronoi diagram of your points (for each point on the plane, find the closest (x_i,y_i) )
• Take a double value to get the Delaunay triangulation : connect (x_i,y_i) and (x_j,y_j) if there is a segment of the dot line, so that (x_i,y_i) and (x_j,y_j) are equidistant (and closer than any other couple).
• On each triangle, find a plane through three vertices. This is the linear interpolation you need.

---

The first two steps in Python are listed below. The regularity of your grid can allow you to speed up the process (it can also ruin 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

Interpolating to irregular 2D data is not so simple. I do not know the true generalization of spline to irregular 2D.

Besides triangulation-based approaches, you can take 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 is very unevenly distributed (dense clusters), you may need to adapt the size of the functions or use an approximation rather than interpolation.

Answer 3

You can use your points as control points on a Bezier (or Bspline) surface, especially if (xi, yi) a rectangle pattern in the XY plane. There is no involvement in this regard.

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

If you want to experiment, this forum post seems to contain simple code in Matlab , and you can use GuIRIT to do the same if you don't have Matlab (although this requires finding the file format of the program).