基于区域互信息的特征级图像配准算法

基于区域互信息的特征级图像配准算法摘要利用互信息法进行图像配准已成为图像处理领域的热点。图像融合的关键是图像配准。图像配准的方法主要有特征点法、曲线法、表面法和矩主轴法等。在此过程中不需要进行图像分割以及特征提取，可以实现图像配准的自动化，而且鲁棒性较强，配准精度高，但是，计算互信息相似度是基于整幅图像的像素灰度，因此计算复杂度较高。最大互信息配准法由于不需要对不同成像模式下的图像灰度间的关系作任何假设，也不需要对图像进行分割或任何预处理，具有自动化程度高的特点。近年来受到了越来越多学者的关注，并且在图像配准领域得到了普遍重视和广泛应用，并被许多图像处理软件包作为标准的配准算法。关键词：区域互信息;特征级图像;配准;算法 目录毕业论文（设计） I基于区域互信息的特征级图像配准算法 I1. 绪论 62. 研究综述 62.1任务背景 62.2过程分析 62.3系统开发技术 83. 基于互信息的图像配准 143.1总体设计思路和图像配准实现 143.2几何变换 153.3互信息计算 153.4灰度级插值技术 183.5优化搜索算法 183.6一维搜索算法 194. 测试和运行 204.1实验环境 204.2环境配置 204.3运行结果 215. 结论 22参考文献 23参考文献 25致谢 26绪论随着现代科技的迅速发展和新型传感器的不断涌现，人们获取图像数据的能力不断提高。由于不同物理特性的传感器所产生的图像不断增多，同一地区往往可以获得大量不同尺度、不同光谱、不同时相的多源图像数据信息。在利用多源图像信息进行数据融合、目标变化检测、目标识别等多源协同处理工作之前，必须进行多源图像配准工作，配准精度的高低直接影响到后续应用效果的好坏。所以，如何对图像进行高精度的配准一直是图像处理领域的热点与重点。医学图像的配准经过了一个较长的发展历程。二十世纪八十年代初，图像配准主要应用在DSA(数字减影血管造影)方面。它针对二维图像采用基于图像灰度的方法，通过检测相关性和灰度值的差异来决定刚性变换的参数。这种方法在今天的三维变换中仍可以借鉴。二十世纪八十年代中后期到九十年代初，不同模态的图像配准开始受到关注，但这时候图像配准和图像融合也主要集中于二维刚性变换。二十世纪九十年代，医学图像配准研究取得了飞速发展，尤其在三维方向和非刚性图形配准上成果非常显著。由此可见，不论在国内还是国外图像配准技术经过多年的研究，图像配准技术己经取得了诸多研究成果。研究综述2.1任务背景图像配准技术从基于特征的配准方法发展到基于统计的配准方法有其突破性的意义。与基于特征的配准方法相比，基于统计的配准方法的突出优点为鲁棒性好、配准精度高、人工干预少。基于统计的配准方法通常是指最大互信息的图像配准方法。基于互信息的图像配准是用两幅图像的联合概率分布与完全独立时的概率分布的广义距离来估计互信息，并作为图像配准的测度。当两幅图像达到最佳配准时，它们的对应像素的灰度互信息应为最大。由于基于互信息的配准对噪声比较敏感，首先，通过滤波和分割等方法对图像进行预处理。然后进行采样、变换、插值、优化从而达到配准的目的。基于互信息的配准技术属于基于像素相似性的方法。它基于图像中所有的像素进行配准，基于互信息的图像配准引入了信息论中的概念，如熵、边缘熵、联合熵和互信息等，可使配准精度达到亚像素级的高精度。2.2过程分析根据任务的功能及目标，图像配准算法可分为以下几步：特征空间、搜索空间、搜索策略与相似性测度四部分组合成医学图像配准技术。如图2.1所示。图2.1算法过程图图像配准的基本框架包括以下四个方面：特征空间、搜索空间、搜索策略和相似性度量。特征空间特种空间的选取对图像配准有着重要的意义，因为特征空间不仅直接关系到图像中的那些特征对配准算法敏感和哪些特征被匹配，而且大体上决定了配准算法的运行速度和鲁棒性性能。在图像的配准中究竟采用何种特征；图像配准算法的运行速度及执行效率；算法对何种特征敏感；配准之后的图像的性能包括鲁棒性，高精度性能等。由此看来在图像配准过程中如果图像特征选择得当，就能够提高图像配准的执行效率，加快运行速度，同时又能够使图像的信息得到充分的表达。2.搜索空间搜索空间是指在配准过程中对图像进程变换的范围及变换的方式。经过上一步的图像特征提取之后，图像配准的下一步操作就是依据两幅图3.搜索策略搜索策略的任务是在搜索空间中找到最优的配准参数，在搜索过程中以相似度量的值作为判优依据。4.相似性度量相似性度量是衡量每次变换结果优劣的准侧，用来对变换结果进行评估，为搜索策略的下一步动作提供依据。相似性度量和特征空间、搜索空间紧密相关，不同的特征空间往往对应不同的相似性度量；而相似性度量的值将直接决定配准变换的选择，以及判断在当前所取的变换模型下图像是否被正确匹配了。通常匹配算法抗干扰的能力是由特征提取和相似性度量共同决定的。常用的相似性度量有互信息、归一化互信息、联合熵。而相似性度量的值将直接决定配准变换的选择，以及判断在当前所取的变换2.3系统开发技术互信息是信息论中的一个基本概念，通常用于描述两个系统间的相关性，或者是一个系统中所包含的另一个系统信息的多少，是两个随机变量A和B之间统计相关性的量度,或是一个变量包含另一个变量的信息量的量度。它可以用熵和来描述以及联合熵（2-1）其中和分别是系统A和B的熵，是A,B的联合熵，表示一直系统A时B的条件熵和一直系统B时A的条件熵。上述各种熵可分别表示为（2-2）（2-3）（2-4）（2-5）（2-6）其中，aA,bB,和分别是系统A和系统B完全独立时的边缘概率密度，是系统A和B的联合概率分布，是已知系统B时A的条件概率分布，是已知系统B时A的条件概率分布，如果联合概率分布密度满足,则随机变量A和B相互独立;如果A和B满足某映射关系T使,则随机变量A和B最大相关。基本Powell算法实现思想是：把整个过程分为若干轮迭代，每一轮迭代都有n个搜索方向，进行n+1次一维搜索。在每一轮搜索中都要先确定一个初始点，从初始点出发沿已知的n个搜索方向一次进行一维搜索，得到一个最好的点，接着再沿与连线方向进行一次一维搜索，得到本轮最好点；然后以改点作为初始点进行新一轮迭代。基本的实现算法如下：（1）给定允许误差、初始点和n个线性无关的方向,,,…,（2）从出发依次沿方向,,,…,进行一维搜索，得,,,…,.再从出发沿与连线方向进行一维搜索，得到（3）如果,停止搜索，得到点;否则置，j=1,2,…,n。在实际应用中，图像之间的关系是很复杂的，变换一般也不是单一的形式，但是为了实验研究，并考虑到误差的可允许性，通常将变换模型分为两大类:线性变换和非线性变换。在图像处理中，经常用到的变换形式主要有刚体变换(rigidtransfonnation)、仿射变换(affinetransformation)、投影变换(Projectivetransformation)、透视变换(perectivetransfo而ation)、多项式变换(transfonnation)等。图2一1显示了几种常见的几何变换。(l)刚体变换。如果第一幅图像变换到第二幅图像后，图像中亮点间的距离仍保持不变，简单的说就是物体或镜头的运动不影响成像物体非让相对尺寸及形状，那么这种变换称为刚体变换，它是平移、旋转和缩放的组合，其数学描述为（2-7）其中，是比例缩放因子，是旋转角，和是平移量。一般最少需要两对点来求解变换参数。(2)仿射变换。仿射变换是更复杂一些的图像形变，它在复杂畸变下仍然能保持很好的数学特性。其数学表达式为（2-8）仿射变换是一种常用的配准变换模型，对于在不同位置的相同视角拍摄同一场景目标而成的两幅图像的配准问题，它能够很好地满足配准要求。仿射变换后直线依然映射为直线，仍然保持平衡关系，其包含六个参数，至少需要三对不在一条直线上的控制点来求解变换参数。(3)二维投影数学描述为（2-9）投影变换与仿射变换都可以保证任意方向的直线变换后仍为直线。但是投影变换的九个自由度(对应于九个系数)可以满足平行四边形的映射，这是具有六个自由度的仿射变换所不能做到的。(4)多项式变换的一般形式可表示为公式（2-10）(5)透视变换通常用于描述三维场景到平面的映射。如果一个目标在三维场景中坐标为，那么它投影到平面上的图像坐标(通过透视变换可描述为（2-11）平移变换：（1）原理分析平移变换是指将图像中所有像素都按照给定的平移量水平、垂直方向移动如图所示，设为原图像上的点，图像水平平移量为,垂直平移量为，则平移后点，坐标变换为从图看出和满足下列关系（2-12）用矩阵表示为（2-13）对4-2的矩阵求逆，可以得到反向变换（2-14）其对应的坐标方程为（2-15）通过式2-15，平移后新图像中的任一点都可以在原图像中找到对应的点。例如，对于新图像中的点（0，0），代入方程组，可以求出在原图像中的对应点为当大于0或大于0时，点不在原图像中，也就是说改点的灰度值在原图像中没有定义，此时可以直接将它的灰度值设置为0或255.相反，原图像右侧的部分像素点没有包含在新图像中，也就是说原图像中有点被移出了显示区域。对图像进行平移变换用到函数se=translate(strel(1),[xy])。上图为用函数se=translate(strel(1),[yx])使图像相对于参考图像在x轴方向向左平移20个像素单位，在y轴方向向上平移30个像素。上图为用函数se=translate(strel(1),[yx])使图像相对于参考图像在x轴方向向左平移70个像素单位，在y轴方向向上平移30个像素。旋转变换：（1）原理分析旋转是指将图像中所有的像素按给定的角度在平面内顺时针或者逆时针沿圆弧路径进行变换。通常约定为正值时做逆时针旋转，为负值时做顺时针旋转。当基准点为原点时，原始点和变换后的角度和坐标关系如图所示。其中r是点到原点的距离，是点与x轴的夹角，是旋转角。利用三角函数，可以计算出旋转前后的坐标满足下列关系。旋转前：旋转后：用矩阵表示为（2-16）对2-16求逆矩阵，可以得到反向变换（2-17）上述旋转是以坐标原点（0，0）为基点进行的，但在实际应用中，往往是以图像中心为基点对图像进行旋转。此时的做法是：先将坐标系平移到图像中心，再进行旋转，然后平移回新的坐标原点。下面推导坐标系平移的转换表达式。将坐标系平移到坐标系处，其中坐标系的原点在坐标系中的坐标为，则它们的坐标变换矩阵表达式为：（2-18）其逆变换转换矩阵表达式为：（2-19）如果图像旋转前中心坐标为,旋转后中心坐标为,则旋转变矩阵的表达式为：（2-20）对图像进行旋转变换主要用函数I=imrotate(J,,'nearest','crop')，参数表示往往是以图像中心为基点对图像进行旋转。此时的做法是：先将坐标系平移到图像中心，再进行旋转，然后平移回新的坐标原点。通常约定为正值时做逆时针旋转，为负值时做顺时针旋转。（2-21）由上式可以看出，信息量随着消息的可能性组合m增多而增多，如果消息只有一种可能性时即事件为必然事件时，那么消息中包含的信息量为零。因此可以看出，可能收到的不同消息越多，对收到哪条消息的不确定性就越大;相反，收到只有一种可能性的消息，不确定性为零，Hartley对消息的度量实际是对不确定性的度量。Hartley度量方法的不足之处是式他所定义信息量是假定所有符号发生的概率相同，但实际情况各符号并不一定都等概发生，为此，Shannon用概率加权来衡量消息出现的可能性，对Hartley的度量方法做出改进。设某一随机过程中有k种可能的情况，每种情况发生的概率分别是，，…，，shannon给出了嫡的如下定义：（2-22）当所有可能的事件均以相等的概率发生时，上式就成了Hartley定义的嫡，并且这时嫡取得最大值，即(2-23)；所以，Hartley嫡是shannon嫡的特殊情形，而Shannon更具有一般性。基于互信息的图像配准3.1总体设计思路和图像配准实现图3.1程序流程图本算法从图像配准的框架入手，着重研究了基于最大互信息的配准方法。为了加快配准速度，针对互信息计算量大和存在局部极值的问题，本文集中于优化策略的研究，在一般优化算法的讨论分析基础之上，提出了改进的优化算法，针对一般Powell法不考虑线性无关问题，本设计采用了改进后的POWELL算法，可以使搜索方向线性无关，共轭程度增加。该程序能够实现平移X、Y和旋转角度R的图像之间的配准，有待改进powell优化算法求取更多的参数，使配准出的图像的精确度更加高。该算法首先获取两幅图像：参考图像ct.jpf和待配准图像convert.jpg.。其中convert.jpg图像是通过translateImage.m生成的。然后运行POWELL算法得到最优的变换参数。其中POWELL算法中包括了坐标的变换，互信息的求取函数（利用最大互信息衡量最优参数），插值算法，还调用了一维搜索等算法用于优化搜索算法。最后，对待配准图像进行坐标变换得到结果。利用POWELL的结果作为参数，再次运行translateImage函数得到最后的图像。3.2几何变换本算法设计应用空间几何变换对图像进行平移变换和旋转变换。se=translate(strel(1),[x,y]);J=imdilate(I,se);其中函数imdilate主要用于对图像I进行膨胀操作，结构体se由strel创建。实现该功能的函数文件为translateImage.m。旋转变换是指图像中所有的像素点（基准点除外）按给定的角度在平面内顺时针或逆时针沿圆弧路径进行交换。通常约定角度为正值时做顺时针旋转，角度为负值时做逆时针旋转。但是在本算法中规定角度为正值时做逆时针旋转，角度为负值时做顺时针旋转。主要步骤如下所示：J=imrotate（J，ang，’bilinear’,’crop’）；该函数功能为逆时针旋转图像J；其中J为需要旋转的图像，ang为要旋转的角度，bilinear为插值方法，本次设计中主要应用双线性插值方法。实现该功能的函数文件为translateImage.m。3.3互信息计算互信息是信息论中的一个基本概念，通常用于描述两个系统间的统计相关性，或者是一个系统中所包含的另一个系统中信息的多少，它可以用信息熵来描述：I(A,B)=H(A)+H(B)-H(A,B)。其中，H(A)和H(B)分别是系统A和B的熵，H(A，B)是它们的联合熵，实现该功能的函数利用一个联合直方图，计算出两幅图像的互信息值。首先建立三个矩阵zeros(256，256)；zeros(1，256)；zeros(1，256)；再进行归一化处理，然后统计出需要计算出互信息的两幅图像的直方图和他们的联合直方图；hab(indexx，indexy)=hab(indexx，indexy)+1；ha(indexx)=ha(indexx)+1；hb(indexy)=ha(indexy)+1；进而计算出两幅图像的信息熵和联合信息熵；P=hab/hsum；Hab=sum(sum(-p(index).*lop(p(index))))；P=ha/hsum；Ha=sum(sum(-p(index).*lop(p(index))))；P=hb/hsum；Hb=sum(sum(-p(index).*lop(p(index))))；最后利用公式mi=Ha+Hb—Hab计算出两幅图像的互信息值。实现该功能的函数文件主要为MMI.m。 该算法在计算互信息时，多次利用直方图求取互信息；ha(indexx)=ha(indexx)+1；hb(indexy)=ha(indexy)+1；直方图是多种数字图像处理技术的基础，能够提供十分有用的图像统计信息。（1）直方图的定义设数字图像的灰度级值，则的直方图可用离散函数表示为k=0,1,…,L-1（3-1）式中，表示第k级灰度值，表灰度图像中灰度级值为的像素个数。为了直观，直方图函数通常用图形表示。在直方图曲线中，水平轴对应度级，纵轴对应于。如图3.3所示：图3.3直方图（2）归一化直方图设图像的归一化直方图函数为，总像素为n，对式（3-1）两边同除以n，则得到的表达式k=0,1,…,L-1（3-2）由于n=,因此上式可以表示为：k=0,,1,…,L-1(3-3)从式3-3可知，就是灰度级的概率分布函数。显然，，即一个归一化的直方图的所有部分之和等于1。在获取互信息时，也需要计算出两幅图像的联合直方图；hab(indexx，indexy)=hab(indexx，indexy)+1；设和均为MN的数字图像，则图像A和B的联合直方图函数可表示为：(3-4)i=1,…,Mj=1,…,N式中，和表示在两幅图像相同位置(i,j)的一对灰度级值，表示同一灰度级值对出现的次数。可见联合直方图就是统计两幅图对应点的灰度级值对出现次数的函数。联合直方图数也可用二维图形表示：x轴表示图像A的灰度值，y轴表示图像B的灰度值。当图像A与图像B完全一致时，图形的所有点都分布在一条斜率为1的直线上。如图3.4所示：图3.4二维联合直方图与归一化直方图类似，对式（3-4）两边同除以全部的灰度级值对出现的次数和n，即可得到归一化的联合直方图函数：（3-5）i=1,2,…,Mj=1,…,N显然，归一化的联合直方图函数就是两幅图像对应点的灰度级值对的联合概率分布函数。从归一化直方图和归一化的联合直方图的定义中，不难看出图像A和图像B的概率分布函数、及它们的联合概率分布函数满足下列关系（3-6）（3-7）3.4灰度级插值技术本次算法中应用了双线性插值法，该方法效果较好，运算时间也比较快，是最邻近插值法的3至4倍。双线性插值方法假定内插点P四周四个点围成的区域内的灰度变化是线性的，从而可以用线性内插法，根据四个临近像素的灰度值，计算出内插点P的灰度值。实现步骤如下所示：dx=(x0-x)；dy=(y0-y)；pixel1=(1-dx)*im(x,y)+dx*im(x+1,y)；pixel2=(1-dx)*im(x,y+1)+dx*im(x+1,y+1)；pixel=(1-dy)*pixel1+dy*pixel2；实现该功能的函数文件为BilinearInterplate.m。3.5优化搜索算法在优化搜索中保持n个搜索方向线性无关至关重要。为了克服基本Powell算法中无法保证搜索方向的线性无关问题，出现了改进的Powell算法，改进的Powell算法优势在变量多时更加明显。改进的Powell算法实现如下：（1）给定允许误差初始点和n个线性无关的方向…，置k=1。（2）置，从出发，依次沿方向，，…，进行一维搜索，得到点，，…，求m,使得令，若则停止计算；否则，进行步骤（3）。（3）求使得令，若则停止计算，得；否则进行步骤（4）。（4）若则令转到步骤（2）利用该算法进行优化搜索，首先设定初始点和允许误差范围，在矩阵中存放D1、D2、D3三个方向矢量。然后从初始点出发，沿方向D1进行搜索，得到最大值fx1，对应点为x1，从x1出发，沿方向D2进行搜索，得到最大值fx2，对应点为x2，从x2出发，沿方向D3进行搜索，得到最大值fx3，对应点为x3。D4=x3-初始点；D4为这一轮产生的新方向矢量。若距离(x[k,n]-x[k,0])<=e(允许误差)则停止计算，得到结果x3。从初始点出发，沿D4方向进行一轮搜索，得到最大值fx4，对应点为x4。初始点换为x4。。同样，若距离(x[k,n]-x[k,0])<=e(允许误差)则停止计算，得到结果x3。判断是否跟换方向矢量。共轭性得到增强，替换选定的方向矢量。否则，共轭性没有的到增强不替换方向向量，从新进行搜索。同样，若距离(x[k,n]-x[k,0])<=e(允许误差)则停止计算，得到结果x3。在此算法中一维搜索采用Brent算法。最后找出最优结果。实现该功能的函数文件为powell.m。3.6一维搜索算法Brent法是一种综合了黄金分割法和抛物线法的搜索方法，巧妙地将两种算法的优点结合起来，具有收敛速度快和鲁棒性好的优点。具体思路如下：首先确定一个包含极小值点的不确定区域，该过程用黄金分割法，然后确定一个极小值点，该过程用抛物线法，若此极小值点落在此不确定区间，则利用该极小值继续进行二次插值；否则放弃该极小值点，改用黄金分割法进行搜索。Brent法在搜索过程中密切关注a、b、u、v、w、x六个点，其中a,b表示包含极小值点的不确定区间[a,b]。u表示最新搜索到的极小值点；w表示上一次搜索到的极小值点；v表示上一次搜索到的w值；x表示目前已经搜索到的最佳极小值点。设目标函数为法f(x),则Brent算法的实现过程如下：（1）给定出初始区间[a,b],精度>0，黄金分割系数cgold=0.381966,k=1。（2）计算v=a+cgold(b-a)，置x=v，w=v；计算置上一次迭代步长e=0.0。（3）计算当前区间中点，则停止搜索，得极小值x,否则转（4）（4）令则采用黄金分割法，转步骤（8）。（5）若（6）（7）若u在[a,b]之外，则用黄金分割法重新求极小变量，则转步骤（8）；若u-a<tol2或b-u<tol2,则用toll代替前面的改变量。（9）计算法f(u),按照各自的定义更新a,b,x,w,v,f(x),f(w),f(v)。置k=k+1，转步骤（3）。实现该功能的函数文件为oneDimSearch.m。测试和运行4.1实验环境本软件的运行环境需要一定的软件以及硬件设备，以便于实现较为理想的运行结果，其所需软件及硬件条件如下：硬件设备：微机一台操作系统：Windowsxpsp3编程语言：MATLAB语言运行环境：MATLAB开发技术：图像配准的相关算法4.2环境配置安装MATLAB7.0系统，需运行系统自带的安装程序setup.exe。运行后，将进入MATLAB安装程序的欢迎对话框(WelcometotheMathWorksInstaller)，用户可以按照安装提示依次操作。MATLAB系统的启动

启动MATLAB系统有3种常见方法：

(1)在Windows桌面，单击任务栏上的“开始”按钮，选择“程序”菜单项，然后选择“MATLAB7.0”程序选项。

(2)运行MATLAB系统启动程序matlab.exe。

(3)利用快捷方式。MATLAB系统的退出

要退出MATLAB系统，也有3种常见方法：

(1)在MATLAB主窗口File菜单中选择ExitMATLAB命令。

(2)在MATLAB命令窗口输入Exit或Quit命令。

(3)单击MATLAB主窗口的“关闭”按钮。4.3运行结果此图像为配准时的标准图像，为待配准图像的参考图像。如图4.1所示：图4.1CT原图像通过调用图像的转换函数，对原图像进行相应的坐标移动、角度的旋转等操作形成待配准图像。如图4.2所示：图4.2待配准图像两幅图像如果已经达到配准，此时互信息值最大，并且两幅图像在几何上达到最佳位置，图4.3为待配准图像经过平移和旋转变换达到互信息值最大时与原图像的配准情况。配准后的图像如图4.3所示：图4.3配准后的CT图像结论多模态图像配准在医学研究中具有极为重要的应用价值，随着科技的发展，医学分析和诊断多利用多模态图像，尤其医学图像配准是现代医学图像处理技术应用的一个重要方面，它是图像融合的第一步。不同传感器所具有的特点不同，为了将多源图像进行有效地组织，发挥各自的优点，我们必须先研究各种传感器所具有的特点，针对这些传感器不同的特点，进行优势互补，帮助我们所获多源图像进行解译。不同模态的图像成像原理不同，而且提供的信息针对性也不同。临床上根据图像提供的信息内容将医学图像分为两大类:一类图像用于描述生理形态的解剖结构的解剖图像，另一类图像用于描述生理代谢或功能信息的功能图像。前者包括计算机X射线断层成像(简称CT)、核磁共振成像(简称MR)、超声成像(简称US)和B型超声检查(简称B超)等。这类图像能准确反映人体器官的解剖结构信息并且有较高的分辨率，但是不能清晰地反映生理的功能信息。功能图像包括正电子发射断层成像(简称PET)、功能磁共振成像(简称FMR)和单光子发射计算机断层成像(简称SPECT)等。这类图像虽然能够提供脏器的功能代谢信息，但是它在反映人体器官的解剖结构信息时分辨率较差。参考文献[1]王菊花,吴增印.多光谱图像无失真压缩的方法研究[J].空间电子技术.2001(Z1)[2]安斌,陈书海,陈华,严卫东,孙煜.纹理特征在多光谱图像分类中的应用[J].激光与红外.2002(03)[3]雷杰,周有喜,吴成柯,李云松,孔繁锵.一种新的高效干涉多光谱图像压缩算法[J].光子学报.2009(06)[4]刘斌,彭嘉雄.基于二通道不可分加性小波的多光谱图像融合[J].光学学报.2007(08)[5]刘娟妮,周诠,李小军,方海.多光谱图像变换域统计量移位鲁棒无损信息隐藏[J].电讯技术.2015(11)[6