基于Harris-SIFT算法和全卷积深度预测的显微镜成像的三维重建研究

1、基于Harris-SIFT算法和全卷积深度预测的显微镜成像的三维重建研究近年来，随着人工智能和机器视觉技术的快速发展，基于图像的三维重建方法越来越受到关注，但在面向显微镜应用领域，其主要应用于平面标本的观测，鲜有针对立体标本成像以及建模的显微设备。研发新一代智能显微设备的同时，开发出面向立体标本的三维模型重建和立体标本几何参数测量的技术具有重要意义。作为机器视觉的热门方向，已有诸多方法和理论用于基于图像的三维重建，其中，运动恢复结构法（Structure from Motion，SfM）1是应用最广泛的经典方法之一。基于SfM设计出了多种三维建模系统，如Wu开发的开源软件Visual SfM2

2、，Moulon等人开发的开源软件OpenMVG3，Schnberger等人开发的开源三维建模系统COLMAP4，5。蒋华强等6针对现有SfM算法重建模型存在的点云稀疏等问题，提出一种利用不同匹配数据进行模型重建的算法。但由于SfM仅计算图像之间成功匹配的特征点的三维信息并将其恢复为三维坐标形成三维点云，而图像中含有的特征点信息相对较少7，因此SfM计算得到的点云模型比较稀疏，重建后的模型精度较低。多视图立体方法（Multi-View Stereo，MVS）8可以从物体多个视角图像计算得到场景的密集三维点云。PMVS（Patch-based MVS）5以SfM重建的稀疏点云为输入，利用图像面邻域

3、信息迭代，使用点云扩张策略进行点云扩张和过滤，最后重建出稠密点云。Schnberger等人9以SfM得到的稀疏重建模型和相机姿态为输入，采用深度图融合的方法来恢复出稠密点云。基于学习的MVS方法如Point-MVSNet10，R-MVSNet11，Cascade-MVSNet12等采用深度图融合的方法恢复较高精度的场景稠密点云也取得不错的结果。以上三维重建算法对拍摄设备要求较高，适用于对大场景成像的纹理等细节丰富的场景，而在实体显微镜应用场景中，由于显微镜工作环境、珠宝等立体标本在显微镜局部放大观测时存在弱纹理、细节不丰富和反光等问题，运用上述算法得到的三维重建效果并不理想。本文提出了一种适用

4、于显微镜应用场景下基于特征提取的多视图立体三维重建算法。以基恩士VHX-6000高级体视显微镜作为实验平台，以集成电路板元器件或珠宝首饰作为实验对象，通过显微镜采集物体不同角度的图像序列（大于10张），对立体标本显微镜图像进行三维重建，取得较好效果。2 相关工作2.1图像采集与预处理显微镜具有高放大率和窄场的特点，对立体标本的图像采集往往只能获得物体较小的局部放大图像，所采集的显微图像包含的纹理信息较少属于弱纹理图像。因此，为了增强图像的细节信息和纹理信息以便后续过程可以提取到图像更多的特征信息，使用直方图均衡化方法增强图像的细节信息以使图像灰度分布更均匀13，使用一种融合光照和彩色信息的灰度

5、图像补偿算法增强图像的纹理信息对图像进行灰度化14。本文的实验设备为基恩士VHX-6000高端数码显微系统，如图1所示。图1基恩士VHX-6000Fig.1KEYENCE VHX-60002.2改进的Harris-SIFT匹配算法特征点提取与匹配是基于图像三维重建算法的关键一步，特征提取与匹配的精度直接影响重建效果。Harris15角点检测算法提取的特征角点比较全面，具有平移和旋转不变性，对参数变化、噪声、图像对比度和光照变化等有很好的鲁棒性，但对图像尺度变化敏感，不具有尺度不变特性。尺度不变特征变换（Scale-invariant Feature Transform， SIFT）16算法对旋

6、转、尺度、亮度等变化具有不变性，对视角转换、噪声、仿射变换有很高的容忍度，但由于算法实现过程中需要通过高斯模糊、降采样等方式构建高斯金字塔，得到的特征点精度不是很高，对边缘光滑的目标无法准确提取特征点，且计算量大，效率低。本文采用一种将Harris与SIFT算法相结合的方法17，首先，用Harris进行图像中的特征角点提取，使用特征角点包含的信息熵值对特征角点进一步筛选；其次，使用马氏距离度量特征点之间的相似性；再次，用随机抽样一致性（RANSAC）方法对匹配模型进行优化，筛选匹配点对，提高匹配的正确率；最后使用显微镜采集的显微图像进行对比实验，验证本算法的匹配准确率。2.2.1Harris角

