Introduction to PyGEL

_J. Andreas Bærentzen, January 2018 (revised February 2021) _

The PyGEL Python library for geometry processing provides Python3 bindings for a subset of the features of the C++ based GEL library.

PyGEL is a Python package consisting of five modules:

• hmesh provides Manifold which is a class that represents polygonal meshes using the halfedge representation. hmesh also provides a slew of functions for manipulating polygonal meshes and the MeshDistance class which makes it simple to compute the distance to a triangle mesh. hmesh also includes the skeleton_to_feq function which makes it straight forward to produce a quad-only mesh from a graph.
• graph contains the Graph class which is used for graphs: i.e. collections of vertices (in 3D) connected by edges. Unlike a Manifold, a Graph does not have to represent a surface. There are also some associated functions which may be useful: in particular, there is the LS_skeletonize and MSLS_skeleton functions which compute a curve skeleton from a Graph and returns the result as a new Graph. The MSLS_skeleton function is multi-scale and hence much faster.
• gl_display provides the Viewer class which makes it simple to visualize meshes and graphs.
• jupyter_display makes it easy to use PyGEL in the context of a Jupyter Notebook. This module contains a function that allows you to create a widget for interactively visualizing a mesh or a graph in a Notebook. The feature is based on the Plotly library and it is possible to export the resulting notebooks to HTML while preserving the interactive 3D graphics in the notebook (as we will see below).
• spatial contains the I3DTree class which is simply a kD-tree specialized for mapping 3D points to integers - typically indices. Of course, scipy.spatial has a more generic class, so this is perhaps not the most important part of PyGEL.

This text provides a brief introduction to PyGEL. If you are planning to use PyGEL, please also cast a glance at the PyGEL Reference Documentation. All features of the API are concisely described in the reference.

PyGEL Installation

PyGEL is on PyPi under the name PyGEL3D. Turned out PyGEL was taken by another project which I am sure is excellent and does something very different, so we still refer to PyGEL as PyGEL, but for the purpose of installation from PyPi, the package is called PyGEL3D. Maybe all you need to do is

> pip install PyGEL3D

Now, you may have to use pip3 (because Python3) and it might be necessary to install with --user. However, it also depends on which platform you are on. GEL and PyGEL are compiled automatically for Windows, MacOS and Linux (specifically Ubuntu 18.04) using GitHub actions every time we push to GitHub. However, it is hard to make sure that the resulting binaries work without fail on every computer, so you may be out of luck with PyPI.

In that case, GEL and PyGEL can easily be compiled using CMake. Clone the github repository and be on your (hopefully merry) way. Instructions are in the README.

Python Distribution and Required Packages

It is recommended that you either use the standard distribution from Python.org or the Anaconda distribution. Anaconda may make it a little easier to install the packages that you need and it also has the advantage that you can install everything as a user with no administrator priviliges. That being said, the regular Python distribution should also work well for PyGEL. If you use the standard distribution, Pip is the tool you would normally use to install new components.

Regardless of distribution, you need to have Plotly (for Python) installed. However, normally you don't have to worry because if you install either from PyPi or a version you built yourself it will ultimately be done with pip and this takes care of dependencies. If you plan to use PyGEL from a Jupyter notebook, you clearly need Jupyter notebook.

PyGEL: Rationale and Design

The GEL library was created by the author (Andreas Bærentzen) with significant contributions from colleagues at the Section for Visual Computing of the Computer Science and Applied Mathematics Department at The Technical University of Denmark. The library has been used for several research projects and as the basis for exercises in a long running course on Geometry Processing. The GEL library is also the recommended platform for solving problems from our book "Guide to Computational Geometry Processing". However, GEL is a C++ library and while that has not been a big limitation in the past, we wanted to offer a geometry processing course in a more general form -- and to undergraduates. In this context, it seemed prudent to switch to Python. While there are several other geometry processing libraries with Python bindings, it turned out to be surprisingly straight forward to use ctypes to develop Python bindings for GEL, and since none of the other options seemed to be a drop in replacement, we decided to create PyGEL.

The scope of PyGEL is a little narrower than the scope of GEL. For instance, PyGEL does not contain voxel grids or a vector matrix library for 2,3, and 4D vectors and matrices. To a large extent this is because C++ libraries often work best if they have few dependencies whereas this appears to be a bit less of a concern with Python where most of the packages we need can be installed with ease using e.g. pip. If you happen to have a strong interest in some of the choices that led to how PyGEL was developed, the process is documented here

