The calculation and source code of the straight line intersecting the space plane

目的：in writing a graphics project,It is common to encounter situations where three-dimensional planes intersect.It is also the foundation of graphics

A space plane is the same as a straight line,There are two main expressions.
1）General algebraic expression： $ax+by+cz+d=0$
2）geometric expression(都是列向量)：$(P,N)$,where is a point on the plane,$N$为平面的法向量.Its translation into algebraic expression is $(P-O{)}^T\ast N=0$（Needs to be transposed to a row vector）
在后续的计算中,在写C++代码中,Generally inclined to geometric expression.Because it is more intuitive,It is also useful for subsequent calculations.Below I will introduce the algorithm for the intersection of two planes based on the geometric expression,Give lastC++代码.

方法一：
Suppose two planes are ：

$$(P_1,N_1)=>(X-P_1{)}^T\ast N_1=0(1)$$

$$(P_2,N_2)=>(X-P_2{)}^T\ast N_2=0(2)$$

There are three situations in which two spatial planes intersect
1：The intersection is a straight line.
2：if the planes are parallel,no intersection.
3：in the case of co-planarity,相交,But not flat.
The following describes the case where the planes intersect as straight lines：

A representation of a straight line where two planes intersect：$O_3$ (起始顶点);$N_3$(方向).
1）First calculate the direction of the line of intersection$N_3$, Because the direction of the line of intersection is found to be perpendicular to bothplane1和plane2The normal vector of the two planes.叉乘$N_1\times N_2$,得到如下：
$$N_3=N_1\times N_2(3)$$

The red is the normal vector of the obtained plane intersection.(图像中的)

2）Need to find a vertex on the line.It can be any vertex,最好的是$P_1$,$P_2$Vertices are projected onto lines$P_1^{\prime },P_2^{\prime }$的中点.为了方便讲述,This blog is only used$P_1$A point projected onto a line$P$,如下图

它是平面$(P_1,N_3)$and the intersection of the line.将$P_1,N_3$写成如下格式：

$$(P_1,N_3)=>(X-P_1{)}^T\ast N_3=0(4)$$

因此用(1)(2)(4)A system of linear equations can be formed as :

组成新的平面.Then use the intersection of three planes to get vertices.如下图示意：
The third normal vector can be expressed as ：
$$\begin{array}{c}(X-P_1{)}^T\ast N_1=0\\ (X-P_2{)}^T\ast N_2=0\\ (X-P_1{)}^T\ast N_3=0\end{array}$$

Unpacking the formula can be obtained as follows：
$$\begin{array}{c}{X}^T\ast N_1=P_1^T\ast N_1\\ X^T\ast N_2=P_2^T\ast N_2\\ X^T\ast N_3=P_1^T\ast N_3\end{array}$$

其中$X$为未知数,它是向量,The vertices on the line can be obtained by the above formula$X^{\ast }$,它的值也就是$P=X^{\ast }$.
The code for this part will be added later.It is also very simple as follows:

#include<Eigen/Eigen>
#include<iostream>
#include<vector>

/* m=[n1, n2, n3] */
Eigen::Matrix3f ConstructMatrix3fFromVectors(Eigen::Vector3f& n1, Eigen::Vector3f& n2, Eigen::Vector3f& n3)
{
    
    Eigen::Matrix3f m;
    m.block<1, 3>(0, 0) = n1;
    m.block<1, 3>(1, 0) = n2;
    m.block<1, 3>(2, 0) = n3;
    return m;
}

float Pnt3d2PlaneDist(Eigen::Vector3f& ppnt, Eigen::Vector3f& pdir, Eigen::Vector3f& pnt)
{
    
	float sd = pdir.norm();
	if (sd < 1e-8)
		return std::numeric_limits<float>::max();
	float sb = std::abs(pdir.dot(ppnt - pnt));
	float d = sb / sd;
	return d;
}

int ComputeLineFromTwoPlanes(Eigen::Vector3f& src_pnt, Eigen::Vector3f& src_dir, 
Eigen::Vector3f& tgt_pnt, Eigen::Vector3f& tgt_dir, Eigen::Vector3f& lpnt, Eigen::Vector3f& ldir)
{
    
    ldir = src_dir.cross(tgt_dir);
    //Determines whether two planes are coplanar,或者平行
    if ((std::abs(ldir[0]) + std::abs(ldir[1]) + std::abs(ldir[2])) < 1e-6)	//Determine whether two planes are parallel
	{
    
	    float dist = Pnt3d2PlaneDist(tgt_pnt, tgt_dir, src_pnt);
	    if (dist < 1e-4)
		    return 0;	//共面
	    else
		    return 1;	//平行面
    }
    //Calculate the vertices projected onto the line
    float pn1 = src_pnt.dot(src_dir);
    float pn2 = tgt_pnt.dot(tgt_dir);
    float pn3 = src_pnt.dot(ldir);
    Eigen::Vector3f b = Eigen::Vector3f(pn1, pn2, pn3);
    Eigen::Matrix3f A = ConstructMatrix3fFromVectors(src_dir, tgt_dir, ldir);
    lpnt = A.inverse()*b;
    return 2;
}

//test adjugate matrix

int main(int argc, char** argv)
{
    
    Eigen::Vector3f sor_c_pnt = Eigen::Vector3f(3.84935, 0.742589, 1.38037);
    Eigen::Vector3f sor_c_dir = Eigen::Vector3f(-0.999087, -0.0426946, -0.00180283);
    Eigen::Vector3f tgt_c_pnt = Eigen::Vector3f(2.73047, 0.188734, 1.95757);
    Eigen::Vector3f tgt_c_dir = Eigen::Vector3f(0.0486091, -0.99881, 0.00407396);
    
    Eigen::Vector3f lpnt, ldir;
    ComputeLineFromTwoPlanes(sor_c_pnt, sor_c_dir, tgt_c_pnt, tgt_c_dir, lpnt, ldir);
    std::cerr <<"lpnt, dir: " << lpnt.transpose() <<", " << ldir.transpose() << std::endl;    
    return 0;
}

