用于曲面配准的非刚体ICP算法最优化步骤

B. Amberg, S. Romdhani, and T. Vetter, “Optimal Step Nonrigid ICP Algorithms for Surface Registration,” in 2007 IEEE Conference on Computer Vision and Pattern Recognition, Minneapolis, MN, USA, 2007, pp. 1–8.

Overview

配准是将template warp到target上。该算法使用局部仿射性进行正则化，假设每个顶点都对应一个仿射变换，并且最小化相邻顶点对应仿射变换的差(局部平滑)。由于一对匹配点无法计算得到一个仿射变换，因此引入正则项作为新的约束。
该文章定义了非刚体ICP框架的最优化步骤，扩展的ICP方法适用于非刚体形变并且保留了ICP原本的收敛性。

构造CostFunction

给定template mesh为$\mathcal{S} = \left( \mathcal{V}, \mathcal{E} \right)$拥有$n$个顶点$\mathcal{V}$和$m$条边$\mathcal{E}$的图。Target $\mathcal{T}$为任意一种中3D空间中可以给其中任意一个点找到最近点的表示。
需要求解的数据是每个template顶点对应的仿射变换矩阵$\bm{X}_i$，完整的形式为

$\bm{X} \coloneqq \begin{bmatrix} \bm{X}_1 & \dots & \bm{X}_n \end{bmatrix} ^T$

距离项

形变后template上的顶点和target之间的距离应该很小，距离项由此构造为

$E_d(\bm{X}) \coloneqq \sum_{\bm{v}_i \in \mathcal{V}} w_i \text{dist}^2(\mathcal{T}, \bm{X}_i \bm{v}_i)$

$\bm{v}_i$为其次向量$\begin{bmatrix} x,y,z,1 \end{bmatrix}^T$，点$\bm{v}$和它在target中最近点的距离为$\text{dist}(\mathcal{T}, \bm{v})$。使用分层边界球方法提高最近点搜索的速度。匹配的可信度为权重$w_i$，当某个点没有对应点时，其权重为$0$。

stiffness项

Stiffness项用来惩罚相邻顶点间仿射变换的不同。使用加权矩阵$\bm{G} \coloneqq \text{diag} (1,1,1,\gamma)$，得到

$E_s(\bm{X}) \coloneqq \sum_{\{i,j\} \in \mathcal{E}} \lVert (\bm{X}_i - \bm{X}_j)\bm{G} \rVert ^2_F$

landmark项

用来初始化和引导配准，对cost function本身影响不大。给定一系列landmark $\mathcal{L} = \{ (\bm{v}_{i_1}, \bm{l}_1), \dots , (\bm{v}_{i_l}, \bm{l}_l) \}$将template顶点映射到target顶点上，

$E_l(\bm{X}) \coloneqq \sum_{\{\bm{v}_i, \bm{l}\} \in \mathcal{L}} \lVert \bm{X}_i \bm{v}_i - \bm{l} \rVert ^2$

完整的CostFunction

$E(\bm{X}) \coloneqq E_d(\bm{X}) + \alpha E_s(\bm{X}) + \beta E_l(\bm{X})$

计算步骤

• 初始化$\bm{X}^0$
• 对每个stiffness $\alpha ^i \in \{ \alpha ^1, \dots \alpha ^n \}$都有$\alpha^i > \alpha^{i+1}$
  • 到$\lVert \bm{X}^j - \bm{X}^{j-1} \rVert < \varepsilon$之前
    • 为$\mathcal{V}(\bm{X}^{j-1})$找到出事对应
    • 将$\bm{X}^j$作为初始对应和$\alpha^i$的最优形变

$\mathcal{V}(\bm{X})$为形变后的点。优化步骤包含2个循环。外循环找到越来越优的形变，由非常强的stiff正则化约束，然后逐步降低。内循环在固定的stiffness下寻找形变，初始对应通过最近点搜索寻找，每一次迭代更新template后重新确定对应关系。

优化方法

