一种显微子图像拼接方法及其大批量图像的拼接方法

1.本发明涉及图像处理领域，具体涉及一种显微子图像拼接方法及其大批量图像的拼接方法。


背景技术：

2.图像拼接是将一组相互重叠的图像序列进行空间匹配对准，经重采样合成形成一幅包含各图像序列信息的宽视角场景的、完整的、高清晰的新图像的技术。该问题涉及的关键技术是如何找到图像的重合部分以及如何将具有重合部分的图像进行拼合，即图像配准和图像融合。图像配准方法可分为基于区域的图像配准和基于特征的图像配准。前者比较两幅图像不同区域的统计信息实现配准，有互相关、互信息、交叉熵等方法，实现简单，不需要提取额外特征，但是计算量比较大，鲁棒性较差，适用场景较少；后者从待配准的图像中提取特定的特征集作为控制结构，利用其之间的对应关系来进行配准，主要有harris、fast、susan、sift、surf等算法，它们在解决图像尺度不变、透视不变性具有一定的优势，并且具有一定的鲁棒性和稳定性，得到广泛应用。
3.然而现有算法具有以下几个缺点：
4.1、图像配准算法经过几十年的发展已经取得了很大的进展，但由于拍摄环境复杂多变，现在还没有一种算法能够解决所有图像的匹配问题。现在的几种方法各有其优缺点，如果能综合利用这些方法的优点将会取得更好的匹配结果。
5.2、现有的图像拼接将研究重点放在了两张图像间的配准，忽略了大视野图像的整体拼接：在面临整体拼接时，采用事先设定的顺序直接与上一张图像进行拼接。这种传统的拼接方法具有简单快速的特点，可以完成粗略的图像拼接；但随着图像数量增加，容易因为累积误差产生拼接效果不佳的情况。
6.现有技术中公开了一种基于改进的harris角点检测拼接超声图像方法的专利，该方法包括s1：图像获取：在ct模拟定位室及治疗室中，先后采集ct及超声图像，本研究所有的图像数据来自三名腹部术后肿瘤患者扫描所获取的体数据，每位患者选取10个超声及cbct的层面图像,s2对s1中获取的图像进行预处理：在做超声图像拼接之前首先将400
×
400的二维超声图像分割成无重叠的4
×
4子图，计算子图中角点响应函数r值分布，为获取自适应阈值做准备,实现了两幅超声拼接图像间的平滑过渡，消除了拼接缝隙，图像拼接的准确性优于传统的harris角点检测的方法。然而上述方法没有结合图片的特征点以及特征向量，图片的拼接准度不够，另外该方法采用传统的拼接路径，即直接与上一张图像进行拼接，无法完成大批量图像的拼接，且容易出现误差累积效应，发生畸变。


技术实现要素：

