最佳缝合线算法（图像融合）

理论根据《图像拼接的改进算法》，据说这个算法可以消除重叠部分运动物体的重影和鬼影问题，所以就编下试试看，反正之前编的那种很老的取平均值法融合、渐入渐出（基于距离）融合、以及改进的三角函数权重融合都只是适合静态图像融合 有重影和鬼影问题  不适合有运动物体的图像融合，所以还是要最佳缝合线算法：看论文上就四步 很清晰也很好懂  结果自己写的时候才发现看起来很容易的也许编起来没那么容易 之前 想得太简单了。

缝合线算法：

function D=bestlinefusion(A,B)
%《图像拼接的改进算法》最佳缝合线算法 图像融合
%先根据之前得到的H矩阵计算重叠区域Rect
[H,W,k]=size(A);
rdata1=-118;
L=W+1+rdata1;
R=W;
n=R-L+1;
%计算得到重叠区域的差值图像  其实我不懂计算差值图像干嘛  只要计算重叠区域的大小就好了 为什么还要计算差值图 后面又没用到
Rect=zeros(H,n);
for i=1:H
    for j=L:R
        Rect(i,j-L+1)=A(i,j)-B(i,j-L+1);
    end
end
Rect=uint8(Rect);%这句要不要呢？
%最终融合图的大小
rdata1=-118;
rdata2=3;
Y=2*W+rdata1+1;
D=zeros(H,Y);
%放路径的矩阵
path=zeros(H,n);
%放强度值 每条路径的强度值strength=color^2+geometry
color=zeros(1,n);
geometry=zeros(1,n);
strength1=zeros(1,n);
strength2=zeros(1,n);
%计算第一行即初始化的强度值
for j=L:R
    y=j-L+1;
    color(y)=A(1,j)-B(1,y);
    if(y==1)
        Bxdao=B(1,y+1)+2*B(2,y+1);
        Bydao=B(2,y)+2*B(2,y+1);
        Aydao=2*A(2,j-1)+A(2,j)+2*A(2,j+1);
        Axdao=A(1,j+1)+2*A(2,j+1)-A(1,j-1)-2*A(2,j-1);
        geometry(y)=(Axdao-Bxdao)*(Aydao-Bydao);
        strength1(y)=color(y)^2+geometry(y);
        path(1,y)=y;
        continue;
    end
    if(j==R)
        Axdao=A(1,j-1)-2*A(2,j-1);
        Aydao=2*A(2,j-1)+A(2,j);
        Bxdao=B(1,y+1)+2*B(2,y+1)-B(1,y-1)-2*B(2,y-1);
        Bydao=2*B(2,y-1)+B(2,y)+2*B(2,y+1);
        geometry(y)=(Axdao-Bxdao)*(Aydao-Bydao);
        strength1(y)=color(y)^2+geometry(y);
        path(1,y)=y;
        continue;
    end
    Axdao=A(1,j+1)+2*A(2,j+1)-A(1,j-1)-2*A(2,j-1);
    Bxdao=B(1,y+1)+2*B(2,y+1)-B(1,y-1)-2*B(2,y-1);
    Aydao=2*A(2,j-1)+A(2,j)+2*A(2,j+1);
    Bydao=2*B(2,y-1)+B(2,y)+2*B(2,y+1);
    geometry(y)=(Axdao-Bxdao)*(Aydao-Bydao);
    strength1(y)=color(y)^2+geometry(y);
    path(1,y)=y;
