Linear geometric algebra rotor estimator for efficient mesh deformation

: The authors solve the problem of estimating the best rotation aligning two sets of corresponding vectors (also known as Wahba's problem or point cloud registration). The proposed method is among the fastest methods reported in recent literatures, moreover it is robust to noise, accurate and simpler than most other methods. It is based on solving the linear equations derived from the formulation of the problem in Euclidean Geometric Algebra. The authors show its efficiency in two applications: the as-rigid-as-possible (ARAP) surface modelling and the more smooth rotation enhanced ARAP mesh animation which is the only method capable of deforming surface modes with quality of tetrahedral models. Mesh deformation is a key technique in games, automated construction and robotics. The ARAP technique along with its improved variants, although have been extensively studied, can still not be achieved efficiently. Linear geometric algebra based rotor solution proposed in this study gives another perspective of the kernel problem. This, however, not only improves the real performance of the three-dimensional mesh deformation, but also provides a brand new computationally efficient solution to the Wahba's problem and point cloud registration, which has been closely related to the automation science and engineering.


Introduction
The concept of mesh deformation arises from the emerging needs of the highly dynamic three-dimensional (3D) computer animation, which has been extensively employed in computer games, visualisation, computer aided design and intelligent manufacturing [1]. Another similar terminology, i.e. the surface flattening has also been widely studied in automated construction [2]. In principle, the target of these approaches is using 2D meshes to construct specific 3D models. As complex 3D models own quite abundant surface details, the deformation process is time consuming. The as-rigidas-possible (ARAP, [3]) models such process by introducing various local rigid transformations. This determines that the mathematical model of ARAP is similar with the form of pointpoint matching, i.e. a classical problem aligning two vector sets (may be of different numbers of points) [4]. This problem is renowned as the Wahba's problem [5] and point cloud registration [6] in aerospace engineering and robotics, respectively. Since the proposal of Wahba's problem in 1965, many robust and computationally efficient algorithms have been proposed. Early efforts have been paid by Shuster, Markley, Mortari who developed different solvers based on various attitude representations like quaternion, rotation matrix, Rodrigues vector and so on [7][8][9]. Recently, Wu et al. studied fast results for online determination of attitude using quaternion formulation of the Wahba's problem [10,11]. Related applications like registration using iterative closest point also highly rely on the efficiency of matching [12].
Historically, representative solutions mainly cover the singlepoint batch processing of point-set alignment. However, it is required in many applications that the solution should be either accurate, robust or continuous, so that the animated meshes can be smoothly placed. A recent method using the gradient descent algorithm (GDA) of the Wahba's problem solves this problem but is not linear [13]. It has also been proposed in [14] that some suboptimal filters can produce continuous results but the accuracy is not optimal. The target of ARAP-based mesh deformation is similar to Wahba's problem but as the amount of meshes are usually huge. Therefore, the practitioners require a new method that is not only accurate, robust, continuous but also, ARAP-friendly. Guided by this requirement, in the paper, the main contributions are: i. A geometric algebra (GA) based formulation of the kernel problem in ARAP has been proposed, which is intuitive, linear and simple. ii. It is able for us to construct continuous Newton algorithm to converge to the ideal solution within only several simple iterations. iii. ARAP and its improved variants are combined with the proposed method in a friendly manner allowing real-time performance of accurate mesh animation.
We thus study the analytical forms of all the operations required. Through derivations we are able to give linear solutions with simple matrix manipulations. It has been compared with representatives of Wahba's problem and point cloud registration and has shown that it owns the same accuracy but less computational burden. The performance of the mesh animation also proves that the accuracy belongs to the optimal ones while the computational efficiency is the highest. This paper is then structured as follows. Section 2 provides the details of the proposed derivations and solution. Section 3 investigates the application of the method to the ARAP and its variants. Section 4 gives experimental results and comparisons and the concluding remarks are drawn in Section 5.

Preliminaries
All the rotation matrices belong to the space of special orthogonal groups, i.e. for 3 × 3 case, the 3D SO(3) := R R ∈ ℝ 3 × 3 , R ⊤ R = I, det (R) = 1 . The GA allows for the geometric product of two 3 × 1 vectors a and b such that in which ⋅ and × are inner (dot) and outer (cross) products. Note that the outer product a × b = ℬ forms a bivector ℬ also in the 3D GA space G 3 . This further gives a normalised version within the scope of the even subalgebra G 3 + , which produces a group of rotor (spin) vectors We use R k to extract the kth element of R.

Proposed work
Given two sets of n points in P = {p j } j = 1 n , Q = {q j } j = 1 n , we try to minimise the following error: where {c j } j = 1 n are scalar weights such that ∑ j = 1 n c j = 1. In that form, the minimisation is non-linear in R. Following some manipulations, we can make it linear by multiplying by R on the right so the constraint RR = 1 is implicit in the energy E(R), although we will need to project R back to the manifold through normalisation. Let F be a column vector of functions which square amounts to the energy E(R) The energy E(R) can be expressed as the matrix product where F j is producing a multivector with the basis {e 1 , e 2 , e 3 , e 123 }. The critical points of the energy E(R) where it is minimised can be found solving ∇E(R) = 0, where the gradient has the following form: where J is the Jacobian matrix of F. Since the energy E(R) is purely quadratic in R, the solution for ∇E(R) = 0 can be found by solving a linear system of equations. An optimal rotor is in the null space of a particular matrix derived from J ⊤ F = 0. An easy way to obtain a solution in the null-space is using one iteration of Newton's method. Due to linearity of the equation where H is independent of R. An optimal solution can then be found by solving a linear system where R 0 is an initial rotor and g(R 0 ) is the gradient of E evaluated at R 0 . Of course, the choice of initial rotor R 0 has a particular effect on finding a local minimum. The energy E(R) is non-convex, its shape is similar to a sinusoidal wave with infinitely many points at minimum energy value, each one at rotor R i = e −(θ + iπ)B for i ∈ ℕ, being B the optimal unique attitude bivector and θ an optimal angle with minimal absolute value. The choice of R 0 leads to find a local solution close to it, the choice R 0 = [1, 0, 0, 0] ⊤ is optimal from the computational point of view and is also desirable because lead us to find solutions close to the identity rotor. Since there are infinite solutions, H is a singular matrix. A simple and efficient way to find its pseudo-inverse is to add a very small regularisation term to E(R). Our strategy is to introduce a new constant rotor R i within a new term ϵ ∥ R − R i ∥ 2 , which we can express as matrix product with Jacobian where I is the 4 × 4 identity matrix. The energy now looks like this where 0 < ϵ < 1 is constant but small (we use ϵ = 10 −6 ). The term ϵ ∥ R − R i ∥ 2 is a regularisation term that helps on finding a solution biased towards R i , so one can think R i as a constant rotor which is pointing R towards a given minimum point. With this regularisation strategy we can not only find the pseudo-inverse of H but also we can make sure to obtain a stable R that is not stuck when the minimal point is close to some far angle such ±π radians. For instance by taking R i as the previous rotor calculated with Newton iteration, we can check that R is at some minimum point, because in the next iteration we should obtain the same R. In the edge cases near to ±π where the R may not be at minimum, the previous solution R i will direct the next solution to a minimum, commonly in one to three iterations. So, regularisation term can be seen as a feedback for reaching a stable minimal R. We will see this more clearly in the algorithm. Keep in mind that we are solving the linear system just once (i.e. we are inverting H once) because the R i term is affecting only the gradient. The gradient and Hessian are now The optimal rotor is as before The Hessian matrix H resembles now a form of Tikhonov regularisation, whose inverse H −1 approaches to the Moore-Penrose pseudo-inverse as ε approaches to zero. We now need to determine how to obtain R i . An stand-alone algorithm can take R i as the previous rotor calculated with Newton iteration. But we also can take R i to be a previous solution in a simulation. We have used both choices in our experiments, with excellent results in both cases. We noticed in a simulation that the previous rotor helps to get the next one preserving the same sense of rotation. We now proceed to lie down our algorithm taking R i from previous rotor calculated with Newton iteration.

Optimal computation of Jacobians:
Let us recall that each term of the sum ∑ j = 1 n c j ∥ Rp j − q j R ∥ 2 can be expressed in matrix F j is producing a multivector with the basis {e 1 , e 2 , e 3 , e 123 }.
Solution to HΔR = g(R 0 ) looks like this where we should normalise the rotor R * for sending it back to the manifold.

Form of the Jacobians:
The Jacobian J j is 4 × 4 which four columns are the directional derivatives of ∂F j ∂e 13 = c j Thus, the leading symmetric matrix H j = J jT J j has a simple form.
The symmetric matrix H j = J jT J j has a simple form in (25) provided that (see (25)) Since H j is symmetric only 10 out of 16 elements need to be actually computed. The system now looks like this Solving just one iteration of the Newton suffice for solving the system, i.e. 19. In the iteration the symmetric matrix H = ∑ j = 1 n H j is constant, but the gradient at right-hand-side depends on R since F j depends on it. We now can proceed to optimise a little more our method. We incorporate the fixed initial guess R 0 = 1 into the gradient Notice that the first row of H j is equal to g j , so g j does not need to be explicitly calculated. Now we introduce the regularisation term F 2 and proceed to replace again R 0 on it, which look like this The final iteration is reduced to Notice that all terms are constant except d which depends on the previous iteration. Note that the matrix ∑ j = 1 n H j + ϵI does not depend on R and its inverse can be precomputed. Since the matrix (∑ j = 1 n H j + ϵI) −1 and the vector ∑ j = 1 n g j are constant, the inner loop amounts to compute a cheap matrix-vector multiplication. Also notice that since the matrix ∑ j = 1 n H j is symmetric, only ten components of J need to be actually computed and since the matrix ∑ j = 1 n H j + ϵI is symmetric its inverse is symmetric and also cheap to compute.

ARAP
Given two meshes P and Q consisting of vertices p i and q i , respectively, and directed edges p i j = p j − p i and q i j = q j − q i , the discrete ARAP energy is defined as where R 1 , …, R m ∈ SO(3) are optimal local rotations, ℰ 1 , …, ℰ m are their corresponding set of 1-ring edges and c i j are weighting coefficients, typically the familiar cotangent weights.
The main idea of this method is breaking the surface into overlapping cells ℰ i and seek for keeping the cells transformations as rigid as possible in the least squares sense. Overlap of the cells is necessary to avoid surface stretching or shearing at the boundary of the cells. The vertices of mesh P are in original position while vertices of mesh Q are the deformed vertices and the matrix R i is the best rigid transformation, in the least squares sense, relating the original and the deformed vertices. This is a non-linear optimisation problem that can be efficiently solved by a simple iterative method that solves two linear least squares sub-problems on each iteration. The first step is to consider the vertices of P and Q constant and obtain the best rigid transformation R i for each cell ℰ i . The second step is to consider the rotations R i constant and computing the optimal deformed vertices q i in the least squares sense. This can be achieved by taking the gradient of E(P, Q) w.r.t. Q and equating the result to zero, which leads to formulate a linear system of equations which can be solved by constraining the position of at least two vertices of Q. The main source of inefficiency in this method is that the first step typically involves solving a series of singular value decomposition (SVD) problems, which is slow even using the optimised solver for 3 × 3 matrices of McAdams et al.. We show how a GA approach can speed up the technique. We first change the energy E(P, Q) for using rotors instead of rotation matrices We first solve the best rotors incrementally first by applying the algorithm developed above. That is an incremental version of our algorithm that only require to solve a linear system for computing the next best rotor.
The second step is computing the optimal vertices q i . Taking the partial derivatives of E(P, Q) w.r.t. q i and equating the result to zero, we obtain Equation (34) can be expressed in matrix form as LQ = C, where L is the symmetric Laplace-Beltrami operator, Q is the column of target positions and C is the right hand side of (34). That is, a Poisson equation. Constraints of the form q i = c i are incorporated into the system by substituting the corresponding variables, i.e. erasing respective rows and columns from L and updating the right-hand side with the values c i . The system is then solved in the least squares sense Mesh deformations are achieved by repositioning the constrained vertices q i = c i , solving the linear subproblem for rotors R i , updating the right-hand-side of (34) and solving the LLS system for q i . Having a good initial guess, the convergence is typically achieved in less than ten iterations (three to four iterations already provides compelling results). The performance improvements over the ARAP surface modelling technique are significant as can be seen in the comparison. The Euclidean rotors are as efficient as quaternions (indeed they are quite the same). Rotors require less storage than rotation matrices (just four numbers) and operations such as scalar multiplication, composition of transformations and addition are more efficient than matrices.

Smooth rotation ARAP (SR-ARAP) energy
The SR-ARAP [15] energy for a smooth map between two 2- The first term in the integral is a membrane energy, and the second term is a bending energy that penalises the difference between rotations. αÂ is a weighting scalar, where Â is the area of the whole surface. Normally, the differential of a mapping of a 2-manifold is a 2 × 2 matrix, which maps tangent vectors from the parametric domain to a tangent plane at a surface point. Here, df is a 3 × 3 matrix, which maps the 3D embedding of the tangent vectors from one surface to another (which is simply the Jacobian matrix of f : ℝ 3 ↦ ℝ 3 ). The discretisation is where R 1 , …, R m ∈ SO(3) are optimal local rotations associated with the vertices; ℰ i is the set of 1-ring edges which are neighbours of ith vertex; w i j are scalar weights; and Â is the triangle mesh area, which is used to make the energy scale invariant. Here the regularisation is conducted in the form of the chordal distance ∥ R i − R j ∥ F 2 between R i and R j . The first term, the membrane energy, is similar to the ARAP discretisation, and the second term, the bending energy, penalises the difference between an edge set rotation and the rotations of neighbouring edge sets. The objective of the membrane term is to lower the distortion of an edge set (resists stretching and shearing), by keeping the map differential close to rigid. The objective of the bending term is to keep the variation in the rotations in an edge set neighbourhood low, such that the neighbourhood would transform as a unit, as much as possible.

Geometric-algebra SR-ARAP (GA SR-ARAP) energy
The GA SR-ARAP energy is shown as follows: Notice that we are expressing rotations using rotors instead of 3 × 3 rotation matrices, also notice that in the second term, the bending energy, we are using the Euclidean (L 2 ) norm instead of the Frobenius norm. Those changes are not affecting the resulting surface deformation in any noticeable way but are very convenient since the resulting optimisation problem can be solved with a much more efficient algorithm as we will show.

Local step
In the local step we consider vertices of P and Q as constants and solve for the best rotation of each cell ℰ i independently by minimising Like in the Wahba problem, in the first term the minimisation is non-linear in R i . Following Perwass, we can make it linear by multiplying by R i on the right This amounts to minimise where the constraint R i R i = 1 is now implicit, although we will need to project R i back to the rotor manifold through normalisation. Let F(R) be the vector of functions which square equals to the energy E(R) = F ⊤ F. The column vector F can be split in two blocks so that where n is the degree of the ith vertex, i.e. n = ℰ i . The energy E(R) can be expressed as the matrix product the factor w i j is absorbing the whole factor (αÂ /n)w i j for the sake of notation simplicity.
Since J 1 and J 2 are constant, the equation ∇E(R) = 0 is linear in R and the solution amounts to solve the null-space of a linear system. As before, we can solve it using only one iteration of Newton method. The Newton increment is given by ΔR = − H −1 ∇E(R).

The Hessian matrix H is
So, the full system looks like this For some initial guess R 0 , that we will choose to be R 0 = 1 as in the previous section. The rest of the work focuses on the computation of Jacobian and Hessian which is in the similar form shown in aforementioned contents. Here we define Notice that the first row of H j is equal to g j , so g j does not need to be calculated. Now the right-hand-side does not depend on R i (see (49)) . The final system to solve is

Global step
The global step is computing the optimal vertices {q i } ∈ Q E SR (P, Q) = min Taking the partial derivatives of E SR (P, Q) w.r.t. q i and equating the result to zero lead us to obtain the linear system of (34) which can be expressed in matrix form as LQ = C, where L is the symmetric Laplace-Beltrami operator, Q is the column of target positions and C is the right hand side of (34). Constraints of the form q i = c i are incorporated into the system by substituting the corresponding variables, i.e. erasing respective rows and columns from L and updating the right-hand side with the values c i . The solution is (L ⊤ L)Q = L ⊤ C, as described previously.

Local relaxation:
In the SR-ARAP method, mesh deformations are obtaining by repositioning the constrained vertices q i = c i , solving local step subproblem for rotors R i , updating the right-hand-side of (34) and solving the global step subproblem system for q i . Unlike the ARAP surface, the optimised rotations in the local step for SR-ARAP are codependent. Since we optimise each rotation independently, while fixing the others, the local step can be considered as a relaxation, thus more than one local iteration should be executed before performing the global step. At least two steps of local relaxations should be done for each global iteration, although for better convergence we should perform more local relaxations for each global iteration.

Real-time SR-ARAP
The local relaxation used to be a serious bottleneck for the SR-ARAP method. However our fast rotor estimation SR-ARAP can be used to compute as many local relaxation steps as needed at the same cost as the simple ARAP method (i.e. without slowing down the solver). Also we show how real-time performance can be achieved on SR-ARAP by adding a simple multiresolution step.

Rotor relaxation:
Recall that rotor estimation for SR-ARAP amounts to solve one linear system where H = ∑ j = 1 n H j + αÂ I , g = ∑ j = 1 n g j and d = ∑ j = 1 n d j . The only term involving the neighbour rotations is ∑ j = 1 n d j . So, for the entire local relaxation we can precompute H −1 and g. So, the computational cost of solving each relaxation step amounts to do a matrix multiplication.

Multiresolution SR-ARAP:
For achieving real-time performance, we optimise the SR-ARAP energy in a low resolution version of the input mesh and then we transfer that solution to the full resolution mesh. To obtain the low resolution mesh, we simplify the mesh using half edge collapses (i.e. the simplified mesh is a triangulation of a subset of the original vertices) while minimising the Quadrics error metric. After we obtained the optimal deformed shape on the simplified mesh, we transfer the optimised rotations to the full resolution mesh and solve the linear system of (34) using as position constraints the optimised vertices of the low resolution mesh. The transference of rotations from the low-res mesh to the highres mesh is equivalent to clustering rotations in the high resolution mesh. The clustering of rotations is based on the connectivity graph of the simplified mesh. Having the a injective map f : Q low → Q hi that maps rotations from low-res mesh Q low to hi-res mesh Q hi we generate the surjective map c: Q hi → Q low which is mapping several (equal) rotations of Q hi to one corresponding rotation of Q low . Our algorithm to construct the surjection c is a simple diffusion process: We start by mapping the rotations from low-res mesh to hi-res mesh using the injection f, then we iteratively copy the assigned rotations to the neighbouring vertices that still does not have a rotation, repeating the process until there are no vertices with neighbours without a rotation assigned. That process is used to precompute the surjection c: Q hi → Q low which is used later for efficient transferring of rotations. Let f −1 denote the inverse image of the injection f. So, the range of f −1 is the entire set of vertices in Q low and the domain of f −1 is the image of f.

Experimental results
In this section, we use some datasets to conduct registration and mesh deformation using our proposed method and other existing algorithms. In the first study, we show the accuracy consistency of the proposed method with representative methods. We use the following model to generate the point cloud data: in which r i is the reference vector and b i stand for the measurement vector in the body frame. We use the ground truth data of R and reference data as shown in Table III in [10]. Here various algorithms for registration are employed for assembled test, they are SVD method by Horn [16], Fast Linear Attitude Estimator (FLAE) [10], GDA [13], Fast Symbolic 3D Registration (FS3R) [11], Fast Analytical 3D Registration (FA3R) [17]. Each of them has been proven to be accurate, robust and computationally efficient. We evaluate various representatives using the loss function value ℒ and computation time on a personal computer of MacBook Pro 2017 with i7-CPU and MATLAB environment. The required parameters of these algorithms are: • Proposed: relative stop criteria 1 × 10 −12 .
With 12 classical cases in [10], we are able to generate the results shown in Table 1. One may observe that the accuracy of compared methods are almost identical. The reason is that these solutions are all based on the framework of (54) so the resulted analytical formulations achieve equivalent mathematical errors.
The developed method has been combined with ARAP and its variants in the aforementioned contents. Based on the energy functions defined using the GA, we are able to implement the proposed algorithm. Here, the libigl-2.1.0 library [18,19] has been invoked for high-efficiency implementation of the mesh deformation via parallel computation through the BLASX library [20]. We use the decimated-knight and cactus models to The model of knight contains 1500 vertices and simulated control points are generated using ground truth rotation. Fig. 1 shows the refined ARAP results of various algorithms. Gradient descent is introduced for classical implementation where the descent factor is set to 5 × 10 −2 . Seen from the figure, it is observed that the GA supervised ARAP can achieve quite accurate deformation after 20 iterations. However, the GDA according to its nature of convergence problem, cannot fully converge to the ground truth value. That is to say in engineering, we may need more iterations or improved Newton methods to converge to a satisfactory result. However, we need to note that this is not computationally efficient and thus should be compared with the proposed GA algorithm.
The cactus model is employed for mesh deformation and the result of GA SR-ARAP is shown in Fig. 2. We can see that under high noise level the algorithm is still capable of achieving good deformation results. This indicates the accuracy of the GA-based point cloud registration. Now we need to compare the computational efficiency of various algorithms. All these algorithms are then implemented in an unoptimised (compilation level) manner guaranteeing the fairness of all the run-time results. The termination criteria of all the algorithms is that the relative difference of the minimised energy of two successive iterations reaches 1 × 10 −8 . The details of the execution time for the cactus model with different numbers of vertices are summarised in Table 2. The GA-based ARAP and variants are more time-efficient than other representatives. We need to point out that although FLAE is analytical and computationally more efficient than the iterative GA algorithm proposed in this paper, it is not so friendly  combining with ARAP. Therefore, the convergence of the FLAEbased ARAP will require more computational burden. For illustration of the proposed method for large numbers of meshes, the Armadillo model has been chosen which contains over 43,000 vertices and 80,000 edges [21]. For ARAP, it is hard to implement real-time animation of the Armadillo model for satisfactory update frequency. Using the proposed GA-ARAP, it is able to speed up the mesh deformation to a large extent. Fig. 3 shows the details of the reference and the deformation result from the proposed GA-ARAP.
The developed algorithm maintains good accuracy and can achieve very high visualisatoin frequency of 60 fps on a typical personal computer with i7-4core CPU. For conventional ARAP, since the computational burden is relatively higher, it can only reach the update rate of 12 fps at most. The good computational behaviour of the proposed method provides the possibility of performing highly real-time mesh deformation such as games and high-speed visualisation.

Conclusion
In this paper, we solve the problem of estimating the best rotation for the alignment of two sets of corresponding 3D points. It is based on solving the linear equations derived from the formulation of the problem in Euclidean GA. The method is fast, robust to noise, accurate, simple and GPU friendly. Applications and experimental validation of mesh deformation are presented for description of the performance of the proposed algorithm. The results show that the slightly losing accuracy of the proposed method leads to advance in execution time consumption. The designed algorithms are ARAP-friendly and may be helpful for real-time mesh deformation in related applications.