具有一定厚度的刚体破碎后会随机形成若干个任意形状的碎块，破碎刚体的复原问题就是将这些碎块按照某种方法进行匹配拼接的过程。两个刚体碎块的匹配就是通过某种方法寻找一个3维刚体变换，使得两个待匹配碎块的曲面的共同部分能够正确地匹配和拼接的过程^[1]。目前，刚体碎块匹配已经在文物复原^[2-4]、医学研究^[5]、对象识别^[6]、遥感图像处理^[7]以及工程计算^[8]等领域得到了较为广泛的应用。

整体来说，碎块的匹配过程可以分为粗匹配和细匹配两个阶段。粗匹配就是将两个碎块的断裂面进行粗略对齐的过程，它能为细匹配带来较好的初值，而细匹配就是指将粗匹配后的两个碎块进一步地对齐，从而完成两个碎块的精确匹配。在众多的碎块匹配算法中，通常采用曲面的特征进行匹配，如曲率^[3]和法向^[9]特征，它们具有旋转和平移不变性，但是离散曲率对噪声敏感，会带来一定的误差。轮廓曲线^[10-11]也是曲面特征描述的一种重要方法，只要正确提取并匹配了刚体碎块断裂面的轮廓曲线，则两个碎块就匹配好了，但是对于轮廓曲线不明显的碎块，该方法并不适合。另外，还可以利用积分不变量^[12]来描述刚体碎块断裂面凹凸不平的几何特征，通过对断裂面上或凸或凹的显著性区域的匹配来完成两个碎块断裂面的匹配。

目前，Besl等人^[13]提出的迭代最近点(ICP)算法也在刚体碎块匹配中得到了较为广泛的应用，该方法通过对两个断裂面上的顶点进行配准来完成两个碎块的断裂面的匹配。ICP算法步骤简单，易于实现，但是在数据量较大时需要的迭代时间比较长，而且对两个待匹配点集的初始位置要求较高，抗噪性也不好。针对ICP算法中存在的这些问题，国内外学者提出了许多改进的ICP算法^[14-19]。比如，Du等人^[20]通过在ICP算法中加入高斯概率模型提出了概率迭代最近点(PICP)算法，该算法不仅在一定程度上提高了配准的精度和速度，而且还具有较强的抗噪性，是一种有效的点云配准算法。Mavridis等人^[21]将模拟退火搜索与标准稀疏ICP算法结合，提出了一种改进的稀疏ICP算法，该算法高效地解决了优化问题，具有较强的抗噪性，能很好地解决两个几何表面的配准。Du等人^[22]提出了尺度迭代最近点算法(SICP)，该算法将一个带有边界的尺度矩阵加入到ICP算法中，解决了不同尺度情况的点云数据的配准问题。虽然这些改进的ICP算法在收敛速度、精度和抗噪性等方面都得到了较大程度的提高，但是对于一些特殊刚体碎块(如秦俑碎块)的匹配效果并不理想。

针对具有厚度且断裂面较原始面明显粗糙的刚体碎块，采用先粗配再细配的方法来完成两个刚体碎块的断裂面匹配。首先对碎块进行曲面分割并提取碎块的断裂面，然后根据积分不变量提取断裂面上的显著性区域并对不同碎块断裂面上的特征区域进行匹配，从而实现碎块的粗匹配。最后，通过加入高斯概率模型、角度约束和动态迭代系数的方式来改进ICP算法，并采用该方法实现两个断裂面的细匹配。

在粗匹配阶段，采用了基于显著性区域的刚体碎块断裂面匹配算法。首先，采用区域生长算法进行曲面分割；然后，根据曲面的粗糙程度提取碎块的断裂面，再根据积分不变量计算断裂面上顶点的凹凸性，进而提取出凹的和凸的显著性区域；最后，通过计算不同碎块的断裂面的显著性区域的相似性对显著性区域进行匹配，从而实现碎块断裂面的粗匹配。

在提取碎块的断裂面之前，首先进行碎块外表面分割，即将碎块的外表面以棱边为界限分割成多张曲面，然后再提取碎块的断裂面。

实验中的碎块数据模型均采用三角网格模型进行表示，这里采用区域生长算法^[23]实现曲面分割，具体方法描述为：首先判断相邻三角面片法矢的夹角是否大于给定阈值，若大于给定阈值则将这两个三角面片分割在两个不同的曲面上，否则分割在同一曲面上，重复该过程直到所有三角面片分割完为止；然后逐一检查分割曲面，若曲面特别小，则将其合并到相邻的曲面中，否则判断该曲面与其相邻曲面的平均法矢夹角，若小于给定阈值则进行合并，否则不合并。

这里主要研究断裂面较原始面粗糙的刚体碎块的匹配问题，因此在曲面分割后可以根据曲面的粗糙程度来提取断裂面。在几何特征上，较为粗糙的断裂面上顶点的法向量变化较原始面的法向量变化大，因此可以根据断裂面上所有顶点的法向量的变化情况来区分断裂面和原始面。

定义曲面$\mathit{\boldsymbol{\varphi }}$上任意顶点$\boldsymbol{p}$的弯曲能量为

$ {e_k}(\boldsymbol{p}) = \frac{1}{k}\sum\limits_{i = 1}^k {\frac{{{{\left\| {{\boldsymbol{n}_p}-{\boldsymbol{n}_{qj}}} \right\|}^2}}}{{{{\left\| {\boldsymbol{p}-{\boldsymbol{q}_j}} \right\|}^2}}}} $

(1)

式中，${{\boldsymbol{q}_j}}$为顶点$\boldsymbol{p}$在曲面$\mathit{\boldsymbol{\varphi }}$上$k$近邻的顶点，${{\boldsymbol{n}_p}}$和${{\boldsymbol{n}_{qj}}}$分别是顶点$\boldsymbol{p}$和${{\boldsymbol{q}_j}}$的法向量，${{{\left\| {\boldsymbol{p}-{\boldsymbol{q}_j}} \right\|}^2}}$表示顶点$\boldsymbol{p}$和${{\boldsymbol{q}_j}}$之间的直线距离。

设置一个以$\boldsymbol{p}$为球心$r$为半径的包围球${B_r}(\boldsymbol{p})$，则点$\boldsymbol{p}$在邻域${B_r}(\boldsymbol{p})$内的局部弯曲能力定义为

