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.
Solution 1:
"words that end in GRY" has made an excellent suggestion in the comments: there are a lot of well-studied techniques for scattered data interpolation. In general, and especially if you had an arbitrary number of points, I would recommend using thin plate splines, which work very well in my experience.
Then again, if you have exactly four points, one might wonder if there is a general interpolation scheme that exactly recovers bilinear interpolation when the points form a rectangle. Well, here's one.
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.
