Spatial Interpolation for Irregular Grid

How would I interpolate to a point P if I have four points around it such that:

Q1 = (x1,y1), Q2 = (x2,y2), Q3 = (x3,y3), Q4 = (x4,y4)

If the coordinates formed a regular 2D grid I would use a bilinear interpolation, but I don't think I can use it for irregular grid points. I have also tried to use the inverse distance weighting and found that the results aren't as accurate as I would like. A simple linear approximation is all I'm looking for. Any help or ideas would be greatly appreciated.

In bilinear interpolation, the interpolating function is of the form $$f(x,y)=a_{xy}xy+b_xx+b_yy+c,$$ which is a quadratic function of $(x,y)$. Let's consider a general quadratic instead, $$f(x,y)=a_{xx}x^2+a_{xy}xy+a_{yy}y^2+b_xx+b_yy+c.$$ The known points give four equations that $f$ must satisfy, \begin{align} f(x_1,y_1)=a_{xx}x_1^2+a_{xy}x_1y_1+a_{yy}y_1^2+b_xx_1+b_yy_1+c&=z_1, \\ f(x_2,y_2)=a_{xx}x_2^2+a_{xy}x_2y_2+a_{yy}y_2^2+b_xx_2+b_yy_2+c&=z_2, \\ f(x_3,y_3)=a_{xx}x_3^2+a_{xy}x_3y_3+a_{yy}y_3^2+b_xx_3+b_yy_3+c&=z_3, \\ f(x_4,y_4)=a_{xx}x_4^2+a_{xy}x_4y_4+a_{yy}y_4^2+b_xx_4+b_yy_4+c&=z_4, \end{align} which can be written as a linear system of equations in the unknown parameters, $$\begin{bmatrix} x_1^2 & x_1y_1 & y_1^2 & x_1 & y_1 & 1 \\ x_2^2 & x_2y_2 & y_2^2 & x_2 & y_2 & 1 \\ x_3^2 & x_3y_3 & y_3^2 & x_3 & y_3 & 1 \\ x_4^2 & x_4y_4 & y_4^2 & x_4 & y_4 & 1 \end{bmatrix} \begin{bmatrix} a_{xx} \\ a_{xy} \\ a_{yy} \\ b_x \\ b_y \\ c \end{bmatrix} = \begin{bmatrix} z_1 \\ z_2 \\ z_3 \\ z_4 \end{bmatrix},$$ or, for short, $$\mathbf X\mathbf a=\mathbf z.$$
We can't solve this yet because it's underdetermined: we have four equations but six unknowns. But of all the possible solutions to this system, we should probably pick the one that has the smallest quadratic coefficients. If the data can be fit by a purely linear model where the quadratic coefficients are zero, then that's the solution we should use. In the rectangular case, even if we do need the $xy$ term, we shouldn't involve the other two quadratic terms because they are redundant with the linear and constant terms; thus we will recover bilinear interpolation. So let's take the solution that minimizes \begin{align} e &= a_{xx}^2+a_{xy}^2+a_{yy}^2 \\ &= \mathbf a^T\mathbf E\mathbf a \end{align} where $\mathbf E$ is the diagonal matrix with entries $[1,1,1,0,0,0]$. Now by the method of Lagrange multipliers, minimizing $\mathbf a^T\mathbf E\mathbf a$ subject to $\mathbf X\mathbf a=\mathbf z$ is equivalent to solving \begin{align} \mathbf E\mathbf a + \mathbf X^T\mathbf\lambda &= \mathbf 0, \\ \mathbf X\mathbf a &= \mathbf z, \end{align} where $\lambda$ is a four-element vector of Lagrange multipliers. Equivalently, $$\begin{bmatrix}\mathbf E & \mathbf X^T \\ \mathbf X & \mathbf 0\end{bmatrix} \begin{bmatrix}\mathbf a \\ \mathbf\lambda\end{bmatrix} = \begin{bmatrix}\mathbf 0 \\ \mathbf z\end{bmatrix},$$ which is a linear system that you can solve as usual.
Anyway, here's an animation where two points have $z=0$ and two have $z=1$. They initially form a square, and then the top edge rotates about its midpoint.
• @Rahul This derivation worked EXTREMELY well for 2D! Now I'm wondering if this is possible to expand to 3D. For example, could the general equation be expanded to $f(x,y,z) = ax^2+by^2+cz^2+dxy+exz+fyz+gx+hy+iz+j$ and E redefined as the diagonal with entries [1,1,1,1,1,1,0,0,0,0]? X would then become dependent on the number of points being used. So for something like interpolating inside of a cubic structure with 8 points, X would be an [8 x 10] matrix – Travis Feb 4 '16 at 0:16