一种多图自动配准及融合拼接方法与流程

本发明涉及机器视觉、图像处理技术领域，尤其是涉及一种多图自动配准及融合拼接方法；

背景技术：

目前，大部分医疗成像设备拍摄的图像都是局部图片，视野受限，全景图像设备价格昂贵，难以推广使用；利用图像拼接技术生成全景图片，成本低、效果好，提供整体目标图像，便于医务人员更好地进行诊断和治疗，在X射线成像、超声检查等方面均可以使用；近年来，拼接技术在国内外取得很大进展，常用配准方法包括变换域匹配、灰度匹配、特征匹配、模板匹配等，通过计算图像间变换参数，将多幅待拼接图像图像变换坐标系，实现拼接；由于医学成像受电子噪声以及射线本身易散射，拍摄位置差异、光照和人体组织固有的复杂因素影响，往往造成图像细节模糊、分辨率低、亮度暗、对比度差的现象出现，常用的配准拼接方法效果并不理想，且多图拼接通常需要人工参与选择，特别是对于医学显微图像，甚至由上百幅图像拼接成一个目标，数据量、工作量极大，且方向、重叠区域位置很难把握，因此，针对医学图像的特殊性，目前的拼接方法很难满足。

本发明针对医学图像的特殊性，提出一种多图自动配准及融合拼接方法，方法通过多尺度Retinex方法增强多幅图像，之后利用改进的相位相关互功率谱快速估计重叠区域，减少配准范围的同时通过冲激函数能量峰值大小实现多图自动排序，接着通过SUFR局部特征提取算法提取待配准区域特征点，改进随机抽样一致算法精确提纯匹配点对，优化模型参数，建立图像变换矩阵实现拼接，最后针对接缝明显的现象，利用NSCT变换算法并制定融合策略进一步对图像融合处理，使得图像接缝过渡平滑自然。

技术实现要素：

为了解决由于医学成像的特殊性，现有配准及融合拼接方法存在的不足的问题，本发明提供一种多图自动配准及融合拼接方法，该方法很大程度改善了因多图手工排序困难、图像细节模糊、分辨率低、噪声高、照度不均等问题对图像配准和拼接的影响，对节省成像设备开销以及医务人员准确诊断提供了帮助，与传统相比有明显优势，在辅助医学成像领域有着良好的应用前景。

本发明解决其技术问题所采用的技术方案是：

一种多图自动配准及融合拼接方法，其方法步骤为:

1.输入未排序的图像序列，通过多尺度Retinex方法增强图像细节，提高对比度；

2.接着图像间用对数极坐标变换改进的相位相关算法互匹配，计算缩放、旋转和平移参数，粗估计图像间重叠区域，并根据冲激函数能量峰值大小制定排序规则,实现多图自动排序；

3.然后通过SUFR算法提取待配准重叠区域特征点，改进随机抽样一致算法精确提纯匹配点对，优化模型参数，建立图像间变换矩阵，依次拼接序列图像；

4.最后利用NSCT变换算法并制定融合策略，通过多尺度分解后的低频、高频子带进一步对拼接图像进行图像融合处理，输出融合后的拼接图像。

本发明有益效果：一种多图配准及拼接融合方法，该方法很大程度改善了因电子噪声和射线本身易散射造成的图像细节模糊、分辨率低、噪声高、照度不均等问题对图像配准和拼接的影响，解决拍摄位置差异问题，对节省成像设备开销以及医务人员准确诊断提供了帮助，与传统相比有明显优势，在辅助医学成像领域有着良好的应用前景。

附图说明

图1为本发明Retinex方法增强前后比较图；其中:(a)原图像；(b)原图直方图；(c)增强后图像；(d)增强后图像直方图。

图2为本发明冲激函数能量分布图；其中:图(a)和(b)为两幅不同位置人体肺部图像，图(c)为图(b)增加噪声、改变照度等处理的结果，图(d)为图(a)和图(b)对数极坐标变换改进的相位相关法冲激函数能量分布情况；图(e)是图(a)和图(c)冲激函数能量分布情。

图3为本发明精确提纯处理匹配图；其中:(a)阈值0.4匹配；(b)阈值0.7匹配；(c)提纯后匹配；(d)光照改变匹配；(e)旋转缩放匹配；(f)模糊匹配。

图4为本发明胯骨自动排序前序列图。

图5为本发明胯骨自动排序后序列图。

图6为融合拼接比较图；其中:(a)直接拼接；(b)平均梯度融合拼接；(c)小波变换融合拼接；(d)本发明拼接。

具体实施方式

一种多图自动配准及融合拼接方法，其具体实施步骤如下：

