Meshes¶
Python imports and toolbox¶
import pygel3d
from pygel3d import hmesh, jupyter_display as jd
import numpy as np
import plotly.graph_objects as go
jd.set_export_mode(True)
Mesh primitives¶
A mesh is made of a collection of vertices and faces. Vertices are associated with 3D coordinates. This is the geometry of the mesh. The faces connect vertices together. This is the topology of the mesh. Mathematically, the topology is the set of open subsets of a set. Intuitively, an open subset can be seen in terms of neighborhoods : a neighborhood of a point is roughly an open subset containing it. By defining open subsets, the topology therefore defines which points are neighbors. When connecting the vertices with faces, these faces define which vertices are neighbors.
vertices = [
[0,0,0], #0
[0,1,0], #1
[0,0,1], #2
[0,1,1], #3
[1,0,0], #4
[1,1,0], #5
[1,0,1], #6
[1,1,1] #7
]
faces = [
[0,3,1], [0,2,3], #x=0
[4,5,7], [4,7,6], #x=1
[0,6,2], [0,4,6], #y=0
[1,3,7], [1,7,5], #y=1
[0,1,5], [0,5,4], #z=0
[2,6,7], [2,7,3], #z=1
]
cube = hmesh.Manifold.from_triangles(vertices, faces)
jd.display(cube, smooth=False)
Many mesh file formats will describe meshes in this way, with a list of vertices and a list of faces. Below the example of the "Wavefront obj" file format.
!cat suzanne_triangles.obj
suzanne = hmesh.load("suzanne_triangles.obj")
jd.display(suzanne)
Triangles are a simple primitive, but many mesh data structures will support faces with a higher valence
bigguy = hmesh.load("bigguy.obj")
jd.display(bigguy, smooth = False)
Note however that other types of faces are usually converted to triangles when it comes to rendering and shading, as more than three points will generally not lie in the same plane.
Mesh circulation with arrays¶
Circulating on a mesh is a classical operation required in many situations. Desired operations are :
- givent a face, list its vertices
- given a face, list the adjacent faces
- given a vertex, list the adjacent faces
- givent a vertex, list the neighboring vertices
However, just storing the mesh as an array of vertices and an array of faces yields inefficient operations.
def face_neighbors(faces, f):
verts = set(faces[f])
neighbors = []
for nf,nface in enumerate(faces):
nverts = set(nface)
if len(verts.intersection(nverts))==2:
neighbors.append(nf)
return neighbors
faces
[[0, 3, 1], [0, 2, 3], [4, 5, 7], [4, 7, 6], [0, 6, 2], [0, 4, 6], [1, 3, 7], [1, 7, 5], [0, 1, 5], [0, 5, 4], [2, 6, 7], [2, 7, 3]]
face_neighbors(faces, 0)
[1, 6, 8]
Note however that when the operation is global, the overall cost can often be mitigated as a loop over the whole mesh is done anyway, and depending on the operation, it is often possible to do everything in a single pass.
def global_neighbors(faces):
edges = {}
neighbors = []
for fi,f in enumerate(faces):
neighbors.append([])
for i in range(len(f)):
edge = tuple(sorted((f[i], f[(i+1)%len(f)]))) ;
if edge in edges:
nf = edges[edge]
neighbors[nf].append(fi)
neighbors[fi].append(nf)
else:
edges[edge] = fi
return neighbors
global_neighbors(faces)
[[1, 6, 8], [0, 4, 11], [3, 7, 9], [2, 5, 10], [1, 5, 10], [3, 4, 9], [0, 7, 11], [6, 2, 8], [0, 7, 9], [8, 2, 5], [4, 3, 11], [10, 6, 1]]
When such neighbors are used intensively, these face neighbors are commonly stored as an additional array with a convention to easily determine the neighboring face along a given edge. Although such neighbors improve the mesh circulation, these are not enough to provide efficiently all the desired neighborhoods. For instance, listing the faces or vertices adjacent to a vertex is still a challenge. When full circulation is desired, the classical data structure is a halfedge data structure.
Halfedge data structure¶
A halfedge data structure stores the mesh from its edges. For each edge of the mesh, two half edges are created, one for each adjacent face. Halfdges are oriented such that :
- no two halfedges in the same face point towards the same vertex
- opposite halfedges on both sides of an edge are oriented in the opposite direction
To circulate on the mesh, each half edge stores :
- a reference to its face
- a reference to the vertex at its tip ;
- a reference to the next halfedge along the face, starting from the vertex it points to ;
- a reference to the opposite halfedge along the same edge ;
- optionally a reference to the previous halfedge along the face
Equipped with this data structure it is then possible to circulate efficiently on the mesh. The data structure used by the module we use here is a halfedge data structure.
h = cube.halfedges()[0]
print(cube.incident_face(h))
print(cube.incident_vertex(h))
print(cube.next_halfedge(h))
print(cube.opposite_halfedge(h))
0 1 1 5
Our module implements circulation wrappers around these halfedges :
print("faces neighboring vertex 0 : ", list(cube.circulate_vertex(0,'f')))
print("vertices neighboring vertex 0 : ", list(cube.circulate_vertex(0,'v')))
print("halfedges emitted from vertex 0 : ", list(cube.circulate_vertex(0,'h')))
h = cube.circulate_vertex(0, 'h')[0]
print("vertex neighboring vertex 0 from halfedge ", h, " is ", cube.incident_vertex(h))
faces neighboring vertex 0 : [9, 5, 4, 1, 0, 8] vertices neighboring vertex 0 : [5, 4, 6, 3, 1, 2] halfedges emitted from vertex 0 : [27, 15, 12, 3, 0, 24] vertex neighboring vertex 0 from halfedge 27 is 5
These circulation wrappers can be implemented just using halfedge circulation
def vertex_neighbors(mesh, v):
#there is no simpler method to get a halfedge emitted at v
h = mesh.circulate_vertex(v, 'h')[0]
neighbors = []
start = h
while True:
neighbors.append(mesh.incident_vertex(h))
h = mesh.next_halfedge(mesh.opposite_halfedge(h))
if h == start:
break
return neighbors
vertex_neighbors(cube, 0)
[5, 2, 1, 3, 6, 4]
Applications¶
First, this data structure offers the same capabilities as the array data structure. Its only problem is that it is a bit more memory consuming, and editing the mesh may require more maintenance to ensure that the halfedges remain properly connected. For instance, computing the volume of the mesh is as simple as :
def mesh_volume(mesh):
volume = 0
mat = np.zeros((3,3))
for f in mesh.faces():
for i,v in enumerate(mesh.circulate_face(f)):
mat[i,:] = mesh.positions()[v]
volume += np.linalg.det(mat) / 6
return volume
print("the cube has volume : ", mesh_volume(cube))
print("suzanne has volume : ", mesh_volume(suzanne))
the cube has volume : 0.9999999999999999 suzanne has volume : 0.6366285731634208
Propagating along the mesh is now also possible.
from queue import SimpleQueue as Queue
def bfs(mesh, start):
Q = Queue()
visited = [False for _ in mesh.vertices()]
distance = [len(visited) for _ in mesh.vertices()]
visited[start] = True
distance[start] = 0
Q.put(start)
while not Q.empty():
v = Q.get()
for nv in mesh.circulate_vertex(v, 'v'):
if not visited[nv]:
visited[nv] = True
distance[nv] = distance[v] + 1
Q.put(nv)
return distance
colors = bfs(suzanne, 0)
jd.display(suzanne, data = colors)
Limitations¶
Apart from being a bit memory consuming, halfedge data structures cannot represent any surface. First, the surface has to be manifold. Mathematically this means each point on the surface has a neighborhood which can be bijectively and continuously (in both directions) deformed into a disk. Boundaries are possible by reducing this conditions to a half disc. This is however not enough. Circulating on half edges requires halfedges to be consistently oriented across the whole mesh. In particular, this means that opposite halfedges have to be oriented in opposite directions, and no two halfedge of the same face points to the same vertex. This not possible on any manifold. A Moëbius strip, a Klein bottle or a Projective plane are classical examples of non orientable manifolds.