Efficient two step optimization for large embedded deformation graph based SLAM

TL;DR

This paper introduces a two-step optimization method for large-scale embedded deformation graph SLAM that maintains constant estimation speed regardless of map size by exploiting graph sparsity.

Contribution

It proposes a novel decoupling approach that reduces computation complexity from quadratic to near constant, enabling scalable large-scale deformable SLAM.

Findings

Achieves near constant computation time as map size increases.

Reduces estimation complexity from O(n^2) to close to O(1).

Validated effectiveness in large-scale deformation scenarios.

Abstract

Embedded deformation nodes based formulation has been widely applied in deformable geometry and graphical problems. Though being promising in stereo (or RGBD) sensor based SLAM applications, it remains challenging to keep constant speed in deformation nodes parameter estimation when model grows larger. In practice, the processing time grows rapidly in accordance with the expansion of maps. In this paper, we propose an approach to decouple nodes of deformation graph in large scale dense deformable SLAM and keep the estimation time to be constant. We observe that only partial deformable nodes in the graph are connected to visible points. Based on this fact, sparsity of original Hessian matrix is utilized to split parameter estimation in two independent steps. With this new technique, we achieve faster parameter estimation with amortized computation complexity reduced from O(n^2) to closing O(1). As a result, the computation cost barely increases as the map keeps growing. Based on our strategy, computational bottleneck in large scale embedded deformation graph based applications will be greatly mitigated. The effectiveness is validated by experiments, featuring large scale deformation scenarios.

Equations

$${\tilde{\mathbf{v}}}_i={\mathbf{R}}_c\sum _{j=1}^m{\omega }_j({\mathbf{v}}_i)[{\mathbf{A}}_j({\mathbf{v}}_i-{\mathbf{g}}_j)+{\mathbf{g}}_j+{\mathbf{t}}_j]+{\mathbf{T}}_c,$$

$${\omega }_j({\mathbf{v}}_i)=1-||{\mathbf{v}}_i-{\mathbf{g}}_j||/d_{max},$$

$$\begin{array}{rl}\underset{{\mathbf{R}}_c,{\mathbf{T}}_c,{\mathbf{A}}_1,{\mathbf{t}}_1...{\mathbf{A}}_m,{\mathbf{t}}_m}{\mathrm{argmin}}& {\omega }_{rot}{E}_{rot}+{\omega }_{reg}{E}_{reg}+{\omega }_{data}{E}_{data}.\end{array}$$

$$E_{rot}=\sum _{j=1}^{m}Rot({\mathbf{A}}_j),$$

$$Rot(\mathbf{A})=({\mathbf{c}}_1^T\cdot {\mathbf{c}}_2{)}^2+({\mathbf{c}}_1^T\cdot {\mathbf{c}}_3{)}^2+({\mathbf{c}}_2^T\cdot {\mathbf{c}}_3{)}^2+$$

$$({\mathbf{c}}_1^T\cdot {\mathbf{c}}_1-1{)}^2+({\mathbf{c}}_2^T\cdot {\mathbf{c}}_2-1{)}^2+({\mathbf{c}}_3^T\cdot {\mathbf{c}}_3-1{)}^2.$$

$$E_{reg}=\sum _{j=1}^m\sum _{k\in \mathbb{N}(j)}{\alpha }_{jk}||{\mathbf{A}}_j({\mathbf{g}}_k-{\mathbf{g}}_j)+{\mathbf{g}}_j+{\mathbf{t}}_j-({\mathbf{g}}_k+{\mathbf{t}}_k)|{|}^2.$$

$$E_{data}=$$

$$\begin{array}{r}\underset{3m\times n}{\mathbf{M}}=\left[\begin{array}{ccc}\cdots & & \cdots \\ \cdots & {\mathbf{e}}_{n,1}& \cdots \\ \cdots & & \cdots \\ \cdots & {\mathbf{e}}_{n,2}& \cdots \\ \cdots & & \cdots \\ \cdots & {\mathbf{e}}_{n,3}& \cdots \\ \cdots & & \cdots \\ \cdots & {\mathbf{e}}_{n,4}& \cdots \\ \cdots & & \cdots \end{array}\right]\end{array}\underset{m\times n}{\mathbf{C}}=\left[\begin{array}{ccc}\cdots & & \cdots \\ \cdots & {\omega }_{\mathbf{N}(n,1)}& \cdots \\ \cdots & & \cdots \\ \cdots & {\omega }_{\mathbf{N}(n,2)}& \cdots \\ \cdots & & \cdots \\ \cdots & {\omega }_{\mathbf{N}(n,3)}& \cdots \\ \cdots & & \cdots \\ \cdots & {\omega }_{\mathbf{N}(n,4)}& \cdots \\ \cdots & & \cdots \end{array}\right],$$

