Python，四面体（scipy.Delaunay）3D点云的一部分答案

【问题标题】：Python, section of a tetrahedralized (scipy.Delaunay) 3D cloud of pointsPython，四面体（scipy.Delaunay）3D点云的一部分
【发布时间】：2014-06-11 12:42:30
【问题描述】：

我想在 3D 空间中绘制船体的“横截面”，船体与平面的交点。

空间由轴X, Y, Z定义，交叉平面平行于XZ由Y = 50定义

首先，我在np.array 中加载了points 的3D 云：

#loading colors

points = np.array([(GEO.XYZRGB(rank, name, X, Y, Z))
                  for rank, name, X, Y, Z in csv.reader(open('colors.csv'))])

点结构为rank, name, X, Y, Z, R, G, B
每个点都由X, Y, Z在3D空间中定义

几个例子：

['2' 'Y2    ' '77.89506204' '87.46909733' '42.72168896' '254' '244' '21']
['3' 'Y4    ' '76.95634543' '83.94271933' '39.48573173' '255' '234' '0']
['4' 'PINKwA' '64.93353667' '59.00840333' '84.71839733' '218' '154' '225']
...

现在，我对点进行了scipy.Delaunay 的四面体化：

# Delaunay triangulation    
tri = scipy.spatial.Delaunay(points[:,[2,3,4]], furthest_site=False)

所以我可以得到所有的vertices（即船体的每个奇异四面体）：

# indices of vertices
indices = tri.simplices

# the vertices for each tetrahedron
vertices = points[indices]

print vertices

我的问题：从这里，我有顶点，我如何找到飞机和船体之间的所有交点？

谢谢

【问题讨论】：

标签： python numpy 3d delaunay

【解决方案1】：

下面我给出 python 代码，给定一组 3d 点和一个平面（由其法线向量和平面上的一个点定义）计算 3d Delaunay 三角剖分（镶嵌）和 Delaunay 边与飞机。

下图显示了单位立方体中与x=0 平面相交的二十个随机点示例的结果（相交点为蓝色）。用于可视化的代码修改自代码in this answer。

要实际计算平面交点，我使用以下代码。基本函数 plane_delaunay_intersection 使用两个辅助函数 - collect_edges 收集 Delaunay 三角剖分的边缘（每条线段只有一个副本），以及 plane_seg_intersection 将线段与平面相交。

代码如下：

from scipy.spatial import Delaunay
import numpy as np

def plane_delaunay_intersection(pts, pln_pt, pln_normal):
    """ 
    Returns the 3d Delaunay triangulation tri of pts and an array of nx3 points that are the intersection
    of tri with the plane defined by the point pln_pt and the normal vector pln_normal.
    """
    tri = Delaunay(points)
    edges = collect_edges(tri)
    res_lst = []
    for (i,j) in edges:
        p0 = pts[i,:]
        p1 = pts[j,:]
        p = plane_seg_intersection(pln_pt, pln_normal, p0, p1)
        if not np.any(np.isnan(p)):
            res_lst.append(p)
    res = np.vstack(res_lst)
    return res, tri 


def collect_edges(tri):
    edges = set()

    def sorted_tuple(a,b):
        return (a,b) if a < b else (b,a)
    # Add edges of tetrahedron (sorted so we don't add an edge twice, even if it comes in reverse order).
    for (i0, i1, i2, i3) in tri.simplices:
        edges.add(sorted_tuple(i0,i1))
        edges.add(sorted_tuple(i0,i2))
        edges.add(sorted_tuple(i0,i3))
        edges.add(sorted_tuple(i1,i2))
        edges.add(sorted_tuple(i1,i3))
        edges.add(sorted_tuple(i2,i3))
    return edges


def plane_seg_intersection(pln_pt, pln_normal, p0, p1):
    t0 = np.dot(p0 - pln_pt, pln_normal)
    t1 = np.dot(p1 - pln_pt, pln_normal)
    if t0*t1 > 0.0:
        return np.array([np.nan, np.nan, np.nan])  # both points on same side of plane

    # Interpolate the points to get the intersection point p.
    denom = (np.abs(t0) + np.abs(t1))
    p = p0 * (np.abs(t1) / denom) + p1 * (np.abs(t0) / denom)
    return p

以下代码用于为上图示例生成输入：

np.random.seed(0)
x = 2.0 * np.random.rand(20) - 1.0
y = 2.0 * np.random.rand(20) - 1.0
z = 2.0 * np.random.rand(20) - 1.0

points = np.vstack([x, y, z]).T
pln_pt = np.array([0,0,0])  # point on plane
pln_normal = np.array([1,0,0])  # normal to plane
inter_pts, tri = plane_delaunay_intersection(points, pln_pt, pln_normal)

【讨论】：