空间平面相交的直线的计算及其源码

目的：在写图形学项目中，常常遇到三维平面相交的情况。它也是图形学的基础

空间平面和直线一样，主要两种表达方式。
1）通用的代数表达方式： ax+by+cz+d=0ax+by+cz+d=0ax+by+cz+d=0
2）几何的表达方式(都是列向量)：(P,N)(P,N)(P,N)，其中为平面的某个点，NNN为平面的法向量。它的转化为代数的表达为(P−O)T∗N=0(P-O)^T*N=0(P−O)T∗N=0（需要转置为行向量）
在后续的计算中，在写C++代码中，一般倾向于几何表达。因为它比较直观，同时也利于后续的计算。下面我将基于几何的表达形式介绍一下两个平面相交的算法，最后给予C++代码。

方法一：
假设两个平面为：

(P1,N1)=>(X−P1)T∗N1=0(1)(P_1,N_1)=>(X-P_1)^T*N_1=0 \space \space \space(1)(P1,N1)=>(X−P1)T∗N1=0 (1)

(P2,N2)=>(X−P2)T∗N2=0(2)(P_2,N_2)=>(X-P_2)^T*N_2=0 \space \space \space(2)(P2,N2)=>(X−P2)T∗N2=0 (2)

两个空间平面相交情况有三种
1：相交为直线。
2：如果平面平行，无相交。
3：者共面的情况下，相交，但是不为平面。
下面介绍一下平面相交为直线的情况：

两个平面的相交的直线的表达：O3O_3O3 (起始顶点)；N3N_3N3(方向)。
1）首先计算交线的方向N3N_3N3, 因为发现交线的方向同时垂直于plane1和plane2两个平面的法向量。叉乘N1×N2N_1\times N_2N1×N2,得到如下：
N3=N1×N2(3)N_3 = N_1\times N_2 \space \space (3)N3=N1×N2 (3)

红色即为得到的平面交线的法向量。(图像中的)

2）需要找到一个在直线上的顶点。它可以是任意的顶点，最好的是P1P_1P1,P2P_2P2顶点投影到直线上P1′,P2′P_1',P_2'P1′,P2′的中点。为了方便讲述，本博客只利用了P1P_1P1投影到直线上的点PPP,如下图

它是平面(P1,N3)(P_1,N_3)(P1,N3)和直线的交点。将P1,N3P_1,N_3P1,N3写成如下格式：

(P1,N3)=>(X−P1)T∗N3=0(4)(P_1,N_3)=>(X-P_1)^T*N_3=0 \space \space \space(4)(P1,N3)=>(X−P1)T∗N3=0 (4)

因此用(1)(2)(4)可以形成一个线性方程组为:

组成新的平面。然后利用三个平面相交得到顶点。如下图示意：
第三个法向量可以表达为：
(X−P1)T∗N1=0(X−P2)T∗N2=0(X−P1)T∗N3=0\begin{gathered} (X-P_1)^T*N_1=0 \\ (X-P_2)^T*N_2=0\\ (X-P_1)^T*N_3=0 \end{gathered} (X−P1)T∗N1=0(X−P2)T∗N2=0(X−P1)T∗N3=0

拆开公式可以得到如下：
XT∗N1=P1T∗N1XT∗N2=P2T∗N2XT∗N3=P1T∗N3\begin{gathered} X^T*N_1=P_1^T*N_1 \\ X^T*N_2=P_2^T*N_2\\ X^T*N_3=P_1^T*N_3 \end{gathered} XT∗N1=P1T∗N1XT∗N2=P2T∗N2XT∗N3=P1T∗N3

其中XXX为未知数，它是向量，通过上述公式可以得到直线上的顶点X∗X^*X∗，它的值也就是P=X∗P=X^*P=X∗。
这部分的代码后续在添加上。也是非常简单如下:

