# highspeedlogic算法仿真-匹配的两幅相邻图像

close all

image1 = i1(:,:,1);

image2 = i2(:,:,1);

image3 = i3(:,:,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%显示要匹配的两幅相邻图像

%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 交替为红色、白色、蓝色和黑色

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

11-02 199

06-13 4万+
09-02 1170
07-14
11-01 3万+
03-24 2362
11-29 3768
09-24 161万+
08-26 1万+
10-28 3万+