CN101702235B - 基于三角剖分的图像配准方法 - Google Patents

Abstract

本发明涉及一种基于三角剖分的图像配准方法，视频图像的获取，并进行A/D转换；输出数字视频图像信号发送到DSP装置进行基于多分辨率方法的去噪处理后进行角点检测摄取特征点；将特征点进行三角形匹配；匹配后处理得出配准图像输出，此方法具有适用范围广、配准精度高等优点，对于视频信号后续的图像融合处理、基于特征级及决策级等人工智能的处理或者人机交互判断决策等工作具有重要意义和实用价值。

Description

基于三角剖分的图像配准方法

技术领域

本发明涉及一种图像处理技术，特别涉及一种基于三角剖分的图像配准方法。

背景技术

由于多源图像是由不同成像传感器或者是由同种成像传感器的不同成像方式获得的多幅图像，图像间可能出现相对平移、旋转和比例缩放等，不能直接进行融合，而必须先进行图像配准，以建立图像间像素-像素的对应关系。图像配准作为多源图像融合的预处理过程，其重要性是显而易见的，它是多源图像融合得以进行的前提条件，图像配准的精度对融合的精度也有显著影响。

常用的图像配准算法本身分有三类：基于像素灰度信息的图像配准算法，如基于互信息的配准算法；基于频率域的图像配准算法，如基于傅立叶变换的配准算法；基于特征的图像配准算法，如基于轮廓的配准算法。从传感器种类分可以分为两类：同类传感器图像配准方法和异类传感器图像配准方法。

同类传感器多源图像是由成像机理相同的传感器获得的，图像的灰度分布差异不大。根据文献报道，对于未配准图像之间的坐标关系为刚体变换(也就是未配准图像之间仅存在旋转和平移)的同类传感器多源图像配准问题，基于互信息的配准算法精度最高，但是基于互信息的图像配准算法中涉及到大量的浮点运算，时间复杂度大，而且在配准参数优化过程中很有可能收敛到局部极值。

异类传感器多源图像是由成像机理不同的传感器获得的，图像间的灰度分布差异很大，一般的基于像素灰度信息的图像配准算法难以适用。基于特征的图像配准算法从输入图像中提取一些显著特征如轮廓，然后对这些显著特征进行匹配，最后估计配准参数，算法鲁棒性高，适于异类传感器多源图像的配准问题。目前常用的算法有基于轮廓的图像配准算法，该算法只是利用输入图像中的闭合轮廓进行配准，首先提取输入图像中的闭合轮廓，然后利用区域不变矩和链码表示的方法对这些闭合轮廓进行匹配，最后根据匹配后的闭合轮廓的质心求取配准参数。

不论是同类传感器图像配准方法还是异类传感器图像配准方法，都存在一些问题，如对于同类传感器的配准算法，目前问题是涉及到大量的浮点运算，时间复杂度大，而且在配准参数优化过程中很有可能收敛到局部极值。对于异类传感器图像配准，在实际配准情况中，有时很难从输入图像提取足够的闭合轮廓。实验表明，匹配的闭合轮廓对的数量会影响配准参数估计的精度。如从输入图像中不能检测到匹配的闭合轮廓对，则这种配准算法就不能适用，在配准过程中只利用输入图像中的闭合轮廓提供的信息，而丢弃了输入图像中的开轮廓所提供的信息，在实际应用中有一定的局限性。

发明内容

本发明是针对现在多源图像配准方法存在实际应用的局限性问题，提出了一种基于三角剖分的图像配准方法，使用了ARM与DSP组合的嵌入式***，该***分别有调用配准策略的方式：启动阶段、用户手动方式和定时更新方式。不论是同类传感器配准还是异类传感器配置，统一利用Delaunay三角剖分作为对比特征进行配准。能够有效的带动图像配准方法在实际中的运用，并有效改善图像配准的性能。

本发明的技术方案为：一种基于三角剖分的图像配准方法，包括如下步骤：

