互信息的圖像配準

1、信息論大作業(yè)基于互信息的圖像配準 班級： 09030901 學號： 2009302311 姓名： 益琛 同組成員：陳升富 黎照1. 引言隨著醫(yī)學、計算機技術及生物工程技術的發(fā)展，醫(yī)學影像學為臨床診斷提供了多種模態(tài)的醫(yī)學圖像，不同的醫(yī)學圖像提供了相關臟器的不同信息：CT(Computed Tomography，電子計算機X 射線斷層掃描)和MRI(Magneticresona nce ima ging，核磁共振成像)以較高的空間分辨率提供了臟器的解剖結構信息。在實際臨床應用中，單一模態(tài)的圖像往往不能提供醫(yī)生所需要的足夠的信息，通常需要將不同模態(tài)的圖像融合在一起，得到更豐富的信息，以便了解病變組織

2、或器官的綜合信息，從而做出準確的診斷或制訂出合適的治療方案。而圖像配準是圖像融合的重要前提，圖像配準是指對一幅圖像進行一定的幾何變換而映射到另一幅圖像中，使得兩幅圖像中的相關點達到空間上的一致。圖像配準主要有兩大類方法，基于灰度的方法和基于特征的方法?；诨叶鹊呐錅史椒ㄖ苯永脠D像的灰度數(shù)據(jù)進行配準，從而避免了因分割而帶來的誤差，因而具有精度較高、魯棒性強、不需要預處理而能實現(xiàn)自動配準的特點。在基于灰度的配準方法中，基于互信息的方法包括互信息和歸一化互信息方法，它們已經被廣泛使用并具有最高的精度。本文使用的是基于互信息的配準方法。2. 圖像配準技術2.1圖像配準技術的數(shù)學定義 數(shù)字圖像可以用一

3、個二維矩陣來表示，如果用、分別表示待配準圖像和參考圖像在點(x,y)處的灰度值，那么圖像、的配準關系可表示為： （1）其中f代表二維的空間幾何變換函數(shù)；g表示一維的灰度變換函數(shù)。配準的主要任務是尋找最佳的空間變換關系f與灰度變換關系g，使兩幅圖像實現(xiàn)最佳對準。其中，空間幾何變換是灰度變換的前提，是實現(xiàn)精準配準的關鍵環(huán)節(jié)。2.2幾何變換空間變換主要解決圖像平面上像素的重新定位問題，式(1)中的空間幾何變換函數(shù)f可用空間變換模型進行描述，常用的空間變換模型有剛體變換、仿射變換、投影變換和非線性變換。剛體變換使得一幅圖像中任意兩點間的距離變換到另一幅圖像中后仍然保持不變；仿射變換使得一幅圖像中的直線

4、經過變換后仍保持直線，并且平行線仍保持平行；投影變換是從三維圖像到二維平面的投影；非線性變換把一條直線變換為一條曲線，一般用代數(shù)多項式來表示。仿射變換是最常用的一種空間變換形式,可以實現(xiàn)圖像的平移、旋轉、按比例縮放等操作，我們在實驗中使用的是此變換模型。仿射變換可以用矩陣形式表示： 當分別取值為、將依次對圖像進行平移、旋轉、按比例縮放操作。2.3 插值技術浮動圖的像素點經過空間變換后，參考圖中對應點的坐標一般來說不是整數(shù)，必須通過插值方法計算該點的灰度值。常用的插值算法有最近鄰插值算法、雙線性插值算法和部分體積插值算法。為了盡可能避免基于互信息配準的局部最優(yōu)問題，本文采用改進PV插值算法。PV

5、插值法是一種專門針對兩幅圖像的聯(lián)合直方圖的更新而設計的插值技術，它并不是真正意義上的插值方法，因為通過此方法并不能計算出反向變換點的灰度值。PV插值法的計算過程如圖1.圖中的(x)為反向變換得到的一個浮點數(shù)點，其四個最近鄰像素點分別為。設參考圖像為r(x),浮動圖像為f(x),則它們的聯(lián)合圖方圖函數(shù)如下。 i=1,2,3,41-dxdxdy1-dy2.4 優(yōu)化算法常用的優(yōu)化算法有：牛頓法、最速下降法、模擬退火法、遺傳算法、單純形法、模式搜索法、Powell法等搜索算法。Powell法不需要對目標函數(shù)進行求導計算，具有收斂速度快、精度高、可靠性好等優(yōu)點，是目前解無約束最優(yōu)化問題十分有效的直接法，