#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);//判断两个平面是否共面，或者平行if ((std::abs(ldir[0]) + std::abs(ldir[1]) + std::abs(ldir[2])) < 1e-6)    //判断两个平面是否平行{float dist = Pnt3d2PlaneDist(tgt_pnt, tgt_dir, src_pnt);if (dist < 1e-4)return 0;  //共面elsereturn 1;   //平行面}//计算投影到直线的顶点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 matrixint 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;
}

运行的结果如下：

方法二：
如果函数的表达为 ax+by+cz+d=0ax+by+cz+d=0ax+by+cz+d=0。两个平面的表达式为
a1x+b1y+c1z+d1=0(5)a_1x+b_1y+c_1z+d_1=0 \space \space \space (5)a1x+b1y+c1z+d1=0 (5)

a2x+b2y+c2z+d2=0(6)a_2x+b_2y+c_2z+d_2=0 \space \space \space (6)a2x+b2y+c2z+d2=0 (6)

它可以写成向量的形式为:

N1T×X=−d1(7)N_1^T \times X=-d_1 \space \space \space (7)N1T×X=−d1 (7)

N2T×X=−d2(8)N_2^T \times X=-d_2 \space \space \space (8)N2T×X=−d2 (8)

1：计算直线方向：使用同样的方法，可以得到N3=N1×N2N_3=N_1 \times N_2N3=N1×N2，
2：在直线上找到一个顶点：因为两个平面都是顶点，因此不会考虑寻找最好的顶点。我们使用youtube上的一个视频，它的介绍如下：

假设其中的一根轴为0(z=0z=0z=0)。然后组成下面的公式。
a1x+b1y+d1=0a2x+b2y+d2=0\begin{gathered} a_1x+b_1y+d_1=0 \\ a_2x+b_2y+d_2=0 \end{gathered} a1x+b1y+d1=0a2x+b2y+d2=0
这样得到的顶点为P=(0,ys,zs)P=(0,y_s,z_s)P=(0,ys,zs)
但是在实际的代码中，对于轴选为000需要考量。我们知道，对于一根轴，如果它几乎没有任何贡献，可能会出现问题。就如果在方程中a1x+b1y+c1z+d1=0(5)a_1x+b_1y+c_1z+d_1=0 \space \space \space (5)a1x+b1y+c1z+d1=0 (5)，如果c1=10000，a1=0.0001,b1=100c_1=10000，a_1=0.0001, b_1=100c1=10000，a1=0.0001,b1=100，那么可能会发现计算的xxx会非常大。这样的交线在xxx轴上系数非常小，表示它跟xxx轴平行。因此在实际的计算交线时候，我们需要通过计算平面交线的主轴。它有法向量N3N_3N3表示。
如果N3N_3N3中xxx系数最大。我们假设交线为P=(0,y,z)P=(0,y,z)P=(0,y,z)，然后带入方程为:

b1y+c1z+d1=0b2y+c2z+d2=0\begin{gathered} b_1y+c_1z+d_1=0 \\ b_2y+c_2z+d_2=0 \end{gathered} b1y+c1z+d1=0b2y+c2z+d2=0

如果N3N_3N3中yyy系数最大。我们假设交线为P=(x,0,z)P=(x,0,z)P=(x,0,z)
a1x+c1z+d1=0a2x+c2z+d2=0\begin{gathered} a_1x+c_1z+d_1=0 \\ a_2x+c_2z+d_2=0 \end{gathered} a1x+c1z+d1=0a2x+c2z+d2=0

如果N3N_3N3中zzz系数最大。我们假设交线为P=(0,y,z)P=(0,y,z)P=(0,y,z)
a1y+b1z+d1=0a2y+b2z+d2=0\begin{gathered} a_1y+b_1z+d_1=0 \\ a_2y+b_2z+d_2=0 \end{gathered} a1y+b1z+d1=0a2y+b2z+d2=0