1)视频图像的获取：将从CCD摄像仪和红外热像仪得到的模拟视频信号转换成为LDVS格式的24位数字像素值，所有的模拟信号所表示的像素经过A/D转换后，得到一帧一帧的数字视频图像信号，将其发送到DSP装置进行基于多分辨率方法的去噪处理；

2)基于多分辨率的分解：采用基于小波框架变换的方法，原始图像与高通滤波器h_n、低通滤波器g_n分别进行滤波，滤波后的结果分别为小波框架在n层的高通系数和n层的低通系数，如果需要继续分解，则需要继续进行插零处理，直到完成既定的分解目标，对于小波框架的高通部分，进行了滤波处理，而低通部分用于进行基于Delaunay的配准策略；

3)角点检测：利用小波框架分解后的低通系数，通过扫描整个低通系数图像的每个象素点完成的，以每个象素点为圆心，以某个事先设定好的r为半径作一个图形学的圆，并在图像的具***置上实现定位，取落在图形学圆上的象素点与圆心处的象素点进行比较，如果所有的图形学圆上的象素点大于或小于圆心出象素点的灰度值到达一个级别t，则认为该点为特征点；

4)三角形匹配：构造一个超级三角形，包含所有特征点，放入三角形链表；将点集中的特征点依次***，在三角形链表中找出其外接圆包含***点的三角形，删除影响三角形的公共边，将***点同影响三角形的全部顶点连接起来，从而完成一个点在Delaunay三角形链表中的***；根据优化准则对局部新形成的三角形进行优化，将形成的三角形放入Delaunay三角形链表；循环执行特征点***，直到所有散点***完毕；

5)配准图像：Delaunay三角形链表中的匹配进行两个刚性特征不变量进行计算，假设l₁，l₂，l₃为三角形的三条边，m₁，m₂，m₃为三角形的三个角：

第一个刚性变换不变量： $V_t=[\frac{l_1}{l_3},\frac{l_2}{l_3},\cos \left(m_3\right)],$

第二个刚性变换不变量：V_m＝[∠m₁，∠m₂，∠m₃]，

然后对他们之间的比值进行比较，当比较后的结果的总和小于一个值时，则可以选择相应的旋转角度参数和平移参数作为最终的配准结果。

所述步骤2)中基于小波框架变换的方法：小波框架的滤波算法程序预定义部分中的小波框架系数的设定所使用的原始小波为DB3小波，其高通滤波系数为：

[0.037828；-0.023849；-0.110624；0.377403；0.852699；0.377404；-0.110624；-0.023849；0.037828]

低通滤波系数为：

[-0.064539；0.040689；0.418092；-0.788486；0.418092；0.040689；-0.064539；0；0]；

对小波系数进行框架化后得到的系数为

其高通滤波系数为：

h₁＝[0.037828；0；-0.023849；0；-0.110624；0；0.377403；0；0.852699；0；0.377404；0；-0.110624；0；-0.023849；0；0.037828]/

低通滤波系数为：

g1＝[-0.064539；0；0.040689；0；0.418092；0；-0.788486；0；0.418092；0；0.040689；0；-0.064539；0；0；0；0]/

以上仅为一次分解时的小波框架系数，如果需要更高层的分解，则需要继续进行插零处理。

本发明的有益效果在于：本发明基于三角剖分的图像配准方法，具有适用范围广、配准精度高等优点，对于视频信号后续的图像融合处理、基于特征级及决策级等人工智能的处理或者人机交互判断决策等工作具有重要意义和实用价值。

附图说明

图1为本发明基于三角剖分的图像配准方法的***硬件结构框图；

图2为本发明基于三角剖分的图像配准方法的软件结构图；

图3为本发明实施例中所用到的CCD图像中实际视频中的一帧；

图4为本发明实施例中所用到的红外图像中实际视频中的一帧；

图5为本发明实施例中最终获得的配准后视频图像的一帧。

具体实施方式