Perhaps, the pivotal feature of GEL, and by extension PyGEL, can be said to be HMesh - the polygonal mesh representation. HMesh is halfedge based. Originally, the halfedge representation was chosen since it is more efficient than the older winged edge representation while also making connectivity information readily available. Furthermore, the halfedge representation makes it trivial to support polygons with arbitrary numbers of edges. From an implementation point of view, many operations are probably harder to implement than for simpler mesh representations, but at this point, the GEL code base seems fairly robust while offering a rich set of mesh manipulation operations.

An important concern (for teaching) was to be able to run PyGEL inside a Jupyter Notebook. It is probably best to use a programming editor or an IDE to develop code, but for presentation purposes, the ability provided by JN to mix code, output from code, and text is very hard to beat. One big stumbling block was that we wanted students to be able to show 3D models inside a Jupyter Notebook. Several frameworks ostensibly make this possible, but because of subtle bugs almost none of these do so in a way that allows the 3D model to survive into a static HTML page generated from the Notebook. One exception turned out to be Plotly. HTML pages generated from ipynb files using the Nbconvert program would actually retain interactive views of 3D models. Thus, the PyGEL package contains a module called jupyter_display with a single function, display, that converts a mesh into an interactive view, embedded in the notebook. One - absolutely maddening - caveat, though, is that LaTeX formulas in the Markdown text remain LaTeX code when you do that. There is some bug which causes MathJax to be unintentionally neutralized by Plotly. Hopefully, that will be resolved, but for now the bug remains.

Loading and viewing a mesh

To use PyGEL you need to import the appropriate module. The code below starts by doing that, then it loads a mesh (i.e. a Manifold) called m, creates a Viewer called viewer, and, finally, it displays m using viewer. Note that the script below could be run from a Jupyter notebook but also whatever editor you might prefer or from an interactive Python shell.

from pygel3d import hmesh, gl_display as gl

m = hmesh.load("../data/bunnygtest.obj")
viewer = gl.Viewer()
viewer.display(m)

Having executed the code above, you should see a window displaying a 3D model of a bunny, and it is now possible to play with the mesh. You can rotate by pressing the left mouse button and dragging the mouse. You can zoom using the right mouse button, and if you hold shift the right mouse button pans instead of zooming. One thing that might be puzzling is that the display function does not return. That is because the viewer is in the same thread of execution as the rest of your script: it is not running in parallel. However, if you want to return the script, all you have to do is to press ESC inside the window. Doing so, you might notice that the image freezes but the window does not go away? That is still entirely as intended. We might want to return to the window after all - either to visualize a different mesh or simply to look a bit more at the one we have. To make the window active, just call event_loop:

viewer.event_loop()

You will have noticed that the window with the bunny came alive again. Hit ESC one more time to exit the viewer. If you really want the window to go away, you can either wait for the entire script to terminate or explicitly delete the object like so:

del viewer

Exploring the Viewer

The display function in Viewer has a number of parameters which were not involved in the simple example above. In the example, display was called only with the mesh m. Instead of m we could have called display with a Graph. The Viewer is happy to display both types of object even at the same time. However, it is much more full fledged as a mesh viewer than a graph viewer.

