基于互信息的图像配准

信息论大作业基于互信息的图像配准班级：金融101学号：2009302311姓名：魏泉1. 引言随着医学、计算机技术及生物工程技术的发展，医学影像学为临床诊断提供了多种模态的医学图像，不同的医学图像提供了相关脏器的不同信息：CT(Computed Tomography ，电子计算机X 射线断层扫描)和MRI(Magneticresona nce ima ging ，核磁共振成像)以较高的空间分辨率提供了脏器的解剖结构信息。

在实际临床应用中，单一模态的图像往往不能提供医生所需要的足够的信息，通常需要将不同模态的图像融合在一起，得到更丰富的信息，以便了解病变组织或器官的综合信息，从而做出准确的诊断或制订出合适的治疗方案。

而图像配准是图像融合的重要前提，图像配准是指对一幅图像进行一定的几何变换而映射到另一幅图像中，使得两幅图像中的相关点达到空间上的一致。

图像配准主要有两大类方法，基于灰度的方法和基于特征的方法。

基于灰度的配准方法直接利用图像的灰度数据进行配准，从而避免了因分割而带来的误差，因而具有精度较高、鲁棒性强、不需要预处理而能实现自动配准的特点。

在基于灰度的配准方法中，基于互信息的方法包括互信息和归一化互信息方法，它们已经被广泛使用并具有最高的精度。

本文使用的是基于互信息的配准方法。