(1)输入未排序的图像序列，通过多尺度Retinex方法增强图像细节，提高对比度；包括如下步骤：

(1.1)将原始图像I(x,y)中像素点分解为两个独立的部分：周围环境亮度决定的入射照度分量L(x,y)和物体表面反射率决定的反射分量R(x,y)，原始图像表示为：I(x,y)＝L(x,y)*R(x,y)；

其中L(x,y)决定像素的动态范围，反映缓慢变化的低频信息，R(x,y)决定像素的内在性质，反映高频细节信息；

(1.2)利用对数分离入射、反射分量，单尺度Retinex的公式表示如下：

上式中(x,y)是像素点的坐标，第i个通道中像素值为I_i(x,y)，其中i＝0,1,2...N，通道个数为N，N＝1表示灰度图像，N＝3表示彩色图像，i取1,2,3分别表示R、G、B三通道，*是卷积符号，高斯环绕函数其中r＝x²+y²，λ为常数，c是环绕尺度,需满足条件∫∫F(x,y)dxdy＝1，其中环绕尺度c值若小，则高斯模板半径小，图像细节增强明显，动态压缩范围大，但颜色失真明显；c值若大，颜色丰富，真实性好，但动态压缩范围和局部细节减弱，多尺度Retinex可以更好地兼顾了色彩保真度、动态压缩范围、图像局部细节等方面的一致性，并消除中低尺度光晕化问题；

(1.3)多尺度Retinex表示为单尺度Retinex输出结果的加权求和，多尺度Retinex公式表示如下：

其中K表示环绕尺度的个数，高斯函数权值为W_k，K＝1时即退化为单尺度Retinex，此时W₁＝1；本发明将高、中、低三个尺度合并，取K＝3，W₁＝W₂＝W₃＝1₃；本发明环绕尺度c的高、中、低三个尺度取值分别为120、80、40；

(1.4)将r_i多(x,y)从对数域向实数域转换可得反射图像R_i(x,y),若图像为彩色图像，则分别在R、G、B三通道计算后合成，若图像为灰度图像，则对灰度值按以上步骤计算即可。

(2)图像间用对数极坐标变换改进的相位相关算法互匹配，计算缩放、旋转和平移参数，粗估计图像间重叠区域，并根据冲激函数能量峰值大小制定排序规则，实现多图自动排序；包括如下步骤：

(2.1)利用对数极坐标变换改进的相位相关算法互匹配图像，可以解决原相位相关算法不适用于图像发生旋转、缩放情况，计算出缩放、旋转和平移参数，提高不同环境的鲁棒性并粗估计出图像间重叠区域，具体步骤如下：

a.设图像I₁(x,y),I₂(x,y)傅里叶变换后为I₁′(u,v),I₂′(u,v)，当I₁，I₂存在平移、旋转和缩放位移时，设λ为缩放比例，θ₀为旋转角度，Δx和Δy为平移距离；

b.将图像所在空间变换到对数极坐标空间，之后可以将图像间的缩放比例λ和旋转角度θ₀表示为对数极坐标下的加减运算，如式(3)所示，其中η₀＝logλ，M₁和M₂分别为I₁′和I₂′的模；

M₂(η,θ)＝M₁(η-η₀,θ-θ₀) (3)

d.对上式(3)运用相位相关算法，由公式(4)计算图像互功率谱，将右边反傅里叶变换，得到冲激函数，其中M^*₂(u,v)为M₂(u,v)复共轭，冲激函数峰值的位置即表示缩放比例和旋转角度，将I₂按照求取的缩放比例和旋转角度的大小反变换后，和I₁再次进行相位相关法计算，得到平移距离(Δx，Δy)；

(2.2)根据相位相关算法及其冲激函数能量峰值大小制定排序规则，实现多图自动排序，解决手动排序易错、耗时的问题，具体步骤如下：

a.根据上步计算两两图像的互功率谱，反傅里叶变换后求得的冲激函数峰值的大小表示重叠区域的多少，其值范围是[0,1]，将其作为相关系数，N幅图则有N-1个相关系数；

b.将图像间的相关系数存入N*N大小的数组中，由于仅求图像间的相关系数，因此对角线不放数值；首尾两幅图各自有一个相邻图像，中间图有两个相邻图像，因此从行列两个方向分别选取相关系数最大的两个，图像均至少与其中一幅相邻，行列求交后对应两个值，则为中间图像，求交后一个值的为首尾图像；

c.对于首尾图像的区分可以根据图像平移量，如平移量为正，则为尾图像，如平移量为负，则为首图像；

d.根据首图可以确定和其相邻的下一幅图，继而逐步确定剩余图像的排序，平移量为负的图像排在右边，平移量为正的图像排在左边。

