泊松融合
2015-11-16 16:02
267 查看
拖拖拉拉快一个月了 这是探索出来的 但感觉效果不是特别好 我目前只能做到这样子了：想把左图的人放进右图的水池中 用泊松融合想实现无缝自然的效果：
左图是直接做mask镶嵌的结果 右图是解了泊松方程后的最终的结果 我所能得到的：
感觉效果不够好 因为这个理想结果是这样：
看这个就比我的自然很多很多 我的应该是有问题的 但我是按照步骤来的 目前不知道错在哪里 在找：[今天经过网友的提醒，找到了，原来是迭代次数问题，我这个效果是只迭代一次的，再迭代几次就能达到理想效果，那篇论文中没有提过这个迭代！害死我了]
A=imread('F:\fisheye\pond.jpg');
src=A;
[ma,na,ka]=size(A);
dst=imread('F:\fisheye\swim.jpg');
se=strel('diamond',10);
B_erode=imerode(dst,se);
Berodelogical=im2bw(B_erode);
imshow(Berodelogical)
B=dst;
[mb,nb,kb]=size(B);
>> for i=1:mb
for j=1:nb
if(Berodelogical(i,j)==1)
Berodelogical(i,j)=0;
else
Berodelogical(i,j)=1;
end
end
end
dstX=100;
dstY=100;
for i=dstY:dstY+mb-1
for j=dstX:dstX+nb-1
ii=i-(dstY-1);
jj=j-(dstX-1);
if(Berodelogical(ii,jj)==1)
A(i,j,1)=B(ii,jj,1);
A(i,j,2)=B(ii,jj,2);
A(i,j,3)=B(ii,jj,3);
end
end
end
ROI=uint8(zeros(mb,nb,3));
for i=dstY:dstY+mb-1
for j=dstX:dstX+nb-1
ii=i-(dstY-1);
jj=j-(dstX-1);
ROI(ii,jj,1)=A(i,j,1);
ROI(ii,jj,2)=A(i,j,2);
ROI(ii,jj,3)=A(i,j,3);
end
end
w=[0,-1,1];
ROIgradienty=imfilter(double(ROI),w,'conv');
ROIgradientx=imfilter(double(ROI),w','conv');
%接下来对梯度求偏导得到融合图像的散度 lap
kernel=[-1,1,0];
lapy=imfilter(ROIgradienty,kernel,'conv');
lapx=imfilter(ROIgradientx,kernel','conv');
lap=lapx+lapy;
%接下来求解系数矩阵
%先将lap四周的点填充为原来的像素值（约束值） 其它的仍然是散度
[m,n,k]=size(lap);
heightRegion=m;
widthRegion=n;
L=n;
J=m;
for i=1:m
for j=1:n
if((i==1&&j==1)||(i==1&&j==n)||(i==m&&j==1)||(i==m&&j==n))
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
else
if(i==1||i==m)
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
else
if(j==1||j==n)
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
end
end
end
end
end
%然后将lap'排成一列 即Ax=b的b矩阵
for i=1:m
for j=1:n
lap1(i,j)=lap(i,j,1);
lap2(i,j)=lap(i,j,2);
lap3(i,j)=lap(i,j,3);
end
end
lapp1=lap1';
lapp2=lap2';
lapp3=lap3';
b1=lapp1(:);
b2=lapp2(:);
b3=lapp3(:);
b=[b1 b2 b3];
%计算系数矩阵Ax=b的A 这里用M
[mm,nn]=size(b);
M=spalloc(mm,mm,6*(mm-2*m-2*n+4));
for i=1:mm
for j=1:mm
if(i==j)
M(i,j)=1;
end
end
end
%上面是边界点的约束值仍然是它本身的像素值 就是乘以1它本身就好
for y=1:heightRegion
for x=1:widthRegion
if x~=1&&y~=1&&y~=heightRegion&&x~=widthRegion
i=(y-1)*L+x;
M(i,i)=-4;
M(i,i-1)=1;
M(i,i+1)=1;
M(i,i-L)=1;
M(i,i+L)=1;
end
end
end
x1=bicg(M,b1,1e-6,400); %用传说中的稀疏矩阵 原来就是为了防止我之前的问题out of memory！当有很多0时的大矩阵
x2=bicg(M,b2,1e-6,400); %这个双共轭梯度解法 来解Mx=b的这个x 400是迭代次数
x3=bicg(M,b3,1e-6,400);
% reshape x to become the pixel intensity of the region
imRegion1= reshape(x1,widthRegion,heightRegion);
imRegion2= reshape(x2,widthRegion,heightRegion);
imRegion3= reshape(x3,widthRegion,heightRegion);
imNew=src;
imNew(dstY:dstY+heightRegion-1,dstX:dstX+widthRegion-1,1)=imRegion1';
imNew(dstY:dstY+heightRegion-1,dstX:dstX+widthRegion-1,2)=imRegion2';
imNew(dstY:dstY+heightRegion-1,dstX:dstX+widthRegion-1,3)=imRegion3';
下面这些是逐步探索的得到上面这个的过程：
想看下泊松融合的效果，据说比多分辨率融合效果要好，因为多分辨率融合会使得融合区域变暗。论文上步骤写得不是特别清楚，看了好几篇都是这样，都写得差不多，介绍下泊松图像剪辑的原理公式那些，然后就开始泊松融合了，这样根本自己编不了 又上网查了下，好像没人写过，但是这个资料貌似还好http://blog.csdn.net/archielau/article/details/8765298这个更好http://blog.csdn.net/baimafujinji/article/details/46787837还讲了泊松融合的物理模型来源
抽象出数学模型用于图像融合中 不过这个讲的其实只能成为图像无缝镶嵌 并不是真正的将泊松用于图像拼接融合中 这个可能还是得自己编写 不过试了下这个泊松图像编辑用于图像镶嵌的程序：
和
将第一幅图无缝镶嵌到第二幅图中 直接用http://blog.csdn.net/baimafujinji/article/details/46787837他的程序 结果：
%图片来自《全景图像拼接关键技术研究》 但这篇论文只讲了用泊松融合来做图片嵌入 其实不是图像拼接
%资料来自于http://blog.csdn.net/baimafujinji/article/details/46787837
%他其实也只用泊松融合来做图片嵌入 并不是完整的图像拼接的程序 这个要自己改写
A=imread('F:\fisheye\bosong1.jpg');
B=imread('F:\fisheye\bosong2.jpg');
[ma,na,ka]=size(A);
[mb,nb,kb]=size(B);
gradient_inner=B(1:mb-2,2:nb-1,:)+B(3:mb,2:nb-1,:)+B(2:mb-1,1:nb-2,:)+B(2:mb-1,3:nb,:)-4*B(2:mb-1,2:nb-1,:);
Lap=[0,1,0;1,-4,1;0,1,0];
I1=conv2(double(B(:,:,1)),double(Lap));
I2=conv2(double(B(:,:,2)),double(Lap));
I3=conv2(double(B(:,:,3)),double(Lap));
gradient_inner(:,:,1)=I1(3:mb,3:nb);
gradient_inner(:,:,2)=I2(3:mb,3:nb);
gradient_inner(:,:,3)=I3(3:mb,3:nb);
dstX=70;
dstY=20;
rebuilt=A(dstY:dstY+mb-1,dstX:dstX+nb-1,:);
for n=1:1000
rebuilt(2:2:mb-1,2:2:nb-1,:)=(rebuilt(1:2:mb-2,2:2:nb-1,:)+rebuilt(3:2:mb,2:2:nb-1,:)+rebuilt(2:2:mb-1,1:2:nb-2,:)+rebuilt(2:2:mb-1,3:2:nb,:)-gradient_inner(1:2:mb-2,1:2:nb-2,:))/4;
rebuilt(3:2:mb-1,3:2:nb-1,:)=(rebuilt(2:2:mb-2,3:2:nb-1,:)+rebuilt(4:2:mb,3:2:nb-1,:)+rebuilt(3:2:mb-1,2:2:nb-2,:)+rebuilt(3:2:mb-1,4:2:nb,:)-gradient_inner(2:2:mb-2,2:2:nb-2,:))/4;
rebuilt(3:2:mb-1,2:2:nb-1,:)=(rebuilt(2:2:mb-2,2:2:nb-1,:)+rebuilt(4:2:mb,2:2:nb-1,:)+rebuilt(3:2:mb-1,1:2:nb-2,:)+rebuilt(3:2:mb-1,3:2:nb,:)-gradient_inner(2:2:mb-2,1:2:nb-2,:))/4;
rebuilt(2:2:mb-1,3:2:nb-1,:)=(rebuilt(1:2:mb-2,3:2:nb-1,:)+rebuilt(3:2:mb,3:2:nb-1,:)+rebuilt(2:2:mb-1,2:2:nb-2,:)+rebuilt(2:2:mb-1,4:2:nb,:)-gradient_inner(1:2:mb-2,2:2:nb-2,:))/4;
end
A(dstY:mb+dstY-1,dstX:nb+dstX-1,:)=rebuilt;
A=uint8(A);
%上面是直接按照博客上的程序 必须迭代插值1000次吗 出来结果不对！镶嵌进去的是一个渐变的晕染的黑月亮一样！
这个不对 不知道是他程序本来就不对还是怎样 我是和他基本一样啊 出来就是不对 得不到他的左图那样
这个是直接镶嵌：
结果这样 就是要让这个边界没有 要自然
按照原理 改成这样：
A=imread('F:\fisheye\bosong1.jpg');
B=imread('F:\fisheye\bosong2.jpg');
[ma,na,ka]=size(A);
[mb,nb,kb]=size(B);
gradient_inner=B(1:mb-2,2:nb-1,:)+B(3:mb,2:nb-1,:)+B(2:mb-1,1:nb-2,:)+B(2:mb-1,3:nb,:)-4*B(2:mb-1,2:nb-1,:);
Lap=[0,1,0;1,-4,1;0,1,0];
I1=conv2(double(B(:,:,1)),double(Lap));
I2=conv2(double(B(:,:,2)),double(Lap));
I3=conv2(double(B(:,:,3)),double(Lap));
gradient_inner(:,:,1)=I1(3:mb,3:nb);
gradient_inner(:,:,2)=I2(3:mb,3:nb);
gradient_inner(:,:,3)=I3(3:mb,3:nb);
dstX=70;
dstY=20;
for i=dstY:dstY+mb-3
for j=dstX:dstX+nb-3
ii=i-(dstY-1);
jj=j-(dstX-1);
A(i,j,1)=(A(i+1,j,1)+A(i,j+1,1)+A(i-1,j,1)+A(i,j-1,1)-gradient_inner(ii,jj,1))/4;
A(i,j,2)=(A(i+1,j,2)+A(i,j+1,2)+A(i-1,j,2)+A(i,j-1,2)-gradient_inner(ii,jj,2))/4;
A(i,j,3)=(A(i+1,j,3)+A(i,j+1,3)+A(i-1,j,3)+A(i,j-1,3)-gradient_inner(ii,jj,3))/4;
end
end
%上面的还是不对 有边界的黑框
这还是不对 是个有边界的不晕染的黑月亮框框！
上面不对 我找到http://blog.csdn.net/hjimce/article/details/45716603这个人写得比较好 原理步骤都讲得很清晰 浅显易懂 所以按照他的步骤一步一步来：
1，得到掩膜 并裁剪出源图中的包含目标的最小矩形 先直接覆盖在A图中：
%http://blog.csdn.net/hjimce/article/details/45716603按照这个资料
%先得到掩膜mask 方便裁剪出一个包含镶嵌图的最小矩形
B=imread('F:\fisheye\bosong2.jpg');
[mb,nb,kb]=size(B);
bmask=im2bw(B);
for i=1:mb
for j=1:nb
if(bmask(i,j)==1)
bottom=i;
end
end
end
for j=1:nb
for i=1:mb
if(bmask(i,j)==1)
right=j;
end
end
end
for i=mb:-1:1
for j=1:nb
if(bmask(i,j)==1)
up=i;
end
end
end
for j=nb:-1:1
for i=1:mb
if(bmask(i,j)==1)
left=j;
end
end
end
wid=-left;
height=bottom-up;
BB=uint8(zeros(height,wid,3));
for i=up:bottom
for j=left:right
ii=i-(up-1);
jj=j-left+1;
BB(ii,jj,1)=B(i,j,1);
BB(ii,jj,2)=B(i,j,2);
BB(ii,jj,3)=B(i,j,3);
end
end
[mbb,nbb,kbb]=size(BB);
dstX=70;
dstY=20;
A=imread('F:\fisheye\bosong1.jpg');
[ma,na,ka]=size(A);
for i=dstY:dstY+mbb-1
for j=dstX:dstX+nbb-1
ii=i-(dstY-1);
jj=j-(dstX-1);
A(i,j,1)=BB(ii,jj,1);
A(i,j,2)=BB(ii,jj,2);
A(i,j,3)=BB(ii,jj,3);
end
end
结果：
掩膜计算得有点不准确 不过这个没太大关系 反正主要关注泊松融合 这个我是直接用二值图判断得到的 误差就随它去了
2，计算BB的x、y方向的梯度场 计算A的x、y方向的梯度场 并将A的要放进月亮的位置的梯度场替换成BB的梯度场得到AA
%计算BB和AA的x方向y方向的梯度场
w=[0,-1,1];
BBx=imfilter(double(BB),w,'conv');
BBy=imfilter(double(BB),w','conv');
%imshow(BBx,[])
AAx=imfilter(double(AA),w,'conv');
AAy=imfilter(double(AA),w','conv');
%imshow(AAx,[])
%对AAx和AAy分别做mask 点乘 即等同于将AA待镶嵌区域的梯度场替换成BB的梯度场
for i=dstY:dstY+mbb-1
for j=dstX:dstX+nbb-1
ii=i-dstY+1;
jj=j-dstX+1;
AAx(i,j,1)=BBx(ii,jj,1);
AAx(i,j,2)=BBx(ii,jj,2);
AAx(i,j,3)=BBx(ii,jj,3);
AAy(i,j,1)=BBy(ii,jj,1);
AAy(i,j,2)=BBy(ii,jj,2);
AAy(i,j,3)=BBy(ii,jj,3);
end
end
结果：
3，计算AA的偏导得到散度 （拉普拉斯坐标）
%接下来对梯度求偏导得到融合图像的散度 lap
kernel=[-1,1,0];
lapx=imfilter(AAx,kernel,'conv');
lapy=imfilter(AAy,kernel','conv');
lap=lapx+lapy;
4，构建系数矩阵 ax=b 其中a是系数矩阵 x是所求的重构的像素点的各个通道值 b是散度
%接下来求解系数矩阵
%先将lap四周的点填充为原来的像素值（约束值） 其它的仍然是散度
[m,n,k]=size(lap);
for i=1:m
for j=1:n
if((i==1&&j==1)||(i==1&&j==n)||(i==m&&j==1)||(i==m&&j==n))
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
else
if(i==1||i==m)
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
else
if(j==1||j==n)
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
end
end
end
end
end
%然后将lap'排成一列 即Ax=b的b矩阵
for i=1:m
for j=1:n
lap1(i,j)=lap(i,j,1);
lap2(i,j)=lap(i,j,2);
lap3(i,j)=lap(i,j,3);
end
end
lapp1=lap1';
lapp2=lap2';
lapp3=lap3';
b1=lapp1(:);
b2=lapp2(:);
b3=lapp3(:);
b=[b1 b2 b3];
%计算系数矩阵Ax=b的A 这里用AAA 先对角线全部填1 然后非周围点就再例外赋值 1 1 -4 1 1
[mm,nn]=size(b);
AAA=int(zeros(mm,mm)); %There is error!!!!!
for i=1:mm
for j=1:mm
if(i==j)
AAA(i,j)=1;
end
end
end
但是到这步就有问题了 因为AAA是25071x25071大小 说超过了MATLAB的最大矩阵范围！！！如果这个问题不能解决 后面都不用做了 更不用说通过Ax=b求出x了 那求不出x 那融合后的每点的值就没办法得到了。 那个博主只说了在OpenCV里怎么泊松融合，因为在我出错的这步，OpenCV里有一种快速解泊松方程的方法，所以OpenCV里不是这样构造A和b来解x的。
然后我重新查了资料，原来pudn上有现成的别人传的泊松融合的程序matlab写的 我下下来 分析了下 运行了下 结果
左图是融合效果 这样看 好像的确是将右图的月亮融合进去了 无缝融合 但我换了幅图试下 结果：
其实下载的程序感官上就像是将左边的飞机处理成右边的飞机的效果 然后放进目标图中 这个融合效果不好 因为它带了背景 其实 这个程序是这样的：
src_img = imread('moon.jpg');
dst_img = imread('mountains.jpg'); %先读取两幅原图 一个是待嵌入图 一个是目标图
img_width = 95; %这个其实就是包含目标图的矩形 包含月亮的矩形 这个是上传文档的作者自己定的
img_height = 90;
src_st_x = 165; %这个就是原图中的点 表示待会儿从哪里开始切出这个矩形
src_st_y = 55;
dst_st_x = 700; %这个就是目标图中的点 表示待会儿从哪里开始放进这个包含月亮的矩形
dst_st_y = 350;
extracted_src_img = src_img(src_st_y:src_st_y+img_height,src_st_x:src_st_x+img_width,:); %截出月亮大小 以及原图中从某点开始的对应月亮区域
extracted_dst_img = dst_img(dst_st_y:dst_st_y+img_height,dst_st_x:dst_st_x+img_width,:);
[rws,cls] = size(extracted_src_img(:,:,1));
dst_img_temp = dst_img;
dst_img_temp(dst_st_y:dst_st_y+img_height,dst_st_x:dst_st_x+img_width,:) = extracted_src_img; %把月亮放进目标图中你之前想放的对应位置
padding_factor = 10;
padded_src_img = pad_image(extracted_src_img,padding_factor); %这个函数是月亮和原图那个小区域扩充边界 上下左右都多填充了10个像素的大小区域 用黑色填充的
padded_dst_img = pad_image(extracted_dst_img,padding_factor); %对切出的矩形和目标图的矩形都这样填充
temp_result = zeros(rws,cls,3);
temp_result(:,:,1) = find_result_channel(padded_src_img(:,:,1),padded_dst_img(:,:,1),padding_factor);
temp_result(:,:,2) = find_result_channel(padded_src_img(:,:,2),padded_dst_img(:,:,2),padding_factor);
temp_result(:,:,3) = find_result_channel(padded_src_img(:,:,3),padded_dst_img(:,:,3),padding_factor);
figure(1); imshow(uint8(temp_result)); %对每个通道计算系数矩阵A和b 求新像素值x 让月亮的颜色和要融合图的一样 ！
dst_img(dst_st_y:dst_st_y+img_height,dst_st_x:dst_st_x+img_width,:) = temp_result; %将处理好的月亮的矩形放进目标图中那个要换掉的矩形的位置
figure(2); imshow(uint8(dst_img));
其中几个子函数是这样的：
function [padded_image] = pad_image(src_img,padding_factor) %其实就是将src_img上下左右都填充padding_factor个像素大小 用黑色填充
rws = size(src_img,1);
cls = size(src_img,2);
temp = zeros(rws,padding_factor);
src_img_temp_1 = zeros(rws,cls + (2*padding_factor),3);
src_img_temp_1(:,:,1) = [temp src_img(:,:,1) temp];
src_img_temp_1(:,:,2) = [temp src_img(:,:,2) temp];
src_img_temp_1(:,:,3) = [temp src_img(:,:,3) temp];
src_img_temp_2 = zeros(rws+ (2*padding_factor),cls+ (2*padding_factor),3);
rws= size(src_img_temp_2,1);
cls = size(src_img_temp_2,2);
temp = zeros(padding_factor,cls);
src_img_temp_2(:,:,1) = [temp; src_img_temp_1(:,:,1); temp];
src_img_temp_2(:,:,2) = [temp; src_img_temp_1(:,:,2); temp];
src_img_temp_2(:,:,3) = [temp; src_img_temp_1(:,:,3); temp];
padded_image = src_img_temp_2;
function result = find_result_channel(src_img,dst_img,padding_factor) %这个其实就是将src_img按照dst_img进行处理 让它们看起来颜色亮度一致 或者说和谐
[rws cls] = size(src_img);
src_img_x = find_x_gradient(src_img,cls);
src_img_y = find_y_gradient(src_img,rws);
dst_img_x = find_x_gradient(dst_img,cls);
dst_img_y = find_y_gradient(dst_img,rws);
new_dst_img_x=dst_img_x;
new_dst_img_y=dst_img_y;
new_dst_img_x(2:rws-1,2:cls-1) = src_img_x(2:rws-1,2:cls-1);
new_dst_img_y(2:rws-1,2:cls-1) = src_img_y(2:rws-1,2:cls-1);
new_dst_img_x_x = find_x_gradient_b_diff(new_dst_img_x,cls);
new_dst_img_y_y = find_y_gradient_b_diff(new_dst_img_y,rws);
divergence_field = new_dst_img_x_x + new_dst_img_y_y;
divergence_field = divergence_field(padding_factor:size(divergence_field,1)-padding_factor-1,padding_factor:size(divergence_field,2)-padding_factor-1);
src_img = src_img(padding_factor+1:size(src_img,1)-padding_factor,padding_factor+1:size(src_img,2)-padding_factor);
dst_img = dst_img(padding_factor+1:size(dst_img,1)-padding_factor,padding_factor+1:size(dst_img,2)-padding_factor);
[rws cls] = size(src_img);
%laplace_sparse_matrix = Create_Laplace_Sparse_Matrix(rws, cls);
laplace_sparse_matrix = sp_laplace(rws, cls);
mask = zeros(rws,cls);
mask(:,1) = 1;
mask(:,cls) = 1;
mask(1,:) = 1;
mask(rws,:) = 1;
mask = logical(mask);
mask = mask(:);
divergence_field = divergence_field(:);
src_img = src_img(:);
dst_img = dst_img(:);
x_dash_dash = dst_img(mask);
A_dash_dash = laplace_sparse_matrix(:,mask);
b_dash_dash = A_dash_dash * x_dash_dash;
b_dash = divergence_field - b_dash_dash;
A_dash = laplace_sparse_matrix(~mask,~mask);
b_dash = b_dash(~mask,:);
x_dash = A_dash\b_dash;
result = zeros(size(divergence_field,1),1);
result(~mask,1) = x_dash;
result(mask,1) = x_dash_dash;
result = reshape(result,rws,cls);
这个函数还有几个子函数：
function [A,N] = sp_laplace(nx,ny) %这个函数是构建系数矩阵
N = nx * ny;
N5 = 5 * N;
irow = zeros(N5, 1);
icol = irow;
NZA = irow;
index = 0;
row = 0;
for j = 1:ny
for k = 1:nx
row = row + 1;
if j > 1
index = index + 1;
NZA (index) = -1.0;
irow(index) = row;
icol(index) = row - nx; % South
end
if k > 1
index = index + 1;
NZA (index) = -1.0;
irow(index) = row;
icol(index) = row - 1; % West
end
index = index + 1;
NZA (index) = 4.0;
irow(index) = row;
icol(index) = row; % P (self)
if k < nx
index = index + 1;
NZA (index) = -1.0;
irow(index) = row;
icol(index) = row + 1; % East
end
if j < ny
index = index + 1;
NZA (index) = -1.0;
irow(index) = row;
icol(index) = row + nx; % North
end
end
end
icol = icol(1:index);
irow = irow(1:index);
NZA = NZA(1:index);
A = sparse (irow, icol, NZA, N, N);
function [Iy] = find_y_gradient_b_diff(I,Ny)
I = double(I);
Iy = I - [I(2:Ny,:); I(1,:)];
function [Iy] = find_y_gradient(I,Ny)
Iy = [I(2:Ny,:); I(1,:)] - I;
function [Ix] = find_x_gradient_b_diff(I,Nx)
I = double(I);
Ix = I - [I(:,2:Nx) I(:,1)];
function [Ix] = find_x_gradient(I,Nx)
Ix =[I(:,2:Nx) I(:,1)]-I; %这是x方向的梯度？这样求梯度 第一列放到最后去 减去原图？
这就是pudn上别人传的用MATLAB写的泊松融合 但我感觉这个不对 我还是觉得应该像http://blog.csdn.net/hjimce/article/details/45716603这个人里的一样：
而不是像上面一样将包含飞机的最小矩形切出来 对比目标图的背景进行处理 让它们好像颜色亮度等和谐 然后放进去 我感觉不是这样！
找到一个更完整的资料http://blog.csdn.net/vsooda/article/details/38823745原来OpenCV3里面已经把泊松融合封装在了cloning这个类里面了！我的版本是2.4.9，所以我就稍微改了下：
#include <stdlib.h>
#include <opencv2/highgui/highgui.hpp>
#include<opencv2/opencv.hpp>
using namespace cv;
void dst(double *mod_diff, double *sineTransform, int h, int w)
{
unsigned long int idx;
Mat temp = Mat(2 * h + 2, 1, CV_32F);
Mat res = Mat(h, 1, CV_32F);
Mat planes[] = { Mat_<float>(temp), Mat::zeros(temp.size(), CV_32F) };
Mat result;
int p = 0;
for (int i = 0; i < w; i++)
{
temp.at<float>(0, 0) = 0.0;
for (int j = 0, r = 1; j < h; j++, r++)
{
idx = j*w + i;
temp.at<float>(r, 0) = (float)mod_diff[idx];
}
temp.at<float>(h + 1, 0) = 0.0;
for (int j = h - 1, r = h + 2; j >= 0; j--, r++)
{
idx = j*w + i;
temp.at<float>(r, 0) = (float)(-1.0 * mod_diff[idx]);
}
merge(planes, 2, result);
dft(result, result, 0, 0);
Mat planes1[] = { Mat::zeros(result.size(), CV_32F), Mat::zeros(result.size(), CV_32F) };
split(result, planes1);
std::complex<double> two_i = std::sqrt(std::complex<double>(-1));
double factor = -2 * imag(two_i);
for (int c = 1, z = 0; c < h + 1; c++, z++)
{
res.at<float>(z, 0) = (float)(planes1[1].at<float>(c, 0) / factor);
}
for (int q = 0, z = 0; q < h; q++, z++)
{
idx = q*w + p;
sineTransform[idx] = res.at<float>(z, 0);
}
p++;
}
}
void idst(double *mod_diff, double *sineTransform, int h, int w)
{
int nn = h + 1;
unsigned long int idx;
dst(mod_diff, sineTransform, h, w);
for (int i = 0; i < h; i++)
for (int j = 0; j < w; j++)
{
idx = i*w + j;
sineTransform[idx] = (double)(2 * sineTransform[idx]) / nn;
}
}
void transpose(double *mat, double *mat_t, int h, int w)
{
Mat tmp = Mat(h, w, CV_32FC1);
unsigned long int idx;
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*(w)+j;
tmp.at<float>(i, j) = (float)mat[idx];
}
}
Mat tmp_t = tmp.t();
for (int i = 0; i < tmp_t.size().height; i++)
for (int j = 0; j < tmp_t.size().width; j++)
{
idx = i*tmp_t.size().width + j;
mat_t[idx] = tmp_t.at<float>(i, j);
}
}
void solve(const Mat &img, double *mod_diff, Mat &result)
{
int w = img.size().width;
int h = img.size().height;
unsigned long int idx, idx1;
double *sineTransform = new double[(h - 2)*(w - 2)];
double *sineTransform_t = new double[(h - 2)*(w - 2)];
double *denom = new double[(h - 2)*(w - 2)];
double *invsineTransform = new double[(h - 2)*(w - 2)];
double *invsineTransform_t = new double[(h - 2)*(w - 2)];
double *img_d = new double[(h)*(w)];
//结果存在img_d
dst(mod_diff, sineTransform, h - 2, w - 2);
transpose(sineTransform, sineTransform_t, h - 2, w - 2);
dst(sineTransform_t, sineTransform, w - 2, h - 2);
transpose(sineTransform, sineTransform_t, w - 2, h - 2);
int cy = 1;
for (int i = 0; i < w - 2; i++, cy++)
{
for (int j = 0, cx = 1; j < h - 2; j++, cx++)
{
idx = j*(w - 2) + i;
denom[idx] = (float)2 * cos(CV_PI*cy / ((double)(w - 1))) - 2 + 2 * cos(CV_PI*cx / ((double)(h - 1))) - 2;
}
}
for (idx = 0; idx < (unsigned)(w - 2)*(h - 2); idx++)
{
sineTransform_t[idx] = sineTransform_t[idx] / denom[idx];
}
idst(sineTransform_t, invsineTransform, h - 2, w - 2);
transpose(invsineTransform, invsineTransform_t, h - 2, w - 2);
idst(invsineTransform_t, invsineTransform, w - 2, h - 2);
transpose(invsineTransform, invsineTransform_t, w - 2, h - 2);
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*w + j;
img_d[idx] = (double)img.at<uchar>(i, j);
}
}
for (int i = 1; i < h - 1; i++)
{
for (int j = 1; j < w - 1; j++)
{
idx = i*w + j;
img_d[idx] = 0.0;
}
}
for (int i = 1, id1 = 0; i < h - 1; i++, id1++)
{
for (int j = 1, id2 = 0; j < w - 1; j++, id2++)
{
idx = i*w + j;
idx1 = id1*(w - 2) + id2;
img_d[idx] = invsineTransform_t[idx1];
}
}
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*w + j;
if (img_d[idx] < 0.0) {
result.at<uchar>(i, j) = 0;
}
else if (img_d[idx] > 255.0)
result.at<uchar>(i, j) = 255;
else {
result.at<uchar>(i, j) = (uchar)img_d[idx];
}
}
}
delete[] sineTransform;
delete[] sineTransform_t;
delete[] denom;
delete[] invsineTransform;
delete[] invsineTransform_t;
delete[] img_d;
}
void poisson_solver(const Mat &img, Mat &gxx, Mat &gyy, Mat &result)
{
int w = img.size().width;
int h = img.size().height;
unsigned long int idx;
Mat lap = Mat(img.size(), CV_32FC1);
lap = gxx + gyy;
Mat bound = img.clone();
//rectangle 外围保持原样，rect内部变为scalar
rectangle(bound, Point(1, 1), Point(img.cols - 2, img.rows - 2), Scalar::all(0), -1);
//rectangle(bound, Point(20, 20), Point(img.cols - 50, img.rows - 50), Scalar::all(0), -1);
double *boundary_point = new double[h*w];
Mat bound_map = cv::Mat(img.size(), CV_8UC1, cv::Scalar(0));
for (int i = 1; i < h - 1; i++)
for (int j = 1; j < w - 1; j++)
{
idx = i*w + j;
boundary_point[idx] = -4 * (int)bound.at<uchar>(i, j) + (int)bound.at<uchar>(i, (j + 1)) + (int)bound.at<uchar>(i, (j - 1))+ (int)bound.at<uchar>(i - 1, j) + (int)bound.at<uchar>(i + 1, j);
bound_map.at<uchar>(i, j) = boundary_point[idx];
}
Mat diff = Mat(h, w, CV_32FC1);
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*w + j;
diff.at<float>(i, j) = (float)(lap.at<float>(i, j) - boundary_point[idx]);
}
}
//diff和lap几乎没什么区别
double *mod_diff = new double[(h - 2)*(w - 2)];
for (int i = 0; i < h - 2; i++)
{
for (int j = 0; j < w - 2; j++)
{
idx = i*(w - 2) + j;
mod_diff[idx] = diff.at<float>(i + 1, j + 1);
}
}
///////////////////////////////////////////////////// Find DST /////////////////////////////////////////////////////
solve(img, mod_diff, result);
delete[] mod_diff;
delete[] boundary_point;
}
void lapx(const Mat &img, Mat &gxx)
{
Mat kernel = Mat::zeros(1, 3, CV_8S);
kernel.at<char>(0, 0) = -1;
kernel.at<char>(0, 1) = 1;
filter2D(img, gxx, CV_32F, kernel);
}
void lapy(const Mat &img, Mat &gyy)
{
Mat kernel = Mat::zeros(3, 1, CV_8S);
kernel.at<char>(0, 0) = -1;
kernel.at<char>(1, 0) = 1;
filter2D(img, gyy, CV_32F, kernel);
}
void getGradientx(const Mat &img, Mat &gx)
{
Mat kernel = Mat::zeros(1, 3, CV_8S);
kernel.at<char>(0, 2) = 1;
kernel.at<char>(0, 1) = -1;
filter2D(img, gx, CV_32F, kernel);
}
void getGradienty(const Mat &img, Mat &gy)
{
Mat kernel = Mat::zeros(3, 1, CV_8S);
kernel.at<char>(2, 0) = 1;
kernel.at<char>(1, 0) = -1;
filter2D(img, gy, CV_32F, kernel);
}
void main()
{
Mat g=imread("bosong.jpg");
Mat I=imread("bosong1.jpg");
Mat gx,gy,sx,sy;
getGradientx(g,gx);
getGradienty(g,gy);
getGradientx(I,sx);
getGradienty(I,sy);
//下面其实就是poisson()函数
Mat fx = Mat(I.size(), CV_32FC3);
Mat fy = Mat(I.size(), CV_32FC3);
fx = gx + sx;
fy = gy + sy;
Mat gxx = Mat(I.size(), CV_32FC3);
Mat gyy = Mat(I.size(), CV_32FC3);
//gxx, gyy 是在x，y方向的laplacian算子
lapx(fx, gxx);
lapy(fy, gyy);
vector<Mat> rgb_channel,rgbx_channel,rgby_channel;
split(gxx, rgbx_channel);
split(gyy, rgby_channel);
split(I, rgb_channel);
Mat rgbchan2=rgb_channel.at(2);
Mat rgbchan1=rgb_channel.at(1);
Mat rgbchan0=rgb_channel.at(0);
Mat rgbxchan2=rgbx_channel.at(2);
Mat rgbychan2=rgby_channel.at(2);
Mat rgbxchan1=rgbx_channel.at(1);
Mat rgbychan1=rgby_channel.at(1);
Mat rgbxchan0=rgbx_channel.at(0);
Mat rgbychan0=rgby_channel.at(0);
Mat output1,output2,output0;
//运行下面这个函数时出错
poisson_solver(rgbchan2,rgbxchan2,rgbychan2,output2);
poisson_solver(rgbchan1,rgbxchan1,rgbychan1,output1);
poisson_solver(rgbchan0,rgbxchan0,rgbychan0,output0);
vector<Mat> output;
output.at(2)=output2;
output.at(1)=output1;
output.at(0)=output0;
Mat result;
merge(output,result);
namedWindow("poissonresult");
imshow("poissonresult",output1);
waitKey(0);
}
结果出错，在运行poisson_solver()时出错：
那个月亮的图我是扩大到和另一张一样大，用0填充的，所以大小一样了，这个错误，，，，，？后来检查了下，这个错误是在solve()最后一个for循环嵌套里，所以百度了下，这个错误很可能是去访问了矩阵中根本不存在的位置(超出了那个矩阵的大小) ，所以我加上了一句if(i>=0 && i<result.cols && j>=0 && j<result.rows)，就是把那两个for循环里的语句放在这个if{}内部，这样就没这个错误了。可是我想显示融合图每个通道的图像时候：
我这里只显示一个通道 可是
出错：显示说这个通道没有实际大小 ：
这个错误又是闹哪样。。。？
继续调试，结合http://blog.csdn.net/vsooda/article/details/38823745和http://blog.csdn.net/hjimce/article/details/45716603如果他们两个人的单独运行 有错误！试着调试结果还是不行 这个报错说我的结果blend是空的 我有试着下VS2015+OpenCV3.0因为想亲自看看那个类 但是本来装好了 结果我打开错误的标示以为是自己装错了 又卸载了 现在又重装！！！哎：
#include<opencv2/opencv.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<opencv2/imgproc/imgproc.hpp>
#include<iostream>
using namespace cv;
using namespace std;
void computeGradientX(const Mat &img,Mat &gx)
{
Mat kernel = Mat::zeros(1, 3, CV_8S);
kernel.at<char>(0,2) = 1;
kernel.at<char>(0,1) = -1;
if(img.channels() == 3)
filter2D(img, gx, CV_32F, kernel);
else if (img.channels() == 1)
{
Mat tmp[3];
for(int chan = 0 ; chan < 3 ; ++chan)
filter2D(img, tmp[chan], CV_32F, kernel);
merge(tmp, 3, gx);
}
}
void computeGradientY( const Mat &img, Mat &gy)
{
Mat kernel = Mat::zeros(3, 1, CV_8S);
kernel.at<char>(2,0) = 1;
kernel.at<char>(1,0) = -1;
if(img.channels() == 3)
filter2D(img, gy, CV_32F, kernel);
else if (img.channels() == 1)
{
Mat tmp[3];
for(int chan = 0 ; chan < 3 ; ++chan)
filter2D(img, tmp[chan], CV_32F, kernel);
merge(tmp, 3, gy);
}
}
void computeLaplacianX(const Mat &img,Mat &laplacianX)
{
Mat kernel = Mat::zeros(1, 3, CV_8S);
kernel.at<char>(0,0) = -1;
kernel.at<char>(0,1) = 1;
filter2D(img, laplacianX, CV_32F, kernel);
}
void computeLaplacianY(const Mat &img,Mat &laplacianY)
{
Mat kernel = Mat::zeros(3, 1, CV_8S);
kernel.at<char>(0,0) = -1;
kernel.at<char>(1,0) = 1;
filter2D(img, laplacianY, CV_32F, kernel);
}
void computeDerivatives(const Mat& destination, const Mat &patch, const Mat &binaryMask,Mat &binaryMaskFloat)
{
//initVariables(destination,binaryMask);//相关变量初始化，没用
Mat destinationGradientX,destinationGradientY,patchGradientX,patchGradientY;
computeGradientX(destination,destinationGradientX);//计算背景图像的x方向梯度
computeGradientY(destination,destinationGradientY);//计算背景图像y方向的梯度
computeGradientX(patch,patchGradientX);//计算ROI区域转换复制到destination一样大小的patch图片x方向梯度
computeGradientY(patch,patchGradientY);//计算y方向梯度
Mat Kernel(Size(3, 3),CV_8UC1);
Kernel.setTo(Scalar(1));
erode(binaryMask, binaryMask, Kernel, Point(-1,-1), 3); //对binaryMask进行腐蚀，去掉边界的毛刺点，让边界曲线平滑一点
binaryMask.convertTo(binaryMaskFloat,CV_32FC1,1.0/255.0);
}
//矩阵点乘，将lhs与rhs点乘得到result，因为有三个通道，估计mat不能实现三通道的矩阵的一次性点乘，所以才有这个函数
void arrayProduct(const cv::Mat& lhs, const cv::Mat& rhs, cv::Mat& result)
{
vector <Mat> lhs_channels;
vector <Mat> result_channels;
split(lhs,lhs_channels);//拆分成3个通道的矩阵
split(result,result_channels);
for(int chan = 0 ; chan < 3 ; ++chan)//三个矩阵进行分别相乘
multiply(lhs_channels[chan],rhs,result_channels[chan]);
merge(result_channels,result);//合成为一个
}
void dst(double *mod_diff, double *sineTransform, int h, int w)
{
unsigned long int idx;
Mat temp = Mat(2 * h + 2, 1, CV_32F);
Mat res = Mat(h, 1, CV_32F);
Mat planes[] = { Mat_<float>(temp), Mat::zeros(temp.size(), CV_32F) };
Mat result;
int p = 0;
for (int i = 0; i < w; i++)
{
temp.at<float>(0, 0) = 0.0;
for (int j = 0, r = 1; j < h; j++, r++)
{
idx = j*w + i;
temp.at<float>(r, 0)=(float)mod_diff[idx];
}
temp.at<float>(h + 1, 0) = 0.0;
for (int j = h - 1, r = h + 2; j >= 0; j--, r++)
{
idx = j*w + i;
temp.at<float>(r, 0) = (float)(-1.0 * mod_diff[idx]);
}
merge(planes, 2, result);
dft(result, result, 0, 0);
Mat planes1[] = { Mat::zeros(result.size(), CV_32F), Mat::zeros(result.size(), CV_32F) };
split(result, planes1);
std::complex<double> two_i = std::sqrt(std::complex<double>(-1));
double factor = -2 * imag(two_i);
for (int c = 1, z = 0; c < h + 1; c++, z++)
res.at<float>(z, 0) = (float)(planes1[1].at<float>(c, 0) / factor);
for (int q = 0, z = 0; q < h; q++, z++)
{
idx = q*w + p;
sineTransform[idx] = res.at<float>(z, 0);
}
p++;
}
}
void idst(double *mod_diff, double *sineTransform, int h, int w)
{
int nn = h + 1;
unsigned long int idx;
dst(mod_diff, sineTransform, h, w);
for (int i = 0; i < h; i++)
for (int j = 0; j < w; j++)
{
idx = i*w + j;
sineTransform[idx] = (double)(2 * sineTransform[idx]) / nn;
}
}
void transpose(double *mat, double *mat_t, int h, int w)
{
Mat tmp = Mat(h, w, CV_32FC1);
unsigned long int idx;
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*(w)+j;
tmp.at<float>(i, j) = (float)mat[idx];
}
}
Mat tmp_t = tmp.t();
for (int i = 0; i < tmp_t.size().height; i++)
for (int j = 0; j < tmp_t.size().width; j++)
{
idx = i*tmp_t.size().width + j;
mat_t[idx] = tmp_t.at<float>(i, j);
}
}
void solve(const Mat &myimg,double *mymod_diff, Mat &myresult)
{
int w = myimg.size().width;
int h = myimg.size().height;
unsigned long int idx, idx1;
double *sineTransform = new double[(h - 2)*(w - 2)];
double *sineTransform_t = new double[(h - 2)*(w - 2)];
double *denom = new double[(h - 2)*(w - 2)];
double *invsineTransform = new double[(h - 2)*(w - 2)];
double *invsineTransform_t = new double[(h - 2)*(w - 2)];
double *img_d = new double[(h)*(w)];
//结果存在img_d
dst(mymod_diff,sineTransform,h-2,w-2);
transpose(sineTransform,sineTransform_t,h-2,w-2);
dst(sineTransform_t,sineTransform,w-2,h-2);
transpose(sineTransform,sineTransform_t,w-2,h-2);
int cy = 1;
for (int i=0;i<w-2;i++,cy++)
{
for (int j=0,cx=1;j<h-2;j++,cx++)
{
idx = j*(w - 2) + i;
denom[idx] = (float)2 * cos(CV_PI*cy / ((double)(w - 1))) - 2 + 2 * cos(CV_PI*cx / ((double)(h - 1))) - 2;
}
}
for (idx = 0; idx < (unsigned)(w - 2)*(h - 2); idx++)
sineTransform_t[idx] = sineTransform_t[idx] / denom[idx];
idst(sineTransform_t, invsineTransform, h - 2, w - 2);
transpose(invsineTransform, invsineTransform_t, h - 2, w - 2);
idst(invsineTransform_t, invsineTransform, w - 2, h - 2);
transpose(invsineTransform, invsineTransform_t, w - 2, h - 2);
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*w + j;
img_d[idx] = (double)myimg.at<uchar>(i, j);
}
}
for (int i = 1; i < h - 1; i++)
{
for (int j = 1; j < w - 1; j++)
{
idx = i*w + j;
img_d[idx] = 0.0;
}
}
for (int i = 1, id1 = 0; i < h - 1; i++, id1++)
{
for (int j = 1, id2 = 0; j < w - 1; j++, id2++)
{
idx = i*w + j;
idx1 = id1*(w - 2) + id2;
img_d[idx] = invsineTransform_t[idx1];
}
}
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
if(i>=0 && i<myresult.cols && j>=0 && j<myresult.rows)
{
idx = i*w + j;
if (img_d[idx] < 0.0)
myresult.at<uchar>(i, j) = 0;
else if (img_d[idx] > 255.0)
myresult.at<uchar>(i, j) = 255;
else
myresult.at<uchar>(i, j) = (uchar)img_d[idx];
}
}
}
delete[] sineTransform;
delete[] sineTransform_t;
delete[] denom;
delete[] invsineTransform;
delete[] invsineTransform_t;
delete[] img_d;
}
//泊松方程求解 输入散度(laplacianX+laplacianY)，及边界点像素即可重建求解
//输入：img为背景图片，laplacianX+laplacianY 为散度
//输出：result重建结果
void poissonSolver(const Mat &img, Mat &laplacianX , Mat &laplacianY, Mat &result)
{
const int w = img.cols;
const int h = img.rows;
Mat lap = Mat(img.size(),CV_32FC1);
lap = laplacianX + laplacianY;//散度
Mat bound = img.clone();
//边界修正，opencv为了方便，直接把图片最外围的像素点排除在外，不参与泊松重建
rectangle(bound,Point(1,1),Point(img.cols-2,img.rows-2),Scalar::all(0),-1);
//Mat boundary_points;
//Laplacian(bound, boundary_points, CV_32F);
//boundary_points = lap - boundary_points;
//Mat mod_diff = boundary_points(Rect(1, 1, w-2, h-2));
//solve(img,mod_diff,result);
//把上面四句和下面的换掉
unsigned long int idx;
double *boundary_point = new double[h*w];
Mat bound_map = cv::Mat(img.size(), CV_8UC1, cv::Scalar(0));
for (int i = 1; i < h - 1; i++)
for (int j = 1; j < w - 1; j++)
{
idx = i*w + j;
boundary_point[idx] = -4 * (int)bound.at<uchar>(i, j) + (int)bound.at<uchar>(i, (j + 1)) + (int)bound.at<uchar>(i, (j - 1))
+ (int)bound.at<uchar>(i - 1, j) + (int)bound.at<uchar>(i + 1, j);
bound_map.at<uchar>(i, j) = boundary_point[idx];
}
Mat diff = Mat(h, w, CV_32FC1);
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*w + j;
diff.at<float>(i, j) = (float)(lap.at<float>(i, j) - boundary_point[idx]);
}
}
double *mod_diff = new double[(h - 2)*(w - 2)];
for (int i = 0; i < h - 2; i++)
{
for (int j = 0; j < w - 2; j++)
{
idx = i*(w - 2) + j;
mod_diff[idx] = diff.at<float>(i + 1, j + 1);
}
}
solve(img, mod_diff, result);
delete[] mod_diff;
delete[] boundary_point;
}
//泊松重建
void poisson(const Mat &destination,const Mat &destinationGradientX,const Mat &destinationGradientY,const Mat &patchGradientX,const Mat &patchGradientY,Mat &outputB,Mat &outputG,Mat &outputR)
{
Mat laplacianX = Mat(destination.size(),CV_32FC3);
Mat laplacianY = Mat(destination.size(),CV_32FC3);
//因为前面已经对destinationGradientX做了固定区域的mask，patchGradientX做了修改区域的mask
laplacianX = destinationGradientX + patchGradientX;//求解整张图片新的梯度场
laplacianY = destinationGradientY + patchGradientY;
computeLaplacianX(laplacianX,laplacianX);//求解梯度的散度 也就是拉普拉坐标
computeLaplacianY(laplacianY,laplacianY);
vector<Mat> rgbx_channel,rgby_channel;
split(laplacianX,rgbx_channel);//通道拆分
split(laplacianY,rgby_channel);
vector<Mat> output1;
split(destination,output1);
//for(int chan =0;chan<3;chan++)
// poissonSolver(output1[chan], rgbx_channel[chan], rgby_channel[chan], output[chan]); //这里有问题说mat里面不许这样？！我改用下面的 大不了多写几句
Mat output1B,output1G,output1R,rgbx_channelB,rgbx_channelG,rgbx_channelR,rgby_channelB,rgby_channelG,rgby_channelR;
output1B=output1.at(0);
output1G=output1.at(1);
output1R=output1.at(2);
rgbx_channelB=rgbx_channel.at(0);
rgbx_channelG=rgbx_channel.at(1);
rgbx_channelR=rgbx_channel.at(2);
rgby_channelB=rgby_channel.at(0);
rgby_channelG=rgby_channel.at(1);
rgby_channelR=rgby_channel.at(2);
poissonSolver(output1B, rgbx_channelB, rgby_channelB, outputB);
poissonSolver(output1G, rgbx_channelG, rgby_channelG, outputG);
poissonSolver(output1R, rgbx_channelR, rgby_channelR, outputR);
//output.at(0)=B; //我这样对吗 平常做右值的可以作为左值吗 ？
//output.at(1)=G;
//output.at(2)=R;
}
//输入：I背景整张图片 wmask背景中待修改的区域的mask
//输出：cloned
void evaluate(const Mat &I,const Mat &wmask,const Mat &patchGradientX,const Mat &patchGradientY, Mat &cloned)
{
bitwise_not(wmask,wmask);//矩阵元素取反操作，这样背景图片保持不变的像素对应的mask值为1
Mat binaryMaskFloatInverted;
wmask.convertTo(binaryMaskFloatInverted,CV_32FC1,1.0/255.0);
//上面已经对patchGradientX做了mask操作 ，这边也对destinationGradientX做mask操作
Mat mydestinationGradientX,mydestinationGradientY;
computeGradientX(I,mydestinationGradientX);//计算背景图像的x方向梯度
computeGradientY(I,mydestinationGradientY);//计算背景图像y方向的梯度
arrayProduct(mydestinationGradientX,binaryMaskFloatInverted, mydestinationGradientX);
arrayProduct(mydestinationGradientY,binaryMaskFloatInverted, mydestinationGradientY);
Mat outputB,outputG,outputR; //改了这里！！！！！！！！！！！！！！！！
poisson(I,mydestinationGradientX,mydestinationGradientY,patchGradientX,patchGradientY,outputB,outputG,outputR);
Mat output[]={outputB,outputG,outputR};
merge(output,3,cloned);
}
//克隆融合外部接口函数
//输入：destination背景图片的整张图片 binaryMask为destination待修改的像素的mask
//patch是由src图片的ROI区域复制过来的图像，其大小与destination相同，patch只有binaryMask区域存的是src的ROI图片
//输出：cloned融合结果整张图片
void normalClone(const Mat &destination, const Mat &patch, const Mat &binaryMask, Mat &cloned)
{
const int w = destination.cols;
const int h = destination.rows;
const int channel = destination.channels();
const int n_elem_in_line = w * channel;
//计算destination在x，y方向的梯度，获得结果为：destinationGradientX destinationGradientY
//计算patch在x，y方向的梯度，获得结果为： patchGradientX patchGradientY
//同时对binaryMask进行边界腐蚀，去除毛刺，让边界变得光滑一点，同时把binaryMask归一化为0~1得binaryMaskFloat
//其实这样归一化后binaryMaskFloat只有数值0,1 1代表即将被修改的像素点
Mat mybinaryMaskFloat,mypatchGradientX,mypatchGradientY;
computeDerivatives(destination,patch,binaryMask,mybinaryMaskFloat);
computeGradientX(patch,mypatchGradientX);//计算ROI区域转换复制到destination一样大小的patch图片x方向梯度
computeGradientY(patch,mypatchGradientY);//计算y方向梯度
//因为patch是一个src包含ROI的最小矩形块图片
//patchGradientX与binaryMaskFloat相乘，这样patchGradientX就只剩下有用的区域了
arrayProduct(mypatchGradientX,mybinaryMaskFloat,mypatchGradientX);
arrayProduct(mypatchGradientY,mybinaryMaskFloat,mypatchGradientY);
evaluate(destination,binaryMask,mypatchGradientX,mypatchGradientY,cloned);
}
//无缝融合，这个函数其实没什么功能，只是为了减少、方便计算，先把要融合的区域裁剪下来，也就是对_src、_dst、_mask 进行裁剪到最小
//输入：_src前景图像 _dst背景图像 _mask前景图像的mask，
//p是用于对应用的，p点指的是在dst中，待融合区域中心点(dst的待融合区域的大小是根据包含熊的矩形的大小确定的)
//输出：_blend融合结果图片
//void seamlessClone(InputArray _src, InputArray _dst, InputArray _mask,Mat &mybinaryMaskFloat, Point p, OutputArray _blend, int flags)
void seamlessClone(InputArray _src, InputArray _dst, InputArray _mask, Point p,OutputArray _blend)
{
const Mat src = _src.getMat();
const Mat dest = _dst.getMat();
const Mat mask = _mask.getMat();
_blend.create(dest.size(), CV_8UC3);//融合后图片的大小肯定跟背景图像一样
Mat blend = _blend.getMat();
int minx = INT_MAX, miny = INT_MAX, maxx = INT_MIN, maxy = INT_MIN;
int h = mask.size().height;
int w = mask.size().width;
Mat gray = Mat(mask.size(),CV_8UC1);
Mat dst_mask = Mat::zeros(dest.size(),CV_8UC1);//背景图像的mask
Mat cs_mask = Mat::zeros(src.size(),CV_8UC3);
Mat cd_mask = Mat::zeros(dest.size(),CV_8UC3);
if(mask.channels() == 3)//如果给定的mask是彩色图 需要转换成单通道灰度图
cvtColor(mask, gray, COLOR_BGR2GRAY );
else
gray = mask;
//计算包含mask的最小矩形，也就是把那只熊包含起来的最小矩形框，这个矩形是位于src的，后面还有一个对应的矩形位于dst 这两个for循环有错误？！！！
for(int i=0;i<h;i++)
{
for(int j=0;j<w;j++)
{
if(gray.at<uchar>(i,j) == 255)
{
minx = std::min(minx,i);
maxx = std::max(maxx,i);
miny = std::min(miny,j);
maxy = std::max(maxy,j);
}
}
}
int lenx = maxx - minx;//计算矩形的宽
int leny = maxy - miny;//计算矩形的高
Mat patch = Mat::zeros(Size(leny, lenx), CV_8UC3);//根据上面的矩形区域，创建一个大小相同矩阵
int minxd = p.y - lenx/2;//计算dst的矩形
int maxxd = p.y + lenx/2;
int minyd = p.x - leny/2;
int maxyd = p.x + leny/2;
CV_Assert(minxd>=0&&minyd>=0&&maxxd<=dest.rows&&maxyd<=dest.cols);
Rect roi_d(minyd,minxd,leny,lenx);//dst 兴趣区域的矩形
Rect roi_s(miny,minx,leny,lenx);//src 兴趣区域矩形
Mat destinationROI = dst_mask(roi_d);
Mat sourceROI = cs_mask(roi_s);
gray(roi_s).copyTo(destinationROI);//
src(roi_s).copyTo(sourceROI,gray(roi_s));
src(roi_s).copyTo(patch, gray(roi_s));//patch为
destinationROI = cd_mask(roi_d);
cs_mask(roi_s).copyTo(destinationROI);//cs_mask为把前景图片的矩形区域图像 转换到背景图片矩形中的图片
//Cloning obj;
//obj.normalClone(dest,cd_mask,dst_mask,blend,flags);
normalClone(dest,cd_mask,dst_mask,blend);
}
void main()
{
//第一步 读取两幅原图 并选定目标图中要融合的的区域
Mat src_primer_img=imread("airplane.jpg");
Mat dst_primer_img=imread("mountains.jpg");
if(!src_primer_img.data||!dst_primer_img.data)
cout<<"damn!wrong in reading!"<<endl;
Mat dst_img=dst_primer_img(Rect(100,300,src_primer_img.cols,src_primer_img.rows));
//imshow("dst",dst_img);
//第二步 对source即待嵌入的图进行腐蚀得到src_img_mask（这里得到的mask结果仍是彩图）
Mat src_img,src_img_mask;
Mat Kernel(Size(3, 3),CV_8UC1);
Kernel.setTo(Scalar(1));
erode(src_primer_img,src_img, Kernel, Point(-1,-1), 3); //对binaryMask进行腐蚀，去掉边界的毛刺点，让边界曲线平滑一点
src_img.convertTo(src_img_mask,CV_32FC1,1.0/255.0);
//computeDerivatives(dst_img,src_img,src_img,src_img_mask);
//imshow("srcmask",src_img_mask);
//imwrite("mask.jpg",src_img_mask);
//第三步
Point p=Point(100,300); //原始大图 要自己定点 中心点 矩形的！！！
Mat blend;
//seamlessClone(src_img,dst_primer_img,src_img_mask,p,blend); 替换成下面的 方便找出从哪里开始错误
Mat src;
src_img.copyTo(src);
Mat dest;
dst_primer_img.copyTo(dest);
Mat mask;
//src_img_mask.copyTo(mask);
mask=Mat::zeros(src_img_mask.size(),CV_32FC1);
//_blend.create(dest.size(), CV_8UC3);//融合后图片的大小肯定跟背景图像一样
dest.copyTo(blend);
int minx = INT_MAX, miny = INT_MAX, maxx = INT_MIN, maxy = INT_MIN;
int h = mask.size().height;
int w = mask.size().width;
Mat gray = Mat(mask.size(),CV_8UC1);
Mat dst_mask = Mat::zeros(dest.size(),CV_8UC1);//背景图像的mask
Mat cs_mask = Mat::zeros(src.size(),CV_8UC3);
Mat cd_mask = Mat::zeros(dest.size(),CV_8UC3);
if(mask.channels() == 3)//如果给定的mask是彩色图 需要转换成单通道灰度图
cvtColor(mask,gray,COLOR_BGR2GRAY );
else
gray = mask;
//imwrite("gray.jpg",gray); //保存显示检查错误 为什么这里显示的是腐蚀后的灰度图 但保存后的图片打开却是个全黑屏呢？奇怪
//imshow("gray",gray);
//计算包含mask的最小矩形，也就是把那只熊包含起来的最小矩形框，这个矩形是位于src的，后面还有一个对应的矩形位于dst 这两个for循环开始出现错误 貌似不是访问像素点越界的问题啊？？if(i>=0&&i<gray.rows&&j>=0&&j<gray.cols)
for(int i=0;i<h;++i)
{
for(int j=0;j<w;++j)
{
if((int)gray.at<float>(i,j)==255)
{
minx = std::min(minx,i);
maxx = std::max(maxx,i);
miny = std::min(miny,j);
maxy = std::max(maxy,j);
}
}
}
int lenx = maxx - minx;//计算矩形的宽
int leny = maxy - miny;//计算矩形的高
Mat patch = Mat::zeros(Size(leny, lenx), CV_8UC3);//根据上面的矩形区域，创建一个大小相同矩阵
int minxd = p.y - lenx/2;//计算dst的矩形
int maxxd = p.y + lenx/2;
int minyd = p.x - leny/2;
int maxyd = p.x + leny/2;
CV_Assert(minxd>=0&&minyd>=0&&maxxd<=dest.rows&&maxyd<=dest.cols);
Rect roi_d(minyd,minxd,leny,lenx);//dst 兴趣区域的矩形
Rect roi_s(miny,minx,leny,lenx);//src 兴趣区域矩形
Mat destinationROI = dst_mask(roi_d);
Mat sourceROI = cs_mask(roi_s);
gray(roi_s).copyTo(destinationROI);
//src(roi_s).copyTo(sourceROI,gray(roi_s)); //这句开始出现错误！！！！！！！！！！！！！！！！！！改成下面这样
src(roi_s).copyTo(sourceROI);
//src(roi_s).copyTo(patch, gray(roi_s));//patch为
src(roi_s).copyTo(patch);
destinationROI = cd_mask(roi_d);
cs_mask(roi_s).copyTo(destinationROI);//cs_mask为把前景图片的矩形区域图像 转换到背景图片矩形中的图片
normalClone(dest,cd_mask,dst_mask,blend);
namedWindow("poissonresult",1);
imshow("poissonresult",blend);
waitKey(0);
}
但是中间结果即待镶嵌图的mask还是挺好的 比我在matlab下自己之前写的得到的要好：
重来
全部重来写的 继续用之前写的MATLAB版本的 就是之前构造系数矩阵太大没写了的 现在接着写 那个原来可以用传说中的稀疏矩阵解决 我以前只听过没真的看过稀疏矩阵 也没去了解过现在知道了：
这是两个原图 我主要是用泊松融合让那只飞机可以无缝融合进右图中 把左图得到的mask放进右图
然后我用我构造系数矩阵然后重新解方程组得到的结果是：
A=imread('F:\fisheye\mountains.jpg');
src=A;
[ma,na,ka]=size(A);
B=imread('F:\fisheye\airplane.jpg');
dst=B;
[mb,nb,kb]=size(B);
se=strel('diamond',10);
B_erode=imerode(B,se);
Berodelogical=im2bw(B_erode);
for i=1:mb
for j=1:nb
if(Berodelogical(i,j)==1)
Berodelogical(i,j)=0;
else
Berodelogical(i,j)=1;
end
end
end
dstX=200;
dstY=100;
for i=dstY:dstY+mb-1
for j=dstX:dstX+nb-1
ii=i-(dstY-1);
jj=j-(dstX-1);
if(Berodelogical(ii,jj)==1)
A(i,j,1)=B(ii,jj,1);
A(i,j,2)=B(ii,jj,2);
A(i,j,3)=B(ii,jj,3);
end
end
end
%上面是把飞机放进目标图里 然后对目标图的新ROI计算梯度 散度
ROI=uint8(zeros(mb,nb,3));
for i=dstY:dstY+mb-1
for j=dstX:dstX+nb-1
ii=i-(dstY-1);
jj=j-(dstX-1);
ROI(ii,jj,1)=A(i,j,1);
ROI(ii,jj,2)=A(i,j,2);
ROI(ii,jj,3)=A(i,j,3);
end
end
w=[0,-1,1];
ROIgradienty=imfilter(double(ROI),w,'conv');
ROIgradientx=imfilter(double(ROI),w','conv');
%接下来对梯度求偏导得到融合图像的散度 lap
kernel=[-1,1,0];
lapy=imfilter(ROIgradienty,kernel,'conv');
lapx=imfilter(ROIgradientx,kernel','conv');
lap=lapx+lapy;
%接下来求解系数矩阵
%先将lap四周的点填充为原来的像素值（约束值） 其它的仍然是散度
[m,n,k]=size(lap);
heightRegion=m;
widthRegion=n;
L=n;
J=m;
for i=1:m
for j=1:n
if((i==1&&j==1)||(i==1&&j==n)||(i==m&&j==1)||(i==m&&j==n))
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
else
if(i==1||i==m)
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
else
if(j==1||j==n)
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
end
end
end
end
end
%然后将lap'排成一列 即Ax=b的b矩阵
for i=1:m
for j=1:n
lap1(i,j)=lap(i,j,1);
lap2(i,j)=lap(i,j,2);
lap3(i,j)=lap(i,j,3);
end
end
lapp1=lap1';
lapp2=lap2';
lapp3=lap3';
b1=lapp1(:);
b2=lapp2(:);
b3=lapp3(:);
b=[b1 b2 b3];
%计算系数矩阵Ax=b的A 这里用M
[mm,nn]=size(b);
M=spalloc(mm,mm,6*(mm-2*m-2*n+4));
for i=1:mm
for j=1:mm
if(i==j)
M(i,j)=1;
end
end
end
%上面是边界点的约束值仍然是它本身的像素值 就是乘以1它本身就好
for y=1:heightRegion
for x=1:widthRegion
if x~=1&&y~=1&&y~=heightRegion&&x~=widthRegion
i=(y-1)*L+x;
M(i,i)=-4;
M(i,i-1)=1;
M(i,i+1)=1;
M(i,i-L)=1;
M(i,i+L)=1;
end
end
end
% use bi-conjugate gradient method to solve the matrix
%x = bicg(M,b',[],400); %稀疏矩阵的专业解法迭代求解 传说中的双共轭梯度法 不太了解 不能用平常的X=A\b’好像
%x=bicg(M,b,1e-6,400); %每一行是一个点的RGB三值 理论上的 但实际不能这样 换成下面
x1=bicg(M,b1,1e-6,400);
x2=bicg(M,b2,1e-6,400);
x3=bicg(M,b3,1e-6,400);
% reshape x to become the pixel intensity of the region
imRegion1= reshape(x1,widthRegion,heightRegion);
imRegion2= reshape(x2,widthRegion,heightRegion);
imRegion3= reshape(x3,widthRegion,heightRegion);
imNew=src;
imNew(dstY:dstY+heightRegion-1,dstX:dstX+widthRegion-1,1)=imRegion1';
imNew(dstY:dstY+heightRegion-1,dstX:dstX+widthRegion-1,2)=imRegion2';
imNew(dstY:dstY+heightRegion-1,dstX:dstX+widthRegion-1,3)=imRegion3';
这个效果不太对！！！我不放飞机进去 放月亮
右图是对月亮做mask后直接相加的结果 左图是后面泊松重建后的效果 忽略两图月亮的毛刺 因为我腐蚀时边界没处理 感觉泊松融合这次还是不够正确啊 可是我检查了下还是按照http://blog.csdn.net/hjimce/article/details/45716603这个人的理论来的啊 后面求解泊松方程时 用了稀疏矩阵 周围的点约束值是它们本身的像素值
其它点就是散度 系数矩阵也对了啊 怎么还是达不到它的效果呢
左图是直接做mask镶嵌的结果 右图是解了泊松方程后的最终的结果 我所能得到的：
感觉效果不够好 因为这个理想结果是这样：
看这个就比我的自然很多很多 我的应该是有问题的 但我是按照步骤来的 目前不知道错在哪里 在找：[今天经过网友的提醒，找到了，原来是迭代次数问题，我这个效果是只迭代一次的，再迭代几次就能达到理想效果，那篇论文中没有提过这个迭代！害死我了]
A=imread('F:\fisheye\pond.jpg');
src=A;
[ma,na,ka]=size(A);
dst=imread('F:\fisheye\swim.jpg');
se=strel('diamond',10);
B_erode=imerode(dst,se);
Berodelogical=im2bw(B_erode);
imshow(Berodelogical)
B=dst;
[mb,nb,kb]=size(B);
>> for i=1:mb
for j=1:nb
if(Berodelogical(i,j)==1)
Berodelogical(i,j)=0;
else
Berodelogical(i,j)=1;
end
end
end
dstX=100;
dstY=100;
for i=dstY:dstY+mb-1
for j=dstX:dstX+nb-1
ii=i-(dstY-1);
jj=j-(dstX-1);
if(Berodelogical(ii,jj)==1)
A(i,j,1)=B(ii,jj,1);
A(i,j,2)=B(ii,jj,2);
A(i,j,3)=B(ii,jj,3);
end
end
end
ROI=uint8(zeros(mb,nb,3));
for i=dstY:dstY+mb-1
for j=dstX:dstX+nb-1
ii=i-(dstY-1);
jj=j-(dstX-1);
ROI(ii,jj,1)=A(i,j,1);
ROI(ii,jj,2)=A(i,j,2);
ROI(ii,jj,3)=A(i,j,3);
end
end
w=[0,-1,1];
ROIgradienty=imfilter(double(ROI),w,'conv');
ROIgradientx=imfilter(double(ROI),w','conv');
%接下来对梯度求偏导得到融合图像的散度 lap
kernel=[-1,1,0];
lapy=imfilter(ROIgradienty,kernel,'conv');
lapx=imfilter(ROIgradientx,kernel','conv');
lap=lapx+lapy;
%接下来求解系数矩阵
%先将lap四周的点填充为原来的像素值（约束值） 其它的仍然是散度
[m,n,k]=size(lap);
heightRegion=m;
widthRegion=n;
L=n;
J=m;
for i=1:m
for j=1:n
if((i==1&&j==1)||(i==1&&j==n)||(i==m&&j==1)||(i==m&&j==n))
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
else
if(i==1||i==m)
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
else
if(j==1||j==n)
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
end
end
end
end
end
%然后将lap'排成一列 即Ax=b的b矩阵
for i=1:m
for j=1:n
lap1(i,j)=lap(i,j,1);
lap2(i,j)=lap(i,j,2);
lap3(i,j)=lap(i,j,3);
end
end
lapp1=lap1';
lapp2=lap2';
lapp3=lap3';
b1=lapp1(:);
b2=lapp2(:);
b3=lapp3(:);
b=[b1 b2 b3];
%计算系数矩阵Ax=b的A 这里用M
[mm,nn]=size(b);
M=spalloc(mm,mm,6*(mm-2*m-2*n+4));
for i=1:mm
for j=1:mm
if(i==j)
M(i,j)=1;
end
end
end
%上面是边界点的约束值仍然是它本身的像素值 就是乘以1它本身就好
for y=1:heightRegion
for x=1:widthRegion
if x~=1&&y~=1&&y~=heightRegion&&x~=widthRegion
i=(y-1)*L+x;
M(i,i)=-4;
M(i,i-1)=1;
M(i,i+1)=1;
M(i,i-L)=1;
M(i,i+L)=1;
end
end
end
x1=bicg(M,b1,1e-6,400); %用传说中的稀疏矩阵 原来就是为了防止我之前的问题out of memory！当有很多0时的大矩阵
x2=bicg(M,b2,1e-6,400); %这个双共轭梯度解法 来解Mx=b的这个x 400是迭代次数
x3=bicg(M,b3,1e-6,400);
% reshape x to become the pixel intensity of the region
imRegion1= reshape(x1,widthRegion,heightRegion);
imRegion2= reshape(x2,widthRegion,heightRegion);
imRegion3= reshape(x3,widthRegion,heightRegion);
imNew=src;
imNew(dstY:dstY+heightRegion-1,dstX:dstX+widthRegion-1,1)=imRegion1';
imNew(dstY:dstY+heightRegion-1,dstX:dstX+widthRegion-1,2)=imRegion2';
imNew(dstY:dstY+heightRegion-1,dstX:dstX+widthRegion-1,3)=imRegion3';
下面这些是逐步探索的得到上面这个的过程：
想看下泊松融合的效果，据说比多分辨率融合效果要好，因为多分辨率融合会使得融合区域变暗。论文上步骤写得不是特别清楚，看了好几篇都是这样，都写得差不多，介绍下泊松图像剪辑的原理公式那些，然后就开始泊松融合了，这样根本自己编不了 又上网查了下，好像没人写过，但是这个资料貌似还好http://blog.csdn.net/archielau/article/details/8765298这个更好http://blog.csdn.net/baimafujinji/article/details/46787837还讲了泊松融合的物理模型来源
抽象出数学模型用于图像融合中 不过这个讲的其实只能成为图像无缝镶嵌 并不是真正的将泊松用于图像拼接融合中 这个可能还是得自己编写 不过试了下这个泊松图像编辑用于图像镶嵌的程序：
和
将第一幅图无缝镶嵌到第二幅图中 直接用http://blog.csdn.net/baimafujinji/article/details/46787837他的程序 结果：
%图片来自《全景图像拼接关键技术研究》 但这篇论文只讲了用泊松融合来做图片嵌入 其实不是图像拼接
%资料来自于http://blog.csdn.net/baimafujinji/article/details/46787837
%他其实也只用泊松融合来做图片嵌入 并不是完整的图像拼接的程序 这个要自己改写
A=imread('F:\fisheye\bosong1.jpg');
B=imread('F:\fisheye\bosong2.jpg');
[ma,na,ka]=size(A);
[mb,nb,kb]=size(B);
gradient_inner=B(1:mb-2,2:nb-1,:)+B(3:mb,2:nb-1,:)+B(2:mb-1,1:nb-2,:)+B(2:mb-1,3:nb,:)-4*B(2:mb-1,2:nb-1,:);
Lap=[0,1,0;1,-4,1;0,1,0];
I1=conv2(double(B(:,:,1)),double(Lap));
I2=conv2(double(B(:,:,2)),double(Lap));
I3=conv2(double(B(:,:,3)),double(Lap));
gradient_inner(:,:,1)=I1(3:mb,3:nb);
gradient_inner(:,:,2)=I2(3:mb,3:nb);
gradient_inner(:,:,3)=I3(3:mb,3:nb);
dstX=70;
dstY=20;
rebuilt=A(dstY:dstY+mb-1,dstX:dstX+nb-1,:);
for n=1:1000
rebuilt(2:2:mb-1,2:2:nb-1,:)=(rebuilt(1:2:mb-2,2:2:nb-1,:)+rebuilt(3:2:mb,2:2:nb-1,:)+rebuilt(2:2:mb-1,1:2:nb-2,:)+rebuilt(2:2:mb-1,3:2:nb,:)-gradient_inner(1:2:mb-2,1:2:nb-2,:))/4;
rebuilt(3:2:mb-1,3:2:nb-1,:)=(rebuilt(2:2:mb-2,3:2:nb-1,:)+rebuilt(4:2:mb,3:2:nb-1,:)+rebuilt(3:2:mb-1,2:2:nb-2,:)+rebuilt(3:2:mb-1,4:2:nb,:)-gradient_inner(2:2:mb-2,2:2:nb-2,:))/4;
rebuilt(3:2:mb-1,2:2:nb-1,:)=(rebuilt(2:2:mb-2,2:2:nb-1,:)+rebuilt(4:2:mb,2:2:nb-1,:)+rebuilt(3:2:mb-1,1:2:nb-2,:)+rebuilt(3:2:mb-1,3:2:nb,:)-gradient_inner(2:2:mb-2,1:2:nb-2,:))/4;
rebuilt(2:2:mb-1,3:2:nb-1,:)=(rebuilt(1:2:mb-2,3:2:nb-1,:)+rebuilt(3:2:mb,3:2:nb-1,:)+rebuilt(2:2:mb-1,2:2:nb-2,:)+rebuilt(2:2:mb-1,4:2:nb,:)-gradient_inner(1:2:mb-2,2:2:nb-2,:))/4;
end
A(dstY:mb+dstY-1,dstX:nb+dstX-1,:)=rebuilt;
A=uint8(A);
%上面是直接按照博客上的程序 必须迭代插值1000次吗 出来结果不对！镶嵌进去的是一个渐变的晕染的黑月亮一样！
这个不对 不知道是他程序本来就不对还是怎样 我是和他基本一样啊 出来就是不对 得不到他的左图那样
这个是直接镶嵌：
结果这样 就是要让这个边界没有 要自然
按照原理 改成这样：
A=imread('F:\fisheye\bosong1.jpg');
B=imread('F:\fisheye\bosong2.jpg');
[ma,na,ka]=size(A);
[mb,nb,kb]=size(B);
gradient_inner=B(1:mb-2,2:nb-1,:)+B(3:mb,2:nb-1,:)+B(2:mb-1,1:nb-2,:)+B(2:mb-1,3:nb,:)-4*B(2:mb-1,2:nb-1,:);
Lap=[0,1,0;1,-4,1;0,1,0];
I1=conv2(double(B(:,:,1)),double(Lap));
I2=conv2(double(B(:,:,2)),double(Lap));
I3=conv2(double(B(:,:,3)),double(Lap));
gradient_inner(:,:,1)=I1(3:mb,3:nb);
gradient_inner(:,:,2)=I2(3:mb,3:nb);
gradient_inner(:,:,3)=I3(3:mb,3:nb);
dstX=70;
dstY=20;
for i=dstY:dstY+mb-3
for j=dstX:dstX+nb-3
ii=i-(dstY-1);
jj=j-(dstX-1);
A(i,j,1)=(A(i+1,j,1)+A(i,j+1,1)+A(i-1,j,1)+A(i,j-1,1)-gradient_inner(ii,jj,1))/4;
A(i,j,2)=(A(i+1,j,2)+A(i,j+1,2)+A(i-1,j,2)+A(i,j-1,2)-gradient_inner(ii,jj,2))/4;
A(i,j,3)=(A(i+1,j,3)+A(i,j+1,3)+A(i-1,j,3)+A(i,j-1,3)-gradient_inner(ii,jj,3))/4;
end
end
%上面的还是不对 有边界的黑框
这还是不对 是个有边界的不晕染的黑月亮框框！
上面不对 我找到http://blog.csdn.net/hjimce/article/details/45716603这个人写得比较好 原理步骤都讲得很清晰 浅显易懂 所以按照他的步骤一步一步来：
1，得到掩膜 并裁剪出源图中的包含目标的最小矩形 先直接覆盖在A图中：
%http://blog.csdn.net/hjimce/article/details/45716603按照这个资料
%先得到掩膜mask 方便裁剪出一个包含镶嵌图的最小矩形
B=imread('F:\fisheye\bosong2.jpg');
[mb,nb,kb]=size(B);
bmask=im2bw(B);
for i=1:mb
for j=1:nb
if(bmask(i,j)==1)
bottom=i;
end
end
end
for j=1:nb
for i=1:mb
if(bmask(i,j)==1)
right=j;
end
end
end
for i=mb:-1:1
for j=1:nb
if(bmask(i,j)==1)
up=i;
end
end
end
for j=nb:-1:1
for i=1:mb
if(bmask(i,j)==1)
left=j;
end
end
end
wid=-left;
height=bottom-up;
BB=uint8(zeros(height,wid,3));
for i=up:bottom
for j=left:right
ii=i-(up-1);
jj=j-left+1;
BB(ii,jj,1)=B(i,j,1);
BB(ii,jj,2)=B(i,j,2);
BB(ii,jj,3)=B(i,j,3);
end
end
[mbb,nbb,kbb]=size(BB);
dstX=70;
dstY=20;
A=imread('F:\fisheye\bosong1.jpg');
[ma,na,ka]=size(A);
for i=dstY:dstY+mbb-1
for j=dstX:dstX+nbb-1
ii=i-(dstY-1);
jj=j-(dstX-1);
A(i,j,1)=BB(ii,jj,1);
A(i,j,2)=BB(ii,jj,2);
A(i,j,3)=BB(ii,jj,3);
end
end
结果：
掩膜计算得有点不准确 不过这个没太大关系 反正主要关注泊松融合 这个我是直接用二值图判断得到的 误差就随它去了
2，计算BB的x、y方向的梯度场 计算A的x、y方向的梯度场 并将A的要放进月亮的位置的梯度场替换成BB的梯度场得到AA
%计算BB和AA的x方向y方向的梯度场
w=[0,-1,1];
BBx=imfilter(double(BB),w,'conv');
BBy=imfilter(double(BB),w','conv');
%imshow(BBx,[])
AAx=imfilter(double(AA),w,'conv');
AAy=imfilter(double(AA),w','conv');
%imshow(AAx,[])
%对AAx和AAy分别做mask 点乘 即等同于将AA待镶嵌区域的梯度场替换成BB的梯度场
for i=dstY:dstY+mbb-1
for j=dstX:dstX+nbb-1
ii=i-dstY+1;
jj=j-dstX+1;
AAx(i,j,1)=BBx(ii,jj,1);
AAx(i,j,2)=BBx(ii,jj,2);
AAx(i,j,3)=BBx(ii,jj,3);
AAy(i,j,1)=BBy(ii,jj,1);
AAy(i,j,2)=BBy(ii,jj,2);
AAy(i,j,3)=BBy(ii,jj,3);
end
end
结果：
3，计算AA的偏导得到散度 （拉普拉斯坐标）
%接下来对梯度求偏导得到融合图像的散度 lap
kernel=[-1,1,0];
lapx=imfilter(AAx,kernel,'conv');
lapy=imfilter(AAy,kernel','conv');
lap=lapx+lapy;
4，构建系数矩阵 ax=b 其中a是系数矩阵 x是所求的重构的像素点的各个通道值 b是散度
%接下来求解系数矩阵
%先将lap四周的点填充为原来的像素值（约束值） 其它的仍然是散度
[m,n,k]=size(lap);
for i=1:m
for j=1:n
if((i==1&&j==1)||(i==1&&j==n)||(i==m&&j==1)||(i==m&&j==n))
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
else
if(i==1||i==m)
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
else
if(j==1||j==n)
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
end
end
end
end
end
%然后将lap'排成一列 即Ax=b的b矩阵
for i=1:m
for j=1:n
lap1(i,j)=lap(i,j,1);
lap2(i,j)=lap(i,j,2);
lap3(i,j)=lap(i,j,3);
end
end
lapp1=lap1';
lapp2=lap2';
lapp3=lap3';
b1=lapp1(:);
b2=lapp2(:);
b3=lapp3(:);
b=[b1 b2 b3];
%计算系数矩阵Ax=b的A 这里用AAA 先对角线全部填1 然后非周围点就再例外赋值 1 1 -4 1 1
[mm,nn]=size(b);
AAA=int(zeros(mm,mm)); %There is error!!!!!
for i=1:mm
for j=1:mm
if(i==j)
AAA(i,j)=1;
end
end
end
但是到这步就有问题了 因为AAA是25071x25071大小 说超过了MATLAB的最大矩阵范围！！！如果这个问题不能解决 后面都不用做了 更不用说通过Ax=b求出x了 那求不出x 那融合后的每点的值就没办法得到了。 那个博主只说了在OpenCV里怎么泊松融合，因为在我出错的这步，OpenCV里有一种快速解泊松方程的方法，所以OpenCV里不是这样构造A和b来解x的。
然后我重新查了资料，原来pudn上有现成的别人传的泊松融合的程序matlab写的 我下下来 分析了下 运行了下 结果
左图是融合效果 这样看 好像的确是将右图的月亮融合进去了 无缝融合 但我换了幅图试下 结果：
其实下载的程序感官上就像是将左边的飞机处理成右边的飞机的效果 然后放进目标图中 这个融合效果不好 因为它带了背景 其实 这个程序是这样的：
src_img = imread('moon.jpg');
dst_img = imread('mountains.jpg'); %先读取两幅原图 一个是待嵌入图 一个是目标图
img_width = 95; %这个其实就是包含目标图的矩形 包含月亮的矩形 这个是上传文档的作者自己定的
img_height = 90;
src_st_x = 165; %这个就是原图中的点 表示待会儿从哪里开始切出这个矩形
src_st_y = 55;
dst_st_x = 700; %这个就是目标图中的点 表示待会儿从哪里开始放进这个包含月亮的矩形
dst_st_y = 350;
extracted_src_img = src_img(src_st_y:src_st_y+img_height,src_st_x:src_st_x+img_width,:); %截出月亮大小 以及原图中从某点开始的对应月亮区域
extracted_dst_img = dst_img(dst_st_y:dst_st_y+img_height,dst_st_x:dst_st_x+img_width,:);
[rws,cls] = size(extracted_src_img(:,:,1));
dst_img_temp = dst_img;
dst_img_temp(dst_st_y:dst_st_y+img_height,dst_st_x:dst_st_x+img_width,:) = extracted_src_img; %把月亮放进目标图中你之前想放的对应位置
padding_factor = 10;
padded_src_img = pad_image(extracted_src_img,padding_factor); %这个函数是月亮和原图那个小区域扩充边界 上下左右都多填充了10个像素的大小区域 用黑色填充的
padded_dst_img = pad_image(extracted_dst_img,padding_factor); %对切出的矩形和目标图的矩形都这样填充
temp_result = zeros(rws,cls,3);
temp_result(:,:,1) = find_result_channel(padded_src_img(:,:,1),padded_dst_img(:,:,1),padding_factor);
temp_result(:,:,2) = find_result_channel(padded_src_img(:,:,2),padded_dst_img(:,:,2),padding_factor);
temp_result(:,:,3) = find_result_channel(padded_src_img(:,:,3),padded_dst_img(:,:,3),padding_factor);
figure(1); imshow(uint8(temp_result)); %对每个通道计算系数矩阵A和b 求新像素值x 让月亮的颜色和要融合图的一样 ！
dst_img(dst_st_y:dst_st_y+img_height,dst_st_x:dst_st_x+img_width,:) = temp_result; %将处理好的月亮的矩形放进目标图中那个要换掉的矩形的位置
figure(2); imshow(uint8(dst_img));
其中几个子函数是这样的：
function [padded_image] = pad_image(src_img,padding_factor) %其实就是将src_img上下左右都填充padding_factor个像素大小 用黑色填充
rws = size(src_img,1);
cls = size(src_img,2);
temp = zeros(rws,padding_factor);
src_img_temp_1 = zeros(rws,cls + (2*padding_factor),3);
src_img_temp_1(:,:,1) = [temp src_img(:,:,1) temp];
src_img_temp_1(:,:,2) = [temp src_img(:,:,2) temp];
src_img_temp_1(:,:,3) = [temp src_img(:,:,3) temp];
src_img_temp_2 = zeros(rws+ (2*padding_factor),cls+ (2*padding_factor),3);
rws= size(src_img_temp_2,1);
cls = size(src_img_temp_2,2);
temp = zeros(padding_factor,cls);
src_img_temp_2(:,:,1) = [temp; src_img_temp_1(:,:,1); temp];
src_img_temp_2(:,:,2) = [temp; src_img_temp_1(:,:,2); temp];
src_img_temp_2(:,:,3) = [temp; src_img_temp_1(:,:,3); temp];
padded_image = src_img_temp_2;
function result = find_result_channel(src_img,dst_img,padding_factor) %这个其实就是将src_img按照dst_img进行处理 让它们看起来颜色亮度一致 或者说和谐
[rws cls] = size(src_img);
src_img_x = find_x_gradient(src_img,cls);
src_img_y = find_y_gradient(src_img,rws);
dst_img_x = find_x_gradient(dst_img,cls);
dst_img_y = find_y_gradient(dst_img,rws);
new_dst_img_x=dst_img_x;
new_dst_img_y=dst_img_y;
new_dst_img_x(2:rws-1,2:cls-1) = src_img_x(2:rws-1,2:cls-1);
new_dst_img_y(2:rws-1,2:cls-1) = src_img_y(2:rws-1,2:cls-1);
new_dst_img_x_x = find_x_gradient_b_diff(new_dst_img_x,cls);
new_dst_img_y_y = find_y_gradient_b_diff(new_dst_img_y,rws);
divergence_field = new_dst_img_x_x + new_dst_img_y_y;
divergence_field = divergence_field(padding_factor:size(divergence_field,1)-padding_factor-1,padding_factor:size(divergence_field,2)-padding_factor-1);
src_img = src_img(padding_factor+1:size(src_img,1)-padding_factor,padding_factor+1:size(src_img,2)-padding_factor);
dst_img = dst_img(padding_factor+1:size(dst_img,1)-padding_factor,padding_factor+1:size(dst_img,2)-padding_factor);
[rws cls] = size(src_img);
%laplace_sparse_matrix = Create_Laplace_Sparse_Matrix(rws, cls);
laplace_sparse_matrix = sp_laplace(rws, cls);
mask = zeros(rws,cls);
mask(:,1) = 1;
mask(:,cls) = 1;
mask(1,:) = 1;
mask(rws,:) = 1;
mask = logical(mask);
mask = mask(:);
divergence_field = divergence_field(:);
src_img = src_img(:);
dst_img = dst_img(:);
x_dash_dash = dst_img(mask);
A_dash_dash = laplace_sparse_matrix(:,mask);
b_dash_dash = A_dash_dash * x_dash_dash;
b_dash = divergence_field - b_dash_dash;
A_dash = laplace_sparse_matrix(~mask,~mask);
b_dash = b_dash(~mask,:);
x_dash = A_dash\b_dash;
result = zeros(size(divergence_field,1),1);
result(~mask,1) = x_dash;
result(mask,1) = x_dash_dash;
result = reshape(result,rws,cls);
这个函数还有几个子函数：
function [A,N] = sp_laplace(nx,ny) %这个函数是构建系数矩阵
N = nx * ny;
N5 = 5 * N;
irow = zeros(N5, 1);
icol = irow;
NZA = irow;
index = 0;
row = 0;
for j = 1:ny
for k = 1:nx
row = row + 1;
if j > 1
index = index + 1;
NZA (index) = -1.0;
irow(index) = row;
icol(index) = row - nx; % South
end
if k > 1
index = index + 1;
NZA (index) = -1.0;
irow(index) = row;
icol(index) = row - 1; % West
end
index = index + 1;
NZA (index) = 4.0;
irow(index) = row;
icol(index) = row; % P (self)
if k < nx
index = index + 1;
NZA (index) = -1.0;
irow(index) = row;
icol(index) = row + 1; % East
end
if j < ny
index = index + 1;
NZA (index) = -1.0;
irow(index) = row;
icol(index) = row + nx; % North
end
end
end
icol = icol(1:index);
irow = irow(1:index);
NZA = NZA(1:index);
A = sparse (irow, icol, NZA, N, N);
function [Iy] = find_y_gradient_b_diff(I,Ny)
I = double(I);
Iy = I - [I(2:Ny,:); I(1,:)];
function [Iy] = find_y_gradient(I,Ny)
Iy = [I(2:Ny,:); I(1,:)] - I;
function [Ix] = find_x_gradient_b_diff(I,Nx)
I = double(I);
Ix = I - [I(:,2:Nx) I(:,1)];
function [Ix] = find_x_gradient(I,Nx)
Ix =[I(:,2:Nx) I(:,1)]-I; %这是x方向的梯度？这样求梯度 第一列放到最后去 减去原图？
这就是pudn上别人传的用MATLAB写的泊松融合 但我感觉这个不对 我还是觉得应该像http://blog.csdn.net/hjimce/article/details/45716603这个人里的一样：
而不是像上面一样将包含飞机的最小矩形切出来 对比目标图的背景进行处理 让它们好像颜色亮度等和谐 然后放进去 我感觉不是这样！
找到一个更完整的资料http://blog.csdn.net/vsooda/article/details/38823745原来OpenCV3里面已经把泊松融合封装在了cloning这个类里面了！我的版本是2.4.9，所以我就稍微改了下：
#include <stdlib.h>
#include <opencv2/highgui/highgui.hpp>
#include<opencv2/opencv.hpp>
using namespace cv;
void dst(double *mod_diff, double *sineTransform, int h, int w)
{
unsigned long int idx;
Mat temp = Mat(2 * h + 2, 1, CV_32F);
Mat res = Mat(h, 1, CV_32F);
Mat planes[] = { Mat_<float>(temp), Mat::zeros(temp.size(), CV_32F) };
Mat result;
int p = 0;
for (int i = 0; i < w; i++)
{
temp.at<float>(0, 0) = 0.0;
for (int j = 0, r = 1; j < h; j++, r++)
{
idx = j*w + i;
temp.at<float>(r, 0) = (float)mod_diff[idx];
}
temp.at<float>(h + 1, 0) = 0.0;
for (int j = h - 1, r = h + 2; j >= 0; j--, r++)
{
idx = j*w + i;
temp.at<float>(r, 0) = (float)(-1.0 * mod_diff[idx]);
}
merge(planes, 2, result);
dft(result, result, 0, 0);
Mat planes1[] = { Mat::zeros(result.size(), CV_32F), Mat::zeros(result.size(), CV_32F) };
split(result, planes1);
std::complex<double> two_i = std::sqrt(std::complex<double>(-1));
double factor = -2 * imag(two_i);
for (int c = 1, z = 0; c < h + 1; c++, z++)
{
res.at<float>(z, 0) = (float)(planes1[1].at<float>(c, 0) / factor);
}
for (int q = 0, z = 0; q < h; q++, z++)
{
idx = q*w + p;
sineTransform[idx] = res.at<float>(z, 0);
}
p++;
}
}
void idst(double *mod_diff, double *sineTransform, int h, int w)
{
int nn = h + 1;
unsigned long int idx;
dst(mod_diff, sineTransform, h, w);
for (int i = 0; i < h; i++)
for (int j = 0; j < w; j++)
{
idx = i*w + j;
sineTransform[idx] = (double)(2 * sineTransform[idx]) / nn;
}
}
void transpose(double *mat, double *mat_t, int h, int w)
{
Mat tmp = Mat(h, w, CV_32FC1);
unsigned long int idx;
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*(w)+j;
tmp.at<float>(i, j) = (float)mat[idx];
}
}
Mat tmp_t = tmp.t();
for (int i = 0; i < tmp_t.size().height; i++)
for (int j = 0; j < tmp_t.size().width; j++)
{
idx = i*tmp_t.size().width + j;
mat_t[idx] = tmp_t.at<float>(i, j);
}
}
void solve(const Mat &img, double *mod_diff, Mat &result)
{
int w = img.size().width;
int h = img.size().height;
unsigned long int idx, idx1;
double *sineTransform = new double[(h - 2)*(w - 2)];
double *sineTransform_t = new double[(h - 2)*(w - 2)];
double *denom = new double[(h - 2)*(w - 2)];
double *invsineTransform = new double[(h - 2)*(w - 2)];
double *invsineTransform_t = new double[(h - 2)*(w - 2)];
double *img_d = new double[(h)*(w)];
//结果存在img_d
dst(mod_diff, sineTransform, h - 2, w - 2);
transpose(sineTransform, sineTransform_t, h - 2, w - 2);
dst(sineTransform_t, sineTransform, w - 2, h - 2);
transpose(sineTransform, sineTransform_t, w - 2, h - 2);
int cy = 1;
for (int i = 0; i < w - 2; i++, cy++)
{
for (int j = 0, cx = 1; j < h - 2; j++, cx++)
{
idx = j*(w - 2) + i;
denom[idx] = (float)2 * cos(CV_PI*cy / ((double)(w - 1))) - 2 + 2 * cos(CV_PI*cx / ((double)(h - 1))) - 2;
}
}
for (idx = 0; idx < (unsigned)(w - 2)*(h - 2); idx++)
{
sineTransform_t[idx] = sineTransform_t[idx] / denom[idx];
}
idst(sineTransform_t, invsineTransform, h - 2, w - 2);
transpose(invsineTransform, invsineTransform_t, h - 2, w - 2);
idst(invsineTransform_t, invsineTransform, w - 2, h - 2);
transpose(invsineTransform, invsineTransform_t, w - 2, h - 2);
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*w + j;
img_d[idx] = (double)img.at<uchar>(i, j);
}
}
for (int i = 1; i < h - 1; i++)
{
for (int j = 1; j < w - 1; j++)
{
idx = i*w + j;
img_d[idx] = 0.0;
}
}
for (int i = 1, id1 = 0; i < h - 1; i++, id1++)
{
for (int j = 1, id2 = 0; j < w - 1; j++, id2++)
{
idx = i*w + j;
idx1 = id1*(w - 2) + id2;
img_d[idx] = invsineTransform_t[idx1];
}
}
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*w + j;
if (img_d[idx] < 0.0) {
result.at<uchar>(i, j) = 0;
}
else if (img_d[idx] > 255.0)
result.at<uchar>(i, j) = 255;
else {
result.at<uchar>(i, j) = (uchar)img_d[idx];
}
}
}
delete[] sineTransform;
delete[] sineTransform_t;
delete[] denom;
delete[] invsineTransform;
delete[] invsineTransform_t;
delete[] img_d;
}
void poisson_solver(const Mat &img, Mat &gxx, Mat &gyy, Mat &result)
{
int w = img.size().width;
int h = img.size().height;
unsigned long int idx;
Mat lap = Mat(img.size(), CV_32FC1);
lap = gxx + gyy;
Mat bound = img.clone();
//rectangle 外围保持原样，rect内部变为scalar
rectangle(bound, Point(1, 1), Point(img.cols - 2, img.rows - 2), Scalar::all(0), -1);
//rectangle(bound, Point(20, 20), Point(img.cols - 50, img.rows - 50), Scalar::all(0), -1);
double *boundary_point = new double[h*w];
Mat bound_map = cv::Mat(img.size(), CV_8UC1, cv::Scalar(0));
for (int i = 1; i < h - 1; i++)
for (int j = 1; j < w - 1; j++)
{
idx = i*w + j;
boundary_point[idx] = -4 * (int)bound.at<uchar>(i, j) + (int)bound.at<uchar>(i, (j + 1)) + (int)bound.at<uchar>(i, (j - 1))+ (int)bound.at<uchar>(i - 1, j) + (int)bound.at<uchar>(i + 1, j);
bound_map.at<uchar>(i, j) = boundary_point[idx];
}
Mat diff = Mat(h, w, CV_32FC1);
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*w + j;
diff.at<float>(i, j) = (float)(lap.at<float>(i, j) - boundary_point[idx]);
}
}
//diff和lap几乎没什么区别
double *mod_diff = new double[(h - 2)*(w - 2)];
for (int i = 0; i < h - 2; i++)
{
for (int j = 0; j < w - 2; j++)
{
idx = i*(w - 2) + j;
mod_diff[idx] = diff.at<float>(i + 1, j + 1);
}
}
///////////////////////////////////////////////////// Find DST /////////////////////////////////////////////////////
solve(img, mod_diff, result);
delete[] mod_diff;
delete[] boundary_point;
}
void lapx(const Mat &img, Mat &gxx)
{
Mat kernel = Mat::zeros(1, 3, CV_8S);
kernel.at<char>(0, 0) = -1;
kernel.at<char>(0, 1) = 1;
filter2D(img, gxx, CV_32F, kernel);
}
void lapy(const Mat &img, Mat &gyy)
{
Mat kernel = Mat::zeros(3, 1, CV_8S);
kernel.at<char>(0, 0) = -1;
kernel.at<char>(1, 0) = 1;
filter2D(img, gyy, CV_32F, kernel);
}
void getGradientx(const Mat &img, Mat &gx)
{
Mat kernel = Mat::zeros(1, 3, CV_8S);
kernel.at<char>(0, 2) = 1;
kernel.at<char>(0, 1) = -1;
filter2D(img, gx, CV_32F, kernel);
}
void getGradienty(const Mat &img, Mat &gy)
{
Mat kernel = Mat::zeros(3, 1, CV_8S);
kernel.at<char>(2, 0) = 1;
kernel.at<char>(1, 0) = -1;
filter2D(img, gy, CV_32F, kernel);
}
void main()
{
Mat g=imread("bosong.jpg");
Mat I=imread("bosong1.jpg");
Mat gx,gy,sx,sy;
getGradientx(g,gx);
getGradienty(g,gy);
getGradientx(I,sx);
getGradienty(I,sy);
//下面其实就是poisson()函数
Mat fx = Mat(I.size(), CV_32FC3);
Mat fy = Mat(I.size(), CV_32FC3);
fx = gx + sx;
fy = gy + sy;
Mat gxx = Mat(I.size(), CV_32FC3);
Mat gyy = Mat(I.size(), CV_32FC3);
//gxx, gyy 是在x，y方向的laplacian算子
lapx(fx, gxx);
lapy(fy, gyy);
vector<Mat> rgb_channel,rgbx_channel,rgby_channel;
split(gxx, rgbx_channel);
split(gyy, rgby_channel);
split(I, rgb_channel);
Mat rgbchan2=rgb_channel.at(2);
Mat rgbchan1=rgb_channel.at(1);
Mat rgbchan0=rgb_channel.at(0);
Mat rgbxchan2=rgbx_channel.at(2);
Mat rgbychan2=rgby_channel.at(2);
Mat rgbxchan1=rgbx_channel.at(1);
Mat rgbychan1=rgby_channel.at(1);
Mat rgbxchan0=rgbx_channel.at(0);
Mat rgbychan0=rgby_channel.at(0);
Mat output1,output2,output0;
//运行下面这个函数时出错
poisson_solver(rgbchan2,rgbxchan2,rgbychan2,output2);
poisson_solver(rgbchan1,rgbxchan1,rgbychan1,output1);
poisson_solver(rgbchan0,rgbxchan0,rgbychan0,output0);
vector<Mat> output;
output.at(2)=output2;
output.at(1)=output1;
output.at(0)=output0;
Mat result;
merge(output,result);
namedWindow("poissonresult");
imshow("poissonresult",output1);
waitKey(0);
}
结果出错，在运行poisson_solver()时出错：
那个月亮的图我是扩大到和另一张一样大，用0填充的，所以大小一样了，这个错误，，，，，？后来检查了下，这个错误是在solve()最后一个for循环嵌套里，所以百度了下，这个错误很可能是去访问了矩阵中根本不存在的位置(超出了那个矩阵的大小) ，所以我加上了一句if(i>=0 && i<result.cols && j>=0 && j<result.rows)，就是把那两个for循环里的语句放在这个if{}内部，这样就没这个错误了。可是我想显示融合图每个通道的图像时候：
我这里只显示一个通道 可是
出错：显示说这个通道没有实际大小 ：
这个错误又是闹哪样。。。？
继续调试，结合http://blog.csdn.net/vsooda/article/details/38823745和http://blog.csdn.net/hjimce/article/details/45716603如果他们两个人的单独运行 有错误！试着调试结果还是不行 这个报错说我的结果blend是空的 我有试着下VS2015+OpenCV3.0因为想亲自看看那个类 但是本来装好了 结果我打开错误的标示以为是自己装错了 又卸载了 现在又重装！！！哎：
#include<opencv2/opencv.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<opencv2/imgproc/imgproc.hpp>
#include<iostream>
using namespace cv;
using namespace std;
void computeGradientX(const Mat &img,Mat &gx)
{
Mat kernel = Mat::zeros(1, 3, CV_8S);
kernel.at<char>(0,2) = 1;
kernel.at<char>(0,1) = -1;
if(img.channels() == 3)
filter2D(img, gx, CV_32F, kernel);
else if (img.channels() == 1)
{
Mat tmp[3];
for(int chan = 0 ; chan < 3 ; ++chan)
filter2D(img, tmp[chan], CV_32F, kernel);
merge(tmp, 3, gx);
}
}
void computeGradientY( const Mat &img, Mat &gy)
{
Mat kernel = Mat::zeros(3, 1, CV_8S);
kernel.at<char>(2,0) = 1;
kernel.at<char>(1,0) = -1;
if(img.channels() == 3)
filter2D(img, gy, CV_32F, kernel);
else if (img.channels() == 1)
{
Mat tmp[3];
for(int chan = 0 ; chan < 3 ; ++chan)
filter2D(img, tmp[chan], CV_32F, kernel);
merge(tmp, 3, gy);
}
}
void computeLaplacianX(const Mat &img,Mat &laplacianX)
{
Mat kernel = Mat::zeros(1, 3, CV_8S);
kernel.at<char>(0,0) = -1;
kernel.at<char>(0,1) = 1;
filter2D(img, laplacianX, CV_32F, kernel);
}
void computeLaplacianY(const Mat &img,Mat &laplacianY)
{
Mat kernel = Mat::zeros(3, 1, CV_8S);
kernel.at<char>(0,0) = -1;
kernel.at<char>(1,0) = 1;
filter2D(img, laplacianY, CV_32F, kernel);
}
void computeDerivatives(const Mat& destination, const Mat &patch, const Mat &binaryMask,Mat &binaryMaskFloat)
{
//initVariables(destination,binaryMask);//相关变量初始化，没用
Mat destinationGradientX,destinationGradientY,patchGradientX,patchGradientY;
computeGradientX(destination,destinationGradientX);//计算背景图像的x方向梯度
computeGradientY(destination,destinationGradientY);//计算背景图像y方向的梯度
computeGradientX(patch,patchGradientX);//计算ROI区域转换复制到destination一样大小的patch图片x方向梯度
computeGradientY(patch,patchGradientY);//计算y方向梯度
Mat Kernel(Size(3, 3),CV_8UC1);
Kernel.setTo(Scalar(1));
erode(binaryMask, binaryMask, Kernel, Point(-1,-1), 3); //对binaryMask进行腐蚀，去掉边界的毛刺点，让边界曲线平滑一点
binaryMask.convertTo(binaryMaskFloat,CV_32FC1,1.0/255.0);
}
//矩阵点乘，将lhs与rhs点乘得到result，因为有三个通道，估计mat不能实现三通道的矩阵的一次性点乘，所以才有这个函数
void arrayProduct(const cv::Mat& lhs, const cv::Mat& rhs, cv::Mat& result)
{
vector <Mat> lhs_channels;
vector <Mat> result_channels;
split(lhs,lhs_channels);//拆分成3个通道的矩阵
split(result,result_channels);
for(int chan = 0 ; chan < 3 ; ++chan)//三个矩阵进行分别相乘
multiply(lhs_channels[chan],rhs,result_channels[chan]);
merge(result_channels,result);//合成为一个
}
void dst(double *mod_diff, double *sineTransform, int h, int w)
{
unsigned long int idx;
Mat temp = Mat(2 * h + 2, 1, CV_32F);
Mat res = Mat(h, 1, CV_32F);
Mat planes[] = { Mat_<float>(temp), Mat::zeros(temp.size(), CV_32F) };
Mat result;
int p = 0;
for (int i = 0; i < w; i++)
{
temp.at<float>(0, 0) = 0.0;
for (int j = 0, r = 1; j < h; j++, r++)
{
idx = j*w + i;
temp.at<float>(r, 0)=(float)mod_diff[idx];
}
temp.at<float>(h + 1, 0) = 0.0;
for (int j = h - 1, r = h + 2; j >= 0; j--, r++)
{
idx = j*w + i;
temp.at<float>(r, 0) = (float)(-1.0 * mod_diff[idx]);
}
merge(planes, 2, result);
dft(result, result, 0, 0);
Mat planes1[] = { Mat::zeros(result.size(), CV_32F), Mat::zeros(result.size(), CV_32F) };
split(result, planes1);
std::complex<double> two_i = std::sqrt(std::complex<double>(-1));
double factor = -2 * imag(two_i);
for (int c = 1, z = 0; c < h + 1; c++, z++)
res.at<float>(z, 0) = (float)(planes1[1].at<float>(c, 0) / factor);
for (int q = 0, z = 0; q < h; q++, z++)
{
idx = q*w + p;
sineTransform[idx] = res.at<float>(z, 0);
}
p++;
}
}
void idst(double *mod_diff, double *sineTransform, int h, int w)
{
int nn = h + 1;
unsigned long int idx;
dst(mod_diff, sineTransform, h, w);
for (int i = 0; i < h; i++)
for (int j = 0; j < w; j++)
{
idx = i*w + j;
sineTransform[idx] = (double)(2 * sineTransform[idx]) / nn;
}
}
void transpose(double *mat, double *mat_t, int h, int w)
{
Mat tmp = Mat(h, w, CV_32FC1);
unsigned long int idx;
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*(w)+j;
tmp.at<float>(i, j) = (float)mat[idx];
}
}
Mat tmp_t = tmp.t();
for (int i = 0; i < tmp_t.size().height; i++)
for (int j = 0; j < tmp_t.size().width; j++)
{
idx = i*tmp_t.size().width + j;
mat_t[idx] = tmp_t.at<float>(i, j);
}
}
void solve(const Mat &myimg,double *mymod_diff, Mat &myresult)
{
int w = myimg.size().width;
int h = myimg.size().height;
unsigned long int idx, idx1;
double *sineTransform = new double[(h - 2)*(w - 2)];
double *sineTransform_t = new double[(h - 2)*(w - 2)];
double *denom = new double[(h - 2)*(w - 2)];
double *invsineTransform = new double[(h - 2)*(w - 2)];
double *invsineTransform_t = new double[(h - 2)*(w - 2)];
double *img_d = new double[(h)*(w)];
//结果存在img_d
dst(mymod_diff,sineTransform,h-2,w-2);
transpose(sineTransform,sineTransform_t,h-2,w-2);
dst(sineTransform_t,sineTransform,w-2,h-2);
transpose(sineTransform,sineTransform_t,w-2,h-2);
int cy = 1;
for (int i=0;i<w-2;i++,cy++)
{
for (int j=0,cx=1;j<h-2;j++,cx++)
{
idx = j*(w - 2) + i;
denom[idx] = (float)2 * cos(CV_PI*cy / ((double)(w - 1))) - 2 + 2 * cos(CV_PI*cx / ((double)(h - 1))) - 2;
}
}
for (idx = 0; idx < (unsigned)(w - 2)*(h - 2); idx++)
sineTransform_t[idx] = sineTransform_t[idx] / denom[idx];
idst(sineTransform_t, invsineTransform, h - 2, w - 2);
transpose(invsineTransform, invsineTransform_t, h - 2, w - 2);
idst(invsineTransform_t, invsineTransform, w - 2, h - 2);
transpose(invsineTransform, invsineTransform_t, w - 2, h - 2);
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*w + j;
img_d[idx] = (double)myimg.at<uchar>(i, j);
}
}
for (int i = 1; i < h - 1; i++)
{
for (int j = 1; j < w - 1; j++)
{
idx = i*w + j;
img_d[idx] = 0.0;
}
}
for (int i = 1, id1 = 0; i < h - 1; i++, id1++)
{
for (int j = 1, id2 = 0; j < w - 1; j++, id2++)
{
idx = i*w + j;
idx1 = id1*(w - 2) + id2;
img_d[idx] = invsineTransform_t[idx1];
}
}
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
if(i>=0 && i<myresult.cols && j>=0 && j<myresult.rows)
{
idx = i*w + j;
if (img_d[idx] < 0.0)
myresult.at<uchar>(i, j) = 0;
else if (img_d[idx] > 255.0)
myresult.at<uchar>(i, j) = 255;
else
myresult.at<uchar>(i, j) = (uchar)img_d[idx];
}
}
}
delete[] sineTransform;
delete[] sineTransform_t;
delete[] denom;
delete[] invsineTransform;
delete[] invsineTransform_t;
delete[] img_d;
}
//泊松方程求解 输入散度(laplacianX+laplacianY)，及边界点像素即可重建求解
//输入：img为背景图片，laplacianX+laplacianY 为散度
//输出：result重建结果
void poissonSolver(const Mat &img, Mat &laplacianX , Mat &laplacianY, Mat &result)
{
const int w = img.cols;
const int h = img.rows;
Mat lap = Mat(img.size(),CV_32FC1);
lap = laplacianX + laplacianY;//散度
Mat bound = img.clone();
//边界修正，opencv为了方便，直接把图片最外围的像素点排除在外，不参与泊松重建
rectangle(bound,Point(1,1),Point(img.cols-2,img.rows-2),Scalar::all(0),-1);
//Mat boundary_points;
//Laplacian(bound, boundary_points, CV_32F);
//boundary_points = lap - boundary_points;
//Mat mod_diff = boundary_points(Rect(1, 1, w-2, h-2));
//solve(img,mod_diff,result);
//把上面四句和下面的换掉
unsigned long int idx;
double *boundary_point = new double[h*w];
Mat bound_map = cv::Mat(img.size(), CV_8UC1, cv::Scalar(0));
for (int i = 1; i < h - 1; i++)
for (int j = 1; j < w - 1; j++)
{
idx = i*w + j;
boundary_point[idx] = -4 * (int)bound.at<uchar>(i, j) + (int)bound.at<uchar>(i, (j + 1)) + (int)bound.at<uchar>(i, (j - 1))
+ (int)bound.at<uchar>(i - 1, j) + (int)bound.at<uchar>(i + 1, j);
bound_map.at<uchar>(i, j) = boundary_point[idx];
}
Mat diff = Mat(h, w, CV_32FC1);
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
idx = i*w + j;
diff.at<float>(i, j) = (float)(lap.at<float>(i, j) - boundary_point[idx]);
}
}
double *mod_diff = new double[(h - 2)*(w - 2)];
for (int i = 0; i < h - 2; i++)
{
for (int j = 0; j < w - 2; j++)
{
idx = i*(w - 2) + j;
mod_diff[idx] = diff.at<float>(i + 1, j + 1);
}
}
solve(img, mod_diff, result);
delete[] mod_diff;
delete[] boundary_point;
}
//泊松重建
void poisson(const Mat &destination,const Mat &destinationGradientX,const Mat &destinationGradientY,const Mat &patchGradientX,const Mat &patchGradientY,Mat &outputB,Mat &outputG,Mat &outputR)
{
Mat laplacianX = Mat(destination.size(),CV_32FC3);
Mat laplacianY = Mat(destination.size(),CV_32FC3);
//因为前面已经对destinationGradientX做了固定区域的mask，patchGradientX做了修改区域的mask
laplacianX = destinationGradientX + patchGradientX;//求解整张图片新的梯度场
laplacianY = destinationGradientY + patchGradientY;
computeLaplacianX(laplacianX,laplacianX);//求解梯度的散度 也就是拉普拉坐标
computeLaplacianY(laplacianY,laplacianY);
vector<Mat> rgbx_channel,rgby_channel;
split(laplacianX,rgbx_channel);//通道拆分
split(laplacianY,rgby_channel);
vector<Mat> output1;
split(destination,output1);
//for(int chan =0;chan<3;chan++)
// poissonSolver(output1[chan], rgbx_channel[chan], rgby_channel[chan], output[chan]); //这里有问题说mat里面不许这样？！我改用下面的 大不了多写几句
Mat output1B,output1G,output1R,rgbx_channelB,rgbx_channelG,rgbx_channelR,rgby_channelB,rgby_channelG,rgby_channelR;
output1B=output1.at(0);
output1G=output1.at(1);
output1R=output1.at(2);
rgbx_channelB=rgbx_channel.at(0);
rgbx_channelG=rgbx_channel.at(1);
rgbx_channelR=rgbx_channel.at(2);
rgby_channelB=rgby_channel.at(0);
rgby_channelG=rgby_channel.at(1);
rgby_channelR=rgby_channel.at(2);
poissonSolver(output1B, rgbx_channelB, rgby_channelB, outputB);
poissonSolver(output1G, rgbx_channelG, rgby_channelG, outputG);
poissonSolver(output1R, rgbx_channelR, rgby_channelR, outputR);
//output.at(0)=B; //我这样对吗 平常做右值的可以作为左值吗 ？
//output.at(1)=G;
//output.at(2)=R;
}
//输入：I背景整张图片 wmask背景中待修改的区域的mask
//输出：cloned
void evaluate(const Mat &I,const Mat &wmask,const Mat &patchGradientX,const Mat &patchGradientY, Mat &cloned)
{
bitwise_not(wmask,wmask);//矩阵元素取反操作，这样背景图片保持不变的像素对应的mask值为1
Mat binaryMaskFloatInverted;
wmask.convertTo(binaryMaskFloatInverted,CV_32FC1,1.0/255.0);
//上面已经对patchGradientX做了mask操作 ，这边也对destinationGradientX做mask操作
Mat mydestinationGradientX,mydestinationGradientY;
computeGradientX(I,mydestinationGradientX);//计算背景图像的x方向梯度
computeGradientY(I,mydestinationGradientY);//计算背景图像y方向的梯度
arrayProduct(mydestinationGradientX,binaryMaskFloatInverted, mydestinationGradientX);
arrayProduct(mydestinationGradientY,binaryMaskFloatInverted, mydestinationGradientY);
Mat outputB,outputG,outputR; //改了这里！！！！！！！！！！！！！！！！
poisson(I,mydestinationGradientX,mydestinationGradientY,patchGradientX,patchGradientY,outputB,outputG,outputR);
Mat output[]={outputB,outputG,outputR};
merge(output,3,cloned);
}
//克隆融合外部接口函数
//输入：destination背景图片的整张图片 binaryMask为destination待修改的像素的mask
//patch是由src图片的ROI区域复制过来的图像，其大小与destination相同，patch只有binaryMask区域存的是src的ROI图片
//输出：cloned融合结果整张图片
void normalClone(const Mat &destination, const Mat &patch, const Mat &binaryMask, Mat &cloned)
{
const int w = destination.cols;
const int h = destination.rows;
const int channel = destination.channels();
const int n_elem_in_line = w * channel;
//计算destination在x，y方向的梯度，获得结果为：destinationGradientX destinationGradientY
//计算patch在x，y方向的梯度，获得结果为： patchGradientX patchGradientY
//同时对binaryMask进行边界腐蚀，去除毛刺，让边界变得光滑一点，同时把binaryMask归一化为0~1得binaryMaskFloat
//其实这样归一化后binaryMaskFloat只有数值0,1 1代表即将被修改的像素点
Mat mybinaryMaskFloat,mypatchGradientX,mypatchGradientY;
computeDerivatives(destination,patch,binaryMask,mybinaryMaskFloat);
computeGradientX(patch,mypatchGradientX);//计算ROI区域转换复制到destination一样大小的patch图片x方向梯度
computeGradientY(patch,mypatchGradientY);//计算y方向梯度
//因为patch是一个src包含ROI的最小矩形块图片
//patchGradientX与binaryMaskFloat相乘，这样patchGradientX就只剩下有用的区域了
arrayProduct(mypatchGradientX,mybinaryMaskFloat,mypatchGradientX);
arrayProduct(mypatchGradientY,mybinaryMaskFloat,mypatchGradientY);
evaluate(destination,binaryMask,mypatchGradientX,mypatchGradientY,cloned);
}
//无缝融合，这个函数其实没什么功能，只是为了减少、方便计算，先把要融合的区域裁剪下来，也就是对_src、_dst、_mask 进行裁剪到最小
//输入：_src前景图像 _dst背景图像 _mask前景图像的mask，
//p是用于对应用的，p点指的是在dst中，待融合区域中心点(dst的待融合区域的大小是根据包含熊的矩形的大小确定的)
//输出：_blend融合结果图片
//void seamlessClone(InputArray _src, InputArray _dst, InputArray _mask,Mat &mybinaryMaskFloat, Point p, OutputArray _blend, int flags)
void seamlessClone(InputArray _src, InputArray _dst, InputArray _mask, Point p,OutputArray _blend)
{
const Mat src = _src.getMat();
const Mat dest = _dst.getMat();
const Mat mask = _mask.getMat();
_blend.create(dest.size(), CV_8UC3);//融合后图片的大小肯定跟背景图像一样
Mat blend = _blend.getMat();
int minx = INT_MAX, miny = INT_MAX, maxx = INT_MIN, maxy = INT_MIN;
int h = mask.size().height;
int w = mask.size().width;
Mat gray = Mat(mask.size(),CV_8UC1);
Mat dst_mask = Mat::zeros(dest.size(),CV_8UC1);//背景图像的mask
Mat cs_mask = Mat::zeros(src.size(),CV_8UC3);
Mat cd_mask = Mat::zeros(dest.size(),CV_8UC3);
if(mask.channels() == 3)//如果给定的mask是彩色图 需要转换成单通道灰度图
cvtColor(mask, gray, COLOR_BGR2GRAY );
else
gray = mask;
//计算包含mask的最小矩形，也就是把那只熊包含起来的最小矩形框，这个矩形是位于src的，后面还有一个对应的矩形位于dst 这两个for循环有错误？！！！
for(int i=0;i<h;i++)
{
for(int j=0;j<w;j++)
{
if(gray.at<uchar>(i,j) == 255)
{
minx = std::min(minx,i);
maxx = std::max(maxx,i);
miny = std::min(miny,j);
maxy = std::max(maxy,j);
}
}
}
int lenx = maxx - minx;//计算矩形的宽
int leny = maxy - miny;//计算矩形的高
Mat patch = Mat::zeros(Size(leny, lenx), CV_8UC3);//根据上面的矩形区域，创建一个大小相同矩阵
int minxd = p.y - lenx/2;//计算dst的矩形
int maxxd = p.y + lenx/2;
int minyd = p.x - leny/2;
int maxyd = p.x + leny/2;
CV_Assert(minxd>=0&&minyd>=0&&maxxd<=dest.rows&&maxyd<=dest.cols);
Rect roi_d(minyd,minxd,leny,lenx);//dst 兴趣区域的矩形
Rect roi_s(miny,minx,leny,lenx);//src 兴趣区域矩形
Mat destinationROI = dst_mask(roi_d);
Mat sourceROI = cs_mask(roi_s);
gray(roi_s).copyTo(destinationROI);//
src(roi_s).copyTo(sourceROI,gray(roi_s));
src(roi_s).copyTo(patch, gray(roi_s));//patch为
destinationROI = cd_mask(roi_d);
cs_mask(roi_s).copyTo(destinationROI);//cs_mask为把前景图片的矩形区域图像 转换到背景图片矩形中的图片
//Cloning obj;
//obj.normalClone(dest,cd_mask,dst_mask,blend,flags);
normalClone(dest,cd_mask,dst_mask,blend);
}
void main()
{
//第一步 读取两幅原图 并选定目标图中要融合的的区域
Mat src_primer_img=imread("airplane.jpg");
Mat dst_primer_img=imread("mountains.jpg");
if(!src_primer_img.data||!dst_primer_img.data)
cout<<"damn!wrong in reading!"<<endl;
Mat dst_img=dst_primer_img(Rect(100,300,src_primer_img.cols,src_primer_img.rows));
//imshow("dst",dst_img);
//第二步 对source即待嵌入的图进行腐蚀得到src_img_mask（这里得到的mask结果仍是彩图）
Mat src_img,src_img_mask;
Mat Kernel(Size(3, 3),CV_8UC1);
Kernel.setTo(Scalar(1));
erode(src_primer_img,src_img, Kernel, Point(-1,-1), 3); //对binaryMask进行腐蚀，去掉边界的毛刺点，让边界曲线平滑一点
src_img.convertTo(src_img_mask,CV_32FC1,1.0/255.0);
//computeDerivatives(dst_img,src_img,src_img,src_img_mask);
//imshow("srcmask",src_img_mask);
//imwrite("mask.jpg",src_img_mask);
//第三步
Point p=Point(100,300); //原始大图 要自己定点 中心点 矩形的！！！
Mat blend;
//seamlessClone(src_img,dst_primer_img,src_img_mask,p,blend); 替换成下面的 方便找出从哪里开始错误
Mat src;
src_img.copyTo(src);
Mat dest;
dst_primer_img.copyTo(dest);
Mat mask;
//src_img_mask.copyTo(mask);
mask=Mat::zeros(src_img_mask.size(),CV_32FC1);
//_blend.create(dest.size(), CV_8UC3);//融合后图片的大小肯定跟背景图像一样
dest.copyTo(blend);
int minx = INT_MAX, miny = INT_MAX, maxx = INT_MIN, maxy = INT_MIN;
int h = mask.size().height;
int w = mask.size().width;
Mat gray = Mat(mask.size(),CV_8UC1);
Mat dst_mask = Mat::zeros(dest.size(),CV_8UC1);//背景图像的mask
Mat cs_mask = Mat::zeros(src.size(),CV_8UC3);
Mat cd_mask = Mat::zeros(dest.size(),CV_8UC3);
if(mask.channels() == 3)//如果给定的mask是彩色图 需要转换成单通道灰度图
cvtColor(mask,gray,COLOR_BGR2GRAY );
else
gray = mask;
//imwrite("gray.jpg",gray); //保存显示检查错误 为什么这里显示的是腐蚀后的灰度图 但保存后的图片打开却是个全黑屏呢？奇怪
//imshow("gray",gray);
//计算包含mask的最小矩形，也就是把那只熊包含起来的最小矩形框，这个矩形是位于src的，后面还有一个对应的矩形位于dst 这两个for循环开始出现错误 貌似不是访问像素点越界的问题啊？？if(i>=0&&i<gray.rows&&j>=0&&j<gray.cols)
for(int i=0;i<h;++i)
{
for(int j=0;j<w;++j)
{
if((int)gray.at<float>(i,j)==255)
{
minx = std::min(minx,i);
maxx = std::max(maxx,i);
miny = std::min(miny,j);
maxy = std::max(maxy,j);
}
}
}
int lenx = maxx - minx;//计算矩形的宽
int leny = maxy - miny;//计算矩形的高
Mat patch = Mat::zeros(Size(leny, lenx), CV_8UC3);//根据上面的矩形区域，创建一个大小相同矩阵
int minxd = p.y - lenx/2;//计算dst的矩形
int maxxd = p.y + lenx/2;
int minyd = p.x - leny/2;
int maxyd = p.x + leny/2;
CV_Assert(minxd>=0&&minyd>=0&&maxxd<=dest.rows&&maxyd<=dest.cols);
Rect roi_d(minyd,minxd,leny,lenx);//dst 兴趣区域的矩形
Rect roi_s(miny,minx,leny,lenx);//src 兴趣区域矩形
Mat destinationROI = dst_mask(roi_d);
Mat sourceROI = cs_mask(roi_s);
gray(roi_s).copyTo(destinationROI);
//src(roi_s).copyTo(sourceROI,gray(roi_s)); //这句开始出现错误！！！！！！！！！！！！！！！！！！改成下面这样
src(roi_s).copyTo(sourceROI);
//src(roi_s).copyTo(patch, gray(roi_s));//patch为
src(roi_s).copyTo(patch);
destinationROI = cd_mask(roi_d);
cs_mask(roi_s).copyTo(destinationROI);//cs_mask为把前景图片的矩形区域图像 转换到背景图片矩形中的图片
normalClone(dest,cd_mask,dst_mask,blend);
namedWindow("poissonresult",1);
imshow("poissonresult",blend);
waitKey(0);
}
但是中间结果即待镶嵌图的mask还是挺好的 比我在matlab下自己之前写的得到的要好：
重来
全部重来写的 继续用之前写的MATLAB版本的 就是之前构造系数矩阵太大没写了的 现在接着写 那个原来可以用传说中的稀疏矩阵解决 我以前只听过没真的看过稀疏矩阵 也没去了解过现在知道了：
这是两个原图 我主要是用泊松融合让那只飞机可以无缝融合进右图中 把左图得到的mask放进右图
然后我用我构造系数矩阵然后重新解方程组得到的结果是：
A=imread('F:\fisheye\mountains.jpg');
src=A;
[ma,na,ka]=size(A);
B=imread('F:\fisheye\airplane.jpg');
dst=B;
[mb,nb,kb]=size(B);
se=strel('diamond',10);
B_erode=imerode(B,se);
Berodelogical=im2bw(B_erode);
for i=1:mb
for j=1:nb
if(Berodelogical(i,j)==1)
Berodelogical(i,j)=0;
else
Berodelogical(i,j)=1;
end
end
end
dstX=200;
dstY=100;
for i=dstY:dstY+mb-1
for j=dstX:dstX+nb-1
ii=i-(dstY-1);
jj=j-(dstX-1);
if(Berodelogical(ii,jj)==1)
A(i,j,1)=B(ii,jj,1);
A(i,j,2)=B(ii,jj,2);
A(i,j,3)=B(ii,jj,3);
end
end
end
%上面是把飞机放进目标图里 然后对目标图的新ROI计算梯度 散度
ROI=uint8(zeros(mb,nb,3));
for i=dstY:dstY+mb-1
for j=dstX:dstX+nb-1
ii=i-(dstY-1);
jj=j-(dstX-1);
ROI(ii,jj,1)=A(i,j,1);
ROI(ii,jj,2)=A(i,j,2);
ROI(ii,jj,3)=A(i,j,3);
end
end
w=[0,-1,1];
ROIgradienty=imfilter(double(ROI),w,'conv');
ROIgradientx=imfilter(double(ROI),w','conv');
%接下来对梯度求偏导得到融合图像的散度 lap
kernel=[-1,1,0];
lapy=imfilter(ROIgradienty,kernel,'conv');
lapx=imfilter(ROIgradientx,kernel','conv');
lap=lapx+lapy;
%接下来求解系数矩阵
%先将lap四周的点填充为原来的像素值（约束值） 其它的仍然是散度
[m,n,k]=size(lap);
heightRegion=m;
widthRegion=n;
L=n;
J=m;
for i=1:m
for j=1:n
if((i==1&&j==1)||(i==1&&j==n)||(i==m&&j==1)||(i==m&&j==n))
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
else
if(i==1||i==m)
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
else
if(j==1||j==n)
lap(i,j,1)=A(i,j,1);
lap(i,j,2)=A(i,j,2);
lap(i,j,3)=A(i,j,3);
end
end
end
end
end
%然后将lap'排成一列 即Ax=b的b矩阵
for i=1:m
for j=1:n
lap1(i,j)=lap(i,j,1);
lap2(i,j)=lap(i,j,2);
lap3(i,j)=lap(i,j,3);
end
end
lapp1=lap1';
lapp2=lap2';
lapp3=lap3';
b1=lapp1(:);
b2=lapp2(:);
b3=lapp3(:);
b=[b1 b2 b3];
%计算系数矩阵Ax=b的A 这里用M
[mm,nn]=size(b);
M=spalloc(mm,mm,6*(mm-2*m-2*n+4));
for i=1:mm
for j=1:mm
if(i==j)
M(i,j)=1;
end
end
end
%上面是边界点的约束值仍然是它本身的像素值 就是乘以1它本身就好
for y=1:heightRegion
for x=1:widthRegion
if x~=1&&y~=1&&y~=heightRegion&&x~=widthRegion
i=(y-1)*L+x;
M(i,i)=-4;
M(i,i-1)=1;
M(i,i+1)=1;
M(i,i-L)=1;
M(i,i+L)=1;
end
end
end
% use bi-conjugate gradient method to solve the matrix
%x = bicg(M,b',[],400); %稀疏矩阵的专业解法迭代求解 传说中的双共轭梯度法 不太了解 不能用平常的X=A\b’好像
%x=bicg(M,b,1e-6,400); %每一行是一个点的RGB三值 理论上的 但实际不能这样 换成下面
x1=bicg(M,b1,1e-6,400);
x2=bicg(M,b2,1e-6,400);
x3=bicg(M,b3,1e-6,400);
% reshape x to become the pixel intensity of the region
imRegion1= reshape(x1,widthRegion,heightRegion);
imRegion2= reshape(x2,widthRegion,heightRegion);
imRegion3= reshape(x3,widthRegion,heightRegion);
imNew=src;
imNew(dstY:dstY+heightRegion-1,dstX:dstX+widthRegion-1,1)=imRegion1';
imNew(dstY:dstY+heightRegion-1,dstX:dstX+widthRegion-1,2)=imRegion2';
imNew(dstY:dstY+heightRegion-1,dstX:dstX+widthRegion-1,3)=imRegion3';
这个效果不太对！！！我不放飞机进去 放月亮
右图是对月亮做mask后直接相加的结果 左图是后面泊松重建后的效果 忽略两图月亮的毛刺 因为我腐蚀时边界没处理 感觉泊松融合这次还是不够正确啊 可是我检查了下还是按照http://blog.csdn.net/hjimce/article/details/45716603这个人的理论来的啊 后面求解泊松方程时 用了稀疏矩阵 周围的点约束值是它们本身的像素值
其它点就是散度 系数矩阵也对了啊 怎么还是达不到它的效果呢
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- stitching.cpp鱼眼图像拼接融合 源码分析
- 柱面投影matlab程序
- 图像拼接(不投影到柱面)(渐入渐出融合) matlab程序
- 最佳缝合线算法（图像融合）
- 多分辨率拼接算法（继最佳缝合线之后）
- 关于文章说明
- Java进阶(十六)使用new Date()和System.currentTimeMillis()获取当前时间戳
- ListView 重用CheckBox出现错乱
- larave5 生成数据库测试数据
- Java进阶(十六)使用new Date()和System.currentTimeMillis()获取当前时间戳
- 还是回文（dp）
- Android开发知识点总结
- 每一天短信发送次数每一分钟发送频率发送限制 hibernate数据库框架
- 【Java基础】增强for循环要注意陷阱
- MmGetSystemRoutineAddress 函数实现代码
- Samba服务器搭建案例
- 轻松把玩HttpClient之配置ssl，采用设置信任自签名证书实现https
- Android（11）——从源码的角度解析setContentView
- centos7的安装及配置
- windows中使用git和开源中国