三维空间 点线面解析
文章目录
- 简介
- 基元类型及表示
- 本篇思想和原则
- 关系求解
- 向量的夹角
- 向量平行或垂直
- 直线
- 平面
- 直线共面
- 距离
- 点点距离
- 点线距离
- 点面距离
- 线线距离
- 线面距离
- 面面距离
- 两个平面求交线
- 三个平面求交点
- 两条直线求交点
- 平面面积
- 平面方程
简介
之前写了3篇室内场景重建的博客, 只是简单介绍了一下方法, 并没有对点线面的具体计算做讨论.
这篇补个坑, 用作交流学习之用.
相关的3篇博客:
- 室内场景重建一 粗略平面拟合
- 室内场景重建二 立方体拟合
- 室内场景重建三 点面分割
基元类型及表示
本篇将重点讨论, 三维空间内, 如下三个基元(Primitive)的代数和几何关系
点(Point), 用一个
1x3的行向量[x,y,z][x, y, z][x,y,z]表示, 分别表示点的三个坐标线(Line), 用一个
2x3的矩阵[[x0,y0,z0],[A,B,C]][[x_0, y_0, z_0], [A, B, C]][[x0,y0,z0],[A,B,C]], 表示, 其中[x0,y0,z0][x_0, y_0, z_0][x0,y0,z0]为该直线上的一个点, [A,B,C][A, B, C][A,B,C]为该直线的方向(不区分正负方向). 该直线的方程如下:
x−x0A=y−y0B=z−z0C\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]为该平面的法向量(不区分正负方向). 该平面的方程如下:本篇指的所有的面, 都是指矩形, 不是不规则的平面. 毕竟室内场景中, 绝大部分的面都是矩形.
Ax+By+Cz−1=0\large Ax + By + Cz -1 = 0 Ax+By+Cz−1=0
本篇思想和原则
- 在计算关系的时候, 将上述三种基元统一表示成同格式的向量来计算
- 尽可能多使用
numpy封装的API, 避免for循环. 一方面是考虑到代码的美观, 另一方面是numpy经过底层加速, 速度会比for循环快不少.
关系求解
至此我们定义了3种基元的表示, 但是可以看到这三种表示是不统一的.
因此我们在计算关系的时候, 采用的最重要的思路就是将他们统一表示成同格式的向量来计算
向量的夹角
为了方便后续计算, 我们先定义求解任意两个向量夹角的函数cosine_between_vectors和arccosine两个函数
两个向量夹角公式各位应该特别熟悉, 我就贴一下不做进一步的证明了
cosine(α,β)=α∗β∣∣α∣∣∗∣∣β∣∣\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 angledef 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°的两个向量, 依然认为是平行的
- 实际应用中, 激光扫描得到的点会有误差, 在编程判断的时候需要体现出来这个误差
因此思路应该是, 给定两个向量vec1,vec2vec_1, vec_2vec1,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_THRESHOLDelif 'vl' == dtype:vec = primitive_2[1]cos = cosine_between_vectors(primitive_1, vec)return np.abs(cos) > COSINE_PARALLEL_THRESHOLDelif 'vp' == dtype:"""由于平面的法向量垂直于平面, 所以若向量与平面平行, 向量应该垂直于法向量"""normal_plane = normal_of_plane(primitive_2)cos = cosine_between_vectors(primitive_1, normal_plane)return np.abs(cos) < COSINE_VERTICAL_THRESHOLDelif '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_THRESHOLDelif '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_THRESHOLDelif '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_THRESHOLDdef 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_THRESHOLDelif 'vl' == dtype:vec = primitive_2[1]cos = cosine_between_vectors(primitive_1, vec)return np.abs(cos) < COSINE_VERTICAL_THRESHOLDelif 'vp' == dtype:"""由于平面的法向量垂直于平面, 所以若向量与平面垂直, 向量应该平行于法向量"""normal_plane = normal_of_plane(primitive_2)cos = cosine_between_vectors(primitive_1, normal_plane)return np.abs(cos) > COSINE_PARALLEL_THRESHOLDelif '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_THRESHOLDelif '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_THRESHOLDelif '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
直线共面
除了平行和垂直, 两条直线的关系还涉及一个共面异面的问题.
给定两条直线L1=[[x1,y1,z1],[A1,B1,C1]],L2=[[x2,y2,z2],[A2,B2,C2]]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]], 则他们共面的充要条件为:
∥x2−x1y2−y1z2−z1A1B1C1A2B2C2∥=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
证明如下:
取两条直线的方向向量m1=[A1,B1,C1],m2=[A2,B2,C2]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=[x2−x1,y2−y1,z2−z1]p = [x_2- x_1, y_2-y_1, z_2-z_1]p=[x2−x1,y2−y1,z2−z1]。
- 充分性， 由于上述的行列式值为0， 说明向量ppp能够用m1,m2m_1, m_2m1,m2线性表示. 而三维空间中，两条向量是不足以表示第三条向量的，这与这个事实矛盾。 因此说明m1,m2m_1, m_2m1,m2只能在二维空间， 说明两者共面
- 必要性，由于两者在同一平面，那么向量ppp能够用m1,m2m_1, m_2m1,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能直接返回一个向量的范数, 默认为二范数
点线距离
这里需要用到一点技巧.
已知直线LLL上一点M1M_1M1, 直线外一点M0M_0M0, 直线的方向向量为sss, 求点M0M_0M0到直线的距离ddd.
这里我们考虑图中平行四边形的面积. 众所周知, 平行四边形的面积SSS, 数值上等于任意两条邻边的矢量积, 也就是向量的叉乘(数量上相等, 实际上向量叉乘的结果仍是向量, 方向由右手螺旋定责决定)
也就是说:
∣∣M0M1×s∣∣==∣∣s×d∣∣∣∣M0M1×s∣∣==∣∣s∣∣×∣∣d∣∣×sin(s,d)∣∣M0M1×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∣∣==∣∣M0M1×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] - pointnormal_line = line_args[1]# np.cross()可以计算两个向量的矢量积, 也就是叉乘return np.linalg.norm(np.cross(MP, normal_line)) / np.linalg.norm(normal_line)
点面距离
中学学过二维空间内的点线距离:
d=∣Ax0+By0+C∣A2+B2\large d = \frac{|Ax_0+By_0+C|}{\sqrt{A^2+B^2}} d=A2+B2∣Ax0+By0+C∣
其实二维空间内的线, 就是三维空间内的面, 因此三维的点面距离为:
d=∣Ax0+By0+Cz0+D∣A2+B2+C2\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=∣C1−C2∣A2+B2\large d = \frac{|C_1 - C_2|}{\sqrt{A^2+B^2}} d=A2+B2∣C1−C2∣
因此, 三维的面面距离为:
d=∣D1−D2∣A2+B2+C2\large d = \frac{|D_1 - D_2|}{\sqrt{A^2+B^2+C^2}} d=A2+B2+C2∣D1−D2∣
但是, 这个方法要求两个平面的A,B,CA, B, CA,B,C相等., 而我是让两个平面的DDD相等, 因此需要做一点变换. 将两个平面方程参数中的A,B,CA, B, CA,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 * kreturn np.abs(plane_tmp[3] - plane_1_args[3]) / np.linalg.norm(normal_of_plane(plane_1_args))
两个平面求交线
已知平面P1=[A1,B1,C1,−1],P2=[A2,B2,C2,−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=[[x0,y0,z0],[A,B,C]]L=[[x_0, y_0, z_0], [A, B, C]]L=[[x0,y0,z0],[A,B,C]].
这里设交线LLL的方向向量为sss, P1,P2P_1, P_2P1,P2的法向量为m1,m2m_1, m_2m1,m2. 由于LLL同时在平面P1,P2P_1, P_2P1,P2内, 因此LLL与两个平面都平行. 也就是说, sss与m1,m2m_1, m_2m1,m2都垂直. 那么怎么由两个已知的向量, 计算得到与他们都垂直的一个新向量?
答案显然是向量矢量积, 也就是前面的叉乘×\times×, 矢量积的结果的方向, 由右手螺旋定责确定, 因此与两个参数向量是垂直的.
那么确定了交线的方向, 怎么确定交线上的一点呢?
也很简单. 设交点为P=[x0,y0,z0]P=[x_0, y_0, z_0]P=[x0,y0,z0], 固定x0=kx_0 = kx0=k, 去求剩余的y0,z0y_0, z_0y0,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和z0k = 0equal_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个平面P1=[A1,B1,C1,−1],P2=[A2,B2,C2,−1],P3=[A3,B3,C3,−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], 求它们的交点Po=[x,y,z]P_o=[x, y, z]Po=[x,y,z]
那么有:
PQ=(111)PQ = \begin{pmatrix} 1 \\ 1 \\ 1\\ \end{pmatrix} PQ=⎝⎛111⎠⎞
其中,
P=(A1B1C1A2B2C2A3B3C3),Q=(xyz)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⎠⎞
很容易求得上述线性方程的解QQQ, 如下:
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|/2planes_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 - ABD = D - B# 四边形面积 == 对角线向量叉乘的模的二分之一return np.linalg.norm(np.cross(AC, BD), axis=1, keepdims=True) / 2
平面方程
前面引入了平面的边界点, 那么问题又来了? 怎么通过一个平面, 去求它的4个边界点?这个问题在几何中是不可解的. 因为几何中的平面本身是无限的, 不存在边界点.
因此实际应用中往往是反过来的. 通过4个边界点, 去求这4点确定的平面的方程.(其实3个点就够了)
假设已知三个点P1=[x1,y1,z1],P2=[x2,y2,z2],P3=[x3,y3,z3]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], 求它们所在的平面Po=[A,B,C,−1]P_o=[A, B, C, -1]Po=[A,B,C,−1].
即求解线性方程组
PQ=(111)PQ = \begin{pmatrix} 1 \\ 1 \\ 1\\ \end{pmatrix} PQ=⎝⎛111⎠⎞
其中,
P=(x1y1z1x2y2z2x3y3z3),Q=(ABC)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 = 0PQ = 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)
三维空间 点线面解析相关推荐
- 三维空间——点线面关系
最基础最重要的概念--叉积,说到叉积就要聊聊行列式. 行列式的代数意义与Cramer法则联系密切,先来个简单的例子, 消除x2得到这样的结果: . 行列式 正是那个分母,其计算和叉积一样. 行列 ...
- JAVA 解析 DXF 文件 点线面圆
一.DXF 文件简介 1.人肉解析 观察几个具有代表性的 dxf 文件,点.文本.线. 使用文本工具直接打开 DXF 文件,可以看到很多字段,这里根据官方文档找规律,找到具有代表性的一些字段如下: 点 ...
- 三维空间长度温度数量_风电叶片模具水循环温度控制机及其智能化控制解析
由于风电叶片基质材料的特殊性,因此对其固化过程的温度变量有严格的要求.传统的叶片表面电阻丝固化加热方式具有动态性能差.跟踪性能低和大滞后等缺点,极易造成叶片基材失效. 针对该问题特别有设计了一种基于P ...
- 三维场景图：用于统一语义、三维空间和相机的结构
三维场景图:用于统一语义.三维空间和相机的结构 3D Scene Graph: A structure for unified semantics, 3D space, and camera 论文链接: ...
- 视觉SLAM开源算法ORB-SLAM3 原理与代码解析
来源:深蓝学院,文稿整理者:何常鑫,审核&修改:刘国庆 本文总结于上交感知与导航研究所科研助理--刘国庆关于[视觉SLAM开源算法ORB-SLAM3 原理与代码解析]的公开课. ORB-SLA ...
- 一文详解LOAM-SLAM原理深度解析
作者丨yc zhang@知乎 来源丨https://zhuanlan.zhihu.com/p/111388877 编辑丨3D视觉工坊 LOAM作为比较古老的激光匹配slam方法,一直以来都霸占着KIT ...
- 卷积神经网络（CNN）数学原理解析
来源:图灵人工智能 作者:Piotr Skalski 编辑:python数据科学 原标题:Gentle Dive into Math Behind Convolutional Neural Netwo ...
- 中美首份8000字长文解析全球热点脑机接口（重磅干货）
来源:硅谷密探 摘要:"我们所想象的一切,都会变为现实." 如果说当今什么技术最接近科幻,那么一定是脑机接口. 脑机接口的研究已经实现了意识打字(1分钟之内平均输入39个字母),还 ...
- 傅立叶变换物理意义解析进阶
1.为什么要进行傅里叶变换,其物理意义是什么? 傅立叶变换是数字信号处理领域一种很重要的算法.要知道傅立叶变换算法的意义,首先要了解傅立叶原理的意义.傅立叶原理表明:任何连续测量的时序或信号,都可以表 ...
最新文章
- 【Xamarin笔记】Events, Protocols and Delegates
- 李开复：AlphaGo 若打败了世界冠军，意味着什么？
- 程序设计入门-C语言基础知识-翁恺-第六周：数组-详细笔记（六）
- 互联网日报 | 6月7日 星期一 | 华为已捐献鸿蒙全部基础能力；芝麻信用7年免押金4000亿；奈雪的茶通过港交所上市聆讯...
- aws rds监控慢sql_使用AWS Database迁移服务进行AWS RDS SQL Server迁移
- 剑指offer(C++)-JZ52：两个链表的第一个公共结点(数据结构-链表)
- c#中的线程Thread
- 使用第三方框架解耦的一种思路—简单工厂模式的运用
- 顺网服务器ip修改工具,一键更换IP工具，修改IP地址 — 活动撸羊毛必备
- 推荐几个高质量图片网站,再也不怕没图装X了
- 笔记本电脑外接显示器
- python基因差异分析_Biopython基因组分析
- PyCharm添加Anaconda中的虚拟环境，Python解释器出现Conda executable is not found（解决方案）
- 系统吞吐量、QPS、并发数、响应时间，以及提高吞吐量的思路
- 杭电1856——并差集
- 手机生产日期查询方法
- 华为/思科已知一个ip查对应mac和交换机接口
- 解释瑞利分布的平方、莱斯分布的平方、高斯分布的平方 服从什么分布？
- 黑马程序员-Java基础加强之枚举
- 第二十八章 Caché 命令大全 TSTART 命令