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

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

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

评论将由博主筛选后显示,对所有人可见 | 还能输入1000个字符 “速评一下”
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页