求解直线与平面的交点
作者:互联网
求解直线与平面的交点
微信公众号:幼儿园的学霸
目录
文章目录
前言
直线与平面的交点求解相关的内容在网上已经有很多资料进行介绍,目前所看到的博文在数学模型建立上都是正确的,但是其编程实现却存在问题,导致只有部分情况下能够正确求出直线与平面的交点,另外一些情况下求出的交点却是错误的.本文对原理进行推导并实现正确的编码.
数学模型推导
已知经过两点
P
1
(
x
1
,
y
1
,
z
1
)
,
P
2
(
x
2
,
y
2
,
z
2
)
P_1(x_1,y_1,z_1),P_2(x_2,y_2,z_2)
P1(x1,y1,z1),P2(x2,y2,z2)的直线与平面
a
x
+
b
y
+
c
z
+
d
=
0
ax+by+cz+d=0
ax+by+cz+d=0相交于点P
,求点P
的坐标.
该问题几何模型如下:
求解该问题的一种方法是建立方程组,联立求解点的坐标.该方法在数学实现上非常直观,但是求方程组的解稍微麻烦.
另一种是通过向量法进行求解.推导过程如下:
设
m
=
P
1
P
→
P
1
P
2
→
,
则
O
P
→
=
O
P
1
→
+
P
1
P
→
=
O
P
1
→
+
m
∗
P
1
P
2
→
设 m = \frac{\overrightarrow {P_1P}}{\overrightarrow{P_1P_2}},则 \\ \overrightarrow {OP} = \overrightarrow {OP_1} + \overrightarrow {P_1P} =\overrightarrow {OP_1} + m*\overrightarrow {P_1P_2}
设m=P1P2
P1P
,则OP
=OP1
+P1P
=OP1
+m∗P1P2
其中:
m
=
P
1
P
→
P
1
P
2
→
=
P
1
D
→
P
1
D
2
→
=
∣
P
1
D
→
∣
∣
P
1
D
2
→
∣
m = \frac{\overrightarrow {P_1P}}{\overrightarrow{P_1P_2}} = \frac{\overrightarrow {P_1D}}{\overrightarrow{P_1D_2}} = \frac{\lvert \overrightarrow {P_1D}\rvert } {\lvert \overrightarrow{P_1D_2} \rvert}
m=P1P2
P1P
=P1D2
P1D
=∣P1D2
∣∣P1D
∣
而
∣
P
1
D
→
∣
\lvert \overrightarrow {P_1D}\rvert
∣P1D
∣为点P1
到平面的距离,计算公式如下:
∣
P
1
D
→
∣
=
∣
a
x
1
+
b
y
1
+
c
z
1
+
d
∣
a
2
+
b
2
+
c
2
\lvert \overrightarrow {P_1D}\rvert = \frac{\lvert {ax_1+by_1+cz_1+d} \rvert} {\sqrt{a^2+b^2+c^2}}
∣P1D
∣=a2+b2+c2
∣ax1+by1+cz1+d∣
另外:
∣
P
1
D
2
→
∣
=
∣
P
1
P
2
→
∣
c
o
s
θ
=
∣
P
1
P
2
→
∣
∗
∣
P
1
P
2
→
⋅
P
1
D
2
→
∣
∣
P
1
P
2
→
∣
∣
P
1
D
2
→
∣
=
∣
P
1
P
2
→
⋅
P
1
D
2
→
∣
P
1
D
2
→
∣
∣
{\lvert \overrightarrow{P_1D_2} \rvert} = {\lvert \overrightarrow{P_1P_2} \rvert}cos{\theta} = {\lvert \overrightarrow{P_1P_2} \rvert} * \frac{\lvert{\overrightarrow{P_1P_2} \cdot \overrightarrow{P_1D_2} } \rvert} {\lvert{\overrightarrow{P_1P_2}} \rvert \lvert{\overrightarrow{P_1D_2}} \rvert} \\ =\lvert{ \overrightarrow{P_1P_2} \cdot \frac{\overrightarrow{P_1D_2}}{\lvert{\overrightarrow{P_1D_2}} \rvert} } \rvert
∣P1D2
∣=∣P1P2
∣cosθ=∣P1P2
∣∗∣P1P2
∣∣P1D2
∣∣P1P2
⋅P1D2
∣=∣P1P2
⋅∣P1D2
∣P1D2
∣
显然,
P
1
D
2
→
∣
P
1
D
2
→
∣
\frac{\overrightarrow{P_1D_2}}{\lvert{\overrightarrow{P_1D_2}} \rvert}
∣P1D2
∣P1D2
表示平面的单位法向量,其用坐标的形式表示如下:
P
1
D
2
→
∣
P
1
D
2
→
∣
=
(
a
,
b
,
c
)
a
2
+
b
2
+
c
2
\frac{\overrightarrow{P_1D_2}}{\lvert{\overrightarrow{P_1D_2}} \rvert} = \frac{(a,b,c)}{\sqrt{a^2+b^2+c^2}}
∣P1D2
∣P1D2
=a2+b2+c2
(a,b,c)
因此:
∣
P
1
D
2
→
∣
=
∣
P
1
P
2
→
⋅
P
1
D
2
→
∣
P
1
D
2
→
∣
∣
=
P
1
D
2
→
∣
P
1
D
2
→
∣
=
∣
a
(
x
2
−
x
1
)
+
b
(
y
2
−
y
1
)
+
c
(
z
2
−
z
1
)
a
2
+
b
2
+
c
2
∣
{\lvert \overrightarrow{P_1D_2} \rvert} = \lvert{ \overrightarrow{P_1P_2} \cdot \frac{\overrightarrow{P_1D_2}}{\lvert{\overrightarrow{P_1D_2}} \rvert} } \rvert = \frac{\overrightarrow{P_1D_2}}{\lvert{\overrightarrow{P_1D_2}} \rvert} \\ = \lvert{ \frac{a(x_2-x_1)+b(y_2-y_1)+c(z_2-z_1)}{\sqrt{a^2+b^2+c^2}} } \rvert
∣P1D2
∣=∣P1P2
⋅∣P1D2
∣P1D2
∣=∣P1D2
∣P1D2
=∣a2+b2+c2
a(x2−x1)+b(y2−y1)+c(z2−z1)∣
综上,可得:
m
=
∣
P
1
D
→
∣
∣
P
1
D
2
→
∣
=
∣
a
x
1
+
b
y
1
+
c
z
1
+
d
a
(
x
2
−
x
1
)
+
b
(
y
2
−
y
1
)
+
c
(
z
2
−
z
1
)
∣
O
P
→
=
O
P
1
→
+
P
1
P
→
=
O
P
1
→
+
m
∗
P
1
P
2
→
或
(
x
,
y
,
z
)
=
(
x
1
,
y
1
,
z
1
)
+
m
∗
(
x
2
−
x
1
,
y
2
−
y
1
,
z
2
−
z
1
)
m = \frac{\lvert \overrightarrow {P_1D}\rvert } {\lvert \overrightarrow{P_1D_2} \rvert} = \lvert{\frac{ax_1+by_1+cz_1+d}{a(x_2-x_1)+b(y_2-y_1)+c(z_2-z_1)} } \rvert \\ \overrightarrow {OP} = \overrightarrow {OP_1} + \overrightarrow {P_1P} =\overrightarrow {OP_1} + m*\overrightarrow {P_1P_2} \\ 或 (x,y,z) = (x_1,y_1,z_1) + m*(x_2-x_1,y_2-y_1,z_2-z_1)
m=∣P1D2
∣∣P1D
∣=∣a(x2−x1)+b(y2−y1)+c(z2−z1)ax1+by1+cz1+d∣OP
=OP1
+P1P
=OP1
+m∗P1P2
或(x,y,z)=(x1,y1,z1)+m∗(x2−x1,y2−y1,z2−z1)
编程实现
/*!
* 求直线与平面的交点
* @param p1
* @param p2 直线经过的两点p1,p2
* @param plane 平面参数 ax+by+cz+d=0
* @param crossP 交点
* @return true--直线与平面相交;false--直线与平面平行(或者位于平面内)
* @author liheng
*/
bool GetLinePlaneCrossP(Eigen::Vector3f p1,Eigen::Vector3f p2,const Eigen::Vector4f& plane,Eigen::Vector3f& crossP)
{
//求分子P1D
auto P1D = plane[0]*p1[0] + plane[1]*p1[1] + plane[2]*p1[2];//ax1+by1+cz1//or采用向量点乘方式
P1D = std::abs( P1D );
auto P2D = plane[0]*p2[0] + plane[1]*p2[1] + plane[2]*p2[2];
P2D = std::abs( P2D );
//保证P2离平面比较近
if( P1D<P2D )
{
std::swap(p1,p2);//点交换
P1D = P2D;//距离交换
}
const Eigen::Vector3f P1P2 = p2-p1;
//求分母P1D2
auto P1D2 = plane[0]*P1P2[0] + plane[1]*P1P2[1] + plane[2]*P1P2[2];
P1D2 = std::abs(P1D2);
if( P1D2<FLT_EPSILON )
return false;//平行
auto m = P1D/P1D2;
crossP = p1 + m*P1P2;
return true;
}
下面函数同样可以
///*!
// * 求直线与平面的交点
// * @param p1
// * @param p2 直线经过的两点p1,p2
// * @param plane 平面参数 ax+by+cz+d=0
// * @param crossP 交点
// * @return true--直线与平面相交;false--直线与平面平行(或者位于平面内)
// * @author liheng
// */
//bool GetLinePlaneCrossP(const Eigen::Vector3f& p1,const Eigen::Vector3f& p2,const Eigen::Vector4f& plane,Eigen::Vector3f& crossP)
//{
// //求分子P1D
// auto P1D = plane[0]*p1[0] + plane[1]*p1[1] + plane[2]*p1[2];//ax1+by1+cz1//or采用向量点乘方式
//
// const Eigen::Vector3f P1P2 = p2-p1;
//
// //求分母P1D2
// auto P1D2 = plane[0]*P1P2[0] + plane[1]*P1P2[1] + plane[2]*P1P2[2];
// if( std::abs(P1D2)<FLT_EPSILON )
// return false;//平行
//
// auto m = P1D/P1D2;
// crossP = p1 - m*P1P2;//注意:此处是 负号- 同时求m时无绝对值符号
//
// return true;
//}
int main(int argc, char ** argv )
{
Eigen::Vector3f p1(0,0,5),p2(0,0,-10),crossP;
Eigen::Vector4f plane(0,0,1,0);//z=0
GetLinePlaneCrossP(p1,p2,plane,crossP);
std::cout<<"crossP:"<<crossP<<std::endl;
return 0;
}
上面的2个函数均可以求出直线与平面交点.不论输入的两点在平面两侧或位于平面同一侧,以及两点到平面的距离排序.
粗略的python版本如下
#!/usr/bin/env python3
#coding=utf-8
#============================#
#Program:main.py
#
#Date:20-6-23
#Author:liheng
#Version:V1.0
#============================#
import numpy as np
x1,y1,z1 = 0,0,10
x2,y2,z2 = 0,0,5
a,b,c,d = 0,0,1,0
# p1 = np.array([x1,y1,z1])#
# p2 = np.array([x2,y2,z2])
# plane_normal = np.array([a,b,c]) #a,b,c,d平面方程系数
#
# P1D = np.vdot(p1,plane_normal)+d
# P1D = np.abs(P1D)
#
# P2D = np.vdot(p2,plane_normal)+d
# P2D = np.abs(P2D)
#
# if P1D<P2D:#点1靠近平面,交换点1 2
# p1,p2 = p2,p1
# P1D = P2D
#
# P1D2 = np.vdot(p2-p1,plane_normal)
# P1D2 = np.abs(P1D2)
#
# n = P1D / P1D2
# p = p1 + n*(p2- p1)#所求交点
# print(p)
//or
p1 = np.array([x1,y1,z1])#
p2 = np.array([x2,y2,z2])
plane_normal = np.array([a,b,c]) #a,b,c,d平面方程系数
P1D = np.vdot(p1,plane_normal)+d
P1D2 = np.vdot(p2-p1,plane_normal)
n = P1D / P1D2
p = p1 - n*(p2- p1)#所求交点
print(p)
参考资料
标签:直线,P1,overrightarrow,rvert,求解,lvert,1D,交点,D2 来源: https://blog.csdn.net/leonardohaig/article/details/113573684