Python，四面体 (scipy.Delaunay) 3D 点云的一部分 [英] Python, section of a tetrahedralized (scipy.Delaunay) 3D cloud of points

问题描述

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

I want to draw a "cross section" of a hull in a 3D space, the intersection of the hull with a plane.

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

The space is defined by axis X, Y, Z, and the crossing plane, parallel to XZ is defined by Y = 50

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

First, I loaded a cloud of 3D points in a np.array :

#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

  并且每个点在 3D 空间中由 X, Y, Z

  and each points is defined in the 3D space by X, Y, Z

  几个例子:

  ['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 四面体化:

  Now, I made a scipy.Delaunay tetrahedralization of the points:

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

  所以我可以得到所有的顶点(即船体的每个奇异四面体):

  So I can get all the vertices (i.e each singular tetrahedron of the hull):

  # indices of vertices
  indices = tri.simplices
  
  # the vertices for each tetrahedron
  vertices = points[indices]
  
  print vertices
  

  我的问题:从这里开始，我有了顶点，如何找到飞机和船体之间的所有交点?

  My question: from here, I have the vertices, how do I find all the points of intersection between the plane and the hull?

  谢谢

  推荐答案

  下面我给出了 python 代码，给定一组 3d 点和一个平面(由其法向量和平面上的一个点定义)计算 3d Delaunay三角剖分(曲面细分)和 Delaunay 边与平面的交点.

  Below I give python code that, given a set of 3d points and a plane (defined by its normal vector and a point on the plane) computes the 3d Delaunay triangulation (tessellation) and the intersection points of the Delaunay edges with the plane.

  下图以单位立方体中与 x=0 平面相交的 20 个随机点为例(相交点为蓝色)显示了结果.用于可视化的代码是从

  The following figure visualizes the result on an example of twenty random points in the unit cube intersecting with the x=0 plane (the intersection points are colored blue). The code used for the visualization is modified from the code in this answer.

  要实际计算平面交点，我使用以下代码.基本函数plane_delaunay_intersection 使用了两个辅助函数——collect_edges 来收集Delaunay 三角剖分的边(每个线段只有一个副本)，以及plane_seg_intersection, 线段与平面相交.

  To actually compute the plane intersection points I use the following code. The basic function plane_delaunay_intersection uses two auxiliary functions - collect_edges to gather the edges of the Delaunay triangulation (only one copy of each segment), and plane_seg_intersection, which intersects a line segment with a plane.

  代码如下:

  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
  

  以下代码用于生成上例图的输入:

  The following code was used to generate the input for the example figure above:

  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)
  

  这篇关于Python，四面体 (scipy.Delaunay) 3D 点云的一部分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！