$ {{\bar e}_{k, r}}(\boldsymbol{p}) = \frac{1}{{\left| {{N_r}(\boldsymbol{p})} \right|}}\sum\limits_{\boldsymbol{q} \in {\boldsymbol{N}_r}(\boldsymbol{p})} {{e_k}(\boldsymbol{q})} $

(2)

式中，${N_r}(\mathit{\boldsymbol{p}})={B_r}(\mathit{\boldsymbol{p}}) \cap \mathit{\boldsymbol{\varphi }}$。

那么，曲面$\mathit{\boldsymbol{\varphi }}$的平均粗糙度为

$ \rho (\mathit{\boldsymbol{\varphi }}) = \frac{1}{{\left| \mathit{\boldsymbol{\varphi }} \right|}}\sum\limits_{p \in \mathit{\boldsymbol{\varphi }}} {{{\bar e}_{k, r}}(\mathit{\boldsymbol{p}})} $

(3)

式中，${\left| \mathit{\boldsymbol{\varphi }} \right|}$表示曲面$\mathit{\boldsymbol{\varphi }}$上的顶点数。若$\rho (\mathit{\boldsymbol{\varphi }})$大于给定阈值$\rho '$，则认为该曲面为断裂面，否则认为该曲面是原始面。

断裂面上的显著性区域是指一个连通的或凹或凸的局部区域，平面区域不予考虑。对断裂面上顶点的凹凸性的判断可以根据体积积分不变量来计算。

3维曲面$\mathit{\boldsymbol{\varphi }}$上任意一点$\boldsymbol{p}$的体积积分不变量的定义为

$ {V_r}(\boldsymbol{p}) = \int_{{B_r}(\boldsymbol{p})} {{I_D}(x){\rm{d}}x} $

(4)

式中，局部邻域$D$是以$\boldsymbol{p}$为球心$r$为半径的球体，记为${B_r}(\boldsymbol{p})$。${{I_D}(x)}$为曲面示性函数，对于曲面上的任意点$x$，当$x$在曲面$\mathit{\boldsymbol{\varphi }}$外侧时，${{I_D}(x)=1}$，当$x$在曲面$\mathit{\boldsymbol{\varphi }}$内侧时，$\int_{{B_r}(p)} {{I_D}(x)=0} $。3维曲面$\mathit{\boldsymbol{\varphi }}$上任意一点$\boldsymbol{p}$的体积积分不变量${V_r}(\boldsymbol{p})$就是球${B_r}(\boldsymbol{p})$在曲面外侧部分的体积，如图 1所示。

这里，${V_r}(\boldsymbol{p})$的值与曲面在点$\boldsymbol{p}$附近的凹凸程度有关，与$r$的大小成正比，且${V_r}(\boldsymbol{p})$不受${B_r}(\boldsymbol{p})$内部噪声点的影响^[24]。因此，${V_r}(\boldsymbol{p})$反映了点$\boldsymbol{p}$邻域曲面的凹凸程度，可以进行不同半径尺度的计算，且对噪声鲁棒。

为了方便描述，定义${S_r}(\boldsymbol{p})$来描述顶点$\boldsymbol{p}$在半径$r$区域内的凹凸性，形式为

$ {S_r}(\boldsymbol{p}) = \frac{3}{{4{\rm{\pi }}{r^3}}}({V_r}(\boldsymbol{p})-\frac{2}{3}{\rm{\pi }}{r^3}) $

(5)

于是有，当$\boldsymbol{p}$为凸顶点时，${V_r}(\boldsymbol{p}) > \frac{2}{3}{\rm{\pi }}{r^3}, {S_r}(\boldsymbol{p}) > 0$；当$\boldsymbol{p}$为凹顶点时，${V_r}(\boldsymbol{p}) < \frac{2}{3}{\rm{\pi }}{r^3}, {S_r}(\boldsymbol{p}) < 0$；当$\boldsymbol{p}$为平面顶点时，${V_r}(\boldsymbol{p})=\frac{2}{3}{\rm{\pi }}{r^3}, {S_r}(\boldsymbol{p})=0$。那么，${S_r}(\boldsymbol{p})$的取值范围为${S_r}(\boldsymbol{p}) \in [-\frac{1}{2}, \frac{1}{2}]$。在实际计算中需要设置一个阈值${\varepsilon _r}$，当${\boldsymbol{S}_r}(\boldsymbol{p}) < - {\varepsilon _r}$时，$\boldsymbol{p}$为凹顶点；当${\boldsymbol{S}_r}(\boldsymbol{p}) > {\varepsilon _r}$时，$\boldsymbol{p}$为凸顶点；当$-{\varepsilon _r} \le {S_r}(\boldsymbol{p}) \le {\varepsilon _r}$时，$\boldsymbol{p}$为平面顶点。

根据以上顶点凹凸性的判断方法，显著性区域的提取算法描述如下：

1)计算断裂面上每个顶点的体积积分不变量${V_r}(\boldsymbol{p})$和${S_r}(\boldsymbol{p})$；

2)根据${S_r}(\boldsymbol{p})$对断裂面上的顶点进行分类，将满足${S_r}(\boldsymbol{p}) > {\varepsilon _r}$的凸顶点记入集合${\boldsymbol{S}_t}$中，满足${S_r}(\boldsymbol{p}) < -{\varepsilon _r}$的凹顶点记入集合${\boldsymbol{S}_a}$中；

3)对集合${\boldsymbol{S}_t}$和${\boldsymbol{S}_a}$中的顶点，采用聚类算法将断裂面划分为多个或凹或凸的显著性区域。若某个显著性区域顶点数目少于给定阈值，可以删掉该区域。

基于显著性区域的粗匹配就是对不同碎块的断裂面上的显著性区域进行匹配。设有两个显著性区域${\boldsymbol{R}_1}$和${\boldsymbol{R}_2}$，并且${\boldsymbol{R}_1}$和${\boldsymbol{R}_2}$分别属于两个不同碎块的断裂面${\mathit{\boldsymbol{\varphi }}_1}$和${\mathit{\boldsymbol{\varphi }}_2}$，若区域${\boldsymbol{R}_1}$和${\boldsymbol{R}_2}$同时满足：类型不同，即一个是凹区域且另一个是凸区域；具有相近的主成分；面积大小接近，则认为它们是相似区域。这里采用主成分分析法(PCA)来比较显著性区域的相似性。