$${\bm{\Lambda}}=\left(\begin{array}[]{ccc}{\mathbf{A}_{1}}&\cdots&{\mathbf{A}_{m}}.\\ \end{array}\right)$$

$${\mathbf{T}}=\left(\begin{array}[]{ccc}{\mathbf{t}_{1}}+{\mathbf{g}_{1}}&\cdots&{\mathbf{t}_{m}}+{\mathbf{g}_{m}}\\ \end{array}\right).$$

$$E_{data}=||{\mathbf{R}}_{\mathrm{c}}\cdot [\mathbf{\Lambda }\cdot \mathbf{M}+\mathbf{T}\cdot \mathbf{C}]+{\mathbf{T}}_{\mathrm{c}}\otimes 1-\tilde{\mathbf{P}}|{|}_F^2.$$

$$E_{data}=||{\mathrm{\mathbf{R}_{c}}}\cdot\left(\begin{array}[]{cc}{\mathrm{\mathbf{\Lambda}}}&{\mathbf{T}}\\ \end{array}\right)\cdot\left(\begin{array}[]{c}{\mathbf{M}}\\ \mathbf{C}\end{array}\right)+{\mathrm{{\mathbf{T}_{c}}}}\otimes\mathbf{1}-\mathbf{\tilde{P}}||_{F}^{2}.$$

$$E_{data}=||{\mathbf{R}}_{\mathrm{c}}[\mathbf{\Pi }\mathbf{\Phi }{]}^T+{\mathbf{T}}_{\mathrm{c}}\otimes 1-\tilde{\mathbf{P}}|{|}_F^2.$$

$$E_{data}$$

$$\displaystyle=||{\mathrm{\mathbf{R}_{c}}}[\left(\begin{array}[]{cc}{\bm{\Pi^{{}^{\prime}}}}&\ {\mathbb{O}}\\ \end{array}\right)\left(\begin{array}[]{c}{\mathbf{\bm{\Phi}}}_{1}\\ {\mathbf{\bm{\Phi}}_{2}}\end{array}\right)]^{T}+{\mathrm{{\mathbf{T}_{c}}}}\otimes\mathbf{1}-\mathbf{\tilde{P}}||_{F}^{2},$$

$$\begin{array}{rl}\underset{{\mathbf{R}}_c,{\mathbf{T}}_c,{\mathbf{\Phi }}_1}{\mathrm{argmin}}& {\omega }_{rot}{\tilde{E}}_{rot}+{\omega }_{reg}{\tilde{E}}_{reg}+{\omega }_{data}{E}_{data}\end{array}$$

$$\begin{array}{rl}\underset{{\mathbf{\Phi }}_2}{\mathrm{argmin}}& {\omega }_{rot}{E}_{rot}+{\omega }_{reg}{E}_{reg}\end{array}$$

$${\color[rgb]{1,0,0}\bm{\mathcal{J}}=\begin{bmatrix}\frac{\partial\mathbf{J}_{1}}{\partial\mathbf{x}_{c}}&\frac{\partial\mathbf{J}_{1}}{\partial\mathbf{x}_{f}}\\ \frac{\partial\mathbf{J}_{2}}{\partial\mathbf{x}_{c}}&\mathbb{O}\end{bmatrix}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\begin{bmatrix}\mathbf{J}_{1c}&\mathbf{J}_{1f}\\ \mathbf{J}_{2c}&\mathbb{O}\end{bmatrix}}$$

$${\color[rgb]{1,0,0}\bm{\mathcal{H}}=\begin{bmatrix}\mathbf{J}_{1c}^{T}\mathbf{J}_{1c}+\mathbf{J}_{2c}^{T}\mathbf{J}_{2c}&\mathbf{J}_{1c}^{T}\mathbf{J}_{1f}\\ \mathbf{J}_{1f}^{T}\mathbf{J}_{1c}&\mathbf{J}_{1f}^{T}\mathbf{J}_{1f}\end{bmatrix}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\begin{bmatrix}\mathbf{\Lambda}_{cc}&\mathbf{\Lambda}_{cf}\\ \mathbf{\Lambda}_{cf}^{T}&\mathbf{\Lambda}_{ff}\end{bmatrix}}$$

