|
CGAL 4.5 - Triangulated Surface Mesh Deformation
|
|
CGAL 4.5 - Triangulated Surface Mesh Deformation
|
This package offers surface mesh deformation algorithms which compute new vertex positions of a surface mesh under positional constraints of some of its vertices, without requiring any additional structure other than the surface mesh itself
This package implements the algorithm described in [5] together with an alternative energy function [3]. The algorithm minimizes a nonlinear deformation energy under positional constraints to preserve rigidity as much as possible. The minimization of the energy relies on solving sparse linear systems and finding closest rotation matrices.
A surface mesh deformation system consists of:
A vertex from the ROI that is not a control vertex is called an unconstrained vertex. These definitions are depicted in Figure 53.1.
In this package, two algorithms are implemented:
Given an edge weighting scheme, both methods iteratively minimize an energy function and produce a different surface mesh at each step until convergence is reached.
Spokes and Rims is the default method proposed by the package. It provides an unconditional convergence while the ARAP method requires the edge weights to be positive. However, the results obtained using the Spokes and Rims method are more dependent on the discretization of the deformed surface (See Figure 53.2).
More details on these algorithms are provided in section Deformation Techniques, Energies and Weighting Schemes.
The deformation methods implemented rely on solving a sparse linear system. The sparse matrix definition depends on the weighting scheme and on the unconstrained and control vertices. The right term depends only on the target positions of the control vertices.
The deformation process is handled by the class Surface_mesh_deformation and the surface mesh is represented as a halfedge data structure that must be a model of the HalfedgeGraph concept.
The class Surface_mesh_deformation provides two groups of functions for the preprocessing (sparse matrix definition) and the deformation (right hand term definition).
The preprocessing consists in computing a factorization of the aforementioned sparse matrix to speed up the linear system resolution. It requires the ROI to be defined. The following conventions are used for the definition of the ROI:
Each time the ROI is modified, the preprocessing function preprocess() must be called. Note that if it is not done, the first deformation step calls this function automatically and has a longer runtime compared to subsequent deformation steps.
The function Surface_mesh_deformation::preprocess() returns true if the factorization is successful, and false otherwise.
Rank deficiency is the main reason for failure. Typical failure cases are:
SurfaceModelingWeights) features too many zeros and breaks the connectivity informationThe choice of the weighting scheme provides a mean to adjust the way the control vertices influences the unconstrained vertices. The defaults provides satisfactory results in general but other weighting schemes may be selected or designed to experiment or improve the results in specific cases.
The deformation of the surface mesh is triggered by the displacement of the control vertices. This is achieved through setting the target positions of the control vertices (directly or by using an affine transformation to be applied to a control vertex or a range of control vertices).
Note that a rotation or a translation of a control vertex is always applied on its last target position set: they are cumulative.
The deformation of the surface mesh happens when calling the function Surface_mesh_deformation::deform(). The number of optimization iterations varies depending on whether the user chooses a fixed number of iterations or a stopping criterion based on the energy variation.
After the call to the deformation function, the input surface mesh is updated and the control vertices are at their target positions and the unconstrained vertices are moved accordingly. The function Surface_mesh_deformation::deform() can be called several times consecutively, in particular if the convergence has not been reached yet (otherwise it has no effect).
Surface_mesh_deformation::preprocess() is first called. The original positions can be updated by calling Surface_mesh_deformation::overwrite_initial_geometry() ( which will also require a new preprocessing step). This behavior is illustrated in Video 1.Two deformation techniques are provided by this package. This section summarizes from the user point of view what is explained in details in the section Deformation Techniques, Energies and Weighting Schemes.
The As-Rigid-As-Possible deformation technique requires the use of a positive weighting scheme to guarantee the correct minimization of the energy. When using the default cotangent weighting scheme, this means that the input surface mesh must be clean. That is, that for all edges in the surface mesh the sum of the angles opposite to the edge in the incident triangles is less that \( \pi \). If this is not the case and the targeted application allows the modification of the surface mesh connectivity, a solution (amongst other) is to bissect (possibly recursively) the problematic edges. See Figure 53.3.
If the mesh connectivity must be preserved, the Spokes and Rims deformation technique is guaranteed to always correctly minimize the energy even if the weights are negative. However, this technique is more dependent on the discretization of the deformed surface (See Figure 53.2).
In this example, the whole surface mesh is used as ROI and a few vertices are added as control vertices. Surface_mesh_deformation::set_target_position() is used for setting the target positions of the control vertices.
File Surface_modeling/all_roi_assign_example.cpp
deform_1.off and deform_2.off. In this example, we use the functions translate() and rotate() on a range of control vertices. Note that the translations and the rotations are defined using a 3D vector type and a quaternion type from the Eigen library.
File Surface_modeling/k_ring_roi_translate_rotate_example.cpp
deform_1.off and deform_2.off. In the previous examples, we used an enriched polyhedron storing an ID in its halfedges and vertices together with the default property maps in the deformation object to access them. In the following example, we show how we can use alternative property maps.
For practical performance however we recommend relying upon the former examples instead, as using a std::map to access indices increases the complexity from constant to logarithmic.
File Surface_modeling/deform_polyhedron_with_custom_pmap_example.cpp
Using a custom weighting scheme for edges is also possible if one provides a model of SurfaceModelingWeights. In this example, the weight of each edge is pre-computed and an internal map is used for storing and accessing them.
Another example is given in the manual page of the concept SurfaceModelingWeights.
File Surface_modeling/custom_weight_for_edges_example.cpp
A plugin for the polyhedron demo is available to test the algorithm. The following video tutorials explain how to use it. When the deformation dock window is open, the picking of control vertices and of the ROI is done by pressing Shift and clicking with the left button of the mouse. The displacement of the vertices is triggered when the Ctrl button is pressed.
This section gives the theoretical background to make the user manual self-contained and at the same time explains where the weights comes in. This allows advanced users of this package to tune the weighting scheme by developing a model of the concept SurfaceModelingWeights used in the class Surface_mesh_deformation.
The Laplacian representation (referred to as Laplace coordinates in [2]) of a vertex in a surface mesh is one way to encode the local neighborhood of a vertex in the surface mesh. In this representation, a vertex \( \mathbf{v}_i \) is associated a 3D vector defined as:
\[ \begin{equation} L(\mathbf{v}_i) = \sum_{\mathbf{v}_j \in N(\mathbf{v}_i)} w_{ij}(\mathbf{v}_i - \mathbf{v}_j), \label{eq:lap_open} \end{equation} \]
where:
The simplest choice for the weights is the uniform scheme where \( w_{ij}=1/|N(\mathbf{v}_i)| \) for each adjacent vertex \(\mathbf{v}_j\). In this case, the Laplacian representation of a vertex is the vector between this vertex and the centroid of its adjacent vertices (Figure 53.6).
In the surface mesh deformation context, a popular choice is the cotangent weight scheme that derives from the discretization of the Laplace operator [4] : Given an edge of the surface mesh, its corresponding cotangent weight is the mean of the cotangents of the angles opposite to the edge. It was shown to produce results that are not biased by the surface mesh of the approximated surface [5].
Considering a surface mesh with \(n\) vertices, it is possible to define its Laplacian representation \(\Delta\) as a \(n \times 3\) matrix:
\[ \begin{equation} \mathbf{L}\mathbf{V} = \Delta, \label{eq:lap_system} \end{equation} \]
where:
This section is an introduction to provide the background for the next two sub-sections describing the algorithms implemented in this package. A system relying only on the approach described below results in non-smooth transitions in the neighborhood of the control vertices. For a survey on different Laplacian-based editing techniques we refer to [1].
The main idea behind Laplacian-based deformation techniques is to preserve the Laplacian representation under deformation constraints. The Laplacian representation of a surface mesh is treated as a representative form of the discretized surface, and the deformation process must follow the deformation constraints while preserving the Laplacian representation as much as possible.
There are different ways to incorporate deformation constraints into the deformation system [1]. This package supports hard constraints, that is, target positions of control vertices are preserved after the deformation.
Given a surface mesh deformation system with a ROI made of \( n \) vertices and \( k \) control vertices, we consider the following linear system:
\[ \begin{equation} \left[ \begin{array}{ccc} \mathbf{L}_f\\ 0 \; \mathbf{I}_c \end{array} \right] \mathbf{V} = \left[ \begin{array}{ccc} {\Delta}_f \\ \mathbf{V}_c \end{array} \right], \label{eq:lap_energy_system} \end{equation} \]
where:
The left-hand side matrix of the system of Eq. \(\eqref{eq:lap_energy_system}\) is a square non-symmetric sparse matrix. To solve the aforementioned system, an appropriate solver (e.g. LU solver) needs to be used. Note that solving this system preserves the Laplacian representation of the surface mesh restricted to the unconstrained vertices while satisfying the deformation constraints.
Given a surface mesh \(M\) with \( n \) vertices \( \{\mathbf{v}_i\} i \in \{1 \dots n \} \) and some deformation constraints, we consider the following energy function:
\[ \begin{equation} \sum_{\mathbf{v}_i \in M} \sum_{\mathbf{v}_j \in N(\mathbf{v}_i)} w_{ij} \left\| (\mathbf{v}'_i - \mathbf{v}'_j) - \mathbf{R}_i(\mathbf{v}_i - \mathbf{v}_j) \right\|^2, \label{eq:arap_energy} \end{equation} \]
where:
An as-rigid-as possible surface mesh deformation [5] is defined by minimizing this energy function under the deformation constraints, i.e. the assigned position \( {v}'_i\) for each vertex \( \mathbf{v}_i\) in the set of control vertices. Defining the one-ring neighborhood of a vertex as its set of adjacent vertices, the intuitive idea behind this energy function is to allow each one-ring neighborhood of vertices to have an individual rotation, and at the same time to prevent shearing by taking advantage of the overlapping of one-ring neighborhoods of adjacent vertices (see Figure 53.7).
There are two unknowns per vertex in Eq. \(\eqref{eq:arap_energy}\): the new positions ( \(\mathbf{v}'_k\)) of the unconstrained vertices and the rotation matrices ( \(\mathbf{R}_i\)). If the energy contribution of each vertex is positive, this boils down to minimizing the energy contribution of each vertex \(\mathbf{v}_i\).
Each such term of the energy is minimized by using a two-step optimization approach (also called local-global approach).
In the first step, the positions of the vertices are considered as fixed so that the rotation matrices are the only unknowns.
For the vertex \(\mathbf{v}_i\), we consider the covariance matrix \(\mathbf{S}_i\):
\[ \begin{equation} \mathbf{S}_i = \sum_{\mathbf{v}_j \in N(\mathbf{v}_i)} w_{ij} (\mathbf{v}_i - \mathbf{v}_j)(\mathbf{v}'_i - \mathbf{v}'_j)^T, \label{eq:cov_matrix} \end{equation} \]
It was shown [6] that minimizing the energy contribution of \(\mathbf{v}_i\) in Eq. \(\eqref{eq:arap_energy}\) is equivalent to maximizing the trace of the matrix \(\mathbf{R}_i \mathbf{S}_i\). \(\mathbf{R}_i \) is the transpose of the unitary matrix in the polar decomposition of \(\mathbf{S}_i\).
In the second step, the rotation matrices are substituted into the partial derivative of Eq. \(\eqref{eq:arap_energy}\) with respect to \(\mathbf{v}'_i\). Assuming the weights are symmetric, setting the derivative to zero results in the following equation:
\[ \begin{equation} \sum_{\mathbf{v}_j \in N(\mathbf{v}_i)} w_{ij}(\mathbf{v}'_i - \mathbf{v}'_j) = \sum_{\mathbf{v}_j \in N(\mathbf{v}_i)} w_{ij} \frac{(\mathbf{R}_i + \mathbf{R}_j)}{2} (\mathbf{v}_i - \mathbf{v}_j). \label{eq:lap_ber} \end{equation} \]
The left-hand side of this equation corresponds to the one of Eq. \(\eqref{eq:lap_open}\), and we can set \(\Delta\) to be the right-hand side. Solving the linear system in Eq. \(\eqref{eq:lap_energy_system}\) gives the new positions of the unconstrained vertices.
This two-step optimization can be applied several times iteratively to obtain a better result.
The original algorithm [5] we described assumes that:
\[ \begin{equation} \sum_{\mathbf{v}_j \in N(\mathbf{v}_i)} (w_{ij} + w_{ji})(\mathbf{v}'_i - \mathbf{v}'_j) = \sum_{\mathbf{v}_j \in N(\mathbf{v}_i)} (w_{ij}\mathbf{R}_i + w_{ji}\mathbf{R}_j)(\mathbf{v}_i - \mathbf{v}_j). \label{eq:lap_ber_asym} \end{equation} \]
A method minimizing another energy function is described next to avoid the latter issue.
The elastic energy function proposed by [3] additionally takes into account all the opposite edges in the facets incident to a vertex. The energy function to minimize becomes:
\[ \begin{equation} \sum_{\mathbf{v}_i \in M} \sum_{(\mathbf{v}_j, \mathbf{v}_k) \in E(\mathbf{v}_i)} w_{jk} \left\| (\mathbf{v}'_j - \mathbf{v}'_k) - \mathbf{R}_i(\mathbf{v}_j - \mathbf{v}_k) \right\|^2, \label{eq:arap_energy_rims} \end{equation} \]
where \(E(\mathbf{v}_i)\) consists of the set of edges incident to \(\mathbf{v}_i\) (the spokes) and the set of edges in the link (the rims) of \(\mathbf{v}_i\) in the surface mesh \(M\) (see Figure 53.8).
The method to get the new positions of the unconstrained vertices is similar to the two-step optimization method explained in As-Rigid-As Possible Deformation. For the first step, the Eq. \(\eqref{eq:cov_matrix}\) is modified to take into account the edges in \(E(\mathbf{v}_i)\):
\[ \begin{equation} \mathbf{S}_i = \sum_{(\mathbf{v}_j, \mathbf{v}_k) \in E(\mathbf{v}_i)} w_{jk} (\mathbf{v}_j - \mathbf{v}_k)(\mathbf{v}'_j - \mathbf{v}'_k)^T, \label{eq:cov_matrix_sr} \end{equation} \]
For the second step, setting partial derivative of Eq. \(\eqref{eq:arap_energy_rims}\) to zero with respect to \(\mathbf{v}_i\) gives the following equation:
\[ \begin{equation} \sum_{\mathbf{v}_j \in N(\mathbf{v}_i)} (w_{ij} + w_{ji})(\mathbf{v}'_i - \mathbf{v}'_j) = \sum_{\mathbf{v}_j \in N(\mathbf{v}_i)} \frac{w_{ij}(\mathbf{R}_i + \mathbf{R}_j + \mathbf{R}_m) + w_{ji}(\mathbf{R}_i + \mathbf{R}_j + \mathbf{R}_n)}{3} (\mathbf{v}_i - \mathbf{v}_j). \label{eq:lap_ber_rims} \end{equation} \]
where \(\mathbf{R}_m\) and \(\mathbf{R}_n\) are the rotation matrices of the vertices \(\mathbf{v}_m\), \(\mathbf{v}_n\) which are the opposite vertices of the edge \(\mathbf{v}_i \mathbf{v}_j\) (see Figure 53.8). Note that if the edge \( \mathbf{v}_i \mathbf{v}_j \) is on the boundary of the surface mesh, then \( w_{ij} \) must be 0 and \( \mathbf{v}_m \) does not exist.
An important property of this approach compared to As-Rigid-As Possible Deformation is that the contribution to the global energy of each vertex is guaranteed to be non-negative when using the cotangent weights [3]. Thus even with negative weights, the minimization of the energy with the iterative method presented is always guaranteed. However, this method is more dependent on the discretization of the deformed surface (See Figure 53.2).
The implementation in this package uses the cotangent weights by default (negative values included) as proposed in [3].
An initial version of this package has been implemented during the 2011 Google Summer of Code by Yin Xu under the guidance of Olga Sorkine and Andreas Fabri. Ilker O. Yaz took over the finalization of the package with the help of Sébastien Loriot for the documentation and the API. The authors are grateful to Gaël Guennebaud for his great help on using the Eigen library and for providing the code to compute the closest rotation.
1.8.4