end
color=zeros(1,n);
geometry=zeros(1,n);
small=0;
%开始扩展 向下一行 从第二行到倒数第二行 最后一行单独拿出来 像第一行一样 因为它的结构差值geometry不好算
for i=2:H-1
    %先把下一行的强度值全部计算出来 到时候需要比较哪三个就拿出哪三个
    for j=L:R
        x=i;
        y=j-L+1;
        color(y)=A(i,j)-B(x,y);
        if(y==1)
            Axdao=2*A(i-1,j+1)+A(i,j+1)+2*A(i+1,j+1)-2*A(i-1,j-1)-A(i,j-1)-2*A(i+1,j-1);
            Bxdao=2*B(x-1,y+1)+B(x,y+1)+2*B(x+1,y+1);
            Aydao=-2*A(i-1,j-1)-A(i-1,j)-2*A(i-1,j+1)+2*A(i+1,j-1)+A(i+1,j)+2*A(i+1,j+1);
            Bydao=-B(x-1,y)-2*B(x-1,y+1)+B(x+1,y)+2*B(x+1,y+1);
            geometry(y)=(Axdao-Bxdao)*(Aydao-Bydao);
            strength2(y)=color(y)^2+geometry(y);
            continue;
        end
        if(j==R)
            Axdao=-2*A(i-1,j-1)-A(i,j-1)-2*A(i+1,j-1);
            Bxdao=2*B(x-1,y+1)+B(x,y+1)+2*B(x+1,y+1)-2*B(x-1,y-1)-B(x,y-1)-2*B(x+1,y-1);
            Aydao=-2*A(i-1,j-1)-A(i-1,j)+2*A(i+1,j-1)+A(i+1,j);
            Bydao=-2*B(x-1,y-1)-B(x-1,y)-2*B(x-1,y+1)+2*B(x+1,y-1)+B(x+1,y)+2*B(x+1,y+1);
            geometry(y)=(Axdao-Bxdao)*(Aydao-Bydao);
            strength2(y)=color(y)^2+geometry(y);
           continue;
        end
        Axdao=2*A(i-1,j+1)+A(i,j+1)+2*A(i+1,j+1)-2*A(i-1,j-1)-A(i,j-1)-2*A(i+1,j-1);
        Bxdao=2*B(x-1,y+1)+B(x,y+1)+2*B(x+1,y+1)-2*B(x-1,y-1)-B(x,y-1)-2*B(x+1,y-1);
        Aydao=-2*A(i-1,j-1)-A(i-1,j)-2*A(i-1,j+1)+2*A(i+1,j-1)+A(i+1,j)+2*A(i+1,j+1);
        Bydao=-2*B(x-1,y-1)-B(x-1,y)-2*B(x-1,y+1)+2*B(x+1,y-1)+B(x+1,y)+2*B(x+1,y+1);
        geometry(y)=(Axdao-Bxdao)*(Aydao-Bydao);
        strength2(y)=color(y)^2+geometry(y);
    end
    for j=1:n
        if(path(i-1,j)==1)
           if(strength2(1)<strength2(2))
              strength1(j)=strength1(j)+strength2(1);
              path(i,j)=1;
           else
              strength1(j)=strength1(j)+strength2(2);
              path(i,j)=2;
           end
        else
            if(path(i-1,j)==n)
              if(strength2(n-1)<strength2(n))
                strength1(j)=strength1(j)+strength2(n-1);
                path(i,j)=n-1;
              else
                strength1(j)=strength1(j)+strength2(n);
                path(i,j)=n;
              end
            else
                if(strength2(path(i-1,j)-1)<strength2(path(i-1,j)))
                    if(strength2(path(i-1,j)-1)<strength2(path(i-1,j)+1))
                        small=strength2(path(i-1,j)-1);
                        path(i,j)=path(i-1,j)-1;
                    else
                        small=strength2(path(i-1,j)+1);
                        path(i,j)=path(i-1,j)+1;
                    end
                else
                    if(strength2(path(i-1,j))<strength2(path(i-1,j)+1))
                        small=strength2(path(i-1,j));
                        path(i,j)=path(i-1,j);
                    else
                        small=strength2(path(i-1,j)+1);
                        path(i,j)=path(i-1,j)+1;
                    end
                end
                strength1(j)=strength1(j)+small;
            end
        end
        small=0;
    end
    strength2=zeros(1,n);
    color=zeros(1,n);
    geometry=zeros(1,n);
end
%单独计算最后一行
i=H;
for j=L:R
        x=i;
        y=j-L+1;
        color(y)=A(i,j)-B(x,y);
        if(y==1)
            Axdao=2*A(i-1,j+1)+A(i,j+1)-2*A(i-1,j-1)-A(i,j-1);
            Aydao=-2*A(i-1,j-1)-A(i-1,j)-2*A(i-1,j+1);
            Bxdao=2*B(x-1,y+1)+B(x,y+1);
            Bydao=-B(x-1,y)-2*B(x-1,y+1);
            continue;
        end
        if(j==R)
            Bxdao=2*B(x-1,y+1)+B(x,y+1)-2*B(x-1,y-1)-B(x,y-1);
            Bydao=-2*B(x-1,y-1)-B(x-1,y)-2*B(x-1,y+1);
            Axdao=-2*A(i-1,j-1)-A(i,j-1);
            Aydao=-2*A(i-1,j-1)-A(i-1,j);
            continue;
        end 
        Axdao=2*A(i-1,j+1)+A(i,j+1)-2*A(i-1,j-1)-A(i,j-1);
        Bxdao=2*B(x-1,y+1)+B(x,y+1)-2*B(x-1,y-1)-B(x,y-1);
        Aydao=-2*A(i-1,j-1)-A(i-1,j)-2*A(i-1,j+1);
        Bydao=-2*B(x-1,y-1)-B(x-1,y)-2*B(x-1,y+1);
        geometry(y)=(Axdao-Bxdao)*(Aydao-Bydao);
        strength2(y)=color(y)^2+geometry(y);
