基于傅立叶梅林变换和最大互信息理论的图像配准方法

摘要

本发明公开了一种基于傅立叶变换和powell算法的图像配准方法，主要解决了基于傅立叶梅林变换配准精度不高和基于powell对初始值敏感的问题提出的。其实现步骤为：(1)输入浮动图像和参考图像(2)低通滤波预处理(3)利用傅立叶梅林变换得到变换参数(4)计算互信息值(5)依据互信息和powell判断是否最优(6)若不是，则继续搜索，搜索时不断重复“空间几何变换—>统计联合直方图—>计算互信息—>最优化判断”的过程(7)搜索到满足精度要求的参数(8)几何变换(9)图像融合。本发明克服了Powell优化算法对初始值敏感的缺点，得到了很好的配准精度。

说明书

技术领域

发明属于图像处理技术领域，更进一步涉及图像配准领域的一种先进行粗匹 配得到初始参数的一种高精度的配准方法。本方法可以用于经过仿射变换的医学 图像和雷达孔径合成SAR图像的配准。

背景技术

近年来，图像配准应用于各种目标识别和图像融合领域，而用进化算法应用 于图像配准领域是一个热点问题。从配准过程来看，最重要的部分就是找到精确 合适的配准参数，利用配准参数对图像进行合适的空间变换，从而得到正确配准 结果。

目前，按照配准过程来说，可以分为基于特征的图像配准和基于灰度的图像 配准，他们的主要区别在于是否包含分割步骤。基于特征的方法包含图像的分割 过程，用于提取图像的特征信息然后进行显著特征的配准，基于灰度的配准方法 无需进行图像的分割与特征提取，直接用图像的统计信息作为配准的相似性度 量。

1994年Chen等人提出了一种基于傅立叶梅林(Fourier-Mellin)变换的图像 配准方法。该方法假设的是两幅图像经过了旋转尺度伸缩和平移变换后的图像。 该方法通过计算两幅图像功率谱的反Fourier变换所对应的峰值位置求取他 们的相对位移，通过对图像幅度谱进行对数-极(Log-Polar)变换，在对数-极空 间用类似于计算相对平移的方法求得旋转角度。大量实验结果发现，对于配准精 度要求不高的情况下，该方法具有很强的鲁棒性，但对于需要进一步配准的精度 的情况下，该方法得到的结果只能作为配准的初始估计值。

通常，浮动图像相对于参考图像是一种仿射变换。当确定了浮动图像与参考 图像的配准参数后，空间变换模型就建立了。对浮动图像进行空间逆变换，就可 得到配准图像。通过计算配准图像和参考图像的联合直方图得到互信息值。互信 息值可通过两幅图像的联合直方图求出，而联合直方图采用部分体积插值法

(PV)求得，而后利用Powell搜索配准参数，由于Powell算法具有依赖初始点 的缺陷，会导致搜索陷入局部极值，或者导致配准失败。

发明内容：

本发明意在克服上述方法的缺点，提出将一种将傅立叶梅林(Fourier-Mellin) 变换和最大互信息理论结合的图像配准方法。

为实现上述目的，本发明具体实现步骤为：

(1)输入两幅待配准图像一幅参考图像一幅浮动图像

(2)用傅立叶变换位移理论求出浮动图像与参考图像的相对平移量

3a)计算两幅图像的离散傅立叶变换；

3b)再通过两幅图像的互功率谱计算相位差；

3c)求互功率谱的逆变换，得到一个脉冲函数，脉冲峰值位置即为相 对平移量Xtran和Ytran。

(3)用傅立叶梅林变换估计浮动图像相对于参考图像的旋转角度

2a)通过二维离散傅立叶变换求出图像幅度谱；

2b)在对数-极(Log-Polar)空间应用傅立叶梅林变换求出旋转角度 theta和缩放因子delta。

(4)对两幅原图像通过低通滤波器进行预处理

(5)设置初始点和初始搜索方向，初始点即设置为(2)(3)步估计的相 对平移量和旋转角度，搜索方向为单位方向

(6)根据上述参数对浮动图像进行空间几何变换(affine)

(7)计算原参考图像与变换后图像的互信息值

7a)使用PV(Partial Volume)插值法计算两幅图像的联合直方图；

7b)求出两幅图像的联合概率分布函数和边缘分布函数；

7c)利用概率分布函数及边缘分布函数求得熵，从而求出归一化互信 息值。

(8)因为互信息关于这三个参数没有具体的表达式，所以无法利用传统优 化算法，这里用方向加速法(Powell)通过直接比较归一化互信息数 值大小移动迭代点求最大值处的三个变换参数值

(9)若达到最大值，则输出配准参数，若没有则返回步骤(6)

