# Chapter 6: Miscellaneous¶

```
import igl
import scipy as sp
import numpy as np
from meshplot import plot, subplot, interact
import os
root_folder = os.getcwd()
```

*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))
```

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)
```