Seamless cloning泊松克隆
原文地址:http://blog.csdn.net/hjimce/article/details/45716603
作者:hjimce
本篇博文主要讲解2004年Siggraph的经典paper:《Poisson Image Editing》,在图像融合领域,融合效果最牛逼的paper。讲这个算法,我没打算讲太多理论的公式,理论的东西,对于大部分数学比较差的人来说看了就头晕。什么散度、拉普拉斯算子、梯度场、泊松方程、泊松方程第一类边界条件(Dirichlet boundary)、泊松方程第二类边界条件(Neumann boundary),如果把这些公式贴上来,估计很多人还没看到算法是怎么实现的,就已经看不下去了。因此我将直接给出离散形式实现方法,算法流程。
开始这个算法前,我需要先讲解一个数学问题:
一、散度计算
现在假设一幅图像为3*3的单通道灰度图像:
![](https://img-blog.csdn.net/20150525102928874?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaGppbWNl/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
我们假设每一点的像素值为V,V(1)表示像素点1的值,那么我们可以定义像素点5的散度的计算公式为:
div(5)=[V(2)+V(4)+V(6)+V(8)]-4*V(5)
说白了就是通过拉普拉斯卷积核,进行卷积,就可以求解散度了。
![](https://img-blog.csdn.net/20150525103751157?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaGppbWNl/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
拉普拉斯卷积核
当然正规的过程应该是先求解像素点5的梯度值,然后在对梯度求导,这样就能得到散度,不过得到的结果其实就是上面的计算公式。
二、泊松重建
OK,现在如果我给定一张图像,那么是不是可以利用拉普拉卷积核,求解每个点的散度(这里需要先说一下,后面用于泊松方程求解散度的时候,应该先求梯度,然后再对梯度求导得到散度,不要直接用卷积核,不然会犯我之前的一个错误,我之前就是直接用卷积核求解散度,导致边界的地方出现了过渡不自然的现象)。
现在反过来,如果我给定每个像素点的散度,我要你求解每个像素点的值,要怎么求取。这便是泊松方程的灵魂了。为了更好的理解重建过程,我现在假设图像的大小是4*4的16个像素点图片,如下:
![](https://img-blog.csdn.net/20150525105137740?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaGppbWNl/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
ok,假设我给你像素点6、7、10、11的散度值div(6)、div(7)、div(10)、div(11),那么我们是不是可以列出如下4个方程:
[V(2)+V(5)+V(7)+V(10)]-4*V(6)=div(6)
[V(3)+V(6)+V(8)+V(11)]-4*V(7)=div(7)
[V(6)+V(9)+V(11)+V(14)]-4*V(10)=div(10)
[V(7)+V(10)+V(12)+V(15)]-4*V(11)=div(11)
这个时候,如果我们只有四个方程,可是里面有16个像素点,也就是说有16个未知数。因此单单靠上面的4个方程,就想把所有的像素值求解出来是不可能的,这样方程有无数多个解。因此我们需要添加约束方程,这个便是泊松重建方程的约束条件了。假设我们添加边界约束条件,也就是说如果我已经知道了上面那副图像最外围一圈的每个像素点的值u,这样我们就可以得到12个约束方程。即:
V(1)=u(1) V(2)=u(2) V(3)=u(3)
V(4)=u(4) V(5)=u(5) V(8)=u(8)
V(9)=u(9) V(12)=u(12) V(13)=u(13)
V(14)=u(14) V(15)=u(15)V(16)=u(16)
上面有12个方程,外加给定的散度4个方程,这样我们有16个方程。这样就可以求解方程组了,这样就能实现通过散度+边界约束条件,实现图像重建。这个便是泊松方程的主要过程。
OK,不管图像多大,如果我们已经知道图片最外一圈的像素值(约束条件),以及其它像素点的散度值,我们就能把这个方程给列出来,构建泊松方程,重建图像。如果到这里你都看懂了,那么我觉得其实已经可以开始写图像融合的算了,是不是觉得算法很简单。说白了就是要求解一个方程组。
因此泊松融合,说的再简单一点,就是构建方程组:
Ax=b
然后通过求解这个方程组得到每个像素点的值。而算法的整个过程可以说是怎么构建方程组的b值 ,而系数矩阵其实是一个系数矩阵,矩阵的每一行有五个非零元素,对应于拉普拉斯算法的卷积核。
三、泊松图像融合
泊松融合可以说是目前融合效果上等的算法,泊松融合对应的文献为《Poisson Image Editing》,这篇文献叫基于泊松方程的图像编辑,没有叫融合,是因为它的神奇功能不仅仅用于简单的融合,还有一大堆的神器功能,当年我看到这篇文献的时候,感觉相当神奇,这个算法唯一的缺点是求解泊松方程需要一定的时间,速度比较慢。我之前自己看这paper把这篇文献的代码写过一遍,然而当时结果会出现偏色现象,所以一直以为自己没有真正看懂这篇paper,而今重新回顾,才发现原来自己的思路没有错,就是一个参数搞错了,计算散度的时候没写对。我们知道,对于一个像素点的散度求解,其实就是拉普拉斯算子滤波的结果:
![](https://img-blog.csdn.net/20150525101221481?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaGppbWNl/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
拉普拉斯算子
因为泊松重建,其实就是求解方程组:
Ax=b
算法的整个过程在于求解系数稀疏矩阵A、及b。只要A、b求出来了,那么我们就可以求解方程组得到x,而x就是我们得到的融合结果的像素颜色值。
因此当时我想当然的以为,b的求解直接用拉普拉斯卷积核对源图像的兴趣区域(ROI),进行卷积就可以得到散度b的值,因为对于给定的一幅图像散度其实就是通过拉普拉斯卷积核进行卷积,得到的结果,就是每个像素点的散度。这种思路本没有错,然而这样会出现边界过渡不自然的现象。
而正确的思路应该是:求解ROI的梯度场Isrc,及背景图像不被修改的像素区域的梯度场Idst。然后通过Isrc+Idst得到整幅待重建图像的梯度场,最后才根据梯度场求解散度。所以千万不要偷懒,不要一步求解散度,要先把待重建图像的梯度场求好,再进行求解散度。
OK,再啰嗦一遍算法的流程,看一下下面的图片,
1、问题描述:
现在假设我们有图像g,,如下图所示:
![](https://img-blog.csdn.net/20150525203216907?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaGppbWNl/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
待克隆图像区域(ROI)
还有一张背景图片S:
![](https://img-blog.csdn.net/20150525204514091?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaGppbWNl/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
背景图片S
现在我们希望把图片g融合粘贴到s中,且实现自然融合的效果:
![](https://img-blog.csdn.net/20150525204820724?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaGppbWNl/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
2、算法实现:
步骤1、计算图像g的梯度场。通过差分的方法,可以求得图像g的梯度场v:
![](https://img-blog.csdn.net/20150525205041364?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaGppbWNl/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
ROI的梯度场
梯度场的求取知道怎么求吧?如果连这都不会,那真的需要把图像最基本的东西好好看一看,说的简单一点就是卷积,我们平时边缘检测的时候,就有用到过计算梯度的模长。在这里我贴一下OpenCV的泊松融合这一步的代码:
- computeGradientX(patch,patchGradientX);
- computeGradientY(patch,patchGradientY);
上面的patch变量,你可以简单的把它理解为图像g,然后计算梯度的函数为:
- void Cloning::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 Cloning::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);
- }
- }
这样我们就可以计算出g的梯度场V(patchGradientX,patchGradientY
)。
步骤2、计算背景图片的梯度场:
- computeGradientX(destination,destinationGradientX);
- computeGradientY(destination,destinationGradientY);
变量destination为背景图像。这样就得到了背景图片的梯度场(destinationGradientX,destinationGradientY),如下图,下图的梯度场我是随便画一画的。
![](https://img-blog.csdn.net/20150526095429996?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaGppbWNl/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
背景图片的梯度场
步骤3、计算融合图像的梯度场。计算完了以后,我们就直接把ROI的梯度场覆盖到S的梯度场上:
- Mat laplacianX = Mat(destination.size(),CV_32FC3);
- Mat laplacianY = Mat(destination.size(),CV_32FC3);
-
-
- laplacianX = destinationGradientX + patchGradientX;
- laplacianY = destinationGradientY + patchGradientY;
上面的代码需要注意的是destinationGradientX、destinationGradientY已经被做了mask操作,所以才能直接相加。具体的mask操作如下:
- arrayProduct(destinationGradientX,binaryMaskFloatInverted, destinationGradientX);
- arrayProduct(destinationGradientY,binaryMaskFloatInverted, destinationGradientY);
-
- void Cloning::arrayProduct(const cv::Mat& lhs, const cv::Mat& rhs, cv::Mat& result) const
- {
- vector <Mat> lhs_channels;
- vector <Mat> result_channels;
-
- split(lhs,lhs_channels);
- split(result,result_channels);
-
- for(int chan = 0 ; chan < 3 ; ++chan)
- multiply(lhs_channels[chan],rhs,result_channels[chan]);
-
- merge(result_channels,result);
- }
上面函数中binaryMaskFloatInverted是一个mask,即Ω区域的值为0,非Ω区域的值为1。
![](https://img-blog.csdn.net/20150526095559500?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaGppbWNl/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
待重建图像的梯度场
总之你只要把背景图片的Ω区域的梯度场直接替换为g的梯度场v就可以了,因此如果你前面想简化计算,其实背景图片Ω区域的梯度场是不需要计算的,因为这一块迟早会被g的梯度场替换掉,你只需要要计算背景图片不被覆盖的区域的梯度场就可以了。这一步就是得到待重建图像的梯度场。
步骤4、求解融合图像的散度。通过步骤3,我们可以得到每个像素点的梯度值,也就是待重建图像的梯度场,因此接着我们需要对梯度求偏导,从而获得散度。
- computeLaplacianX(laplacianX,laplacianX);
- computeLaplacianY(laplacianY,laplacianY);
其相关调用函数:
- void Cloning::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);
- }
其实:
- computeLaplacianX(laplacianX,laplacianX);
- computeLaplacianY(laplacianY,laplacianY);
这两句代码就是对梯度(laplacianX,laplacianY)在x和y方向上求偏导。因此最后散度的计算为:
- lap = laplacianX + laplacianY;
步骤5、求解系数矩阵。OK,第4步我们已经把散度计算完毕,回顾一下前面的泊松重建方程,Ax=b,b便是散度,因此接着我们需要只要构建系数矩阵,还有约束方程就ok了,这一步因为opencv的源码是用了泊松方程的快速求解的方法,它没有直接按我们的一般理解去求A,然后x=A-1*b。因为泊松方程有快速的求解方法,如果直接用求解A,然后求A得逆矩阵,那计算真不是一般的大。假如待重建图像的大小是1000*1000的,那么系数矩阵的大小就是(1000*1000)X(1000*1000)的方阵。虽然A最后是稀疏矩阵,但是这么庞大的矩阵,搞起来也要崩溃啊,其实也不是很慢,差不多也就几十秒钟的时间,计算机的计算速度感觉还是挺快的。这一步我贴一下其它的代码,因为opencv没有直接构建A矩阵,它是用了泊松方程的快速求解算法进行求解的,求解算法里面有正弦、余弦函数,因此我猜它是用FFT方法求解泊松方程的,具体我没有细看,也没有必要细看,因为这个求解方程不该是我们关注的重点,我们需要关注的是怎么构建这个方程。所以我还是得讲一下普通的解法,系数矩阵A到底是个什么玩意。
其实矩阵A,我前面已经提到过了。矩阵A的对角线的元素为-4,然后每行有对应的其它4个非零元素,其值为1,因为我们拉普拉斯卷积核的时候,就是这样搞的。还有一点我们图像边界像素点的值应该为1。为了简单理解,我现在回到博文最开始的部分:
如果一幅图像,除了边界像素点之外,上面3*3图像的边界像素点为1、2、3、4、6、7、8、9。其它像素点的散度(上图中的像素5)我都已经知道了。那么我就可以列出泊松方程:
[V(2)+V(4)+V(6)+V(8)]-4*V(5)=div(5)
然后如果在把一幅图像的边界像素点的像素值告诉你,那么你就可以求解泊松方程了,假设约束点的值为u。以上面3*3的图像为例,最后系数矩阵A的构造为:
![](https://img-blog.csdn.net/20150526102823772?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaGppbWNl/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
然后最后列出Ax=b的结果为:
![](https://img-blog.csdn.net/20150526103249415?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaGppbWNl/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
这样分别求解三个通道的方程,我们就可以获得每个点的像素R,G,B值了。再啰嗦一遍上面系数矩阵A的特点,图像最外围一圈的边界的对角线元素之为1,因为这些点是约束方程,其它的非边界点就直接根据拉普拉斯的卷积核就可以了。到了这里我觉得我应经讲的没法再详细了,就这样吧。opencv的源码如果你看不懂,建议看一下这个:http://eric-yuan.me/poisson-blending-2/
最后贴一下这个算法的神器融合效果:
![](https://img-blog.csdn.net/20150526104623862?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaGppbWNl/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
这篇博文只是讲了《Poisson Image Editing》第一个功能,普通无缝融合功能。后面将继续讲解其它神器功能的实现,敬请期待。
作者:hjimce 联系qq:1393852684 更多资源请关注我的博客:http://blog.csdn.net/hjimce 原创文章,转载请保留本行信息,如有错误,欢迎指正。
最后把opencv的完整版普通融合的代码贴在这里,供大家学习:
-
-
-
-
- void cv::seamlessClone(InputArray _src, InputArray _dst, InputArray _mask, Point p, OutputArray _blend, int flags)
- {
- 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);
- Mat cs_mask = Mat::zeros(src.size(),CV_8UC3);
- Mat cd_mask = Mat::zeros(dest.size(),CV_8UC3);
-
- if(mask.channels() == 3)
- cvtColor(mask, gray, COLOR_BGR2GRAY );
- else
- gray = mask;
-
- 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;
- 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);
- Rect roi_s(miny,minx,leny,lenx);
-
- 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));
-
- destinationROI = cd_mask(roi_d);
- cs_mask(roi_s).copyTo(destinationROI);
-
-
- Cloning obj;
- obj.normalClone(dest,cd_mask,dst_mask,blend,flags);
-
- }
-
-
-
-
- void Cloning::normalClone(const Mat &destination, const Mat &patch, const Mat &binaryMask, Mat &cloned, int flag)
- {
- const int w = destination.cols;
- const int h = destination.rows;
- const int channel = destination.channels();
- const int n_elem_in_line = w * channel;
-
-
-
-
- computeDerivatives(destination,patch,binaryMask);
-
-
-
- arrayProduct(patchGradientX,binaryMaskFloat, patchGradientX);
- arrayProduct(patchGradientY,binaryMaskFloat, patchGradientY);
-
-
- evaluate(destination,binaryMask,cloned);
- }
- void Cloning::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 Cloning::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 Cloning::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 Cloning::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 Cloning::dst(const Mat& src, Mat& dest, bool invert)
- {
- Mat temp = Mat::zeros(src.rows, 2 * src.cols + 2, CV_32F);
-
- int flag = invert ? DFT_ROWS + DFT_SCALE + DFT_INVERSE: DFT_ROWS;
-
- src.copyTo(temp(Rect(1,0, src.cols, src.rows)));
-
- for(int j = 0 ; j < src.rows ; ++j)
- {
- float * tempLinePtr = temp.ptr<float>(j);
- const float * srcLinePtr = src.ptr<float>(j);
- for(int i = 0 ; i < src.cols ; ++i)
- {
- tempLinePtr[src.cols + 2 + i] = - srcLinePtr[src.cols - 1 - i];
- }
- }
-
- Mat planes[] = {temp, Mat::zeros(temp.size(), CV_32F)};
- Mat complex;
-
- merge(planes, 2, complex);
- dft(complex, complex, flag);
- split(complex, planes);
- temp = Mat::zeros(src.cols, 2 * src.rows + 2, CV_32F);
-
- for(int j = 0 ; j < src.cols ; ++j)
- {
- float * tempLinePtr = temp.ptr<float>(j);
- for(int i = 0 ; i < src.rows ; ++i)
- {
- float val = planes[1].ptr<float>(i)[j + 1];
- tempLinePtr[i + 1] = val;
- tempLinePtr[temp.cols - 1 - i] = - val;
- }
- }
-
- Mat planes2[] = {temp, Mat::zeros(temp.size(), CV_32F)};
-
- merge(planes2, 2, complex);
- dft(complex, complex, flag);
- split(complex, planes2);
-
- temp = planes2[1].t();
- dest = Mat::zeros(src.size(), CV_32F);
- temp(Rect( 0, 1, src.cols, src.rows)).copyTo(dest);
- }
-
- void Cloning::idst(const Mat& src, Mat& dest)
- {
- dst(src, dest, true);
- }
-
-
-
-
-
- void Cloning::solve(const Mat &img, Mat& mod_diff, Mat &result)
- {
-
- const int w = img.cols;
- const int h = img.rows;
-
-
- Mat res;
- dst(mod_diff, res);
-
- for(int j = 0 ; j < h-2; j++)
- {
- float * resLinePtr = res.ptr<float>(j);
- for(int i = 0 ; i < w-2; i++)
- {
- resLinePtr[i] /= (filter_X[i] + filter_Y[j] - 4);
- }
- }
-
- idst(res, mod_diff);
-
- unsigned char * resLinePtr = result.ptr<unsigned char>(0);
- const unsigned char * imgLinePtr = img.ptr<unsigned char>(0);
- const float * interpLinePtr = NULL;
-
-
- for(int i = 0 ; i < w ; ++i)
- result.ptr<unsigned char>(0)[i] = img.ptr<unsigned char>(0)[i];
-
- for(int j = 1 ; j < h-1 ; ++j)
- {
- resLinePtr = result.ptr<unsigned char>(j);
- imgLinePtr = img.ptr<unsigned char>(j);
- interpLinePtr = mod_diff.ptr<float>(j-1);
-
-
- resLinePtr[0] = imgLinePtr[0];
-
- for(int i = 1 ; i < w-1 ; ++i)
- {
-
-
- float value = interpLinePtr[i-1];
- if(value < 0.)
- resLinePtr[i] = 0;
- else if (value > 255.0)
- resLinePtr[i] = 255;
- else
- resLinePtr[i] = static_cast<unsigned char>(value);
- }
-
-
- resLinePtr[w-1] = imgLinePtr[w-1];
- }
-
-
- resLinePtr = result.ptr<unsigned char>(h-1);
- imgLinePtr = img.ptr<unsigned char>(h-1);
- for(int i = 0 ; i < w ; ++i)
- resLinePtr[i] = imgLinePtr[i];
- }
-
-
-
- void Cloning::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();
-
- 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);
- }
-
- void Cloning::initVariables(const Mat &destination, const Mat &binaryMask)
- {
- destinationGradientX = Mat(destination.size(),CV_32FC3);
- destinationGradientY = Mat(destination.size(),CV_32FC3);
- patchGradientX = Mat(destination.size(),CV_32FC3);
- patchGradientY = Mat(destination.size(),CV_32FC3);
-
- binaryMaskFloat = Mat(binaryMask.size(),CV_32FC1);
- binaryMaskFloatInverted = Mat(binaryMask.size(),CV_32FC1);
-
-
- const int w = destination.cols;
- filter_X.resize(w - 2);
- for(int i = 0 ; i < w-2 ; ++i)
- filter_X[i] = 2.0f * std::cos(static_cast<float>(CV_PI) * (i + 1) / (w - 1));
-
- const int h = destination.rows;
- filter_Y.resize(h - 2);
- for(int j = 0 ; j < h - 2 ; ++j)
- filter_Y[j] = 2.0f * std::cos(static_cast<float>(CV_PI) * (j + 1) / (h - 1));
- }
-
- void Cloning::computeDerivatives(const Mat& destination, const Mat &patch, const Mat &binaryMask)
- {
- initVariables(destination,binaryMask);
-
- computeGradientX(destination,destinationGradientX);
- computeGradientY(destination,destinationGradientY);
-
- computeGradientX(patch,patchGradientX);
- computeGradientY(patch,patchGradientY);
-
- Mat Kernel(Size(3, 3), CV_8UC1);
- Kernel.setTo(Scalar(1));
- erode(binaryMask, binaryMask, Kernel, Point(-1,-1), 3);
-
- binaryMask.convertTo(binaryMaskFloat,CV_32FC1,1.0/255.0);
- }
-
- void Cloning::scalarProduct(Mat mat, float r, float g, float b)
- {
- vector <Mat> channels;
- split(mat,channels);
- multiply(channels[2],r,channels[2]);
- multiply(channels[1],g,channels[1]);
- multiply(channels[0],b,channels[0]);
- merge(channels,mat);
- }
-
- void Cloning::arrayProduct(const cv::Mat& lhs, const cv::Mat& rhs, cv::Mat& result) const
- {
- vector <Mat> lhs_channels;
- vector <Mat> result_channels;
-
- split(lhs,lhs_channels);
- split(result,result_channels);
-
- for(int chan = 0 ; chan < 3 ; ++chan)
- multiply(lhs_channels[chan],rhs,result_channels[chan]);
-
- merge(result_channels,result);
- }
-
- void Cloning::poisson(const Mat &destination)
- {
- Mat laplacianX = Mat(destination.size(),CV_32FC3);
- Mat laplacianY = Mat(destination.size(),CV_32FC3);
-
-
- laplacianX = destinationGradientX + patchGradientX;
- laplacianY = destinationGradientY + patchGradientY;
-
- computeLaplacianX(laplacianX,laplacianX);
- computeLaplacianY(laplacianY,laplacianY);
-
- split(laplacianX,rgbx_channel);
- split(laplacianY,rgby_channel);
-
- split(destination,output);
-
- for(int chan = 0 ; chan < 3 ; ++chan)
- {
- poissonSolver(output[chan], rgbx_channel[chan], rgby_channel[chan], output[chan]);
- }
- }
-
-
- void Cloning::evaluate(const Mat &I, const Mat &wmask, const Mat &cloned)
- {
- bitwise_not(wmask,wmask);
-
- wmask.convertTo(binaryMaskFloatInverted,CV_32FC1,1.0/255.0);
-
- arrayProduct(destinationGradientX,binaryMaskFloatInverted, destinationGradientX);
- arrayProduct(destinationGradientY,binaryMaskFloatInverted, destinationGradientY);
-
- poisson(I);
-
- merge(output,cloned);
- }
上面的代码对应于opencv的这个算法的第一个功能“Normal Cloning”,后面还有五大神奇的功能。其具体功能选项如下:
具体使用看文献《Poisson Image Editing》
* 1- Normal Cloning
* 2- Mixed Cloning 1与2的区别见文献图片6
* 3- Monochrome Transfer 细节风格转换文献中的图片5,这个就像paper《Style Transfer for Headshot Portraits》一样的功能
* 4- Color Change 文献图片11
* 5- Illumination change 文献图片10
* 6- Texture Flattening 文献图片9
参考文献:
1、Opencv3.0
2、《Poisson Image Editing》