如图1所示基于三角剖分的图像配准方法的***硬件结构框图，采用了两个成像传感器，一个是AVT工业CCD摄像头F-032B/C，另外一个是Infrared Solutions公司可提供一款基于非冷却辐射热计技术的热成像仪：IR-160型成像仪，其可输出160×120像素NTSC或PAL视频输出信号；传感器种类选择开关和DSP图像融合芯片均采用了TI公司的TMS320C6000；A/D多路采集芯片采用24位A/D转换器CS5381；ARM芯片采用了ARM公司生产的ARM940T。

如图2所示基于三角剖分的图像配准方法的软件结构图，软件架构中安排的配准过程分别通过三种方式实现调用。1)实现了***启动初始化阶段的自动图像配准过程，该调用是在冷启动、热启动过程中自动加载的；2)实现了定时调用图像配准的功能，该调用是通过***中设置的时钟寄存器溢出后的软件标志自动实现调用的，该时钟寄存器的初值决定了定时图像配准程序调用的间隔时间，一般情况下，该时间设置为300秒/次；3)实现了用户触发方式的调用图像配准程序，该调用是在基于ARM的嵌入式操作***的功能界面上通过一个按钮触发完成的，这通常是在用户发现配准后的简单融合图像出现较大偏差时加以补偿的动作。

在ARM芯片上完成嵌入式控制台程序的启动、在线运行及业务逻辑的控制：嵌入式控制台程序中包括了“载入”、“图像A”、“图像B”、“图像配准”、“选项”等功能按钮以及图像数据的显示界面，该显示功能可以使用的支撑平台(包括但不限于)：OpenGL、DirectX等；“载入”按钮主要完成图像数据的自动载入显示过程；“图像配准”按钮实现前述的图像配准程序的调用；“选项”具体给出的是各种图像配准参数以及扩展的图像配准方法的设定，本发明仅给出了“基于Delaunay三角剖分的图像配准方法”；图像数据的显示界面仅是调用了OpenGL、DirectX等支撑平台的API函数实现了视图窗口，并在该视图窗口中显示配准后的图像；“图像A”、“图像B”按钮分别是在视图窗口中显示两异类传感器的动态图像。

基于三角剖分的图像配准方法的DSP模块上进行了视频图像的采集和输入、视频图像的去噪以及视频图像的输出。其具体描述如下：

图像的采集和输入：视频图像的获取是通过A/D转换器，首先将从CCD摄像仪(F-032B/C)和红外热像仪(IR-160)得到的模拟视频信号转换成为LDVS格式的24位数字像素值。将所有的CCD获得模拟信号所表示的像素经过A/D转换后，终将得到一帧一帧的数字视频图像信号，由于F-032B/C的分辨率是640×480，而IR-160的分辨率是160×120，因此需要在DSP侧对CCD信号进行重采样处理，使之与红外传感器得到的160×120能够完全匹配，这里需要说明的是，经过手动调整，是能够使得两个摄像头得到的视频图像达到完全配准的，然后根据这个角度对其进行机械固定。得到数字视频图像后，将其发送到DSP装置进行基于多分辨率方法的去噪处理。

数字图像信号的每一帧在DSP装置上进行基于多分辨率的数字图像信号的去噪处理，该处理过程分为基于多分辨率的分解，这里采用基于小波框架变换的方法：

小波框架的滤波算法程序预定义部分中的小波框架系数的设定所使用的原始小波为DB3小波，其高通滤波系数为：

[0.037828；-0.023849；-0.110624；0.377403；0.852699；0.377404；-0.110624；-0.023849；0.037828]

低通滤波系数为：

[-0.064539；0.040689；0.418092；-0.788486；0.418092；0.040689；-0.064539；0；0]；

对小波系数进行框架化后得到的系数为

其高通滤波系数为：

h1＝[0.037828；0；-0.023849；0；-0.110624；0；0.377403；0；0.852699；0；0.377404；0；-0.110624；0；-0.023849；0；0.037828]/

低通滤波系数为：

g1＝[-0.064539；0；0.040689；0；0.418092；0；-0.788486；0；0.418092；0；0.040689；0；-0.064539；0；0；0；0]/

以上仅为一次分解时的小波框架系数，如果需要更高层的分解，则需要继续进行插零处理。

