Skip to content

Chapter 6: Miscellaneous

PyPI version buildwheels

import igl
import scipy as sp
import numpy as np
from meshplot import plot, subplot, interact

import os
root_folder = os.getcwd()

Libigl contains a wide variety of geometry processing tools and functions for dealing with meshes and the linear algebra related to them: far too many to discuss in this introductory tutorial. We’ve pulled out a couple of the interesting functions in this chapter to highlight.

Mesh Statistics

Libigl contains various mesh statistics, including face angles, face areas and the detection of singular vertices, which are vertices with more or less than 6 neighbours in triangulations or 4 in quadrangulations.

The example computes these quantities and does a basic statistic analysis that allows to estimate the isometry and regularity of a mesh:

v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "horse_quad.obj"))

## Count the number of irregular vertices, the border is ignored
irregular = igl.is_irregular_vertex(v, f) 
v_count = v.shape[0]
irregular_v_count = np.sum(irregular)
irregular_ratio = irregular_v_count / v_count

print("Irregular vertices: \n%d/%d (%.2f%%)\n"%(irregular_v_count, v_count, irregular_ratio * 100))

## Compute areas, min, max and standard deviation
area = igl.doublearea(v, f) / 2.0

area_avg = np.mean(area)
area_min = np.min(area) / area_avg
area_max = np.max(area) / area_avg
area_ns = (area - area_avg) / area_avg
area_sigma = np.sqrt(np.mean(np.square(area_ns)))

print("Areas (Min/Max)/Avg_Area Sigma: \n%.2f/%.2f (%.2f)\n"%(area_min, area_max, area_sigma))

## Compute per face angles, min, max and standard deviation
angles = igl.internal_angles(v, f)
angles = 360.0 * (angles / (2 * np.pi))

angle_avg = np.mean(angles)
angle_min = np.min(angles)
angle_max = np.max(angles)
angle_ns = angles - angle_avg
angle_sigma = np.sqrt(np.mean(np.square(angle_ns)))

print("Angles in degrees (Min/Max) Sigma: \n%.2f/%.2f (%.2f)\n"%(angle_min, angle_max, angle_sigma))

Irregular vertices: 
659/2400 (27.46%)

Areas (Min/Max)/Avg_Area Sigma: 
0.01/6.26 (0.88)

Angles in degrees (Min/Max) Sigma: 
2.36/171.79 (25.02)

The first row contains the number and percentage of irregular vertices, which is particularly important for quadrilateral meshes when they are used to define subdivision surfaces: every singular point will result in a point of the surface that is only C^1.

The second row reports the area of the minimal element, maximal element and the standard deviation. These numbers are normalized by the mean area, so in the example above 5.33 max area means that the biggest face is 5 times larger than the average face. An ideal isotropic mesh would have both min and max area close to 1.

The third row measures the face angles, which should be close to 60 degrees (90 for quads) in a perfectly regular triangulation. For FEM purposes, the closer the angles are to 60 degrees the more stable will the optimization be. In this case, it is clear that the mesh is of bad quality and it will probably result in artifacts if used for solving PDEs.

Subdivision surfaces

Given a coarse mesh (aka cage) with vertices V and faces F, one can createa higher-resolution mesh with more vertices and faces by subdividing every face. That is, each coarse triangle in the input is replaced by many smaller triangles. Libigl has three different methods for subdividing a triangle mesh.

An “in plane” subdivision method will not change the point set or carrier surface of the mesh. New vertices are added on the planes of existing triangles and vertices surviving from the original mesh are not moved.

By adding new faces, a subdivision algorithm changes the combinatorics of the mesh. The change in combinatorics and the formula for positioning the high-resolution vertices is called the “subdivision rule”.

For example, in the in plane subdivision method of igl.upsample, vertices are added at the midpoint of every edge: v_{ab} = \frac{1}{2}(v_a + v_b) and each triangle (i_a,i_b,i_c) is replaced with four triangles: (i_a,i_{ab},i_{ca}), (i_b,i_{bc},i_{ab}), (i_{ab},i_{bc},i_{ca}), and (i_{bc},i_{c},i_{ca}). This process may be applied recursively, resulting in a finer and finer mesh.

The subdivision method of igl.loop is not in plane. The vertices of the refined mesh are moved to weight combinations of their neighbors: the mesh is smoothed as it is refined (Loop, 1987). This and other smooth subdivision methods can be understood as generalizations of spline curves to surfaces. In particular the Loop subdivision method will converge to a C^1 surface as we consider the limit of recursive applications of subdivision. Away from “irregular” or “extraordinary” vertices (vertices of the original cage with valence not equal to 6), the surface is C^2. The combinatorics (connectivity and number of faces) of igl.loop and igl.upsample are identical: the only difference is that the vertices have been smoothed in igl.loop.

Finally, libigl also implements a form of in plane “false barycentric subdivision” in igl.false_barycentric_subdivision. This method simply adds the barycenter of every triangle as a new vertex v_{abc} and replaces each triangle with three triangles (i_a,i_b,i_{abc}), (i_b,i_c,i_{abc}), and (i_c,i_a,i_{abc}). In contrast to igl.upsample, this method will create triangles with smaller and smaller internal angles and new vertices will sample the carrier surfaces with extreme bias.

