3D Linear Interpolation

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.


An example 3D mesh that shows an upward sloping warp

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.).

A 3D plane defined by 4 points, and is subdivided into two regions. All points have a known or measured Z "height"
A 3D plane defined by 4 points, and is subdivided into two regions. All points have a known or measured Z "height" 

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.