$${\color[rgb]{1,0,0}\begin{bmatrix}\mathbf{\Lambda}_{cc}&\mathbf{\Lambda}_{cf}\\ \mathbf{\Lambda}_{cf}^{T}&\mathbf{\Lambda}_{ff}\end{bmatrix}\begin{bmatrix}\mathbf{x}_{c}\\ \mathbf{x}_{f}\end{bmatrix}=\begin{bmatrix}\mathbf{J}_{1c}^{T}\mathbf{F}+\mathbf{J}_{2c}^{T}\mathbf{F}\\ \mathbf{J}_{1f}^{T}\mathbf{F}\end{bmatrix}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\begin{bmatrix}\mathbf{y}_{c}\\ \mathbf{y}_{f}\end{bmatrix}}$$

$${\color[rgb]{1,0,0}\begin{bmatrix}\mathbf{\Lambda}_{cc}-\mathbf{\Lambda}_{cf}\mathbf{\Lambda}_{ff}^{-1}\mathbf{\Lambda}_{cf}^{T}&\mathbb{O}\\ \mathbf{\Lambda}_{cf}^{T}&\mathbf{\Lambda}_{ff}\end{bmatrix}\begin{bmatrix}\mathbf{x}_{c}\\ \mathbf{x}_{f}\end{bmatrix}=\begin{bmatrix}\mathbf{y}_{c}-\mathbf{\Lambda}_{cf}\mathbf{\Lambda}_{ff}^{-1}\mathbf{y}_{f}\\ \mathbf{y}_{f}\end{bmatrix}}$$

Full text

I INTRODUCTION

Simultaneous Localization and Mapping (SLAM) is a technique widely used in robotics for pose estimation and environment mapping. While SLAM in the rigid scenario is quite mature, SLAM with deformations is still new to the community. Several attempts in this direction include RGBD camera based deformable human body reconstruction and RGB camera based soft-tissue SLAM in Minimally Invasive Surgery (MIS). The key challenge in developing deformable SLAM is how to define and formulate the deformation.

As an early stage, researchers assume deformation scenario is technically similar to SLAM in dynamic environment. By assuming a large amount of target features to be static, approaches like [1] and [2] were proposed, adopted conventional extended Kalman filter (EKF) SLAM and Parallel Tracking and Mapping (PTAM) by applying some thresholds to separate rigid and non-rigid feature points. More recent work ORB-SLAM [3] have been proposed and modified in [4] [5] [6] [7]. The above works show that extending traditional rigid SLAM to non-rigid scenarios is possible. However, map reconstruction based on the rigid assumption inevitably leads to gaps within rigid patches from different perspectives. The more extensive deformation of the environment, the larger gaps exist in the recovered map.

A better solution is to model the deformation. Recognizing map deformation modeling is vital and indeed many researches turn to the computer vision society for efficient deformation description. One of the most well-known and well-studied techniques is the Embedded Deformation (ED) [8], which was initially proposed for designing smooth character motions in cartoons. ED graph makes use of a sparsely interconnected scattered nodes with attributes of local rigid transformation. By mixing these local rigid transformations with weights, a global deformation is simulated in a discrete but smooth form. ED graph has been widely applied in RGBD and stereo SLAM [9] [10] and 3D human reconstructions [11] [12] (discrete hierarchical warp grid but very similar to ED graph) [13] [14] [15]. The ED graph based SLAM formulation is able to describe the deformation of the environment in addition to camera pose, which blossoms with the application like [11] [12] [13] [15] for human motion reconstruction and [10] for large scale SLAM in MIS.

In the meantime, for more challenging monocular equipments, Finite Element Method (FEM) [16] [17] and Structure from Template (SfT) [18] [19] were proposed to simulate deformation. Both methods use a predefined grid or triangular mesh template to describe the soft-tissue. However, to the best of our knowledge, it is hard to apply these methods when the map is incrementally built and no complete real-time implementation has demonstrated how it will effectively be applied in large scenarios.

Overall, ED graph based formulation is the most applicable and widely used approach for modeling deformations. Numerous real-time dense SLAM systems have validated its effectiveness, particularly in batch scenarios, with fast sequential or parallel implementation. However, ED graph based formulation also comes with disadvantages and limitations. One major issue in ED graph is that when a new observation is incorporated, the number of nodes increases dramatically, posing heavy computational burden. Little attention has been paid in the field of truncated signed distance function based 3D human reconstruction because the target size, as well as map extent, are predefined. State estimation and map updating are all confined within a volume. In more general cases, however, as reported in [10], when reconstructing geometry without a predefined volume, the environment size is unbounded due to the sharp growth of the graph, which eventually implies an amortized $O(n^{2})$ complexity with respect to the number of nodes in the graph. In a word, optimizing an expanding ED graph in an unconstrained space significantly limits the performance of the system.

