Once built, you can execute the assignment from inside the build/ by running
+on a given mesh:
+
+
./deformation [path to mesh.obj]
+
+
+
When the mesh is blue, the system is in “place handles” mode. Click on the mesh
+to select vertex locations for control point handles. After pressing space to switch to
+deformation mode, drag the handles. Pressing m will switch between different
+deformation methods.
+
+
+
+
+
+
Background
+
+
In this assignment we explore smooth deformation of an existing shape. Shape
+deformation has many
+applications in geometry process; we will focus on interactive handle-based
+deformation. In this setup, the user repositions a sparse set of points and the
+goal is to propagate the transformation at these “handles” to the rest of the
+shape. To be interactive we should aim for computing the new deformation of the
+shape at around 30 frames per
+second.
+
+
Shape deformation is the transformation from its rest shape to a
+new/current/deformed shape. If the position of a point on some 3D rest shape
+ given by \(\hat{\x} ∈ \R³\) then we will write that the unknown position
+on the deformed shape is given by \(\x ∈ \R³\). We can write this point’s
+displacement vector as \(\x - \hat{\x} =: \d ∈ \R³\).
+
+
The propagation of the handles’ deformation can be thought of in two
+complementary ways:
as a shape optimization problem, where we try to define a new shape that
+ retains the details of the old shape but fulfills the handle constraints.
+
+
+
In the following discussion we will take advantage of the ability to switch
+between thinking of the unknowns as the positions of the deformed shape (\(\x\))
+and displacements (\(\d\)). These views are equivalent, but often one or the
+other provides a better intuitive understanding.
+
+
Continuity
+
+
We will limit ourselves to continuous deformations of shapes. That is, the
+shape will not tear, crack or change its topological features.
+
+
If we represent our shape discretely as a triangle
+mesh (e.g., with rest vertices
+in \(\hat{\V}
+∈ \R^{n × 3}\) and faces in \(F ∈ \{1,…,n\}^{m × 3}\), then we can trivially
+ensure a continuous deformation by determining new vertex positions \(\V\). The
+topology (connectivity) of the mesh (\(F\)) will not change.
+
+
Generic Distortion Minimization
+
+
A rest surface \(\hat{S}\)
+immersed in \(\R³\) can
+be described as a mapping \(\hat{\x}\) from some 2D parametric domain \(Ω\). For any
+parameters \(u\) and \(v\), \(\hat{\x}\) describes the 3D position:
+
+
\[
+\hat{\x}(u,v) ∈ \R³.
+\]
+
+
Similarly the deformed surface can be represented as a position function \(\x:
+Ω → \R³\). The displacement vector field is thus a function taking the
+difference: \(\d(u,v) = \x(u,v) - \hat{\x}(u,v)\).
+
+
For the handle-based deformation problem we would like to find a new surface
+(defined by \(\x\)) that:
+
+
+
adds as little distortion as possible to the shape, and
+
satisfies the users constraints at selected handle positions.
+
+
+
We can cast this as an energy optimization problem. Suppose we have
+energy functional
+\(E(\x)\) that
+measures the amount of distortion between the new shape (\(\x\)) and the rest
+shape \(\hat{x}\), then we could optimize for the best possible shape \(\x\) by
+minimizing \(E\):
+
+
\[
+\min_\x E(\x).
+\]
+
+
To ensure that the user’s \(k\) handle points are interpolated we add the
+constraints:
+
+
\[
+\text{ subject to } \x(u_i, v_i) = \g_i \ ∀ i = \{1, … , k\},
+\]
+where \(\g_i\) is the position of the $i$th control point handle.
+
+
While the constraints are straightforward, we have many choices for how to
+formulate the energy function \(E\). A natural choice is to measure distortion in
+an egalitarian way by integrating a local measure of distortion at all points
+on the surface:
+
+
\[
+\min_\x ∫_Ω ‖ e(\x) ‖² \ dA \quad \text{ subject to } \x(u_i, v_i) = \g_i \ ∀ i = \{1, … , k\},
+\]
+
+
where \(e\) is a vector- or scalar- valued function measuring local (unsquared)
+distortion. We will now consider different choices for \(e\).
+
+
Linear Methods
+
+
If we assume that the deformation between the rest shape given by \(\hat{\x}\)
+and the new shape given by \(\x\) is small then we can measure the distortion
+of the deformation in terms of the smoothness of the displacement field. This
+simplest methods will integrate the magnitude of derivatives of the
+displacement field (\(\d\)): if the displacement field has large variations or
+sudden changes then it is inducing a lot of distortion.
+
+
Gradient-based energy
+
+
Let us first consider minimizing the integral of squared variation of the
+displacement field:
+
+
\[
+\min_\d ∫_Ω ‖ ∇\d ‖_F^2 \ dA \quad \text{ subject to } \d_i =
+\g_i-\hat{\x}_i \ ∀ i = \{1, … , k\},
+\]
+where
+\(∇\d = \left(\begin{array}{ccc}
+\frac{∂d^x}{∂u} & \frac{∂d^y}{∂u} & \frac{∂d^z}{∂u} \\
+\frac{∂d^x}{∂v} & \frac{∂d^y}{∂v} & \frac{∂d^z}{∂v}
+\end{array}\right)\)
+the Jacobian
+matrix of the displacement field \(\d\)
+
+
+
Deformation Gradient
+
+
If \(\I ∈ \R^{3 × 3}\) is the identity matrix, then the quantity \(\F := \I +
+∇\d\) is referred to as the deformation
+gradient
+in the mechanics
+community.
+
+
+
This is simply the familiar Dirichlet
+energy applied to each
+coordinate function of the displacement field independently.
+
+
We can discretize this over our triangle mesh surface the same way we have in
+smoothing and parameterization assignments:
+
+
\[
+\min_{\D} \tr{\D^\transpose \L \D} \quad \text{subject to } \D_\text{handles} =
+\g_\text{handles} - \hat{\V}_\text{handles},
+\]
+where the rows of \(\g_\text{handles} ∈ \R^{k × 3}\) contains the new positions
+of the \(k\) control point handles.
+
+
While easy to implement, this method suffers from a couple immediate problems:
+
+
+
it is not smooth at constraints, and
+
then influence of handles dies off too quickly.
+
+
+
+
+Two regions of selected points (purple) are constrained. The gradient-based
+energy produces a continuous displacement but is not smooth at the inner
+constrained points (sharp crease). The Laplacian-based energy produces a
+displacement with continuous positions and
+derivatives.
+
+
+
By minimizing the Dirichlet energy, each coordinate of “displacement field” is
+a harmonic function.
+Intuitively (however abstractly) we can think of each function as diffusing
+the user’s constraints as if they were
+heat values. As such, each
+function diffuses quickly to an average, slowly varying value over most of the
+domain. As a displacement field a constant value for each coordinate function
+would mean a translation: most of the shape is simply translated.
+
+
+
+
+
+
The gradient operator (\(∇\)) is a linear
+operator. We can alternatively view
+our minimization above in terms of the unknown positions \(\x\):
+
+
\[
+\min_\d ∫_Ω ‖ ∇\d ‖_F^2 \ dA ⇒
+\min_\x ∫_Ω ‖ ∇(\x - \hat{\x}) ‖_F^2 \ dA ⇒
+\min_\x ∫_Ω ‖ \underbrace{∇\x}_\text{after} -
+\underbrace{∇\hat{\x}}_\text{before} ‖_F^2 \ dA.
+\]
+If we think of the gradient of the position function \(∇\x\) (with respect to the
+underlying parameterization \(u,v\)) as a local geometric feature
+descriptor then this
+energy can be re-understood as measuring the difference in this feature before
+and after the deformation. This is very sensible as we are trying to measure
+distortion. We would expect that a low-distortion deformation would correspond
+with a small change to local features.
+
+
Unfortunately, the gradient of the position function \(\x\) is a poor,
+first-order local feature.
+
+
Laplacian-based energy
+
+
If we model distortion as the change in a local feature descriptor, then a
+natural local and relative descriptor would be one that compared the position
+of some point on the shape to the average of its local neighborhood. We have
+studied an operator that computes this in the smoothing assignment. The
+Laplace(-Beltrami) operator
+can be derived as taking exactly the difference of a functions value at a point
+and the average (i.e., centroid) of
+an infinitesimal region around that point:
(see, e.g., “Differential coordinates for local mesh morphing and deformation”
+[Alexa et al. 2003] and expanded upon in “Laplacian Surface Editing”
+[Sorkine et al. 2004]).
+
+
When applied to the embedding function \(\x\) the Laplace operator computes the
+difference in position between a point and its local neighborhood. This
+vector points in the
+normal direction and its
+magnitude corresponds inversely with how flat the surface is
+locally.
+
+
Let’s replace the gradient feature above with this second-order feature
+descriptor and massage our optimization problem back in terms of
+displacements:
+
+
\[
+\min_\x ∫_Ω ‖ \underbrace{∆\x}_\text{after} -
+\underbrace{∆\hat{\x}}_\text{before} ‖^2 \ dA ⇒
+\min_\x ∫_Ω ‖ ∆(\x - \hat{\x}) ‖^2 \ dA ⇒
+\min_\d ∫_Ω ‖ ∆\d ‖^2 \ dA.
+\]
+
+
Just as we can show that harmonic functions (\(∆\d = 0\)) minimize the Dirichlet
+energy, we can use calculus of
+variations apply
+Green’s identitytwice to
+show that minimizers of the squared-Laplacian energy are _bi-_harmonic
+functions (\(∆∆\d = 0\) or \(∆²\d = 0\)). Obviously all harmonic functions are also
+biharmonic functions. This implies that the space of biharmonic functions is
+strictly larger. In particular, this will allow use to ensure continuity of
+first derivatives across constrained values at handles.
+
+
+
+
+
+
The fact that each coordinate of the displacement field \(\d\) will be a
+bi-harmonic function is not so elucidating, but by minimizing the
+squared-Laplacian energy integrated over the domain, we can say that it is
+as-harmonic-as-possible. Harmonic functions include constant functions (i.e.,
+translations) but also any displacement that bends in one direction by an
+equal and opposite amount as it bends in the other direction: \(∆\d = 0 → ∂\d/∂u
++ ∂\d/∂v = 0 → ∂\d/∂u = -∂\d/∂v\).
+
+
To discretize this energy we can make use of our discrete Laplacian operator
+\(\L\). This matrix computes the locally integrated Laplacian of a given
+function specified by per-vertex values \(\f\). Now we would like to integrate
+the square of the point-wise Laplacian. We can approximate the point-wise
+Laplacian by the local integral
+average of the Laplacian
+\(\M^{-1}\L \f\). Integrating this over the mesh we have our complete approximate
+of the energy:
+
+
\[
+\begin{align}
+∫_Ω ‖∆\d‖² \;dA &≈ \tr{ \D^\transpose \L^\transpose \M^{-\transpose} \M \M^{-1} \L
+\D}\\
+&= \tr{ \D^\transpose \underbrace{\L^\transpose \M^{-1} \L}_{\Q} \D },
+\end{align}
+\]
+where \(\M ∈ \R^{n × n}\) is the mass-matrix for the given mesh and \(\Q ∈
+\R^{n×n}\) can be thought of as the bi-Laplacian matrix.
+
+
+
\(k\)-harmonic
+
+
The logical continuation of harmonic and biharmonic deformation is to consider
+triharmonic (\(∆³\d = 0\)) and tetraharmonic (\(∆⁴\d = 0\)) and so on. It’s
+straightforward to implement these, though there are diminishing returns and
+increasing costs.
+
+
+
Precomputation
+
+
With out loss of generality, assume that the rows of the unknown displacements
+\(\D\) have been sorted so that displacements corresponding to handle vertices
+are in the bottom part:
If we don’t change which vertices are handles, but only change the positions
+of the selected handles, then only \(\D_\text{h}\) changes above. In particular,
+the matrix \(\Q_\text{u,u}\) is unchanged. Therefore, we can
+prefactorize it so that
+applying its inverse is fast (igl::min_quad_with_fixed does this for you).
+
+
+
Actually the entire term \(\Q_\text{u,u}^{-1} \Q_\text{u,h} =: \W ∈ \R^{(n-k)
+× k}\) does not change. The columns of \(\W\) reveal how unknown displacements
+respond to each handle point’s displacement (see also “An intuitive framework
+for real-time freeform modeling” [Botsch & Kobbelt 2004]). This provides a
+gateway to the relationship with linear blend skinning and automatic
+(biharmonic) weighting functions (see “Bounded Biharmonic Weights for
+Real-Time Deformation” [Jacobson et al. 2011]).
+
+
+
Trouble in paradise
+
+
Biharmonic displacements work well for small deformations and deformations that
+do not imply a large rotation. However, the Laplacian of the position function
+\(∆\x\) as a feature descriptor is not rotation invariant. This problem is true
+for all linear differential features including the gradient of the embedding
+function \(∇\x\) considered above (see “On linear variational surface deformation
+methods” [Botsch & Sorkine 2008]).
+
+
+
+
+
+
This means that if the user transforms all of the handle locations by a rigid
+transformation \(\T\) these energies will not measure zero for a displacement
+equivalent to applying the rigid transformation \(\T\) to the entire shape. We
+would like global rotation invariance, but we would also like this property
+to apply locally so that parts of the shape can easily rotate.
+
+
As-rigid-as-possible
+
+
In the scenario where each handle
+\(i\) are perfectly transformed by a rigid
+transformation\(\x_i =
+\Rot \hat{\x}_i + \t\), where \(\Rot ∈ SO(3) ⊂ \R^{3×3}\) is a rotation matrix and
+\(\t∈\R^3\) is a translation vector. If an
+oracle could only tell us this
+particular rigid transformation then we could repair
+the gradient-based energy above by pre-rotating the rest shape by this
+transformation:
where the translation vector \(\t\) falls out because a translation has constant
+gradient.
+
+
We do not know the rotation \(\Rot\) ahead of time, but we could be as generous
+as possible and use the “best” rotation \(\Rot ← \argmin_\Rot ∫_Ω ‖ ∇\x - \Rot
+∇\hat{\x} ‖² \;dA\):
Optimizing this energy will ensure global rotation invariance. To ensure
+local rotation invariance, we can replace \(\Rot ∈ SO(3)\) with a spatially
+varying function\(\Rot : Ω → SO(3)\) that outputs a “best” rotation for any
+point on the shape (see “A simple geometric model for elastic deformations”
+[Chao et al. 2010]). In this way, the optimal rotation will be locally rigid
+everywhere, or as-rigid-as-possible (ARAP).
+
+
+
+
+
+
+
For embedded solid shapes, we can take the rest shape given by \(\hat{\x}\) as
+the parameterization \(Ω\), so that \(∇\hat{\x} = \I\). This allows us to rewrite
+the as-rigid-as-possible energy as the square of the difference between the
+deformation gradient and the closest rotation:
This form provides a bridge between the as-rigid-as-possible energy common in
+geometry processing to corotated linear elasticity used in
+mechanics/physically-based simulation (made explicit in “A simple geometric
+model for elastic deformations” [Chao et al. 2010]). See Section 3.4 of “FEM
+Simulation of 3D Deformable Solids” [Sifakis 2012] for a graphics-mechanics
+perspective.
+
+
+
Discrete as-rigid-as-possible energy
+
+
For a triangle mesh with displacing vertices, the gradient of the embedding
+function is constant inside each triangle. In this way we can write the raw
+gradient energy above as a double sum over all half-edges \(ij\) of all faces \(f\)
+in the mesh:
+
+
\[
+½ ∫_Ω ‖ ∇ \x - ∇\hat{\x}‖² \;dA = ½ ∑\limits_{f ∈ F} ∑\limits_{ ij ∈ f} c_{ij} ‖
+(\v_i-\v_j) - (\hat{\v}_i-\hat{\v}_j)‖²,
+\]
+where \(c_{ij}\) is cotangent of the angle opposite half-edge \(ij\).
+
+
To inject localized best fit rotations, we will assign an unknown rotation
+matrix \(\Rot_k\) to each vertex \(k\) of the mesh and accounts for a third of the
+energy integrated over incident triangles:
+\[
+½ ∫_Ω ‖ ∇ \x - \Rot ∇\hat{\x}‖² \;dA =
+⅙ ∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)}
+c_{ij} ‖ (\v_i-\v_j) - \Rot_k (\hat{\v}_i-\hat{\v}_j)‖²,
+\]
+where \(F(k)\) is the set of all faces incident on the \(k\)-th vertex.
+
+
+
Where do rotations live?
+
+
We have assigned a rotation for each vertex: there are \(n\) rotation matrices
+as auxiliary degrees of freedom. This is in contrast to assigning
+rotations per face (as in “A Local/Global Approach to Mesh Parameterization”
+[Liu et al. 2008]). Per-face–or more generally per-element–rotations work
+well for codimension zero
+objects (triangle meshes in \(\R²\) or tetrahedral meshes in \(\R³\)). But for
+triangle mesh surfaces in \(\R³\) (i.e., codimension 1), per-face rotations
+would lead to “crumpling”) because bending along edges would not be
+measured.
+
+
+
Optimization
+
+
The simplest method for optimizing the ARAP energy is by alternating between
+
+
+
finding the optimal rotations \(\Rot_k\) assuming the vertex positions \(\V\) are
+ fixed, and
+
finding the optimal vertex positions \(\V\) assuming all rotations \(\Rot_k\) are
+ fixed.
+
+
+
Each rotation \(\Rot_k\) only affects the local energy and doesn’t interact with
+the other rotations. So each can be optimized locally. In contrast, the
+mesh vertex positions \(\V\) depend on each other requiring a global solve. In
+the geometry processing literature, this is known as a local-global
+optimization (see “As-rigid-as-possible surface modeling” [Sorkine & Alexa
+2007]). It is also known as “alternating block coordinate descent” because we
+have separated the variables into disjoint sets (\(\V,\Rot_1,…,\Rot_n\)) and taking
+the optimal descent direction for each independently.
+
+
Observing the discrete energy above we can see that the energy is quadratic in
+\(\V\) and quadratic in each \(\Rot_k\). Let’s start by separating the terms that are
+quadratic and linear in \(\V\):
if we stack the rotation matrices \(\Rot_k\) into large matrix \(\Rot ∈ \R^{3n ×
+3}\) then we can write this energy in matrix form as:
+
+
\[
+\tr{ \V^\transpose \L \V } + \tr{ \V^\transpose \K \Rot },
+\]
+where \(\L ∈ \R^{n × n}\) is the familiar cotangent discrete Laplacian matrix and
+\(\K ∈ \R^{n × 3n}\) sparse matrix containing cotangents multiplied against
+differences across edges in the rest mesh (e.g., \(\hat{\v}_i - \hat{\v}_j\)).
+
+
+
I’m so confused. What’s in the \(\K\) matrix?
+
+
Let’s take it slow. The \(\K\) matrix is represents the
+bilinear form that combines unknown
+vertex positions and unknown rotations. We have identified above that we can
+write this in summation or matrix form:
+\[
+⅙ ∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)} c_{ij} (\v_i-\v_j)^\transpose \Rot_k (\hat{\v}_i-\hat{\v}_j)
+= \tr{ \V^\transpose \K \Rot },
+\]
+but how did we get here?
+
+
Let’s start with the summation form. The constants of this formula are the
+\(c_{ij}\) terms and the \((\hat{\v}_i-\hat{\v}_j)\) terms. Since these always
+appear together, let us merge them into weighted edge difference vectors
+\(c_{ij} (\hat{\v}_i-\hat{\v}_j) =: \hat{\e}_{ij} ∈ \R³\):
+
+
\[
+⅙ ∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)} \underbrace{(\v_i-\v_j)^\transpose
+\Rot_k \hat{\e}_{ij}}_{∈\R},
+\]
+the inner term in the summation is an inner
+product; that is, a
+scalar. Let’s expose this
+by expanding the matrix-vector products of the inner-product:
+\[
+⅙ ∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)} ∑\limits_{α=1}^3 ∑\limits_{β=1}^3
+ (v_i^α - v_j^α)R_k^{αβ}\hat{e}_{ij}^β.
+\]
+
+
If our mesh is stored as a vertex list and face list, it’s not easy/efficient
+to loop over per-vertex rotations (outer sum) and then over all half-edges of
+incident faces (second sum). Instead, let’s rearrange these sums to loop over
+all faces first, then the half-edges of that face, and then over all
+per-vertex rotations that involve this half-edge:
+\[
+⅙ ∑\limits_{f=1}^m ∑\limits_{ ij ∈ E(f)} \quad ∑\limits_{k | ij ∈F(k)} \quad ∑\limits_{α=1}^3 ∑\limits_{β=1}^3
+ (v_i^α - v_j^α)R_k^{αβ}\hat{e}_{ij}^β,
+\]
+where the third sum is over all rotations \(k\) such that the half-edge \(ij\)
+belongs to the half-edges of the faces incident on the \(k\)-th vertex: \(k | ij
+∈F(k)\). Well, this means \(k\) can either be \(i\) or \(j\) or the third vertex of
+the \(f\)-th face.
+
+
Now let’s turn our attention back to the
+summand. The terms indexed by \(α\)
+never mix. That is, we never add/multiply \(v_i^α\), \(v_j^γ\), and \(R_k^{δβ}\)
+unless \(α=γ=δ\). This implies that we can write this summation in matrix form
+as:
+\[
+\V₁^\transpose \K₁ \Rot₁ +
+\V₂^\transpose \K₂ \Rot₂ +
+\V₃^\transpose \K₃ \Rot₃,
+\]
+where \(\V_α ∈ \R^n\) is \(α\)-th column of \(\V\), \(\Rot_α ∈ \R^{3n}\) is the \(α\)-th column of
+\(\Rot\) and \(\K₁,\K₂,\K₃ ∈ \R^{n × 3n}\) are sparse matrices.
+
+
Further, the constant term \(\hat{e}^β_{ij}\) in the summation acts the same on
+\((v_i^α - v_j^α)R_k^{αβ}\) for any \(α\) value. This implies that \(\K₁=\K₂=\K₃\),
+so we can reduce the matrix form to:
+\[
+\tr{\V \K \Rot}.
+\]
+
+
Finally, we can answer what is in each entry \(K_{vw}\), where \(v\) chooses the
+row and \(w=3k+β\) chooses the column of \(\K\). Treating our nested summation as
+nested for-loops, we can increment entries in \(\K\)
+\[
+\begin{align}
+K_{i\ 3k+β} & \mathrel{+}= \hat{e}_{ij}^β, \\
+K_{j\ 3k+β} & \mathrel{+}= -\hat{e}_{ij}^β, \\
+\end{align}
+\]
+for each half-edge encountered.
+
+
+
Local step
+
+
Minimizing this energy with respect \(\Rot\) corresponds to minimizing:
+
+
\[
+\tr{ \underbrace{\V^\transpose \K}_{\C^\transpose} \Rot },
+\]
+where \(\C ∈ \R^{3n × 3}\) stacks weighted covariance matrices \(\C_k ∈ \R^{3 ×
+3}\) for each region covered by the corresponding rotation \(\Rot_k\). We have
+seen this problem before in the registration assignment. For each \(\C_k\),
+\(\Rot_k\) will be the closest rotation matrix solved via singular value
+decomposition.
+
+
Global step
+
+
Minimizing the energy above with respect to \(\V\) corresponds to solving a
+Dirichlet energy-like minimization problem:
where \(\K \Rot =: \B ∈ \R^{n × 3}\) is a matrix of rotated vertex gradients.
+Adding the handle constraints to the corresponding rows of \(\V\) this is easily
+minimized by setting all partial derivatives with respect to the unknowns in
+\(\V\) equal to zero (as in the linear methods above).
+
+
Implementation
+
+
In order to facilitate interactive deformation we would like our local and
+global iterations to be computed as quickly as possible. Since the quadratic
+form in the global step is the same regardless of the current rotations or
+current handle positions, we can
+prefactorize it (again,
+as above). The matrix \(\K\) also does not depend on the rotations, current
+positions or handle positions, so we can pre-build this matrix. This way
+computing \(\C\) and \(\B\) is a simple matrix-matrix multiplication.
+
+
+
Note: When constructing \(\K\) it’s easiest to iterate over all
+half-edges in the mesh (by iterating over all faces and then each of the
+three edges). Each half-edge \(ij\)contributes terms tying \(\v_i,\v_j\) to
+each of the (three) rotations \(\Rot_k\) that apply against their difference
+(see “Fast Automatic Skinning Transformations” [Jacobson et al. 2012]).
+
+
+
Still trouble in paradise
+
+
+
+
+
+
The as-rigid-as-possible deformation method for surfaces described above has a
+number of remaining problems:
+
+
+
like its gradient-based energy ancestor the deformation is not smooth at
+ constraints;
+
the energy punishes bending of the surface–which is good–but does so in
+ a way that diminishes as the mesh becomes higher and higher resolution, in
+ otherwords, the discrete energy is mesh-resolution dependent; and
+
the energy is biased by the original combinatorics of the mesh (even in
+ flat regions).
+
+
+
Tasks
+
+
Blacklist
+
+
+
igl::arap
+
igl::arap_linear_block
+
igl::covariance_scatter_matrix
+
igl::harmonic
+
+
+
Whitelist
+
+
+
igl::polar_svd3x3 (or your previous assignment’s closest_rotation)
+
igl::min_quad_with_fixed
+
igl::cotmatrix_entries
+
igl::cotmatrix (or your previous implementation)
+
igl::massmatrix (or your previous implementation)
+
+
+
src/biharmonic_precompute.cpp
+
+
Precompute data needed to efficiently solve for a biharmonic deformation given
+a mesh with vertices V and faces F and a list of selected vertices as
+indices b into V. The output should be a prefacorized system using the
+data struct employed by igl::min_quad_with_fixed.
+
+
src/biharmonic_solve.cpp
+
+
Given precomputation data and a list of handle displacements determine
+displacements for all vertices in the mesh.
+
+
src/arap_precompute.cpp
+
+
Precompute data needed to efficiently conduct local-global iterations for an
+arap deformation. This includes the data struct employed by
+igl::min_quad_with_fixed to solve the global step" and constructing the
+bilinear form K that mixes rotation degrees of freedom with unknown
+positions for preparing the covariance matrices of the local step and the
+linear term of the global step.
+
+
src/arap_single_iteration.cpp
+
+
Given precomputed data (data and K), handle positionsbc and current
+positions of all vertices U, conduct a single iteration of the local-global
+solver for minimizing the as-rigid-as-possible energy. Output the positions
+of all vertices of the mesh (by overwriting U).
+
+
+
diff --git a/README.md b/README.md
index 5222e7c..6a673e3 100644
--- a/README.md
+++ b/README.md
@@ -17,80 +17,682 @@ on a given mesh:
./deformation [path to mesh.obj]
+When the mesh is blue, the system is in "place handles" mode. Click on the mesh
+to select vertex locations for control point handles. After pressing space to switch to
+deformation mode, drag the handles. Pressing `m` will switch between different
+deformation methods.
+
+
+
## Background
+In this assignment we explore smooth deformation of an existing shape. [Shape
+deformation](https://en.wikipedia.org/wiki/Deformation_(mechanics)) has many
+applications in geometry process; we will focus on interactive **_handle-based
+deformation_**. In this setup, the user repositions a sparse set of points and the
+goal is to propagate the transformation at these "handles" to the rest of the
+shape. To be interactive we should aim for computing the new deformation of the
+shape at around 30 [frames per
+second](https://en.wikipedia.org/wiki/Frame_rate).
+
+Shape deformation is the transformation from its _rest shape_ to a
+new/current/deformed shape. If the position of a point on some 3D rest shape
+ given by $\hat{\x} ∈ \R³$ then we will write that the unknown position
+on the deformed shape is given by $\x ∈ \R³$. We can write this point's
+displacement vector as $\x - \hat{\x} =: \d ∈ \R³$.
+
+The propagation of the handles' deformation can be thought of in two
+complementary ways:
+
+ 1. as a [scattered data
+ interpolation](https://en.wikipedia.org/wiki/Multivariate_interpolation#Irregular_grid_.28scattered_data.29)
+ problem, where handles provide sparse samples of an unknown [displacement
+ field](https://en.wikipedia.org/wiki/Displacement_field_(mechanics)), or
+ 2. as a shape optimization problem, where we try to define a _new_ shape that
+ retains the details of the old shape but fulfills the handle constraints.
+
+In the following discussion we will take advantage of the ability to switch
+between thinking of the unknowns as the positions of the deformed shape ($\x$)
+and displacements ($\d$). These views are equivalent, but often one or the
+other provides a better intuitive understanding.
+
+#### Continuity
+
+We will limit ourselves to _continuous_ deformations of shapes. That is, the
+shape will not tear, crack or change its topological features.
+
+If we represent our shape discretely as a [triangle
+mesh](https://en.wikipedia.org/wiki/Triangle_mesh) (e.g., with _rest_ vertices
+in $\hat{\V}
+∈ \R^{n × 3}$ and faces in $F ∈ \{1,…,n\}^{m × 3}$, then we can _trivially_
+ensure a continuous deformation by determining new vertex positions $\V$. The
+topology (connectivity) of the mesh ($F$) will not change.
+
+#### Generic Distortion Minimization
+
+A rest surface $\hat{S}$
+[immersed](https://en.wikipedia.org/wiki/Immersion_(mathematics)) in $\R³$ can
+be described as a mapping $\hat{\x}$ from _some_ 2D parametric domain $Ω$. For any
+parameters $u$ and $v$, $\hat{\x}$ describes the 3D position:
+
+\\[
+\hat{\x}(u,v) ∈ \R³.
+\\]
+
+Similarly the deformed surface can be represented as a position function $\x:
+Ω → \R³$. The displacement vector field is thus a function taking the
+difference: $\d(u,v) = \x(u,v) - \hat{\x}(u,v)$.
+
+For the handle-based deformation problem we would like to find a new surface
+(defined by $\x$) that:
+
+ 1. adds as little _distortion_ as possible to the shape, and
+ 2. satisfies the users constraints at selected handle positions.
+
+We can cast this as an energy optimization problem. Suppose we have
+energy [functional](https://en.wikipedia.org/wiki/Functional_(mathematics))
+$E(\x)$ that
+measures the amount of distortion between the new shape ($\x$) and the rest
+shape $\hat{x}$, then we could optimize for the best possible shape $\x$ by
+minimizing $E$:
+
+\\[
+\min_\x E(\x).
+\\]
+
+To ensure that the user's $k$ handle points are interpolated we add the
+constraints:
+
+\\[
+\text{ subject to } \x(u_i, v_i) = \g_i \ ∀ i = \{1, … , k\},
+\\]
+where $\g_i$ is the position of the $i$th control point handle.
+
+While the constraints are straightforward, we have many choices for how to
+formulate the energy function $E$. A natural choice is to measure distortion in
+an egalitarian way by integrating a _local_ measure of distortion at all points
+on the surface:
+
+\\[
+\min_\x ∫_Ω ‖ e(\x) ‖² \ dA \quad \text{ subject to } \x(u_i, v_i) = \g_i \ ∀ i = \{1, … , k\},
+\\]
+
+where $e$ is a vector- or scalar- valued function measuring local (unsquared)
+distortion. We will now consider different choices for $e$.
+
### Linear Methods
-Measure distortion with respect to small displacements
+If we assume that the deformation between the rest shape given by $\hat{\x}$
+and the new shape given by $\x$ is _small_ then we can measure the distortion
+of the deformation in terms of the smoothness of the displacement field. This
+simplest methods will integrate the magnitude of derivatives of the
+displacement field ($\d$): if the displacement field has large variations or
+sudden changes then it is inducing a lot of distortion.
#### Gradient-based energy
-- Dependent on global coordinate system
+Let us first consider minimizing the integral of squared variation of the
+displacement field:
+
+\\[
+\min_\d ∫_Ω ‖ ∇\d ‖_F^2 \ dA \quad \text{ subject to } \d_i =
+\g_i-\hat{\x}_i \ ∀ i = \{1, … , k\},
+\\]
+where
+$∇\d = \left(\begin{array}{ccc}
+\frac{∂d^x}{∂u} & \frac{∂d^y}{∂u} & \frac{∂d^z}{∂u} \\
+\frac{∂d^x}{∂v} & \frac{∂d^y}{∂v} & \frac{∂d^z}{∂v}
+\end{array}\right)$
+the [Jacobian](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant)
+matrix of the displacement field $\d$
+
+> ##### Deformation Gradient
+>
+> If $\I ∈ \R^{3 × 3}$ is the identity matrix, then the quantity $\F := \I +
+> ∇\d$ is referred to as the [deformation
+> gradient](https://en.wikipedia.org/wiki/Finite_strain_theory#Deformation_gradient_tensor)
+> in the [mechanics](https://en.wikipedia.org/wiki/Continuum_mechanics)
+> community.
+
+This is simply the familiar [Dirichlet
+energy](https://en.wikipedia.org/wiki/Dirichlet's_energy) applied to each
+coordinate function of the displacement field _independently_.
+
+We can discretize this over our triangle mesh surface the same way we have in
+smoothing and parameterization assignments:
+
+\\[
+\min_{\D} \tr{\D^\transpose \L \D} \quad \text{subject to } \D_\text{handles} =
+\g_\text{handles} - \hat{\V}_\text{handles},
+\\]
+where the rows of $\g_\text{handles} ∈ \R^{k × 3}$ contains the new positions
+of the $k$ control point handles.
+
+While easy to implement, this method suffers from a couple immediate problems:
+
+ 1. it is not smooth at constraints, and
+ 2. then _influence_ of handles dies off too quickly.
+
+
+
+By minimizing the Dirichlet energy, each coordinate of "displacement field" is
+a [harmonic function](https://en.wikipedia.org/wiki/Harmonic_function).
+Intuitively (however abstractly) we can think of each function as _diffusing_
+the user's constraints as if they were
+[heat](https://en.wikipedia.org/wiki/Heat_equation) values. As such, each
+function diffuses quickly to an average, slowly varying value over most of the
+domain. As a displacement field a constant value for each coordinate function
+would mean a translation: most of the shape is simply translated.
+
+
+
+The gradient operator ($∇$) is a [_linear_
+operator](https://en.wikipedia.org/wiki/Linear_map). We can alternatively view
+our minimization above in terms of the unknown positions $\x$:
+
+\\[
+\min_\d ∫_Ω ‖ ∇\d ‖_F^2 \ dA ⇒
+\min_\x ∫_Ω ‖ ∇(\x - \hat{\x}) ‖_F^2 \ dA ⇒
+\min_\x ∫_Ω ‖ \underbrace{∇\x}_\text{after} -
+\underbrace{∇\hat{\x}}_\text{before} ‖_F^2 \ dA.
+\\]
+If we think of the gradient of the position function $∇\x$ (with respect to the
+underlying parameterization $u,v$) as a local geometric [feature
+descriptor](https://en.wikipedia.org/wiki/Feature_(computer_vision)) then this
+energy can be re-understood as measuring the difference in this feature before
+and after the deformation. This is very sensible as we are trying to measure
+distortion. We would expect that a low-distortion deformation would correspond
+with a small change to local features.
+
+Unfortunately, the gradient of the position function $\x$ is a _poor_,
+first-order local feature.
+
+#### Laplacian-based energy
+
+If we model distortion as the change in a local feature descriptor, then a
+natural local and _relative_ descriptor would be one that compared the position
+of some point on the shape to the average of its local neighborhood. We have
+studied an operator that computes this in the smoothing assignment. The
+[Laplace(-Beltrami) operator](https://en.wikipedia.org/wiki/Laplace_operator)
+can be derived as taking exactly the difference of a functions value at a point
+and the average (i.e., [centroid](https://en.wikipedia.org/wiki/Centroid)) of
+an infinitesimal region around that point:
+
+\\[
+∆ f(\x) = \lim_{|B(\x)| → 0} \frac{1}{|B(\x))|} ∫_{B(\x)} f(\z) \;d\z - f(\x)
+\\]
+
+(see, e.g., "Differential coordinates for local mesh morphing and deformation"
+[Alexa et al. 2003] and expanded upon in "Laplacian Surface Editing"
+[Sorkine et al. 2004]).
+
+When applied to the embedding function $\x$ the Laplace operator computes the
+difference in _position_ between a point and its local neighborhood. This
+_vector_ points in the
+[normal](https://en.wikipedia.org/wiki/Normal_(geometry)) direction and its
+magnitude corresponds inversely with [how flat the surface is
+locally](https://en.wikipedia.org/wiki/Mean_curvature).
+
+Let's replace the gradient feature above with this _second-order_ feature
+descriptor and massage our optimization problem back in terms of
+displacements:
+
+\\[
+\min_\x ∫_Ω ‖ \underbrace{∆\x}_\text{after} -
+\underbrace{∆\hat{\x}}_\text{before} ‖^2 \ dA ⇒
+\min_\x ∫_Ω ‖ ∆(\x - \hat{\x}) ‖^2 \ dA ⇒
+\min_\d ∫_Ω ‖ ∆\d ‖^2 \ dA.
+\\]
+
+Just as we can show that harmonic functions ($∆\d = 0$) minimize the Dirichlet
+energy, we can use [calculus of
+variations](https://en.wikipedia.org/wiki/Calculus_of_variations) apply
+[Green's identity](https://en.wikipedia.org/wiki/Green's_identities) _twice_ to
+show that minimizers of the squared-Laplacian energy are _bi-_harmonic
+functions ($∆∆\d = 0$ or $∆²\d = 0$). Obviously all harmonic functions are also
+biharmonic functions. This implies that the space of biharmonic functions is
+strictly _larger_. In particular, this will allow use to ensure continuity of
+first derivatives across constrained values at handles.
+
+
+
+The fact that each coordinate of the displacement field $\d$ will be a
+bi-harmonic function is not so elucidating, but by minimizing the
+squared-Laplacian energy integrated over the domain, we can say that it is
+as-_harmonic_-as-possible. Harmonic functions include constant functions (i.e.,
+translations) but also any displacement that _bends_ in one direction by an
+equal and opposite amount as it bends in the other direction: $∆\d = 0 → ∂\d/∂u
++ ∂\d/∂v = 0 → ∂\d/∂u = -∂\d/∂v$.
+
+To discretize this energy we can make use of our discrete Laplacian operator
+$\L$. This matrix computes the locally _integrated_ Laplacian of a given
+function specified by per-vertex values $\f$. Now we would like to integrate
+the square of the _point-wise_ Laplacian. We can approximate the point-wise
+Laplacian by the local [integral
+average](https://en.wikipedia.org/wiki/Mean_of_a_function) of the Laplacian
+$\M^{-1}\L \f$. Integrating this over the mesh we have our complete approximate
+of the energy:
+
+\\[
+\begin{align}
+∫_Ω ‖∆\d‖² \;dA &≈ \tr{ \D^\transpose \L^\transpose \M^{-\transpose} \M \M^{-1} \L
+\D}\\
+&= \tr{ \D^\transpose \underbrace{\L^\transpose \M^{-1} \L}_{\Q} \D },
+\end{align}
+\\]
+where $\M ∈ \R^{n × n}$ is the mass-matrix for the given mesh and $\Q ∈
+\R^{n×n}$ can be thought of as the bi-Laplacian matrix.
+
+> #### $k$-harmonic
+>
+> The logical continuation of harmonic and biharmonic deformation is to consider
+> _triharmonic_ ($∆³\d = 0$) and _tetraharmonic_ ($∆⁴\d = 0$) and so on. It's
+> straightforward to implement these, though there are diminishing returns and
+> increasing costs.
+
+##### Precomputation
+
+With out loss of generality, assume that the rows of the unknown displacements
+$\D$ have been sorted so that displacements corresponding to handle vertices
+are in the bottom part:
+
+\\[
+\D = \left(\begin{array}{c}
+\D_\text{u} \\
+\D_\text{h}
+\end{array}
+\right)
+\\]
+
+Since the displacements at handles are _known_ before the optimization, we can
+separate the knowns and unknowns in the energy:
+
+\\[
+\min_{\D_\text{u}}
+\tr{(\D_\text{u}^\transpose \ \D_\text{h}^\transpose)
+\left(\begin{array}{cc}
+\Q_\text{u,u} & \Q_\text{u,h} \\
+\Q_\text{h,u} & \Q_\text{h,h}
+\end{array}\right)
+\left(\begin{array}{c}
+ \D_\text{u} \\
+ \D_\text{h}
+\end{array}
+\right)} \\
+\min_{\D_\text{u}}
+\tr{\D_\text{u}^\transpose \Q_\text{u,u} \D_\text{u} +
+2 \D_\text{u}^\transpose \Q_\text{u,h} \D_\text{h} +
+\underbrace{\D_\text{h}^\transpose \Q_\text{h,h}
+\D_\text{h}}_\text{constant}} \\
+\min_{\D_\text{u}}
+\tr{
+\D_\text{u}^\transpose \Q_\text{u,u} \D_\text{u} +
+2 \D_\text{u}^\transpose \Q_\text{u,h} \D_\text{h}}
+\\]
+where $\Q_\text{u,u} ∈ \R^{(n-k) × (n-k)}$ is the quadratic coefficients matrix
+corresponding to the unknown displacements.
+
+This quadratic optimization problem may solved by setting all partial
+derivatives with respect to degrees of freedom in $\D_\text{u}$ to zero:
+
+\\[
+2 \Q_\text{u,u} \D_\text{u} + 2 \Q_\text{u,h} \D_\text{h}
+= 0 → \D_\text{u} = \Q_\text{u,u}^{-1} \Q_\text{u,h} \D_\text{h}
+\\]
+
+If we don't change _which_ vertices are handles, but only change the positions
+of the selected handles, then only $\D_\text{h}$ changes above. In particular,
+the matrix $\Q_\text{u,u}$ is unchanged. Therefore, we can
+[prefactorize](https://en.wikipedia.org/wiki/Cholesky_decomposition) it so that
+_applying its inverse_ is fast (`igl::min_quad_with_fixed` does this for you).
+
+> Actually the entire term $\Q_\text{u,u}^{-1} \Q_\text{u,h} =: \W ∈ \R^{(n-k)
+> × k}$ does not change. The columns of $\W$ reveal how unknown displacements
+> respond to each handle point's displacement (see also "An intuitive framework
+> for real-time freeform modeling" [Botsch & Kobbelt 2004]). This provides a
+> gateway to the relationship with linear blend skinning and automatic
+> (biharmonic) weighting functions (see "Bounded Biharmonic Weights for
+> Real-Time Deformation" [Jacobson et al. 2011]).
+
+##### Trouble in paradise
+
+Biharmonic displacements work well for small deformations and deformations that
+do not imply a large rotation. However, the Laplacian of the position function
+$∆\x$ as a feature descriptor is _not_ rotation invariant. This problem is true
+for all linear differential features including the gradient of the embedding
+function $∇\x$ considered above (see "On linear variational surface deformation
+methods" [Botsch & Sorkine 2008]).
+
+
+
+This means that if the user transforms all of the handle locations by a rigid
+transformation $\T$ these energies will _not_ measure zero for a displacement
+equivalent to applying the rigid transformation $\T$ to the entire shape. We
+would like _global_ rotation invariance, but we would also like this property
+to apply _locally_ so that parts of the shape can easily rotate.
+
+### As-rigid-as-possible
+
+In the scenario where each handle
+$i$ are perfectly transformed by a [rigid
+transformation](https://en.wikipedia.org/wiki/Rigid_transformation) $\x_i =
+\Rot \hat{\x}_i + \t$, where $\Rot ∈ SO(3) ⊂ \R^{3×3}$ is a rotation matrix and
+$\t∈\R^3$ is a translation vector. If an
+[oracle](https://en.wikipedia.org/wiki/Oracle) could only tell us this
+particular rigid transformation then we could repair
+the gradient-based energy above by pre-rotating the rest shape by this
+transformation:
+
+\\[
+\begin{align}
+∫_Ω ‖ ∇ \x - ∇(\Rot \hat{\x} + \t) ‖² \;dA
+ &= ∫_Ω ‖ ∇ \x - ∇(\Rot \hat{\x}) - ∇\t‖² \;dA \\
+ &= ∫_Ω ‖ ∇ \x - \Rot ∇\hat{\x} ‖² \;dA,
+\end{align}
+\\]
+
+where the translation vector $\t$ falls out because a translation has constant
+gradient.
+
+We do not know the rotation $\Rot$ ahead of time, but we could be as generous
+as possible and use the "best" rotation $\Rot ← \argmin_\Rot ∫_Ω ‖ ∇\x - \Rot
+∇\hat{\x} ‖² \;dA$:
+
+\\[
+∫_Ω \left\|∇\x - \left( \argmin_\Rot ∫_Ω ‖ ∇\x - \Rot ∇\hat{\x} ‖² \;dA \right)∇\hat{\x}\right\|² \;dA.
+\\]
+
+If we treat $\Rot$ as a degree of freedom along with the unknown positions
+$\x$, we can unify this into an optimization over $\x$ and $\Rot$:
+
+\\[
+\min_{\x,\Rot∈SO(3)} ∫_Ω \left\|∇\x - \Rot ∇\hat{\x}\right\|² \;dA.
+\\]
+
+Optimizing this energy will ensure _global_ rotation invariance. To ensure
+_local_ rotation invariance, we can replace $\Rot ∈ SO(3)$ with a spatially
+varying _function_ $\Rot : Ω → SO(3)$ that outputs a "best" rotation for any
+point on the shape (see "A simple geometric model for elastic deformations"
+[Chao et al. 2010]). In this way, the optimal rotation will be locally rigid
+everywhere, or _as-rigid-as-possible_ (ARAP).
+
+
+
+> For embedded solid shapes, we can take the rest shape given by $\hat{\x}$ as
+> the parameterization $Ω$, so that $∇\hat{\x} = \I$. This allows us to rewrite
+> the as-rigid-as-possible energy as the square of the difference between the
+> [deformation gradient](#deformationgradient) and the closest rotation:
+>
+> \\[
+> ∫_Ω \left\|∇\x - \Rot ∇\hat{\x}\right\|² \;dA \\
+> ∫_Ω \left\|(∇\x + \I - \I) - \Rot \I \right\|² \;dA \\
+> ∫_Ω \left\|(\I + ∇\x - ∇\hat{x}) - \Rot \right\|² \;dA \\
+> ∫_Ω \left\|(\I + ∇\d) - \Rot \right\|² \;dA \\
+> ∫_Ω \left\|\F - \Rot \right\|² \;dA \\
+> \\]
+>
+> This form provides a bridge between the as-rigid-as-possible energy common in
+> geometry processing to _corotated linear elasticity_ used in
+> mechanics/physically-based simulation (made explicit in "A simple geometric
+> model for elastic deformations" [Chao et al. 2010]). See Section 3.4 of "FEM
+> Simulation of 3D Deformable Solids" [Sifakis 2012] for a graphics-mechanics
+> perspective.
+>
+
+#### Discrete as-rigid-as-possible energy
+
+For a triangle mesh with displacing vertices, the gradient of the embedding
+function is constant inside each triangle. In this way we can write the raw
+gradient energy above as a double sum over all half-edges $ij$ of all faces $f$
+in the mesh:
+
+\\[
+½ ∫_Ω ‖ ∇ \x - ∇\hat{\x}‖² \;dA = ½ ∑\limits_{f ∈ F} ∑\limits_{ ij ∈ f} c_{ij} ‖
+(\v_i-\v_j) - (\hat{\v}_i-\hat{\v}_j)‖²,
+\\]
+where $c_{ij}$ is cotangent of the angle opposite half-edge $ij$.
+
+To inject localized best fit rotations, we will assign an unknown rotation
+matrix $\Rot_k$ to each vertex $k$ of the mesh and accounts for a third of the
+energy integrated over incident triangles:
+\\[
+½ ∫_Ω ‖ ∇ \x - \Rot ∇\hat{\x}‖² \;dA =
+⅙ ∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)}
+c_{ij} ‖ (\v_i-\v_j) - \Rot_k (\hat{\v}_i-\hat{\v}_j)‖²,
+\\]
+where $F(k)$ is the set of all faces incident on the $k$-th vertex.
+
+> ##### Where do rotations _live_?
+> We have assigned a rotation for each vertex: there are $n$ rotation matrices
+> as auxiliary degrees of freedom. This is in contrast to assigning
+> rotations per face (as in "A Local/Global Approach to Mesh Parameterization"
+> [Liu et al. 2008]). Per-face--or more generally per-element--rotations work
+> well for [codimension](https://en.wikipedia.org/wiki/Codimension) zero
+> objects (triangle meshes in $\R²$ or tetrahedral meshes in $\R³$). But for
+> triangle mesh surfaces in $\R³$ (i.e., codimension 1), per-face rotations
+> would lead to "crumpling") because _bending_ along edges would not be
+> measured.
+
+#### Optimization
+
+The simplest method for optimizing the ARAP energy is by alternating between
+
+ 1. finding the optimal rotations $\Rot_k$ assuming the vertex positions $\V$ are
+ fixed, and
+ 2. finding the optimal vertex positions $\V$ assuming all rotations $\Rot_k$ are
+ fixed.
+
+Each rotation $\Rot_k$ only affects the local energy and doesn't interact with
+the _other_ rotations. So each can be optimized _locally_. In contrast, the
+mesh vertex positions $\V$ depend on each other requiring a _global_ solve. In
+the geometry processing literature, this is known as a local-global
+optimization (see "As-rigid-as-possible surface modeling" [Sorkine & Alexa
+2007]). It is also known as "alternating block coordinate descent" because we
+have separated the variables into disjoint sets ($\V,\Rot_1,…,\Rot_n$) and taking
+the optimal descent direction for each independently.
+
+Observing the discrete energy above we can see that the energy is quadratic in
+$\V$ and quadratic in each $\Rot_k$. Let's start by separating the terms that are
+quadratic and linear in $\V$:
+
+\\[
+⅙ \underbrace{∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)} c_{ij} (\v_i-\v_j)^\transpose(\v_i-\v_j)}_\text{quadratic}
++
+⅙ \underbrace{∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)} c_{ij} (\v_i-\v_j)^\transpose \Rot_k (\hat{\v}_i-\hat{\v}_j)}_\text{linear}
+\\]
+
+if we stack the rotation matrices $\Rot_k$ into large matrix $\Rot ∈ \R^{3n ×
+3}$ then we can write this energy in matrix form as:
+
+\\[
+\tr{ \V^\transpose \L \V } + \tr{ \V^\transpose \K \Rot },
+\\]
+where $\L ∈ \R^{n × n}$ is the familiar cotangent discrete Laplacian matrix and
+$\K ∈ \R^{n × 3n}$ sparse matrix containing cotangents multiplied against
+differences across edges in the rest mesh (e.g., $\hat{\v}_i - \hat{\v}_j$).
+
+> ###### I'm so confused. What's in the $\K$ matrix?
+>
+> Let's take it slow. The $\K$ matrix is represents the
+> [bilinear form](https://en.wikipedia.org/wiki/Bilinear_form) that combines unknown
+> vertex positions and unknown rotations. We have identified above that we can
+> write this in summation or matrix form:
+> \\[
+> ⅙ ∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)} c_{ij} (\v_i-\v_j)^\transpose \Rot_k (\hat{\v}_i-\hat{\v}_j)
+> = \tr{ \V^\transpose \K \Rot },
+> \\]
+> but how did we get here?
+>
+> Let's start with the summation form. The _constants_ of this formula are the
+> $c_{ij}$ terms and the $(\hat{\v}_i-\hat{\v}_j)$ terms. Since these always
+> appear together, let us merge them into weighted edge difference vectors
+> $c_{ij} (\hat{\v}_i-\hat{\v}_j) =: \hat{\e}_{ij} ∈ \R³$:
+>
+> \\[
+> ⅙ ∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)} \underbrace{(\v_i-\v_j)^\transpose
+> \Rot_k \hat{\e}_{ij}}_{∈\R},
+> \\]
+> the inner term in the summation is an [inner
+> product](https://en.wikipedia.org/wiki/Inner_product_space); that is, a
+> [scalar](https://en.wikipedia.org/wiki/Scalar_(mathematics)). Let's expose this
+> by expanding the matrix-vector products of the inner-product:
+> \\[
+> ⅙ ∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)} ∑\limits_{α=1}^3 ∑\limits_{β=1}^3
+> (v_i^α - v_j^α)R_k^{αβ}\hat{e}_{ij}^β.
+> \\]
+>
+> If our mesh is stored as a vertex list and face list, it's not easy/efficient
+> to loop over per-vertex rotations (outer sum) and then over all half-edges of
+> incident faces (second sum). Instead, let's rearrange these sums to loop over
+> all faces first, then the half-edges of that face, and then over all
+> per-vertex rotations that involve this half-edge:
+> \\[
+> ⅙ ∑\limits_{f=1}^m ∑\limits_{ ij ∈ E(f)} \quad ∑\limits_{k | ij ∈F(k)} \quad ∑\limits_{α=1}^3 ∑\limits_{β=1}^3
+> (v_i^α - v_j^α)R_k^{αβ}\hat{e}_{ij}^β,
+> \\]
+> where the third sum is over all rotations $k$ such that the half-edge $ij$
+> belongs to the half-edges of the faces incident on the $k$-th vertex: $k | ij
+> ∈F(k)$. Well, this means $k$ can either be $i$ or $j$ or the third vertex of
+> the $f$-th face.
+>
+> Now let's turn our attention back to the
+> [summand](https://en.wiktionary.org/wiki/summand). The terms indexed by $α$
+> never _mix_. That is, we never add/multiply $v_i^α$, $v_j^γ$, and $R_k^{δβ}$
+> unless $α=γ=δ$. This implies that we can write this summation in matrix form
+> as:
+> \\[
+> \V₁^\transpose \K₁ \Rot₁ +
+> \V₂^\transpose \K₂ \Rot₂ +
+> \V₃^\transpose \K₃ \Rot₃,
+> \\]
+> where $\V_α ∈ \R^n$ is $α$-th column of $\V$, $\Rot_α ∈ \R^{3n}$ is the $α$-th column of
+> $\Rot$ and $\K₁,\K₂,\K₃ ∈ \R^{n × 3n}$ are sparse matrices.
+>
+> _Further_, the constant term $\hat{e}^β_{ij}$ in the summation acts the same on
+> $(v_i^α - v_j^α)R_k^{αβ}$ for any $α$ value. This implies that $\K₁=\K₂=\K₃$,
+> so we can reduce the matrix form to:
+> \\[
+> \tr{\V \K \Rot}.
+> \\]
+>
+> Finally, we can answer what is in each entry $K_{vw}$, where $v$ chooses the
+> row and $w=3k+β$ chooses the column of $\K$. Treating our nested summation as
+> nested for-loops, we can increment entries in $\K$
+> \\[
+> \begin{align}
+> K_{i\ 3k+β} & \mathrel{+}= \hat{e}_{ij}^β, \\
+> K_{j\ 3k+β} & \mathrel{+}= -\hat{e}_{ij}^β, \\
+> \end{align}
+> \\]
+> for each half-edge encountered.
-- Not smooth at constraints
-- Each coordinate of "displacement field" is a harmonic function (picture of
- arrows and picture of x-coordinate)
- - Each coordinate is a heat diffusion toward some average
- - As-_constant_-as-possible
-#### Laplacian-based energy
-- Each coordinate of "displacement field" is a biharmonic function
- - larger space of functions than harmonic
- - includes functions that are smooth at the constraints
- - As-_harmonic_-as-possible (not simply harmonic since we have to satisfy more
- constraints): locally harmonic roughly corresponds to locally
- "fair"/non-ripply (technically, harmonic means bending in one direction
- equals opposite amount of bending in the other direction. [Saddle]()
- shaped)
+##### Local step
-- Corresponds to higher-order interpolation à la Hermite splines
+Minimizing this energy with respect $\Rot$ corresponds to minimizing:
-Dependent on local coordinate system, but still confused by large rotations.
+\\[
+\tr{ \underbrace{\V^\transpose \K}_{\C^\transpose} \Rot },
+\\]
+where $\C ∈ \R^{3n × 3}$ stacks weighted covariance matrices $\C_k ∈ \R^{3 ×
+3}$ for each region _covered_ by the corresponding rotation $\Rot_k$. We have
+seen this problem before in the registration assignment. For each $\C_k$,
+$\Rot_k$ will be the closest rotation matrix solved via [singular value
+decomposition](https://en.wikipedia.org/wiki/Singular_value_decomposition).
-#### $k$-harmonic
+##### Global step
-The logical continuation of harmonic and biharmonic deformation is to consider
-triharmonic and tetraharmonic and so on. It's straightforward to extend these,
-though there are diminishing returns and increasing costs.
+Minimizing the energy above with respect to $\V$ corresponds to solving a
+Dirichlet energy-like minimization problem:
-### Non-Linear Methods
+\\[
+\min_\V \tr{ \V^\transpose \L \V } + \tr{ \V^\transpose \B }
+\\]
-Take a step back. Why didn't we like the gradient based energy? Global
-coordinate frame, non-smoothness at constraints. Put aside non-smoothness at
-constraints. Let's try to tackle large rotations.
+where $\K \Rot =: \B ∈ \R^{n × 3}$ is a matrix of rotated vertex gradients.
+Adding the handle constraints to the corresponding rows of $\V$ this is easily
+minimized by setting all partial derivatives with respect to the unknowns in
+$\V$ equal to zero (as in the linear methods above).
-**Deformation gradient** for 2D shapes
+##### Implementation
-#### As-rigid-as-possible energy
+In order to facilitate interactive deformation we would like our local and
+global iterations to be computed as quickly as possible. Since the quadratic
+form in the global step is the _same_ regardless of the current rotations or
+current handle positions, we can
+[prefactorize](https://en.wikipedia.org/wiki/Cholesky_decomposition) it (again,
+as above). The matrix $\K$ also does not depend on the rotations, current
+positions or handle positions, so we can pre-build this matrix. This way
+computing $\C$ and $\B$ is a simple matrix-matrix multiplication.
-Warm-up: For 2D shapes
+> **Note:** When constructing $\K$ it's easiest to iterate over _all_
+> half-edges in the mesh (by iterating over all faces and then each of the
+> three edges). Each half-edge $ij$ _contributes_ terms tying $\v_i,\v_j$ to
+> _each_ of the (three) rotations $\Rot_k$ that apply against their difference
+> (see "Fast Automatic Skinning Transformations" [Jacobson et al. 2012]).
-> ##### Volumes
+##### Still trouble in paradise
-##### Surfaces
+
-Overlapping integration regions
+The as-rigid-as-possible deformation method for surfaces described above has a
+number of remaining problems:
+
+ 1. like its gradient-based energy ancestor the deformation is not smooth at
+ constraints;
+ 2. the energy punishes _bending_ of the surface--which is good--but does so in
+ a way that diminishes as the mesh becomes higher and higher resolution, in
+ otherwords, the discrete energy is mesh-resolution dependent; and
+ 3. the energy is biased by the original combinatorics of the mesh (even in
+ flat regions).
## Tasks
### Blacklist
- - `igl::harmonic`
- `igl::arap`
+ - `igl::arap_linear_block`
- `igl::covariance_scatter_matrix`
+ - `igl::harmonic`
### Whitelist
- - `igl::fit_rigid` (or your previous assignment's `closest_rotation`)
+ - `igl::polar_svd3x3` (or your previous assignment's `closest_rotation`)
+ - `igl::min_quad_with_fixed`
+ - `igl::cotmatrix_entries`
+ - `igl::cotmatrix` (or your previous implementation)
+ - `igl::massmatrix` (or your previous implementation)
+
+### `src/biharmonic_precompute.cpp`
+
+Precompute data needed to efficiently solve for a biharmonic deformation given
+a mesh with vertices `V` and faces `F` and a list of selected vertices as
+indices `b` into `V`. The output should be a prefacorized system using the
+`data` struct employed by `igl::min_quad_with_fixed`.
+
+### `src/biharmonic_solve.cpp`
-### k-harmonic precomputation
+Given precomputation `data` and a list of handle _displacements_ determine
+_displacements_ for all vertices in the mesh.
-### k-harmonic deformation
+### `src/arap_precompute.cpp`
-### covariance scatter matrix
+Precompute data needed to efficiently conduct local-global iterations for an
+arap deformation. This includes the `data` struct employed by
+`igl::min_quad_with_fixed` to solve the global step" and constructing the
+bilinear form `K` that mixes rotation degrees of freedom with unknown
+positions for preparing the covariance matrices of the local step and the
+linear term of the global step.
-### arap precompuation
+### `src/arap_single_iteration.cpp`
-### arap deformation
+Given precomputed data (`data` and `K`), handle _positions_ `bc` and current
+positions of all vertices `U`, conduct a _single_ iteration of the local-global
+solver for minimizing the as-rigid-as-possible energy. Output the _positions_
+of all vertices of the mesh (by overwriting `U`).
diff --git a/images/bump-k-harmonic.jpg b/images/bump-k-harmonic.jpg
new file mode 100644
index 0000000..83b1246
Binary files /dev/null and b/images/bump-k-harmonic.jpg differ
diff --git a/images/bunny-biharmonic.gif b/images/bunny-biharmonic.gif
new file mode 100644
index 0000000..14a8a4e
Binary files /dev/null and b/images/bunny-biharmonic.gif differ
diff --git a/images/bunny-harmonic.gif b/images/bunny-harmonic.gif
new file mode 100644
index 0000000..bed1fa9
Binary files /dev/null and b/images/bunny-harmonic.gif differ
diff --git a/images/knight-arap-high-vs-low-resolution.gif b/images/knight-arap-high-vs-low-resolution.gif
new file mode 100644
index 0000000..302b3f6
Binary files /dev/null and b/images/knight-arap-high-vs-low-resolution.gif differ
diff --git a/images/knight-arap-large-rotation.gif b/images/knight-arap-large-rotation.gif
new file mode 100644
index 0000000..3b21405
Binary files /dev/null and b/images/knight-arap-large-rotation.gif differ
diff --git a/images/knight-biharmonic-large-rotation.gif b/images/knight-biharmonic-large-rotation.gif
new file mode 100644
index 0000000..e0e10db
Binary files /dev/null and b/images/knight-biharmonic-large-rotation.gif differ
diff --git a/images/knight-deformation.gif b/images/knight-deformation.gif
new file mode 100644
index 0000000..a8d0f81
Binary files /dev/null and b/images/knight-deformation.gif differ
diff --git a/include/arap_precompute.h b/include/arap_precompute.h
new file mode 100644
index 0000000..77a81e9
--- /dev/null
+++ b/include/arap_precompute.h
@@ -0,0 +1,33 @@
+#ifndef ARAP_PRECOMPUTE_H
+#define ARAP_PRECOMPUTE_H
+#include
+#include
+
+namespace igl
+{
+ template struct min_quad_with_fixed_data;
+}
+
+// Precompute data needed to efficiently conduct local-global iterations for an
+// arap deformation. This includes the `data` struct employed by
+// `igl::min_quad_with_fixed` to solve the global step" and constructing the
+// bi-linear form `K` that mixes rotation degrees of freedom with unknown
+// positions for preparing the covariance matrices of the local step and the
+// linear term of the global step.
+//
+// Inputs:
+// V #V by dim vertex positions
+// F #F by simplex-size list of element indices
+// b #b indices into V of handle vertices
+// Outputs:
+// data pre-factorized system matrix etc. (see `igl::min_quad_with_fixed`)
+// K #R*dim by #V
+void arap_precompute(
+ const Eigen::MatrixXd & V,
+ const Eigen::MatrixXi & F,
+ const Eigen::VectorXi & b,
+ igl::min_quad_with_fixed_data & data,
+ Eigen::SparseMatrix & K);
+
+#endif
+
diff --git a/include/arap_single_iteration.h b/include/arap_single_iteration.h
new file mode 100644
index 0000000..687eb5e
--- /dev/null
+++ b/include/arap_single_iteration.h
@@ -0,0 +1,29 @@
+#ifndef ARAP_SINGLE_ITERATION_H
+#define ARAP_SINGLE_ITERATION_H
+#include
+#include
+
+namespace igl
+{
+ template struct min_quad_with_fixed_data;
+}
+
+// Given precomputed data (`data` and `K`), handle _positions_ `bc` and current
+// positions of all vertices `U`, conduct a _single_ iteration of the
+// local-global solver for minimizing the as-rigid-as-possible energy. Output
+// the _positions_ of all vertices of the mesh (by overwriting `U`).
+//
+// Inputs:
+// data pre-factorized system matrix etc. (see `arap_precompute`
+// K pre-constructed bi-linear term of energy combining rotations and
+// positions
+// U #V by dim list of current positions (see output)
+// Outputs:
+// U #V by dim list of new positions (see input)
+void arap_single_iteration(
+ const igl::min_quad_with_fixed_data & data,
+ const Eigen::SparseMatrix & K,
+ const Eigen::MatrixXd & bc,
+ Eigen::MatrixXd & U);
+
+#endif
diff --git a/include/biharmonic_precompute.h b/include/biharmonic_precompute.h
new file mode 100644
index 0000000..c466e34
--- /dev/null
+++ b/include/biharmonic_precompute.h
@@ -0,0 +1,29 @@
+#ifndef BIHARMONIC_PRECOMPUTE_H
+#define BIHARMONIC_PRECOMPUTE_H
+#include
+#include
+
+
+namespace igl
+{
+ template struct min_quad_with_fixed_data;
+}
+
+// Precompute data needed to efficiently solve for a biharmonic deformation
+// given a mesh with vertices `V` and faces `F` and a list of selected vertices
+// as indices `b` into `V`. The output should be a prefacorized system using
+// the `data` struct employed by `igl::min_quad_with_fixed`.
+//
+// Inputs:
+// V #V by dim vertex positions
+// F #F by simplex-size list of element indices
+// b #b indices into V of handle vertices
+// Outputs:
+// data pre-factorized system matrix etc. (see `igl::min_quad_with_fixed`)
+void biharmonic_precompute(
+ const Eigen::MatrixXd & V,
+ const Eigen::MatrixXi & F,
+ const Eigen::VectorXi & b,
+ igl::min_quad_with_fixed_data & data);
+
+#endif
diff --git a/include/biharmonic_solve.h b/include/biharmonic_solve.h
new file mode 100644
index 0000000..e2018ed
--- /dev/null
+++ b/include/biharmonic_solve.h
@@ -0,0 +1,25 @@
+#ifndef BIHARMONIC_SOLVE_H
+#define BIHARMONIC_SOLVE_H
+#include
+
+namespace igl
+{
+ template struct min_quad_with_fixed_data;
+}
+
+// Given precomputation data and a list of handle _displacements_ determine
+// _displacements_ for all vertices in the mesh.
+//
+// Inputs:
+// data pre-factorized system matrix etc. (see `biharmonic_precompute`
+// bc #b by 3 list of displacements for each handle vertex
+// Outputs:
+// D #V by 3 list of displacements
+void biharmonic_solve(
+ const igl::min_quad_with_fixed_data & data,
+ const Eigen::MatrixXd & bc,
+ Eigen::MatrixXd & D);
+
+
+#endif
+
diff --git a/main.cpp b/main.cpp
index 5d6487a..737b6b0 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,20 +1,23 @@
+#include "biharmonic_precompute.h"
+#include "biharmonic_solve.h"
+#include "arap_precompute.h"
+#include "arap_single_iteration.h"
+#include
#include
#include
#include
#include
#include
-#include
-#include
#include
#include
-#include
#include
#include
+// Undoable
struct State
{
- Eigen::MatrixXd CV;
- Eigen::MatrixXd CU;
+ // Rest and transformed control points
+ Eigen::MatrixXd CV, CU;
bool placing_handles = true;
} s;
@@ -50,85 +53,113 @@ int main(int argc, char *argv[])
Eigen::MatrixXd V,U;
Eigen::MatrixXi F;
long sel = -1;
- Eigen::RowVector3f last;
+ Eigen::RowVector3f last_mouse;
+ igl::min_quad_with_fixed_data biharmonic_data, arap_data;
+ Eigen::SparseMatrix arap_K;
// Load input meshes
igl::read_triangle_mesh(
- (argc>1?argv[1]:"../shared/data/bunny.off"),V,F);
+ (argc>1?argv[1]:"../shared/data/decimated-knight.off"),V,F);
U = V;
- // Load data into MatrixXd rather than VectorXd for simpler `smooth` API
- // Just use y-coordinates as data to be smoothed
- // Create a libigl Viewer object to toggle between point cloud and mesh
igl::viewer::Viewer viewer;
std::cout<bool
{
- last = Eigen::RowVector3f(
- viewer.current_mouse_x,
- viewer.core.viewport(3) - viewer.current_mouse_y,
- 0);
+ last_mouse = Eigen::RowVector3f(
+ viewer.current_mouse_x,viewer.core.viewport(3)-viewer.current_mouse_y,0);
if(s.placing_handles)
{
+ // Find closest point on mesh to mouse position
int fid;
- Eigen::Vector3f bc;
- // Cast a ray in the view direction starting from the mouse position
+ Eigen::Vector3f bary;
if(igl::unproject_onto_mesh(
- last.head(2),
+ last_mouse.head(2),
viewer.core.view * viewer.core.model,
viewer.core.proj,
viewer.core.viewport,
V, F,
- fid, bc))
+ fid, bary))
{
- push_undo();
- s.CV.conservativeResize(s.CV.rows()+1,3);
- s.CV.row(s.CV.rows()-1) =
- V.row(F(fid,0))*bc(0)+ V.row(F(fid,1))*bc(1)+ V.row(F(fid,2))*bc(2);
- update();
- return true;
+ long c;
+ bary.maxCoeff(&c);
+ Eigen::RowVector3d new_c = V.row(F(fid,c));
+ if(s.CV.size()==0 || (s.CV.rowwise()-new_c).rowwise().norm().minCoeff() > 0)
+ {
+ push_undo();
+ s.CV.conservativeResize(s.CV.rows()+1,3);
+ // Snap to closest vertex on hit face
+ s.CV.row(s.CV.rows()-1) = new_c;
+ update();
+ return true;
+ }
}
}else
{
- // Get closest control point
+ // Move closest control point
Eigen::MatrixXf CP;
igl::project(
Eigen::MatrixXf(s.CU.cast()),
viewer.core.view * viewer.core.model,
viewer.core.proj, viewer.core.viewport, CP);
- Eigen::VectorXf D = (CP.rowwise()-last).rowwise().norm();
+ Eigen::VectorXf D = (CP.rowwise()-last_mouse).rowwise().norm();
sel = (D.minCoeff(&sel) < 30)?sel:-1;
if(sel != -1)
{
- last(2) = CP(sel,2);
+ last_mouse(2) = CP(sel,2);
push_undo();
update();
return true;
@@ -141,16 +172,25 @@ int main(int argc, char *argv[])
{
if(sel!=-1)
{
- Eigen::RowVector3f drag(
+ Eigen::RowVector3f drag_mouse(
viewer.current_mouse_x,
viewer.core.viewport(3) - viewer.current_mouse_y,
- last(2));
+ last_mouse(2));
Eigen::RowVector3f drag_scene,last_scene;
- igl::unproject(drag, viewer.core.view * viewer.core.model, viewer.core.proj,viewer.core.viewport,drag_scene);
- igl::unproject(last, viewer.core.view * viewer.core.model, viewer.core.proj,viewer.core.viewport,last_scene);
+ igl::unproject(
+ drag_mouse,
+ viewer.core.view*viewer.core.model,
+ viewer.core.proj,
+ viewer.core.viewport,
+ drag_scene);
+ igl::unproject(
+ last_mouse,
+ viewer.core.view*viewer.core.model,
+ viewer.core.proj,
+ viewer.core.viewport,
+ last_scene);
s.CU.row(sel) += (drag_scene-last_scene).cast();
- last = drag;
-
+ last_mouse = drag_mouse;
update();
return true;
}
@@ -166,6 +206,12 @@ int main(int argc, char *argv[])
{
switch(key)
{
+ case 'M':
+ case 'm':
+ {
+ method = (Method)(((int)(method)+1)%((int)(NUM_METHODS)));
+ break;
+ }
case 'R':
case 'r':
{
@@ -173,12 +219,24 @@ int main(int argc, char *argv[])
s.CU = s.CV;
break;
}
+ case 'U':
+ case 'u':
+ {
+ // Just trigger an update
+ break;
+ }
case ' ':
push_undo();
s.placing_handles ^= 1;
- if(!s.placing_handles)
+ if(!s.placing_handles && s.CV.rows()>0)
{
+ // Switching to deformation mode
s.CU = s.CV;
+ Eigen::VectorXi b;
+ igl::snap_points(s.CV,V,b);
+ // PRECOMPUTATION FOR DEFORMATION
+ biharmonic_precompute(V,F,b,biharmonic_data);
+ arap_precompute(V,F,b,arap_data,arap_K);
}
break;
default:
@@ -188,32 +246,33 @@ int main(int argc, char *argv[])
return true;
};
+ // Special callback for handling undo
viewer.callback_key_down =
[&](igl::viewer::Viewer &, unsigned char key, int mod)->bool
{
- switch(key)
+ if(key == 'Z' && (mod & GLFW_MOD_SUPER))
{
- case 'Z':
- if(mod & GLFW_MOD_SUPER)
- {
- if(mod & GLFW_MOD_SHIFT)
- {
- redo();
- }else
- {
- undo();
- }
- }
- update();
- break;
+ (mod & GLFW_MOD_SHIFT) ? redo() : undo();
+ update();
+ return true;
+ }
+ return false;
+ };
+ viewer.callback_pre_draw =
+ [&](igl::viewer::Viewer &)->bool
+ {
+ if(viewer.core.is_animating && !s.placing_handles && method == ARAP)
+ {
+ arap_single_iteration(arap_data,arap_K,s.CU,U);
+ update();
}
return false;
};
viewer.data.set_mesh(V,F);
viewer.core.show_lines = false;
viewer.core.is_animating = true;
+ viewer.data.face_based = true;
update();
viewer.launch();
-
return EXIT_SUCCESS;
}
diff --git a/shared b/shared
index 8e5599c..9e70205 160000
--- a/shared
+++ b/shared
@@ -1 +1 @@
-Subproject commit 8e5599cfa7eb07492bb17682971373a85ff88c01
+Subproject commit 9e7020587d5ca3fc8a5e18866973e916dbadd5f0
diff --git a/src/arap_precompute.cpp b/src/arap_precompute.cpp
new file mode 100644
index 0000000..75b5b00
--- /dev/null
+++ b/src/arap_precompute.cpp
@@ -0,0 +1,41 @@
+#include "arap_precompute.h"
+#include
+#include
+typedef Eigen::Triplet T;
+
+void arap_precompute(
+ const Eigen::MatrixXd & V,
+ const Eigen::MatrixXi & F,
+ const Eigen::VectorXi & b,
+ igl::min_quad_with_fixed_data & data,
+ Eigen::SparseMatrix & K)
+{
+ Eigen::SparseMatrix L;
+ igl::cotmatrix(V, F, L);
+
+ igl::min_quad_with_fixed_precompute(L, b, Eigen::SparseMatrix(), false, data);
+
+ int nV = V.rows();
+ std::vector tripletList;
+
+ for(int f = 0; f < F.rows(); f++) {
+ for(int a = 0; a < 3; a++) {
+ int i = F(f, (a + 1) % 3);
+ int j = F(f, (a + 2) % 3);
+
+ Eigen::Vector3d diff = V.row(i) - V.row(j);
+ Eigen::Vector3d eij = L.coeff(i, j) * diff / 6.0;
+
+ for(int ki = 0; ki < 3; ki++) {
+ int k = F(f, (a + ki) % 3);
+ for(int B = 0; B < 3; B++) {
+ tripletList.push_back(T(i, 3 * k + B, eij[B])); // Duplicates will be summed
+ tripletList.push_back(T(j, 3 * k + B, -eij[B]));
+ }
+ }
+ }
+ }
+
+ K.resize(nV, 3 * nV);
+ K.setFromTriplets(tripletList.begin(), tripletList.end());
+}
diff --git a/src/arap_single_iteration.cpp b/src/arap_single_iteration.cpp
new file mode 100644
index 0000000..de406f6
--- /dev/null
+++ b/src/arap_single_iteration.cpp
@@ -0,0 +1,24 @@
+#include "arap_single_iteration.h"
+
+#include
+#include
+
+void arap_single_iteration(
+ const igl::min_quad_with_fixed_data & data,
+ const Eigen::SparseMatrix & K,
+ const Eigen::MatrixXd & bc,
+ Eigen::MatrixXd & U)
+{
+ Eigen::MatrixXd C = K.transpose() * U;
+ Eigen::MatrixXd R(C.rows(), C.cols());
+
+ for(int k = 0; k < data.n; k++) {
+ Eigen::Matrix3d C_k = C.block(k * 3, 0, 3, 3);
+ Eigen::Matrix3d R_k;
+ C_k.normalize();
+ igl::polar_svd3x3(C_k, R_k);
+ R.block(k * 3, 0, 3, 3) = R_k;
+ }
+
+ igl::min_quad_with_fixed_solve(data, K*R, bc, Eigen::VectorXd(), U);
+}
diff --git a/src/biharmonic_precompute.cpp b/src/biharmonic_precompute.cpp
new file mode 100644
index 0000000..c6ee405
--- /dev/null
+++ b/src/biharmonic_precompute.cpp
@@ -0,0 +1,25 @@
+#include "biharmonic_precompute.h"
+
+#include
+#include
+#include
+
+void biharmonic_precompute(
+ const Eigen::MatrixXd & V,
+ const Eigen::MatrixXi & F,
+ const Eigen::VectorXi & b,
+ igl::min_quad_with_fixed_data & data)
+{
+ Eigen::SparseMatrix L;
+ igl::cotmatrix(V, F, L);
+
+ Eigen::SparseMatrix M;
+ igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_DEFAULT, M);
+
+ Eigen::SimplicialLDLT> solver(M);
+ Eigen::SparseMatrix M_inv_L = solver.solve(L).sparseView();
+ Eigen::SparseMatrix Q = L.transpose() * M_inv_L; // L^T M^-1 L
+
+ igl::min_quad_with_fixed_precompute(Q, b, Eigen::SparseMatrix(), false, data);
+}
+
diff --git a/src/biharmonic_solve.cpp b/src/biharmonic_solve.cpp
new file mode 100644
index 0000000..1a0d2ec
--- /dev/null
+++ b/src/biharmonic_solve.cpp
@@ -0,0 +1,11 @@
+#include "biharmonic_solve.h"
+#include
+
+void biharmonic_solve(
+ const igl::min_quad_with_fixed_data & data,
+ const Eigen::MatrixXd & bc,
+ Eigen::MatrixXd & D)
+{
+ igl::min_quad_with_fixed_solve(data, Eigen::VectorXd::Zero(data.n), bc, Eigen::VectorXd(), D);
+}
+
+
+
+
+
+
+
+