分解过程是原始图像(或每层的低通图像)与高通滤波器h_n、低通滤波器g_n分别进行滤波，滤波后的结果分别为小波框架在n层的高通系数和n层的低通系数。如果需要继续分解，则使用n层的低通系数再次进行这部分运算，依次类推，直到完成既定的分解目标。

对于小波框架的低通部分和高通部分，后续分别进行了不同的处理。对于高通部分，本算法进行了滤波处理，而低通部分用于进行基于Delaunay的配准策略。

基于小波框架的滤波算法的具体步骤为：首先对原图像进行小波框架的多分辨率分解，建立基于混合高斯模型。

其中混合高斯模型为： $P\left(\mathrm{\φ}\right|\mathrm{\ψ})=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{w}_i{g}_i\left(\mathrm{\φ}\right)---\left(1\right)$

其中w_i是混合权重，且 $\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{w}_i=1,$ g_i(·)是混合元素，其中i＝1，...，M。每一个元素的分布都为一高斯分布，其均值为μ_i，方差为∑_i：

$g_i\left(\mathrm{\φ}\right)=\frac{1}{{\left(2\mathrm{\π}\right)}^{d/2}{\left|{\mathrm{\Σ}}_i\right|}^{1/2}}\exp \{-\frac{1}{2}{(\mathrm{\φ}-{\mathrm{\μ}}_i)}^T{\mathrm{\Σ}}_i^{-1}(\mathrm{\φ}-{\mathrm{\μ}}_i)\}---\left(2\right)$ 其中d是采样数据的维数。

本实施例中使用了5×5窗口的系数对混合高斯模型进行估计，该估计过程是一个不完全数据的估计问题，本发明采用两次迭代的方法分别求取两个最优值：

第一步计算条件概率密度：

$g_{k,i,l}\left[I_i\left(l\right)\right]=\frac{\frac{{\mathrm{\λ}}_{k,i}}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{k,i}^2}}\exp (-\frac{{(I_i\left(l\right)-{\mathrm{\α}}_{i}s\left(l\right)-{\mathrm{\β}}_i)}^2}{2{\mathrm{\σ}}_{k,i}^2})}{\underset{p=1}{\overset{K}{\mathrm{\Σ}}}\frac{{\mathrm{\λ}}_{p,i}}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{p,i}^2}}\exp (-\frac{{(I_i\left(l\right)-{\mathrm{\α}}_{i}s\left(l\right)-{\mathrm{\β}}_i)}^2}{2{\mathrm{\σ}}_{p,i}^2})}---\left(3\right)$

这一步是计算在参数以及原始数据(图像)给定情况下，混合高斯分布第k项的条件概率密度。

第二步更新参数α_i，α_i的值在{-1，0，1}中选择以使(5)式的值最大：

$Q=-\frac{1}{2}\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[\ln \left({\mathrm{\σ}}_{k,i}^2\right)+\frac{{(I_i\left(l\right)-{{\mathrm{\α}}^{\′}}_{i}S\left(l\right)-{\mathrm{\β}}_i)}^2}{2{\mathrm{\σ}}_{k,i}^2}]\·g_{k,i,l}\left[I_i\left(l\right)\right]---\left(4\right)$

第三步重新计算条件概率密度分布g_k，i，l，更新真实场景S(l)：

$S^{\′}\left(l\right)=\frac{\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}(I_i\left(l\right)-{\mathrm{\β}}_i){{\mathrm{\α}}^{\′}}_i\frac{g_{k,i,l}\left(I_i\left(l\right)\right)}{{\mathrm{\σ}}_{k,i}^2}}{\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{{\mathrm{\α}}^{\′}}_i^2\frac{g_{k,i,l}\left(I_i\left(l\right)\right)}{{\mathrm{\σ}}_{k,i}^2}}---\left(5\right)$

