基于互信息的图像配准

1、 下面代码MATLAB编程实现了基于互信息的图像配准。采取的是改进的Powell搜索策略算法，在4个方向上采用的是黄金分割法。当基准图像和待配准图像的互信息达到最大时，两幅图像配准的效果最好，同时获得此时图像的变换系数，然后进行空间变换，获得图像融合图。运行结果一下图中图4-1为十字架参考图像，图4-2为十字架浮动图像，图4-3为获得空间变换参数后的变换图像，图4-4为参考图像和浮动图像的配准后图像 图4-1参考图像 图4-2 浮动图像 图4-3变换图像 图4-4 配准效果图 运行结果二 在下图中， 图4-5 为一幅名为圣保罗教堂的钥匙图像，是参考图像，图4-6为浮动图像，该图像与图4-5相比

2、，X轴和Y轴方向上发生一定的移动，并且发生一定的旋转。图4-7为获得空间变换参数后经过空间变换的图像，图4-8为配准效果图 图4-5 参考图像 图 4-6 浮动图像 图 4-6 变换图像 图4-7 配准效果图 运行结果三下图中图4-9为一幅人脸图像，作为参考图。图4-10为另一幅人脸图像，二者的神情是有一定差异的，作为浮动图像，将这两张图像进行配准。图4-11是进行优化策略获得配准参数，对图4-10经过空间变换获得的图像。图4-12为参考图像和浮动图像的配准效果图。 图4-9人脸参考图像 图4-10 人脸浮动图像 图4-11 人脸变换图像 图4-12 配准效果图 表4-1为程序运行时，获得的图

3、像匹配参数以及互信息值和图像配准运行时间。 表4-1 配准参数表 图像标号X轴偏移量 Y轴偏移量 旋转角度 互信息值图像配准时间(s)4-1与4-20.46625-0.507017.190950.6760131.4-5与4-60.529442.842413.550351.8618301.4-8与4-94.06384-14.094147.596911.3424686. 从互信息大小来看，图4-5和图4-6的互信息值最大，得到的图像配准的效果最好，但是我们不能把互信息的大小作为评判三队图像中那对图像配准的好坏，因为带配准的图像不论图和平移旋转，两幅图像的相互关系是一定的即互信息的大小是一定的，它们

4、所包含的相同信息是有限度的，所以，只能与本身图像或者参考图像相比，来判断图像配准的好坏。从配准时间来看，可以发现，图像越复杂，图像配准花费的时间越长。比如对于十字架图像的配准来说，图像比较简单，花费的时间最少；对于人脸图像来说，图像相对来说比较复杂，花费的时间最长。从图像配准的整体时间来看，基于Harris角点的图像配准花费的时间才几秒，而基于互信息的图像处理花费的时间远远大于几秒，这是基于互信息图像配准的缺点。具体上对本方法分析，一方面是互信息计算涉及到大量的浮点运算，另一方面是采用的是Powell搜索策略，在搜索的每个方向上，用黄金分割法搜索每一小步，都需要计算并比较互信息，花费很长时间。

5、从图像配准的效果来看，尤其是对于刚体变换，获得到的空间变换参数是很精确的图像的配准效果还是不错的。本文使用的方法获得参数是3个，对于二维变换，还能满足要求，但是对于三维图像，匹配的效果很差，不能满足要求。其代码如下：function Q=Peizhun(I,J)I=imread(cross1.jpg);J=imread(cross2.jpg);if isrgb(I) I1 = rgb2gray(I);else I1=I;endif isrgb(J) J1= rgb2gray(J);else J1=J;endI=double(I1);J=double(J1);ticRegistrationPar

6、ameters=Powell(I,J);tocElapsedTime=toc;x=RegistrationParameters(1);y=RegistrationParameters(2);ang=RegistrationParameters(3);MI_Value=RegistrationParameters(4);RegistrationResult=sprintf(X,Y,Angle=%.5f%.5f%.5f,x,y,ang)MI_Value=sprintf(MI_Value=%.4f,MI_Value)X=round(x);Y=round(y);Ang=round(ang);se =

