基于PCA和FFT的青铜器三维碎片虚拟拼接方法与流程

本发明涉及青铜器文物虚拟复原方法，尤其涉及青铜器三维碎片点云特征匹配和拼接方法，具体为基于PCA和FFT的青铜器三维碎片虚拟拼接方法。

背景技术：

青铜器是我国古代主要的生产和祭祀用品，代表了先进的技术和艺术，是重要的文物的组成部分。但是，随着时间的流逝，这些年代久远的历史文化遗产，在经受了若干年的自然风化、侵蚀和人为破坏之后，多数已经变得残缺不全，成为若干个形状任意的子物体。

目前，主要是依靠手工对考古发现的大量破碎青铜器文物进行复原。这样，不仅在操作过程中容易损坏，而且青铜器修复质量极大地依赖技术人员个人的知识和经验，拼接速度缓慢，以至于大量碎片长期堆积，加之人工操作误拼误修在所难免，其损失不可估量。实践证明利用计算机技术辅助文物复原技术，既可降低文物修复难度，加快其复原速度，也可避免人工修复过程对文物的二次损害。因此，利用计算机虚拟复原技术辅助实现破碎青铜器的虚拟拼接是必然选择。

目前，文物虚拟拼接方法主要分为两类：(1)针对已知几何形状进行匹配，主要是基于对称轴进行匹配，但由于青铜器碎片往往存在形变，所以利用该方法匹配其容易产生误差。(2)针对未知几何形状的碎片，如陶瓷和壁画碎片等，主要是利用碎片彼此的轮廓曲线、光度或纹理信息的匹配程度为度量进行虚拟拼接，而青铜器碎片边缘轮廓往往不规则，且碎片的纹理常常被锈迹所覆盖，因此基于轮廓、光度和纹理拼接难度较大。

新近有文献“Pairwise matching of 3D fragments using fast fourier transform”利用快速傅里叶变换拼接石器取得成功，证明该方法对断面缺损和噪声影响鲁棒性好，但与之相比，青铜器碎片边缘薄，缺乏断裂面特征，所以，此方法并不能直接适用。

为此，需要有一种方法来专门解决边缘薄、易腐蚀、有缺损的青铜器碎片虚拟拼接的问题。

技术实现要素：

本发明为了解决青铜器碎片拼接中边缘薄、特征点少、边缘被腐蚀甚至缺损的问题，提供了一种基于PCA和FFT的青铜器三维碎片虚拟拼接方法。

本发明是采用如下的技术方案实现的：基于PCA和FFT的青铜器三维碎片虚拟拼接方法，包括以下步骤：

(1)用精度在0.03mm以上三维扫描仪分别采集两青铜器碎片的三维点云数据，使用geomagic studio软件对三维点云数据提取碎片轮廓点形成碎片轮廓并以1mm进行重采样；

(2)对待拼接的两个碎片轮廓分别进行分段：计算第一个碎片轮廓的曲率其中，r＝(x,y,z)，x,y,z分别为轮廓点的方向坐标，r′和r″分别为r的一阶导数和二阶导数，||·||³为矩阵二范数的三次幂，以曲率局部极大值点作为角点，以每个角点及其周围m个点构成一个待匹配特征曲线段；对于第二个碎片轮廓，以轮廓上每个点及其周围m个点构成一个待匹配特征曲线段，则待匹配特征曲线段的长度均为m+1个点；

(3)利用PCA为每个待匹配特征曲线段建立一个局部坐标系，以第一个碎片轮廓r₁上的第i个待匹配特征曲线段r_1i(i小于或等于所提取的角点个数)为例，设点q_j(j＝1,…,m+1)为待匹配特征曲线段r_1i上的点，则该待匹配特征曲线段的协方差矩阵为：其中，设λ₁,λ₂,λ₃(λ₁≥λ₂≥λ₃)和e₁′,e₂′,e₃′分别为协方差矩阵对应的特征值和特征向量，则局部坐标系的三个垂直的基向量分别为e₁′,e₂′,e₃′；

(4)把待匹配特征曲线段从全局坐标系转换到局部坐标系：全局坐标系为{o,e₁,e₂,e₃}，局部坐标系为{o′,e₁′,e₂′,e₃′}，局部坐标系的基向量e′_j在全局坐标系下的坐标是(x_j,y_j,z_j)，j＝1,2,3，则从全局坐标系转换到局部坐标系的过渡矩阵为：A＝(e′₁,e′₂,e′₃)，若局部坐标系原点在全局坐标系下的坐标为(x₀,y₀,z₀)，待匹配特征曲线段上某点在全局坐标系坐标为(x,y,z)，则其在局部坐标系下的坐标为(x′,y′,z′)；