ov, of = igl.read_triangle_mesh(os.path.join(root_folder, "data", "decimated-knight.off"))
uv, uf = igl.upsample(ov, of)
lv, lf = igl.loop(ov, of)

p = subplot(ov, of, shading={"wireframe": True}, s=[1, 3, 0])
subplot(uv, uf, shading={"wireframe": True}, s=[1, 3, 1], data=p)
subplot(lv, lf, shading={"wireframe": True}, s=[1, 3, 2], data=p)
p

# @interact(mode=['Coarse','Upsample', 'Loop'])
# def switch(mode):
#     if mode == "Coarse":
#         plot(ov, of, shading={"wireframe": True}, plot=p)
#     if mode == "Upsample":
#         plot(uv, uf, shading={"wireframe": True}, plot=p)
#         #p.update_object(vertices=uv, faces=uf)
#     if mode == "Loop":
#         plot(lv, lf, shading={"wireframe": True}, plot=p)

Data smoothing

A noisy function f defined on a surface \Omega can be smoothed using an energy minimization that balances a smoothing term E_S with a quadratic fitting term:

u = \operatorname{argmin}_u \alpha E_S(u) + (1-\alpha)\int_\Omega ||u-f||^2 dx

The parameter \alpha determines how aggressively the function is smoothed.

A classical choice for the smoothness energy is the Laplacian energy of the function with zero Neumann boundary conditions, which is a form of the biharmonic energy. It is constructed using the cotangent Laplacian L and the mass matrix M: QL = L'*(M\L). Because of the implicit zero Neumann boundary conditions however, the function behavior is significantly warped at the boundary if f does not have zero normal gradient at the boundary.

In (Stein, 2017) it is suggested to use the Biharmonic energy with natural Hessian boundary conditions instead, which corresponds to the hessian energy with the matrix QH = H'*(M2\H), where H is a finite element Hessian and M2 is a stacked mass matrix. The matrices H and QH are implemented in libigl as igl.hessian and igl.hessian_energy respectively.

In the following example the differences between the Laplacian energy with zero Neumann boundary conditions and the Hessian energy can be clearly seen: whereas the zero Neumann boundary condition in the third image bias the isolines of the function to be perpendicular to the boundary, the Hessian energy gives an unbiased result.

The following example shows a function on the beetle mesh, the function with added noise, the result of smoothing with the Laplacian energy and zero Neumann boundary conditions, and the result of smoothing with the Hessian energy.

v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "beetle.off"))
e = igl.edges(f)

# Constructing an exact function to smooth
z_exact = v[0:, 2] + 0.5 * v[0:, 1] + v[0:, 1] * v[0:, 1] + v[0:, 2] * v[0:, 2] * v[0:, 2]

# Make the exact function noisy
s = 0.2 * (np.max(z_exact) - np.min(z_exact))
np.random.seed(5)
z_noisy = z_exact + s * np.random.rand(*z_exact.shape)

# Constructing the squared Laplacian and squared Hessian energy
l = igl.cotmatrix(v, f)
m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_BARYCENTRIC)

m_inv_l = sp.sparse.linalg.spsolve(m, l)
ql = l.T @ m_inv_l
qh = igl.hessian_energy(v, f)

# Solve to find Laplacian-smoothed and Hessian-smoothed solutions
al = 8e-4;
zl = sp.sparse.linalg.spsolve(al * ql + (1 - al) * m, al * m.dot(z_noisy))
ah = 5e-6;
zh = sp.sparse.linalg.spsolve(ah * qh + (1 - ah) * m, ah * m.dot(z_noisy))

# Calculate isolines
ilx_v, ilx_e = igl.isolines(v, f, z_exact, 30)
iln_v, iln_e = igl.isolines(v, f, z_noisy, 30)
ill_v, ill_e = igl.isolines(v, f, zl, 30)
ilh_v, ilh_e = igl.isolines(v, f, zh, 30)


#TODO add edges to subplot
p = subplot(v, f, z_exact, s=[2, 2, 0])
# p.view.add_edges(ilx_v, ilx_e)

subplot(v, f, z_noisy, s=[2, 2, 1], data=p)
# p.view.add_edges(iln_v, iln_e)

subplot(v, f, zl, s=[2, 2, 2], data=p)
# p.view.add_edges(ill_v, ill_e)

subplot(v, f, zh, s=[2, 2, 3], data=p)
# p.view.add_edges(ilh_v, ilh_e)
p

# e_id = p.add_edges(ilx_v, ilx_e)
# @interact(mode=['Original', 'Noisy', 'Biharmonic smoothing (0-Neumann)', 'Biharmonic smoothing (Natural Hessian)'])
# def switch(mode):
#     global e_id
#     p.remove_object(e_id)
#     if mode == "Original":
#         p.update_object(colors=z_exact)
#         e_id = p.add_edges(ilx_v, ilx_e)
#     if mode == "Noisy":
#         p.update_object(colors=z_noisy)
#         e_id = p.add_edges(iln_v, iln_e)
#     if mode == "Biharmonic smoothing (0-Neumann)":
#         p.update_object(colors=zl)
#         e_id = p.add_edges(ill_v, ill_e)
#     if mode == "Biharmonic smoothing (Natural Hessian)":
#         p.update_object(colors=zh)
#         e_id = p.add_edges(ilh_v, ilh_e)

