互信息的图像配准

1、信息论大作业基于互信息的图像配准 班级： 09030901 学号： 2009302311 姓名： 益琛 同组成员：陈升富 黎照1. 引言随着医学、计算机技术及生物工程技术的发展，医学影像学为临床诊断提供了多种模态的医学图像，不同的医学图像提供了相关脏器的不同信息：CT(Computed Tomography，电子计算机X 射线断层扫描)和MRI(Magneticresona nce ima ging，核磁共振成像)以较高的空间分辨率提供了脏器的解剖结构信息。在实际临床应用中，单一模态的图像往往不能提供医生所需要的足够的信息，通常需要将不同模态的图像融合在一起，得到更丰富的信息，以便了解病变组织

2、或器官的综合信息，从而做出准确的诊断或制订出合适的治疗方案。而图像配准是图像融合的重要前提，图像配准是指对一幅图像进行一定的几何变换而映射到另一幅图像中，使得两幅图像中的相关点达到空间上的一致。图像配准主要有两大类方法，基于灰度的方法和基于特征的方法。基于灰度的配准方法直接利用图像的灰度数据进行配准，从而避免了因分割而带来的误差，因而具有精度较高、鲁棒性强、不需要预处理而能实现自动配准的特点。在基于灰度的配准方法中，基于互信息的方法包括互信息和归一化互信息方法，它们已经被广泛使用并具有最高的精度。本文使用的是基于互信息的配准方法。2. 图像配准技术2.1图像配准技术的数学定义 数字图像可以用一

3、个二维矩阵来表示，如果用、分别表示待配准图像和参考图像在点(x,y)处的灰度值，那么图像、的配准关系可表示为： （1）其中f代表二维的空间几何变换函数；g表示一维的灰度变换函数。配准的主要任务是寻找最佳的空间变换关系f与灰度变换关系g，使两幅图像实现最佳对准。其中，空间几何变换是灰度变换的前提，是实现精准配准的关键环节。2.2几何变换空间变换主要解决图像平面上像素的重新定位问题，式(1)中的空间几何变换函数f可用空间变换模型进行描述，常用的空间变换模型有刚体变换、仿射变换、投影变换和非线性变换。刚体变换使得一幅图像中任意两点间的距离变换到另一幅图像中后仍然保持不变；仿射变换使得一幅图像中的直线

4、经过变换后仍保持直线，并且平行线仍保持平行；投影变换是从三维图像到二维平面的投影；非线性变换把一条直线变换为一条曲线，一般用代数多项式来表示。仿射变换是最常用的一种空间变换形式,可以实现图像的平移、旋转、按比例缩放等操作，我们在实验中使用的是此变换模型。仿射变换可以用矩阵形式表示： 当分别取值为、将依次对图像进行平移、旋转、按比例缩放操作。2.3 插值技术浮动图的像素点经过空间变换后，参考图中对应点的坐标一般来说不是整数，必须通过插值方法计算该点的灰度值。常用的插值算法有最近邻插值算法、双线性插值算法和部分体积插值算法。为了尽可能避免基于互信息配准的局部最优问题，本文采用改进PV插值算法。PV

5、插值法是一种专门针对两幅图像的联合直方图的更新而设计的插值技术，它并不是真正意义上的插值方法，因为通过此方法并不能计算出反向变换点的灰度值。PV插值法的计算过程如图1.图中的(x)为反向变换得到的一个浮点数点，其四个最近邻像素点分别为。设参考图像为r(x),浮动图像为f(x),则它们的联合图方图函数如下。 i=1,2,3,41-dxdxdy1-dy2.4 优化算法常用的优化算法有：牛顿法、最速下降法、模拟退火法、遗传算法、单纯形法、模式搜索法、Powell法等搜索算法。Powell法不需要对目标函数进行求导计算，具有收敛速度快、精度高、可靠性好等优点，是目前解无约束最优化问题十分有效的直接法，