运行的结果如下：

方法二：
If the expression of the function is $ax+by+cz+d=0$. The expressions for the two planes are
$$a_{1}x+b_{1}y+c_{1}z+d_1=0(5)$$

$$a_{2}x+b_{2}y+c_{2}z+d_2=0(6)$$

It can be written in the form of a vector as :

$$N_1^T\times X=-d_1(7)$$

$$N_2^T\times X=-d_2(8)$$

1：Calculate the straight line direction：使用同样的方法,可以得到$N_3=N_1\times N_2$,
2：Find a vertex on the line：Because both planes are vertices,So finding the best vertex is not considered.我们使用youtube上的一个视频,它的介绍如下：

Suppose one of the axes is 0($z=0$).Then make up the following formula.
$$\begin{array}{c}{a}_{1}x+b_{1}y+d_1=0\\ a_{2}x+b_{2}y+d_2=0\end{array}$$
The resulting vertex is $P=(0,y_s,z_s)$
但是在实际的代码中,Selected for axis$0$need to be considered.我们知道,for an axis,If it contributes almost nothing,可能会出现问题.Just if in the equation$a_{1}x+b_{1}y+c_{1}z+d_1=0(5)$,如果$c_1=10000,a_1=0.0001,b_1=100$,Then it might be possible to find computational$x$会非常大.Such an intersection exists$x$The on-axis coefficients are very small,means it follows$x$轴平行.Therefore, in the actual calculation of the intersection,We need to calculate the main axis of the plane intersection.It has normal vectors$N_3$表示.
如果$N_3$中$x$系数最大.We assume that the intersection line is $P=(0,y,z)$,Then enter the equation as :

$$\begin{array}{c}{b}_{1}y+c_{1}z+d_1=0\\ b_{2}y+c_{2}z+d_2=0\end{array}$$

如果$N_3$中$y$系数最大.We assume that the intersection line is $P=(x,0,z)$
$$\begin{array}{c}{a}_{1}x+c_{1}z+d_1=0\\ a_{2}x+c_{2}z+d_2=0\end{array}$$

如果$N_3$中$z$系数最大.We assume that the intersection line is $P=(0,y,z)$
$$\begin{array}{c}{a}_{1}y+b_{1}z+d_1=0\\ a_{2}y+b_{2}z+d_2=0\end{array}$$

具体实现代码如下：

int GetIntersectionLineFromTwoPlanesFunction(Eigen::Vector4f& sor_fun, Eigen::Vector4f& tgt_fn, 
		Eigen::Vector3f& lpnt3d, Eigen::Vector3f& ldir3d)
	{
    
		Eigen::Vector3f sorn = Eigen::Vector3f(sor_fun[0], sor_fun[1], sor_fun[2]);
		Eigen::Vector3f tgtn = Eigen::Vector3f(tgt_fn[0], tgt_fn[1], tgt_fn[2]);
		Eigen::Vector3f ln = sorn.cross(tgtn);
		ldir3d = ln;
		//Judge the three axes,Find the largest axis
		if ((std::abs(ln[0]) + std::abs(ln[1]) + std::abs(ln[2])) < 1e-6)	//Determine whether two planes are parallel
		{
    
			Eigen::Vector3f sm_pnt;
			SamplePnt3dFromPlaneFnc(sor_fun, sm_pnt);
			float dist = Pnt3d2PlaneDist(tgt_fn, sm_pnt);
			if (dist < 1e-4)
				return 0;	//共面
			else
				return 1;	//平行面
		}
		int max_cn = FindMaxCoordFromVec3D(ln);
		float d1, d2;
		d1 = sor_fun[3];
		d2 = tgt_fn[3];
		if (max_cn == 0)
		{
    
			lpnt3d[0] = 0.0f;
			lpnt3d[1] = (d2*sorn[2] - d1*tgtn[2]) / ln[0];	//ln[0]为b_1,c_1和b_2,c_2的矩阵det|A|
			lpnt3d[2] = (d1*tgtn[1] - d2*sorn[1]) / ln[0];
			return 2;
		}
		else if (max_cn == 1)
		{
    
			lpnt3d[0] = (d1*tgtn[2] - d2*sorn[2]) / ln[1]; //ln[1]为a_1,c_1和a_2,c_2的矩阵det|A|
			lpnt3d[1] = 0.0f;
			lpnt3d[2] = (d2*sorn[0] - d1*tgtn[0]) / ln[1];
			return 3;
		}
		else
		{
    
			lpnt3d[0] = (d2*sorn[1] - d1*tgtn[1]) / ln[2]; //ln[2]为a_1,b_1和a_2,b_2的矩阵det|A|
			lpnt3d[1] = (d1*tgtn[0] - d2*sorn[0]) / ln[2];
			lpnt3d[2] = 0.0f;
			return 4;
		}
	}

FindMaxCoordFromVec3D(ln);This is to find the longest axis.

如果上述2,3Case planes are parallel or coplanar：
First get a point on the plane,Then calculate the distance from this vertex to the plane.因此如下：
SamplePnt3dFromPlaneFnc(sor_fun, sm_pnt);This is going to any point on the plane
Pnt3d2PlaneDist(tgt_fn, sm_pnt); 计算点到平面的距离.It is used to determine whether two faces are coplanar.

There is also a better way to intersect planes,But I haven't fully understood it yet,The link to paste it is below：
http://tbirdal.blogspot.com/2016/10/a-better-approach-to-plane-intersection.html