CN101901481B - 一种图像拼接方法 - Google Patents

Abstract

本发明提供一种图像拼接方法，具体是一种二维宽景超声图像拼接方法，应用于实时获取或已存储的超声视频图像序列，利用相邻两帧图像的高度相关性，即配准图和活动图，检测活动图像中关键特征点，估计活动图和配准图之间的运动向量，用于完成相邻图像的特征点之间的运动匹配和全局运动参数估计，即得到一个旋转角度和位移量，进而完成两幅图像的配准过程；根据得到的运动参数将当前图像映射到配准图上或者宽景图像的坐标系，进而实现已有的宽景图和配准后的图像之间的融合。如果宽景图像超过当前显示设备窗口大小，将图像整体向后移动若干个像素点，直至宽景图像中新拼接的部分完全显示出来；否则，处理下一幅图像。

Description

一种图像拼接方法

技术领域

本发明提供一种图像拼接方法，尤其涉及一种将多幅二维图像按时间顺序进行配准、融合、拼接以及生成全景图像处理的方法。

背景技术

图像拼接(image mosaic)是计算机视觉热门的研究领域，已经成为多媒体、医学图像处理和计算机图形学领域中的热点问题。图像拼接问题可以定义成对一系列空间重叠的图像进行特征提取、运动参数估计及图像融合和增强的处理过程，最终生成一个无缝的、高分辨率的宽视角图像。

图像拼接技术主要分为三个主要步骤：图像预处理、图像配准、图像融合与边界平滑，图像预处理主要指图像配准前，将图像进行噪声抑制、纹理和对比度的增强以及直方图的归一化等预处理，使配准图和活动图不存在明显的差异。图像预处理主要是为图像配准做准备，让图像质量能够满足图像配准的要求：如果图像质量不理想时不经过图像预处理就进行图像配准，很容易造成两幅图像中的一些子区域之间或特征点之间的误匹配。图像配准主要指基于配准图和活动图中的图像关键特征或灰度信息，寻找最佳的匹配特征点或子区域，搜索每个子区域对或特征点对之间运动向量，最终估计出两幅图像之间的全局的、线性或非线性的运动变换参数。

图像拼接的成功与否关键是图像的配准的效果。然而，通常的待配准的医学图像之间，不同目标区域可能存在多种非线性变换，或者存在大面积的无明显特征区域（如均匀纹理区域或同色区域等），这些情况大大增加了图像匹配的难度和挑战。因此，一个好的图像配准算法应该能够在多种情况下准确找到图像间的对应信息和特征点，将图像对齐。

图像融合指在完成图像匹配以后，对图像进行拼接、缝合，并对缝合的边界进行平滑处理，让缝合边界区域自然过渡。

图像配准主要基于活动图和配准图中的匹配特征或灰度信息，寻找最佳匹配的特征点对或图像子区域对，计算每个特征点对或子区域对间的运动向量，并估计两幅图像之间的全局的线性或非线性的运动变换参数。大多数超声图像拼接技术采用最小二乘法原则估计全局的运动变换参数，如US5782766,CN1839760A，CN101455576A等都采用最小二乘法则来估计全局运动参数，而且由于处理时间上的限制，参数估计和特征点匹配都是通过一次迭代完成。

应用SAD相关的方法对提取的特征点或子区域进行匹配的准确率较低，而且即使在较为理想的匹配结果中，有一部分特征点或子区域会落在在局部的运动组织上时，这样特征点或子区域便成为所谓的孤立点（outlier）,会造成运动估计的误差；这主要是由于局部的组织运动区域和背景有不同的运动参数，而基于灰度值的SAD匹配结果再直接应用最小二乘法估计运动变换参数的方法无法克服孤立点所造成的影响。

大部分医学图像拼接、配准方法都基于特征点或子区域的匹配来完成运动估计，找到匹配的特征点对或子区域对之后，再应用贝叶斯估计来完成对运动参数的计算，这样的配准过程通常无法收敛到最优解而且耗时较长，一个较为直观的办法是采用迭代近邻点算法（Iterative Closest Point Method）来提高配准精度，简称ICP算法，这样使得运动参数估计得以收敛到一个状态空间的最佳点。然而，ICP算法的有效性非常依赖于运动参数的初始值，为了有效解决这一问题，该技术方案考虑应用图像拼接中的传统Ransac算法来搜索一个较为联想的运动参数的初始值。

现有图像拼接技术都是在配准图的相应的特征点或子区域的一个更大的图像领域内，搜索与活动图相匹配的特征点和子区域来实现特征点匹配算法的，该匹配方法的缺点是搜索过程非常耗时，而应用启

发性搜索只有可能找到一个局部最优解，如CN101455576A应用爬山法来解决特征点匹配算法的耗时问题。CN1839760A指出，由于噪音和组织运动等因素，会造成特征点的运动偏移量的误差，但其并没有给出有效的解决办法。运动的目标区域内的特征点通常不参与运动参数的估计，这是由于运动区域的特征点本身通常不是刚性运动，当目标区域运动速度较快时，可以考虑对运动区域进行有效的检测并单独进行配准。

