1　引言

光学显微镜由于其特殊的作用广泛应用于工业检测、科学研究和医学分析等领域，但传统光学显微镜视场受限于放大倍率，在高放大倍率下只能观察到样本的部分内容，这给使用者带来许多不便。通过图像拼接技术可将同一显微切片的图像合成具备较宽视野和较高分辨率的图像，从而解决显微镜视场受限于放大倍率的问题。显微图像通常存在大量复杂且相似的微观结构，对配准算法鲁棒性要求较高。随着显微相机成像分辨率的提高，基于特征的拼接方法提取到大量特征进行匹配，计算量极大，因此，研究高效的显微图像拼接算法具有重要意义。

基于特征的拼接方法是图像拼接应用最广泛的方法之一，特征包括图像中角点、线、轮廓、封闭区域等^［1］。其中，以加速稳健特征（SURF）算法^［2］为基础，设计出了多种满足不同应用场景的拼接算法，如赵立杰等^［3］针对活性污泥显微图像表征信息有限等问题，提出通过特征检测的方式获取距离矩阵，再采用Floyd算法优化拼接路径实现多图拼接。贾仁庆等^［4］为了解决荧光图像和明场图像位置不匹配的问题，提出通过SURF算法采集图像特征后，利用粒子群优化算法对图像进行配准，有效提高浮游植物鉴定的精度。张寅等^［5］为了解决遥感图像检测目标小等难题，提出了先通过特征检测算法建立高低层级特征融合结构，利用特征融合模块增强小目标的特征信息，从而实现小目标检测。Liu等^［6］针对无序显微图像拼接耗时长的问题，提出通过特征检测建立基于合并排序的图像拼接模型，在不遍历所有图像的情况下，直接对匹配后的图像进行拼接。

显微图像存在大量相似结构是造成特征点匹配效率低的主要原因，为了避开特征点的检测，部分学者提出选择线特征作为检测对象^［7］进行图像配准、通过模板匹配进行初始运动估计再应用光流法进一步匹配^［8］、边缘检测结合线性搜索法配准图像^［9-10］等方法拼接显微图像。同时，通过深度学习方式^［11-12］进行图像拼接，需要大量数据和高性能图形处理器（GPU）训练数据模型。因此，显微设备无法完美地部署深度学习算法。开源图像处理软件OpenPano、AutoStitch、ImageJ提供了图像拼接功能插件^［13］，这为学者们所提的拼接算法提供比较对象。

上述改进的拼接算法和图像拼接系统中，图像采集与处理独立，图像拼接是后处理过程。当图像分辨率较高时，难以实现快速配准，因此，拼接时间和效果成为了衡量拼接算法的标准。本文提出一种适用于显微镜应用场景下基于载物台运动信息的图像拼接算法。以配备电动载物台的光学显微镜作为实验平台，通过载物台运动信息改进特征匹配算法。以特征丰富且微观结构相似的脊髓横切片、双子叶植物茎切片和特征稀疏的植物花粉作为实验对象，通过显微镜采集不同重叠区域的图像进行拼接，该实验结果优于AutoStitch软件拼接的结果。

2　相关工作

2.1　载物台运动信息与成像间的关系

显微镜载物台XYZ轴分别通过丝杆与步进电机组合控制，XY轴控制载物台水平和竖直平面方向的移动，Z轴控制聚焦工作，如图1所示。图像拼接时，载物台以“弓”字形移动，图像采集顺序如图2所示，灰色区域为图像重叠区域，区域的大小由电动载物台移动的距离控制。

图 1. 运动平台

Fig. 1. Motion platform

下载图片 查看所有图片

图 2. 拼接路径

Fig. 2. Stitching path

下载图片 查看所有图片

显微镜的视场大小与CCD相机靶面尺寸和放大倍率有着密切关系，视场大小F可表示为

式中：A为靶面尺寸；Γ为光学放大倍率。像元尺寸P可表示为

式中：σ为图像分辨率。

载物台需要丝杆与步进电机搭配完成移动，载物台移动距离与成像间的关系S可表示为

式中：N为丝杆旋转圈数。综合上式，可得

2.2　显微图像结构相似度

