Retrieving indices from Mesh3D alpha shapes

Hello all,

I’m following this tutorial on finding the alpha shapes of a 3D point cloud in order to create a non-convex hull of the point cloud. I have successfully implemented the tutorial. However, I would now like to save the alpha shape data for post-processing in other software.

Here’s the snippet from the tutorial where the mesh3d object with alpha shapes is created:

``````simplexes = Mesh3d(alphahull =10.0,
name = '',
x =x,
y= y,
z= z
)
``````

Is it possible to somehow get inside that simplexes object so that I can access the indices of the points that make-up the alpha shapes?

So far I have found this post which seems to indicate what I’m asking isn’t possible at the time it was written. If so, is there some other way to find the 3D alpha shapes of a point-cloud so that I can save them? It’s actually not necessary that I plot the data. I’m just interested in finding the non-convex boundary of a 3D point-cloud and saving the indices/mesh of that boundary.

Thanks in advance for any feedback.

Hi all,

I too am interested in finding the indices of the points on the mesh boundary. I was wondering if there was any movement on this?

nwp: did you find another solution to this problem?

Thanks

@arthurPy

Plotly Mesh3d defines and displays an alpha-shape, but cannot return its simplicial structure.
To get it you can use `pyvista` (`pip install pyvista`), that can call the function `Delaunay_3d` with a prescribed
`alpha` (pyvista alpha = 1/alphahull; alphahull is an attribute of the Plotly Mesh3d). This function returns a mesh
and from `mesh.cells` and `mesh.offset` you can extract the simplices of dimension 0, 1, 2, 3 that define the simplicial structure of the alpha shape. In this case you can plot the corresponding alpha shape as a Mesh3d from points and triangles (i.e. providing x, y, z, i, j, k to define a an instance of go.Mesh3d). Example https://plot.ly/~empet/15549/.

Thanks @empet for your reply! And thanks for the example. I have been playing with this notebook, but I’m not sure I fully understand.

What I am looking for are the indices of the points on the boundary of the alphahull.

For example. If we define points at the vertices of a cube and one additional point (the 0th entry of pts) in the interior of the cube:
`pts = np.asarray([[15, 15, 15], [10, 10, 10], [20, 10, 10], [20, 20, 10], [10, 20, 10], [10, 10, 20], [10, 20, 20], [20, 20, 20], [10, 10, 20], [20, 10, 20]])`

Is it possible to define a function that takes the mesh as an input and returns the indices of the points on boundary:
`exterior_indices(mesh) = [1,2,3,4,5,6,7]`
and the indices of the points not on the boundary `interior_indices(mesh) = [0]`

Possibly this is what the function `edge_trace` is for in your example?
Also when I run `mesh.cells` or `mesh,offset` I get an `AttributeError`.

Thanks for the help.

@arthurPy

The function `alphashape3d()` in the posted Jupyter notebook defines an alpha shape calling the pyvista `Delaunay_3d` with a prescribed parameter alpha.
The lines of code inside this function body process the output of the `Delaunay_3d`.

To be able to inspect the contents of the associated elements (`mesh.offset`, for example), I’ll paste the lines of code in `alphashape3d()` in a new version below, as a main part of code. Hence to run it copy, please, all the functions in that notebook, and as a main part of code the following
lines:

``````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=0.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!!!!

fig1 = go.Figure(data =[scatt, alphashape])
title = f'Alpha Shape of a 3D Point Cloud<br>alpha={alpha}'

fig1.update_layout(title_text=title, title_x=0.5, width=800, height=800)
fig1.show()
``````

fig1 with an opacity less than 1 in Mesh3d displays the general structure of the generated alpha complex, with
2129 triangular faces:

``````np.array(faces).shape
(2129, 3)
``````

The following lines of code extract the faces corresponding to the alpha complex boundary, i.e. what you need:

``````boundary_mesh = mesh.extract_geometry()
boundary_faces = boundary_mesh.faces.reshape((-1,4))[:, 1:]

my_mesh = get_mesh(points3d,boundary_faces, opacity=1)
figb = go.Figure(data=[scatt, my_mesh])
figb.update_layout(title_text=title, title_x=0.5, width=800, height=800)
figb.show()
``````

Check the number of `boundary_faces`:

``````boundary_faces.shape
(509, 3)
``````

Hence the number of boundary faces is 4x less than the number of alpha complex faces.

Note that:

1. the function ‘edge_trace()’ can be used to display the isolated edges i.e. the edges that are not triangle sides in faces.
You can call it when the edges extracted from the `Delaunay_3d` is not an empty list.
2. The pyvista data structures are a bit different from the usual 3d meshes, read from a Wavefront `obj` file, `ply`, `off` or `stl` format files .

To see what information you can read from `boundary_mesh`
just type

``````help(boundary_mesh)
``````

Thanks so much for this.

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],
)

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

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.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.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()

``````

Welcome to Plotly forum!
It seems that it’s not clear how a mesh is defined. It consists in two arrays of shape (nr_points, 3), respectively
(nr_faces, 3) .The former is an array of floats, representing the mesh points, while the latter, an array of ints, representing the mesh faces (triangles).
Hence in your code `boundary_faces` is an array of ints. For example, displaying the first 5 rows we get:

``````print(boundary_faces)

array([[139,  49, 138],
[ 66,  65,  24],
[ 65,  24, 138],
[ 24, 138,  66],
[ 63,  51,  57]], dtype=int64)
``````

The first row tell us that the first boundary face has as vertices the points `points3d[139], points3d[49], points3d[138]`. Now its clear why you got plots on different scales: you did not use the coordinates of points but the point indices.
To get the right plot, just comment out the line of code `x1, y1, z1 = boundary_faces.T`,
``````boundary_points= points3d[boundary_faces]