Linear, Bi-linear or better I have been writing some code to do some interpolation of 2D data on an irregular grid. So far what I have done is:


*

*Triangulate the known points using Delaunay.

*Find the vertices of the triangles which the points to interpolate lie with-in.

*Use the Barycentric transformation to interpolate the points.  
Visually:

The black lines - convex hull.
The cyan solid lines - delaunay triangulation. The dotted cyan lines connect the points to interpolate with the vertices of the triangles they lie with-in. 
So the question is how would I go about doing bi-linear interpolation instead, or maybe there is a better approach?   
Please note I am not a mathematician so forgive me if I have use any terminology incorrectly.   
 A: You could use 
Cubic Hermite Splines to interpolate the values at the vertices. The interpolated points will be linear in the 1st order (connect along the edges), but not in the 2nd order since the surface through each triangle's vertices is underconstrained.
Wikipedia gives a 2d case. For the 3d case you need to compute the weighted sum of 10 control points. The weight is given by the Bernstein function, generalized for a barycentric coordinate zyx:
Babc = (6 * z^a * x^b * y^c) / (a!*b!*c!).
for x+y+z=1 and a+b+c = 3. 
So for instance at point A(z=1) so B300 = 6* 1^3 * 0^0 * 0^0 / 6 = 1. So trivially at the control point of vertex A, the value is 1* the value of A.
Generate 6 additional control points, by using the normal at the vertices to compute two slopes along each edge. The normal can be computed from all neighbours of every vertex. These control points (for instance A+dA/3, where dA is the slope at A along the edge to B) are multiplied by Bernstein function B210 and so on, analogous to the 2d case on wiki.
Finally you need a middle control point multiplied by B111. This is the parameter that is underconstrained. It's like a tetrahedron added to the spline, the apex of which is the control point. If it's the average of A,B,C then the tetrahedron will be flat. This control point won't influence the value at the edges, but it will influence the slope, so if you want the result to be 'smooth' you'll have to figure out this control point relative to the neighbouring triangles.