虽然点匹配会在优化过程中不断更新，但是在每一次循环的时候是固定的。在对应固定的情况下，cost function是一个可以最小化的稀疏二次系统，重写为

$\bar{E}(\bm{X}) \coloneqq \bar{E}_d(\bm{X}) + \alpha E_s(\bm{X}) + \beta E_l(\bm{X})$

距离项

根据重写后的cost function可以看出只有距离项需要重写，因为此时已经确定了具体的点匹配。

$\begin{gathered} \bar{E}_d (\bm{X}) \coloneqq \sum_{\bm{v}_i \in \mathcal{V}} w_i \lVert \bm{X}_i \bm{v}_i - \bm{u}_i \rVert ^2 \\ = \begin{Vmatrix} (\bm{W} \otimes \bm{I}_3) \left(\begin{bmatrix} \bm{X}_1 \\ & \ddots \\ & & \bm{X}_n \end{bmatrix} \begin{bmatrix} \bm{v}_1 \\ \vdots \\ \bm{v}_n \end{bmatrix} - \begin{bmatrix} \bm{u}_1 \\ \vdots \\ \bm{u}_n \end{bmatrix} \right) \end{Vmatrix} ^2 \end{gathered}$

其中$\bm{W} \coloneqq \text{diag}(w_1, \dots, w_n)$，$\otimes$是Kronecker乘积。但是这种形式不容易微分，因此在重新改写一下，

$\begin{aligned} \bar{E}_d(\bm{X}) &= \lVert \bm{W} (\bm{DX} - \bm{U}) \rVert ^2_F \\ \bm{D} & \coloneqq \begin{bmatrix} \bm{v}_1^T \\ & \bm{v}_2^T \\ && \ddots \\ &&& \bm{v}_n^T \end{bmatrix} \\ \bm{U} & \coloneqq \begin{bmatrix} \bm{u}_1 , \dots, \bm{u}_n \end{bmatrix}^T \end{aligned}$

stiffness项

引入node-arc关系矩阵，矩阵由mesh的拓扑图直接定义。每个arc(edge)一行，每个node(vertex)一列。如果边$r$连接顶点$(i,j)$，行$r$中的非$0$元素为$\bm{M}_{ri}=-1$和$\bm{M}_{rj}=1$，有

$E_s(\bm{X}) = \lVert (\bm{M} \otimes \bm{G}) \bm{X} \rVert^2_F$

landmark项

仿照距离项构造$\bm{D}_L$和$\bm{U}_L$

$E_l(\bm{X}) = \lVert \bm{D}_L \bm{X} - \bm{U}_L \rVert ^2_F$

完整的CostFunction

得到一个二次函数

$\begin{aligned} \bar{E}(\bm{X}) &= \begin{Vmatrix} \begin{bmatrix} \alpha \bm{M} \otimes \bm{G} \\ \bm{WD} \\ \beta \bm{D}_L \end{bmatrix} \bm{X} - \begin{bmatrix} \bm{0} \\ \bm{WU} \\ \bm{U}_L \end{bmatrix} \end{Vmatrix} ^2_F \\ &= \lVert \bm{AX} - \bm{B} \rVert ^2_F \end{aligned}$

如果$\bm{A}$的秩小于$4n$，这个问题的约束不够，是病态的问题。该文认为仿射变换的自由度为12。文中有证明$\bm{A}$的列秩为$4n$。

数据缺失和鲁棒性

当template中的顶点$\bm{v}_i$，对应到target中的缺失数据时，权重$w_i$设为$0$。
检测没有对应顶点有3个测试，满足其一则去掉对应$(\bm{X}_i \bm{v}_i, \bm{u}_i)$：

• $\bm{u}_i$在target mesh的边界上
• $\bm{X}_i \bm{v}_i$和$\bm{u}_i$处的mesh法线夹角大于阈值
• 线段$\bm{X}_i \bm{v}_i$到$\bm{u}_i$与形变后的template相交

优化过程中找到的新的对应会使残差变大，以前到优化方法都无法解决这个问题，文中提出的优化过程可以解决这个问题。