本文详细探讨了三维空间中点、线、面的表示方法，并介绍了如何计算它们之间的关系，如向量夹角、平行与垂直判断、共面性、距离计算以及交点查找。此外，还提供了相应的Python函数实现，适用于室内场景重建等领域的几何计算。
文章目录
简介
之前写了3篇室内场景重建的博客, 只是简单介绍了一下方法, 并没有对点线面的具体计算做讨论.
这篇补个坑, 用作交流学习之用.
相关的3篇博客:
基元类型及表示
本篇将重点讨论, 三维空间内, 如下三个基元(Primitive)的代数和几何关系
-
点(Point), 用一个
1x3的行向量 [ x , y , z ] [x, y, z] [x,y,z]表示, 分别表示点的三个坐标 -
线(Line), 用一个
2x3的矩阵 [ [ x 0 , y 0 , z 0 ] , [ A , B , C ] ] [[x_0, y_0, z_0], [A, B, C]] [[x0,y0,z0],[A,B,C]], 表示, 其中 [ x 0 , y 0 , z 0 ] [x_0, y_0, z_0] [x0,y0,z0]为该直线上的一个点, [ A , B , C ] [A, B, C] [A,B,C]为该直线的方向(不区分正负方向). 该直线的方程如下:
x − x 0 A = y − y 0 B = z − z 0 C \large \frac{x-x_0}{A} = \frac{y-y_0}{B} = \frac{z-z_0}{C} Ax−x0=By−y0=Cz−z0 -
面(Plane), 用一个
1x4的行向量 [ A , B , C , − 1 ] [A, B, C, -1] [A,B,C,−1]表示, 其中 [ A , B , C ] [A, B, C] [A,B,C]为该平面的法向量(不区分正负方向). 该平面的方程如下:本篇指的所有的面, 都是指矩形, 不是不规则的平面. 毕竟室内场景中, 绝大部分的面都是矩形.
A x + B y + C z − 1 = 0 \large Ax + By + Cz -1 = 0 Ax+By+Cz−1=0
本篇思想和原则
- 在计算关系的时候, 将上述三种基元统一表示成同格式的向量来计算
- 尽可能多使用
numpy封装的API, 避免for循环. 一方面是考虑到代码的美观, 另一方面是numpy经过底层加速, 速度会比for循环快不少.
关系求解
至此我们定义了3种基元的表示, 但是可以看到这三种表示是不统一的.
因此我们在计算关系的时候, 采用的最重要的思路就是将他们统一表示成同格式的向量来计算
向量的夹角
为了方便后续计算, 我们先定义求解任意两个向量夹角的函数cosine_between_vectors和arccosine两个函数
两个向量夹角公式各位应该特别熟悉, 我就贴一下不做进一步的证明了
c
o
s
i
n
e
(
α
,
β
)
=
α
∗
β
∣
∣
α
∣
∣
∗
∣
∣
β
∣
∣
\large cosine(\alpha, \beta) = \frac{\alpha * \beta}{||\alpha||*||\beta||}
cosine(α,β)=∣∣α∣∣∗∣∣β∣∣α∗β
def arccosine(cos, dtype='degree'):
"""
返回一个角度
@param cos:
@param dtype:
@return: np.float
"""
assert dtype in ['rad', 'degree'], 'Wrong value of param {dtype}'
angle = np.arccos(cos)
if 'degree' == dtype:
return np.degrees(angle)
return angle
def cosine_between_vectors(vec_1, vec_2):
"""
两个向量之间的夹角余弦
@param vec_1:
@param vec_2:
@return: float
"""
if not type(vec_1) == np.array:
vec_1 = np.array(vec_1)
assert 1 == len(vec_1.shape), 'Wrong shape of param {vec_1}' # 向量1必须是1xM的形状
if not type(vec_2) == np.array:
vec_2 = np.array(vec_2)
assert 1 == len(vec_2.shape), 'Wrong shape of param {vec_2}' # 向量2必须是1xM的形状
assert vec_1.shape[0] == vec_2.shape[0], '{vec_1} is not the same shape with {vec_2}'
normal_1_len = np.linalg.norm(vec_1) # 向量1的模
normal_2_len = np.linalg.norm(vec_2) # 向量2的模
return vec_1.dot(vec_2) / (normal_1_len * normal_2_len)
向量平行或垂直
有了上面求向量夹角的函数, 可以判断两个向量是否平行或垂直.
注意:
- 由于本篇中忽略了向量的正负方向, 因此即使夹角为180°的两个向量, 依然认为是平行的
- 实际应用中, 激光扫描得到的点会有误差, 在编程判断的时候需要体现出来这个误差
因此思路应该是, 给定两个向量 v e c 1 , v e c 2 vec_1, vec_2 vec1,vec2, 需要计算两者的夹角余弦的绝对值, 如果绝对值非常接近1, 说明两者平行(再次强调本文所有方向都忽略正反), 如果非常接近0, 说明两者垂直.
这里我们定义两个阈值COSINE_PARALLEL_THRESHOLD和COSINE_VERTICAL_THRESHOLD.
那么进一步封装, 本篇处理的不仅仅是向量, 而是3类基元(点, 线, 面). 除了点不存在平行垂直的概念, 剩下的向量, 直线, 平面都可以两两组合判断平行或者垂直. 一共有向量和向量, 向量和直线, 向量和平面, 直线和直线, 直线和平面, 平面和平面六种组合, 其中除了向量和向量可以直接做计算, 剩下的情况都需要做一点处理才能计算.
做处理的思路也很简单, 将参与计算的两个基元, 转化成两个向量做计算.
直线
原则上直线也是一种向量, 因此这种情况和上面向量和向量的流程是一样的, 唯一需要注意的是, 向量是一个行向量, 而本篇中的直线是一个矩阵, 不能直接将行向量和矩阵做处理, 应该将直线的方向向量, 也就是矩阵的第二行拿出来, 这就转化成了两个向量
平面
同样, 考虑将平面转化成一个向量参与计算, 那么首先想到的应该是平面的法向量
- 如果这个向量和平面平行, 那么这个向量应该垂直于平面的法向量
- 反过来, 如果向量和平面垂直, 那么向量和法向量应该平行.
这里唯一需要注意的是, 涉及到平面的垂直和平行的关系会变化.
那么我们很容易得到下面的is_parallel()和is_vertical()两个函数:
COSINE_PARALLEL_THRESHOLD = 0.999 # 两条向量夹角的cos值的绝对值超过这个值即认为这两条向量平行
COSINE_VERTICAL_THRESHOLD = 0.001 # 两条向量夹角的cos值的绝对值小于这个值即认为这两条向量垂直
def is_parallel(primitive_1, primitive_2, dtype=None):
"""
根据两个基元的类型, 判断两个基元是否平行
@param primitive_1:
@param primitive_2:
@param dtype: 表示基元1和基元2的类型, 由以下表示
- 'v' : 表示向量, 由向量的方向参数[x, y, z]表示
- 'l' : 表示直线, 由直线的方程参数[[x0, y0, z0], [A, B, C]]表示
- 'p' : 表示平面, 由平面的方程参数[A, B, C, -1]表示
可按下两两组合
组合\意义 基元1 基元2
- 'vv' : 向量 向量
- 'vl' : 向量 直线
- 'vp' : 向量 平面
- 'll' : 直线 直线
- 'lp' : 直线 平面
- 'pp' : 平面 平面
@return: bool
"""
if not type(primitive_1) == np.array:
primitive_1 = np.array(primitive_1)
if not type(primitive_2) == np.array:
primitive_2 = np.array(primitive_2)
assert dtype in ['vv', 'vl', 'vp', 'll', 'lp', 'pp'], 'wrong value of param {dtype}'
if 'vv' == dtype:
cos = cosine_between_vectors(primitive_1, primitive_2)
return np.abs(cos) > COSINE_PARALLEL_THRESHOLD
elif 'vl' == dtype:
vec = primitive_2[1]
cos = cosine_between_vectors(primitive_1, vec)
return np.abs(cos) > COSINE_PARALLEL_THRESHOLD
elif 'vp' == dtype:
"""
由于平面的法向量垂直于平面, 所以若向量与平面平行, 向量应该垂直于法向量
"""
normal_plane = normal_of_plane(primitive_2)
cos = cosine_between_vectors(primitive_1, normal_plane)
return np.abs(cos) < COSINE_VERTICAL_THRESHOLD
elif 'll' == dtype:
vec_1, vec_2 = primitive_1[1], primitive_2[1]
cos = cosine_between_vectors(vec_1, vec_2)
return np.abs(cos) > COSINE_PARALLEL_THRESHOLD
elif 'lp' == dtype:
"""
由于平面的法向量垂直于平面, 所以若线与平面平行, 线应该垂直于法向量
"""
vec = primitive_1[1]
normal_plane = normal_of_plane(primitive_2)
cos = cosine_between_vectors(vec, normal_plane)
return np.abs(cos) < COSINE_VERTICAL_THRESHOLD
elif 'pp' == dtype:
"""
由于平面的法向量垂直于平面, 所以若两个平面平行, 则两条法线也平行
"""
normal_plane_1 = normal_of_plane(primitive_1)
normal_plane_2 = normal_of_plane(primitive_2)
cos = cosine_between_vectors(normal_plane_1, normal_plane_2)
return np.abs(cos) > COSINE_PARALLEL_THRESHOLD
def is_vertical(primitive_1, primitive_2, dtype):
"""
根据两个基元的类型, 判断两个基元是否垂直
@param primitive_1:
@param primitive_2:
@param dtype: 表示基元1和基元2的类型, 由以下表示
- 'v' : 表示向量, 由向量的方向参数[x, y, z]表示
- 'l' : 表示直线, 由直线的方程参数[[x0, y0, z0], [A, B, C]]表示
- 'p' : 表示平面, 由平面的方程参数[A, B, C, -1]表示
可按下两两组合
组合\意义 基元1 基元2
- 'vv' : 向量 向量
- 'vl' : 向量 直线
- 'vp' : 向量 平面
- 'll' : 直线 直线
- 'lp' : 直线 平面
- 'pp' : 平面 平面
@return: bool
"""
if not type(primitive_1) == np.array:
primitive_1 = np.array(primitive_1)
if not type(primitive_2) == np.array:
primitive_2 = np.array(primitive_2)
assert dtype in ['vv', 'vl', 'vp', 'll', 'lp', 'pp'], 'wrong value of param {dtype}'
if 'vv' == dtype:
cos = cosine_between_vectors(primitive_1, primitive_2)
return np.abs(cos) < COSINE_VERTICAL_THRESHOLD
elif 'vl' == dtype:
vec = primitive_2[1]
cos = cosine_between_vectors(primitive_1, vec)
return np.abs(cos) < COSINE_VERTICAL_THRESHOLD
elif 'vp' == dtype:
"""
由于平面的法向量垂直于平面, 所以若向量与平面垂直, 向量应该平行于法向量
"""
normal_plane = normal_of_plane(primitive_2)
cos = cosine_between_vectors(primitive_1, normal_plane)
return np.abs(cos) > COSINE_PARALLEL_THRESHOLD
elif 'll' == dtype:
vec_1, vec_2 = primitive_1[1], primitive_2[1]
cos = cosine_between_vectors(vec_1, vec_2)
return np.abs(cos) < COSINE_VERTICAL_THRESHOLD
elif 'lp' == dtype:
"""
由于平面的法向量垂直于平面, 所以若线与平面垂直, 线应该平行于法向量
"""
vec = primitive_1[1]
normal_plane = normal_of_plane(primitive_2)
cos = cosine_between_vectors(vec, normal_plane)
return np.abs(cos) > COSINE_PARALLEL_THRESHOLD
elif 'pp' == dtype:
"""
由于平面的法向量垂直于平面, 所以若两个平面垂直, 则两条法线也垂直
"""
normal_plane_1 = normal_of_plane(primitive_1)
normal_plane_2 = normal_of_plane(primitive_2)
cos = cosine_between_vectors(normal_plane_1, normal_plane_2)
return np.abs(cos) < COSINE_VERTICAL_THRESHOLD
直线共面
除了平行和垂直, 两条直线的关系还涉及一个共面异面的问题.
给定两条直线
L
1
=
[
[
x
1
,
y
1
,
z
1
]
,
[
A
1
,
B
1
,
C
1
]
]
,
L
2
=
[
[
x
2
,
y
2
,
z
2
]
,
[
A
2
,
B
2
,
C
2
]
]
L_1=[[x_1, y_1, z_1], [A_1, B_1, C_1]], L_2=[[x_2, y_2, z_2], [A_2, B_2, C_2]]
L1=[[x1,y1,z1],[A1,B1,C1]],L2=[[x2,y2,z2],[A2,B2,C2]], 则他们共面的充要条件为:
∥
x
2
−
x
1
y
2
−
y
1
z
2
−
z
1
A
1
B
1
C
1
A
2
B
2
C
2
∥
=
0
\begin{Vmatrix} \large x_2- x_1 & \large y_2-y_1 & \large z_2-z_1 \\ \large A_1 & \large B_1 & \large C_1 \\ \large A_2 & \large B_2 & \large C_2 \\ \end{Vmatrix} \large = 0
∥∥∥∥∥∥x2−x1A1A2y2−y1B1B2z2−z1C1C2∥∥∥∥∥∥=0
证明如下:
取两条直线的方向向量 m 1 = [ A 1 , B 1 , C 1 ] , m 2 = [ A 2 , B 2 , C 2 ] m_1 = [A_1, B_1, C_1], m_2 = [A_2, B_2, C_2] m1=[A1,B1,C1],m2=[A2,B2,C2], 再取直线上两点连成的向量 p = [ x 2 − x 1 , y 2 − y 1 , z 2 − z 1 ] p = [x_2- x_1, y_2-y_1, z_2-z_1] p=[x2−x1,y2−y1,z2−z1]。
- 充分性， 由于上述的行列式值为0， 说明向量 p p p能够用 m 1 , m 2 m_1, m_2 m1,m2线性表示. 而三维空间中，两条向量是不足以表示第三条向量的，这与这个事实矛盾。 因此说明 m 1 , m 2 m_1, m_2 m1,m2只能在二维空间， 说明两者共面
- 必要性，由于两者在同一平面，那么向量 p p p能够用 m 1 , m 2 m_1, m_2 m1,m2线性表示， 因此该行列式值为0
函数实现如下：
def is_in_same_plane(line_1_args, line_2_args):
"""
判断两条直线是否共面
直线共面的充要条件:
|x1-x2 y1-y2 z1-z2|
| A1 B1 C1 | == 0
| A2 B2 C2 |
@param line_1_args: [[x1, y1, z1], [A1, B1, C1]]
@param line_2_args: [[x2, y2, z2], [A2, B2, C2]]
@return: bool
"""
if not type(line_1_args) == np.array:
line_1_args = np.array(line_1_args)
assert 2 == len(line_1_args.shape), 'Wrong shape of param {line_1_args}'
if not type(line_2_args) == np.array:
line_2_args = np.array(line_2_args)
assert 2 == len(line_2_args.shape), 'Wrong shape of param {line_2_args}'
martix = np.array([
line_1_args[0] - line_2_args[0],
line_1_args[1],
line_2_args[1]
])
return 0 == np.abs(np.linalg.det(martix))
距离
这是一个大头。分析一下， 点， 线， 面， 任意组合都能计算距离。因此一共有点点距离， 点线距离， 点面距离， 线线距离， 线面距离， 面面距离6种组合。
点点距离
这个不用多说, 实现如下:
def distance_point_to_point(point_1, point_2):
"""
点到点的距离
@param point_1:
@param point_2:
@return: float
"""
if not type(point_1) == np.array:
point_1 = np.array(point_1)
assert 1 == len(point_1.shape), 'Wrong shape of param {point_1}'
if not type(point_2) == np.array:
point_2 = np.array(point_2)
assert 1 == len(point_2.shape), 'Wrong shape of param {point_2}'
return np.linalg.norm(point_1 - point_2) # norm能直接返回一个向量的范数, 默认为二范数
点线距离
这里需要用到一点技巧.
已知直线 L L L上一点 M 1 M_1 M1, 直线外一点 M 0 M_0 M0, 直线的方向向量为 s s s, 求点 M 0 M_0 M0到直线的距离 d d d.
这里我们考虑图中平行四边形的面积. 众所周知, 平行四边形的面积 S S S, 数值上等于任意两条邻边的矢量积, 也就是向量的叉乘(数量上相等, 实际上向量叉乘的结果仍是向量, 方向由右手螺旋定责决定)
也就是说:
∣
∣
M
0
M
1
×
s
∣
∣
=
=
∣
∣
s
×
d
∣
∣
∣
∣
M
0
M
1
×
s
∣
∣
=
=
∣
∣
s
∣
∣
×
∣
∣
d
∣
∣
×
s
i
n
(
s
,
d
)
∣
∣
M
0
M
1
×
s
∣
∣
=
=
∣
∣
s
∣
∣
×
∣
∣
d
∣
∣
×
1
\large ||M_0M_1 \times s || == ||s \times d|| \\ \large ||M_0M_1 \times s || == ||s|| \times ||d|| \times sin(s, d) \\ \large ||M_0M_1 \times s || == ||s|| \times ||d|| \times 1 \\
∣∣M0M1×s∣∣==∣∣s×d∣∣∣∣M0M1×s∣∣==∣∣s∣∣×∣∣d∣∣×sin(s,d)∣∣M0M1×s∣∣==∣∣s∣∣×∣∣d∣∣×1
那么很容易得出: (其中
×
\times
×表示向量叉乘, 不是一般的乘法)
∣
∣
d
∣
∣
=
=
∣
∣
M
0
M
1
×
s
∣
∣
∣
∣
s
∣
∣
\large ||d|| == \frac{||M_0M_1 \times s||}{||s||} \\
∣∣d∣∣==∣∣s∣∣∣∣M0M1×s∣∣
def distance_point_to_line(point, line_args):
"""
点到直线的距离
|MP x N|/|N|
@param point: [x, y, z]
@param line_args: [[x0, y0, z0], [A, B, C]]
"""
if not type(point) == np.array:
point = np.array(point)
assert 1 == len(point.shape), 'Wrong shape of {point}'
if not type(line_args) == np.array:
line_args = np.array(line_args)
assert 2 == len(line_args.shape), 'Wrong shape of param {line_args}'
MP = line_args[0] - point
normal_line = line_args[1]
# np.cross()可以计算两个向量的矢量积, 也就是叉乘
return np.linalg.norm(np.cross(MP, normal_line)) / np.linalg.norm(normal_line)
点面距离
中学学过二维空间内的点线距离:
d
=
∣
A
x
0
+
B
y
0
+
C
∣
A
2
+
B
2
\large d = \frac{|Ax_0+By_0+C|}{\sqrt{A^2+B^2}}
d=A2+B2∣Ax0+By0+C∣
其实二维空间内的线, 就是三维空间内的面, 因此三维的点面距离为:
d
=
∣
A
x
0
+
B
y
0
+
C
z
0
+
D
∣
A
2
+
B
2
+
C
2
\large d = \frac{|Ax_0+By_0+Cz_0 + D|}{\sqrt{A^2+B^2+C^2}}
d=A2+B2+C2∣Ax0+By0+Cz0+D∣
def distance_point_to_plane(point, plane_args):
"""
点到面的距离
@param point:
@param plane_args:
@return: float
"""
if not type(point) == np.array:
point = np.array(point)
assert 1 == len(point.shape), 'Wrong shape of param {point}'
if not type(plane_args) == np.array:
plane_args = np.array(plane_args)
assert 1 == len(plane_args.shape), 'Wrong shape of param {plane_args}'
a = np.hstack([point, 1])
b = np.linalg.norm(normal_of_plane(plane_args))
return np.abs(a.dot(plane_args)) / b
线线距离
乍看很难, 其实很简单.
首先, 验证这两条线是否平行, 不平行不存在线线距离这么一说.
然和, 只要求其中一条线上的某一点, 到另一条直线的距离即可.
def distance_line_to_line(line_1_args, line_2_args):
"""
直线到直线的距离
需要两条直线平行
该距离 == 其中一条直线上的点, 到另一条直线的距离
@param line_1_args: [[x10, y10, z10], [A1, B1, C1]]
@param line_2_args: [[x20, y20, z20], [A2, B2, C2]]
"""
if not type(line_1_args) == np.array:
line_1_args = np.array(line_1_args)
assert 2 == len(line_1_args.shape), 'Wrong shape of param {line_1_args}'
if not type(line_2_args) == np.array:
line_2_args = np.array(line_2_args)
assert 2 == len(line_2_args.shape), 'Wrong shape of param {line_2_argse}'
assert is_parallel(line_1_args, line_2_args, dtype='ll'), \
'Abort! {line_1_args} and {line_2_args} are not parallel.'
return distance_point_to_line(line_1_args[0], line_2_args)
线面距离
同上, 验证平行, 然后求点到面的距离.
def distance_line_to_plane(line_args, plane_args):
"""
直线到平面的距离
需要直线和平面平行
确认平行后, 直接返回点[x0, y0, z0]到平面plane_args的距离即可
@param line_args: [[x0, y0, z0], [A, B, C]]
@param plane_args: [A, B, C, D, -1]
@return: float
"""
if not type(line_args) == np.array:
line_args = np.array(line_args)
assert 2 == len(line_args.shape), 'Wrong shape of param {line_args}'
if not type(plane_args) == np.array:
plane_args = np.array(plane_args)
assert 1 == len(plane_args.shape), 'Wrong shape of param {plane_args}'
assert is_parallel(line_args, plane_args, dtype='lp'), \
'Abort! line {line_args} and plane {plane_args} are not parallel.'
return distance_point_to_plane(line_args[0], plane_args)
面面距离
前面说过, 三维的面, 就是二维的线. 而二维空间中, 直线之间的距离公式为:
d
=
∣
C
1
−
C
2
∣
A
2
+
B
2
\large d = \frac{|C_1 - C_2|}{\sqrt{A^2+B^2}}
d=A2+B2∣C1−C2∣
因此, 三维的面面距离为:
d
=
∣
D
1
−
D
2
∣
A
2
+
B
2
+
C
2
\large d = \frac{|D_1 - D_2|}{\sqrt{A^2+B^2+C^2}}
d=A2+B2+C2∣D1−D2∣
但是, 这个方法要求两个平面的
A
,
B
,
C
A, B, C
A,B,C相等., 而我是让两个平面的
D
D
D相等, 因此需要做一点变换. 将两个平面方程参数中的
A
,
B
,
C
A, B, C
A,B,C缩放到相等, 才能继续计算.
同理, 需要验证两个平面平行.
def distance_plane_to_plane(plane_1_args, plane_2_args):
"""
计算两个平面之间的距离
要求两个平面平行
平行平面距离公式 |D1-D2|/sqrt(A^2+B^2+C^2)
@param plane_1_args: [A1, B1, C1, -1]
@param plane_2_args: [A2, B2, C2, -1]
@return: float
"""
if not type(plane_1_args) == np.array:
plane_1_args = np.array(plane_1_args)
assert 1 == len(plane_1_args.shape), 'Wrong shape of param {plane_1_args}'
if not type(plane_2_args) == np.array:
plane_2_args = np.array(plane_2_args)
assert 1 == len(plane_2_args.shape), 'Wrong shape of param {plane_2_args}'
assert is_parallel(plane_1_args, plane_2_args, dtype='pp'), \
'Abort! {plane_1_args} and {plane_2_args} are not parallel.'
# 这里要求两个平面的A,B,C相同, 因此先放大两者的A, B, C到相等
k = plane_1_args[0] / plane_2_args[0]
plane_tmp = plane_2_args * k
return np.abs(plane_tmp[3] - plane_1_args[3]) / np.linalg.norm(normal_of_plane(plane_1_args))
两个平面求交线
已知平面 P 1 = [ A 1 , B 1 , C 1 , − 1 ] , P 2 = [ A 2 , B 2 , C 2 , − 1 ] P_1 = [A_1, B_1, C_1, -1], P_2=[A_2, B_2, C_2, -1] P1=[A1,B1,C1,−1],P2=[A2,B2,C2,−1], 确定他们的交线 L = [ [ x 0 , y 0 , z 0 ] , [ A , B , C ] ] L=[[x_0, y_0, z_0], [A, B, C]] L=[[x0,y0,z0],[A,B,C]].
这里设交线 L L L的方向向量为 s s s, P 1 , P 2 P_1, P_2 P1,P2的法向量为 m 1 , m 2 m_1, m_2 m1,m2. 由于 L L L同时在平面 P 1 , P 2 P_1, P_2 P1,P2内, 因此 L L L与两个平面都平行. 也就是说, s s s与 m 1 , m 2 m_1, m_2 m1,m2都垂直. 那么怎么由两个已知的向量, 计算得到与他们都垂直的一个新向量?
答案显然是向量矢量积, 也就是前面的叉乘 × \times ×, 矢量积的结果的方向, 由右手螺旋定责确定, 因此与两个参数向量是垂直的.
那么确定了交线的方向, 怎么确定交线上的一点呢?
也很简单. 设交点为 P = [ x 0 , y 0 , z 0 ] P=[x_0, y_0, z_0] P=[x0,y0,z0], 固定 x 0 = k x_0 = k x0=k, 去求剩余的 y 0 , z 0 y_0, z_0 y0,z0. 两个未知数, 两个方程. 不难.
def inter_line_of_two_planes(plane_1_args, plane_2_args):
"""
求两个平面交线的方程
将两个平面的法向量叉乘, 得到交线的方向向量
再在交线上任取一点, 即可确定交线的方程
平面由[A, B, C, -1]表示组成
@param plane_1_args: [A1, B1, C1, -1]
@param plane_2_args: [A2, B2, C2, -1]
@return: np.ndarray
"""
assert 1 == len(plane_1_args.shape), 'Wrong shape of param {plane_1_args}'
assert 1 == len(plane_2_args.shape), 'Wrong shape of param {plane_2_args}'
plane_normal_1 = normal_of_plane(plane_1_args)
plane_normal_2 = normal_of_plane(plane_2_args)
# 交线垂直于两个平面的法向量, 因此可以用两者的叉乘表示
line_normal = np.cross(plane_normal_1, plane_normal_2)
# 固定x0 = k, 求y0和z0
k = 0
equal_1 = np.array([plane_1_args[1], plane_1_args[2]])
equal_2 = np.array([plane_2_args[1], plane_2_args[2]])
P = np.array([equal_1, equal_2])
b = np.array([1 - plane_1_args[0] * k, 1 - plane_2_args[0] * k])
solve = np.linalg.solve(P, b)
return np.array([
[k, solve[0], solve[1]], # 交线上某个点的坐标
line_normal # 交线的方向
])
三个平面求交点
这就更简单了, 已知3个平面 P 1 = [ A 1 , B 1 , C 1 , − 1 ] , P 2 = [ A 2 , B 2 , C 2 , − 1 ] , P 3 = [ A 3 , B 3 , C 3 , − 1 ] P_1=[A_1, B_1, C_1, -1], P_2=[A_2, B_2, C_2, -1], P_3=[A_3, B_3, C_3, -1] P1=[A1,B1,C1,−1],P2=[A2,B2,C2,−1],P3=[A3,B3,C3,−1], 求它们的交点 P o = [ x , y , z ] P_o=[x, y, z] Po=[x,y,z]
那么有:
P Q = ( 1 1 1 ) PQ = \begin{pmatrix} 1 \\ 1 \\ 1\\ \end{pmatrix} PQ=⎝⎛111⎠⎞
其中,
P
=
(
A
1
B
1
C
1
A
2
B
2
C
2
A
3
B
3
C
3
)
,
Q
=
(
x
y
z
)
P= \begin{pmatrix} A_1 & B_1 & C_1 \\ A_2 & B_2 & C_2 \\ A_3 & B_3 & C_3 \\ \end{pmatrix}, Q = \begin{pmatrix} x \\ y \\ z \\ \end{pmatrix}
P=⎝⎛A1A2A3B1B2B3C1C2C3⎠⎞,Q=⎝⎛xyz⎠⎞
很容易求得上述线性方程的解
Q
Q
Q, 如下:
def inter_point_of_three_planes(plane_1_args, plane_2_args, plane_3_args):
"""
求三个平面的交点
平面由[A, B, C, -1]表示组成
@param plane_1_args: [A1, B1, C1, -1]
@param plane_2_args: [A2, B2, C2, -1]
@param plane_3_args: [A3, B3, C3, -1]
@return: np.ndarray
"""
assert 1 == len(plane_1_args.shape), 'Wrong shape of param {plane_1_args}'
assert 1 == len(plane_2_args.shape), 'Wrong shape of param {plane_2_args}'
assert 1 == len(plane_3_args.shape), 'Wrong shape of param {plane_3_args}'
P = normals_of_planes([plane_1_args, plane_2_args, plane_3_args])
ones = np.array([1, 1, 1])
return np.linalg.solve(P, ones)
两条直线求交点
这里我是直接手动联立两个方程求解的, 得到:
def inter_point_of_two_lines(line_1_args, line_2_args):
"""
求两条直线的交点
要求两条直线共面且不平行
@param line_1_args:
@param line_2_args:
@return: np.ndarray
"""
if not type(line_1_args) == np.array:
line_1_args = np.array(line_1_args)
assert 2 == len(line_1_args.shape), 'Wrong shape of param {line_1_args}'
if not type(line_2_args) == np.array:
line_2_args = np.array(line_2_args)
assert 2 == len(line_2_args.shape), 'Wrong shape of param {line_2_args}'
assert is_in_same_plane(line_1_args, line_2_args), \
'Abort! Line_1 and line_2 are not in the same plane.'
assert not is_parallel(line_1_args, line_2_args, dtype='ll'), \
'Abort! Line_1 and line_2 are parallel.'
(x1, y1, z1), (A1, B1, C1) = line_1_args[0], line_1_args[1]
(x2, y2, z2), (A2, B2, C2) = line_2_args[0], line_2_args[1]
x = (A2 * x1 - A1 * x2) / (A2 - A1)
y = (B2 * y1 - B1 * y2) / (B2 - B1)
z = (C2 * z1 - C1 * z2) / (C2 - C1)
return np.array([x, y, z])
平面面积
前面说过, 这里的平面指的是矩形. 鉴于矩形也是一种平行四边形. 因此我们只需要知道平面的两条邻边向量即可. 而很不幸, 本篇定义的平面表示中, 不存在这样的两个向量.
因此这里需要引入其他的方法. 利用平面的4个边界点
利用这些边界点, 很容易求的平面的一组邻边向量, 直接叉乘即可.
也可以采用对角线叉乘。
def areas_of_planes(planes_bounds):
"""
计算所有平面的面积
单个平面面积 == |ACxBD|/2
planes_bounds需要保证:
图形上相邻的点, 在数组上也相邻
图形上处于对角的点, 在数组上索引相差2
@param planes_bounds: 平面的四个边界点
@return : np.ndarray
"""
if not type(planes_bounds) == np.array:
planes_bounds = np.array(planes_bounds)
assert 3 == len(planes_bounds.shape), 'Wrong shape of param {planes_bounds}'
# shape: (n,), (n,) (n,) (n,)
A, B, C, D = planes_bounds[:, 0], planes_bounds[:, 1], planes_bounds[:, 2], planes_bounds[:, 3]
"""
对角线向量
"""
AC = C - A
BD = D - B
# 四边形面积 == 对角线向量叉乘的模的二分之一
return np.linalg.norm(np.cross(AC, BD), axis=1, keepdims=True) / 2
平面方程
前面引入了平面的边界点, 那么问题又来了? 怎么通过一个平面, 去求它的4个边界点?这个问题在几何中是不可解的. 因为几何中的平面本身是无限的, 不存在边界点.
因此实际应用中往往是反过来的. 通过4个边界点, 去求这4点确定的平面的方程.(其实3个点就够了)
假设已知三个点 P 1 = [ x 1 , y 1 , z 1 ] , P 2 = [ x 2 , y 2 , z 2 ] , P 3 = [ x 3 , y 3 , z 3 ] P_1=[x_1, y_1, z_1], P_2=[x_2, y_2, z_2], P_3=[x_3, y_3, z_3] P1=[x1,y1,z1],P2=[x2,y2,z2],P3=[x3,y3,z3], 求它们所在的平面 P o = [ A , B , C , − 1 ] P_o=[A, B, C, -1] Po=[A,B,C,−1].
即求解线性方程组
P
Q
=
(
1
1
1
)
PQ = \begin{pmatrix} 1 \\ 1 \\ 1\\ \end{pmatrix}
PQ=⎝⎛111⎠⎞
其中,
P
=
(
x
1
y
1
z
1
x
2
y
2
z
2
x
3
y
3
z
3
)
,
Q
=
(
A
B
C
)
P= \begin{pmatrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ x_3 & y_3 & z_3 \\ \end{pmatrix}, Q = \begin{pmatrix} A \\ B \\ C \\ \end{pmatrix}
P=⎝⎛x1x2x3y1y2y3z1z2z3⎠⎞,Q=⎝⎛ABC⎠⎞
def args_of_planes(planes_bounds):
"""
通过每个平面的四个顶点的前三个解析一个平面的参数
Ax+By+Cz-1 = 0
PQ = 0
其中P = bound[0:3]
每个平面返回np.array([A,B,C,-1])
@param planes_bounds:
@return: np.array
"""
if not type(planes_bounds) == np.array:
planes_bounds = np.array(planes_bounds)
assert 3 == len(planes_bounds.shape), 'Wrong shape of param {planes_bounds}'
ans = []
ones = np.array([1, 1, 1])
for bound in planes_bounds:
P = np.array([bound[0], bound[1], bound[2]])
ans.append(np.hstack([np.linalg.solve(P, ones), [-1]]))
return np.array(ans, dtype=np.float)
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