Even though little is known in speeding up ED graph based SLAM systems, numerous publications are dedicated to the graph based optimization in rigid scenarios. A typical graph optimization problem models robot poses (or landmark positions) as nodes in the graph, while the edges encode the relative measurement between connected nodes [20]. Graph sparsification is the most widely applied technique to marginalize subsets of nodes [21] [22] [23]. The key process in this topic is to sparsify edges and marginalize nodes based on indicators like divergence [24]. We will show in this paper that marginalization methods are of great value to enhance efficiency in ED based SLAM. In this paper, we will mitigate the $O(n^{2})$ computational bottleneck encountered in the expanding environment, for ED graph based SLAM as reported in [10]. We analyze the spatial relationship between ED graph and observation and reveal the inherent sparsity pattern of the Hessian matrix. With this discovery, we classify ED nodes into points relevant (PR) and points irrelevant (PI) nodes and propose a decoupled optimization strategy.

In this work we convert the ED graph based formulation into a matrix form, which unveils a sparsity pattern that can be further exploited by restructuring the Jacobian matrix to benefit efficient marginalization. Based on this insight, we develop a method which split ED graph optimization into two steps and a lossy decoupled optimization strategy is proposed. The accuracy and effectiveness are tested on both the ED based geometry deformation as well as ED based SLAM application. The proposed strategy is validated on the MIS-SLAM framework [10] with and without the proposed strategy. For more details on the MIS-SLAM framework, please refer to the flowchart of in [10].

II Revisit general ED graph based SLAM

A ED graph, as other graph representations, comprises nodes and edges. See Fig. 1 for an example. We use the term ‘node’ to describe a ED node, and ‘vertex’ to describe a feature point. Each ED node $j$ defines a local rigid motion modeled by ${\mathbf{A}_{j},\mathbf{t}_{j}}$, where $\mathbf{A}_{j}$ $\in\mathbb{R}^{3\times 3}$ is an affine matrix and $\mathbf{t}_{j}$ $\in\mathbb{R}^{3}$ is a translation vector. Given a vertex $\mathbf{v}_{i}$, neighbouring ED nodes define a deformation field indicating how $\mathbf{v}_{i}$ is deformed to target position $\mathbf{\tilde{v}}_{i}$. For a mesh based geometry, shape deformation is equivalent to translate vertices in the ED graph.

In practice, to avoid over-parameterization, the ED node $\mathbf{g}_{j}$ are deliberately positioned to minimize the size of parameters. For example in [8] [11] [13] [25], ED nodes are generated by uniformly downsampling the original model. The deformation between the point $\mathbf{v}_{i}$ and $\tilde{\mathbf{v}}_{i}$ is modeled as a transformation in the following form:

[TABLE]

where $\mathbf{R}_{c}$ and $\mathbf{T}_{c}$ are global rotation and translation relating to camera motion, $m$ is the number of neighboring nodes. The global camera pose $[\mathbf{R}_{c},\mathbf{T}_{c}]$ is not considered in the original ED work [8], however added later on in [11] [13] [25]. $\omega_{j}(\mathbf{v}_{i})$ is a weight that quantifies the contribution of the ED node $j$ to the vertex $\mathbf{v}_{i}$. We choose the number of neighboring nodes to $m=4$ and define the weight:

[TABLE]

where $d_{max}$ is the maximum distance from the vertex $\mathbf{v}_{i}$ to the neighbouring $k+1$ nearest ED node.

Eq. (1) shows how one vertex is transformed with the ED graph. Conversely, if the source ($\mathbf{v}_{i}$) and target ($\tilde{\mathbf{v}}_{i}$) point pairs are given, the parameters of ED graph ($\mathbf{A}_{j}$ and $\mathbf{t}_{j}$) can also be inferred. Thus, the ED graph ED based SLAM formulation is to infer the ED graph parameter from a cluster of arbitrary source and target points pairs. Then, ED graph is applied to deform the whole model. To the best of our knowledge, all ED graph based SLAM formulation consist of at least three terms: Rotation constraint, regularization constraint and a penalty term on the distance of deformed source and target point pairs. Some methods introduce more innovating terms, like visual hull term in [13], key features in [12] which are beyond the scope. A ED graph based SLAM is to minimize:

[TABLE]

The major key point pairs matching term is $E_{data}$, while two soft constraints $E_{rot}$ and $E_{reg}$ regulate deformation. $E_{rot}$ hauls the affine matrix close to $SO(3)$ by minimizing the column vectors $\mathbf{c}_{1}$, $\mathbf{c}_{2}$ and $\mathbf{c}_{3}$ in:

[TABLE]

[TABLE]