end
for j=1:n
        if(path(i-1,j)==1)
           if(strength2(1)<strength2(2))
              strength1(j)=strength1(j)+strength2(1);
              path(i,j)=1;
           else
              strength1(j)=strength1(j)+strength2(2);
              path(i,j)=2;
           end
        else
            if(path(i-1,j)==n)
              if(strength2(n-1)<strength2(n))
                strength1(j)=strength1(j)+strength2(n-1);
                path(i,j)=n-1;
              else
                strength1(j)=strength1(j)+strength2(n);
                path(i,j)=n;
              end
            else
                if(strength2(path(i-1,j)-1)<strength2(path(i-1,j)))
                    if(strength2(path(i-1,j)-1)<strength2(path(i-1,j)+1))
                        small=strength2(path(i-1,j)-1);
                        path(i,j)=path(i-1,j)-1;
                    else
                        small=strength2(path(i-1,j)+1);
                        path(i,j)=path(i-1,j)+1;
                    end
                else
                    if(strength2(path(i-1,j))<strength2(path(i-1,j)+1))
                        small=strength2(path(i-1,j));
                        path(i,j)=path(i-1,j);
                    else
                        small=strength2(path(i-1,j)+1);
                        path(i,j)=path(i-1,j)+1;
                    end
                end
                strength1(j)=strength1(j)+small;
            end
        end
        small=0;
end        
%比较strength1里放的每条路径的强度值的总和 谁最小 就选path中对应的那一列的路径
[minzhi,minth]=min(strength1);
mypath=path(:,minth);
%mypath放的就是最佳缝合线选出的路径 这条路径坐标是参考图A 右边是目标图B
for i=1:H
    for j=1:mypath(i)+L-1
        D(i,j,1)=A(i,j,1);
        D(i,j,2)=A(i,j,2);
        D(i,j,3)=A(i,j,3);
    end
    for j=mypath(i)+L-1:Y-1
         x=i;
         y=j-L+1;
         D(i,j,1)=B(x,y,1);
         D(i,j,2)=B(x,y,2);
         D(i,j,3)=B(x,y,3);
    end
end

 D=uint8(D);

还是以那两幅图为例：

结果：

差值图像Rect：

可是最终得到的融合后的图像D怎么是下面这个鬼样子：这就是选出来的最佳缝合线？？戳瞎我的狗眼啊！这哪里最佳了？？！！

貌似理论部分没编错啊 哪里有问题呢  还是运行最佳缝合线本来就应该是这样？为了验证下  在另一篇也讲最佳缝合线算法的论文里 《全景图像拼接关键技术研究》下载了两张图片，看找出来最佳缝合线是否和这个作者的近似或者一样

原图：通过放进OPENCV求匹配对 有32对匹配对 然后将这些匹配对输出给MATLAB计算H变换矩阵 我用了30对坐标 得到H矩阵 再进行H矩阵变换 后：然后用下渐入渐出融合 结果：用下改进的三角函数权重的融合 结果果然改进的肉眼看不出来效果  但它们两个都有重影 比较严重的重影！用下刚刚编的最佳缝合线算法看看 这是差值图像然后这是最佳缝合线得到的D：可以看到用最佳缝合线后 没有重影了  那棵树那里以及那棵大树旁边两棵树没有重影了  可以和渐入渐出以及改进的那个对比下   这个已经没有重影了   原来最佳缝合线还是有用的喔   

我想把这幅图中找到的最佳缝合线画出来 ，其实我放大后看得到那条最佳缝，我想想怎么画出来：

L=W+1+rdata1;
imshow(D)
hold on;
for i=1:137
     D(i,L+mypath(i),1)=255;
     D(i,L+mypath(i),2)=255;
     D(i,L+mypath(i),3)=255;
end
hold off;

这条白线就是用上面的最佳缝合线程序找到的最佳缝合线  感觉好像一道白色的闪电 上面看不太出来  用黑色画下：

这样就看得清了 这就是找到的最佳缝合线  的确是避开了那些重影物体  我是这样来写的 其实就是和上面的H矩阵不一样 这个计算出来平移距离是76  竖直方向其实有5个像素的平移量  但我懒得改上面的了   直接只改了rdata1=-76 所以和那个论文中找的最佳缝合线有点点差别  不过没多大关系  因为我指望着像《图像拼接的改进算法》里说的那样找到最佳缝合线后还要进行多分辨率拼接的   就像之前我编渐入渐出时初步得到H矩阵变换后的结果也是有点没对齐的，但是融合时用渐入渐出后它就完完全全对齐了   这个也是一样的   我猜想后续的多分辨率拼接也可以做到这样的效果   虽然现在得到的最佳缝合线有点没对齐的样子  因为我没用竖直方向的平移量rdata2=5，其实加上一样的   但太麻烦了因为我再看一次程序 然后重新计算 我下次再去改！现在这样就行了 竖直方向没对齐的我就依靠后续的多分辨率算法好了   但当竖直方向的平移量rdata2很大时以及有旋转关系时 也要重新计算  不能只算一个rdata1   绝对不能这样  会影响后续多分辨率的效果！