Algorithm for Line-in-Plane intersection in 3D

Simplify notation a bit by calling the point on the plane $P,$ the first point on the line $Q,$ and the difference between the second point and first point on the line $v$ (so $v$ is the vector from the first point to the second point).

The equation of the plane is $n\cdot x + d = 0.$ Since $P$ is on the plane, $d$ must equal $-n\cdot P.$

The equation of the line is $Q + tv$, where $t$ is any real number. The value of $t$ where the line and plane intersect must satisfy $$ n\cdot (Q + tv) + d = 0. $$

Solving this equation for $t$ yields $$ t = \frac{-d - n\cdot Q}{n\cdot v}. $$

Use homogeneous coordinates and the Plücker matrix of the line. If the two points that define the line are $\mathbf A$ and $\mathbf B$, this matrix is $L=\mathbf A\mathbf B^T-\mathbf B\mathbf A^T$. If $\mathbf\pi$ is the vector that represents a plane, its intersection with the line is $$L\mathbf\pi = (\mathbf B^T\mathbf\pi)\mathbf A-(\mathbf A^T\mathbf\pi)\mathbf B.$$ That is, the intersection of a line defined by two points and a plane can be compute with two dot products, two scalar products, a vector addition and then three divisions for the dehomogenization. If the last coordinate of $L\mathbf\pi$ is zero, then the line is parallel to the plane; if the entire product vanishes, then the line lies on the plane.

In this case, $\mathbf A=(x_3,y_3,z_3,1)^T$, $\mathbf B=(x_4,y_4,z_4,1)^T$ and $\mathbf\pi=(x_1,y_1,z_1,-x_1x_2-y_1y_2-z_1z_2)^T$. The elements of the latter are the coefficients of the point-normal equation of the plane.