close all
i1 = imread('10031.bmp');
image1 = i1(:,:,1);
i2 = imread('10036.bmp');
image2 = i2(:,:,1);
i3 = imread('10041.bmp');
image3 = i3(:,:,1);
%image1 = imread('D:\matlab\work\10037.bmp');
%image2 = imread('D:\matlab\work\10041.bmp');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%显示要匹配的两幅相邻图像
%figure(1);
%imshow(image1);
%figure(2);
%imshow(image2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
height = size(image1,1);%基准图像高
width = size(image1,2);%基准图像宽
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%取匹配模板
Pdy = 10;%匹配块在y方向上的偏移量
newHeight = (double(uint8(height/10))*10-Pdy);
newWidth = (newHeight/10);
Pdx = double(uint8((width-newWidth)/2));%匹配块在x方向上的偏移量
RoiSorce = zeros(newHeight,newWidth);%新匹配块数组
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%取匹配块并赋初值
for j=(Pdy/2+1):newHeight+(Pdy/2)
for i=1:newWidth
RoiSorce((j-Pdy/2),i) = image1(j,(Pdx+i));
end
end
%figure(3),
%imshow(uint8(RoiSorce));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SSDvector = zeros(1,(width-newWidth));
%粗匹配 未考虑旋转
minSSD = 10000;
sumSSD = 0;
RecordI = 0;%记录X位移量
RecordJ = 0;%记录Y位移量
%改成二重循环
%K = 1;
q=0;%控制是否进行粗匹配算法，因为算法运行时间长，图像不改变不需要重新计算
while q==1
for z=1:(height-newHeight)
for k=1:(width-newWidth)
sumSSD = 0;
for j=1:newHeight
for i=1:newWidth
temp1 = double(image2(z-1+j,k-1+i));
temp2 = double(RoiSorce(j,i));
temp = temp1-temp2;
m = temp*temp;
sumSSD = sumSSD+m;
end
end
sumSSD = sumSSD/(height*width);
if sumSSD<minSSD
minSSD = sumSSD;
RecordI = k;
RecordJ = z;
end
%SSDvector(K) = sumSSD;
%K=K+1;
end
end
q=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%显示粗匹配的区域 未考虑旋转
%for j=1:newHeight
%for i=1:newWidth
%RoiTar(j,i) = image2(RecordJ+j,RecordI+i);
%end
%end
%figure(4),
%imshow(uint8(RoiTar));
%figure(5)
%plot(SSDvector);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%考虑旋转分量
%test US new 31&36
RecordI = 159;
RecordJ = 4;
minSSD = 59.0408;
%end of test
Ox = double(uint8(newWidth/2));%匹配块中心点x坐标
Oy = newHeight/2;%匹配块中心点y坐标
RecordT = 0;%记录角度偏移量
%考虑旋转角度不大于正负2.50度
for t=-2.50:0.1:2.50%由右向左旋转
tg = tan(t/180*pi);
sumSSD = 0;
for j=1:newHeight
xR = (Oy-j)*tg;
if xR<0
xR = double(uint16(-xR));
xR = -xR-1;
else
xR = double(uint16(xR));
end
for i=1:newWidth
temp1 = double(image2(RecordJ+j,RecordI+i+xR));
temp2 = double(RoiSorce(j,i));
temp = temp1-temp2;
m = temp*temp;
sumSSD = sumSSD+m;
end
end
sumSSD = sumSSD/(height*width);
if sumSSD<minSSD
minSSD = sumSSD;
RecordT = t;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%精细匹配算法
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w=0;%判断是否调用精细匹配
dx = Pdx-RecordI;
dy = (Pdy/2+1)-RecordJ;
%刚体变换仿射矩阵参数
m = [cos(RecordT/180*pi);sin(RecordT/180*pi);-dx;-sin(RecordT/180*pi);cos(RecordT/180*pi);-dy];
while w==1
mm = zeros(6,1);
mmm = zeros(6,1);
A = zeros(6,6);
AA = zeros(6,6);
deta = eye(size(A));
em = zeros(6,1);
bb = zeros(6,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jj = height-2;%防止越界
ii = width-2;%防止越界
l = 0.001;%l不变化时优选值
e1 = 0;
e2 = 0;
q=0;%迭代次数
p=0;%迭代次数
while q==0
e1 = 0;
p=0;
for j=(Pdy/2+1):newHeight+(Pdy/2)
for i=(Pdx+1):(Pdx+newWidth)
%计算匹配点x1，y1 双线性插值
x1 = m(1)*i+m(2)*j+m(3);%有可能越界，怎么处理？
y1 = m(4)*i+m(5)*j+m(6);%有可能越界，怎么处理？
if x1<1
x1 = 1;
end
if x1>=width
x1 = ii;
end
if y1<1
y1 = 1;
end
if y1>=height
y1 = jj;
end
x11 = uint8(x1);%x1取整
y11 = uint8(y1);%y1取整
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%做双线性插值
x111 = double(x11);
y111 = double(y11);
b = x1-x111;
a = y1-y111;
t1 = a*double(image2(uint16(y111),uint16(x111+1)))+(1-a)*double(image2(uint16(y111+1),uint16(x111+1)));
t2 = a*double(image2(uint16(y111+1),uint16(x111)))+(1-a)*double(image2(uint16(y111),uint16(x111)));
Iold = double(image1(uint16(j),uint16(i)));
Inew = b*t1+(1-b)*t2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%求灰度插值，梯度值Gx，Gy
e = Inew-Iold;
e1 = e*e+e1;%最小平方差函数E
x22 = uint16(x1+1);%x1取整
y22 = uint16(y1+1);%y1取整
x222 = double(x22);
y222 = double(y11);
%%%%%%%%%%%%%%%%%%
%做双线性插值
t1 = a*double(image2(uint16(y222),uint16(x222+1)))+(1-a)*double(image2(uint16(y222+1),uint16(x222+1)));
t2 = a*double(image2(uint16(y222+1),uint16(x222)))+(1-a)*double(image2(uint16(y222),uint16(x222)));
Inew = b*t1+(1-b)*t2;
Gx = e-Inew;
x222 = double(x11);
y222 = double(y22);
t1 = a*double(image2(uint16(y222),uint16(x222+1)))+(1-a)*double(image2(uint16(y222+1),uint16(x222+1)));
t2 = a*double(image2(uint16(y222+1),uint16(x222)))+(1-a)*double(image2(uint16(y222),uint16(x222)));
Inew = b*t1+(1-b)*t2;
Gy = e-Inew;
%%%%%%%%%%%%%%%%%%
%求e对仿射系数的偏导数
em(1) = Gx*x1;
em(2) = Gx*y1;
em(3) = Gx;
em(4) = Gy*x1;
em(5) = Gy*y1;
em(6) = Gy;
%求Hessian矩阵
for jj=1:6
for ii=1:6
A(jj,ii) = em(jj)*em(ii)+A(jj,ii);
end
end
%求加权梯度向量
for ii=1:6
bb(ii) = -e*em(ii)-bb(ii);
end
end
end
while p==0
%求改进后的Hessian矩阵
AA = (1+deta*l)*A;
%求m增量
mm = deta/AA*bb;
if abs(e1-e2)<0.00001
p=1;
q=1;
else
e2 = 0;
mmm = m+mm;
for j=(Pdy/2+1):newHeight+(Pdy/2)
for i=(Pdx+1):(Pdx+newWidth)
%计算匹配点x1，y1 双线性插值
x1 = mmm(1)*i+mmm(2)*j+mmm(3);%有可能越界，怎么处理？
y1 = mmm(4)*i+mmm(5)*j+mmm(6);%有可能越界，怎么处理？
if x1<1
x1 = 1;
end
if x1>=width
x1 = ii;
end
if y1<1
y1 = 1;
end
if y1>=height
y1 = jj;
end
x11 = uint16(x1);%x1取整
y11 = uint16(y1);%y1取整
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%做双线性插值
x111 = do