(10)对最优配准参数变换后的图像和原参考图像进行图像融合，获得配 准融合图像

本发明首先利用两幅图像功率谱的反Fourier变换所对应的峰值位置求取两 幅图像的相对位移，通过对图像幅度谱进行对数-极(Log-Polar)坐标变换，在 对数-极空间用类似于计算相对平移的方法求得到旋转角度，然后利用求得的旋 转角度和平移位移来建立空间变换模型，以最大互信息作为相似性度量，通过 Powell算法进行寻优，得到最佳配准参数，对浮动图像进行空间逆变换，最后将 变换后的浮动图像和原参考图像进行图像融合，就可得到配准融合图像。与现有 技术相比，本发明具有以下优点：

第一：本发明在应用傅立叶梅林(Fourier-Mellin)变换进行图像配准之后再 应用基于最大互信息理论的Powell优化算法进行精度匹配，使得配准精度大大 提高，克服了傅立叶梅林变换对于高精度要求的配准不具有鲁棒性的缺点。

第二：本发明在应用Powell算法时，由于在基本的Powell算法中每次都会 无条件的使用新的搜索方向代替原来旧的搜索方向，而本发明在替换搜索方向时 会考虑线性无关的问题，当搜索方向线性无关时能够保证每轮迭代中以搜索方向 为列的列式不为零，因此这些方向是线性无关的，而且随着迭代的增加搜索方向 的共轭程度逐渐增强。

第三：本发明应用了傅立叶梅林(Fourier-Mellin)变换配准的结果参数作为 Powell搜索的初始值，克服了Powell算法对初始值敏感的缺点，配准效果不佳 的缺点，使得配准精度大大提高。

附图说明：

图1为发明的流程图

图2为本发明应用于合成孔径雷达SAR图像的配准结果。

图3为本发明应用于合成孔径雷达SAR图像的配准结果。

图4为本发明应用于医学图像的配准结果。

具体实施方式：

下面结合附图1对本发明的步骤做进一步描述。

步骤1，输入两幅图像，一幅参考图像和一幅浮动图像，这两幅图像是具有 仿射变换关系的两幅图像。

步骤2，对这两幅图像分别进行低通滤波预处理。

低通滤波处理意在对图像进行简单的去噪处理。

步骤3，用傅立叶变换位移求出浮动图像与参考图像的相对平移量。

设f₂(x,y)、f₁(x,y)为两幅图像，f₂(x,y)为f₁(x,y)在x和y方向分别 平移x₀和y₀后的图像，即

f₂(x,y)＝f₁(x-x₀,y-y₀)   (1)

若f₁和f₂对应的傅立叶变换分别为F₁(u,v)和F₂(u,v)，则他们之间有如下关 系：$F_2(u,v)=F_1(u,v)e^{-j({\mathrm{ux}}_0+{\mathrm{vy}}_0)}---\left(2\right)$

f₁(x,y)和f₂(x,y)的互功率谱设为G(u,v)，则

$G(u,v)=\frac{F_1(u,v)F_2^{*}(u,v)}{\left|F_1(u,v)F_2^{*}(u,v)\right|}=e^{j({\mathrm{ux}}_0+{\mathrm{vy}}_0)}---\left(3\right)$

其中F₂^*表示F₂的复共轭。

可以看出，互功率谱的相位等价于图像间的相位差，通过对式(3)进行傅 立叶逆变换，在(x,y)空间将形成一个脉冲函数，该脉冲函数在其他各处的函数 值均为零，只有在脉冲峰值处的函数值不为零，脉冲峰值位置即为参考图像和浮 动图像之间的平移量x₀和y₀。

因此，本发明可以通过对式(3)进行傅立叶逆变换，在(x,y)空间将形成一 个脉冲函数，找到脉冲峰值位置，该脉冲峰值位置即为浮动图像与参考图像的相 对平移量x₀和y₀。

步骤4，用傅立叶梅林变换估计浮动图像的旋转角度和缩放因子。

考虑被配准的两幅图像s(x,y)和r(x,y)，其中r(x,y)表示浮动图像，s(x,y) 表示参考图像，其中s(x,y)是r(x,y)经过平移、旋转和尺度缩放变换后的图像， 即：

s(x,y)＝r[(σ(xcosα+ysinα)-x₀，σ(-xsinα+ycosα)-y0)]   (4)

式中，x₀和y₀为浮动图像相对参考图像的平移量，α为浮动图像相对参考图像的 旋转角度，σ为缩放因子。

那么s(x,y)和r(x,y)对应的傅立叶Fourier变换S(u,v)和R(u,v)之间满足：

