Problem description:
Given four points $P_i$ with coordinates $(x_i, y_i, z_i)$ find the $z$-value at point $C$ with known $(x_c, y_c)$ that lies within the quadrilateral formed by the $P_i$s.
I am aware that the correct way to do it is the one described here or here. I came up with what I though is an alternative solution which however yields different (false?) results. My approach was the following:
- Construct a line (any line) that passes through C and intersects line $P_1P_2$ at, say, $P_{12}$ and $P_3P_4$ at $P_{34}$.
- Perform a 1D linear interpolation along $P_1P_2$ to find the value at $P_{12}$
- Perform a 2nd 1D linear interpolation along $P_3P_4$ to find the value at $P_{34}$
- Finally, perform a 3rd 1D linear interpolation along $P_{12}P_{34}$ to find the value at $C$
- Is there an intuitive explanation as to why my approach is off?
- Why is it that from all the possible surfaces that include points $P_i$, the one resulting from the equations at the provided links is the correct one?
Before we start, it is a good idea to refresh first.
Let's define a linear interpolation function $L$, $$\bbox{ L(x , v_0, v_1) = (1 - x) v_0 + x v_1 } \tag{1a}\label{NA1a}$$ The notation may look a bit odd, but don't let it confuse you. I just wanted to write it as a function in a form where both the interpolating coordinate ($x$) and the values ($v_0$ at $x = 0$, and $v_1$ at $x = 1$) are parameters: $$\bbox{ \begin{aligned} L(0, v_0 , v_1) &= v_0 \\ L(1, v_0 , v_1) &= v_1 \\ \end{aligned} } \tag{1b}\label{NA1b}$$
In this form, bilinear interpolation is defined as $$\bbox{ B(x, y, v_{00}, v_{01}, v_{10}, v_{11} ) = L\bigr(y, L(x, v_{00}, v_{10}), L(x, v_{01}, v_{11}) \bigr) } \tag{2a}\label{NA2a}$$ where $$\bbox{ \begin{aligned} B(0, 0, v_{00}, v_{01}, v_{10}, v_{11} ) &= v_{00} \\ B(0, 1, v_{00}, v_{01}, v_{10}, v_{11} ) &= v_{01} \\ B(1, 0, v_{00}, v_{01}, v_{10}, v_{11} ) &= v_{10} \\ B(1, 1, v_{00}, v_{01}, v_{10}, v_{11} ) &= v_{11} \\ \end{aligned} } \tag{2b}\label{NA2b}$$ If we expand $\eqref{NA2a}$ and collect the coefficients, we get $$\bbox{ \begin{aligned} B(x, y, v_{00}, v_{01}, v_{10}, v_{11}) &= v_{00} \\ \; &+ (v_{10} - v_{00}) x \\ \; &+ (v_{01} - v_{00}) y \\ \; &+ (v_{11} - v_{10} - v_{01} + v_{00}) x y \\ \end{aligned} } \tag{2c}\label{NA2c}$$ The important thing to note here is that although $B$ is bilinear, it is not a plane. (As an example, play with $v_{00} = v_{01} = v_{10} = 0$, $v_{11} = 1$.)
Let's look at OP's problem.
There are three applications of bilinear interpolation here: $$\bbox{\begin{aligned} x(u, v) &= B(u, v, x_{00}, x_{01}, x_{10}, x_{11}) \\ y(u, v) &= B(u, v, y_{00}, y_{01}, y_{10}, y_{11}) \\ z(u, v) &= B(u, v, z_{00}, z_{01}, z_{10}, z_{11}) \\ \end{aligned}}$$ with $0 \le u, v \le 1$, and $z$ being the value OP is obviously interpolating between the four points, so that when $x = x_{00}$ and $y = y_{00}$, then $u = 0$ and $v = 0$ and therefore $z = z_{00}$. Similarly for the three other points.
The key thing to notice is that the interpolation is linear in $(u, v)$, not necessarily in $(x, y)$.
Consider a quadrilateral with $(x_{00}, y_{00}) = (0, 0)$, $(x_{01}, y_{01}) = (0, 1)$, $(x_{10}, y_{10}) = (0, -1)$, and $(x_{11}, y_{11}) = (1, 0)$. Now, $$\bbox{ \begin{cases} x = u v \\ y = v - u \end{cases} \iff \begin{cases} u = -\frac{y}{2} + \sqrt{\frac{y^2}{4} + x} \\ v = \frac{y}{2} + \sqrt{\frac{y^2}{4} + x} \end{cases} }$$ within this quadrilateral, and $$\bbox{ z(x,y) = z_{00} + ( z_{00} - z_{01} - z_{10} + z_{11} ) x + \frac{z_{01} - z_{10}}{2} y - (2 z_{00} - z_{01} - z_{10})\sqrt{\frac{y^2}{4} + x} }$$ See the last summand? It means that for this particular case, $z(x,y)$ is not a bilinear function, and contains that nasty $\sqrt{y^2/4 + x}$ term.
Because that quadrilateral is still a simple convex quadrilateral, we cannot assume that $z(x, y)$ is bilinear for quadrilaterals in general. (If we were to explore this further, we'd find out that $z(x,y)$ is bilinear only when $x_{11} = x_{01} + x_{10} - x_{00}$ and $y_{11} = y_{01} + y_{10} - y_{00}$, i.e. the quadrilateral is actually a rectangle.)
This, in turn, means that OP is trying to fit a proverbial square root peg into a proverbial bilinear hole. It won't work.
Bilinear interpolation is easy to understand from a geometric point.
The easiest way is to consider a simple (non-self-intersecting) quadrilateral $ABCD$, where we interpolate using bilinear coordinates $(u,v)$, $0 \le u \le 1$ and $0 \le v \le 1$.
For example, let $A = \bigr( x(0,0), y(0,0) \bigr)$, $B = \bigr( x(1,0), y(1,0) \bigr)$, $C = \bigr( x(1,1), y(1,1) \bigr)$, and $D = \bigr( x(0,1), y(0,1) \bigr )$.
To find the point $(x, y)$ that corresponds to $(u, v)$, we
Find point $E = A + v (D - A)$ on $\overline{AD}$
(I.e., split $\overline{AD}$ as $v:(1-v)$, and name that point $E$)
Find point $F = B + v ( C - B)$ on $\overline{BC}$
(I.e., split $\overline{BC}$ as $v:(1-v)$, and name that point $F$)
Find the desired point $G = E + u (F - E)$ on $\overline{EF}$
(I.e., the desired point splits $\overline{EF}$ as $u:(1-u)$)
Perhaps this was what OP was after? I'm not sure.