泊松融合-matlab实现

泊松图像融合原理很容易理解但是实现的时候有很多坑，在找资料的时候发现这方面的代码很少，虽然修改多次但代码仍有很多不足的地方。

如果想了解泊松融合原理可以看一看本站其他一些文章写的很好。

废话不说上代码：

main.m

clc,clear
a=double(imread('2.jpg'));
b=double(imread('1.jpg'));
[row_a,col_a,g_a]=size(a);%目标文件
[row_b,col_b,g_b]=size(b);%背景文件
mask_A=zeros(row_a,col_a);
imshow(uint8(b));
flat=1;
while flat==1
    [y_start,x_start]=ginput(1);
    x_end=x_start+row_a-1;
    y_end=y_start+col_a-1;
    flat=0;
    if x_end>row_b||y_end>col_b
        disp('Array out of line');
        flat=1;
    end
end
mask_B=ones(row_b,col_b);
mask_B(x_start:x_end,y_start:y_end)=0;
global number
number=1;
%%
[A,B]=create_AB(a(:,:,1),mask_A,b(:,:,1),mask_B);
x_R= A\B;
%%
[~,B]=create_AB(a(:,:,2),mask_A,b(:,:,2),mask_B);
x_G= A\B;
%%
[~,B]=create_AB(a(:,:,3),mask_A,b(:,:,3),mask_B);
x_B= A\B;
%%
[hxb,lxb,l]=find(mask_B==0);%他是按照一列列进行遍历的
s=[hxb';lxb'];
[s1,id1]=sort(s(1,:));
s_z=[s1;s(2,id1)];
hxb=s_z(1,:);
lxb=s_z(2,:);
for i=1:length(l)
    b(hxb(i),lxb(i),1)=x_R(i);
    b(hxb(i),lxb(i),2)=x_G(i);
    b(hxb(i),lxb(i),3)=x_B(i);
end
imshow(uint8(b))

create_AB.m

function [A,B]=create_AB(image_A,mask_A,image_B,mask_B)
%假设插入的图片为方形
%image_A为插入图片
%mask_A为背景图片的模子，待插入部分为0，不插入部分为1
global number
[hxb,lxb,~]=find(mask_B==0);
B=zeros(length(hxb),1);%创建B
s=[hxb';lxb'];
[s1,id1]=sort(s(1,:));
s_z=[s1;s(2,id1)];%将行下标和列下标的顺序整理为行优先
hxb=s_z(1,:);
lxb=s_z(2,:);
min_x=min(hxb);
min_y=min(lxb);
max_x=max(hxb);
max_y=max(lxb);
k=2;
for i=1:length(hxb)
    mask_B(hxb(i),lxb(i))=k;
    k=k+1;
end
if number==1
A=sparse(zeros(length(hxb),length(hxb)));
end
if number~=1
    A=0;
end
k=1;%k为正在填写A矩阵的第k行
for i=min_x:max_x
    for j=min_y:max_y
        if i==min_x||i==max_x||j==min_y||j==max_y
            if number==1%如果第一次运行该求A函数
                A(k,k)=1;
            end
            B(k)=image_B(i,j);%边界点B为其像素值
            k=k+1;
        end
        if i~=min_x && i~=max_x&&j~=min_y&&j~=max_y
            if number ==1
                A(k,mask_B(i-1,j)-1)=1;
                
                A(k,mask_B(i,j-1)-1)=1;
                
                A(k,mask_B(i,j+1)-1)=1;
                
                A(k,mask_B(i+1,j)-1)=1;
                
                A(k,mask_B(i,j)-1)=-4;
            end
            k=k+1;
        end
    end
end
%%
%求散度,将B补充完全
[A_hxb,A_lxb,~]=find(mask_A==0);
min_Ax=min(A_hxb);
min_Ay=min(A_lxb);
max_Ax=max(A_hxb);
max_Ay=max(A_lxb);
k=1;
for i=min_Ax:max_Ax
    for j=min_Ay:max_Ay
        if i==min_Ax||i==max_Ax||j==min_Ay||j==max_Ay
            k=k+1;
            continue;
        end
        B(k)=image_A(i-1,j)+image_A(i,j-1)+image_A(i,j+1)+image_A(i+1,j)-4*image_A(i,j);%非边界点B为散度
        k=k+1;
    end
end
number=number+1;
end

首先说一说遇到的坑1.导入图片的时候要将unit8型转化为doule型，在写入插入后图片时要转化为unit8型2.AX=B中X为要插入的点的像素值，主要也是要求它，其中A为行数为length（插入像素点）列数也为 length（插入像素点）的矩阵，A分为两部分，第一部分为该行所代表的像素点为边界点时，这时此行所对应的B为该边界点的像素值，第二部分为改行所代表的像素点为非边界点时，此时B所对应的为该点的散度值。3.用[hxb,lxb,~]=find(mask==0)查找要插入部分和待插部分时，返回的行下表和列下标也就是第一个和第二个返回值，它是以列优先排列的需要将他做一个变换为行优先。

在说一说代码的需要改进之处1.首先他只能融合方形的图片，因为我在编写代码时优点偷懒只想完成老师留的作业没有考虑到     代码的可扩展型。2.如你们从上面看到的那样我将A,B放在一个函数中求了，如果是三通道图片，A将计算三次（事实上只需计算一次即可），好在我进行了一点小小的改进，设置了一个全局计数器用来统计 create_AB运行了几次也就是A被计算了几次，如果运行一次以上我就不在运行计算A的代码 了，这样大量的计算资源不会被浪费，即使这样程序还是有很多处地方应去改善。

  ps：如果做实际用途的话可以使用opencv中现成的API，融合效果与执行效率比我这个要好

融合效果：

插入图片：

背景图片：

选取插入的坐标值

融合效果：