7、点检测Harrs角点检测算法是利用移动的局部窗口在图像中计算灰度变化值，若在任意方向上移动都存在较大的灰度变化就认为该窗口中存在角点。通过设置合理的阈值来决定检测到的角点数。Harris角点检测变化如图2所示。图2Harris角点检测变化Fig.2Change of Harris corner detection窗口由像素点(x,y)移动(u,v)产生的灰度变化E可以表示为：E(x,y)=x,yW(x,y)I(x+u,y+v)I(x,y)2（1）其中：I是图像灰度函数，W是窗口函数，将式（1）用泰勒级数化简可得：E=u,vMu,vT（2）其中：M=W(x,y)I2xIxIyIxIyI2y，为自

8、相关矩阵；Ix，Iy分别为像素点在x方向和y方向上的梯度值。为了在显微图像中可能存在的弱纹理区域检测到尽可能多的特征点，对图像逐像素点计算灰度变化。使自相关矩阵M的权重x,yW(x,y)=1，则有M=I2xIxIyIxIyI2y=H(x,y)，表示图像每个像素点在x，y方向的二阶偏导数和x、y方向的导数的Hessian矩阵。从图像中提取到的角点可表示为：R(x,y)=detH(x,y)trH(x,y)2（3）其中：R(x,y)为角点响应函数；det(H)为矩阵H的行列式；tr(H)为矩阵H的迹；为常数（取值范围0.040.06），在本文实验中函数R的局部极大值点即为特征角点，此时取=0.04。

9、通过式（3）初步检测出图像特征点。但由于采用逐一计算每一个像素点的局部窗口内灰度变化，因此检测结果仅反映图像局部变化情况，并会将背景中噪声点或纹理不清晰的像素点判断为特征点。熵值18的大小与局部区域包含的信息量成正比，背景中噪声点或纹理不清晰的像素点包含信息量较少导致计算出的熵值低，熵值大的特征点比熵值小的特征点匹配的准确率高。用以下两式进行熵计算：p=ijf(i,j)/N2（4）Hf=i=1nj=1npij(1p)ij（5）其中：pij为图像某像素处的灰度值与周围像素的灰度分布的综合特征，N为图像尺度，Hf为图像二维熵，n,i,j为像素灰度值。分别对图像中初步检测出的特征角点利用式（5）在以