显微图像存在大量相似结构是造成特征点对误匹配的主要原因之一。SURF算法对显微图像特征检测后每个特征点用64维描述子表示，相似结构造成特征点间描述子区别不大，在特征点匹配环节极易误匹配，从而影响单应性矩阵求解的准确性。暴力匹配（BF）算法通过一对多的方式计算描述子最近欧氏距离，得到最近欧氏距离的特征点后再反向交叉匹配，若仍然匹配到相同特征点则判定为匹配对，在特征点相似的情况下极易造成误匹配。快速最近邻搜索（Flann）算法利用了k-d树或k-means树算法来对特征进行聚类建模搜索到最近邻与次临近的特征点，描述子相似的情况下也容易造成误匹配。由此可知，显微图像结构的相似度影响着特征匹配的效率，用特征点间的欧氏距离表示结构相似度，欧氏距离越小特征点间越相似，具体可表示为

式中：d(pn,pm)为特征点pn与特征点pm的欧氏距离；xi、yi分别为pn和pm的描述子。则在一张图像中特征点pi的相似度可表示为

式中：Si为相似度；M为检测到的特征点总数；N为阈值下与pi相似的特征点数。

2.3　预测拼接模型

特征点的提取和匹配是基于特征拼接算法的关键一步，其精度直接影响单应性矩阵的准确性，特征点匹配是拼接过程中耗时最长的过程，优化特征点匹配环节是提高拼接效率的关键。

提出一种基于载物台运动信息预测匹配点的拼接算法。通过SURF算法在相邻两张图像重叠区域中提取特征点，图像的重叠区域大小由丝杆旋转圈数决定。通过载物台运动的方向及距离预测待匹配点的位置，筛选出以预测位置为中心固定区域内最小欧氏距离的特征点，将其作为实际匹配点，减少了预测区域外相似特征点的影响。利用载物台运动信息计算匹配特征点对的斜率粗筛选误匹配特征点对，用随机抽样一致性（RANSAC）^［14］算法对匹配模型进行优化并计算单应性矩阵。最后，通过单应性矩阵变换图像完成拼接。

预测拼接模型如图3所示，i、j分别为两张相邻图像，D为载物台移动时相邻图像间重叠区域的像元数，“+”表示通过载物台运动信息预测的待匹配特征点位置，两黑点连线表示实际匹配特征点对，由于机械误差和使用环境的影响，电动载物台在运动过程中并不能保持线性，其存在运动误差，虚线框表示预测匹配点的区域。

图 3. 预测拼接模型

Fig. 3. Predictive stitching model

下载图片 查看所有图片

图3中，图像i中的特征点通过预测拼接模型在图像j中匹配到最优配对的特征点。由于仅在预测区域内搜索匹配特征点，极大缩减了匹配特征点的计算量，减少了预期区域外相似特征点引起的误匹配。

2.3.1　特征检测

采用SURF算法对图像重叠区域进行特征检测与描述，SURF算法通过使用积分图像来计算Hessian矩阵的行列式，从而检测图像中的关键点，并使用Haar小波响应来计算关键点的描述子。与尺度不变特征转换（SIFT）算法相比，SURF算法具有更快的计算速度和更好的尺度不变性能^［15］。

2.3.2　改进特征匹配

图像间的重叠区域由载物台预设运动轨迹和机械定位误差决定，如图4所示，“+”表示预测匹配点坐标，连线为实际配对特征点。重叠区域中的特征点pn从第i平面图像移动到第j个平面图像的特征点pm，ΔD为图像i中特征点pn在图像j中以预测匹配特征点坐标pn'为中心的正方形区域边长。则预测匹配特征点的坐标可表示为

式中：Si为图像i的大小，图像i和图像j大小一致；Pi为图像i中的特征点。计算图像i中每个特征点Pi与其在图像j中预测区域内所有特征点Pj的欧氏距离，d(Pi,Pj)可表示为

式中：xi、yi分别为Pi和Pj的描述子。当其欧氏距离最小时，两特征点为最相似的特征点对。此时

式中：MIN（·）为距离最小值函数。

图 4. 预测匹配特征点

Fig. 4. Predict matching feature points

下载图片 查看所有图片

由于显微镜电动载物的精密结构，其在运动过程中误差范围可控。因此，ΔD的值影响计算效率，ΔD越小，区域内特征点的数量越少，计算效率越高，但设定值过小时会导致区域内无特征点。

2.3.3　斜率粗筛选和RANSAC精匹配