Regularization is to ensure smoothness of the deformed surface, the motion field generated from ED nodes should be consistent; otherwise, the deformed shape will be fragmented and distorted. This is a general practice named ‘as-rigid-as-possible’ and defined in following form:

[TABLE]

Similar to [8], $\alpha_{jk}$ is the weight quantifying the influences of the neighboring ED nodes. $\mathbb{N}(j)$ is the set of all connected nodes to node $j$. Normally, the number of neighboring nodes is 4.

To solve geometrical model to frame registration, ‘back-projection’ formulation is adopted as a substitution to iterative closest point (ICP). Readers may refer to [13] [11] [25]. For simplicity, we introduce the basic source and target key point pairs described by [8]. Let $\mathbf{\tilde{v}}_{i}$ and $\mathbf{\hat{v}}_{i}$ be the pairs of source-target points (defined in Eq. (1)). Normally predefined in interactive phase, these key points define how model is to be deformed. These key points define a data term which minimizes a distance as:

[TABLE]

III Efficient two step optimization

III-A Matrix form of ED graph deformation

To fully exploit the structure we rewrite the data term Eq. (7) in a matrix form for the convenience of sparsity analysis. Let us consider a group of predefined key source points ${\mathbf{P}}=[\mathbf{v}_{1}...\mathbf{v}_{n}]$ and the key target points ${\mathbf{\tilde{P}}}=[\tilde{\mathbf{v}}_{1}...\tilde{\mathbf{v}}_{n}]$ in the matrix form. According to Eq. (1), each point ${\mathbf{v}_{i}}$ is deformed by its 4 neighboring nodes. Thus we define two matrices $\mathbf{M}$ and $\mathbf{C}$ (Eq. (8)):

[TABLE]

where $\mathbf{e}_{n,k}=\omega_{\mathbf{N}(n,k)}*(\mathbf{v}_{i}-\mathbf{g}_{\mathbf{N}(n,k)})$ and $\mathbf{N}(i,j)(j=1,2,3,4)$ is the set of neighboring node to point $\mathbf{v}_{i}$. In matrices $\mathbf{M}$ and $\mathbf{C}$, note that non-zero elements are not aligned. Each column only has 4 non-zero elements (neighboring nodes). The sum of each column in matrix $\mathbf{C}$ is 1 due to the point to node topology. In Eq. (7), a source point $\mathbf{v}_{i}$ is transformed by its 4 neighboring node which yields 4 non-zero elements in $\mathbf{M}$ and $\mathbf{C}$. The sum of all weight $\omega_{j}(\mathbf{v}_{i})$ is 1. Note that different source points have different topology, thus the location of non-zero elements are not aligned well in each column. We arrange the parameters of ED nodes $\mathbf{A}_{i}$ and $\mathbf{t}_{i}$ in the following form:

[TABLE]

[TABLE]

Then Eq. (7) takes the following form:

[TABLE]

where $\otimes$ is the kronecker product. $\mathbf{1}$ is $1\times n$ vector of ones. And $||\cdot||_{F}^{2}$ is the Frobenius norm. Eq. (11) can be written compactly in the following form:

[TABLE]

If we further define ${\bm{\Pi}}=[{\mathbf{M}}^{T}\ {\mathbf{C}}^{T}]$ and ${\mathbf{\bm{\Phi}}}=[{\mathbf{\bm{\Lambda}}}\ {\mathbf{T}}]^{T}$, Eq. (12) can be formed as follows:

[TABLE]

The property of Jacobian of $E_{data}$ is determined by ${\bm{\Pi}}$.

III-B Sparsity in ED graph formulation

It is natural to solve Eq. (13) in batch. As the number of vertices increases, the dimension of Jacobian relating to state $\mathbf{\bm{\Phi}}$ increase dramatically. Fortunately, we can explore the fact that only partial ED nodes are connected to currently visible points, due to the limited field of view (FOV) of the camera. See Fig. 2 as an example of constant FOV of depth image. The model keeps expanding while the target depth remains in small size. A typical ED node and target depth relationship is illustrated in Fig. 2(d); 2/3 of the nodes are not within depth FOV resulting no contribution to $E_{data}$. Fig. 3 shows a typical Jacobian of the cost function. In $E_{data}$ block, the shadow region indicates nodes connected to points (PR nodes) while zero block shows the nodes (PI nodes) connecting to inactive points. In this paper, we will exploit the sparity of PR nodes in the zero bocks.

The same sparsity also applies to Eq. (8). By rearranging matrix ${\bm{\Pi}}$ from Fig. 3(a) to Fig. 3(b), we achieve a new Jacobian with zero block. Using this new matrix, Eq. (13) writes:

[TABLE]