设显著性区域$\boldsymbol{R}$包含$n$个顶点$\{ {\boldsymbol{p}_1}, {\boldsymbol{p}_2}, \cdots, {\boldsymbol{p}_n}\} $，顶点的坐标为$({x_i}, {y_i}, {z_i})$，则$\boldsymbol{R}$的质心为

$ \bar p = \frac{1}{n}\sum\limits_{i = 1}^n {{\boldsymbol{p}_i}} $

(6)

$\boldsymbol{R}$的3×3的协方差矩阵为

$ \boldsymbol{M} = \sum\limits_{i = 1}^n [ ({\boldsymbol{p}_i} - \boldsymbol{\bar p}){({\boldsymbol{p}_i} - \boldsymbol{\bar p})^{\rm{T}}}] $

(7)

假设计算得到的$\boldsymbol{M}$特征值为$\lambda _0^R \ge \lambda _1^R \ge \lambda _2^R$，则将其定义为区域$\boldsymbol{R}$的主成分，其对应的特征向量为$\boldsymbol{n}_0^R, \boldsymbol{n}_1^R, \boldsymbol{n}_2^R$，将其定义为区域$\boldsymbol{R}$的主方向。

那么，定义各向异性特征为

$ A(\boldsymbol{R}) = {\left| {\frac{{\lambda _1^R}}{{\lambda _0^R}}} \right|^{\frac{1}{2}}} $

(8)

对于两个显著性区域${\boldsymbol{R}_1}$和${\boldsymbol{R}_2}$，且${\boldsymbol{R}_1}$和${\boldsymbol{R}_2}$分别属于两个不同碎块的断裂面${\mathit{\boldsymbol{\varphi }}_1}$和${\mathit{\boldsymbol{\varphi }}_2}$，定义

$ SA({\boldsymbol{R}_1}, {\boldsymbol{R}_2}) = \left| {\frac{{A({\boldsymbol{R}_1})-A({\boldsymbol{R}_2})}}{{A({\boldsymbol{R}_1}) + A({\boldsymbol{R}_2})}}} \right| $

(9)

若$SA ({\boldsymbol{R}_1}, {\boldsymbol{R}_2}) < {\varepsilon _a}$，则认为${\boldsymbol{R}_1}$和${\boldsymbol{R}_2}$具有相近的主成分，${\varepsilon _a}$为提前设定的阈值。

设${\boldsymbol{R}_1}$和${\boldsymbol{R}_2}$的面积分别为$N ({\boldsymbol{R}_1})$和$N ({\boldsymbol{R}_2})$，则定义${\varepsilon _a}$为

$ SN({\boldsymbol{R}_1}, {\boldsymbol{R}_2}) = \left| {\frac{{N({\boldsymbol{R}_1})-N({\boldsymbol{R}_2})}}{{N({\boldsymbol{R}_1}) + N({\boldsymbol{R}_2})}}} \right| $

(10)

若$SN ({\boldsymbol{R}_1}, {\boldsymbol{R}_2}) < {\varepsilon _n}$，则认为${\boldsymbol{R}_1}$和${\boldsymbol{R}_2}$的面积大小相近，${\varepsilon _n}$为提前设定的阈值。

根据以上对显著性区域的特征的计算结果，再采用穷举法寻找最优匹配，对每一个可能的匹配计算其3维变换，即可完成刚体碎块断裂面的粗匹配。

在细匹配阶段，采用改进的ICP算法来实现刚体碎块断裂面上顶点的进一步匹配。这里通过在ICP算法中引入高斯概率模型、角度约束和动态迭代系数来对其改进，以进一步提高碎块匹配的速度、精度和抗噪性等，达到两个刚体碎块的精确匹配。

给定两个待配准的点集$M=\{ {m_i}\} _{i=1}^{{N_M}}$和$\boldsymbol{D}=\{ {\boldsymbol{d}_j}\} _{j=1}^{{N_D}}$，${N_M}$和${N_D}$分别表示点集$\boldsymbol{M}$和$\boldsymbol{D}$的点数，若它们之间满足3维刚体变换，则它们之间的距离可描述为

$ \begin{array}{l} \mathop {\min }\limits_{\boldsymbol{R}, \boldsymbol{t}, j \in \{ 1, 2, \cdots, {N_D}\} } (\sum\limits_{i = 1}^{{N_M}} {\left\| {(\boldsymbol{R}{\boldsymbol{m}_i} + \boldsymbol{t})-{\boldsymbol{d}_j}} \right\|_2^2} )\\ \;\;\;\;{\rm{s}}.\;{\rm{t}}.\;\;\;\;{\boldsymbol{R}^{\rm{T}}}\boldsymbol{R} = {\boldsymbol{I}_n}, \det (\boldsymbol{R}) = 1 \end{array} $

(11)

式中，$\boldsymbol{R}$为旋转矩阵，$\boldsymbol{t}$为平移矩阵。

在每次迭代过程中，ICP算法都要在$\boldsymbol{D}$中通过寻找距离$\boldsymbol{M}$最近的点来建立它们之间的相关性，以实现刚体变换。那么，ICP算法的基本步骤如下：

1)根据第$k$-1步的已知刚体变换${\boldsymbol{R}_{k-1}}$和${\boldsymbol{t}_{k-1}}$，首先将点集$\boldsymbol{M}$进行${\boldsymbol{R}_{k-1}}{\boldsymbol{m}_i} + {\boldsymbol{t}_{k-1}}$变换，然后再建立两个点集之间的相关性${c_k}(i)$，其数学描述为

$ {c_k}(i) = \mathop {\arg \;\min }\limits_{j \in \{ 1, 2, \cdots, {N_D}\} } \left\| {({\boldsymbol{R}_{k-1}}{\boldsymbol{m}_i} + {\boldsymbol{t}_{k-1}})-{\boldsymbol{d}_j}} \right\|_2^2 $

(12)

式中，$i=1, 2, \cdots, {N_M}$。

2)计算点集$\boldsymbol{M}$和$\boldsymbol{D}$的刚体变换，其数学描述为