6、应用相当广泛，所以我们在实验中采用该算法。Powell算法实现如下：(1) 给定允许误差，初始点和n个线性无关的方向 ，置k=1.(2) 置，从出发，依次沿方向进行一维搜索，得到点。再从出发，沿方向作一维搜索点，得到点。（3）若，则停止搜索，得到点；否则，置 返回步骤（2）.2.4.2 Powell算法中的一维搜索算法brent方法。Brent法思路：开始时利用黄金分割法确定一个较小的包含极小点的不确定区间，然后利用抛物线法获得一个极小点，若此极小点落在此不确定区间，则利用该极小点继续进行二次插值；否则放弃该点，改用黄金分割法搜索。算法中密切关注a,b,u,v,w,x这六个点，其中a,b表示包

7、含极小点的不确定区间；u表示最新搜索到的极小点；w表示上一次搜索到的极小点；v表示上一次的w值；x表示当前已搜到的最佳极小点。算法实现步骤如下（设目标函数为f(x)）：(1)给定初始区间，精度要求，黄金分割系数（2）计算，置；计算，置；置上一次迭代步长。（3）计算当前区间中点，若，则停止搜索，的极小值,否则转（4）。（4）令，若，则采用黄金分割法，转（8）。（5）若，则采用黄金分割法，转（8）。（6）过三点构造抛物线函数，计算（7）若在之外，则用黄金分割法重新求极小点，转步骤（8）；若u相对于的改变量大于上一次的改变量，则转步骤（8）；若或，则用代替前面的改变量。（8）按黄金分割法确定点,且u

8、在区间和中长度较大的一个；若u相对于x的改变量小于，则用代替前面的改变量。（9）计算，按照各自的定义更新。置，转步骤（3）。3 基于互信息的图像配准方法3.1 互信息的计算互信息是信息理论中的一个基本概念，通常用于描述两个系统间的统计相关性，或者是一个系统中所包含的另一个系统中信息的多少，它可以用熵来描述： （2） 其中，和分别是系统A和B的熵，是它们的联合熵，依次定义如下： (3) (4) (5)其中和分别是系统A和B完全独立时的的概率分布。是系统A和B的联合概率分布。令图像A和B的互信息为，将式（3），（4），（5），分别代入式（2），即可得到图像互信息的计算公式： 3.2 配准方法首先根

9、据两幅图像的基本情况预设一个初始参数，其中为裁剪旋转角的图像2 行的第一个索引。为裁剪旋转角的图像2 列的第一个索引，为旋转的角度，为比例因子。然后按照给定的初始参数对图像2 进行变换，并计算图像1 和图像2 的互信息，然后利用最优化工具箱中的fminsearch 函数在附近寻找使图像1 和图像2 互信息最大的点，直至搜索到满足精度要求的参数；最后输出配准参数。4. 图像配准的实现4.1配准流程首先对参考图像和浮动图像按照给定的初始点使用PV插值法统计联合直方图并计算互信息值；然后利用POWELL算法依据最大互信息理论判断所得参数是否最优，若不是，则继续搜索较优参数，在搜索时会不断重复“空间几

10、何变换（affine）-统计联合直方图（PV插值法）-计算互信息值-最优化判断”的过程，直至搜索到满足精度要求的参数；最后输出配准参数。输入参考图像输入浮动图像设置初始点和初始搜素方向空间几何变换计算互信息值最优化否是输出配准参数4.2.所用到的M文件及其源代码4.2.1 ImageRegistration.mfunction varargout = ImageRegistration(varargin)gui_Singleton = 1;gui_State = struct('gui_Name', mfilename, . 'gui_Singleton', g

11、ui_Singleton, . 'gui_OpeningFcn', ImageRegistration_OpeningFcn, . 'gui_OutputFcn', ImageRegistration_OutputFcn, . 'gui_LayoutFcn', , . 'gui_Callback', );if nargin && ischar(varargin1) gui_State.gui_Callback = str2func(varargin1);end if nargout varargout1:nargo

12、ut = gui_mainfcn(gui_State, varargin:);else gui_mainfcn(gui_State, varargin:);end addpath(pwd);function ImageRegistration_OpeningFcn(hObject, eventdata, handles, varargin)handles.output = hObject; guidata(hObject, handles); function varargout = ImageRegistration_OutputFcn(hObject, eventdata, handles

13、)varargout1 = handles.output; function pushbutton1_Callback(hObject, eventdata, handles)global I; %调用OpenImage.m读入参考图像并获取文件名、图像大小%filename ,pathname=uigetfile('*.jpg''*.bmp''*.bmp','Ñ¡ÔñͼƬ');str=pathname filename;I=imread(str)

