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

Hi @LuisCastilloV98,

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,
and insert instead:

boundary_points= points3d[boundary_faces]
x1, y1, z1 = boundary_points.T