泊松融合
来源：互联网 发布：女用催情药 知乎 编辑：程序博客网 时间：2024/05/01 15:02
拖拖拉拉快一个月了 这是探索出来的 但感觉效果不是特别好 我目前只能做到这样子了：想把左图的人放进右图的水池中 用泊松融合想实现无缝自然的效果：
左图是直接做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
%他其实也只用泊松融合来做图片嵌入 并不是完整的图像拼接的程序 这个要自己改写
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()时出错：
出错：显示说这个通道没有实际大小 ：
继续调试，结合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版本的 就是之前构造系数矩阵太大没写了的 现在接着写 那个原来可以用传说中的稀疏矩阵解决 我以前只听过没真的看过稀疏矩阵 也没去了解过现在知道了：
然后我用我构造系数矩阵然后重新解方程组得到的结果是：
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';
这个效果不太对！！！我不放飞机进去 放月亮
- 泊松图像融合(泊松融合)
- 泊松融合
- 泊松融合
- 泊松融合
- 泊松融合(Poisson blend)
- 泊松融合-Possion Blending
- 泊松图像融合与地形合成
- 图像算法---泊松融合算法集锦
- 泊松融合-opencv代码实现
- python实现泊松图像融合
- 泊松图像融合算法 C++ 实现
- 融合
- 融合
- 融合
- 融合
- 图像的泊松(Poisson)编辑、泊松融合完全详解（1）
- 图像的泊松(Poisson)编辑、泊松融合完全详解（2）
- 图像的泊松(Poisson)编辑、泊松融合完全详解
- 在 CentOS 6.7 上配置 Cisco AnyConnect VPN
- 提升页面转化率100招
- 自定义布局
- 堆栈
- jquery操作select下拉框的多种方法(选中,取值,赋值等) 转载
- 泊松融合
- 技术进步与管理创新
- 【Hibernate】Hibernate多表查询后，将查询结果转化为实体
- centos7的安装及配置
- setTimeout
- MmGetSystemRoutineAddress 函数实现代码
- Java进阶(十六)使用new Date()和System.currentTimeMillis()获取当前时间戳
- java使用map去重复
- 网络中的一些专业名词解释