14、;axes(handles.axes1);imshow(I);handles.data=I;guidata(hObject,handles);figure(1);imshow(handles.data); function pushbutton3_Callback(hObject, eventdata, handles)handles.Old_I=handles.data;handles.Old_J=handles.data2; I,J=GLPF(handles);handles.data=I;handles.data2=J;guidata(hObject,handles); ticRegis

15、trationParameters=Powell(handles);tocElapsedTime=toc; handles.RegistrationParameters=RegistrationParameters;y=RegistrationParameters(1);x=RegistrationParameters(2);ang=RegistrationParameters(3);MI_Value=RegistrationParameters(4);RegistrationResult=sprintf('X,Y,Angle=%.5f%.5f%.5f',x,y,ang);MI

16、_Value=sprintf('MI_Value=%.4f',MI_Value);ElapsedTime=sprintf('Elapsed Time=%.3f',ElapsedTime);axes(handles.axes3)RegistrationImage=Register(handles);imshow(RegistrationImage)set(handles.text1,'string',RegistrationResult);set(handles.text2,'string',MI_Value);set(handle

17、s.text3,'string',ElapsedTime); function pushbutton2_Callback(hObject, eventdata, handles)global J;filename ,pathname=uigetfile('*.jpg''*.bmp''*.bmp','Ñ¡ÔñͼƬ');str=pathname filename;J=imread(str);axes(handles.axes2);ims

18、how(J);handles.data2=J;guidata(hObject,handles);figure(2);imshow(handles.data2);4.2.2 Powell.mfunction RegistrationParameters=Powell(handles)e=0.1;X0=0 0 0;D=1 0 0;0 1 0;0 0 1;num=0;while(num<200) num=num+1; d1=D(1,:);%d1Ϊ¾ØÕóDµÄµÚÒ»

19、08;У¬³õʼËÑË÷·½Ïò X1,fX1=OneDimSearch(X0,d1,handles); d2=D(2,:);%d2Ϊ¾ØÕóDµÄµÚ¶þÐУ¬³õʼËÑË÷·½

20、Ïò X2,fX2=OneDimSearch(X1,d2,handles); d3=D(3,:);%d3Ϊ¾ØÕóDµÄµÚÈýÐУ¬³õʼËÑË÷·½Ïò£¬ÈýάËÑË÷¹

21、02;ÓÐÈý¸ö·½Ïò X3,fX3=OneDimSearch(X2,d3,handles); fX0=PV(X0(1),X0(2),-X0(3),handles); Diff=fX1-fX0 fX2-fX1 fX3-fX2; maxDiff,m=max(Diff);%maxº¯ÊýµÄÓ÷¨£¬·µ»ØmaxdiffÎ&

22、#170;ÏòÁ¿DiffµÄ×î´óÔªËØ£¬mΪ×î´óÔªËصÄÐòºÅ d4=X3-X0; temp1=X3-X0; Conditon1=sum(temp1.*temp1); if Conditon1<=e break end X4,fX4,landa=One

23、DimSearch(X0,d4,handles); X0=X4; temp2=X4-X3; Conditon2=sum(temp2.*temp2); if Conditon2<=e X3=X4; break end temp3=sqrt(fX4-fX0)/(maxDiff+eps); if(abs(landa)>temp3) D(4,:)=d4; for i=m:3 D(i,:)=D(i+1,:); end endendRegistrationParameters(1)=-X3(1);RegistrationParameters(2)=-X3(2);RegistrationPara