(3)对之前粗估计的图像间重叠区域通过SUFR局部特征提取算法提取待配准区域特征点，改进随机抽样一致算法，精确提纯匹配点对，优化模型参数，根据参考图像建立图像间目标模型变换矩阵，依次拼接序列图像；包括如下步骤：

(3.1)构建尺度空间，求解局部特征点：使用Hessian矩阵行列式在尺度空间寻找局部极大值，建立9×9尺度空间大小的盒式滤波模板对图像滤波，简化高斯二阶微分模板；滤波后的Hessian矩阵行列式近似为:

det(H)＝D_xx×D_yy-(0.9×D_xy)² (5)

D_xx、D_xy、D_yy为滤波模板和图像卷积后的结果；选择不同尺度Hessian检测算子的盒式滤波器构建多尺度空间金字塔图像，寻找不同尺度的特征点；

(3.2)计算主方向：特征点作中心，设s是特征点所在的尺度值，计算半径6s邻域内的点在x和y方向的Haar小波响应，同时设置高斯权重系数；为了求取主方向，以步长0.2弧度左右旋转π/3范围扇形滑动窗口，对范围内的Haar小波响应值累加，得到新的矢量，整个圆形区域遍历后主方向为最长矢量对应的方向，即最大Haar小波响应累加值对应的方向；

(3.3)生成特征描述算子：将坐标轴旋转至主方向，围绕特征点，取边长为20s的方形区域，划分成4×4的子窗口，每个子窗口有5s×5s个像元，利用尺度是2s的小波对子窗口图像计算响应值，采样后得到垂直主方向的dx和沿主方向的dy，以特征点作为中心进行高斯加权计算，计算各个子窗口的响应值，共4×4个子窗口，特征描述算子形成4×4×4＝64维特征矢量后再对向量进行归一化处理；

(3.4)特征匹配：两幅图中计算I₁图中某特征点和I₂图所有特征点的欧氏距离，当最近和次近距离比小于设定阈值0.6时，则特征点匹配；

(3.5)特征点提纯：具体步骤如下：

a.从参考图像I₁中选取匹配点所在坐标的最大、最小值，将包含匹配特征点区域平分成M*M块，本发明取M＝6，去掉块中不包含匹配点的块；

b.随机选取I₁中9个不同块，并从9个块中各随机选取一个匹配点，和与之对应的变换图像I₂中的9个相应匹配点组成样本集，由其中8对点作为模型参数计算临时模型基本矩阵h；如下式所示，其中h₁₁，…，h₃₂为8个自由度参数；

c.检测第9对匹配点是否是否为临时模型基本矩阵的支撑集，若不是，则重新选择9对匹配点，重复b.步骤；若是，则将临时模型矩阵h认为是候选模型矩阵h，并检测所有匹配点是否在候选模型矩阵h上；

d.设共有P对在候选模型矩阵h上，即支撑集为P，检查P是否大于指定阈值，若是，则得到目标模型矩阵h，若否，则转到b.步骤，重新选择9对匹配特征点；

e.通过P对匹配特征点优化目标模型矩阵h的参数，设P对匹配特征点中每一对(x_i,y_i)和对应点(x_i′,y_i′)之间的误差为e，e＝I₂(x_i′,y_i′)-I₁(x_i,y_i)，迭代计算匹配点对的误差函数E的最小值，更新矩阵h，其中当E小于指定阈值时，如阈值＝0.08，得到最终的目标模型矩阵H，根据H变换图像I₂完成配准；

f.根据参考图像建立图像间变换矩阵，依次拼接序列图像；

上述步骤中，从不同的匹配点块中选择匹配点，可以使得构建矩阵的匹配点分布更为均匀，得到的矩阵更加稳定；虽然随机抽样一致性算法本身确定模型矩阵参数时，所需匹配点对和样本集合中的点对数量相同，即只需要8个匹配点即可，但改进方法中，多选第9个匹配点方便验证是否属于临时模型，可以及时重新选择随机样本集，而不用像原方法那样每选一个随机样本集都要寻找其对应支撑集，因此新方法大大节省了时间开销，提高了目标模型建立的效率。

(4)利用NSCT变换算法并制定融合策略，通过多尺度分解后的低频、高频子带进一步对拼接图像进行图像融合处理，输出融合后的拼接图像；

由于被拼接图像照度不均等问题普遍存在，直接拼接易出现明显接缝，NSCT算法由非采样金字塔滤波器组和非采样方向滤波器组构成，在保持平移不变性的同时具有良好的方向选择性、时频局域性和多分辨率的特性，将其用于拼接图像融合中对接缝、伪影有较好的消除效果，可以有效保留边缘和细节，使得融合后的图像清晰，细节丰富，融合过程及规则定义包括如下步骤：