${\mathrm{\β}}_i=\frac{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}(I_i\left(l\right)-{\mathrm{\α}}_i^{\′}{S}^{\′}\left(l\right))\frac{g_{k,i,l}\left(I_i\left(l\right)\right)}{{\mathrm{\σ}}_{k,i}^2}}{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\frac{g_{k,i,l}\left(I_i\left(l\right)\right)}{{\mathrm{\σ}}_{k,i}^2}}---\left(6\right)$

第四步重新计算条件概率密度分布g_k，i，l、真实场景S(l)和传感器偏移量β_i，更新噪声的模型参数λ_k，i、σ_k，i ²和β_i：

${{\mathrm{\λ}}^{\′}}_{k,i}=\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{g}_{k,i,l}\left(I_i\left(l\right)\right),k=1,K,K;i=1,K,q---\left(7\right)$

${{\mathrm{\σ}}^{\′}}_{k,i}^2=\frac{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{(I_i\left(l\right)-{{\mathrm{\α}}^{\′}}_{i}s\left(l\right)-{\mathrm{\β}}_i)}^2{g}_{k,i,l}\left(I_i\left(l\right)\right)}{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{g}_{k,i,l}\left(I_i\left(l\right)\right)},k=1,K,K;i=1,K,q---\left(8\right)$

第五步，使用新参数重复第一到第四步。

上述叠代过程中，当参数收敛到某一个特定的范围，确定融合图像为：

$S^{\′}\left(l\right)=\frac{\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{I}_i\left(l\right){{\mathrm{\α}}^{\′}}_i^2\frac{g_{k,i,l}\left(I_i\left(l\right)\right)\·C_i\left(l\right)}{{\mathrm{\σ}}_{k,i}^2}}{\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{{\mathrm{\α}}^{\′}}_i^2\frac{g_{k,i,l}\left(I_i\left(l\right)\right)}{{\mathrm{\σ}}_{k,i}^2}}---\left(9\right)$

其中C_i是边缘特征影响因子， $\underset{i=1}{\overset{q}{\mathrm{\Σ}}}{C}_i=1.$

叠代求出公式(1)和公式(2)中的w_i，期望值μ_i以及方差∑_i。小波框架系数由期望值μ_i进行替换。

基于三角剖分的图像配准方法的DSP模块上执行累积分割检测特征提取技术，其具体描述如下：

利用小波框架分解后的低通系数。检测方法通过扫描整个低通系数图像的每个象素点完成的，以每个象素点为圆心，以某个事先设定好的r为半径作一个图形学的圆，并在图像的具***置上实现定位，取落在图形学圆上的象素点与圆心处的象素点进行比较，如果所有的图形学圆上的象素点大于或小于圆心出象素点的灰度值到达一个级别t，则认为该点为特征点；在该实施例中，本发明使用了3作为半径，这意味着在图形学圆上有12个象素。经过大量的实验验证，9个象素值应该是最佳技术个数，也就是说这12个象素中有9个象素比圆心处的象素的灰度值大于t，则可认为是特征点。

基于三角剖分的图像配准方法的DSP模块上执行随机增量法，其具体描述如下：

按照随机的顺序依次***点集中的点，在整个过程中都要维护并更新一个与当前点集对应的Delaunay三角剖分。本发明采用简化算法实现Delaunay三角剖分，其具体步骤为：1、初始化：构造一个超级三角形，包含所有散点，放入三角形链表。2、将点集中的散点依次***，在三角形链表中找出其外接圆包含***点的三角形(称为该点的影响三角形)，删除影响三角形的公共边，将***点同影响三角形的全部顶点连接起来，从而完成一个点在Delaunay三角形链表中的***。3、根据优化准则对局部新形成的三角形进行优化。将形成的三角形放入Delaunay三角形链表。4、循环执行上述第2步，直到所有散点***完毕。

三角最优匹配程序使用了两个指标对叠代过程中的Delaunay三角形进行比对，当比较后的结果的总和小于一个值时，则可以选择相应的旋转角度参数和平移参数作为最终的配准结果；本发明中匹配过程所使用的优化算法采用了自适应模拟退火的算法。

这两个指标分别表示了两个特征不变量，假设l₁，l₂，l₃为三角形的三条边，m₁，m₂，m₃为三角形的三个角：