7.本发明提供一种显微子图像拼接方法，解决显微图像快速且准确通过特征匹配进行两两拼接的问题。
8.本发明又一目的在于提供一种大批量显微子图像拼接方法，解决大批量显微图像
拼接因累计误差与特征匹配不佳而导致拼接效果不佳的问题。
9.为了达到上述技术效果，本发明的技术方案如下：
10.一种显微子图像拼接方法，包括以下步骤：
11.s1：获取原始显微子图像数据集，对原始显微子图像数据做预处理；
12.s2：根据步骤s1所述的已预处理的原始显微子图像数据集，对原始显微子图像进行两两拼接，其中两两拼接包括以下步骤：
13.s21：分别提取两幅图像的特征点以及特征点其所对应的特征向量；
14.s22：根据步骤s21所得的特征点以及特征向量对两幅图像进行图像配准；
15.s23：根据s22所得的配准信息对两幅图像进行图像融合拼接；
16.本发明利用图像的特征点以及特征点其所对应的特征向量，能够为后续图像配准以及图像融合拼接提供有效特征信息，实现更快的特征匹配速度；以两两拼接为基本拼接环节，结合树型拼接路径完成大批量原始显微子图像的拼接，能够快速完整复原显微图像。
17.进一步地，所述步骤s1中，所述对原始显微子图像数据做预处理为对原始显微子图像进行两倍下采样，得到尺寸为原图像四分之一的图像。
18.针对大尺寸图像，直接提取特征点，占用计算空间大且用时长，因而需要对图像进行下采样。
19.进一步地，所述步骤s21中，所述提取图像的特征点，包括以下步骤：
20.s211：设单张图像中的像素点坐标为(x,y)，对坐标为(x,y)的像素点在x方向上进行二阶偏导得到在y方向上进行二阶偏导得到在x方向上进行一阶偏导后在y方向上再进行一阶偏导得到i
xy
(x,y)，利用所得到的和i
xy
(x,y)，构建自相关矩阵m(x,y)：
[0021][0022]
根据自相关矩阵m(x,y)，结合角点响应函数计算得到每一个像素点的r值，其中角点响应函数为
[0023][0024]
det(m)表示矩阵m(x,y)的行列式，tr(m)表示矩阵m对角线之和，α和β为矩阵m的特征值，分别代表x方向和y方向的梯度，det(m)＝α
·
β，tr(m)＝α+β，ε为一极小里量；
[0025]
s212：图像按照宽n等分，高m等分，平均分成n
×
m块，对于每一块图像提取该块中所有像素点坐标(x,y)，并按r值的大小,将像素点坐标(x,y)进行从大到小排序，根据排序结果，保留该块图像r值在前80％的像素点作为候选角点，其余的像素点剔除；
[0026]
s213：选用一个k
×
k的窗口矩阵，按照步长为1，沿x,y的轴向顺序处理整幅图像，其中k值设置不能大于图像的长和宽，对于窗口矩阵所框选的图像部分，在该图像中保留框选范围内r值最大的候选角点，并将其余的候选角点剔除，完成整幅图像的角点检测，并将所检测出的角点作为图像的特征点；
[0027]
s214：根据步骤ⅰ所得的特征点的邻域像素，计算该点的梯度方向作为特征点向，其梯度计算公式为：
[0028][0029][0030]
其中，l为该特征点的尺度空间值，(x,y)为该特征点坐标，m(x,y)表示该特征点的梯度幅值，θ(x,y)表示该特征点的梯度方向。
[0031]
步骤s211参照经典harris角点检测算法的思想，同时经典harris角点检测算法采用的经典角点响应函数为：
[0032]
r＝det(m)-ktr2(m)
[0033]
角点位置梯度α、β都较大，r值随着也较大，即图像分块的目的除为了实现并行处理之外，还能够将子块分次调入内存处理，增强算法的适用性，避免过大图像处理不了的问题，图像分块方式一般存在两种形式：一种是已知分块的数量n*m份，即将原图像的宽高分别分成n和m等份；另一种就是已知子块的宽高，然后从原点开始切分。
[0034]
所用窗口矩阵为最大值滤波器，保留r值最大的点，一定程度上减少角点的数目且使角点分布的更加均匀,有效避免角点集聚现象的产生。
[0035]
k值为常数项，范围在0.04-0.06，但k值的选取易造成所检测角点的点位置偏差，本发明步骤s211采用的角点响应函数，避免了k值的选取，能够减小所检测角点的点位置偏差，ε可防止分母为零时出现异常。
[0036]
进一步地，所述步骤s21中，所述提取图像的特征点其所对应的特征向量，包括以下步骤：
[0037]
s215：将坐标轴按照特征点的特征方向进行旋转，再将特征点附近的邻域划分为4
×
4个子区域，在每个子区域内计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值，形成一个种子点，每个子区域的梯度方向直方图将0到360度划分为8个方向范围，每个范围为45度，这样，每个种子点共有8个方向的梯度强度信息，那么一共有4
×4×
8＝128个数据，最终每个角点均对应形成一个128维的sift特征向量h＝(h1,h2,
…
,hi,
…
,h
128
),hi表示一个数据；
[0038]
s216：对所得到的特征向量的长度做归一化处理得到归一化向量：
[0039]
。
[0040]
为了保证特征向量具有旋转不变性，需要以特征点为中心，将特征点附近邻域内图像梯度的位置和方向旋转一个方向角θ，即将原图像x轴转到与主方向相同的方向。
[0041]
在求出4*4*8的128维特征向量后，此时sift特征向量已经去除了尺度变化、旋转等几何变形因素的影响。而图像的对比度变化相当于每个像素点乘上一个因子，光照变化是每个像素点加上一个值，但这些对图像归一化的梯度没有影响。因此将特征向量的长度归一化，则可以进一步去除光照变化的影响。
[0042]
进一步地，所述对两幅图像进行图像配准，包括以下步骤：
[0043]
s221：设x＝{x1,x2,
…
,xi,
…
}为图像1的特征向量集合，其中xi对应特征点的特征向量y＝{y1,y2,
…
,yj,
…
}为图像2的特征点集合，其中yj对应特征点的特
征向量i＝{
…
,(xi,yj),
…
}为匹配点集合，特征点xi和特征点yj的距离表示为d(xi,yj)，采用欧氏距离方法结合特征点的特征向量计算：
[0044][0045]
利用双向-k最近邻匹配方法进行特征点的筛选与匹配；
[0046]
s222：采用随机抽样一致算法ransac再次筛选并去除匹配点集合i1的错误匹配点，完成最终的匹配点集合i；
[0047]
双向-k最近邻匹配方法，对匹配点进行二次筛选，消除误匹配。
[0048]
进一步地，所述双向-k最近邻匹配方法，包括以下步骤：
[0049]
s2211：设正向匹配阈值为r1，其中r1取值范围为0.5到0.8，反向匹配阈值为r2，其中r2＝0.9，对于图像1的某个特征点xi，利用欧氏距离方法，将图像2对应的所有特征点均与特征点xi计算距离大小，查找到与xi距离最近的两个特征点y
k1
和特征点y
k2
，若则我们将(xi,y
k1
)视为待定的匹配点对，加入匹配点候选集合i0；
[0050]
s2212：对于图像2的某个特征点yj，利用欧氏距离方法，将图像1对应的所有特征点均与特征点yj计算距离大小，查找到与yj距离最近的两个特征点x
k1
和特征点x
k2
，若且该点对(x
k1
,yj)属于匹配点候选集合i0，则我们将其视为正确的匹配点对，加入匹配点集合i1。
[0051]
进一步地，所述采用随机抽样一致算法ransac的计算步骤如下：
[0052]
s2221：从匹配点集合i1中随机抽出不共线的4个样本点作为内点并加入内点集li，计算出变换矩阵hi，记为模型mi；
[0053]
s2222：将匹配点集合i1中的所有剩余点代入模型mi中，若计算结果满足模型mi，则视为内点并加入内点集li；
[0054]
s2223：记录当前内点集li元素个数，并作为最优内点集l，继续迭代重复步骤s2221和s2222得到新的模型mj和内点集lj,若新的内点集lj元素个数大于最优内点集l元素个数，则将内点集lj作为新的最优内点集l，并将模型更新为模型mj,并更新迭代次数k；若新的内点集lj元素个数小于等于最优内点集l元素个数，则维持最优内点集l不变；
[0055]
s2224：待迭代次数k完成后将最优内点集l作为最终的匹配点集合i。
[0056]
进一步地，所述步骤s2223中更新迭代次数k采用以下公式：
[0057][0058]
其中z为所需成功率，范围在0到1之间，k是求解模型需要最少的点的个数，与最初的内点数一致；p是当前内点集内点数量与总匹配点数量的比值。
[0059]
迭代次数初始设定为无穷大，每次更新模型参数估计时，用当前的“内点”比值当成p(内点的概率)来估算出迭代次数。
[0060]
每次计算模型使用n个点的情况下，选取的点至少有一个外点的情况(采样失败的概率)就是：1-pk；
[0061]
那么能采样到正确的k个点去计算出正确模型的概率是:
[0062]
z＝1-(1-pk)k[0063]
z是所需成功率设置，比如想要达到99％的成功率，则z＝0.99。
[0064]
进一步地，所述对两幅图像进行图像融合拼接，包括以下步骤：
[0065]
s231：根据匹配点集合i，设图像1的匹配点坐标为(x,y)，则对应图像2的匹配点坐标为(x
′
,y
′
),根据仿射变换求出变换矩阵h，公式如下：
[0066][0067]
其中a
11
,a
12
,a
21
,a
22
是尺度、旋转角度和裁切角度，a
x
和ay是图像变换参数；
[0068]
s232：利用所求的变换矩阵h将图像1映射到图像2所在的坐标系中，此时图像1和图像2位于同一坐标系下，存在重叠区域，根据重叠区域的对应坐标采用加权平均融合方法得到最终拼接图像，其中加权平均融合方法采用以下公式：
[0069][0070]
式中：f1(x,y)和f2(x,y)表示待拼接的两幅图像，f(x,y)表示拼接完成的图像，w1和w2表示权值，且xi为当前像素点的横坐标，x
l
为重叠区域的左边界，xr为重叠区域的右边界，运算过程中w1由1变到0，w2由0变到1，从而实现两幅图像的平滑过渡。
[0071]
显微镜拍摄图像时，通过移动工作台，从左到右，从上到下，按照顺序拍摄，最终得到的图像集也是有序的，每一行拍摄到的图像为一组，有多组。由于显微镜的图像只有平移没有旋转，对平移图像融合效果好的、速度快的加权平均融合算法
[0072]
一种大批量显微子图像拼接方法，其特征在于，以两两拼接作为基本拼接环节，结合树型拼接路径对大批量原始显微子图像进行拼接复原原始显微图像，所采用树型拼接路径拼接是横向路径拼接与纵向路径拼接的交替进行，横向路径拼接是对每一行不同的显微子图像进行两两拼接，纵向路径拼接是对每一列的不同显微子图像进行两两拼接，每一次的横向路径拼接或者纵向路径拼接中的若干次拼接设置为并行处理，每一次的横向路径拼接或者纵向路径拼接结束显微子图像的尺寸均同步增大，通过重复横向拼接路径与纵向拼接路径的交替作用将原始显微图像完整复原。
[0073]
与现有技术相比，本发明技术方案的有益效果是：
[0074]
一种显微子图像拼接方法利用图像所检测的角点作为图像的特征点，进而得到特征点的特征向量，为后续图像配准以及图像融合拼接提供有效特征信息，实现更快的特征匹配速度；采用双向-k最近邻匹配，能有效减少错配点，使得到的配准图更加准确，提高配准精度；一种大批量显微子图像拼接方法以两两拼接为基本拼接环节，结合树型拼接路径完成大批量原始显微子图像的拼接，能够快速完整复原显微图像。
附图说明
[0075]
图1为一种显微子图像拼接方法的流程图；
[0076]
图2为实施例1所使用的大批量原始显微子图像集ⅰ；
[0077]
图3为实施例1采用树型路径进行拼接的过程示意图；
[0078]
图4为实施例1使用一种大批量显微子图像拼接方法所得到的原始显微图像；
[0079]
图5为实施例2所使用的大批量原始显微子图像集ⅱ；
[0080]
图6为实施例2采用树型路径进行拼接的过程示意图；
[0081]
图7为实施例2使用一种大批量显微子图像拼接方法所得到的原始显微图像；
[0082]
图8为实施例3所使用的大批量原始显微子图像集ⅲ；
[0083]
图9为实施例3采用树型路径进行拼接的过程示意图；
[0084]
图10为实施例3使用一种大批量显微子图像拼接方法所得到的原始显微图像。
具体实施方式
[0085]
附图仅用于示例性说明，不能理解为对本专利的限制；
[0086]
对于本领域技术人员来说，附图中某些公知说明可能省略是可以理解的。
[0087]
下面结合附图和实施例对本发明的技术方案做进一步的说明。
[0088]
实施例1：
[0089]
本实施例使用大批量原始显微子图像数据集ⅰ，如图2所示，共4张
×
4组。
[0090]
一种大批量显微子图像拼接方法，如图1所示，包括以下步骤：
[0091]
s1：获取原始显微子图像数据集，对原始显微子图像进行两倍下采样，得到尺寸为原图像四分之一的图像；
[0092]
s2：根据步骤s1所述的已预处理的原始显微子图像数据集，对原始显微子图像进行两两拼接，其中两两拼接包括以下步骤：
[0093]
s21：分别提取两幅图像的特征点以及特征点其所对应的特征向量；
[0094]
s211：设单张图像中的像素点坐标为(x,y)，对坐标为(x,y)的像素点在x方向上进行二阶偏导得到在y方向上进行二阶偏导得到在x方向上进行一阶偏导后在y方向上再进行一阶偏导得到i
xy
(x,y)，利用所得到的和i
xy
(x,y)，构建自相关矩阵m(x,y)：
[0095][0096]
根据自相关矩阵m(x,y)，结合改进的角点响应函数计算得到每一个像素点的r值，其中改进的角点响应函数为
[0097][0098]
det(m)表示矩阵m(x,y)的行列式，tr(m)表示矩阵m对角线之和，α和β为矩阵m的特征值，分别代表x方向和y方向的梯度，det(m)＝α
·
β，tr(m)＝α+β，ε为一极小里量；
[0099]
s212：图像按照宽n等分，高m等分，平均分成n
×
m块，对于每一块图像提取该块中所有像素点坐标(x,y)，并按r值的大小,将像素点坐标(x,y)进行从大到小排序，根据排序结果，保留该块图像r值在前80％的像素点作为候选角点，其余的像素点剔除；
[0100]
s213：选用一个k
×
k的窗口矩阵，按照步长为1，沿x,y的轴向顺序处理整幅图像，其中k值设置不能大于图像的长和宽，对于窗口矩阵所框选的图像部分，在该图像中保留框
选范围内r值最大的候选角点，并将其余的候选角点剔除，完成整幅图像的角点检测，并将所检测出的角点作为图像的特征点；
[0101]
s214：根据步骤ⅰ所得的特征点的邻域像素，计算该点的梯度方向作为特征点向，其梯度计算公式为：
[0102][0103][0104]
其中，l为该特征点的尺度空间值，(x,y)为该特征点坐标，m(x,y)表示该特征点的梯度幅值，θ(x,y)表示该特征点的梯度方向。
[0105]
s215：将坐标轴按照特征点的特征方向进行旋转，再将特征点附近的邻域划分为4
×
4个子区域，在每个子区域内计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值，形成一个种子点，每个子区域的梯度方向直方图将0到360度划分为8个方向范围，每个范围为45度，这样，每个种子点共有8个方向的梯度强度信息，那么一共有4
×4×
8＝128个数据，最终每个角点均对应形成一个128维的sift特征向量h＝(h1,h2,
…
,hi,
…
,h
128
),hi表示一个数据；
[0106]
s216：对所得到的特征向量的长度做归一化处理得到归一化向量：
[0107]
。
[0108]
s22：根据步骤s21所得的特征点以及特征向量对两幅图像进行图像配准；
[0109]
s221：设x＝{x1,x2,
…
,xi,
…
}为图像1的特征向量集合，其中xi对应特征点的特征向量y＝{y1,y2,
…
,yj,
…
}为图像2的特征点集合，其中yj对应特征点的特征向量i＝{
…
,(xi,yj),
…
}为匹配点集合，特征点xi和特征点yj的距离表示为d(xi,yj)，采用欧氏距离方法结合特征点的特征向量计算：
[0110][0111]
s2211：设正向匹配阈值为r1，其中r1取值范围为0.5到0.8，反向匹配阈值为r2，其中r2＝0.9，对于图像1的某个特征点xi，利用欧氏距离方法，将图像2对应的所有特征点均与特征点xi计算距离大小，查找到与xi距离最近的两个特征点y
k1
和特征点y
k2
，若则我们将(xi,y
k1
)视为待定的匹配点对，加入匹配点候选集合i0；
[0112]
s2212：对于图像2的某个特征点yj，利用欧氏距离方法，将图像1对应的所有特征点均与特征点yj计算距离大小，查找到与yj距离最近的两个特征点x
k1
和特征点x
k2
，若且该点对(x
k1
,yj)属于匹配点候选集合i0，则我们将其视为正确的匹配点对，加入匹配点集合i1。
[0113]
s222：采用随机抽样一致算法ransac再次筛选并去除匹配点集合i1的错误匹配点，完成最终的匹配点集合i；
[0114]
s2221：从匹配点集合i1中随机抽出不共线的4个样本点作为内点并加入内点集li，计算出变换矩阵hi，记为模型mi；
[0115]
s2222：将匹配点集合i1中的所有剩余点代入模型mi中，若计算结果满足模型mi，则视为内点并加入内点集li；
[0116]
s2223：记录当前内点集li元素个数，并作为最优内点集l，继续迭代重复步骤s2221和s2222得到新的模型mj和内点集lj,若新的内点集lj元素个数大于最优内点集l元素个数，则将内点集lj作为新的最优内点集l，并将模型更新为模型mj,并更新迭代次数k；若新的内点集lj元素个数小于等于最优内点集l元素个数，则维持最优内点集l不变；其中更新迭代次数k采用以下公式：
[0117][0118]
其中z为所需成功率，范围在0到1之间，k是求解模型需要最少的点的个数，与最初的内点数一致；p是当前内点集内点数量与总匹配点数量的比值。
[0119]
s2224：待迭代次数k完成后将最优内点集l作为最终的匹配点集合i。
[0120]
s23：根据s22所得的配准信息对两幅图像进行图像融合拼接；
[0121]
s231：根据匹配点集合i，设图像1的匹配点坐标为(x,y)，则对应图像2的匹配点坐标为(x
′
,y
′
),根据仿射变换求出变换矩阵h，公式如下：
[0122][0123]
其中a
11
,a
12
,a
21
,a
22
是尺度、旋转角度和裁切角度，a
x
和ay是图像变换参数；
[0124]
s232：利用所求的变换矩阵h将图像1映射到图像2所在的坐标系中，此时图像1和图像2位于同一坐标系下，存在重叠区域，根据重叠区域的对应坐标采用加权平均融合方法得到最终拼接图像，其中加权平均融合方法采用以下公式：
[0125][0126]
式中：f1(x,y)和f2(x,y)表示待拼接的两幅图像，f(x,y)表示拼接完成的图像，w1和w2表示权值，且xi为当前像素点的横坐标，x
l
为重叠区域的左边界，xr为重叠区域的右边界，运算过程中w1由1变到0，w2由0变到1，从而实现两幅图像的平滑过渡。
[0127]
s3：以步骤s2所述的两两拼接作为基本拼接环节，结合树型拼接路径对大批量原始显微子图像进行拼接复原原始显微图像，具体过程如下，如图3所示：
[0128]
横向路径进行两两拼接一次，纵向路径进行两两拼接一次，横向路径进行两两拼接一次，纵向路径进行两两拼接一次，最终将原始显微图像完整复原，如图4所示。
[0129]
实施例2：
[0130]
本实施例使用大批量原始显微子图像数据集ⅱ，，如图5所示，共5张
×
4组。
[0131]
步骤s1和s2与实施例1相同：
[0132]
s3：以步骤s2所述的两两拼接作为基本拼接环节，结合树型拼接路径对大批量原始显微子图像进行拼接复原原始显微图像，具体过程如下，如图6所示：
[0133]
横向路径进行两两拼接一次，纵向路径进行两两拼接一次，纵向路径进行两两拼接一次，横向路径进行两两拼接一次，横向路径进行两两拼接一次，最终将原始显微图像完整复原，如图7所示。
[0134]
实施例3：
[0135]
本实施例使用大批量原始显微子图像数据集ⅲ，，如图8所示，共5张
×
4组。
[0136]
步骤s1和s2与实施例1相同：
[0137]
s3：以步骤s2所述的两两拼接作为基本拼接环节，结合树型拼接路径对大批量原始显微子图像进行拼接复原原始显微图像，具体过程如下，如图9所示：
[0138]
横向路径进行两两拼接一次，横向路径进行两两拼接一次，纵向路径进行两两拼接一次，纵向路径进行两两拼接一次，横向路径进行两两拼接一次，最终将原始显微图像完整复原，如图10所示。
[0139]
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。