(5)将两个碎片轮廓的待匹配特征曲线段投影到x-z平面和y-z平面，分别提取它们在同一平面的特征投影曲线段的特征函数：式中，g＝1或2分别表示第一个碎片轮廓和第二个碎片轮廓，u_gk表示第g个碎片轮廓的投影曲线是由k个点构成的点集，u_gk的生成方法是：在x-z平面的x轴上和y-z平面的y轴上分别根据实验选取一个区间，将两区间都分为N+1部分(N为2的整数次幂)，划分区间的点集即为{u_gk}(k＝0,…N-1)；对于点集{u_gk}中的每个点，在其特征投影曲线段上都有距离该点最近的点p_l(l＝1,…M)，M为与u_gk中每个点邻近的点的个数；定义式中d(u_gk,p_l)为u_gk中的点和p_l间的距离，w_l为反距离值；z_pl是点p_l的z轴坐标；

(6)计算特征投影曲线段的特征函数的傅里叶系数，a_k＝Re(FFT(f(u_gk)))，b_k＝Im(FFT(f(u_gk)))，FFT(f(u_gk))是特征函数f(u_gk)的快速傅里叶变换，Re(FFT(f(u_gk))和Im(FFT(f(u_gk))分别为FFT(f(u_gk))的实部和虚部；

(7)计算特征投影曲线段的傅里叶光谱能量：其中，a_k和b_k为傅里叶系数；

(8)计算特征投影曲线段间的总能量差异值dd，人工选取阈值T，若dd≤T，则该两特征投影曲线段为粗匹配对；

(9)将两个碎片轮廓分别通过平移旋转变换到与其粗匹配特征曲线段对应的局部坐标系中，并将轮廓投影到x-y平面上，计算其中一个轮廓上任意一点P与另一轮廓上离它最近的点Q之间的欧氏距离d_PQ，自定义阈值t，若d_PQ≤t则属于匹配的点数，计算两轮廓可匹配且是连续匹配的点数，即为两碎片轮廓的公共长度，计算每种可能的匹配的轮廓公共长度，选取公共长度最大的匹配对为最终匹配结果，实现碎片轮廓的精匹配。

本发明用到的基本方法是主成分分析法(PCA)和快速傅里叶变换(Fast Fourier Transform，FFT)，其中，前者是一种数据简化技术，主要用于在新的坐标系统中减少数据集的维数并保留最大的特征，在本发明中用该方法实现了坐标系的变换；而后者是将空间域信号转换成易于分析的频域信号的经典方法，在本发明中通过FFT将点云碎片的轮廓用频域信号表示，便于特征匹配。本发明不仅将PCA和FFT两种方法相结合，而且提出将两个待匹配的轮廓特征曲线段先投影到两个坐标平面进行粗匹配，再将粗匹配结果投影到第三个坐标平面进行精匹配，不仅保证了匹配精度，而且还有利于降低运算量。

上述基于PCA和FFT的青铜器三维碎片虚拟拼接方法，特征曲线段选取点数越少，特征就越不明显，可能的匹配就越多，精匹配计算量就越大，且易造成误匹配，而点数越大，精匹配速度虽然快了，但是对公共长度要求较高，可能会找不到对应的匹配，所以，本发明每个待匹配特征曲线段长度(即：前述的m+1)的取值范围优选为14-16，不仅可以获得较好的匹配结果，还能保证具有较快的运行速度。

上述基于PCA和FFT的青铜器三维碎片虚拟拼接方法，提取特征函数时在x-z平面和y-z平面选取的区间分别为[-5,5]，[-1,1]，M＝4，N＝4或8，可以获得较好的拼接效果。

上述基于PCA和FFT的青铜器三维碎片虚拟拼接方法，特征曲线段间的总能量差异值的合适阈值T选取区间为0.03-0.05。精匹配时自定义阈值t选取区间为0.5-1.5，可以获得较好的拼接效果。

附图2、图3为青铜器碎片虚拟拼接的实例，其中，图2为青铜器碎片点云图像，图3为本方法的青铜器碎片点云拼接效果图像。

本发明将碎片三维轮廓匹配问题转换成特征曲线段匹配，将空间曲线匹配转换成三视图匹配，并进一步分解为二次分级匹配(包括粗匹配和精匹配)，解决了青铜器碎片拼接中边缘薄、特征点少、边缘被腐蚀甚至缺损的问题。

附图说明

图1为本发明的流程图。

图2为青铜器碎片点云图像。

图3为青铜器碎片点云拼接效果图像。

具体实施方式

基于PCA和FFT的青铜器三维碎片虚拟拼接方法，包括以下步骤：

(1)用精度在0.03mm以上三维扫描仪分别采集两青铜器碎片的三维点云数据，使用geomagic studio软件对三维点云数据提取碎片轮廓点形成碎片轮廓并以1mm进行重采样；

(2)对待拼接的两个碎片轮廓分别进行分段：计算第一个碎片轮廓的曲率其中，r＝(x,y,z)，x,y,z分别为轮廓点的方向坐标，r′和r″分别为r的一阶导数和二阶导数，||·||³为矩阵二范数的三次幂，以曲率局部极大值点作为角点；以每个角点及其周围m个点构成一个待匹配特征曲线段；对于第二个轮廓，以轮廓上每个点及其周围m个点构成一个待匹配特征曲线段，则两个曲线段的长度均为m+1个点；

(3)利用PCA为每个待匹配特征曲线段建立一个局部坐标系，以第一个轮廓r₁上的第i个特征曲线段r_1i(i小于或等于上所提取的角点个数)为例，设点q_j(j＝1,…,m+1)为待匹配特征曲线段r_1i上的点，则该特征曲线段的协方差矩阵为：其中，设λ₁,λ₂,λ₃(λ₁≥λ₂≥λ₃)和e₁′,e₂′,e₃′分别为协方差矩阵对应的特征值和特征向量，则局部坐标系的三个垂直的基向量分别为e₁′,e₂′,e₃′；

(4)完成从全局坐标系到局部坐标系的转换：全局坐标系为{o,e₁,e₂,e₃}，局部坐标系为{o′,e₁′,e₂′,e₃′}，局部坐标系的基向量e′_j在全局坐标系下的坐标是(x_j,y_j,z_j)，j＝1,2,3，则从全局坐标系转换到局部坐标系的过渡矩阵为：A＝(e′₁,e′₂,e′₃)，若局部坐标系原点在全局坐标系下的坐标为(x₀,y₀,z₀)，待匹配特征曲线段上某点在全局坐标系和局部坐标系下的坐标分别为(x,y,z)和(x′,y′,z′)；

(5)将两个碎片轮廓的待匹配特征曲线段投影到x-z平面和y-z平面，分别提取它们在同一平面的特征投影曲线段的特征函数：其中,u_gk表示第g个轮廓投影曲线是由k个点构成的点集，g＝1或2分别表示第一个碎片轮廓和第二个碎片轮廓，u_gk的生成方法是:在x-z平面的x轴上和y-z平面的y轴上分别根据实验选取一个区间，将该区间分为N+1部分(N为2的整数次幂)，划分区间的点集即为{u_gk}(k＝0,…N-1)；对于点集{u_gk}中的每个点，在其特征投影曲线段上都有距离该点最近的点p_l(l＝1,…M),M为与u_gk中每个点邻近的点的个数；定义式中d(u_gk,p_l)为u_gk的点和p_l间的距离，w_l为反距离值，z_pl是点p_l的z轴坐标；

(6)计算特征投影曲线段的特征函数的傅里叶系数a_k＝Re(FFT(f(u_gk)))，b_k＝Im(FFT(f(u_gk)))，FFT(f(u_gk))是特征函数f(u_gk)的快速傅里叶变换，Re(FFT(f(u_gk))和Im(FFT(f(u_gk))分别为FFT(f(u_gk))的实部和虚部；

(7)计算特征投影曲线段的傅里叶光谱能量：其中，a_k和b_k为傅里叶系数；

(8)计算特征投影曲线段间的总能量差异值dd，人工选取合适阈值T，若dd≤T，则该两特征投影曲线段为粗匹配对；

(9)将两个碎片轮廓分别通过平移旋转变换到与其粗匹配特征曲线段对应的局部坐标系中，并将轮廓投影到x-y平面上，计算其中一个轮廓上任意一点P与另一轮廓上离它最近的点Q之间的距离d_PQ，自定义阈值t，若d_PQ≤t则属于匹配的点数，计算两轮廓可匹配且是连续匹配的点数，即为两轮廓的公共长度。计算每种可能的匹配的轮廓公共长度，选取公共长度最大的匹配对为最终匹配结果。

上述基于PCA和快速傅里叶变换的青铜器三维碎片拼接方法，轮廓分段选取14-16个点为一段。

上述基于PCA和快速傅里叶变换的青铜器三维碎片虚拟拼接方法，提取特征函数时在x-z平面和y-z平面选取的区间分别为[-5,5]，[-1,1]，M＝4，N＝4或8，可以获得较好的拼接效果。

上述基于PCA和快速傅里叶变换的青铜器三维碎片虚拟拼接方法，特征曲线段间的总能量差异值的合适阈值T选取区间为0.03-0.05。精匹配时自定义阈值t选取区间为0.5-1.5，可以获得较好的拼接效果。