I’m not sure if it executes the code wrong, but when I execute it the points of “boundary_faces” are on another scale than the original figure.
import numpy as np
import plotly.graph_objs as go
import pyvista as pv #pip install pyvista
def pcloud(points3d, marker_size=1.5, marker_color='#454F8C'):
#define the trace representing the point cloud
points3d=np.asarray(points3d)
if points3d.ndim != 2 or points3d.shape[1] != 3:
raise ValueError('your data is not a 3D point cloud')
return go.Scatter3d(
name='',
mode='markers',
x=points3d[:,0],
y=points3d[:,1],
z=points3d[:,2],
marker_size=marker_size,
marker_color=marker_color
)
def get_mesh(points3d, faces, color='lightblue', opacity=1):
# define the Mesh3d representing the alpha-shape
points3d=np.asarray(points3d)
if points3d.ndim != 2 or points3d.shape[1] != 3:
raise ValueError('your data is not a 3D point cloud')
faces = np.asarray(faces)
i, j, k = np.asarray(faces).T
return go.Mesh3d(
color=color,
opacity=opacity,
i=i, j=j, k=k,
x =points3d[:,0],
y= points3d[:,1],
z= points3d[:,2],
flatshading=True
)
def alphashape3d(points3d, alpha=1):
#extract 0, 1, 2, 3-simplices of the alpha shape constructed from points3d
# Here alpha =1/alphahull, where alphahull is a property of a Plotly alpha shape
cloud = pv.PolyData(points3d)
mesh = cloud.delaunay_3d(alpha=alpha)
unconnected_points3d = [] #isolated 0-simplices
edges = [] # isolated edges, 1-simplices
faces = [] # triangles that are not faces of some tetrahedra
tetrahedra = [] # 3-simplices
for k in mesh.offset:
length = mesh.cells[k]
if length == 2:
edges.append(list(mesh.cells[k+1: k+length+1]))
elif length ==3:
faces.append(list(mesh.cells[k+1: k+length+1]))
elif length == 4:
tetrahedra.append(list(mesh.cells[k+1: k+length+1]))
elif length == 1:
unconnected_points3d.append(mesh.cells[k+1])
else:
raise ValueError('For a 3d point cloud only unconneted points3d, edges,\
triangles and tetrahedra are set up')
return unconnected_points3d, edges, faces, tetrahedra
def get_tetrahedra_faces(tetrahedra, faces):
# extract tetrahedra faces and append them to the existing faces
for t in tetrahedra:
faces.extend([[t[0], t[1], t[2]],
[t[0], t[2], t[3]],
[t[0], t[3], t[1]],
[t[1], t[2], t[3]]])
return faces
points3d = np.loadtxt(np.DataSource().open('https://raw.githubusercontent.com/empet/Plotly-plots/master/Data/data-file.txt'))
x, y, z = points3d.T
scatt = pcloud(points3d, marker_color='#454F8C', marker_size=3)
cloud = pv.PolyData(points3d) # set up the pyvista point cloud structure
#alphaparameter for alpha shape definition
alpha=5
#Extract the total simplicial structure of the alpha shape defined via pyvista
mesh = cloud.delaunay_3d(alpha=alpha)
#and select its simplexes of dimension 0, 1, 2, 3:
unconnected_points3d = [] #isolated 0-simplices
edges = [] # isolated edges, 1-simplices
faces = [] # triangles that are not faces of some tetrahedra
tetrahedra = [] # 3-simplices
for k in mesh.offset: #HERE WE CAN ACCESS mesh.offset
length = mesh.cells[k]
if length == 2:
edges.append(list(mesh.cells[k+1: k+length+1]))
elif length ==3:
faces.append(list(mesh.cells[k+1: k+length+1]))
elif length == 4:
tetrahedra.append(list(mesh.cells[k+1: k+length+1]))
elif length == 1:
unconnected_points3d.append(mesh.cells[k+1])
else:
raise ValueError('For a 3d point cloud only unconneted points3d, edges,\
triangles and tetrahedra are set up')
#get tetrahedra faces and append them to the previous list of isolated faces
faces = get_tetrahedra_faces(tetrahedra, faces)
alphashape = get_mesh(points3d, faces, opacity=1) #OR SET A LOWER OPACITY!!!!
boundary_mesh = mesh.extract_geometry()
boundary_faces = boundary_mesh.faces.reshape((-1,4))[:, 1:]
my_mesh = get_mesh(points3d,boundary_faces, opacity=1)
x1, y1, z1 = boundary_faces.T
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, c='r', marker='o')
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.show()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x1, y1, z1, c='r', marker='o')
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.show()