algorithm - spline surface interpolation -

let's have n number of points defining surface on z-axis

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

now want able approximate f(x,y). looking algorithm linear , spline approximation. example algorithms or @ least pointers great.

this vague description of approach make linear approximation.

1. determine voronoi diagram of points (for every point in plane, find nearest (x_i,y_i))
2. take dual of delaunay triangulation: connect (x_i,y_i) , (x_j,y_j) if there line segment of points (x_i,y_i) , (x_j,y_j) equidistant (and closer other pair).
3. on each triangle, find plane through 3 vertices. linear interpolation need.

---

the following implements first 2 steps in python. regularity of grid may allow speed things (it may mess 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 np     matrix = np.array( [ [p[0],p[1],p[0]**2+p[1]**2,1] p in triangle+point ] )     return is_ccw(triangle) == ( np.linalg.det(matrix) > 0 )  triangulation = set() """ triangle in delaunay triangulation if , if circumscribing circle contains no other point.  implementation o(n^4).  faster methods @ http://en.wikipedia.org/wiki/delaunay_triangulation """ triangle in itertools.combinations(points,3):     triangle_empty = true     point in points:         if point in triangle:             next         if circumcircle_contains_point(triangle,point):             triangle_empty = false             break     if triangle_empty:         triangulation.add(triangle)