Interpolation over triangles to get a linear approximation in two dimensions conflicts with the rectangular layout over a grid characteristic of most data sets, although the problem is readily resolved by splitting the rectangles into a pair of triangles.
Another approach to interpolating rectangular data is to use the cartesian product of coordinatewise linear interpolation. It has the further advantage that the corners of the coordinate rectangles can be the four interpolation points which the method requires.
Consider the Vandermonde determinant
The weight matrix is a Kronecker product
This time the interpolation formula
z = YT M-1 X uses the function values at the corners of the rectangle,
|A - B x + C y - D x y + E z = 0||(40)|
The pole in this formula is nice because it is capable of accomodating a saddle point in the contours, although the formula will still have problems with maxima and minima.