综上所述，现有图像拼接技术无法解决特征点或子区域的匹配精度差问题，也无法有效地处理运动目标对刚性配准参数估计造成的影响，以及特征点和子区域匹配搜索算法的耗时问题。

发明内容

本发明提供一种图像拼接方法，本发明考虑到应用视频处理中的运动检测技术，在每个特征点上，计算伪运动向量（pseudo motionvector）用来补偿已经估计到的运动向量，这样的处理不仅可以缩小在特征点匹配过程中的搜索范围，而且可以更快速到找到最优匹配，从而为直接地解决图像中由于感兴趣目标区域运动而产生的图像配准误差。本发明提出应用运动补偿方案以及运动目标区域的滤除方法，更快速到找到最优的匹配的运动向量，并且可以有效地减轻运动目标区域所带来的图像配准误差。

本发明为解决上述技术问题所采用的技术方案为：

一种图像拼接方法，图像序列为，I₁，...,I₁，...,I_n，其时间间隔为Δ，包括以下步骤：

A．建立二维坐标，将第一帧图像设成配准图f和拼接图X₀，并设置活动图g=I₂和处理步长为Δ=1及i=2；

B．在活动图g上寻找一系列特征点：F={(x_j，y_j)|j=1,…,m}，

其中(x_j，y_j)为第j个特征点在活动图中的位置即坐标值；

C．基于特征点集F，迭代搜索初始的最优的子集使得初始运动参数估计

最佳，其中α₀为运动变换中的旋转角度,

和

分别是在X轴方向和Y轴方向的位移参数；

D．对每个特征点，估计在配准图上所对应的特征点的位置：{(x*_j，y*_j)∈f|j=1，...,m}，并迭代求解一组最佳的运动参数；

E．基于估计到的运动参数T=（α,C^X，C^Y）,将当前活动图Ii拼接到当前拼接图，再将活动图其设置成配准图,若拼接的图像宽度大于显示图像窗口，则将图像中所有像素点向后平移至图像新拼接的部分或区域可以完全显示出来；

F．i=i+Δ，如果第i帧图像存在，则返回步骤B;否则将拼接图像左移若干个像素点，直至图像新拼接的部分显示出来，并输出拼接图像。

所述B步骤包括以下步骤：

B1.对当前活动图像g计算其水平方向梯度值和垂直方向梯度值：g_x和g_y，核函数分别为：

$h_x=\left[\begin{array}{ccc}1& 0& -1\\ 2& 0& -2\\ 1& 0& -1\end{array}\right]$

$h_y=\left[\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right]$

并求解梯度幅值图像g_e=|g_x|+|g_y|；

B2.为了将得到的幅值图像g_e做二值化，在[0,255]范围内线性搜索一个阈值ftg，满足以下度量f_b值最小；

$f_b=\frac{s\·t{({\stackrel{\‾}{z}}_1-{\stackrel{\‾}{z}}_2)}^2}{s+t}$

${\stackrel{\‾}{z}}_1=\underset{g_e(x,y)>\mathrm{ftg}}{\mathrm{\Σ}}{g}_e(x,y)/\left|\right\{(x,y)|g_e(x,y)>\mathrm{ftg}\}|$

${\stackrel{\‾}{z}}_2=\underset{g_e(x,y)\≤\mathrm{ftg}}{\mathrm{\Σ}}{g}_e(x,y)/\left|\right\{(x,y)|g_e(x,y)\≤\mathrm{ftg}\}|$

s=|{(x,y)|g_e(x,y)>ftg}|

t=|{(x,y)|g_e(x,y)≤ftg}|

B3.根据得到的阈值ftg计算二值图像g_b，然后对二值图像g_b进行平滑滤波，滤波窗口大小为w₁×w₁，w₁默认值为5，滤波窗口内每个系数值均为1，得到的滤波图像为一个非二值图像为g_F；

B4.将活动图像均分成m个不重叠的子区域，每个区域是个w₂×w₂，例如针对512x512大小的图像，w₂默认值可以取16；

