[Math] 3D point projection on 2d plane using normal vector and projection vector

Question

I want to make a program for finding a 3D point projection on 2D plane. Input would be coordinates of point $P (x,y,z)$, normal vector of projection plane and projection vector / direction vector $(s_x,s_y,s_z)$ of points onto plane. I want to calculate $x,y$ coordinates of 3D point on 2D plane but don't know how to do this.

Answer 1

You haven’t provided enough information for a complete solution: you need to choose a coordinate system for the plane. You’ve also underspecified the plane itself—where is it in relation to the origin?—but that’s a minor issue. Here’s a sketch of how you might construct the complete projection: Split the problem up into a series of simpler transformations: there’s a 3-D parallel projection onto the plane, a change of coordinate system to one in which the projection plane is a coordinate plane, and a reduction to 2-D. It might also be useful to do the coordinate system transformations in several pieces.

Since this is a parallel projection, we can w.l.o.g. assume that the plane passes through the origin. If it doesn’t, the actual projection differs from this one by a translation, which can be accounted for later. For the same reason, we can assume that the origin of the plane’s coordinate system coincides with the 3-D origin. Let $\mathbf s=(s_x,s_y,s_z)^T$ be the projection direction and $\mathbf n=(n_x,n_y,n_z)^T$ be a normal to the plane.

A common way to construct the 3-D projection is to switch to a coordinate system in which two of the coordinate axes are on the plane, and the third parallels the projection direction. If $\mathtt B$ is the matrix that represents the mapping from this coordinate system to the standard one (with the projection direction last), the projection is given by $\mathtt P = \mathtt B \operatorname{diag}(1,1,0) \, \mathtt B^{-1}$. This involves an unavoidable matrix inversion, but there’s another way to construct $\mathtt P$. Instead, switch to homogeneous coordinates and use a general formula for a central projection: if $\mathbf V$ is the homogeneous coordinate vector of the viewpoint and $\mathbf P$ the homogeneous vector that represents the plane, then the projection is $$(\mathbf V^T\mathbf P)\,I_4-\mathbf V\mathbf P^T.$$ For your projection, $\mathbf P = [n_x:n_y:n_z:d]$ for some $d$ that’s a function of the distance of the plane from the origin, and $\mathbf V = [s_x:s_y:s_z:0]$, the point at infinity that represents the projection direction. Plugging these values into the above formula results in $$\begin{bmatrix} \mathbf V^T\mathbf P-n_x s_x & -n_y s_x & -n_z s_x & -d s_x \\ -n_x s_y & \mathbf V^T\mathbf P-n_y s_y & -n_z s_y & -d s_y \\ -n_x s_z & -n_y s_z & \mathbf V^T\mathbf P-n_z s_z & -d s_z \\ 0&0&0& \mathbf V^T\mathbf P \end{bmatrix}.$$ We’re assuming that the plane passes through the origin so that $d=0$, which allows this matrix to be reduced to the $3\times3$ matrix $$\mathtt P = I_3 - {\mathbf s \mathbf n^T \over \mathbf s^T \mathbf n} = \frac1{\mathbf s^T \mathbf n} \begin{bmatrix} \mathbf s^T \mathbf n - n_x s_x & -n_y s_x & -n_z s_x \\ -n_x s_y & \mathbf s^T \mathbf n - n_y s_y & -n_z s_y \\ -n_x s_z & -n_y s_z & \mathbf s^T \mathbf n - n_z s_z \end{bmatrix}.$$ The quantity $\mathbf s^T \mathbf n$ is the dot product of the plane normal and the projection direction vector. Note that it’s not necessary for either of these to be unit vectors. Normalization is taken care of by the division by this dot product. You don’t really have to construct a matrix for this projection, though, because $$\mathtt P\mathbf v = \mathbf v - {\mathbf n^T \mathbf v \over \mathbf n^T \mathbf s}\mathbf s,$$ so you really only need to compute a pair of dot products, multiply a vector by a scalar and add two vectors. Moreover, if you’re transforming many points, one of these dot product only needs to be computed once.

The next step is to change to a coordinate system in which the first two axes lie on the plane. If you choose a pair of orthogonal unit vectors $\mathbf u_1$ and $\mathbf u_2$ to represent the desired directions of the 2-D $x$- and $y$-axes, then this transformation is performed by the orthogonal matrix $$\mathtt R = \begin{bmatrix} \mathbf u_1^T \\ \mathbf u_2^T \\ \mathbf n^T/\|\mathbf n\| \end{bmatrix},$$ i.e., the matrix with these three vectors as its rows. There’s no need to actually perform this normalization of the last row of $\mathtt R$ since you’re not going to use it, anyway. $\mathtt R$ might not be a rotation because your choice of the two unit vectors could form a left-handed basis of $\mathbb R^3$, but it’s guaranteed to be orthogonal.

You can now shift to 2-D by dropping the last coordinate, which can be represented by dropping the last row of $\mathtt R$. The complete 3-D to 2-D transformation is then the product $$\mathtt M = \begin{bmatrix} \mathbf u_1^T \\ \mathbf u_2^T \end{bmatrix} \mathtt P.$$ The resulting coordinates of the 3-D point $\mathbf v$ are thus the dot products of $\mathtt P\mathbf v$ with $\mathbf u_1$ and $\mathbf u_2$, respectively: $$\mathbf u_i^T\mathtt P\mathbf v = \mathbf u_i^T\mathbf v - {\mathbf n^T \mathbf v \over \mathbf n^T \mathbf s} \mathbf u_i^T \mathbf s = \left(\mathbf u_i - {\mathbf u_i^T \mathbf s \over \mathbf n^T \mathbf s}\mathbf n \right)^T \mathbf v.\tag{*}$$ You can, of course, assemble the parenthesized quantities into the rows of a $3\times2$ transformation matrix, if needed. For an orthogonal projection, $\mathbf s$ is parallel to $\mathbf n$, in which case this expression reduces to the dot product of $\mathbf v$ with a basis vector, as one would expect.

You might also need to make a final translational adjustment if the 2-D origin isn’t the image of the original 3-D origin.

For some suggestions on how to choose an orthonormal coordinate system for the plane, see, for example, this question, this one, this answer (for a simple method that uses a reflection instead of a rotation) and related questions.

---

As a final note, if you use the first method I mentioned above to construct $\mathtt P = \mathtt B\operatorname{diag}(1,1,0)\,\mathtt B^{-1}$ and choose for $\mathtt B$ the matrix $\begin{bmatrix}\mathbf u_1 & \mathbf u_2 & \mathbf s\end{bmatrix}$, with the $\mathbf u_i$ as above, you’ll find that the final 3-D to 2-D transformation matrix is just the first two rows of $\mathtt B^{-1}$. The expression in (*) gives you an explicit formula for them. I’ll leave it to you to determine whether computing the transformation with the explicit formula or inverting $\mathtt B$ is more efficient for your application.