$ ({\boldsymbol{R}_k}, {\boldsymbol{t}_k}) = \mathop {\arg \;\min }\limits_{\boldsymbol{R}_k^{\rm{T}}{\boldsymbol{R}_k} = {\boldsymbol{I}_n}, \det (\boldsymbol{R}) = 1} (\sum\limits_{i = 1}^{{N_M}} {\left\| {{\boldsymbol{R}_k}{\boldsymbol{m}_i} + {\boldsymbol{t}_k}-{\boldsymbol{d}_{{c_k}(i)}}} \right\|_2^2} ) $

(13)

重复步骤1)2)直到达到迭代终止条件为止。由于ICP算法是一个局部收敛算法，因此旋转矩阵和平移矩阵的初值非常重要。好的初值不仅能保证算法收敛到全局最小值，而且能大大提高算法的执行效率，初值的选取可参考文献[25-26]的方法。

ICP算法的改进主要分3个步骤进行，首先在ICP算法中加入高斯概率模型^[21]以提高算法的抗噪性，然后加入角度约束来解决因旋转角的变化过大而引起的配准效果不佳的问题，最后加入动态迭代系数来提高算法的收敛速度。

若点集$\boldsymbol{M}=\{ {\boldsymbol{m}_i}\} _{i=1}^{{N_M}}$和$\; \boldsymbol{D}=\{ {\boldsymbol{d}_j}\} _{j=1}^{{N_D}}$含噪声，则可将这两个点集的配准问题转化为概率密度估计问题。根据全概率定理，点集$\boldsymbol{M}$的概率公式可表示为

$ p(\boldsymbol{m}) = \sum\limits_{i = 1}^{{N_M}} {p({\boldsymbol{d}_j})p(\boldsymbol{m}\left| {{\boldsymbol{d}_j}} \right.)} $

(14)

如果点集$\boldsymbol{M}$和$\boldsymbol{D}$都满足高斯概率模型，那么它们之间的关系可描述为

$ p(\boldsymbol{m}\left| \boldsymbol{d} \right.) = \frac{1}{{{{(2{\rm{\pi }}{\sigma ^2})}^{n/2}}}}{\exp ^{-\frac{{\left\| {\boldsymbol{T}(\boldsymbol{m})-\boldsymbol{d}} \right\|_2^2}}{{2{\sigma ^2}}}}} $

(15)

式中，$n$表示点集的维数，${{\sigma ^2}}$表示高斯概率模型的方差，$\boldsymbol{T}(\boldsymbol{m})$表示$\boldsymbol{m}$的刚体变换，$\left| {\left. {\boldsymbol{T}(\boldsymbol{m})-\boldsymbol{d}} \right\|_2^2} \right.$表示点集$\boldsymbol{T}(\boldsymbol{m})$和$\boldsymbol{d}$之间的距离。

一般来说，对于3维扫描仪器获取的初始数据，都是杂乱的3维点云数据，需要对其进行模型简化，使其在数据量上适当降低并满足均匀分布。因此对于点集$\boldsymbol{D}$，有

$ p(\boldsymbol{d}) = \frac{1}{{{N_D}}} $

(16)

那么，似然函数的目标函数可定义为

$ \begin{array}{l} F(\boldsymbol{m};\boldsymbol{R}, \boldsymbol{t}, {\sigma ^2}) =-\sum\limits_{i = 1}^{{N_M}} {{\rm{I}}{{\rm{n}}_e}} (p({\boldsymbol{m}_i}\left| {\boldsymbol{R}, \boldsymbol{t}, {\sigma ^2}} \right.)) = \\ \;\;\;-\sum\limits_{i = 1}^{{N_M}} {{\rm{I}}{{\rm{n}}_e}} (\sum\limits_{i = 1}^{{N_D}} {p({\boldsymbol{m}_j})p({\boldsymbol{m}_i}\left| {{\boldsymbol{d}_j}, } \right.\boldsymbol{R}, \boldsymbol{t}, {\sigma ^2})} ) \end{array} $

(17)

首先，建立求解刚体变换的新目标函数为

$ \begin{array}{l} M({{\mathit{\boldsymbol{\tilde R}}}_k}, {{\mathit{\boldsymbol{\tilde t}}}_k}, {{\mathit{\boldsymbol{\tilde \sigma }}}_k}) = \sum\limits_{i = 1}^{{N_M}} {\sum\limits_{j = 1}^{{N_D}} {(p({\mathit{\boldsymbol{d}}_j}\left| {{\mathit{\boldsymbol{m}}_i}, {\mathit{\boldsymbol{R}}_{k - 1}}, {\mathit{\boldsymbol{t}}_{k - 1}}, \sigma _{k - 1}^2} \right.)} } \\ \;\;\;\;\;\;\;\;\;\;\;\;{\rm{log(}}p({\mathit{\boldsymbol{d}}_j})p({\mathit{\boldsymbol{m}}_i}\left| {{\mathit{\boldsymbol{d}}_j}, {{\mathit{\boldsymbol{\tilde R}}}_k}, {{\mathit{\boldsymbol{\tilde t}}}_k}, \tilde \sigma _k^2} \right.){\rm{)}}) \end{array} $

(18)

令$p_i^{k -1}=p ({\mathit{\boldsymbol{d}}_{{c_k}(i)}}\left| {{\mathit{\boldsymbol{m}}_i}\left| {{\mathit{\boldsymbol{d}}_j}, {{\mathit{\boldsymbol{\tilde R}}}_k}, {{\mathit{\boldsymbol{\tilde t}}}_k}, \tilde \sigma _k^2} \right.)} \right.$。由于ICP算法中两个点集的相关性是一对一的，因此式(18)可简化为

$ \begin{array}{l} \mathit{\boldsymbol{M}}({{\mathit{\boldsymbol{\tilde R}}}_k}, {{\mathit{\boldsymbol{\tilde t}}}_k}, \mathit{\tilde \sigma }_k^2) = - \sum\limits_{i = 1}^{{N_M}} {\mathit{\boldsymbol{p}}_i^{k - 1}} \times \frac{{\left\| {{{\mathit{\boldsymbol{\tilde R}}}_k}{\mathit{\boldsymbol{m}}_i} + {{\mathit{\boldsymbol{\tilde t}}}_k} - {\mathit{\boldsymbol{d}}_{{c_k}(i)}}} \right\|_2^2}}{{2{{\mathit{\tilde \sigma }}^2}}} - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_{i = 1}^{{N_M}} {\mathit{\boldsymbol{p}}_i^{k - 1}} \times \log ({N_D}{(2{\rm{ \mathsf{ π} }}\tilde \sigma _k^2)^{\frac{n}{2}}}) \end{array} $

