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.
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.
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.
Here is the code:
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)
