Interpolation is the act of making an educated guess based on other known information. In the contents of this article, I'm going to specifically be looking at the derivation of a process to make interpolations based on a grid of 3D points. So eventually, given an input $X$ and $Y$ coordinate, we can find the interpolated Z height. This method is implemented in the Voltfolio G-Code Leveler which is mostly made to aid in PCB Leveling.

So we can start off by assuming a 2D array of points evenly spaced by $\Delta{X}$ and $\Delta{Y}$. The data at each point in the array corresponds to a measured $Z$ height. Four points $P$, $Q$, $R_I$ and $R_{II}$ define a rectangular region, that is off-set from the reference origin (IE if $P$ is $(0,0,Z)$ Then $X_{offset}$ and $Y_{offset}$ will both be 0.).

First, we have to confirm that our given point falls in the region we are looking at. Sure this seems pretty trivial, but if implemented in a program, this will have to be done. Assume a 2D arbitrary point $(X,Y)$. In order for this point to fall in the region of concern, the following inequalities must be true:

$$X_{offset} \leq X \leq (X_{offset}+\Delta{X})\\

Y_{offset} \leq Y \leq (Y_{offset}+\Delta{Y})$$

We now know that our point is in the rectangular region, however we must now find out if the point is in sub region I or II. This is important because a plane can only be constrained by 3 points. The interpolation will always uses points $P$ and $Q$, but we must determine if $R_I$ or $R_{II}$ is to be used. The point is located in sub region I if the following inequality is true. This will indicate that $R_I$ must be used for the interpolation.

$$\frac{X}{\Delta{X}} > \frac{Y}{\Delta{Y}}$$

If the above inequality is not satisfied, then the point is in sub region II, and satisfies the following relation. $R_{II}$ must be used.

$$\frac{X}{\Delta{X}} \leq \frac{Y}{\Delta{Y}}$$

Now that we have located our point, we can go on to determining Z. Let us define the following:

$$P=(X_1,Y_1,Z_1)\\

Q=(X_2,Y_2,Z_2)\\

R=(X_3,Y_3,Z_3)$$

From these three points, two vectors can be defined:

$$\vec{PQ} = <X_2-X_1, Y_2-Y_1,Z_2-Z_1>\\

\vec{PR} = <X_3-X_1, Y_3-Y_1,Z_3-Z_1>$$

Now it is possible to find the the normal vector formed by $\vec{PQ}$ and $\vec{PR}$. This can be done by carrying out the cross product:

$$\vec{n} = \vec{PQ} \times \vec{PR} =\begin{vmatrix}

\vec{i} & \vec{j} & \vec{k} \\

X_2-X_1 & Y_2-Y_1 & Z_2-Z_1 \\

X_3-X_1 & Y_3-Y_1 & Z_3-Z_1

\end{vmatrix}$$

The value of this determinate is found to be:

$$\vec{n} = (Y_2-Y_1)(Z_3-Z_1)\vec{i} + (Z_2-Z_1)(X_3-X_1)\vec{j} + (X_2-X_1)(Y_3-Y_1)\vec{k} - (Y_2-Y_1)(X_3-X_1)\vec{k} - (Z_2-Z_1)(Y_3-Y_1)\vec{i} - (X_2-X_1)(Z_3-Z_1)\vec{j}$$

It can keep things much cleaner from here on out if we develop a short hand for the $\vec{i}$, $\vec{j}$ and $\vec{k}$ parts of the above equation:

$$L = [(Y_2-Y_1)(Z_3-Z_1)-(Z_2-Z_1)(Y_3-Y_1)]\vec{i}\\

M = [(Z_2-Z_1)(X_3-X_1)-(X_2-X_1)(Z_3-Z_1)]\vec{j}\\

N = [(X_2-X_1)(Y_3-Y_1)-(Y_2-Y_1)(X_3-X_1)]\vec{i}$$

A 3D plane can now be defined given the normal vector and a point on the plane, which is assumed to be our point $(X,Y,Z)$ where $X$ and $Y$ are known values:

$$L(X-X_1)+M(Y-Y_1)+N(Z-Z_1)=0$$

Finally, solving for $Z$, we find our solution:

$$Z=\frac{-L(X-X_1)-M(Y-Y_1)}{N} - Z_1$$

Phew. Now that wasn't so bad. It's only a little bit of Calculus III, but nothing too unmanageable (or in my case, not too un-recallable). This general equation can now be used for many different interpolation purposes. Specifically for the PCB leveling problem, it allows a small amount of points to be measured in an array. Then when you need to know the $Z$ height at a point somewhere in between two recorded values, this equation can be applied and a linear extrapolated value can be found.