2. 图像配准技术2.1图像配准技术的数学定义 数字图像可以用一个二维矩阵来表示，如果用),(1y x I、),(2y x I 分别表示待配准图像和参考图像在点(x,y)处的灰度值，那么图像I 1、I 2的配准关系可表示为：))),(((),(12y x f g y x I I= （1）其中f 代表二维的空间几何变换函数；g 表示一维的灰度变换函数。

配准的主要任务是寻找最佳的空间变换关系f 与灰度变换关系g ，使两幅图像实现最佳对准。

其中，空间几何变换是灰度变换的前提，是实现精准配准的关键环节。

2.2几何变换空间变换主要解决图像平面上像素的重新定位问题，式(1)中的空间几何变换函数f 可用空间变换模型进行描述，常用的空间变换模型有刚体变换、仿射变换、投影变换和非线性变换。

刚体变换使得一幅图像中任意两点间的距离变换到另一幅图像中后仍然保持不变；仿射变换使得一幅图像中的直线经过变换后仍保持直线，并且平行线仍保持平行；投影变换是从三维图像到二维平面的投影；非线性变换把一条直线变换为一条曲线，一般用代数多项式来表示。

仿射变换是最常用的一种空间变换形式,可以实现图像的平移、旋转、按比例缩放等操作，我们在实验中使用的是此变换模型。

仿射变换可以用矩阵形式表示：1[x 1y 1]=0[x 0y 1]T =0[x 0y 1]111221223132001t t t t t t ⎛⎫⎪ ⎪ ⎪⎝⎭当T 分别取值为1000101xyt t ⎛⎫⎪ ⎪ ⎪⎝⎭、cos sin 0sin cos 0001θθθθ-⎛⎫ ⎪ ⎪ ⎪⎝⎭、000001xy s s ⎛⎫⎪⎪ ⎪⎝⎭将依次对图像进行平移、旋转、按比例缩放操作。

2.3 插值技术浮动图的像素点经过空间变换后，参考图中对应点的坐标一般来说不是整数，必须通过插值方法计算该点的灰度值。

常用的插值算法有最近邻插值算法、双线性插值算法和部分体积插值算法。

为了尽可能避免基于互信息配准的局部最优问题，本文采用改进PV 插值算法。

PV 插值法是一种专门针对两幅图像的联合直方图的更新而设计的插值技术，它并不是真正意义上的插值方法，因为通过此方法并不能计算出反向变换点的灰度值。

PV 插值法的计算过程如图1.图中的T α(x)为反向变换得到的一个浮点数点，其四个最近邻像素点分别为n n n n 4321,,,。

设参考图像为r(x),浮动图像为f(x),则它们的联合图方图函数h rf 如下。

1=∑iiw i=1,2,3,4w n h i i rf f x r i =+∀))(),(()1(*)1(1dy dx w --=)1(*2dy dx w -=dy dx w *)1(4-=dy dx w *3=2.4.1常用的优化算法有：牛顿法、最速下降法、模拟退火法、遗传算法、单纯形法、模式搜索法、Powell 法等搜索算法。

Powell 法不需要对目标函数进行求导计算，具有收敛速度快、精度高、可靠性好等优点，是目前解无约束最优化问题十分有效的直接法，应用相当广泛，所以我们在实验中采用该算法。

Powell 算法实现如下：(1)给定允许误差)0(<εε，初始点x )0(和n 个线性无关的方向dddn ),1()2,1()1,1(,...,, ，置k=1.(2)置x xk k )1()0,(-=，，从x k )0,(出发，依次沿方向dddn k k k ),()2,()1,(,...,,进行一维搜索，得到点xxx n k k k ),()2,()1,(,...,,。

再从x n k ),(出发，沿方向xxdk n k n k )0,(),()1,(-=+作一维搜索点，得到点x k )(。

（3）若ε<--x xk k )1()(，则停止搜索，得到点x k )(；否则，置n j ddj k j k ,...,1,)1,(),1(==++ 1+=k k返回步骤（2）.2.4.2 Powell 算法中的一维搜索算法——brent 方法。

Brent 法思路：开始时利用黄金分割法确定一个较小的包含极小点的不确定区间，然后利用抛物线法获得一个极小点，若此极小点落在此不确定区间，则利用该极小点继续进行二次插值；否则放弃该点，改用黄金分割法搜索。

算法中密切关注a,b,u,v,w,x 这六个点，其中a,b 表示包含极小点的不确定区间][b a ,；u 表示最新搜索到的极小点；w 表示上一次搜索到的极小点；v 表示上一次的w 值；x 表示当前已搜到的最佳极小点。

算法实现步骤如下（设目标函数为f(x)）：(1)给定初始区间][b a ,，精度要求0>ε，黄金分割系数1,381966.0==k cgold （2）计算)(*a b cgold a v -+=，置v w v x ==,；计算)(v f ，置)()(),()(v f x f v f w f ==；置上一次迭代步长0.0=e 。

（3）计算当前区间中点2/)(b a xm +=，若ε<-xm x ，则停止搜索，的极小值x ,否则转（4）。

（4）令1*22,*1tol tol x tol ==ε，若1tol e <，则采用黄金分割法，转（8）。

（5）若)()()(x f v f w f ====，则采用黄金分割法，转（8）。

（6）过三点))(,()),(,()),(,(v f v w f w x f x 构造抛物线函数，计算))()()(())()()(())()(()())()(()(*2/122v f x f w x w f x f v x v f x f w x w f x f v x x u -----------= （7）若u 在][b a ,之外，则用黄金分割法重新求极小点，转步骤（8）；若u 相对于x 的改变量大于上一次的改变量，则转步骤（8）；若2t o l a u <-或2tol u b <-，则用1tol 代替前面的改变量。

（8）按黄金分割法确定点u ,且u 在区间][x a ,和][b x ,中长度较大的一个；若u 相对于x 的改变量小于1tol ，则用1tol 代替前面的改变量。

（9）计算)(u f ，按照各自的定义更新)(),(),(,,,,,v f w f x f v w x b a 。

置1+=k k ，转步骤（3）。

3 基于互信息的图像配准方法3.1 互信息的计算互信息是信息理论中的一个基本概念，通常用于描述两个系统间的统计相关性，或者是一个系统中所包含的另一个系统中信息的多少，它可以用熵来描述：),()()(),(B A H B H A H B A I -+= （2） 其中，)(A H 和)(B H 分别是系统A 和B 的熵，),(B A H 是它们的联合熵，依次定义如下：)(log )()(a a A H P p A aA∑-= (3))(log )()(b b B H P p B bB∑-= (4)),(log ),(),(,b a b a B A H P p AB ba AB∑-= (5)其中)(,,a B b A a pA∈∈和)(b P B 分别是系统A 和B 完全独立时的的概率分布。

),(b a PAB是系统A 和B 的联合概率分布。

令图像A 和B 的互信息为),(B A I ，将式（3），（4），（5），分别代入式（2），即可得到图像互信息的计算公式：),(log ),()(log )()(log )(),(2,22b a b a a a a a B A I P P P P P P AB ba AB B bB A aA ∑∑∑+--=3.2 配准方法首先根据两幅图像的基本情况预设一个初始参数X 0，其中)1(0X 为裁剪 旋转)3(0X 角的图像2 行的第一个索引。

)2(0X 为裁剪旋转)3(0X 角的图像2X为旋转的角度，)4(0X为比例因子。

然后按照给定的列的第一个索引，)3(初始参数对图像2 进行变换，并计算图像1 和图像2 的互信息，然后利用最优化工具箱中的fminsearch 函数在X0附近寻找使图像1 和图像2 互信息最大的点，直至搜索到满足精度要求的参数；最后输出配准参数。

4. 图像配准的实现4.1配准流程首先对参考图像和浮动图像按照给定的初始点使用PV插值法统计联合直方图并计算互信息值；然后利用POWELL算法依据最大互信息理论判断所得参数是否最优，若不是，则继续搜索较优参数，在搜索时会不断重复“空间几何变换（affine）-统计联合直方图（PV插值法）-计算互信息值-最优化判断”的过程，直至搜索到满足精度要求的参数；最后输出配准参数。

4.2.所用到的M文件及其源代码4.2.1． ImageRegistration.mfunction varargout = ImageRegistration(varargin)gui_Singleton = 1;gui_State = struct('gui_Name', mfilename, ...'gui_Singleton', gui_Singleton, ...'gui_OpeningFcn', @ImageRegistration_OpeningFcn, ...'gui_OutputFcn', @ImageRegistration_OutputFcn, ...'gui_LayoutFcn', [], ...'gui_Callback', []);if nargin && ischar(varargin{1})gui_State.gui_Callback = str2func(varargin{1});endif nargout[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});elsegui_mainfcn(gui_State, varargin{:});endaddpath(pwd);function ImageRegistration_OpeningFcn(hObject, eventdata, handles, varargin)handles.output = hObject;guidata(hObject, handles);function varargout = ImageRegistration_OutputFcn(hObject, eventdata, handles)varargout{1} = handles.output;function pushbutton1_Callback(hObject, eventdata, handles)global I;%%%调用OpenImage.m读入参考图像并获取文件名、图像大小%%%[filename ,pathname]=uigetfile({'*.jpg';'*.bmp';'*.bmp'},'Ñ¡ÔñͼƬ'); str=[pathname filename];I=imread(str);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);ticRegistrationParameters=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_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(handles.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);imshow(J);handles.data2=J;guidata(hObject,handles);figure(2);imshow(handles.data2);4.2.2Powell.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µÄµÚÒ»ÐУ¬³õʼËÑË÷·½Ïò[X1,fX1]=OneDimSearch(X0,d1,handles);d2=D(2,:);%d2Ϊ¾ØÕóDµÄµÚ¶þÐУ¬³õʼËÑË÷·½Ïò[X2,fX2]=OneDimSearch(X1,d2,handles);d3=D(3,:);%d3Ϊ¾ØÕóDµÄµÚÈýÐУ¬³õʼËÑË÷·½Ïò£¬ÈýάËÑË÷¹ÊÓÐÈý¸ö·½Ïò [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ΪÏòÁ¿DiffµÄ×î´óÔªËØ£¬mΪ×î´óÔªËصÄÐòºÅd4=X3-X0;temp1=X3-X0;Conditon1=sum(temp1.*temp1);if Conditon1<=ebreakend[X4,fX4,landa]=OneDimSearch(X0,d4,handles);X0=X4;temp2=X4-X3;Conditon2=sum(temp2.*temp2);if Conditon2<=eX3=X4;breakendtemp3=sqrt((fX4-fX0)/(maxDiff+eps));if(abs(landa)>temp3)D(4,:)=d4;for i=m:3D(i,:)=D(i+1,:);endendendRegistrationParameters(1)=-X3(1);RegistrationParameters(2)=-X3(2);RegistrationParameters(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>btemp=a;a=b;b=temp;endv=a+cgold*(b-a);w=v;x=v;e=0.0;fx=Fx(x,X,direction);fv=fx;fw=fx;for iter=1:IterTimesxm=0.5*(a+b);if abs(x-xm)<=Epsilon*2-0.5*(b-a)breakendif abs(e)>Epsilonr=(x-w)*(fx-fv);q=(x-v)*(fx-fw);p=(x-v)*(q-(x-w)*r;q=2*(q-r);if q>0p=-p;endq=abs(q);etemp=e;e=d;if not(abs(p)>=abs(0.5*q*etemp)||p<=q*(a-x)||p>=q*(b-x)) d=p/q;u=x+d;if u-a<Epsilon*2||b-u<Epsilon*2d=MySign(Epsilon,xm-x);endelseif x>=xme=a-x;elsee=b-x;endd=cgold*e;endelseif x>=xm;e=a-x;elsee=b-x;endd=cgold*e;endif abs(d)>=Epsilonu=x+d;elseu=x+MySign(Epsilon,d);endfu=Fx(u,X,direction);if fu<=fxif u>=xa=x;elseb=x;endv=w;fv=fw;w=x;fw=fx;x=u;fx=fu;elseif u<xa=u;elseb=u;endif fu<=fw||w==xv=w;fv=fw;w=u;fw=fu;elseif fu<=fv||v==x||v==w v=u;fv=fu;endendendendlanda=x;Y=X+x*direction;fY=Fx(x,X,direction);function[mySign]=MySign(a,b)if b>0Result=abs(a);elseResulr=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=zeros(1,256);if max(max(a))~=min(min(a))a=(a-min(min(a)))/(max(max(a))-min(min(a)));elsea=zeros(M,N);endif max(max(b))~=min(min(b))b=(b-min(min(b)))/(max(max(b))-min(min(b)));elseb=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;tform=maketform('affine',T);coordinate_x=zeros(width,height);coordinate_y=zeros(width,height);for i=1:widthfor j=1:heightcoordinate_x(i,j)=i;endendfor i=1:widthfor j=1:heightcoordinate_y(i,j)=j;endend[w z]=tforminv(tform,coordinate_x,coordinate_y);for i=1:widthfor j=1:heightsource_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;elsem=fix(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_a2,index_b)+(1-dx)*dy;hab(index_a3,index_b)=hab(index_a3,index_b)+dx*dy;endendendhabsum=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).*log2(pb(index))));mi=Ha+Hb-Hab;4.2.4 Register.mfunction[RegistrationImage]=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,height));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:widthfor j=1:heighttx(i,j)=i;endendfor i=1:widthfor j=1:heightty(i,j)=j;endend[w z]=tforminv(tform,tx,ty);for i=1:widthfor j=1:heightsource_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);elseif(source_x/double(uint16(source_x))==1.0&&(source_y/double(uint16(so urce_y)==1.0)))new_J(i,j)=J(int16(source_x),int16(source_y));elsea=double(uint16(source_x));b=double(uint16(source_y));x11=double(J(a,b));x12=double(J(a,b+1));x21=double(J(a+1,b));x22=double(J(a+1,b+1));new_J(i,j)=uint8((b+1-source_y)*((source_x-a)*x21+(a+1-source_x)*x11) +(source_y-b)*((source_x-a)*x22+(a+1-source_x)*x12));endendendendJ=new_J;I=uint8(I);J=uint8(J);RegistrationImage=uint8(J);5.配准结果参考图像浮动图像配准结果变换参数及最大互信息值X,Y,Angle=[-1.35741][1.39364][9.3372]My Value=[1.9907]Elasped time=[473.458]6.配准结果从配准结果可以看出，选用powell算法及一维brent法的情况下。