B5.在每个子区域搜索出一个特征点集合F_i={(x_j，y_j)|j=1，...,m}使其满足该点的在图像g_F(x，y）上的灰度值最大,并且g_F(x，y）>Th，阈值Th默认值为5。

所述C步骤包括以下步骤：

C1.估计初始迭代参数，

如果|β–γ|<ε=0.0001,T₀=T_i-1，E=F，，β和γ分别为前两帧运动旋转的角度，否则，执行以下步骤；

C2.在前特征值排在前59位的特征点作为子集，在这一子集中随机选择4个特征点为特征集合F₁；

C3.对F₁其中的每一像素点(x，y)，将其旋转为β角度以后，在配准图以相应的点为中心的一个区域中搜索和其匹配的、并使匹配值SAD值最小的像素点，匹配的模板窗口大小为w₃×w₃，w₃默认值为5，搜索区域为[x–|Prev_C^X|,x+2|Prev_C^X|]×[y–|Prev_C^Y|,y+|Prev_C^Y|]，这里假设探头的运动方向为从左向右，Prev_C^X和Prev_C^Y是前一帧的最终运动参数估计值；

C4.对F₁中的每个像素点，在配准图f中都找到相应的匹配点，然后应用最小二乘法估计参数 $T_R=({\mathrm{\&alpha;}}_R,C_R^x,C_R^Y):$

$\left(\begin{array}{c}x*\\ y*\end{array}\right)=\left(\begin{array}{cc}\cos \left({\mathrm{\&alpha;}}_R\right)& -\sin \left({\mathrm{\&alpha;}}_R\right)\\ \sin \left({\mathrm{\&alpha;}}_R\right)& \cos \left({\mathrm{\&alpha;}}_R\right)\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{C}_R^X\\ C_R^Y\end{array}\right)$

C5.在配准图f中找到相应的特征点集F*=T_R(F)，以F中的每个像素点(x，y）为中心的模板和以F*中对应的映射像素点(x*,y*）为中心的模板之间的Normalized Correlation的值NC，模板大小为w3×w3:

$\mathrm{NC}=\frac{{\mathrm{\&Sigma;}}_{k=1}^K(a_k-\stackrel{\&OverBar;}{a})(b_k-\stackrel{\&OverBar;}{b})}{\sqrt{{\mathrm{\&Sigma;}}_{k=1}^K{(a_k-\stackrel{\&OverBar;}{a})}^2}\sqrt{{\mathrm{\&Sigma;}}_{k=1}^K{(b_k-\stackrel{\&OverBar;}{b})}^2}}$

a和b为像素点在活动图g上和配准图f上的灰度值;

C6.根据NC的值，再将F分成两组A={(x，y)|NC((x，y),(x*,y*))>ftg2}和B={(x，y)|NC((x，y),(x*,y*))≤ftg₂}，ftg₂为一阈值，默认为0.7；如果|A|/|F|>0.5，T₀=T_R，E=A，进入步骤C7；否则重复C2－C5；

C7.输出E和T₀。

所述D步骤包括以下步骤：

D1.将特征点集E映射到配准图f上，得到新的特征点集合E*=T(E);

D2.对E中的每一像素点(x，y)在配准图中相应的匹配点(x*,y*)=T(x，y)∈E*，执行基于光流的运动补偿,

x_M=x*-g_Δt(x,y)/g_x(x,y);

y_M=y*-g_Δt(x,y)/g_y(x,y);

g_Δt(x,y)=f(x*,y*)-g(x,y);

；

D3.在以(x_M,y_M)中心的子区域,[x_M-w₄,x_M+2×w₄]×[y_M-w₅,y_M+w₅]，搜索和(x，y)更为匹配的像素点(x₂,y₂)使匹配的NormalizedCorrelation的值最小，模板匹配的窗口大小为w₃×w₃，w₄和w₅的默认值为3；

D4.找出E中子集E₁={(x，y)|NC（(x，y)，(x*,y*））≤1-ftg₃}，ftg₃默认值可以为0.3；在配准图上有相应的匹配点集E2；

D5.如果集合|E₁|>10，基于集合E₁和E₂应用M-estimator估计运动参数T=（α,C^X，C^Y），否则直接跳到步骤D6，这里Err_k是第k个特征点的匹配误差；

D6.重复子步骤D1至D5共5次，如果最终检测的运动向量大于设定阈值，跳到步骤D1，否则输出T。

所述E步骤包括以下步骤：

E1.根据运动变换T到活动图g映射成f*=T_i*(g)，图像f*的灰度值可以通过双线性插值的到;

E2.令f_M=f*,将当前拼接图X_i-1与f_M进行融合，如果当前像素点的既在前M-1帧出现过，又在f_M出现过，进行基于最近邻的图像融合，即：

$f(x,y)=\frac{\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{w}_i{f}_i(x,y)}{\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{w}_i}$

$w_i=\exp (-\underset{P\&Element;R\left(N\right)}{\mathrm{\&Sigma;}}|f_i(x,y)-f_i(x_P,y_P)\left|\right)$

R(N)是像素点(x,y)的3×3邻域；

E3.对新生成的拼接图分别做一次Laplacian边缘增强和各项同性滤波。

所述F步骤包括以下步骤：

F1.如果最终检测的运动向量小于或等于设定阈值,Δ=Δ+1；

F2.如果最终检测的运动向量大于设定阈值：如果Δ>1，则Δ=Δ-1并且i=i-1，返回步骤B,否则终止图像拼接，提示调整探头运动速度；

F3.令i=i+Δ，如果第i帧图像存在并且拼接后图像宽度小于感兴趣区域的窗口宽度，则返回步骤B;否则，左移图像若干个像素点，直至图像新拼接的部分全部落在感兴趣区域的窗口，并输出图像。

本发明提供了一种更快速有效的图像配准方法，提高配准的精确度。

附图说明

图1是本发明实施例二维宽景超声图像成像系统硬件结构图。

具体实施方式

下面根据附图和实施例对本发明作进一步详细说明:

图1是该发明所应用的超声成像系统的结构组件图：用户通过输入设备配置超声扫描所需的参数及宽景成像的配置参数，用户输入设备主要和计算机中央处理器连接，可以对图像存储区域12中的数据进行图像显示、处理,存储及打印等操作。超声图像和宽景图像都可以通过显示器13显示给用户；计算机中央处理器和超声系统的主处理器14相连接，并可通过其对超声成像过程中的参数配置、及模式选择进行操作。探头16在受检器官的皮肤表面进行运动扫描，用以获取感兴趣区域15的超声图像，而目标区域15可能需要一幅大视野的图像来显示整个器官或组织，比如对一个较大器官的尺寸的测量。超声成像系统可将图像序列22拼接成一个宽景图像24，然后将目标区域15的完整组织结构信息显示给用户。

发射控制单元14产生的发射波形经过探头16向目标区域15发射超声波，接收到的超声回波，经过波束合成器17的延时控制、通道合成等处理后，接收控制器18则将射频信号转换成视频信号进一步处理，处理后的信号经过扫描转换单元20处理后转换成超声视频序列22；主处理器14可以通过传输控制电路19设置、调整和监测超声系统的工作参数，并用以产生电流驱动探头工作。

以上描述简单介绍了宽景超声成像的设备，这和大多数医用超声系统体系结构没有差别，而基于超声视频的宽景图像的拼接方法基本不依赖于超声设备。以下描述从超声图像序列22生成宽景图像24的过程：

1.从图像存储区域中提取图像帧序列,I₁，...,I_i，...,I_n，将配准图和活动图设成：f=I₁，g=I₂；再取宽景图像X₁=f,且有i=2，步长为Δ；

2.在第i幅图像上，也就是当前活动图，寻找一系列特征点：F={(x_j,y_j)|j=1,...,m}；

3.基于特征点集F,迭代找出初始的最优的子集

使得初始运动参数估计

最佳；

4.基于特征点子集E和初始运动变换T₀，迭代搜索E中每个特征点的运动到配准图上的相应像素点，即集合E*={(x*_j，y*_j)∈I_i|j=1，...,m_E}，基于E和E*应用M-estimator迭代估计最佳的运动参数T_i=（α,C^X,C^Y）；

5.基于T_i=（α,C^X,C^Y）,把活动图g融合到宽景图X_i-1；

6.i=i+Δ如果第i帧图像存在，则返回步骤2；否则结束退出。

上述处理程序包含三大部分，图像子区域的划分和关键特征点的检测，活动图像和配准图像之间的配准，活动图与宽景图像的融合。图像配准的关键特征点通常决定了图像配准的效果，好的特征点会克服图像中的噪音等图像质量差的问题，本发明基于图像边缘提取和分割技术，采用了一个鲁棒性较强的特征提取的方法检测关键特征点集合。一种实现鲁棒性较强的特征点检测包含以下几个子步骤：

(1)对当前活动图像20计算其水平方向梯度值和垂直方向梯度值：g_x和g_y，核函数分别为：

$h_x=\left[\begin{array}{ccc}1& 0& -1\\ 2& 0& -2\\ 1& 0& -1\end{array}\right]$

$h_y=\left[\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right]$

并求解梯度幅值图像g_e=|g_x|+|g_y|

(2)为了将得到的幅值图像g_e做二值化，在[0,255]范围内线性搜索一个阈值ftg,满足以下度量f_b值最小；

$f_b=\frac{s\&CenterDot;t{({\stackrel{\&OverBar;}{z}}_1-{\stackrel{\&OverBar;}{z}}_2)}^2}{s+t}$

${\stackrel{\&OverBar;}{z}}_1=\underset{g_e(x,y)>\mathrm{ftg}}{\mathrm{\&Sigma;}}{g}_e(x,y)/\left|\right\{(x,y)|g_e(x,y)>\mathrm{ftg}\}|$

${\stackrel{\&OverBar;}{z}}_2=\underset{g_e(x,y)\&le;\mathrm{ftg}}{\mathrm{\&Sigma;}}{g}_e(x,y)/\left|\right\{(x,y)|g_e(x,y)\&le;\mathrm{ftg}\}|$

s=|{(x,y)|g_e(x,y)>ftg}|

t=|{(x,y)|g_e(x,y)≤ftg}|

(3)根据得到的阈值ftg计算二值图像g_b，然后对二值图像g_b进行平滑滤波，滤波窗口大小为w₁×w₁（w₁默认值为5），滤波窗口内每个系数值均为1，得到的滤波图像为一个非二值图像为g_F

(4)将活动图像均分成m个不重叠的子区域，每个区域是个w₂×w₂，例如针对512×512大小的图像，w₂默认值可以取16；

(5)在每个子区域搜索出一个特征点集合F_i={(x_j，y_j)|j=1，...,m}使其满足该点在图像g_F(x,y）上的灰度值最大，并且g_F(x,y）>Th，阈值Th默认值为5。

图像序列U放置在一段存储区域（内存或硬盘），按图像的生成的时间先后顺序排放，并按同样顺序读取。本发明的特征点提取中二值化算法采用了一种动态阈值的方法，即将活动图中的所有像素的灰度值作为样本进行分类，当寻找一个最优的阈值使得类内距离之和与类间距之比达到最小。基于二值化图像即Edge图像计算每个像素点周围像素点灰度值之和，作为效果特征来检测每个图像子区域的特征点，这样的特征点提取办法具有更强的鲁棒性。

通常来说，图像特征点的数目对图像配准精度影响很大，特征点越多，最后配准的误差越小，但会使配准过程更耗时。即使有足够的特征点，往往在配准图中搜索对应的特征点精度也不是很理想，并且需要在一个很大区域搜索与之相匹配的像素点，换句话说，大部分的特征点匹配过程非常依赖于初始的运动变换参数，为了解决这一技术难题我们应用计算机视觉中的Ransac算法来同时估计有效的特征点集合和初始运动变换参数，找到一个最优的初始运动变换参数和一个有效的特征子集。应用Ransac算法来估计初始的最优特征点子集和初始的运动参数,其包含以下处理过程(1)-(7):

(1)估计初始迭代参数，

如果|β–γ|<ε=0.0001,T₀=T_i-1，E=F,(β和γ分别为前两帧运动旋转的角度)否则，执行以下步骤：

(2)在前特征值排在前59位的特征点作为子集，在这一子集中随机选择4个特征点为特征集合F₁；

(3)对F₁其中的每一像素点(x,y)，将其旋转为β角度以后，在配准图以相应的点为中心的一个区域中搜索和其匹配的、并使匹配值SAD值最小的像素点，匹配的模板窗口大小为w₃×w₃（w₃默认值为5），搜索区域为[x–|Prev_C^X|,x+2|Prev_C^X|]×[y–|Prev_C^Y，y+|Prev_C^Y|]，这里假设探头的运动方向为从左向右，Prev_C^X和Prev_C^Y是前一帧的最终运动参数估计值；

(4)F₁中的每个像素点中找到在配准图中找到相应的匹配点，应用最小二乘法估计参数 $T_R=({\mathrm{\&alpha;}}_R,C_R^x,C_R^Y)$

$\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}\cos \left(\mathrm{\&alpha;}\right)& -\sin \left(\mathrm{\&alpha;}\right)\\ \sin \left(\mathrm{\&alpha;}\right)& \cos \left(\mathrm{\&alpha;}\right)\end{array}\right)+\left(\begin{array}{c}{C}_0^X\\ C_0^Y\end{array}\right)---\left(5\right)$

(5)在配准图R中找到相应的特征点集F*=T_R(F)以F中的每个像素点(x,y）为中心的模板和以F*中对应的映射像素点(x*,y*）为中心的模板之间的Normalized Correlation的值NC，模板大小为w3×w3:

$\mathrm{NC}=\frac{{\mathrm{\&Sigma;}}_{k=1}^K(a_k-\stackrel{\&OverBar;}{a})(b_k-\stackrel{\&OverBar;}{b})}{\sqrt{{\mathrm{\&Sigma;}}_{k=1}^K{(a_k-\stackrel{\&OverBar;}{a})}^2}\sqrt{{\mathrm{\&Sigma;}}_{k=1}^K{(b_k-\stackrel{\&OverBar;}{b})}^2}}---\left(6\right)$

a和b为像素点在活动图g上和配准图f上的灰度值;

(6)接着按NC的值将F分成两组A={(x，y)|NC（(x，y)，(x*,y*））>ftg2}和B={(x，y)|NC（(x，y)，(x*,y*））≤ftg₂}，ftg₂为一阈值，默认为0.7；如果|A|/|F|>0.7，T₀=T_R，E=A，进入步骤（7）；否则重复(2)－(5)；

(7)输出E和T₀；

RANSAC（Random Sample Consensus）是随机抽样一致性算法的缩写。它可以应用于任何一个基于数据集的模型估计问题。给定一个数据集，如果数据集中有一定数目的奇异点（Outliers），在图像配准过程中，由配准算法检测到的图像特征点往往包含大量的奇异点，这是由于图像获取过程中所导致的图像质量下降和图像噪音所引起的；而传统的模型参数的估计方法无法去除这些奇异的特征点对参数估计的影响，解决这一问题可以用随机抽样一致性算法。Ransac算法是一种随机的优化算法，因此，即使每次运算求出的结果可能尽不相同，但每一次随机迭代都有可能给出一个更合理的结果，因此提高迭代次数会改善模型估计的效果。

Ransac实际上是一种随机采样的方法，其目标是找到一个最小采样数目的最优数据子集来估计模型参数，以上步骤A-G是Ransac方法在特征点匹配算法中的一种实现：E中点的数目越多，参数估计的越准确。

在得到比较理想的初始特征点子集E和初始运动参数估计T0之后，我们应用迭代近邻点算法来进一步求解更为精确的图像配准所需的运动变换参数T：设T=T₀，然后执行ICP即迭代近邻点算法的以下处理(1)-(7)：

(1)将E映射到配准图f上得到新的映射点集E*=T（E);

对F中的每一像素点(x，y)在配准图中相应的匹配点(x*,y*)=T(x，y)∈E*执行基于光流的运动补偿,

x_M=x*-g_Δt(x,y)/g_x(x,y);

y_M=y*-g_Δt(x,y)/g_y(x,y);

g_Δt(x,y)=f(x*,y*)-g(x,y);

；

(2)在以(x_M，y_M)中心的子区域[x_M-w₄,x_M+2×w₄]×[y_M-w₅,y_M+w₅]，搜索和(x，y)更为匹配的像素点(x₂,y₂)使匹配的NormalizedCorrelation的值最小，模板匹配的窗口大小为w₃×w₃（w₄和w₅的默认值为3）；

(3)找出E中子集E₁={(x，y)|NC（(x，y)，(x*,y*））≤1-ftg₃}，ftg₃默认值可以为0.3；在配准图上有相应的匹配点集E₂；

(4)如果集合|E₁|>10，基于集合E₁和E₂应用M-estimator估计运动参数T=（α,C^X，C^Y），否则直接跳到步骤（5），这里Err_k是第k个特征点的匹配误差；

(5)重复子步骤(1)至(4)共5次，如果最终检测的运动向量非常大，跳到步骤（1），否则输出T。

超声图像拼接技术的特征点匹配算法是通过在配准图上的相应的像素点的一个更大的邻域里，搜索与活动图上的特征点相匹配的特征点，这样的搜索过程非常耗时，而应用启发性搜索一般情况下只能找到一个局部最优解，应用迭代近邻点算法虽然可以解决局部最优解的问题，但是搜索本身的耗时问题没有完全解决，且有可能搜索时间变得更长；但基于每个像素点上的灰度信息计算一个伪运动向量可以用来补偿已经估计到的运动向量，这样的处理可以缩小在特征点匹配过程中的搜索邻域的范围，从而解决特征点匹配的耗时问题；如果感兴趣区域存在运动组织，需要滤除运动目标区域内的特征点，可以考虑应用相邻两帧或多帧图像相减来检测的运动物体区域，，图像相减的结果也可以应用形态学方法进处理，例如Top-hat算法和Bottom-hat算法，把“运动物体区域”中的关键特征点滤除。这样的处理可以消除超声图像噪音和组织运动等因素对图像配准所造成的误差。

图像配准后，可以根据得到的活动图和配准图之间的运动变换参数，将活动图g=I_i拼接到当前的宽景图像X_i-1，完成图像融合部分。实际上，该步骤是将活动图根据运动参数T将活动图融合到当前配准图f上的一个过程：可按时间顺序将活动图g拼接到现有的宽景图X_i-1上，生成新的宽景图像X_i：假设当前活动图为g=I_i和配准图f=I_i-1，那么T_i-1是配准过程得到的配准图f和活动图g之间的刚性运动变换；特别地，活动图经过配准后的生成图像是经过若干次迭代完成的，即图像T_i*(I_i)，且T_i*=T_i·T_i-1...T₀，再将已有的宽景图像X_i-1和配准后的图像T_i*(I_i)进行融合成新的宽景图X_i。全景图中的任何一个像素点可能会在序列中的若干个帧图像中都出现，因此在计算全景图中该点的值时，可以充分利用此点在视频帧序列中所有出现的值。考虑到运动组织和噪音的存在以及成像参数的差异性等原因的影响，可以对这些值进行简单的处理，如一下步骤中的子步骤(2)。

(1)按运动变换T到活动图g映射成f*=T_i*(g)，图像f*的灰度值可以通过双线性插值的到;

(2)令f_M=f*,将当前拼接图X_i-1与f_M进行融合，如果当前像素点的既在前M-1帧出现过，又在f_M出现过，进行基于最近邻的图像融合，即：

$f(x,y)=\frac{\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{w}_i{f}_i(x,y)}{\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{w}_i}---\left(7\right)$

$w_i=\exp (-\underset{P\&Element;R\left(N\right)}{\mathrm{\&Sigma;}}|f_i(x,y)-f_i(x_P,y_P)\left|\right)$

R(N)是像素点(x，y)的3×3邻域；

(3)对新生成的拼接图分别做一次Laplacian边缘增强和各项同性滤波。

当探头运动过快或过慢时，都会影响视频图像相邻帧的重叠的信息量，直接影响图像配准的速度和效果，适当地调节配准步长Δ会提高图像配准和拼接的效率：如果新生成的宽景图Xi的宽度大于图像显示窗口，则将该图像中所有像素点向后平移若干个像素点，直至图像新拼接的部分可以完全显示出来；i=i+Δ如果第i帧图像存在，则返回步骤2，宽景图像拼接步长调整的流程如下:

(1)如果最终检测的运动向量非常小，Δ=Δ+1；

(2)如果最终检测的运动向量非常大：如果Δ>1，则Δ=Δ-1并且i=i-1，返回步骤II,否则终止图像拼接，提示调整探头运动速度；

(3)令i=i+Δ，如果第i帧图像存在并且拼接后图像宽度小于感兴趣区域的窗口宽度，则返回步骤II;否则，左移图像若干个像素点，直至图像新拼接的部分全部落在感兴趣区域的窗口，并输出图像。

本领域技术人员不脱离本发明的实质和精神，可以有多种变形方案实现本发明，以上所述仅为本发明较佳可行的实施例而已，并非因此局限本发明的权利范围，凡运用本发明说明书及附图内容所作的等效结构变化，均包含于本发明的权利范围之内。

Claims (5)

1.一种图像拼接方法，其图像序列为：I₁，...，I_i，...，I_n，其时间间隔为Δ，包括以下步骤：

A.建立二维坐标，将第一帧图像设成配准图f和拼接图X₀，并设置活动图g＝I₂和处理步长为Δ＝1及i＝2；

B.在活动图g上寻找一系列特征点F＝{(x_j，y_j)|j＝1，...，m}，其中(x_j，y_j)为第j个特征点在活动图中的位置坐标值；

C.基于特征点集F，迭代找出初始的最优的子集

基于E估计的运动参数

为在一个有效的特征点子集

上的最佳匹配效果，其中α₀为运动变换中的旋转角度，

和

分别是在轴X方向和Y轴方向的位移参数；

D.对每个特征点，估计在配准图上所对应的特征点的位置：{(x^* _j，y^* _j)∈f|j＝1，...，m}，其中(x^* _j，y^* _j)为第j个特征点在配准图中的位置坐标值，并求解一组最佳的运动参数T；

E.基于求解的运动参数，将当前活动图I_i拼接到当前拼接图，再将该活动图设置成配准图，若拼接的图像宽度大于显示图像窗口，则将图像中所有像素点向后平移至图像新拼接的部分或区域可以完全显示出来；

F.i＝i+Δ如果第i帧图像存在，则返回步骤B；否则将拼接图像左移若干个像素点，直至图像新拼接的部分显示出来，并输出拼接图像，

其特征在于：所述D步骤中的求解方式为迭代求解，所述D步骤包括以下步骤：

D1.令T＝T₀，将E映射到配准图f上得到新的映射点集E^*＝T(E)；

D2.对E中的每一像素点(x，y)在配准图中相应的匹配点(x^*，y^*)＝T(x，y)∈E^*执行基于光流的运动补偿：

x_M＝x^*-g_Δt(x，y)/g_x(x，y)；

y_M＝y^*-g_Δt(x，y)/g_y(x，y)；

g_Δt(x，y)＝f(x^*，y^*)-g(x，y)；

D3.在以(x_M，y_M)中心的子区域，[x_M-w₄，x_M+2×w₄]×[y_M-w₅，y_M+w₅]，搜索和(x，y)更为匹配的像素点(x₂，y₂)使匹配的Normalized Correlation的值最小，模板匹配的窗口大小为w₃×w₃，w₄和w₅的默认值为3；

D4.找出E中子集E₁＝{(x，y)|NC((x，y)，(x^*，y^*))≤1-ftg₃}，ftg₃默认值可以为0.3；在配准图上有相应的匹配点集E₂；

D5.如果集合|E₁|＞10，基于集合E₁和E₂应用M-estimator估计运动参数T＝(α，C^X，C^Y)，否则直接跳到步骤D6，这里Err_k是第k个特征点的匹配误差；

D6.重复子步骤D1至D5共5次，如果最终检测的运动向量大于设定阈值，跳到步骤D1，否则输出T。

2.根据权利要求1所述的一种图像拼接方法，其特征在于所述B步骤包括以下步骤：

B1.对当前活动图像计算其水平方向梯度值和垂直方向梯度值：g_x和g_y，核函数分别为：

$h_x=\left[\begin{array}{ccc}1& 0& -1\\ 2& 0& -2\\ 1& 0& -1\end{array}\right]$

$h_y=\left[\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right]$

并求解梯度幅值图像g_e＝|g_x|+|g_y|；

B2.将得到的幅值图像g_e做二值化，在[0，255]范围内线性搜索一个阈值ftg，满足以下度量f_b值最小；

$f_b=\frac{s\&CenterDot;t{({\stackrel{\&OverBar;}{z}}_1+{\stackrel{\&OverBar;}{z}}_2)}^2}{s+t}$

${\stackrel{\&OverBar;}{z}}_1=\underset{g_e(x,y)>\mathrm{ftg}}{\mathrm{\&Sigma;}}{g}_e(x,y)/\left|\right\{(x,y)|g_e(x,y)>\mathrm{ftg}\}|$

${\stackrel{\&OverBar;}{z}}_2=\underset{g_e(x,y)\&le;\mathrm{ftg}}{\mathrm{\&Sigma;}}{g}_e(x,y)/\left|\right\{(x,y)|g_e(x,y)\&le;\mathrm{ftg}\}|$

s＝|{(x，y)|g_e(x，y)＞ftg}|

t＝|{(x，y)|g_e(x，y)≤ftg}|

B3.根据得到的阈值ftg计算二值图像g_b，然后对二值图像g_b进行平滑滤波，滤波窗口大小为w₁×w₁，w₁默认值为5，滤波窗口内每个系数值均为1，得到的滤波图像为一个非二值图像g_F；

B4.将活动图像均分成m个不重叠的子区域，每个区域是个w₂×w₂，针对512×512大小的图像，w₂默认值取16；

B5.在每个子区域搜索出一个特征点集合F_i＝{(x_j，y_j)|j＝1，...，m}使其满足该点的在图像g_F(x，y)上的灰度值最大，并且g_F(x，y)＞Th，阈值Th默认值为5。

3.根据权利要求1所述的一种图像拼接方法，其特征在于所述C步骤包括以下步骤：

C1.估计初始迭代参数，

如果|β-γ|＜ε＝0.0001，T₀＝T_i-1，E＝F，β和γ分别为前两帧运动旋转的角度，否则，执行以下步骤；

C2.将前特征值排在前59位的特征点作为子集，在这一子集中随机选择4个特征点为特征集合F₁；

C3.对F₁其中的每一像素点(x，y)，将其旋转为β角度以后，在配准图以相应的点为中心的一个区域中搜索和其匹配的、并使匹配值SAD值最小的像素点，匹配的模板窗口大小为w₃×w₃，w₃默认值为5，搜索区域为[x-|Prev_C^X|，x+2|Prev_C^X|]×[y-|Prev_C^Y|，y+|Prev_C^Y|]，这里探头的运动方向为从左向右，Prev_C^X和Prev_C^Y是前一帧的最终运动参数估计值；

C4.对F₁中的每个像素点，在配准图f中都找到相应的匹配点，然后应用最小二乘法估计参数 $T_R=({\mathrm{\&alpha;}}_R,C_R^X,C_R^Y):$

$\left(\begin{array}{c}{x}^{*}\\ y^{*}\end{array}\right)=\left(\begin{array}{cc}\cos \left({\mathrm{\&alpha;}}_R\right)& -\sin \left({\mathrm{\&alpha;}}_R\right)\\ \sin \left({\mathrm{\&alpha;}}_R\right)& \cos \left({\mathrm{\&alpha;}}_R\right)\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{C}_R^X\\ C_R^Y\end{array}\right)$

C5.在配准图f中找到相应的特征点集F^*＝T_R(F)，以F中的每个像素点(x，y)为中心的模板和以F^*中对应的映射像素点(x^*，y^*)为中心的模板之间的NormalizedCorrelation的值NC，模板大小为w3×w3：

$\mathrm{NC}=\frac{{\mathrm{\&Sigma;}}_{k=1}^K(a_k-\stackrel{\&OverBar;}{a})(b_k-\stackrel{\&OverBar;}{b})}{\sqrt{{\mathrm{\&Sigma;}}_{k=1}^K{(a_k-\stackrel{\&OverBar;}{a})}^2}\sqrt{{\mathrm{\&Sigma;}}_{k=1}^K{(b_k-\stackrel{\&OverBar;}{b})}^2}}$

a和b为像素点在活动图g上和配准图f上的灰度值；

C6.根据得到NC的值，将特征点集合F分成两组A＝{(x，y)|NC((x，y)，(x^*，y^*))＞ftg2}和B＝{(x，y)|NC((x，y)，(x^*，y^*))≤ftg₂}，ftg₂为一阈值，默认为0.7；如果|A|/|F|＞0.5，T₀＝T_R，E＝A，进入步骤C7；否则重复C2-C5；

C7.输出E和T₀。

4.根据权利要求1所述的一种图像拼接方法，其特征在于所述E步骤包括以下步骤：

E1.按运动变换T_i-1，把活动图g映射成f^*＝T_i ^*(g)，图像f^*的灰度值通过双线性插值得到；

E2.令f_M＝f^*，将当前拼接图X_i-1与f_M进行融合，如当前像素点在前M-1帧出现过，又在f_M出现过，进行基于最近邻的图像融合，即：

$f(x,y)=\frac{\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{w}_i{f}_i(x,y)}{\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{w}_i}$

$w_i=\exp (-\underset{P\&Element;R\left(N\right)}{\mathrm{\&Sigma;}}|f_i(x,y)-f_i(x_P,y_P)\left|\right)$

R(N)是像素点(x，y)的3×3邻域；

E3.对新生成的拼接图分别做一次Laplacian边缘增强和各项同性滤波。

5.根据权利要求1所述的一种图像拼接方法，其特征在于所述F步骤包括以下步骤：

F1.如果最终检测的运动向量小于或等于设定阈值，Δ＝Δ+1；

F2.如果最终检测的运动向量大于设定阈值：如果Δ＞1，则Δ＝Δ-1并且i＝i-1，返回步骤B，否则终止图像拼接，提示调整探头运动速度；

F3.令i＝i+Δ，如果第i帧图像存在并且拼接后图像宽度小于感兴趣区域的窗口宽度，则返回步骤B；否则，左移图像若干个像素点，直至图像新拼接的部分全部落在感兴趣区域的窗口，并输出图像。

Images