(19)

式中，只有$ -\sum\limits_{i=1}^{{N_M}} {\mathit{\boldsymbol{p}}_i^{k -1}} \times \frac{{\left\| {{{\mathit{\boldsymbol{\tilde R}}}_k}{\mathit{\boldsymbol{m}}_i} + {{\mathit{\boldsymbol{\tilde t}}}_k} -{\mathit{\boldsymbol{d}}_{{c_k}(i)}}} \right\|_2^2}}{{2{{\mathit{\tilde \sigma }}^2}}}$与旋转矩阵和平移矢量有关，因此最大化式(19)的问题就等价于关于旋转矩阵和平移向量最大化式(19)的第1部分的问题，因此求解目标进一步转换为

$ \begin{array}{l} ({\mathit{\boldsymbol{R}}_k}, {\mathit{\boldsymbol{t}}_k}) = \mathop {\arg \;\min }\limits_{\mathit{\boldsymbol{\tilde R}}_k^{\rm{T}}{{\mathit{\boldsymbol{\tilde R}}}_k} = {\mathit{\boldsymbol{I}}_n}, \det (\mathit{\boldsymbol{\tilde R}}) = 1, {{\mathit{\boldsymbol{\tilde t}}}_k}} \sum\limits_{i = 1}^{{N_M}} {\mathit{\boldsymbol{p}}_i^{k - 1}} \times \\ \;\;\;\;\;\;\;\;\;\;\frac{{\left\| {{{\mathit{\boldsymbol{\tilde R}}}_k}{\mathit{\boldsymbol{m}}_i} + {{\mathit{\boldsymbol{\tilde t}}}_k} - {\mathit{\boldsymbol{d}}_{{c_k}(i)}}} \right\|_2^2}}{{2{{\mathit{\tilde \sigma }}^2}}} \end{array} $

(20)

角度约束是指为刚体变换中的旋转角设置边界(即上界和下界)，其目的是为了降低角度变化过大引起的配准失败问题。

首先，定义旋转矩阵$\boldsymbol{R}={\boldsymbol{R}_x}{\boldsymbol{R}_y}{\boldsymbol{R}_z}$。${\boldsymbol{R}_x}, {\boldsymbol{R}_y}, {\boldsymbol{R}_z}$分别为

$ \begin{array}{l} {\boldsymbol{R}_x} = \left[{\begin{array}{*{20}{c}} 1&0&0\\ 0&{\cos {\theta _x}}&{-\sin {\theta _x}}\\ 0&{\sin {\theta _x}}&{\cos {\theta _x}} \end{array}} \right]\\ {\boldsymbol{R}_y} = \left[{\begin{array}{*{20}{c}} {\cos {\theta _y}}&0&{\sin {\theta _y}}\\ 0&1&0\\ {-\sin {\theta _y}}&0&{\cos {\theta _y}} \end{array}} \right]\\ {\boldsymbol{R}_z} = \left[{\begin{array}{*{20}{c}} {\cos {\theta _z}}&{-\sin {\theta _z}}&0\\ {\sin {\theta _z}}&{\cos {\theta _z}}&0\\ 0&0&1 \end{array}} \right] \end{array} $

${\theta _x}$, ${\theta _y}$, ${\theta _z}$是旋转角^[27]，${\theta _xb}$, ${\theta _yb}$, ${\theta _zb}$是旋转角的均值，$\Delta {\theta _x}, \Delta {\theta _y}, \Delta {\theta _z}$是旋转角的偏差，${\theta _{xb}}-\Delta {\theta _x}, {\theta _{yb}}-\Delta {\theta _y}, {\theta _{zb}}-\Delta {\theta _z}$是旋转角的下界，${\theta _{xb}} + \Delta {\theta _x}, {\theta _{yb}} + \Delta {\theta _y}$和${\theta _{zb}} + \Delta {\theta _z}$是旋转角的上界。

这里，采用PCA方法来计算旋转角的主成分矢量。对于3维点集的刚体配准，如果这两个点集能够正确配准，那它们必然具有相似的点分布。因此当它们正确配准时，旋转角的3对主成分矢量都接近0°或180°。所以旋转角可以通过旋转角的3对主成分矢量来估计。

假设点集$\boldsymbol{M}$和$\boldsymbol{D}$的主成分分别为：${\boldsymbol{V}_M}=\{ {\boldsymbol{V}_{M1}}, {\boldsymbol{V}_{M2}}, {\boldsymbol{V}_{M3}}\} $和${\boldsymbol{V}_D}=\{ {\boldsymbol{V}_{D1}}, {\boldsymbol{V}_{D2}}, {\boldsymbol{V}_{D3}}\} $，${\boldsymbol{V}_{Mi}}$和${\boldsymbol{V}_{Di}}$是3×1的矢量。那么点集$\boldsymbol{D}$中有8组主成分满足条件

$ \begin{array}{l} \boldsymbol{V}_D^1 = \{ {\boldsymbol{V}_{D1}}, {\boldsymbol{V}_{D2}}, {\boldsymbol{V}_{D3}}\} \\ \boldsymbol{V}_D^2 = \{ {\boldsymbol{V}_{D1}}, {\boldsymbol{V}_{D2}}, -{\boldsymbol{V}_{D3}}\} \\ \boldsymbol{V}_D^3 = \{ {\boldsymbol{V}_{D1}}, -{\boldsymbol{V}_{D2}}, {\boldsymbol{V}_{D3}}\} \\ \boldsymbol{V}_D^4 = \{ {\boldsymbol{V}_{D1}}, -{\boldsymbol{V}_{D2}}, - {\boldsymbol{V}_{D3}}\} \\ \boldsymbol{V}_D^5 = \{ - {\boldsymbol{V}_{D1}}, {\boldsymbol{V}_{D2}}, {\boldsymbol{V}_{D3}}\} \\ \boldsymbol{V}_D^6 = \{ - {\boldsymbol{V}_{D1}}, {\boldsymbol{V}_{D2}}, - {\boldsymbol{V}_{D3}}\} \\ \boldsymbol{V}_D^7 = \{ - {\boldsymbol{V}_{D1}}, - {\boldsymbol{V}_{D2}}, {\boldsymbol{V}_{D3}}\} \\ \boldsymbol{V}_D^8 = \{ - {\boldsymbol{V}_{D1}}, - {\boldsymbol{V}_{D2}}, - {\boldsymbol{V}_{D3}}\} \end{array} $