10、特征角点为中心的nn邻域内计算图像块的局部熵，统计每幅图像中特征点熵的均值，取熵均值的80%作为阈值，删除低于阈值的特征点，对特征点进行提纯。对于图像弱纹理区域统计的熵均值较低，可以删除较少的点进而在弱纹理区域保留较多的特征点，从而提高算法在弱纹理区域匹配的正确率，同时提高算法的效率。2.2.2生成SIFT特征描述子Harris算法提取到特征角点，还需要采用SIFT算法来描述并建立特征向量。SIFT通过利用特征角点邻域像素的梯度方向分布特性为每个关键点指定方向参数，使算子具备旋转不变性。m(x,y)=(L(x+1,y)L(x1,y)2+(L(x,y+1)L(x,y1)2（6）(x,y)=arc

11、tanL(x,y+1)L(x,y1)L(x+1,y)L(x1,y)（7）其中：m(x,y)为(x,y)处梯度的幅度，(x,y)为(x,y)处的梯度方向，L(x,y)为所有关键点的尺度。在特征点周围nn邻域内以尺度大小建立一个大采样区域，用式（6）、（7）计算窗口内像素点的梯度方向和大小，然后在该邻域内构建一个小采样区域，并用梯度直方图统计该小采样区域中每个“格子”内像素的梯度方向和大小，形成种子点，从而构建出特征描述子。在本文实验中，SIFT首先将采样区域旋转至特征点方向上，然后在其周围1616的邻域内计算梯度方向和大小。然后在每个44的采样窗口中计算8个方向的梯度直方图，绘制每个方向的累加值

12、，即可形成一个种子点。如图3所示，每个特征角点有44个种子点，每个种子点有8个方向，即每个特征角点由448=128维的描述子组成。图3SIFT描述特征点Fig.3SIFT description of feature points2.2.3特征匹配为了提高对显微弱纹理图像的匹配精度，用马氏距离19替代原始SIFT算法中的欧氏距离作为两幅图像中特征角点之间的的相似性判定度量。欧式距离将特征点的不同属性之间的差别等同看待，有时不能满足实际要求。马氏距离与欧式距离的不同之处在于：它考虑到特征点各种特性之间的联系；并且能对两个样本集的相似性进行有效度量，不受尺度变化及变量之间相关性的影响；它还可以排除

13、变量之间的相关性的干扰。马氏距离可以表示为：S=EXE(X)YE(Y)（8）DM(X,Y)=(XY)TS1(XY)（9）其中：X，Y为两个特征向量，S为两组特征向量的协方差矩阵，E为均值。从式（9）中可以发现，马氏距离计算是建立在总体样本的基础上，即计算的是两个向量间的整体距离。相比欧式距离，马氏距离多了一个协方差矩阵的逆矩阵，能更准确的度量特征点之间的相似性。虽然增加了计算量，但马氏距离把方差归一化，使得特征之间的相似关系更加合理，匹配精度更高。为进一步提高算法在弱纹理显微图像上的匹配准确率，需要剔除掉误匹配点，因此采用随机抽样一致性（RANSAC）20算法替代原SIFT算法中的比例阈值法，

14、选择一组局域熵值高的特征点作为初始模型，用初始模型去测试所有的其它数据，如果某个点适用于估计的模型则认为它也是局内点，否则是局外点；然后，用所有假设的局内点去重新估计模型。算法通过K次迭代，可精确估计匹配模型参数。2.3实体显微镜场景下物体的三维重建2.3.1稀疏点云重建利用基恩士VHX-6000可变倍镜头可将智能显微镜镜头调整到一个固定的角度，然后上下、前后、左右等多方向移动载物台，从不同视角拍摄立体样本获得图像序列。本文采用如文献4中的增量SfM方法重建物体的三维结构。增量SfM是具有迭代重建组件的顺序处理方法。首先，用本文算法对图像进行特征点提取与匹配，采用包含正确匹配特征点对最多的两视

15、图作为重建种子，通过它们之间的匹配点对估计出本征矩阵，由本征矩阵采用SVD分解法计算出相机位姿。然后，依次将下一最佳匹配图像Inext（匹配点对数量最多）添加进模型迭代，重复估计相机位姿并对相机进行标定；同时将场景中匹配的点通过三角化来生成三维点，并计算出相机的运动轨迹和场景中物体的三维信息。最后，对所有已生成的三维点和已估计出的相机位姿使用光束平差法（Bundle Adjustment，BA）进行优化。BA结束，得到相机估计参数和场景几何信息，即稀疏的3D点云。2.3.1.1有序图像输入由于增量SfM需要对所有图像依次进行匹配，用无序图像进行模型重建时由于图像的无序导致图像位姿的先验条件未知

16、，在匹配环节存在多张图片之间的交叉匹配和多次交叉验证匹配点，所以计算量大，实时性差。为提高图像的匹配效率，采用有序图像序列作为输入。通常相邻拍摄的两张图像所包含的共同特征信息最多，也更加容易正确匹配。在采集图像时根据显微镜载物台移动的轨迹以及图像拍摄时间对所采集的物体不同视角图像添加标签，将图像根据载物台运行轨迹和时间的先后排列成一组顺序固定的图像序列（类似视频中的帧）。显微镜载物台工作距离较小，拍摄的相邻图片中两两之间物体信息的重合很多，包含的特征点信息最多；在特征匹配时使相邻的图片先进行特征匹配。在SfM算法选好初始的两张匹配图象时默认把与这两张图像相邻的图像作为下一最佳匹配图像Inext

17、进行重建，依次将与上一张图像邻近的图像作为下一最佳输入进行迭代，这样所有图像只需添加匹配一次即可完成所有图像的重建。2.3.1.2三角化SfM在图像正确匹配和估计得到相机参数之后，就可以将图像场景之间正确匹配的特征点通过三角化方法恢复为空间三维点，进而得到三维点云。三角化的方法如图4所示。图4三角化Fig.4Triangulation在图4中，I1，I2为两幅图像，P1，P2分别为两幅图像上匹配的特征点；O1，O2分别为两幅图像拍摄视角时的相机中心，将每幅图像的相机中心与其特征点的连线延长在空间中交于点P，则P点为图像上的点所恢复的世界坐标系下的三维点。得到三维点之后通过BA调整，过滤离群值，

18、估计得到物体的三维点云。2.3.2稠密点云重建SfM重建得到的点云模型为稀疏点云，只有大致的模型结构，缺少了很多细节部分，所以还需要将稀疏点云进行稠密化。多视图立体（MVS）是一种将SfM输出的稀疏点云稠密化的方法，MVS重建利用SfM算法输出的稀疏点云、相机参数及相机位姿，然后融合图像的深度图获得物体的稠密三维点云。但对于显微镜工作场景下拍摄的物体局部存在弱纹理的图像，由于表面灰度信息有限，传统MVS方法图像深度信息估计误差较大。因此，针对上述在显微镜应用场景中的物体稠密重建问题，利用CNN网络与MVS结合的方法来估计和预测采集的显微图像的深度图，进而使用深度图进行稠密重建。为保证模型收敛，

19、本文采用文献21提出的基于深度残差网络（ResNet-50）22-23的全卷积神经网络对采集的显微图像序列进行深度估计和预测。在ResNet-50的基础上设置了反卷积层和上采样层逐步增加每层的输出大小，并在反卷积层输出后增加了激活函数（ReLu）24，对输入图像进行同样比例大小的深度预测。通过比较全卷积神经网络预测的深度与MVS所得深度的绝对差值并设置阈值，将大于阈值的预测深度融合到MVS估计的深度图中得到改善的深度图，然后利用SfM得到的稀疏点云和相机的位姿信息与改善后的深度图恢复出物体的稠密点云。用于预测图像深度的全卷积网络结构如图5所示。为了增加系统的实时性并减少计算时内存的占用量，去掉

20、了原ResNet-50中的平均池化层和全连接层部分，仅保留ResNet-50 的卷积层结构来提取图像中的特征信息，从而可以输出尺寸为77，2 048通道的特征图，减少了计算参数，加速模型的训练，提高预测效率。图5深度预测网络结构Fig.5Depth prediction network structure本文的研究旨在恢复显微镜应用场景中所拍摄物体的近似真实的三维结构。重建物体近似真实的三维结构是以正确估计图像深度为前提的。所以，为了尽可能正确地预测图像的深度信息，约束预测的图像深度值和真实值之间的误差。结合文献21采用平均平方误差作为损失函数，损失函数L定义为：L(y)=i=1mj=1n(1

21、yijtij+)2mn（10）其中：在图像(i,j)处预测的深度值为yij，在图像(i,j)处真实深度值为tij，图像包含的像素点数为mn，为常数。L将预测深度图像素点的预测值和真实值之间的平均误差定义为损失，该损失函数对预测的深度值y惩罚较严格，可以很好地减少显微镜载物台背景区域带来的深度误差影响，进一步约束了前景物体的深度误差。3 实验结果与分析本文的实验平台为基恩士VHX-6000高性能实体显微镜系统，以电子电路板上元件为实体标本采集的立体标本图像序列来验证本算法的可行性。实验PC平台配置为：系统为Ubuntu20.04.3LTS，CPU为英特尔Xeon E5-1630 V3，主频为3.

22、70 GHz，内存为16 G，显卡为英伟达Quadro K2200，显存为4 G，硬盘为2 T固态硬盘。如图1所示的基恩士VHX-6000智能显微镜可以通过人为转动手柄灵活调整显微镜物镜的观测角度，以物镜垂直载物台为原点可以任意左右倾斜090，并通过上下、旋转以及水平任意角度移动载物台实现对立体标本的多视角观测。在本文实验中将显微镜物镜向左调整倾斜60，放大倍数设置为20倍，然后移动载物台采集了36张电路板上的电容的多视角图像序列。对图像根据载物台移动轨迹和采集时间进行排序，使其成为有序图像序列。采集的图像如图6所示。图6电容的多视角图像Fig.6Multi-view image of cap

23、acitance随机从拍摄的36张图像中选取2张图像进行图像匹配实验，选取的图片如图7所示。将选取的2张图像经过本文的图像灰度化补偿算法进行预处理，然后与Harris和SIFT匹配算法经行实验比较，实验数据如表1所示，三种匹配算法的实验结果图如图8所示。图7选取的待匹配图像Fig.7Selected image to be matched表1不同匹配算法实验数据Tab.1Experimental data of different matching algorithmsMatching algorithmImg_L feature_pointImg_R feature_pointMatchin

24、g_pointMatching_accuracyTime/sHarris1 4031 3721 05780.61%6.02SIFT1 7861 63172391.53%11.64Ours85092739294.82%4.61图8三种匹配算法实验结果Fig.8Experimental results of three matching algorithms可以发现，尽管其他算法能提取出较多的特征点，但本文算法在针对显微图像中弱纹理区域的匹配上具有明显优势。如图8所示，本文算法在对电容壁和电容顶部弱纹理区域的特征点匹配中的误匹配较少，而其他算法存在明显的误匹配现象。如表1，本算法也具有更好的性能指

25、标，原因是本文算法在用Harris提取出特征角点之后又利用熵值法对特征点进一步筛选，减少了后续特征点匹配过程的计算量，减少了运算时间；筛选出的包含丰富图像特征信息的点，以及在匹配环节采用马氏距离代替欧氏距离，可更好地计算弱纹理区域特征点之间的相似性，提高了匹配的准确率。匹配准确率相比SIFT算法提高了3.29%，时间比Harris算法减少了23.42%。本文在稀疏重建环节改进了特征提取与匹配，将使用无序图像进行增量SfM重建改为使用有序图像重建。对图6所示的36张图像进行稀疏重建并与SfM算法和OpenMVG算法进行实验比较，实验结果如图9所示，实验数据如表2所示。图9稀疏点云重建Fig.9S

26、parse point cloud reconstruction表2SfM稀疏重建实验数据Tab.2Experimental data of SfM sparse reconstructionAlgorithmNum_image_matchesTime/minNum_point_cloudOpenMVG3 10710.83913 257SFM3 05711.14614 088Ours4 2819.72416 440可以看到改进后的SfM重建结果具有更好的结构完整性，本文算法在电容壁处的弱纹理区域重建后的点云比较密集，本文算法相比SfM算法整体点云数量增加16.69%，相比OpenMVG算法重建

27、时间节省了10.29%。在使用本文深度预测网络模型进行深度预测前先对模型进行训练。本文使用基恩士-VHX6000超景深显微镜拍摄的宝石、电子元器件等不同颜色、形状、材质的立体标本的多视角显微图像构建的本地图像数据集，随机从中取出部分图像构成训练集和验证集，对深度预测网络进行监督训练。基恩士-VHX6000超景深显微镜具有3D测量功能可以同步扫描立体标本的深度信息生成深度图，即在拍摄立体标本多视角图像时可以同步扫描同一视角的深度图作为模型训练的输入图像，图像对应的深度图作为训练标签。模型的前半部分即卷积层部分使用ResNet50在公开数据集ImageNet25上预训练好的权重参数作为初始化参数，

28、模型训练时用AdaDelta优化训练效率26，通过批处理的随机梯度下降法优化网络参数，本文训练中batch size设置为4，迭代过程中采用了“早停”机制，即在迭代过程中观察每次迭代输出预测值与真实值的误差，当在验证集上预测值与真实值的误差不再明显变化时停止训练，此时损失函数L在0.217附近小幅波动、没有明显变化，说明模型收敛，完成本文模型的训练。图10展示了部分预测结果，从预测的深度图可以看出立体标本“电容器”的具体形状，深度变化均匀，很好地预测了立体标本的深度。但不足之处在于，在场景背景处与“电容器”壁身弱纹理处存在噪声以及部分缺失，这在之后用MVS深度图进行修正并利用修正后的深度图恢复

29、三维结构时进行了弥补。图10深度预测结果Fig.10Depth prediction result通过以SfM稀疏重建阶段获得的物体稀疏点云和估计的场景几何参数为基础，用本文全卷积神经网络方法预测的图像深度和MVS估计的深度图相融合，以增强显微镜场景下对弱纹理区域的表征效果，重建出物体的稠密点云。将本文的稠密重建方法的重建效果与开源的三维建模系统COLMAP中的MVS方法相比较，估计的深度图如图11所示，重建后的密集点云如图12所示。图11深度图Fig.11Depth map图12稠密重建后的点云Fig.12Densely reconstructed point cloud在深度图估计方面，C

30、OLMAP中MVS方法估计出的深度图会受显微镜聚焦位置的影响，与场景中物体真实结构深度信息不一致。本文方法用全卷积神经网络预测的图像深度对MVS方法估计的深度图进行矫正，得到的深度图更加精确。在稠密点云方面，可以明显看到本文方法比COLMAP得到的稠密点云具有更好的结构完整性，在图12中电容周身弱纹理区域的重建结构更加完整，重建得到的点云数量明显多于COLMAP。稠密重建实验数据如表3所示。表3稠密重建实验数据Tab.3Dense reconstruction experimental dataAlgorithmRMSERELLOG10Time/minNum_point_cloudSSIMCOLMAP2.3960.1210.07312.335560 0120.64Ours1.1270.0860.0429.724735 0160.89用均方根误差（VRMSE）、平均相对误差（VREL）和平均log10误差（VLOG10）这3种常见的预测评价指标验证本文使用的全卷积神经网络模型预测深度的精确性，表示为：VRMSE=i=1N(yidi)2N（11）VREL=i=1N|yidi|yiN（12）VLOG10=i=1N|log10yilog10di|N（13）其中：yi为第i个像素点的真实深度，di为i个像素点的深度预测值，N为图像像素点总数。可以发