当预测区域内存在多个相似的特征点时，通过预测匹配点的方法计算得到的特征点对仍存在误匹配的情况。为了提高匹配的准确性，引入特征点对斜率粗筛选和RANSAC精匹配的机制。通过载物台运动信息计算出特征点对在同一坐标内的斜率，粗筛选出与载物台运动趋势不同的斜率。同时，准确匹配的特征点对斜率大致相同，计算出斜率分布曲线，将斜率出现次数较低的特征点对剔除如图5所示，从而实现斜率粗筛选匹配对。最后，再采用RANSAC算法对特征点对进行精确匹配，得到更准确的匹配结果。同一坐标系下特征点对的斜率k_i可表示为

式中：xi,yi、xj,yj为特征点的坐标；w为图像的宽度。

图 5. 特征点对斜率粗筛选

Fig. 5. Feature point pair slope coarse screening

下载图片 查看所有图片

RANSAC算法从粗筛选后的数据集中随机选取4对不共线的特征点对，计算出单应性矩阵，并用该矩阵对数据进行变换。通过计算变换后的点和对应真实点之间的误差，判断每个特征点是否为内点。根据内点的数量和误差等因素，评估单应性矩阵的质量，再决定继续迭代和停止迭代的时刻。最终，选择内点数量最多的特征点对计算最终的单应性矩阵，从而实现特征点精确匹配^［14］。

2.4　图像配准

待配准图像通过单应性矩阵变换后，以参考图像为背景，在参考图像坐标系中插入变换后的待配准图像，根据插入位置和尺寸进行拼接。设待配准图像B上的像素坐标为(x,y)，则经过单应性矩阵H变换后的配准图像为B'，B'的像素坐标为(x',y')，则有

式中：h0h1h3h4为尺度及旋转变化量；h2为水平方向位移量；h5为垂直方向位移量；h6h7为水平与垂直方向变化量。显微镜成像仅存在平移关系，无旋转与尺度变换关系，因此，在变换图像时h0h1h3h4可以等于0，式（12）可变换为

以参考图像为中心，待配准图像有4种拼接方式。拼接方向决定了参考图像与待配准图像间的相对位置，如图6所示。当stitching direction为1和3时，即待配准图像在相对于参考图像左侧与上侧做变换。由于参考图像左上角坐标为(0,0)，因此，B'图像的部分像素点坐标为负值超出参考图像的坐标系，仅在参考图像坐标系中显示部分内容，如图7所示。

图 6. 拼接模式

Fig. 6. Stitching mode

下载图片 查看所有图片

图 7. Stitching direction为1的配准图像

Fig. 7. Registering image with stitching direction of 1

下载图片 查看所有图片

当stitching direction为1和3时，需移动配准图像B'完全到参考图像的坐标系中，移动图像B'的计算公式为

式中：dx、dy分别为图像B'顶点坐标超出X轴坐标系最小的x值和超出Y轴坐标系最小的y值。

2.5　改进图像融合

在显微图像拼接过程中，图像之间的不连续性、曝光不同等因素可能导致拼接缝隙的出现，而加权平均融合算法可以在保持拼接效果的同时尽可能地减少这些缝隙。此外，该算法还能够在融合速度和质量之间找到一个平衡点，可以兼顾融合效果和融合时间。加权平均融合算法通过为重叠区域的像素值分配一定的权重来计算两张图像的像素值。在图像重叠部分，权重会渐变，这种方法简单而有效，可以产生更平滑的拼接图像。如图8（a）所示，加权平均融合算法定义为

式中：I为融合后的图像；I1、I2分别为全景图像和配准图像；a1、a2分别为两张图像融合时重叠区域像素点的权值，且满足a1+a2=1。在图8（a）中有

式中：x为像素点坐标；Xl、Xr分别为重叠区域左边界和右边界。

图 8. 图像融合。（a）两张图像融合；（b）多张图融合

Fig. 8. Image fusion. (a) Fusion of two images; (b) fusion of multiple images

下载图片 查看所有图片

当进行多图拼接时，配准图像I2与参考图像接壤，还与全景图像I1存在重叠部分，如图8（b）所示。此时，传统加权平均融合算法无法兼顾两部分重叠区域。用角度θ表示重叠区域每列像素点：

式中：θ的取值范围为0,π/2，利用三角函数取权重可表示为

且满足

式中：Yu、Yb分别为重叠区域上边界、下边界；且满足a3+a4=1。

利用角度表示像素点，从而通过三角函数权重渐变的特性实现重叠区域边缘平滑过度，有效地提高图像融合的质量。

