基于互信息的图像配准Word文档格式.doc
《基于互信息的图像配准Word文档格式.doc》由会员分享，可在线阅读，更多相关《基于互信息的图像配准Word文档格式.doc（17页珍藏版）》请在冰点文库上搜索。
dx
dy
1-dy
2.4优化算法
2.4.1常用的优化算法有：
牛顿法、最速下降法、模拟退火法、遗传算法、单纯形法、模式搜索法、Powell法等搜索算法。
Powell法不需要对目标函数进行求导计算，具有收敛速度快、精度高、可靠性好等优点，是目前解无约束最优化问题十分有效的直接法，应用相当广泛，所以我们在实验中采用该算法。
Powell算法实现如下：
（1）给定允许误差，初始点和n个线性无关的方向，置k=1.
（2）置，，从出发，依次沿方向进行一维搜索，得到点。
再从出发，沿方向作一维搜索点，得到点。
（3）若，则停止搜索，得到点；
否则，置
返回步骤
（2）.
2.4.2Powell算法中的一维搜索算法——brent方法。
Brent法思路：
开始时利用黄金分割法确定一个较小的包含极小点的不确定区间，然后利用抛物线法获得一个极小点，若此极小点落在此不确定区间，则利用该极小点继续进行二次插值；
否则放弃该点，改用黄金分割法搜索。
算法中密切关注a,b,u,v,w,x这六个点，其中a,b表示包含极小点的不确定区间；
u表示最新搜索到的极小点；
w表示上一次搜索到的极小点；
v表示上一次的w值；
x表示当前已搜到的最佳极小点。
算法实现步骤如下（设目标函数为f（x））：
（1）给定初始区间，精度要求，黄金分割系数
（2）计算，置；
计算，置；
置上一次迭代步长。
（3）计算当前区间中点，若，则停止搜索，的极小值,否则转（4）。
（4）令，若，则采用黄金分割法，转（8）。
（5）若，则采用黄金分割法，转（8）。
（6）过三点构造抛物线函数，计算
（7）若在之外，则用黄金分割法重新求极小点，转步骤（8）；
若u相对于的改变量大于上一次的改变量，则转步骤（8）；
若或，则用代替前面的改变量。
（8）按黄金分割法确定点,且u在区间和中长度较大的一个；
若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配准方法
首先根据两幅图像的基本情况预设一个初始参数，其中为裁剪
旋转角的图像2行的第一个索引。
为裁剪旋转角的图像2列的第一个索引，为旋转的角度，为比例因子。
然后按照给定的初始参数对图像2进行变换，并计算图像1和图像2的互信息，然后利用最优化
工具箱中的fminsearch函数在附近寻找使图像1和图像2互信息最大的
点，直至搜索到满足精度要求的参数；
最后输出配准参数。
4.图像配准的实现
4.1配准流程
首先对参考图像和浮动图像按照给定的初始点使用PV插值法统计联合直方图并计算互信息值；
然后利用POWELL算法依据最大互信息理论判断所得参数是否最优，若不是，则继续搜索较优参数，在搜索时会不断重复“空间几何变换（affine）-统计联合直方图（PV插值法）-计算互信息值-最优化判断”的过程，直至搜索到满足精度要求的参数；
输入参考图像
输入浮动图像
设置初始点和初始搜素方向
空间几何变换
计算互信息值
最优化
否
是
输出配准参数
4.2.所用到的M文件及其源代码
4.2.1．ImageRegistration.m
functionvarargout=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'
[]）;
ifnargin&
&
ischar（varargin{1}）
gui_State.gui_Callback=str2func（varargin{1}）;
end
ifnargout
[varargout{1:
nargout}]=gui_mainfcn（gui_State,varargin{:
}）;
else
gui_mainfcn（gui_State,varargin{:
addpath（pwd）;
functionImageRegistration_OpeningFcn（hObject,eventdata,handles,varargin）
handles.output=hObject;
guidata（hObject,handles）;
functionvarargout=ImageRegistration_OutputFcn（hObject,eventdata,handles）
varargout{1}=handles.output;
functionpushbutton1_Callback（hObject,eventdata,handles）
globalI;
%%%调用OpenImage.m读入参考图像并获取文件名、图像大小%%%
[filename,pathname]=uigetfile（{'
*.jpg'
;
'
*.bmp'
},'
Ñ
¡
Ô
ñ
Í
¼
Æ
¬
）;
str=[pathnamefilename];
I=imread（str）;
axes（handles.axes1）;
imshow（I）;
handles.data=I;
guidata（hObject,handles）;
figure
（1）;
imshow（handles.data）;
functionpushbutton3_Callback（hObject,eventdata,handles）
handles.Old_I=handles.data;
handles.Old_J=handles.data2;
[I,J]=GLPF（handles）;
handles.data2=J;
tic
RegistrationParameters=Powell（handles）;
toc
ElapsedTime=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（'
ElapsedTime=[%.3f]'
ElapsedTime）;
axes（handles.axes3）
[RegistrationImage]=Register（handles）;
imshow（RegistrationImage）
set（handles.text1,'
string'
RegistrationResult）;
set（handles.text2,'
set（handles.text3,'
functionpushbutton2_Callback（hObject,eventdata,handles）
globalJ;
J=imread（str）;
axes（handles.axes2）;
imshow（J）;
figure
（2）;
imshow（handles.data2）;
4.2.2Powell.m
function[RegistrationParameters]=Powell（handles）
e=0.1;
X0=[000];
D=[100;
010;
001];
num=0;
while（num<
200）
num=num+1;
d1=D（1,:
%d1Î
ª
¾
Ø
Õ
ó
Dµ
Ä
µ
Ú
Ò
»
Ð
£
³
õ
Ê
Ë
÷
·
½
Ï
ò
[X1,fX1]=OneDimSearch（X0,d1,handles）;
d2=D（2,:
%d2Î
¶
þ
[X2,fX2]=OneDimSearch（X1,d2,handles）;
d3=D（3,:
%d3Î
È
ý
Î
¹
Ó
¸
ö
[X3,fX3]=OneDimSearch（X2,d3,handles）;
fX0=PV（X0
（1）,X0
（2）,-X0（3）,handles）;
Diff=[fX1-fX0fX2-fX1fX3-fX2];
[maxDiff,m]=max（Diff）;
%maxº
¯
Ã
¨
maxdiffÎ
Á
¿
Diffµ
×
î
´
mÎ
º
Å
d4=X3-X0;
temp1=X3-X0;
Conditon1=sum（temp1.*temp1）;
ifConditon1<
=e
break
end
[X4,fX4,landa]=OneDimSearch（X0,d4,handles）;
X0=X4;
temp2=X4-X3;
Conditon2=sum（temp2.*temp2）;
ifConditon2<
X3=X4;
temp3=sqrt（（fX4-fX0）/（maxDiff+eps））;
if（abs（landa）>
temp3）
D（4,:
）=d4;
fori=m:
3
D（i,:
）=D（i+1,:
end
RegistrationParameters
（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;
ifa>
b
temp=a;
a=b;
b=temp;
v=a+cgold*（b-a）;
w=v;
x=v;
e=0.0;
fx=Fx（x,X,direction）;
fv=fx;
fw=fx;
foriter=1:
IterTimes
xm=0.5*（a+b）;
ifabs（x-xm）<
=Epsilon*2-0.5*（b-a）
ifabs（e）>
Epsilon
r=（x-w）*（fx-fv）;
q=（x-v）*（fx-fw）;
p=（x-v）*（q-（x-w）*r;
q=2*（q-r）;
ifq>
p=-p;
q=abs（q）;
etemp=e;
e=d;
ifnot
（abs（p）>
=abs（0.5*q*etemp）||p<
=q*（a-x）||p>
=q*（b-x））
d=p/q;
u=x+d;
ifu-a<
Epsilon*2||b-u<
Epsilon*2
d=MySign（Epsilon,xm-x）;
end
else
ifx>
=xm
e=a-x;
else
e=b-x;
d=cgold*e;
else
ifx>
=xm;
e=a-x;
e=b-x;
d=cgold*e;
ifabs（d）>
=Epsilon
u=x+d;
u=x+MySign（Epsilon,d）;
fu=Fx（u,X,direction）;
iffu<
=fx
ifu>
=x
a=x;
b=x;
v=w;
fv=fw;
w=x;
fw=fx;
x=u;
fx=fu;
ifu<
x
a=u;
b=u;
iffu<
=fw||w==x
v=w;
fv=fw;
w=u;
fw=fu;
iffu<
=fv||v==x||v==w
v=u;
fv=fu;
landa=x;
Y=X+x*direction;
fY=Fx（x,X,direction）;
function[mySign]=MySign（a,b）
ifb>
Result=abs（a）;
Resulr=abs（a）;
mySign=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.3PV.m
function[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）;
ifmax（max（a））~=min（min（a））
a=（a-min（min（a）））/（max（max（a））-min（min（a）））;
a=zeros（M,N）;
ifmax（max（b））~=min（min（b））
b=（b-min（min（b）））/（max（max（b））-min（min（b）））;
b=zeros（M,N）;
a=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=[100;
xy1];
t2=[100;
-u-v1];
t3=[cos（rad）-sin（rad）0;
sin（rad）cos（rad）0;
t4=[100;
uv1];
T=t2*t3*t4*t1;
tform=maketform（'
affine'
T）;
coordinate_x=zeros（width,height）;
coordinate_y=zeros（width,height）;
fori=1:
width
forj=1:
height
coordinate_x（i,j）=i;
coordinate_y（i,j）=j;
[wz]=tforminv（tform,coordinate_x,coordinate_y）;
source_x=w（i,j）;
s
升级会员 