(4.1)设步骤(3)中图I₁和I₂经矩阵H变换后的重叠区域为A(x,y)和B(x,y)，对其分别做NSCT变换，进行多尺度分解，得到系列子带和其中l＝1,2,3...,L，L是分解的最大层数，k＝1,2,3...,2ⁿ，2ⁿ是每层分解的方向数，分解后的低频子带为

(4.2)对于低频子带，图像变化较为缓慢，采用平均融合方法，

(4.3)对于高频子带，反映了人眼比较敏感的边缘、细节特征，采用基于区域方差显著性的高频子带加权融合方法，令区域方差匹配度为M，匹配度阈值为T，一般取0.5≤T≤1，融合权值为w，其中图像A和B的区域方差显著性表示为G(A)和G(B)；

若匹配度M大于或等于阈值T，则子带和根据G(A)、G(B)的大小采用加权平均融合策略如下式所示；

若匹配度M小于阈值T，则选择区域方差显著性大的高频子带输出，即G(A)≥G(B)时，G(A)＜G(B)时，

(4.4)重构图像，输出融合后的拼接图像。

案列分析：

1.图1(a)-(d)依次为人体肺部X射线成像后的原图像、原图直方图、增强后图像和增强后图像直方图，相比增强前图像，多尺度Retinex算法增强后图像细节有更好的描述，图像信息丰富，边缘清晰，且得到的灰度直方图相比原图有更好的灰度分布；

2.表1为图2图像定量质量评价结果，增强改进后，图像质量在对比度、信息和清晰度上有更好的计算值，更有利于提高医学图像质量。

表1质量评价

图2中，图(a)和(b)为两幅不同位置人体肺部图像，图(c)为图(b)增加噪声、改变照度等处理的结果，图(d)为图(a)和图(b)对数极坐标变换改进的相位相关法冲激函数能量分布情况，图(e)是图(a)和图(c)冲激函数能量分布情况，图(d)、(e)中最大峰值的大小表示两图重叠区域多少，最大峰值的位置表示旋转缩放及平移运动情况；图(e)中由于噪声、照度等影响，出现分散分布的小峰值，但大峰值仍很显著，可以据大峰值参数信息粗估计图像重叠区域位置，缩小了后续检测的特征点提取范围，提高效率，之后即可根据两两图像的峰值大小及制定的排序规则实现多图自动排序。

图3为通过SUFR局部特征提取算法提取待配准区域特征点，(a)图和(b)图分别为特征点的欧氏距离阈值在0.4和0.7的情况下未经精确提纯处理的匹配结果，当阈值较小时，匹配较好，但匹配点较少，不利于后期建立参数模型，当阈值较大时，由于图像本身分辨率低，导致误匹配现象明显；(c)图为通过改进的随机抽样一致算法精确匹配，去除错误点后，欧氏距离阈值0.7时的匹配结果；(d)-(f)图分别为第二幅图加强光照、顺时针旋转10度，以及增加高斯模糊后的匹配效果，可以看出，提纯后的SURF有较好抗复杂环境的能力。

图4中(a)-(f)为6幅局部胯骨排序前图像，曝光度不同，部分图像对比度差；表2、表3为图4中(a)-(f)6幅图像根据制定的排序规则，找出首尾以及中间序列图像的过程。表2的内容为通过改进的相位相关算法两两互匹配后相关系数值，选取每行每列最大两个相关系数值后，放入表3中，表3数值下方括号内的值表示图像间平移量大小，其中相交后仅有一个值的为(b)、(c)两幅图，根据规则，(b)为首图，(c)为尾图，其余为中间图像，根据首图(b)可判断出(e)为(b)的后续图像，接着判断出(f)为(e)的后续图像，依次排序。如图5所示，最终排序结果为(b)-(e)-(f)-(a)-(d)-(c)，解决了多图手工排序困难、出错率高的问题。

表2各图间相关系数

表3行列求交及平移量

图6中(a)图为直接将排序后的多图进行拼接的结果，可以看出接缝痕迹明显，(b)和(c)分别为常用的平均梯度融合拼接、小波变换融合拼接结果，平均梯度方法抗干扰能力差，融合清晰度不够；基于小波分解的Mallat算法低维效果较好，但不能稀疏表示高维特性，且缺乏平移不变性，接缝过度仍不够理想，(d)图为本发明拼接结果，由于拼接前通过改进的多尺度Retinex方法增强图像，改善了图像的对比度、突出了细节内容，利用NSCT变换算法，制定融合策略，通过多尺度分解后的低频、高频子带进一步对拼接图像进行融合处理，使得图像间接缝过渡平滑自然。