|S(u,v)|＝σ^-2R[σ^-1(ucosα+vsinα),σ^-1(ucosα+vsinα)]|   (5)

其中|·|表示求幅度谱，从上式可以看出，旋转角度α和缩放因子σ可以和平移量 进行分离计算，频谱幅度仅与旋转角度和缩放因子有关，而与平移量无关，因此， 可以通过图像幅度谱求出旋转角度和缩放因子。

定义：r_p1(θ,logρ)＝r_p(θ,ρ)   (6)

s_p1(θ,logρ)＝s_p(θ,ρ)   (7)

其中r_p和s_p分别是r(x,y)和s(x,y)在极坐标系(θ，ρ)中的幅度谱，即是 r_p＝|R(u,v)|，s_p＝|S(u,v)|，其中θ＝arctan(v/u)，α为旋转 角度，σ为缩放因子，那么很容易的得出

s_p1(θ,logρ)＝r_p1(θ-α,logρ-logσ)   (8)

或者s_p1(θ,λ)＝r_p1(θ-α,λ-κ)   (9)

其中λ＝logρ,κ＝logσ

可以看出，通过上述变换，式(9)变换为与式(1)相同的形式，这样就在 对数-极(Log-Polar)空间应用Fourier变换位移，按照式(2)和式(3)求得α和κ。 如果对数的底取为e，那么

σ＝e^κ   (10)

这样就求出了旋转角度α和比例因子σ，根据求出的α和σ对图像s(x,y)进 行反变换得到图像s₁(x,y)，然后通过式(2)和式(3)在图像空间计算出s₁(x,y) 和r(x,y)之间的平移量x₀和y₀。式(9)被称为傅立叶梅林(Fourier-Mellin)变 换。

步骤5，设置初始点和初始方向，初始方向即为上述步骤中通过傅立叶变换 得到的相对平移量x₀和y₀和旋转角度α。这三个参数即为初始配准参数。

步骤6，对浮动图像进行空间几何变换。

在图像配准过程中，空间变换模型采用仿射变换描述。对于两幅图像，首先 建立统一的坐标系，定义图像的中心为坐标原点，这与通常定义的图像的左上角 为坐标原点不同，所以需要对浮动图像进行坐标平移，浮动图像要经过水平移动、 垂直移动和绕原点旋转的变换去和参考图像进行匹配，那么浮动图像的任意点 O_F到参考图像的坐标点O_R的仿射变换可以描述为：

O_R＝t₁(c,d)×t₂(θ)×t₃(a,b)×t₄(tx,ty)×O_F   (11)

其中，O_F为浮动图像中任一点坐标，O_R是参考图像中任一点坐标。O_F经过仿射 变换到R坐标系中。t₁(c,d)、t₂(θ)、t₃(a,b)、t₄(tx,ty)都是3×3的矩阵。 t₁(c,d)是将浮动图像原点从图像左上角移动到图像中心，移动位移为(c,d)； t₂(θ)将浮动图像绕着新原点(图像中心)旋转，旋转角度是θ；t₃(a,b)是将旋转 后的图像再将坐标原点平移回左上角，移动位移为(a,b)；t₄(tx,ty)是浮动图像 沿着水平、垂直方向分别作tx,ty平移。由于浮动图像和参考图像分辨率相同， 所以t₁(c,d)和t₃(a,b)中的坐标原点平移参数相同，且图像尺寸一旦确定，则平 移参数的值也确定。

因此就剩下寻找使得归一化信息达到最大值时的tx、ty、θ三个参数，即: f_max＝maxf(tx,ty,θ)，属于多参数优化问题。

根据设置的初始参数对浮动图像做空间几何变化得到变换后图像。

步骤7，使用PV(Partial Volume)插值法计算联合直方图，从而计算原参考 图像与变换后图像的互信息值。

参考图像R和浮动图像F的概率分布函数是统计各自的灰度值在整幅图像中 出现的频数，联合概率分布函数是统计两幅图像在相同位置上一对灰度值出现的 频数。本文采用PV(Partial Volume)插值法求归一化的联合直方图，再求出联 合概率分布函数。

设浮动图像F上任意一点的坐标b(i,j)经过反向变换到参考图像的坐标系 内，得到新点T_r(b)，其邻近像素点分别是n₁,n₂,n₃,n₄。为了避免参考图像引入 原本不存在的灰度值，所以采用PV插值利用T_r(b)邻近的四个点计算联合灰度直 方图，即修改联合直方图中四个灰度值对：

h(F(b),R(n₁))、h(F(b),R(n₂))、h(F(b),R(n₃))、h(F(b),R(n₄))。各个修正量由图 1中ω₁到ω₄所示的面积确定，且四点增量之和为1，PV插值法更重要的一点是， 联合直方图中每个灰度对的值是许多小幅度变化的小数的和，不再以1增加，从 而得到光滑的目标函数，有利于优化搜索。