6、應用相當廣泛，所以我們在實驗中采用該算法。Powell算法實現(xiàn)如下：(1) 給定允許誤差，初始點和n個線性無關的方向 ，置k=1.(2) 置，從出發(fā)，依次沿方向進行一維搜索，得到點。再從出發(fā)，沿方向作一維搜索點，得到點。（3）若，則停止搜索，得到點；否則，置 返回步驟（2）.2.4.2 Powell算法中的一維搜索算法brent方法。Brent法思路：開始時利用黃金分割法確定一個較小的包含極小點的不確定區(qū)間，然后利用拋物線法獲得一個極小點，若此極小點落在此不確定區(qū)間，則利用該極小點繼續(xù)進行二次插值；否則放棄該點，改用黃金分割法搜索。算法中密切關注a,b,u,v,w,x這六個點，其中a,b表示包

7、含極小點的不確定區(qū)間；u表示最新搜索到的極小點；w表示上一次搜索到的極小點；v表示上一次的w值；x表示當前已搜到的最佳極小點。算法實現(xiàn)步驟如下（設目標函數(shù)為f(x)）：(1)給定初始區(qū)間，精度要求，黃金分割系數(shù)（2）計算，置；計算，置；置上一次迭代步長。（3）計算當前區(qū)間中點，若，則停止搜索，的極小值,否則轉（4）。（4）令，若，則采用黃金分割法，轉（8）。（5）若，則采用黃金分割法，轉（8）。（6）過三點構造拋物線函數(shù)，計算（7）若在之外，則用黃金分割法重新求極小點，轉步驟（8）；若u相對于的改變量大于上一次的改變量，則轉步驟（8）；若或，則用代替前面的改變量。（8）按黃金分割法確定點,且u

8、在區(qū)間和中長度較大的一個；若u相對于x的改變量小于，則用代替前面的改變量。（9）計算，按照各自的定義更新。置，轉步驟（3）。3 基于互信息的圖像配準方法3.1 互信息的計算互信息是信息理論中的一個基本概念，通常用于描述兩個系統(tǒng)間的統(tǒng)計相關性，或者是一個系統(tǒng)中所包含的另一個系統(tǒng)中信息的多少，它可以用熵來描述： （2） 其中，和分別是系統(tǒng)A和B的熵，是它們的聯(lián)合熵，依次定義如下： (3) (4) (5)其中和分別是系統(tǒng)A和B完全獨立時的的概率分布。是系統(tǒng)A和B的聯(lián)合概率分布。令圖像A和B的互信息為，將式（3），（4），（5），分別代入式（2），即可得到圖像互信息的計算公式： 3.2 配準方法首先根

9、據(jù)兩幅圖像的基本情況預設一個初始參數(shù)，其中為裁剪旋轉角的圖像2 行的第一個索引。為裁剪旋轉角的圖像2 列的第一個索引，為旋轉的角度，為比例因子。然后按照給定的初始參數(shù)對圖像2 進行變換，并計算圖像1 和圖像2 的互信息，然后利用最優(yōu)化工具箱中的fminsearch 函數(shù)在附近尋找使圖像1 和圖像2 互信息最大的點，直至搜索到滿足精度要求的參數(shù)；最后輸出配準參數(shù)。4. 圖像配準的實現(xiàn)4.1配準流程首先對參考圖像和浮動圖像按照給定的初始點使用PV插值法統(tǒng)計聯(lián)合直方圖并計算互信息值；然后利用POWELL算法依據(jù)最大互信息理論判斷所得參數(shù)是否最優(yōu)，若不是，則繼續(xù)搜索較優(yōu)參數(shù)，在搜索時會不斷重復“空間幾

10、何變換（affine）-統(tǒng)計聯(lián)合直方圖（PV插值法）-計算互信息值-最優(yōu)化判斷”的過程，直至搜索到滿足精度要求的參數(shù)；最后輸出配準參數(shù)。輸入參考圖像輸入浮動圖像設置初始點和初始搜素方向空間幾何變換計算互信息值最優(yōu)化否是輸出配準參數(shù)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