具体实现代码如下：

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;//判断三根轴，找到最大的一根轴if ((std::abs(ln[0]) + std::abs(ln[1]) + std::abs(ln[2])) < 1e-6) //判断两个平面是否平行{Eigen::Vector3f sm_pnt;SamplePnt3dFromPlaneFnc(sor_fun, sm_pnt);float dist = Pnt3d2PlaneDist(tgt_fn, sm_pnt);if (dist < 1e-4)return 0; //共面elsereturn 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);这是找最长的一根轴。

如果上述2,3情况平面平行或者共面：
先要在平面上取到一点，然后计算这个顶点到平面的距离即可。因此如下：
SamplePnt3dFromPlaneFnc(sor_fun, sm_pnt);这是在平面上任意去一点
Pnt3d2PlaneDist(tgt_fn, sm_pnt); 计算点到平面的距离。它用于判断两个面是否共面。

还有一个更加好的平面相交的方法，但是自己暂且没有完全看懂，贴上它的链接如下：
http://tbirdal.blogspot.com/2016/10/a-better-approach-to-plane-intersection.html

空间平面相交的直线的计算及其源码相关推荐

两直线平行交叉相乘_计算几何算法5. 直线、线段和平面相交（2D和3D）
直线和线段相交平面相交直线-平面相交两平面相交三个平面相交实现 intersect2D_2Segments() inSegment() intersect3D_SegmentPlane() ...
c++ 圆上任意点坐标计算_线性代数总结第三章向量代数与几何计算（空间平面和直线）...
我的公众号"每日晴天",可关注领取我的笔记pdf版哦~ -------------------------------------------------------------- ...
平面/空间杆系结构有限元编程计算（MATLAB）
GitHub链接:MATLAB源码博主QQ:915339719,有问题可以随时交流程序简介本程序可以对绝大多数的平面/空间杆系结构进行静力学计算.自振频率计算,以及平面杆系结构的稳定性计算,具有 ...
python计算点到直线的距离_Python求平面内点到直线距离的实现
近期遇到个问题,需要计算平面内点到直线的距离,发现数学知识都还给老师了,度娘后找到计算方法,特此记录. 点到直线的计算公式: 通过公式推导,得到信息: a:直线斜率 b:固定值-1 c:直线截距b 转 ...
笛卡尔空间轨迹规划（直线、圆弧）
目录毕设(5)-笛卡尔空间轨迹规划(直线.圆弧) 直线轨迹规划圆弧轨迹规划 Matlab代码验证毕设中用到了很多代码,其中一部分我通过看书和看论文学习并实现的代码,会通过Gitee仓库分享出来, ...
二维平面内两直线交点计算
对于几何关系的计算,两点确定一条直线,两直线的交点的计算经常见到.最简单matlab编程方法是通过两点列出直线方程,使用matlab工具解方程.但是计算时发现速度较为缓慢. 表达式定义空间直线的表达 ...
如何判断两个平面相交_数学提高平面与平面垂直的判定方法是什么
一般地,两个平面相交,如果它们所成的二面角是直二面角,就说这两个平面互相垂直.一个平面过另一个平面的垂线,则这两个平面垂直. 平面与平面垂直的判定方法 1.定义法:如果两个平面所成的二面角为90°,那 ...
python 求平面两点距离_Python求平面内点到直线距离的实现
近期遇到个问题,需要计算平面内点到直线的距离,发现数学知识都还给老师了,度娘后找到计算方法,特此记录. 点到直线的计算公式: 通过公式推导,得到信息: A:直线斜率 B:固定值-1 C:直线截距b 转 ...
高中数学必修2：平面解析几何之直线与圆、圆与圆的位置关系
今天分享关于高中数学必修2平面解析几何中的直线与圆.圆与圆的位置关系知识点,分别从三个方面讲解,用5个经典习题进行解答全过程. 1.直线与圆的位置关系 2.圆与圆的位置关系 2．求圆的弦长的常用方法 ...

空间平面相交的直线的计算及其源码

目的：在写图形学项目中，常常遇到三维平面相交的情况。它也是图形学的基础

空间平面相交的直线的计算及其源码相关推荐

最新文章

热门文章