24、meters(3)=-X3(3);RegistrationParameters(4)=fX3;function Y,fY,landa=OneDimSearch(X,direction,handles)%һάËÑË÷²ÉÓÃbrent·½·¨a=-5;b=5;Epsilon=0.0001;cgold=0.381966;IterTimes=200;if a>b temp=a; a=b; b=temp;endv=a+cgold*(b

25、-a);w=v;x=v;e=0.0;fx=Fx(x,X,direction);fv=fx;fw=fx;for iter=1:IterTimes xm=0.5*(a+b); if abs(x-xm)<=Epsilon*2-0.5*(b-a) break end if abs(e)>Epsilon r=(x-w)*(fx-fv); q=(x-v)*(fx-fw); p=(x-v)*(q-(x-w)*r; q=2*(q-r); if q>0 p=-p; end q=abs(q); etemp=e; e=d; if not (abs(p)>=abs(0.5*q*etemp)|p

26、<=q*(a-x)|p>=q*(b-x) d=p/q; u=x+d; if u-a<Epsilon*2|b-u<Epsilon*2 d=MySign(Epsilon,xm-x); end else if x>=xm e=a-x; else e=b-x; end d=cgold*e; end else if x>=xm; e=a-x; else e=b-x; end d=cgold*e; end if abs(d)>=Epsilon u=x+d; else u=x+MySign(Epsilon,d); end fu=Fx(u,X,direction);

27、if fu<=fx if u>=x a=x; else b=x; end v=w; fv=fw; w=x; fw=fx; x=u; fx=fu; else if u<x a=u; else b=u; end if fu<=fw|w=x v=w; fv=fw; w=u; fw=fu; else if fu<=fv|v=x|v=w v=u; fv=fu; end end endendlanda=x;Y=X+x*direction;fY=Fx(x,X,direction);functionmySign=MySign(a,b)if b>0 Result=abs(a)

28、;else Resulr=abs(a);endmySign=Result;function fx=Fx(x,X,direction,handles)fx=-PV(X(1)+direction(1)*x,X(2)+direction(2)*x,-(X(3)+direction(3)*x),handles); 4.2.3 PV.mfunction mi=PV(x,y,ang,handles)a=handles.data;a=double(a);b=handles.data2;b=double(b);M,N=size(a);hab=zeros(256,256);ha=zeros(1,256);hb=

29、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(max(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

30、 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;tform=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

31、 j=1:height coordinate_y(i,j)=j; endendw z=tforminv(tform,coordinate_x,coordinate_y);for i=1:width for j=1:height source_x=w(i,j); source_y=z(i,j); if(source_x>width-1|source_y>height-1|double(uint16(source_x)<=1|double(uint16(source_y)<=1) hab(a(1,1),a(1,1)=hab(a(1,1),a(1,1)+1; else m=f

32、ix(source_x); n=fix(source_y); index_b=b(i,j); index_a0=a(m,n); index_a1=a(m+1,n); index_a2=a(m,n+1); index_a3=a(m+1,n+1); dx=source_x-m; dy=source_y-n; hab(index_a0,index_b)=hab(index_a0,index_b)+(1-dx)*(1-dy); hab(index_a1,index_b)=hab(index_a1,index_b)+dx*(1-dy); hab(index_a2,index_b)=hab(index_a

33、2,index_b)+(1-dx)*dy; hab(index_a3,index_b)=hab(index_a3,index_b)+dx*dy; end endendhabsum=sum(sum(hab);index=find(hab=0);pab=hab/habsum;Hab=sum(sum(-pab(index).*log2(pab(index);pa=sum(pab,2);index=find(pa=0);Ha=sum(sum(-pa(index).*log2(pa(index);pb=sum(pab,1);index=find(pb=0);Hb=sum(sum(-pb(index).*

34、log2(pb(index);mi=Ha+Hb-Hab;4.2.4 Register.mfunctionRegistrationImage=Register(handles)I=handles.data;J=handles.data2;y=handles.RegistrationParameters(1);x=handles.RegistrationParameters(2);ang=-handles.RegistrationParameters(3); nrows,ncols=size(J);width=nrows;height=ncols;new_J=uint8(zeros(width,h

35、eight); a=(width-1)/2;c=a;b=(height-1)/2;d=b; rad=pi/180*ang;t1=1 0 0;0 1 0;x y 1;t2=1 0 0;0 1 0;-a -b 1;t3=cos(rad) -sin(rad) 0;sin(rad) cos(rad) 0;0 0 1;t4=1 0 0;0 1 0;c d 1;T=t2*t3*t4*t1;tform=maketform('affine',T);tx=zeros(width,height);ty=zeros(width,height);for i=1:width for j=1:height tx(i,j)=i; endendfor i=1:width for j=1:height ty(i,j)=j; endendw z=tforminv(tform,tx,ty);for i=1:width for j=1:height source_x=w(i,j); source_y=z(i,j); if(source_x>width-1|source_y>height-1|double(uint16(source_x)<=0|double(uint16(source_y)<=0) new_J(i,j)=J(1,1); else