问四面体(scipy.Delaunay)三维点云的截面

EN

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)