第一个刚性变换不变量：

$V_t=[\frac{l_1}{l_3},\frac{l_2}{l_3},\cos \left(m_3\right)]$ (公式1)

第二个刚性变换不变量：

V_m＝[∠m₁，∠m₂，∠m₃](公式2)

本发明中红外图像的小波框架低通系数中存在24个Delaunay三角，在可见光图像的低通系数中存在18个Delaunay三角。这样一共有24×18(432)个匹配结果，这些匹配结果都进行两个刚性特征不变量进行计算，然后对他们之间的比值进行比较，选取比值最大的那一组匹配。匹配后的图像在小波框架领域内进行简单的象素级融合，本发明采用了基于传统小波的融合方法，以表示配准后的效果。

本发明方法的整个图像融合过程中需要用到两幅视频图像的输入，如图3和图4为CCD图像和红外图像，图5为最终的配准结果。

Claims (2)

1.一种基于三角剖分的图像配准方法，其特征在于，方法包括如下步骤：

1)视频图像的获取：将从CCD摄像仪和红外热像仪得到的模拟视频信号转换成为LDVS格式的24位数字像素值，所有的模拟视频信号所表示的像素经过A/D转换后，得到一帧一帧的数字视频图像信号，将其发送到DSP装置进行基于多分辨率方法的去噪处理；

2)基于多分辨率的分解：采用基于小波框架变换的方法，原始图像与高通滤波器h_n、低通滤波器g_n分别进行滤波，滤波后的结果分别为小波框架在n层的高通系数和n层的低通系数，如果需要继续分解，则使用n层的低通系数再次进行这部分运算，依次类推，直到完成既定的分解目标；对于小波框架的高通部分，进行了滤波处理，而低通部分用于进行基于Delaunay的配准策略；

3)角点检测：利用小波框架分解后的低通系数；通过扫描整个低通系数图像的每个像素点完成的，具体为以每个像素点为圆心，以某个事先设定好的r为半径作一个图形学的圆，并在图像的具***置上实现定位，取落在图形学圆上的像素点与圆心处的像素点进行比较，如果所有的图形学圆上的像素点大于或小于圆心处像素点的灰度值到达一个级别t，则认为该点为特征点；

4)三角形匹配：构造一个超级三角形，包含所有特征点，放入三角形链表；将点集中的特征点依次***，在三角形链表中找出其外接圆包含***点的三角形，删除影响三角形的公共边，将***点同影响三角形的全部顶点连接起来，从而完成一个点在Delaunay三角形链表中的***；根据优化准则对局部新形成的三角形进行优化，将形成的三角形放入Delaunay三角形链表；循环执行特征点***，直到所有散点***完毕；

5)配准图像：Delaunay三角形链表中的匹配使用了两个刚性变换不变量进行计算，假设l₁，l₂，l₃为三角形的三条边，m₁，m₂，m₃为三角形的三个角：

第一个刚性变换不变量：

第二个刚性变换不变量：V_m＝[∠m₁，∠m₂，∠m₃]，

然后对他们之间的比值进行比较，当比较后的结果的总和小于一个值时，则可以选择相应的旋转角度参数和平移参数作为最终的配准结果。

2.根据权利要求1所述的基于三角剖分的图像配准方法，其特征在于，所述步骤2)中基于小波框架变换的方法：小波框架的滤波算法程序预定义部分中的小波框架系数的设定所使用的原始小波为DB3小波，其高通滤波系数为：

[0.037828；-0.023849；-0.110624；0.377403；0.852699；0.377404；-0.110624；-0.023849；0.037828]

低通滤波系数为：

[-0.064539；0.040689；0.418092；-0.788486；0.418092；0.040689；-0.064539；0；0]；

对小波系数进行框架化后得到的系数为：

高通滤波系数

低通滤波系数

以上仅为一次分解时的小波框架系数，如果需要更高层的分解，则使用n层的低通系数再次进行这部分运算，依次类推，直到完成既定的分解目标。

Images