2.6　多图拼接

在显微图像多图拼接模式中，预设载物台运动轨迹为“弓”字形，如图2所示。从全景图像中截取与待配准图像尺寸相同的图像作为参考图像，将待配准图像与截取的图像进行两两拼接，将其拼接结果按Pi位置粘贴至全景图像如图9所示。

图 9. 多图拼接

Fig. 9. Multi-image stitching

下载图片 查看所有图片

由于显微图像间仅存在平移关系，无旋转与尺度变换关系，在多图拼接模式中，采用帧到帧的拼接模式。stitching direction为2时的算法流程如图10所示，有效避免传统“先采集后拼接”模式中因局部失配造成的影响。求解单应性矩阵中8个未知数至少需要4对优秀特征点对，若不满足则重新采集待配准图像。

图 10. Stitching direction为2的算法流程

Fig. 10. Algorithm flow with stitching direction of 2

下载图片 查看所有图片

Pi表示参考图像与配准图像拼接后的结果粘贴到全景图像的左上角坐标。当stitching direction为2时，Pi可表示为

式中：Pi,x、Pi,y分别为第i次粘贴拼接结果图像左上角坐标的x值、y值；Pi-1,x为第i-1次粘贴拼接结果图像左上角坐标的x值；Wi-1为第i-1次两两拼接图像的宽；w、h分别为参考图像的宽度和高度；Hi-1为第i-1次全景图像的高度。

同理，当stitching direction为1时，Pi可表示为

式中：Wi为第i次两两拼接图像的宽。

当stitching direction为3时，分为两种情况，当载物台从左向右运动再向上运动时，Pi可表示为

式中：Hi为第i次两两拼接图像的高；Pi-1,y为第i-1次粘贴拼接结果图像左上角坐标的y值。当载物台从右向左运动再向上运动时，Pi可表示为

当stitching direction为4时，分为两种情况，当载物台从左向右运动再向下运动时，Pi可表示为

式中：Wi-1为第i-1次最终拼接图像的宽度。当载物台从右向左运动再向下运动时，Pi可表示为

3　分析与讨论

显微图像拼接系统在Windows 10操作系统下，使用OpenCV 4.4.0作为开发环境，Visual Studio 2017作为编程平台。PC配置为Intel（R）Xeon（R）CPU E5-1630 v3处理器，3.70 GHz主频，16 GB内存，Quadro K2200 4GB NVIDIA显卡，ST1000DM003希捷机械硬盘。通过对分辨率为1600×1200的200组脊髓横切片、木本双子叶植物茎横切片以及特征稀疏的植物花粉显微图像进行拼接，与BF和Flann匹配算法比较性能，以及多图拼接结果与AutoStitch软件拼接结果相比较。

3.1　显微图像结构相似度评价指标

以特征点间的相似性作为显微图像结构相似度的评价指标。通过计算图像中某个特征点与其余特征点的欧氏距离来衡量相似度，当欧氏距离小于阈值时，判定为相似特征点。如表1所示，在相同阈值下，显微图像相比自然场景图像具有更多的相似特征点，即更多的相似结构。图11中用菱形标出了与标“+”处特征点相似的特征点。在图11（a）脊髓显微结构中存在大量相似的点结构，这些点结构难以区分，导致特征相似度最高达到21.15%。

图 11. 图像结构相似性。（a）脊髓切片；（b）双子叶切片；（c）自然场景

Fig. 11. Similarity of image structure. (a) Spinal cord; (b) dicotyledons; (c) natural scene

下载图片 查看所有图片

3.2　匹配点预测坐标及实际匹配坐标

图12为显微镜放大倍率为200倍时，电动载物台X轴移动控制下植物花粉成像重叠区域为20%和30%的情况。在这些情况下，3组特征点的坐标预测值和实际匹配值都得到展示。通过判断匹配特征点与预测区域内待匹配特征点间的最小欧氏距离，确定实际匹配点。由于机械精度误差及显微镜设备周围环境的影响，实际匹配坐标与预测匹配坐标存在误差。X轴误差在0.19%~0.25%范围内，Y轴误差在0.16%~0.33%范围内，如表2所示。

图 12. 不同重叠区域预测坐标与实际坐标。（a）20%图像重叠区域；（b）30%图像重叠区域

Fig. 12. Predicted and actual coordinates of different overlapping areas. (a) 20% image overlap area; (b) 30% image overlap area

