Chapter 4: Parametrization¶
import igl
import scipy as sp
import numpy as np
from meshplot import plot, subplot, interact
import os
root_folder = os.getcwd()
In computer graphics, we denote as surface parametrization a map from the surface to \(\mathbf{R}^2\). It is usually encoded by a new set of 2D coordinates for each vertex of the mesh (and possibly also by a new set of faces in one to one correspondence with the faces of the original surface). Note that this definition is the inverse of the classical differential geometry definition.
A parametrization has many applications, ranging from texture mapping to surface remeshing. Many algorithms have been proposed, and they can be broadly divided in four families:
-
Single patch, fixed boundary: these algorithm can parametrize a disk-like part of the surface given fixed 2D positions for its boundary. These algorithms are efficient and simple, but they usually produce high-distortion maps due to the fixed boundary.
-
Single patch, free boundary: these algorithms let the boundary deform freely, greatly reducing the map distortion. Care should be taken to prevent the border to self-intersect.
-
Global parametrization: these algorithms work on meshes with arbitrary genus. They initially cut the mesh in multiple patches that can be separately parametrized. The generated maps are discontinuous on the cuts (often referred as seams).
-
Global seamless parametrization: these are global parametrization algorithm that hides the seams, making the parametrization “continuous”, under specific assumptions that we will discuss later.
Harmonic parametrization¶
Harmonic parametrization (Eck, 2005) is a single patch, fixed boundary parametrization algorithm that computes the 2D coordinates of the flattened mesh as two harmonic functions.
The algorithm is divided in 3 steps:
- Detection of the boundary vertices
- Map the boundary vertices to a circle
- Compute two harmonic functions (one for u and one for the v coordinate). The harmonic functions use the fixed vertices on the circle as boundary constraints.
The algorithm is coded with libigl in the following example. bnd contains the indices of the boundary vertices, bnd_uv their position on the UV plane, and “1” denotes that we want to compute an harmonic function (2 will be for biharmonic, 3 for triharmonic, etc.). Note that each of the three
functions is designed to be reusable in other parametrization algorithms.
The UV coordinates are then used to apply a procedural checkerboard texture to the mesh.
v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "camelhead.off"))
## Find the open boundary
bnd = igl.boundary_loop(f)
## Map the boundary to a circle, preserving edge proportions
bnd_uv = igl.map_vertices_to_circle(v, bnd)
## Harmonic parametrization for the internal vertices
uv = igl.harmonic(v, f, bnd, bnd_uv, 1)
v_p = np.hstack([uv, np.zeros((uv.shape[0],1))])
p = subplot(v, f, uv=uv, shading={"wireframe": False, "flat": False}, s=[1, 2, 0])
subplot(v_p, f, uv=uv, shading={"wireframe": True, "flat": False}, s=[1, 2, 1], data=p)
# @interact(mode=['3D','2D'])
# def switch(mode):
# if mode == "3D":
# plot(v, f, uv=uv, shading={"wireframe": False, "flat": False}, plot=p)
# if mode == "2D":
# plot(v_p, f, uv=uv, shading={"wireframe": True, "flat": False}, plot=p)
Least squares conformal maps¶
Least squares conformal maps parametrization (Levy, 2002) minimizes the conformal (angular) distortion of the parametrization. Differently from harmonic parametrization, it does not need to have a fixed boundary.
LSCM minimizes the following energy:
\[ E_{LSCM}(\mathbf{u},\mathbf{v}) = \int_X \frac{1}{2}| \nabla \mathbf{u}^{\perp} - \nabla \mathbf{v} |^2 dA \]
which can be rewritten in matrix form as (Mullen, 2008):
\[ E_{LSCM}(\mathbf{u},\mathbf{v}) = \frac{1}{2} [\mathbf{u},\mathbf{v}]^t (L_c - 2A) [\mathbf{u},\mathbf{v}] \]
where L_c is the cotangent Laplacian matrix and A is a matrix such that [\mathbf{u},\mathbf{v}]^t A [\mathbf{u},\mathbf{v}] is equal to the vector area of the mesh.
Using libigl, this matrix energy can be written in a few lines of code. The
cotangent matrix can be computed using igl.cotmatrix. Note that we want to apply the Laplacian matrix to the u and v coordinates at the same time, thus we need to extend it taking the left Kronecker product with a 2x2 identity matrix. The area matrix is computed with igl.vector_area_matrix.
The final energy matrix is L_{flat} - 2A. Note that in this case we do not need to fix the boundary. To remove the null space of the energy and make the minimum unique, it is sufficient to fix two arbitrary vertices to two arbitrary positions. The full source code is provided in the following LSCM parametrization example.
v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "camelhead.off"))
# Fix two points on the boundary
b = np.array([2, 1])
bnd = igl.boundary_loop(f)
b[0] = bnd[0]
b[1] = bnd[int(bnd.size / 2)]
bc = np.array([[0.0, 0.0], [1.0, 0.0]])
# LSCM parametrization
_, uv = igl.lscm(v, f, b, bc)
p = subplot(v, f, uv=uv, shading={"wireframe": False, "flat": False}, s=[1, 2, 0])
subplot(uv, f, uv=uv, shading={"wireframe": False, "flat": False}, s=[1, 2, 1], data=p)
# @interact(mode=['3D','2D'])
# def switch(mode):
# if mode == "3D":
# plot(v, f, uv=uv, shading={"wireframe": False, "flat": False}, plot=p)
# if mode == "2D":
# plot(uv, f, uv=uv, shading={"wireframe": True, "flat": False}, plot=p)
As-rigid-as-possible parametrization¶
As-rigid-as-possible parametrization (Liu, 2008) is a powerful single-patch, non-linear algorithm to compute a parametrization that strives to preserve distances (and thus angles). The idea is very similar to ARAP surface deformation: each triangle is mapped to the plane trying to preserve its original shape, up to a rigid rotation.
The algorithm can be implemented reusing the functions discussed in the
deformation chapter: igl.ARAP and arap.solve. The only
difference is that the optimization has to be done in 2D instead of 3D and that
we need to compute a starting point. While for 3D deformation the optimization
is bootstrapped with the original mesh, this is not the case for ARAP
parametrization since the starting point must be a 2D mesh.
In the following example, we initialize the optimization with harmonic parametrization. Similarly to LSCM, the boundary is free to deform to minimize the distortion.
v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "camelhead.off"))
## Find the open boundary
bnd = igl.boundary_loop(f)
## Map the boundary to a circle, preserving edge proportions
bnd_uv = igl.map_vertices_to_circle(v, bnd)
## Harmonic parametrization for the internal vertices
uv = igl.harmonic(v, f, bnd, bnd_uv, 1)
arap = igl.ARAP(v, f, 2, np.zeros(0))
uva = arap.solve(np.zeros((0, 0)), uv)
p = subplot(v, f, uv=uva, shading={"wireframe": False, "flat": False}, s=[1, 2, 0])
p = subplot(uva, f, uv=uva, shading={"wireframe": False, "flat": False}, s=[1, 2, 1], data=p)
# @interact(mode=['3D','2D'])
# def switch(mode):
# if mode == "3D":
# plot(v, f, uv=uva, shading={"wireframe": False, "flat": False}, plot=p)
# if mode == "2D":
# plot(uva, f, uv=uva, shading={"wireframe": True, "flat": False}, plot=p)
Planarization¶
A quad mesh can be transformed in a planar quad mesh with Shape-Up (Bouaziz, 2012), a local/global approach that uses the global step to enforce surface continuity and the local step to enforce planarity.
The following example planarizes a quad mesh until it satisfies a user-given planarity threshold. The colors represent the planarity of the quads.
# Load a quad mesh generated by a conjugate field
vqc, fqc, _ = igl.read_off(os.path.join(root_folder, "data", "inspired_mesh_quads_Conjugate.off"))
# Convert it to a triangle mesh
fqc_tri = np.zeros((fqc.shape[0] * 2, 3), dtype="int64")
fqc_tri[:fqc.shape[0]] = fqc[:, :3]
fqc_tri[fqc.shape[0]:, 0] = fqc[:, 2]
fqc_tri[fqc.shape[0]:, 1] = fqc[:, 3]
fqc_tri[fqc.shape[0]:, 2] = fqc[:, 0]
# Planarize it
vqc_p = igl.planarize_quad_mesh(vqc, fqc, 10, 0.005)
# Calculate a color to each quad that corresponds to its planarity
planarity = igl.quad_planarity(vqc, fqc)
planarity_p = igl.quad_planarity(vqc_p, fqc)
c = np.concatenate([planarity, planarity])
c_p = np.concatenate([planarity_p, planarity_p])
c_range = [min(np.min(c), np.min(c_p)), max(np.max(c), np.max(c_p))]
p = subplot(vqc, fqc_tri, c, shading={"normalize": c_range}, s=[1, 2, 0])
subplot(vqc_p, fqc_tri, c_p, shading={"normalize": c_range}, s=[1, 2, 1], data=p)
# @interact(mode=['Curved','Planar'])
# def switch(mode):
# if mode == "Curved":
# p.update_object(colors=c)
# if mode == "Planar":
# p.update_object(vertices=vqc_p, colors=c_p)