为了使点集$\boldsymbol{M}$和$\boldsymbol{D}$配准成功，旋转矩阵为${\boldsymbol{R}^j}$必须满足：${\mathit{\boldsymbol{R}}^j}{\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{M}}}=\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{D}}^j$，所以${\mathit{\boldsymbol{R}}^j}\mathit{\boldsymbol{=V}}_\mathit{\boldsymbol{D}}^j\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{M}}^{ -1}, j \in \{ 1, 2, \cdots, 8\}$。旋转矩阵${\boldsymbol{R}^j}$中能够使点集$\boldsymbol{M}$和$\boldsymbol{D}$最佳配准的旋转矩阵${\boldsymbol{R}_{opt}}$即为最佳旋转矩阵，其具体计算方法为

$ {\boldsymbol{R}_{opt}} = \left[ {\begin{array}{*{20}{l}} {\;\;\;\cos ({\theta _{yb}})\cos ({\theta _{zb}})}&{ - \cos ({\theta _{yb}})\sin ({\theta _{zb}})}&{\;\;\;\sin ({\theta _{yb}})}\\ {\begin{array}{*{20}{l}} {\;\;\;\sin ({\theta _{xb}}){\rm{sin}}({\theta _{yb}})\cos ({\theta _{zb}}) + }\\ {\;\;\;{\rm{cos}}({\theta _{xb}})\sin ({\theta _{zb}})} \end{array}}&{\begin{array}{*{20}{l}} { - \sin ({\theta _{xb}})\sin ({\theta _{yb}})\sin ({\theta _{zb}}) + }\\ {\;\;\;\cos ({\theta _{xb}})\cos ({\theta _{zb}})} \end{array}}&{ - \sin ({\theta _{xb}})\cos ({\theta _{yb}})}\\ {\begin{array}{*{20}{l}} { - \cos ({\theta _{xb}})\sin ({\theta _{yb}})\cos ({\theta _{zb}}) + }\\ {\;\;\;\sin ({\theta _{xb}})\sin ({\theta _{zb}})} \end{array}}&{\begin{array}{*{20}{l}} {\;\;\;\cos ({\theta _{xb}})\sin ({\theta _{yb}})\sin ({\theta _{zb}}) + }\\ {\;\;\;\sin ({\theta _{xb}})\sin ({\theta _{zb}})} \end{array}}&{\;\;\;\cos ({\theta _{xb}})\cos ({\theta _{yb}})} \end{array}} \right] $

(21)

由此，旋转角的均值${\theta _{xb}}, {\theta _{yb}}, {\theta _{zb}}$都可以通过式(21)计算得到。旋转角的偏差$\Delta {\theta _x}, \Delta {\theta _y}, \Delta {\theta _z}$的值取决于点集的相似程度，具体由实验情况来给出，为了计算方便可以设置为$\Delta {\theta _x}=\Delta {\theta _y}=\Delta {\theta _z}=10^\circ $。

那么，点集$\boldsymbol{M}$和$\boldsymbol{D}$的配准问题可进一步描述为