下载图片 查看所有图片

3.3　所提算法与BF、Flann匹配算法性能比较

选取200组不同类别、大小为1600×1200的显微图像进行实验，分别比较所提匹配算法与BF和Flann匹配算法的性能。以匹配率和匹配耗时作为衡量标准，列出其中3组实验结果如表3所示。当图像特征点较少时三者匹配耗时区别不大，但所提匹配算法获得更多的优秀匹配对。随着图像检测到的特征点数增加，所提匹配算法在匹配率和匹配耗时上优于BF和Flann匹配算法，匹配率提高7.95%~26.52%。

图13为3组不同的显微图像通过3种匹配算法后得到的初匹配和优秀匹配结果。由初匹配结果可知，所提匹配算法有效减弱显微结构相似性带来的误匹配影响。该算法仅在预测区域内筛选待匹配的特征点，缩小了匹配点的范围，减少了特征点匹配过程的计算量，从而减少运算时间。将初匹配结果通过相同的阈值的Low’s算法获得优秀匹配结果，所提匹配算法获得更多的优秀匹配对，提高了求解单应性矩阵的准确度，从而提高图像配准效率。

图 13. 不同匹配算法性能比较实验结果。（a）BF初匹配；（b）Flann初匹配；（c）所提算法初匹配；（d）BF优秀匹配；（e）Flann优秀匹配；（f）所提算法优秀匹配

Fig. 13. Experimental results of performance comparison of different matching algorithms. (a) BF first match; (b) Flann first match; (c) proposed algorithm first match; (d) BF good match; (e) Flann good match; (f) proposed algorithm good match

下载图片 查看所有图片

3.4　图像拼接与融合

显微镜采集图像时因曝光率不同，会造成拼接图像的亮度不一致如图14（a）所示，采用所提的改进融合算法对拼接后的图像进行加权融合，使图像拼接处更平滑自然如图14（b）所示。

图 14. 融合图像。（a）融合图像前；（b）融合图像后

Fig. 14. Fusion images. (a) Before fusing image; (b) after fusing image

下载图片 查看所有图片

实验对比两组放大倍率不同的显微图像分别采取AutoStitch开源软件拼接和所提预测拼接的效果。如图15、图16所示，二者拼接效果相差不大，但AutoStitch拼接图像边缘出现变形。如表4所示，所提算法拼接耗时优于AutoStitch软件，这是因为该软件采用SIFT提取特征点和128维向量描述特征点，增加了特征匹配计算量，以及采取多波段融合拼接图像。

图 15. 3×3图像拼接。（a）拼接原图像；（b）AutoStitch拼接；（c）所提算法拼接

Fig. 15. 3×3 images stitching. (a) Stitched original images; (b) AutoStitch stitching; (c) proposed method stitching

下载图片 查看所有图片

图 16. 4×5图像拼接。（a）AutoStitch拼接；（b）所提算法拼接

Fig. 16. 4×5 images stitching. (a) AutoStitch stitching; (b) proposed method stitching

下载图片 查看所有图片

如图17（a）所示，两套丝杆与步进电机的组合分别控制载物台X轴与Y轴移动距离。如图17（b）所示，在多图拼接模式下，显微镜载物台按预设轨迹运动时显示实时拼接过程部分图像。

图 17. 显微图像拼接系统。（a）电控载物台；（b）（c）实时拼接图像部分过程

Fig. 17. Microscopic image stitching system. (a) Motorized stage; (b) (c) part of the process of stitching images in real time

下载图片 查看所有图片

4　结论

针对显微图像存在大量相似结构引起匹配效率不高，以及传统显微图像拼接系统中图像采集和拼接互相独立，抗干扰性不强等问题，提出了一种结合载物台运动信息的显微图像预测拼接算法。充分利用高精度电动载物台运动信息预测待匹配特征点位置，降低因显微图像相似结构造成的误匹配影响。利用匹配特征点对的斜率快速粗筛选误匹配点对，减少RANSAC算法精匹配迭代次数。最后，通过改进的加权平均融合消除多图拼接缝隙平滑图像。实验对比结果表明：所提算法降低了显微图像特征点误匹配率，有效利用载物台运动信息加速了显微图像拼接过程。当显微图像分辨率为1600×1200且特征丰富时，系统的拼接速度达到2 frame·s^-1，这为进一步研究显微镜实时拼接算法提供依据。