当浮动图像与参考图像完全配准时，即没有任何平移或者旋转，图像的联合 直方图的点集中在对角线上，此时归一化互信息值最大；当浮动图像沿着垂直方 向依次移动的话，其联合直方图的点越来越偏离对角线，分散的分布在对角线周 围。

两幅图像的最大互信息计算公式为：I(R,F)＝H(R)+H(F)-H(R,F)， Studholme等人的研究表明最大互信息对待配准的两幅图像间的重叠比较敏感， 因此对最大互信息公式进行简单扩展，得到归一化互信息计算公式：

$(R,F)=\frac{H\left(R\right)+H\left(F\right)}{H(R,F)}---\left(12\right)$

归一化互信息作为两幅图像的相似性测度能更好的反映配准的变化，减少对图像 重叠区域的敏感度，提高配准精度。上式中H(R)是参考图像的熵，H(F)是浮动 图像的熵；H(R，F)是参考图像和浮动图像的联合熵、需要先求出联合概率分布 函数，然后再分别求出关于R，F的边缘分布函数，最后再求出各自的熵，带入 归一化互信息公式中求得归一化互信息值。

步骤8，根据最大互信息理论利用Powell算法寻优。

(1)给定允许误差ε＞0初始点X⁰和n个线性无关的方向

d^(1,1)，d^(1,2)…d^(1,n)置k＝1。

(2)置x^(k,0)＝x^(k-1)，从x^(k,0)出发，依次沿方向d^(k,1)，d^(k,2)…d^(k,n) 进行一维搜索，得到点x^(k,1)，x^(k,2),…,x^(k,n)，求m,使得 f(x^(k,m-1))-f(x^(k,m))＝max{f(x^(k,j-1))-f(x^(k,j))} 令d^(k,n+1)＝x^(k,n)-x^(k,0)若||x^(k,n)-x^(k,0)||≤ε则停止计算；否 则，进行第(3)步。

(3)求λ_n+1使得f(x^(k,0))+λ_n+1d^(k,n+1)＝minf(x^(k,0)+λd^(k,n+1))令 x^(k+1,0)＝x^(k)＝x^(k,0)+λ_n+1d^(k,n+1)若||x^(k)-x^(k-1)||≤ε 则停止计算，得到x^(k)否则进行步骤(4)

(4)若$\left|{\lambda }_{n+1}\right|>[\frac{f(x^{(k,0)}-f\left(x^{(k+\mathrm{1,0})}\right))}{f\left(x^{(k,m-1)}\right)-f\left(x^{(k,m)}\right)}{]}^{\frac{1}{2}}$则令

d^(k+1,j)＝d^(k,j),j＝1,…,m-1

d^(k+1,j)＝d^(k,j+1),j＝m,…,n

k＝k+1

否则，令$\left(\begin{array}{c}{d}^{(k+1,j)}=d^{(k,j)},j=1,...,n\\ k=k+1\end{array}\right)$转到步骤(2)

由于互信息关于这三个参数没有具体的函数表达式，所以无法利用传统的优 化算法，因此，本文采用方向加速法(Powell)通过直接比较归一化互信息数值 的大小移动迭代点求出最大值处的三个参数值。利用上述Powell寻优方法对配 准参数进行迭代寻优，其中一维搜索运用黄金分割法。

步骤9，若达到最大值，则输出配准参数，若没有则返回步骤6。

步骤10，对应用最优配准参数变换后的浮动图像和原参考图像进行图像融 合，获得配准融合图像。

图像融合过程中，为了可以清晰的看到配准的结果与误差。我们采用将配准 结果图中的灰度重新进行赋值处理，参考图像的灰度值总体赋值成原来的百分之 三十，变换后的浮动图像的灰度值赋值成原来的百分之七十。若是灰度值为255 或者0，则直接赋值为原灰度值。这样一来，可以使得配准的结果明显清晰的显 示，配准误差也可以直观看出。

本发明的效果可以通过下列进一步仿真说明：

本发明图像配准方法对两组合成孔径雷达SAR图像和一组医学图像进行配 准仿真，从配准结果可以直观看出配准精度的提高和以及良好的效果。

对合成孔径雷达SAR图像的仿真结果如图2和图3所示，可以直观看出将 仿射变换过的浮动图像通过了旋转和平移后与参考图像准确的配准了，并且完整 的融合。

对一组医学图像仿真结果如图4所示，可以直观看出，其配准图像和参考图 像完美的配准并融合了。

通过上述实验结果仿真，验证了本发明的可行性与准确性。