$ \left\{ \begin{array}{l} ({\mathit{\boldsymbol{R}}_k}, {\mathit{\boldsymbol{t}}_k}) = \mathop {\arg \;\min }\limits_{\mathit{\boldsymbol{\tilde R}}_k^{\rm{T}}{{\mathit{\boldsymbol{\tilde R}}}_k} = {\mathit{\boldsymbol{I}}_n}, \det (\mathit{\boldsymbol{\tilde R}}) = 1, {{\mathit{\boldsymbol{\tilde t}}}_k}} \sum\limits_{i = 1}^{{N_M}} {\mathit{\boldsymbol{p}}_i^{k - 1}} \times \\ \;\;\;\;\;\;\;\;\;\;\frac{{\left\| {{{\mathit{\boldsymbol{\tilde R}}}_k}{\mathit{\boldsymbol{m}}_i} + {{\mathit{\boldsymbol{\tilde t}}}_k} - {\mathit{\boldsymbol{d}}_{{c_k}(i)}}} \right\|_2^2}}{{2{{\mathit{\tilde \sigma }}^2}}}\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;{\mathit{\boldsymbol{R}}^{\rm{T}}}\mathit{\boldsymbol{R}} = {\mathit{\boldsymbol{I}}_n}, \det (\mathit{\boldsymbol{R}}) = 1\\ \;\;\;\;\;\;\;{\theta _x} \in [{\theta _{xb}}-\Delta {\theta _x}, {\theta _{xb}} + \Delta {\theta _x}]\\ \;\;\;\;\;\;\;{\theta _y} \in [{\theta _{yb}}-\Delta {\theta _y}, {\theta _{yb}} + \Delta {\theta _y}]\\ \;\;\;\;\;\;\;{\theta _z} \in [{\theta _{zb}}-\Delta {\theta _z}, {\theta _{zb}} + \Delta {\theta _z}] \end{array} \right. $

(22)

式(22)是一个不等式约束的优化问题，可以采用基于定区间内目标函数单调性的封闭算法^[26]来解决。

定义均方根(RMS)、方差为

$ {\rm{RMS = (}}\sum\limits_{i = 1}^{{N_M}} {{\mathit{\boldsymbol{p}}_i}} \times \left\| {{\mathit{\boldsymbol{R}}_k}{\mathit{\boldsymbol{m}}_i} + {\mathit{\boldsymbol{t}}_k} - {\mathit{\boldsymbol{d}}_{{c_k}(i)}}} \right\|_2^2{{\rm{)}}^{\frac{1}{2}}} $

(23)

$ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _k^2 = \sum\limits_{i = 1}^{{N_M}} {\mathit{\boldsymbol{p}}_i^{k - 1}} \times \frac{{\left\| {{\mathit{\boldsymbol{R}}_k}{\mathit{\boldsymbol{m}}_i} + {\mathit{\boldsymbol{t}}_k} - {\mathit{\boldsymbol{d}}_{{c_k}(i)}}} \right\|_2^2}}{n} $

(24)

那么，高斯概率方差的更新为

$ \sigma _k^2 = \left\{ \begin{array}{l} \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _{k-1}^2/\lambda \;\;\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _{k-1}^2/\lambda > \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _k^2\\ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _k^2\;\;\;\;\;\;\;\;\;\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _{k-1}^2/\lambda < \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _k^2 \end{array} \right. $

(25)

式中，$\lambda $称作退火系数。当$\lambda $=1时，方差不会更新，就是所谓的ICP算法。当$\lambda $取值较大时，就可以有效地提高迭代速度，但是容易陷入局部最小值。这里，一般$\lambda $的取值在1到2之间，其取值由具体的实验来决定。

因此，后验高斯概率为

$ \begin{array}{l} p({\mathit{\boldsymbol{d}}_{{c_k}(i)}}\left| {{\mathit{\boldsymbol{m}}_i}\left| {{\mathit{\boldsymbol{R}}_k}, {\mathit{\boldsymbol{t}}_k}, \sigma _k^2} \right.} \right.) = \\ \frac{{1/{{(2{\rm{ \mathsf{ π} }}\sigma _k^2)}^{\frac{n}{2}}}{{\exp }^{ - \frac{{\left\| {{\mathit{\boldsymbol{R}}_k}{\mathit{\boldsymbol{m}}_i} + {\mathit{\boldsymbol{t}}_k} - {\mathit{\boldsymbol{d}}_{{c_k}(i)}}} \right\|_2^2}}{{2\sigma _k^2}}}}}}{{\sum\limits_{i = 1}^{{N_M}} 1 /{{(2{\rm{ \mathsf{ π} }}\sigma _k^2)}^{\frac{n}{2}}}{{\exp }^{ - \frac{{\left\| {{\mathit{\boldsymbol{R}}_k}{\mathit{\boldsymbol{m}}_i} + {\mathit{\boldsymbol{t}}_k} - {\mathit{\boldsymbol{d}}_{{c_k}(i)}}} \right\|_2^2}}{{2\sigma _k^2}}}}}} \end{array} $

(26)

动态迭代系数是一个可以自动调整刚体变换参数的整数值，它可以在不影响ICP算法的配准精度和收敛方向的情况下，大大减少算法的迭代次数、提高算法的收敛速度。

在ICP算法中加入动态迭代系数$h$的步骤如下：

1)计算刚体变换矢量$\boldsymbol{q}={[\boldsymbol{R}\left| \boldsymbol{t} \right.]^{\rm{T}}}$以及${\boldsymbol{q}_{k-1}}$的相邻两次迭代的变化量$\Delta {\boldsymbol{q}_{k-1}}$；

2)用刚体变换矢量$\Delta {\boldsymbol{q}_{k-1}}$更新基本ICP算法中的${\boldsymbol{M}_k}$共$h$次，即执行${\boldsymbol{M}_k}=\; \Delta {\boldsymbol{q}_{k-1}}{\boldsymbol{M}_{k-1}}$共$h$次。ICP算法通过${\boldsymbol{M}_k}=\; {\boldsymbol{q}_{k-1}}{\boldsymbol{M}_0}$来计算新的点集，而加入动态迭代系数$h$后的ICP算法则通过${\boldsymbol{M}_k}=\; \Delta {\boldsymbol{q}_{k-1}}{\boldsymbol{M}_{k-1}}$来计算新的点集。

通常，动态迭代系数$h$取大于等于0的整数。当$h$不同时，算法的收敛速度也不同，如图 2所示^[28]。随着动态迭代系数$h$的增大，收敛速度也会越来越快，但是当$h$增大到一定的程度，收敛曲线会出现震荡，可能不再收敛。

改进的ICP算法可描述如下：

1)给定刚体变换初值$\boldsymbol{q}={[{\boldsymbol{R}_0}, {\boldsymbol{t}_0}]^{\rm{T}}}$，${\boldsymbol{R}_0}$为初始旋转矩阵，${\boldsymbol{t}_0}$为初始平移矢量，令${\boldsymbol{M}_0}={\boldsymbol{R}_0}\boldsymbol{M} + {\boldsymbol{t}_0}$，退火系数$\lambda \in (1, 2]$，初始概率$p_i^0=\frac{1}{{{N_x}}}, i=1, 2, \cdots, {N_x}$，迭代次数$k$=0，动态迭代系数$h$=0；

2)估计旋转角度${\theta _x}, {\theta _y}, {\theta _z}$的边界，即${\theta _x} \in [{\theta _{xb}}-\Delta {\theta _x}, {\theta _{xb}} + \Delta {\theta _x}]$, ${\theta _y} \in [{\theta _{yb}}-\Delta {\theta _y}, {\theta _{yb}} + \Delta {\theta _y}]$和${\theta _z} \in [{\theta _{zb}}-\Delta {\theta _z}, {\theta _{zb}} + \Delta {\theta _z}]$；

3)令$k$=$k$+1，并建立相关性${c_k}(i) \in \{ 1, 2, \cdots, {N_D}\} $；

4)利用奇异值分解(SVD)的方法来计算旋转矩阵${\boldsymbol{R}_k}$和平移矢量${\boldsymbol{t}_k}$；

5)计算最优旋转角${\theta _{xo}}, {\theta _{yo}}, {\theta _{zo}}$，最优旋转矩阵${\boldsymbol{R}_{{\rm{opt}}}}={\boldsymbol{R}_{xo}}{\boldsymbol{R}_{yo}}{\boldsymbol{R}_{zo}}$，平移矢量${\mathit{\boldsymbol{t}}_k}=\frac{1}{{{N_M}}}(-{\mathit{\boldsymbol{R}}_k}\sum\limits_{i=1}^{{N_M}} {{m_i}} + \sum\limits_{i=1}^{{N_M}} {{\mathit{\boldsymbol{d}}_{{c_k}(i)}}})$，以及${\boldsymbol{q}_k}$的相邻两次迭代的变化量$\Delta {\boldsymbol{q}_k}$；

6)计算${\boldsymbol{M}_k}={\boldsymbol{R}_k}\boldsymbol{M} + {\boldsymbol{t}_k}$；

7)判断均方根误差RMS，若${\rm{RM}}{{\rm{S}}_k}-{\rm{RM}}{{\rm{S}}_{k-1}} > 0$，则执行$h$=$h$+1操作，否则执行$h$=0操作；