7、translate(strel(1),X Y);J2 = imdilate(J1,se);J3 = imrotate(J2,Ang,nearest,crop);%FusionImage,RegistrationImage=Fusion(x,y,ang,I,J);%imshow(FusionImage)FusionImage=Fus(I1,J3);figure(1);imshow(I1);title(Base image);figure(2);imshow(J1);title(Float image);figure(3);imshow(J3);title(Trans image);figure(

8、4);imshow(FusionImage);title(Registration image);function FusionImage=Fus(I1,J3)%获取原图像及配准参数%对浮动图像空间变换及插值，以便得到配准后的输出图像nrows,ncols=size(J3);width=nrows;height=ncols;%new_J=uint8(zeros(width,height); %产生“融合”后的效果图 I=double(I1)/255; J=double(J3)/255; IJ=double(zeros(width,height); for m=1:width for n=1:h

9、eight if I(m,n)0.999 | J(m,n)0.999 IJ(m,n)=0.8; elseif I(m,n)=0 |J(m,n)=0 IJ(m,n)=0.01; else IJ(m,n)=(I(m,n)*0.3+J(m,n)*0.7); end end end IJ=IJ*255; IJ=uint8(IJ); FusionImage=IJ;function fx=Fx(x,X,direction,I,J)fx=-PV(X(1)+direction(1)*x,X(2)+direction(2)*x,-(X(3)+direction(3)*x),I,J);function mi=PV

10、(x,y,ang,I,J)% x 水平平移像素 y -垂直平移像素 ang -旋转角度 利用PV插值法算出两幅图像的联合直方图的值，结果返回联合直方图数组的值%a=double(I);b=double(J);M,N=size(a);hab=zeros(256,256);ha=zeros(1,256);hb=zeros(1,256);if max(max(a)=min(min(a) a=(a-min(min(a)/(max(max(a)-min(min(a);else a=zeros(M,N);endif max(max(b)-min(min(b) b=(b-min(min(b)/(max(ma

11、x(b)-min(min(b);else b=zeros(M,N);enda=double(int16(a*255)+1;b=double(int16(b*255)+1;%反射变换空间几何变换%width,height=size(b);u=(width-1)/2;v=(height-1)/2;rad=pi/180*ang;t1=1 0 0;0 1 0;x y 1;t2=1 0 0;0 1 0;-u -v 1;t3=cos(rad) -sin(rad) 0;sin(rad) cos(rad) 0;0 0 1;t4=1 0 0;0 1 0;u v 1;T=t2*t3*t4*t1; %几何变换tfo

12、rm=maketform(affine,T);coordinate_x=zeros(width,height);coordinate_y=zeros(width,height);for i=1:width for j=1:height coordinate_x(i,j)=i; endendfor i=1:width for j=1:height coordinate_y(i,j)=j; endendw z=tforminv(tform,coordinate_x,coordinate_y); %反向空间变换%利用PV插值法统计联合直方图%for i=1:width for j=1:height

13、source_x=w(i,j); source_y=z(i,j); if (source_xwidth-1 | source_yheight-1 | double(uint16(source_x)=1 | double(uint16(source_y)=e) if FcFd Fc=Fd; a=c; b=b; c=d; d=a+0.618*(b-a); Fd=Fx(d,X,direction,I,J); else Fd=Fc; a=a; b=d; d=c; c=a+0.382*(b-a); Fc=Fx(c,X,direction,I,J); endendY=X+(b+a)/2)*directio

14、n;landa=(b+a)/2;fY=-Fx(landa,X,direction,I,J);function RegistrationParameters=Powell(I,J)e=0.1;X0=-5 -3 -4;D=1 0 0;0 1 0;0 0 1;while(1) d1=D(1,:); %提取矩阵中第一行的元素 X1,fX1,landa=OneDimSearch(X0,d1,I,J); d2=D(2,:); %提取矩阵中第二行的元素 X2,fX2,landa=OneDimSearch(X1,d2,I,J); d3=D(3,:); %提取矩阵中第三行的元素 X3,fX3,landa=One