where ${\mathbf{\bm{\Phi}}}_{1}$ contains $\mathbf{A}_{j}$ and $\mathbf{t}_{j}$ of PI node set and ${\mathbf{\bm{\Phi}}_{2}}$ contains $\mathbf{A}_{j}$ and $\mathbf{t}_{j}$ of PR node set. $\bm{\Pi^{{}^{\prime}}}$ is the subset of $\bm{\Pi}$ relating to PI nodes (the shadow region in Fig. 2).

III-C Lossy Two level optimization

Explained in Section III-B, the size of the PR nodes is almost constant due to the limited size of depth image in expanding scenario. For instance, the point cloud generated from Hamlyn dataset [26] (grabbed from monitor) is $320\times 240=76800$ at most. As the model grows, the total number of nodes in ED graph is increasing but the number of nodes is almost constant. Taking advantage of Eq. (14), the optimization can be divided into two levels: the optimization of PR nodes ${\mathbf{\bm{\Phi}}}_{1}$ and the optimization of the rest PI nodes ${\mathbf{\bm{\Phi}}_{2}}$.

In specific, we first optimize (${\mathbf{\bm{\Phi}}}_{1}$, ${\mathrm{\mathbf{R}_{c}}}$ and ${\mathrm{{\mathbf{T}_{c}}}}$) by fixing ${\mathbf{\bm{\Phi}}_{2}}$ in (Level I) optimization, to obtain an estimation of ${\mathbf{\bm{\Phi}}}_{1}$, ${\mathrm{\mathbf{R}_{c}}}$ and ${\mathrm{{\mathbf{T}_{c}}}}$. Then the value of the parameters obtained from level I, will be fixed in Level II, together with the two soft constraints $E_{rot}$ and $E_{reg}$ to optimize the parameter $\Phi_{2}$. We explicitly enforce this idea by formulating following energy function:

[TABLE]

[TABLE]

Eq. (15) and Eq. (16) are the Level I and Level II energy functions, where $\tilde{E}_{rot}$ and $\tilde{E}_{reg}$ are the curtailed energy function of $E_{rot}$ and $E_{reg}$ containing $\mathrm{\bm{\Phi_{1}}}$. In other words, the size of Eq. (15) is only related to the size of PR nodes. Therefore, the computational complexity in Level I is reduced from $O(n^{2})$ to constant $O(1)$ thanks to the constant size of ${\mathbf{\bm{\Phi}}}_{1}$. Although optimizing Level II is still $O(n^{2})$, considering the scale of data term $E_{data}$ is far larger than the rest, the computational cost in Level II is relatively low. Note that the new strategy keeps the time consuming step Level I constant while Level II still $O(n^{2})$. But the size of Level II is almost negligible comparing with Level I.

III-D Connection with Marginalization and Information Loss

In this section, we will draw the connection of the proposed two level optimization method with an exact marginalization based method. The analysis will show that the information loss is very low, illustrating the feasibility of the decoupled optimization Eq. (15) and Eq. (16).

When generating Eq. (14), the Jacobian shown in Fig. 3 is re-ordered by classifying $[\mathbf{A}_{1},\mathbf{t}_{1}...\mathbf{A}_{m},\mathbf{t}_{m}]$ into PR nodes ${\mathbf{\bm{\Phi}}}_{1}$ and PI nodes ${\mathbf{\bm{\Phi}}_{2}}$. The state in cost function $[\mathbf{R}_{c},\mathbf{T}_{c},\mathbf{A}_{1},\mathbf{t}_{1}...\mathbf{A}_{m},\mathbf{t}_{m}]$ are classified as $\mathbf{x}_{c}(\mathbf{R}_{c},\mathbf{T}_{c},\mathbf{A}_{1},\mathbf{t}_{1}...\mathbf{A}_{k},\mathbf{t}_{k})$ and $\mathbf{x}_{f}(\mathbf{A}_{k+1},\mathbf{t}_{k+1}...\mathbf{A}_{m},\mathbf{t}_{m})$. Fig. 3 shows the Jacobian in the new order. The first two term are combined due to their sparsity because $E_{rot}$ and $E_{reg}$ are constraints between ED nodes, i.e. unrelated to feature points. The only full block in Fig. 3 is $E_{data}$ w.r.t $\mathrm{\bm{\Phi_{1}}}$ (shadow region), in specific $\frac{\partial\mathbf{J}_{2}}{\partial\mathbf{x}_{c}}$. Let us write down the Jacobian and Hessian as,

[TABLE]

[TABLE]

The density of $\frac{\partial\mathbf{J}_{2}}{\partial\mathbf{x}_{c}}$ is the only dense block in Hessian $\bm{\mathcal{H}}$. Let us use Schur complement [20] from classic rigid SLAM and separate the optimization as follows:

[TABLE]

[TABLE]

By using Schur complement, we can solve $\mathbf{x}_{c}$ independent of $\mathbf{x}_{f}$. The computation of $\mathbf{x}_{c}$ (including $\mathrm{\mathbf{R}_{c}}$, $\mathrm{\mathbf{T}_{c}}$ and ${\mathbf{\bm{\Phi}}}_{1}$) is constant (explained in Section III-C). After solving $\mathbf{x}_{c}$, the optimization of $\mathbf{x}_{f}$ is quite cheap as the $\bm{\Lambda}_{ff}=(\frac{\mathbf{J}_{1}}{\partial\mathbf{x}_{f}})^{T}(\frac{\mathbf{J}_{1}}{\partial\mathbf{x}_{f}})$ is only related to $E_{rot}$ and $E_{reg}$. The sparsity of Hessian and small number of nodes ($n\gg m$) makes the time of solving $\mathbf{x}_{f}$ much less.

Fig. 4 is the Jacobian of $E_{rot}$ and $E_{reg}$ relating to all nodes. The first term $E_{rot}$ is the sum of error of affine transformation (Eq. (4)) making the Jacobian strictly diagonal. The second term $E_{reg}$ defines the transformation error among ED nodes (Eq. (6)). The major part $\mathbf{A}_{j}(\mathbf{g}_{k}-\mathbf{g}_{j})+\mathbf{g}_{j}+\mathbf{t}_{j}$ is also within one node $j$ except the very last $-(\mathbf{g}_{k}+\mathbf{t}_{k})$. The last variable $\mathbf{t}_{k}$ makes the Jacobian not strictly diagonal (please refer to Fig. 4(b)). We come up with an idea that by ignoring Jacobian of $\mathbf{t}_{k}$, and reordering two diagonal Jacobians, $[\frac{\partial{\mathcal{J}}_{1}}{\partial\mathbf{x}_{c}}\frac{\partial{\mathcal{J}}_{1}}{\partial\mathbf{x}_{f}}]$ becomes diagonal and $\bm{\Lambda}_{cf}=(\frac{\partial{\mathcal{J}}_{1}}{\partial\mathbf{x}_{c}})^{T}(\frac{\partial{\mathcal{J}}_{1}}{\partial\mathbf{x}_{f}})=0$. In this case, Eq. (20) can be separated into two equations, with respect to $\mathbf{x}_{c}$ and $\mathbf{x}_{f}$ respectively, which also correspond to the Level I and Level II of the proposed method. In a word, our two level lossy method can be considered as an approximation of Eq. (21), under the assumption that $\bm{\Lambda}_{cf}=0$ roughly stands.

Fig. 5 visualizes the feasibility of the lossy decoupled optimization approach in geometry. PR nodes (green) are the only nodes connected to visible points and contribute to $E_{data}$. All PI nodes (purple) merely share edges with PR nodes and are passively deformed according to the behaviors of PR nodes. Equivalently, the inter-nodes relations in the Jacobian of $E_{reg}$ (Fig. 4(b)) shows these connections (Fig. 4(b)). In view of this, our lossy decoupled optimization approach first optimize PR nodes and then estimate PI nodes.

In conclusion, solving energy function Eq. (3) by ignoring Jacobian of $\mathbf{t}_{k}$ is equivalent to the proposed decoupled optimization in Eq. (15) and Eq. (16).

The information loss of the proposed approach is relatively low. Fig. 5 illustrates the connection between PR nodes and PI nodes is weak on the boundary, in contrast with the dense connections among PR nodes. It also demonstrates how the connection between PR and PI nodes are removed in Level I, and the connection between $\frac{\partial\mathbf{J}_{1}}{\partial\mathbf{x}_{c}}$ and $\frac{\partial\mathbf{J}_{1}}{\partial\mathbf{x}_{f}}$ are removed resulting in $\bm{\Lambda}_{cf}=(\frac{\partial\mathbf{J}_{1}}{\partial\mathbf{x}_{c}})^{T}(\frac{\partial\mathbf{J}_{1}}{\partial\mathbf{x}_{f}})=\mathbb{O}$. The information between two PR nodes is strong while that among the PI nodes is weak. The weak information is neglected in Level I, attributing to relatively low information loss in optimization process.

IV Results and discussion

