cgal/Surface_modeling/doc/Surface_modeling/Surface_modeling.txt

442 lines
23 KiB
Plaintext

namespace CGAL {
/*!
\mainpage User Manual
\anchor Chapter_SurfaceModeling
\cgalAutoToc
\authors Olga Sorkine-Hornung, Yin Xu and Ilker %O. Yaz
\image latex main_image_suggestion_4_resized.png
\image html main_image_suggestion_4_resized.png
\section Surface_Modeling_Intro Introduction
A surface mesh deformation algorithm provides new positions of the vertices 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 \cgalCite{Sorkine2007AsRigidAs} together with an
alternative energy function \cgalCite{Chao2010SimpleGeomModel}.
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. More details
are provided in \ref Surface_Modeling_Overview.
\section Surface_Modeling_Def Definitions
A <em>surface mesh deformation system</em> consists of:
- a triangulated surface mesh (<em>surface mesh</em> in the following),
- a set of vertices defining the region to deform (called the <em>region-of-interest</em> and abbreviated <em>ROI</em>),
- a subset of vertices from the ROI that the user wants to move (called the <em>control vertices</em>),
- a target position for each control vertex (defining the <em>deformation constraints</em>).
A vertex from the ROI that is not a control vertex is called a <em>free vertex</em>. These definitions are illustrated on \cgalFigureRef{Simple_roi}.
\cgalFigureBegin{Simple_roi, roi_example.png}
The ROI features the green vertices (the free vertices) and the red vertices (the control vertices).
(Left) The initial surface mesh; (Right) A target position is defined for each control vertex,
the coordinates of the free vertices are updated by the deformation algorithm.
\cgalFigureEnd
In this package, two slightly different algorithms are implemented:
- The <em>As-Rigid-As-Possible</em> (<em>ARAP</em>) method described in \cgalCite{Sorkine2007AsRigidAs}
- The <em>Spokes and Rims</em> method \cgalCite{Chao2010SimpleGeomModel} that is a slightly modified version of ARAP that is more
stable but produces asymmetric results (See \cgalFigureRef{Arap_spokes_comparison}).
\cgalFigureBegin{Arap_spokes_comparison, arap_spokes_comparison.png}
Comparison between \ref Surface_Modeling_Overview_ARAP and \ref Surface_Modeling_Overview_ARAP_Rims.
On a square with spikes mesh, the ROI consists of the green vertices. The control vertices are the red ones.
We translate the control vertices along the normal of the plane and observe the result produced by \ref
Surface_Modeling_Overview_ARAP (center) and \ref Surface_Modeling_Overview_ARAP_Rims (right) from the same view point. The
latter version does not produce a symmetric result.
\cgalFigureEnd
Both methods are iterative and produce a different mesh at each step until convergence is reached.
For more details on these algorithms, see section \ref Surface_Modeling_Overview.
\section Surface_Modeling_API User Interface Description
The deformation process is handled by the class `Deform_mesh` and the
surface mesh is represented as a halfedge data structure that must
be a model of the `HalfedgeGraph` concept.
The class `Deform_mesh` provides two groups of functions for the preprocessing and the deformation.
<!--
We present the API by considering the user point of view during a deformation session
(see \cgalFigureRef{Flow_of_API} for an illustration).
\cgalFigureBegin{Flow_of_API, usage_flow_of_API.png}
The deformation work flow.
\cgalFigureEnd
-->
\subsection Preprocess_section Preprocessing
The preprocessing requires the ROI to be set (both control and free vertices).
The following conventions are used for the definition of the ROI:
- a vertex added as a control vertex is added to the ROI
- a control vertex removed from the ROI is no longer considered as a control vertex
- a control vertex that is removed is not removed from the ROI
Once the free and control vertices are set, the function `Deform_mesh::preprocess()` should be called to speed up solving the system.
Note that if it is not done, the first deformation step will call this function automatically and will have a longer runtime
compared to subsequent deformation steps.
This function also extends the ROI by adding as <em>fixed control vertex</em> any vertex not in the ROI that is edge-adjacent
to a vertex of the ROI. The target position of a fixed control vertex is its current position and cannot be changed.
The function `Deform_mesh::preprocess()` returns `true` if the factorization of the linear system in Eq.
\f$\eqref{eq:lap_energy_system}\f$ is successful, and `false` otherwise.
Rank deficiency is the main reason for failure. Typical failure cases are:
- All the vertices are in the ROI and no control vertices are set
- The weighting scheme used features too many zero and breaks the connectivity information
\note The ROI does not need to be a connected component of the surface mesh graph. However for better performances
it is better to use an individual deformation object for each connected component.
\subsection Deform_section Deformation
The deformation of the surface mesh is triggered by the displacement of the control vertices.
This is done by setting the target position 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 mesh happens when calling the function `Deform_mesh::deform()`. The number of optimization iterations
varies depending on whether the user chose a fixed number of iterations or a stopping criteria based on the energy variation.
After the call to the deformation function, the input mesh is updated and the control vertices are at
their target positions and the free vertices are moved accordingly.
The function `Deform_mesh::deform()` can be called as many time as desired.
<!-- The current deformed mesh is only used as the initial guess. -->
\warning Vertices can be inserted into or removed from the ROI and the set of control vertices at any time.
In particular, any vertex which is no longer inside the ROI will be assigned to its original position when
`Deform_mesh::preprocess()` is first called.
The original positions can be updated by calling `Deform_mesh::overwrite_original_positions()` (
which will also require a new preprocessing step).
This behavior is illustrated in the video \ref Surface_Modeling_Demo_Tutorial.
\subsection Surface_modeling_examples Examples
\subsubsection Example_1 Using the Whole Mesh as Region-of-Interest
In this example, the whole mesh is used as ROI and a few vertices are added as control vertices.
`Deform_mesh::set_target_position()` is used for setting the target positions of the control vertices.
\cgalExample{Surface_modeling/all_roi_assign_example.cpp}
\cgalFigureBegin{Example_1_results, example_1_results.png}
Deformation results when running example \ref Example_1 : `deform_1.off` and `deform_2.off`.
\cgalFigureEnd
\subsubsection Example_2 Using Affine Transformation on Range of Vertices
In this example, we use the functions `translate()` and `rotate()` on a range of control vertices.
\cgalExample{Surface_modeling/k_ring_roi_translate_rotate_example.cpp}
\cgalFigureBegin{Example_2_results, example_2_results.png}
Deformation results when running example \ref Example_2 : `deform_1.off` and `deform_2.off`.
\cgalFigureEnd
\subsubsection Example_3 Using an Enriched Polyhedron with Ids
In this example, we use an <i>enriched</i> polyhedron storing an ID in its halfedges and vertices
and use it to store indices. The main advantage is to decrease the complexity of accessing
associated indices to edges and vertices from logarithmic to constant.
\cgalExample{Surface_modeling/enriched_polyhedron_with_custom_pmap_example.cpp}
\subsubsection Example_4 Using a Custom Edge Weighting Scheme
Using a custom weighting scheme for edges is also possible if one provides a model of `SurfaceModelingWeightCalculator`.
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 on the manual page of the concept `::SurfaceModelingWeightCalculator`.
\cgalExample{Surface_modeling/custom_weight_for_edges_example.cpp}
\section Surface_Modeling_Demo How to Use the Demo
\todo update these videos when the reivew is converging
A plugin for the polyhedron demo is available to test the algorithm. The following video tutorials explain how to use it.
\subsection Surface_Modeling_Demo_Showcase An Introduction and Real World Example
@htmlonly
<center>
<iframe width="420" height="315" src="http://www.youtube.com/embed/WvC47fVSPGU" frameborder="0" allowfullscreen></iframe>
<iframe width="420" height="315" src="http://www.youtube.com/embed/0oTi5aKdcmg" frameborder="0" allowfullscreen></iframe>
</center>
@endhtmlonly
\subsection Surface_Modeling_Demo_Tutorial Demonstration for Deleting Control Vertices / ROI and Overwrite
@htmlonly
<center>
<iframe width="420" height="315" src="http://www.youtube.com/embed/xPAG_AOLvU0" frameborder="0" allowfullscreen></iframe>
</center>
@endhtmlonly
\section Surface_Modeling_Overview Overview of the Method
The algorithm implemented is a so-called <em>Laplacian-based deformation</em> technique. We first give an overview of
the principal of these methods.
\subsection Surface_Modeling_Overview_Laplacian Laplacian Representation
The <em>Laplacian coordinates</em> (also called <em>Laplace coordinates</em> \cgalCite{botsch2010polygon})
of a vertex in a surface mesh are one way to <em>encode</em> the local neighborhood of a vertex in the mesh.
The Laplacian coordinates of a vertex \f$ \mathbf{v}_i \f$ is a 3D vector defined as:
\f[
\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}
\f]
where:
- \f$N(\mathbf{v}_i)\f$ denotes the set of vertices edge-adjacent to \f$\mathbf{v}_i\f$;
-\f$w_{ij}\f$ denotes a weight between \f$\mathbf{v}_i\f$ and \f$\mathbf{v}_j\f$.
The simplest choice for the weights is the uniform scheme where \f$ w_{ij}=1/|N(\mathbf{v}_i)| \f$ for each adjacent vertex
\f$\mathbf{v}_j\f$. In this case, the Laplacian coordinates of a vertex is the vector between this vertex
and the centroid of its adjacent vertices (\cgalFigureRef{Simple_laplacian}).
The cotangent weight scheme \cgalCite{Pinkall1993Cotangent} is more independent from the discretization of the original
surface: Given an edge of the surface mesh, its corresponding cotangent weight is the mean of the
cotangents of the angles opposite to the edge.
\cgalFigureBegin{Simple_laplacian,simple_mesh_with_laplacian.png}
Laplacian coordinates of \f$ v_i \f$ with uniform weights: the red square vertex is the centroid of the vertices adjacent to \f$ v_i \f$. The Laplacian coordinates \f$ \delta \f$ are represented as a blue vector.
\cgalFigureEnd
Considering the whole surface mesh, it is possible to write a linear system describing the Laplacian coordinates of all its vertices:
\f[
\begin{equation}
\mathbf{L}\mathbf{V} = \Delta,
\label{eq:lap_system}
\end{equation}
\f]
where:
- \f$n\f$ denotes the number of vertices in the surface mesh,
- \f$\mathbf{L}\f$ is a \f$n \times n\f$ sparse matrix, called the <em>Laplacian matrix</em>. Its elements \f$ m_{ij} \f$, \f$i,j \in \{1 \dots n\} \f$ are defined as follows:
- \f$ m_{ii} \f$ = \f$ \sum_{\mathbf{v}_j \in N(\mathbf{v}_i)}w_{ij} \f$,
- \f$ m_{ij} = -w_{ij} \f$ if \f$ \mathbf{v}_j \in N(\mathbf{v_i}) \f$,
- \f$ 0 \f$ otherwise.
- \f$\mathbf{V}\f$ is a \f$n \times 3\f$ matrix made of the %Cartesian coordinates of the vertices;
- \f$\Delta\f$ is a \f$n \times 3\f$ matrix called the Laplacian coordinates of the surface mesh.
\subsection Surface_Modeling_Overview_Laplacian_Deformation Naive Laplacian Deformation
The main idea behind Laplacian-based deformation techniques is to preserve the Laplacian coordinates
under deformation constraints. Laplacian coordinates of the surface mesh are treated as a representative form of the discretized surface,
and the deformation process must follow the deformation constraints while preserving the Laplacian coordinates as much as possible.
There are different ways to incorporate deformation constraints into the deformation system \cgalCite{Botsch2008OnLinearVariational}.
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 \f$ n \f$ vertices and \f$ k \f$ control vertices, we consider the following linear system:
\f[
\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}
\f]
where:
- \f$\mathbf{V}\f$ is a \f$n \times 3\f$ matrix. It is the unknown of the system that represents the coordinates of the vertices after deformation. The system is build such that the \f$ k \f$ last row corresponds to the control vertices.
- \f$\mathbf{L}_f\f$ denotes the Laplacian matrix of the free vertices. It is a \f$n-k \times n \f$ matrix as defined in Eq. \f$\eqref{eq:lap_system}\f$, but removing the rows corresponding to the control vertices.
- \f$\mathbf{I}_c\f$ is a \f$k \times k\f$ identity matrix.
- \f${\Delta}_f\f$ denotes the Laplacian coordinates of the free vertices as defined in Eq. \f$\eqref{eq:lap_system}\f$, but removing the rows corresponding to the control vertices.
- \f$\mathbf{V}_c\f$ is a \f$k \times 3\f$ matrix containing the %Cartesian coordinates of the target positions of the control vertices.
The left-hand side matrix of the system of Eq.\f$\eqref{eq:lap_energy_system}\f$ 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 coordinates of the surface mesh restricted to the free vertices while satisfying the deformation constraints.
\note This section is an introduction to provide the background for the two next sections describing the algorithms implemented
in this package. A system relying only on the approach described above results in non-smooth transitions in the neighborhood of
the control vertices. For a survey on different Laplacian-based editing techniques please refer to
\cgalCite{Botsch2008OnLinearVariational}.
\subsection Surface_Modeling_Overview_ARAP As-Rigid-As Possible Deformation
Given a surface mesh \f$M\f$ with \f$ n \f$ vertices \f$ \{\mathbf{v}_i\} i \in \{1 \dots n \} \f$ and some deformation
constraints, we consider the following energy function:
\f[
\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}
\f]
where:
- \f$\mathbf{R}_i\f$ is a \f$ 3 \times 3 \f$ rotation matrix for the vertex \f$\mathbf{v}_i\f$
- \f$w_{ij}\f$ denotes a weight
- \f$N(\mathbf{v}_i)\f$ denotes the set of vertices adjacent to \f$\mathbf{v}_i\f$ in \f$M\f$
- \f$N(\mathbf{v}'_i)\f$ deontes a new position of the vertex \f$N(\mathbf{v}_i)\f$ after a given deformation
An as-rigid-as possible surface mesh deformation \cgalCite{Sorkine2007AsRigidAs} is defined by minimizing
this energy function under the deformation constraints, that is fixing \f$ {v}'_i\f$ for each vertex \f$ \mathbf{v}_i\f$
in the set of control vertices. The intuitive idea behind this energy function is to allow each one-ring neighborhood of
vertices to have an individual rotation, but at the same time prevents shearing by taking advantage of the
overlapping of one-ring neighborhood of adjacent vertices (see \cgalFigureRef{Overlapping_cells}).
\cgalFigureBegin{Overlapping_cells,overlapping_cells.png}
Overlaps of one-ring neighborhood of vertices. The one-ring neighborhood of four vertices
are drawn with different colors, the corresponding vertex is colored accordingly.
Two adjacent vertices are in the one-ring neighborhood of each other.
\cgalFigureEnd
There are two unknowns per vertex in the Eq. \f$\eqref{eq:arap_energy}\f$, the new positions (\f$\mathbf{v}'_k\f$) of the free
vertices and the rotation matrices (\f$\mathbf{R}_i\f$).
If the energy contribution of each vertex is positive, this boils down to minimizing the energy contribution of each vertex \f$\mathbf{v}_i\f$.
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 unknown.
For the vertex \f$\mathbf{v}_i\f$, we consider the covariance matrix \f$\mathbf{S}_i\f$:
\f[
\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}
\f]
It has been proven \cgalCite{Sorkine2009LeastSquaresRigid} that minimizing the energy contribution of
\f$\mathbf{v}_i\f$ in Eq. \f$\eqref{eq:arap_energy}\f$ is equivalent to maximizing the trace of the matrix
\f$\mathbf{R}_i \mathbf{S}_i\f$. \f$\mathbf{R}_i \f$ is the transpose of the unitary matrix in the polar decomposition
of \f$\mathbf{S}_i\f$.
In the second step, the rotation matrices are substituted into the partial derivative of Eq.\f$\eqref{eq:arap_energy}\f$
with respect to \f$\mathbf{v}'_i\f$. Assuming the weights are symmetric, setting the derivative to zero results in the
following equation:
\f[
\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}
\f]
The left-hand side of this equation corresponds to the one of Eq.\f$\eqref{eq:lap_open}\f$,
and we can set \f$\Delta\f$ to be the right-hand side.
Solving the linear system in Eq. \f$\eqref{eq:lap_energy_system}\f$ gives the new positions of the free vertices.
This two-step optimization can be applied several times iteratively to obtain a better result.
\note The matrix built with the Laplacian matrix of the free vertices in the left-hand side of
Eq. \f$\eqref{eq:lap_energy_system}\f$ depends only on the initial mesh structure and on which vertices
are control vertices. Once the control vertices are set, we can use a direct solver to factorize
the sparse matrix in Eq. \f$\eqref{eq:lap_energy_system}\f$, and reuse this factorization during
each iteration of the optimization procedure.
The original algorithm \cgalCite{Sorkine2007AsRigidAs} we described assumes that:
- the weight between two vertices is symmetric. In order to support asymmetric weights in our implementation,
we slightly changed Eq. \f$\eqref{eq:lap_ber}\f$ to:
\f[
\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}
\f]
- The energy contribution of each vertex is positive. If the weight between two vertices is always positive, this is always the case.
However when using the cotangent weighting scheme (the default in our implementation), if the sum of the angle opposite to an edge is greater than \f$ \pi \f$,
its cotangent weight is negative. As a workaround for bad quality meshes, we eliminate those negative weights
by setting them to zero.
A slightly modified version is described in the next section to avoid the latter issue.
\subsection Surface_Modeling_Overview_ARAP_Rims Spokes and Rims Version
The modified energy function proposed by \cgalCite{Chao2010SimpleGeomModel} takes also into account
all the opposite edges in the facets incident to a vertex. The energy function to minimize becomes:
\f[
\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}
\f]
where \f$E(\mathbf{v}_i)\f$ consists of the set of edges incident to \f$\mathbf{v}_i\f$ (the <em>spokes</em>) and
the set of edges in the link (the <em>rims</em>) of \f$\mathbf{v}_i\f$ in the surface mesh \f$M\f$
(see \cgalFigureRef{Spoke_and_rim_edges}).
\cgalFigureBegin{Spoke_and_rim_edges, spoke_and_rim_edges_2.png}
The vertices \f$ \mathbf{v}_n\f$ and \f$ \mathbf{v}_m\f$ are the opposite vertices to the edge \f$
(\mathbf{v}_i,\mathbf{v}_j)\f$.
\cgalFigureEnd
The method to get the new positions of the free vertices is similar to the two-step optimization
method explained in \ref Surface_Modeling_Overview_ARAP.
For the first step, the Eq. \f$\eqref{eq:cov_matrix}\f$ is modified to take into account the edges in \f$E(\mathbf{v}_i)\f$:
\f[
\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}
\f]
For the second step, setting partial derivative of Eq. \f$\eqref{eq:arap_energy_rims}\f$ to zero with
respect to \f$\mathbf{v}_i\f$ gives the following equation:
\f[
\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}
\f]
where \f$\mathbf{R}_m\f$ and \f$\mathbf{R}_n\f$ are the rotation matrices of the vertices \f$\mathbf{v}_m\f$,
\f$\mathbf{v}_n\f$ which are the opposite vertices of the edge \f$(\mathbf{v}_i, \mathbf{v}_j)\f$
(see \cgalFigureRef{Spoke_and_rim_edges}). Note that if the edge \f$ (\mathbf{v}_i \mathbf{v}_j) \f$ is on
the boundary of the surface mesh, then \f$ w_{ij} \f$ must be 0 and \f$ \mathbf{v}_m
\f$ does not exist.
An important property of this approach compared to \ref Surface_Modeling_Overview_ARAP is that the contribution to the global energy
of each vertex is guaranteed to be non-negative when using the cotangent weights \cgalCite{Chao2010SimpleGeomModel}.
We are using cotangent weights by default (negative values included) as proposed in \cgalCite{Chao2010SimpleGeomModel}.
A drawback of this method is that the results can be asymmetric as shown on \cgalFigureRef{Arap_spokes_comparison}.
\section Surface_Modeling_History Design and Implementation History
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 optimally using the Eigen library.
*/
} /* namespace CGAL */