Deformation - Laplacian editing, ARAP

Path of the code: [05_laplacian_deformation]
The objective of this scene is to model interactive laplacian deformation similarily to the following two articles

Current state of the code

In the current state of the program a planar grid surface is displayed. And an interactive vertex selection system is proposed (with no effect on the surface so far).
The blue vertices are supposed model static vertex position that should be preserved during deformation.
The red vertices are supposed to model target vertices that the user can moves.
To translate the target vertices, press shift key and drag and drop the mouse. A white sphere indicates where the target vertex has been translated.
Another mode "Select constraint" is proposed where these constraints can be modified.

Laplacian surface deformation

The deformation of the surface can be expressed as the set of position \(q_i\) minimizing the following energy
This energy can be reformulated in a matrix form with unknown \(q=(q_0,\cdots,q_{N-1})\)
Note that \(M\) is a large sparse matrix of size \( (N+N_c) \times N\), where \(N_c\) is the number of constraints.
> Your goal is to build the matrix and vector associated to this energy, and solve the linear system in order to compute the new surface.

Solving linear system using Eigen

Eigen is C++ library for efficient matrix manipulation and linear system solver. You can use it to encode your matrix and vector associated to the least square problem. This library (composed only of header files) is already included with the code library (in third_party/eigen).
Note that you can use dense matrix solver, but to fully take advantage of the sparse structure of the matrix, you can directly use the Sparse matrix and solver proposed by Eigen.
If you are not familiar with Eigen, the following code shows a basic example of sparse system solver for least square problem (note that several solvers are proposed in Eigen, and you may experiment with them) 
#include <iostream>
#include <Eigen/Sparse>

int main()
{
    // Build matrix
    Eigen::SparseMatrix<float> M;
    
    
    int N = 50;
    M.resize(N,N);
    // Note: A sparse matrix doesn't store explicitly 0-value entries

    // Reserve space to store non-zero coefficients
    //   ex. N rows, 10 non-zero coeffs per rows (only a guess for efficiency purpose)
    M.reserve(Eigen::VectorXi::Constant(N,10));

    // Fill the matrix ... M.coeffRef(i,j) = value
    for(int k=0; k<N; ++k) {
        M.coeffRef(k,k)=1.0f; M.coeffRef(k,(k+1)%N)=-0.3f;
    }
    
    // Compress the matrix representation
    M.makeCompressed();

    // Factorization for a solver (here Conj. Gradient)
    Eigen::LeastSquaresConjugateGradient< Eigen::SparseMatrix<float> > solver;
    solver.compute(M);
    solver.setTolerance(1e-6f);

    // Build RHS
    Eigen::VectorXf b;
    b = Eigen::VectorXf::Random(N);

    // Solve linear system
    Eigen::VectorXf x = solver.solve(b);
    // or solver.solveWithGuess(b, x0) 
    //   faster with good guess

    std::cout <<"Check solution : "
       << (M*x-b).norm() << std::endl;   
    std::cout <<"N iterations : "
       << solver.iterations() << std::endl;   

    return 0;
}

Solving for Laplacian deformation

Notes:

As-Rigid-As-Possible

The As-Rigid-As-Possible method adds to the previous energy formulation the use of a matrix \(R(q_i)\): rotation corresponding to the optimal rigid transform between \(p_i\) (and its neighborhood) and \(q_i\) (and its neighborhood).
This energy is solved using an iterative process interleaving between theses two steps
> Implement the ARAP deformation and oberve that the rotation of the surface automatically adapts to the constraints.