Marching Tetrahedra

Often 3D data is captured as scalar field defined over space f(\mathbf{x}) : \mathcal{R}^3 \rightarrow \mathcal{R}. Lurking within this field, iso-surfaces of the scalar field are often salient geometric objects. The iso-surface at value v is composed of all points \mathbf{x} in \mathcal{R}^3 such that f(\mathbf{x}) = v. A core problem in geometry processing is to extract an iso-surface as a triangle mesh for further mesh-based processing or visualization. This is referred to as iso-contouring.

“Marching Tetrahedra” (Treece, 1999) is a famous method for iso-contouring tri-linear functions f on a 3D simplicial complex (aka a tet mesh). The core idea of this method is to contour the iso-surface passing through each cell (if it does at all) with a predefined topology (aka connectivity) chosen from a look up tabledepending on the function values at each vertex of the cell. The method iterates (“marches”) over all cells (“tetrahedra”) in the complex and stitches together the final mesh.

In libigl, igl.marching_tets constructs a triangle mesh (v,f) approximating the iso-level set for the value isovalue from an input scalar field s sampled at the vertices of a tet mesh locations (tv, tt):

v, f = igl.marching_tets(tv, tt, s, isovalue)

tv = np.load(os.path.join(root_folder, "data", "marching_cube_tv.npy"))
tt = np.load(os.path.join(root_folder, "data", "marching_cube_tt.npy"))
s = np.linalg.norm(tv, axis=1)

svs = []
sfs = []
for i in np.linspace(0.05, 0.75, 15):
    sv, sf, _, _ = igl.marching_tets(tv, tt, s, i)
    svs.append(sv)
    sfs.append(sf)


i = 0
for t in [3, 8, 11]:
    if i == 0:
        p = subplot(svs[t], sfs[t], s = [1, 3, i])
    else:
        subplot(svs[t], sfs[t], s = [1, 3, i], data=p)
    i += 1

p


# @interact(t=(0, 14))
# def update(t=0):
#     global oid
#     p.remove_object(oid)
#     oid = p.add_mesh(svs[t], sfs[t])    

# p = plot(sv, sf, return_plot=True)
# oid = 0

# @interact(t=(0, 14))
# def update(t=0):
#     global oid
#     p.remove_object(oid)
#     oid = p.add_mesh(svs[t], sfs[t])

References

  1. Christian Schüller, Ladislav Kavan, Daniele Panozzo, Olga Sorkine-Hornung. Locally Injective Mappings, 2013. ↩

  2. Qingnan Zhou, Eitan Grinspun, Denis Zorin. Mesh Arrangements for Solid Geometry, 2016 ↩

  3. J Andreas Baerentzen and Henrik Aanaes. Signed distance computation using the angle weighted pseudonormal, 2005. ↩

  4. Sofien Bouaziz, Mario Deuss, Yuliy Schwartzburg, Thibaut Weise, Mark Pauly Shape-Up: Shaping Discrete Geometry with Projections, 2012 ↩

  5. Akash Garg, Alec Jacobson, Eitan Grinspun. Computational Design of Reconfigurables, 2016 ↩

  6. Hugues Hoppe. Progressive Meshes, 1996 ↩

  7. Alec Jacobson, Ladislav Kavan, and Olga Sorkine. Robust Inside-Outside Segmentation using Generalized Winding Numbers, 2013. ↩

  8. Charles Loop. Smooth Subdivision Surfaces Based on Triangles, 1987. ↩

  9. W.E. Lorensen and Harvey E. Cline. Marching cubes: A high resolution 3d surface construction algorithm, 1987. ↩

  10. Michael Rabinovich, Roi Poranne, Daniele Panozzo, Olga Sorkine-Hornung. Scalable Locally Injective Mappings, 2016. ↩

  11. William J. Schroeder, William E. Lorensen, and Steve Linthicum. Implicit Modeling of Swept Surfaces and Volumes, 1994. ↩

  12. Kenshi Takayama, Alec Jacobson, Ladislav Kavan, Olga Sorkine-Hornung. A Simple Method for Correcting Facet Orientations in Polygon Meshes Based on Ray Casting, 2014. ↩

  13. G.M. Treece, R.W. Prager, and A.H.Gee Regularised marching tetrahedra: improved iso-surface extraction, 1999. ↩

  14. Keenan Crane, Clarisse Weischedel, and Max Wardetzky. Geodesics in Heat: A New Approach to Computing Distance Based on Heat Flow, 2013. ↩

  15. Alexander I. Bobenko and Boris A. Springborn. A discrete Laplace-Beltrami operator for simplicial surfaces, 2005. ↩

  16. Zhongshi Jiang, Scott Schaefer, Daniele Panozzo. SCAF: Simplicial Complex Augmentation Framework for Bijective Maps, 2017 ↩