Our goal is to have both qualitative assessments and quantitative comparisons between the original ED graph optimization and the proposed method. For qualitative comparison, we show the sacrificed accuracy has few impacts on final reconstructed map. A dolphin model is downloaded from turbosquid (https://www.turbosquid.com). w.r.t quantitative test in SLAM, both methods are compared on a synthetic dataset and datasets from [26], where we choose three in-vivo stereo videos with deformation and rigid scope movement. The experiments are conducted on the same hardware and sorfware setting of MIS-SLAM [25]. The module of state estimation in MIS-SLAM is modified to the proposed approach. Note that an iterative solver, i.e. the preconditioned conjugate gradient method, is employed to solve the resulting linear systems, as it provides a way of parallel computing on GPU.

IV-A Qualitative ED deformation comparisons

Fig. 6 shows a qualitative comparison. It is not aiming as a proof of the superiority of the proposed method over the original ED graph method. [8] has already claimed real-time implementation on CPU as well as very nice results. Aiming at speeding up deformable SLAM application, the qualitative result of our lossy decoupled approach is not comparable to the batch estimation of ED. However, we want to illustrate that the proposed method can achieve a similar result and the difference is not visible to the naked eye or difficult to make out. Fig. 6 confirms that the deformed shape performed by our approach do actually have similar result. It looks similar to the original ED graph partially due to simple topology of dolphin. Other complex models like human results in a visible but not very apparent difference. Although the decoupled optimization works well, the results can be much worse than ED when the nodes are too sparse. The deformation is dependent on PR nodes and the insufficiency of PR nodes (or nodes in conventional one step ED) causes vertex deformation that the expected target is not reached.

IV-B Time complexity comparisons

Fig. 7 show a tiny one-step simulation and the result. Tractable time consumption remains small because PR node does not change.

We compare the original MIS-SLAM [25] with the improved version. Fig. 8 illustrates the running time for three Hamlyn datasets (model 6, 20 and 21). In all scenarios decoupled optimization yields better efficiency than batch processing especially in case of long term process (the last dataset in Fig. 8). In the first few steps, the robot is steady and ED graph is not expanding significantly. This attributes to the similar processing time in the first few hundred steps. When the robot moves, the ED graph expands intensively and processing time increases abruptly in state optimization. In view of this, by limiting the node graph, decoupled optimization keeps time consumption stable due to the constant PR node scale. The attached video shows the range of movement in model 6 is much smaller than model 20 and 21. That’s the main reason the proposed method does not improve model 6 too much. As the environment gets larger, our approach keeps much lower time consumption.

IV-C Accuracy comparisons

The lossy decoupled optimization strategy inevitably attributes the loss of accuracy. Section IV-A shows the quality of the deformed map is well preserved in ED deformation process. Moreover, we compare the proposed method and the original one on the same parameters and weights of terms in SLAM application (MIS-SLAM). Different from arbitrary key points matching in ED deformation formulation (Eq. (7)), in SLAM application the $E_{data}$ is in the form of model-to-depth scan matching like point-to-plane ICP. For quantitative validation, we measure the point-to-point distance of deforming map and target scan. For direct validation to ground truth, three synthetic data (heart, right kidney, and stomach) are generated by deliberately deforming models from CT scanned phantom. As a compliment, the three laparoscopy datasets from Hamlyn are tested, but only the back-projection error in each iteration is available since there’s no ground truth. In the batch approach, the average distance of back-projection registration of the three simulation scenarios is 0.18mm (model 6), 0.13mm (model 20) and 0.12mm (model 21). While dataset with ground truth (Hamlyn dataset 10/11) achieves 0.08mm, 0.21mm (Average errors). With our decoupled optimization approach, we achieve 0.31mm (model 6), 0.26mm (model 20), 0.22mm (model 21) and 0.14mm, 0.29mm errors. The attached video shows no big difference in terms of structure and texture.

We also test the average error. On top of the in-vivo dataset, synthetic data are generated. Please refer to the attached video for more details. Roughly speaking, the proposed method sacrifices $32\%$ accuracy in exchange for $50\%$ speed gains.

V Conclusion

We propose a decoupled approach for ED graph optimization that reduces computational complexity from $O(n^{2})$ to near $O(1)$. The decoupled optimization structure achieves faster computation in scenario of expanding environment. Our strategy sacrifices small amount of accuracy in exchange for near-constant processing speed. The constant computation complexity of the lossy strategy should have great potential in ED graph based SLAM in unbounded map scenario.

The node marginalization strategy in this paper, however, is straightforward, in the sense that it only classifies nodes based on the node-vertex connectivity. It’s reasonable because different from the pose graph, ED graph is parallelized in GPU since the time consumption requirement is more strict. However, it remains to be of great interest to test if pose graph pruning method like Kullback–Liebler divergence outperforms the proposed work while remains acceptable consumption in GPU environment.