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 = 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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
e = Inew-Iold;
e2 = e*e+e2;%最小平方差函数E
end
end
if (e1-e2)>0
m=mmm;
l = l/10;
p=1;
q=0;
else
l = l*10;
p=0;
end
end
end
end
w=0;
%AM = m
%BM = [cos(RecordT/180*pi);sin(RecordT/180*pi);-dx;-sin(RecordT/180*pi);cos(RecordT/180*pi);-dy]
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%显示匹配的区域 考虑刚体变换模型 全图像显示
M1 = [cos(RecordT/180*pi),sin(RecordT/180*pi),-dx;-sin(RecordT/180*pi),cos(RecordT/180*pi),-dy;0,0,1];
T = eye(size(M1));
T = T/M1;%T是M1的逆矩阵
%根据精细匹配的结果进行拼接
OffsetY = 100;
OffsetX = 50;
%将第一幅图像放在（OffsetY，OffsetX）位置处
newRoiTar = zeros(2*height,2*width);
for j=1:height
for i=1:width
newRoiTar(j+OffsetY,i+OffsetX) = image1(j,i);
end
end
d = 0.0;
for j=1:height
for i=1:width
C = [i;j;1];
C = T*C;
xNew = C(1);%将图像2上的x坐标映射到图像1上
yNew = C(2);%将图像2上的y坐标映射到图像1上
k = yNew+OffsetY;
wid = width-T(1,3);
if xNew>=1 && yNew>=1 && yNew<=height && xNew<=width
d = i/wid;
a = double(image1(uint16(yNew),uint16(xNew)));
b = double(image2(j,i));
if a==0
d = 1.00;
elseif b==0
d = 0.00;
end
temp = (1-d)*a+d*b;
newRoiTar(uint16(k),uint16(xNew+OffsetX)) = temp;
else
d = 1.00;
newRoiTar(uint16(k),uint16(xNew+OffsetX)) = d*double(image2(j,i));
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%拼接第三幅图像
RecordI = 162;
RecordJ = 4;
RecordT = 0;
dx = Pdx-RecordI;
dy = (Pdy/2+1)-RecordJ;
M2 = [cos(RecordT/180*pi),sin(RecordT/180*pi),-dx;-sin(RecordT/180*pi),cos(RecordT/180*pi),-dy;0,0,1];
T = eye(size(M2*M1));
T = T/(M2*M1);%T是M1的逆矩阵
d = 0.0;
for j=1:height
for i=1:width
C = [i;j;1];
C = T*C;
xNew = C(1);%将图像3上的x坐标映射到图像1上
yNew = C(2);%将图像3上的y坐标映射到图像1上
k = yNew+OffsetY;
wid = width-T(1,3);
if xNew>=1 && yNew>=1 && yNew<=height && xNew<=width
d = i/wid;
a = double(image1(uint16(yNew),uint16(xNew)));
b = double(image3(j,i));
if a==0
d = 1.00;
temp = (1-d)*a+d*b;
elseif b==0
d = 0.00;
temp = (1-d)*a+d*b;
else
temp = (1-d)*a+d*b;
end
newRoiTar(uint16(k),uint16(xNew+OffsetX)) = temp;
else
%可能存在非空点，应增加判断
d = 1.00;
%if newRoiTar(uint16(k),uint16(xNew+OffsetX))==0
newRoiTar(uint16(k),uint16(xNew+OffsetX)) = d*double(image3(j,i));
%else
%newRoiTar(uint16(k),uint16(xNew+OffsetX)) = (newRoiTar(uint16(k),uint16(xNew+OffsetX))+d*double(image3(j,i)))/2;
%end
end
end
end
figure(7),
imshow(uint8(newRoiTar));
figure(8),
imshow(uint8(newRoiTar));
colorbar
%增加颜色显示
colormap(copper)
%rgbplot(hot)%显示颜色变化曲线
%colormapeditor%编辑颜色
%G=hot;
%hsv 色彩饱和值（以红色开始和结束）
%hot 从黑到红到黄到白
%cool 青蓝和洋红的色度
%pink 粉红的彩色度
%gray 线性灰度
%bone 带一点蓝色的灰度
%jet hsv的一种变形（以蓝色开始和结束）
%copper 线性铜色度 好的
%prim 三棱镜。交替为红色、橘黄色、黄色、绿色和天蓝色 ?????
%flag 交替为红色、白色、蓝色和黑色
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
