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

p = subplot(v, f, z_exact, s=[2, 2, 0])

subplot(v, f, z_noisy, s=[2, 2, 1], data=p)

subplot(v, f, zl, s=[2, 2, 2], data=p)

subplot(v, f, zh, s=[2, 2, 3], data=p)
p

# @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)