• mode : a single character that determines how the mesh is visualized:
  • 'w' - wireframe, which is the default visualization mode because we so often want to see the polygons.
  • 'i' - isophote rendering. This will show some funny lines which are lines on the surface such that there is the same angle between the normal and the light source. In other words, isophotes are curves of even intensity of shading.
  • 'g' - glazed (this is supposed to look a bit like ceramics),
  • 's' - scalar field. Using this mode, you can draw a scalar field on the surface if you provide a scalar value for each vertex.
  • 'l' - line field. This is the same sort of visualization as "scalar field" but you need to provide a 3D vector for each vertex.
  • 'n' - normal shading. Fairly boring, but sometimes this is what you want`
  • 'x' - renders the mesh transparent in a constant color. This is probably only useful if you want to show a mesh and a graph together, but then it is very useful.
• smooth : if True (which is default) we use vertex normals. Otherwise, face normals.
• bg_col : background color. The default is dark grey [0.3,0.3,0.3]
• data : per vertex data for visualization. This is either a scalar or vector field. It is ignored unless one of these two visualization modes is selected. The default here is None.
• reset_view : if False view is as left in the previous display call. If True, the view is reset to the default. This can be useful if you make changes to a mesh and then want to return to the view you had, and the default is indeed False.
• once : if True we immediately exit the event loop and return. However, the window stays and if the event loop is called from this or any other viewer, the window will still be responsive. Default is False.

Jupyter Caveat

If you are using PyGEL from within a Jupyter Notebook you may find it absolutely intolerable to use the gl_display.Viewer. At least I find that the Jupyter widgets for visualizing meshes are far more useful in Jupyter and some of my students have insisted that the gl_display.Viewer is broken. It is not, but it can feel that way when it stops the execution of the cells in your notebook.

Visualizing in Jupyter

So let us talk about the jupyter_display module. Sometimes, we do want the ability to show a mesh in a Jupyter Notebook. To do so, we need a different visualization tool. As now mentioned a few times, such a tool has been developed. It is based on the visualization library called Plotly which works quite well in conjunction with Jupyter Notebooks. Displaying a mesh with jupyter_display is as easy as with gl_display but the features are fewer and a bit different. In the following example, jupyter_display is used to visualize a mesh using wireframe (default) and flat shading. Perhaps, it behooves me to also explain the set_export_mode function. This function must be called if you want to be able to export the notebook to HTML with the interactive visualization widgets preserved! Also, it appears that if you do not call it, display must be the last thing you do in a cell. I hate these gotchas, but it is a rather tall software stack we are working on, and it gets a bit wobbly.

from pygel3d import jupyter_display as jd
jd.set_export_mode(True)
jd.display(m, smooth=False)

while that works well, we will often need to associate data with the mesh. jupyter_display supports only scalar fields. In the example below, the vertex x coordinates are used as a scalar field.

pos = m.positions()
jd.display(m, data=pos[:,0], wireframe=False)

Unfortunately, the Plot.ly based visualization is not quite as flexible as the one based on OpenGL, so these two examples exhaust the features of jupyter_display.

Working with meshes

The point of having a mesh representation is to be able to work with the geometry. We need to be able to visit the faces, edges, and vertices of the mesh and make queries about both geometry and connectivity. As our first example, assume we want to find the smallest and the biggest x coordinate of any vertex in the mesh. The code below loops over all vertices and finds the minimum and the maximum x coordinate which are subsequently printed out.

min_x_coord = 1e32
max_x_coord = -1e32

pos = m.positions()
for v in m.vertices():
    min_x_coord = min(min_x_coord, pos[v][0])
    max_x_coord = max(max_x_coord, pos[v][0])

print(min_x_coord, max_x_coord)

-0.0947041 0.0610337

At this point, you might be asking yourself what pos and v are? Let us explain that starting with v. v is a vertex index which means that it is an integer we can use to refer to a vertex. For instance, if we have all our vertex positions in an array then v would be used to index into that array, referring to a specific vertex position.

In fact, pos is precisely an array that contains all of the vertex positions -- and it is not a copy. Any changes that you make to pos are directly reflected in the mesh. To illustrate the power of that, let us randomize all vertex positions.

from random import random
from numpy import array

pos_backup = array(pos)
for v in m.vertices():
    pos[v] = [random(), random(), random()]
jd.display(m)

pos[:] = pos_backup # so let us restore those positions while we remember ...

Note that we save the vertex positions before nuking them. pos_backup = array(pos) creates a new array with the contents of pos and assigns it to pos_backup, thereby making a copy of the positions. Afterwards, we use the slice notation pos[:] = pos_backup to copy back the vertex positions. If we had just typed pos = pos_backup then pos would simply have been another reference to the pos_backup array, but the actual vertex positions had not been changed. Of course, the last bit of information is more about Python than PyGEL, but coming from C++, it was initially confusing to me that in Python a=b does not copy the contents of b into a but simply makes a a new name for the data in b.

Moving on.

To give just one example of iteration over the faces of a mesh, consider

avg_area = 0
F = m.faces()
for f in F:
    avg_area += m.area(f)
avg_area /= len(F)
print("Average area : ", avg_area)

Average area :  8.34321070732137e-06

Above, we exploited that m.faces() is just an iterable container with all the face indices. We can do the exact same thing with halfedges, but at this point, you have probably gotten the idea.

Circulation

Very frequently, we want to visit all the vertices that are neighbors of a given vertex or all the halfedges that emanate from that vertex or all of the faces which are incident on the vertex. Such queries can be carried out using the circulate_vertex function. Given a mesh, m, and a vertex v, we can find the neighbors of v using:

N = m.circulate_vertex(v,mode='v')

The '' indicates that we want the neighbor vertices, so when the function returns, N contains an iterable sequence of the neighboring vertices. If we had given mode='f' or mode='h', we had obtained the incident vertices or emanating halfedges. mode= can be omitted, by the way, and defaults to 'v'.

In the code example below, we smooth the mesh by repeatedly assigning the average of the neighboring vertices to each vertex. This is a crude way of performing smoothing, but it works for demonstration purposes.

from numpy import zeros

for _ in range(0,50):
    avg_pos = zeros(pos.shape)
    for v in m.vertices():
        N = m.circulate_vertex(v,'v')
        for vn in N:
            avg_pos[v] += pos[vn]
        avg_pos[v] /= len(N)
    pos[:] = avg_pos
jd.display(m)

Circulation does not only work for vertices, we can also circulate around faces. For instance,

H = m.circulate_face(f, 'h')

produces an iterable sequence, H, containing all the halfedges that make up the boundary of face f. In a similar way, we can get the vertices of f and the adjacent faces.

Other mesh operations

There is much we have not covered above. The Manifold class has many member functions that operate on individual halfedges, faces, or vertices and which can be used to refine or simplify the mesh. There are also several functions outside of the Manifold class which can be used to manipulate the entire mesh. The reader is strongle encouraged to peruse the reference documentation (PyGEL Reference Documentation) to get a better overview of the features.

Below is one more simple example that shows the use of a few more functions. This time, we load the original bunny mesh, its holes are closed, the mesh is triangulated and finally simplified to 5%. We cleanup to remove vertices that are no longer used, and, finally, show the mesh in the Notebook. Also, the numbers of vertices before and after simplification are printed.

bunny = hmesh.load("../data/bunny.obj")
print("vertices before simplification :", bunny.no_allocated_vertices())
hmesh.close_holes(bunny)
hmesh.triangulate(bunny)
hmesh.quadric_simplify(bunny, 0.05)
bunny.cleanup()
print("vertices after simplification :", bunny.no_allocated_vertices())
jd.display(bunny)

vertices before simplification : 34834
vertices after simplification : 1742

Clearly, the only practical thing that has been achieved by the script above is to approximately reproduce the reduced Stanford Bunny mesh that we use in the other examples.

Annotating meshes

There is one more feature of gl_display.Viewer which is arguably important and which we have not covered. The viewer allows users to select annotation points by ctrl-clicking on the mesh. Simply spawning a viewer as shown in the example below is sufficient to try out this feature. In the example, the annotation points are then printed. If we change the position of an annotation point, it will also move inside the viewer.

pos[:] = pos_backup
viewer = gl.Viewer()
viewer.display(m)
ap = viewer.annotation_points()
for p in ap:
    print(p)
del viewer # or you will have a stale window lying on the desktop

[-0.05125621  0.0946497   0.04424643]
[-0.0061235   0.10561395  0.04395739]
[-0.00760656  0.07759839  0.05759882]

Finding points in space

In geometry processing, we frequently have to deal with collections of irregularly placed points in space. To facilitate queries on this type of data, we have exposed a kD tree class in PyGEL. A kD tree allows us to search for the point closest to a given query point much faster (with asymptotic complexity $O(log(N))$ rather than $O(N)$) than if we simply look through the point list. In the example below, we insert all vertices from our mesh in the I3DTree data structure and then use closest_point to locate the vertex closest to the origin. I have not timed it against the SciPy alternative, so I cannot say if it is faster or slower, but it is less flexible and maybe slightly simpler to use.

from pygel3d import spatial

tree = spatial.I3DTree()

for v in m.vertices():
    tree.insert(pos[v], v)
tree.build()
k,v = tree.closest_point([0,0,0], 1.0)
print("key = ", k, " idx = ", v)

key =  [0.00472706, 0.0340662, 0.00304373]  idx =  1805

Computing distance fields

Another frequently used representation for geometry is distance fields. A distance field is simply a function that maps a point in space to the distance from that point to the closest point on a given surface. The MeshDistance allows precisely for the computation of such signed distances form arbitrary points in space to the closest point on a Manifold.

mdist = hmesh.MeshDistance(m)
print("Distance from origin :", mdist.signed_distance([0,0,0]))

Distance from origin : [0.03452703]

Skeletonization

So far we have not touched upon the use of the Graph class which is one of the newest additions to PyGEL (and GEL). The reason we included this feature was mainly to provide easy access to our skeletonization algorithm Bærentzen and Rotenberg. Since we have been sticking with the bunny so far, let us turn it into skeleton as a final trick. The skeleton of the bunny looks a bit weird in isolation, but you should be able to guess which edges correspond to which features...

from pygel3d import graph

g = graph.from_mesh(bunny)
s = graph.LS_skeleton(g)
jd.display(s)

Computed 901 separators
Found 372 separators
Packed 18 separators
Finding separators: 4.59714
Packing separators: 0.119717