8)判断动态迭代系数$h$，若$h$>0，则通过执行${\boldsymbol{M}_k}=\Delta {\boldsymbol{q}_k}({\boldsymbol{M}_k})\; $共$h$次来更新点集${\boldsymbol{M}_k}$；

9)计算两组点集的方差，并更新高斯概率的方差；

10)计算后验高斯概率，直到满足$\left| {{\rm{RM}}{{\rm{S}}_k}-{\rm{RM}}{{\rm{S}}_{k-1}}} \right| < \varepsilon $或$k > {\rm{Ste}}{{\rm{p}}_{\max }}$算法才终止，否则转到步骤3)，这里$\varepsilon $和${\rm{Ste}}{{\rm{p}}_{\max }}$是预先设置的阈值。

实验中的所有数据均为采用3维扫描仪获取的刚体碎块数据，采用三角网格模型表示，且待匹配碎块满足简单的3维刚体变换。这里进行了两组实验来完成不同类型的碎块匹配。实验1采用的是石头和灰泥碎块数据模型^[29]，是网络上的公共数据模型；实验2采用的是在秦始皇陵扫描的秦俑碎块数据模型，该模型含有较大的噪声，数据比较复杂，匹配难度较大。具体的匹配过程如下：首先采用第1节介绍的基于显著性区域的粗匹配算法将碎块进行粗略对齐，然后采用第2节提出的改进的ICP算法将碎块的断裂面进一步精确匹配，从而完成碎块的最终匹配。

采用了4组石头碎块和灰泥碎块用于匹配实验。首先采用基于显著性区域的粗匹配算法将碎块进行粗略对齐，然后再用改进的ICP算法将碎块的断裂面进行细匹配。实验中选取的旋转角偏差为$\Delta {\theta _x}=\Delta {\theta _y}=\Delta {\theta _z}=10^\circ $，动态迭代系数$h$=4，其匹配结果如图 3所示。

从图 3的匹配结果来看，改进的ICP算法能很好地完成碎块的匹配。第1组的石头碎块匹配属于完全匹配类型，第2组的石头碎块匹配属于部分匹配类型，第3组的灰泥碎块匹配属于完全匹配类型，第4组的灰泥碎块匹配属于部分匹配类型。为了进一步验证改进ICP算法的性能，在基于显著性区域粗匹配算法的基础上，再分别采用ICP算法和PICP算法^[21]进行细匹配。ICP算法、PICP算法和改进的ICP算法等3种算法的匹配点数、RMS、迭代次数以及匹配耗时等匹配参数如表 1所示。

从表 1的匹配结果来看，ICP算法、PICP算法以及改进的ICP算法均能完成碎块的细匹配，但是PICP算法的匹配精度和速度明显优于ICP算法，而改进的ICP算法的匹配精度和速度又明显优于PICP算法，特别是在匹配速度方面的提高尤其明显。因此，改进的ICP算法是一种高精度、快速的匹配算法。

两组实验数据采用3D扫描仪获取的秦俑碎块数据，采用3维网格模型表示，如图 4所示。

图 4中第1组的实验数据是秦俑的两个腿部碎块数据模型(图 4(a))，这两个碎块相对保存完整，没有明显的损坏或缺失。图 4中第2组的实验数据是秦俑的两个胸部碎块数据模型(图 4(b))，这两个碎块在拼接的断裂面处存在一定的细小碎块的缺失，不能完全光滑拼接，因此匹配难度较大。该实验中选取的旋转角偏差为$\Delta {\theta _x}=\Delta {\theta _y}=\Delta {\theta _z}=10^\circ $，动态迭代系数$h$=3。

与实验1的碎块匹配流程类似，实验2也是首先采用基于显著性区域的粗匹配算法将碎块进行粗略对齐，然后再采用ICP算法、PICP算法以及改进的ICP算法将碎块的断裂面进一步细匹配，从而完成碎块的最终匹配。对于图 4的两组秦俑碎块，这3种ICP算法的匹配结果和参数(包括匹配的点数、RMS、迭代次数和耗时)分别如图 5和表 2所示。

从图 5和表 2的匹配结果来看，ICP算法不能很好地完成秦俑碎块的匹配，存在较大的匹配误差，而PICP算法和改进的ICP算法能较好地实现秦俑碎块的匹配，并且改进的ICP算法在匹配精度和匹配速度方面比PICP算法又有了很大程度的提高。这是由于PICP算法能够克服噪声对匹配结果的干扰，而改进的ICP算法不仅能够克服噪声对匹配结果的干扰，还由于角度约束和动态迭代系数的加入而大大提高了匹配的精度和收敛速度。因此，改进的ICP算法是一种高精度、快速、抗噪的刚体碎块匹配算法。

破碎刚体复原问题是计算机图形学、计算机视觉和模式识别中一个有待解决的问题。目前，国内外研究的较多是厚度不计的碎片拼接，对于具有一定厚度的刚体碎块的拼接研究较少。针对具有厚度的刚体碎块数据模型，文中采用了先粗配再细配的方法来完成碎块的断裂面匹配。在粗匹配阶段，首先采用区域生长算法进行曲面分割，将碎块外表面分割成一张张的曲面，然后根据曲面的粗糙程度提取出碎块的断裂面，再根据积分不变量提取断裂面上的显著性区域，最后通过计算不同碎块的断裂面的显著性区域的相似性对显著性区域进行匹配，从而实现碎块断裂面的粗略对齐。在细匹配阶段，提出了一种改进的ICP算法来对两个刚体碎块的断裂面进行进一步精确对齐。该算法是在ICP算法的基础上加入高斯概率模型、角度约束和动态迭代系数的方式来改进ICP算法的，使得算法在匹配精度、速度和抗噪性等方面的性能都得到了很大的提高，能够将两个断裂面进行更加精确的对齐。该算法不仅能对公共数据集中的刚体碎块进行完美匹配，而且还能对复杂的秦俑碎块数据进行较好的完全匹配和部分匹配，是一种非常有效的刚体碎块断裂面匹配算法，具有更广阔的应用领域和应用前景。在今后的研究中，要进一步考虑秦俑碎块在部分缺失、风化、二次断裂、受潮和变形等方面的问题，针对更加复杂的情况提出更加高效、精确的秦俑碎块匹配算法，使其能够实现秦俑等特殊刚体碎块数据模型的完美匹配，最终实现自动的虚拟秦俑复原。