Seamless cloning泊松克隆

Seamless cloning泊松克隆


原文地址：http://blog.csdn.net/hjimce/article/details/45716603


作者：hjimce


本篇博文主要讲解2004年Siggraph的经典paper：《Poisson Image Editing》，在图像融合领域，融合效果最牛逼的paper。讲这个算法，我没打算讲太多理论的公式，理论的东西，对于大部分数学比较差的人来说看了就头晕。什么散度、拉普拉斯算子、梯度场、泊松方程、泊松方程第一类边界条件(Dirichlet boundary)、泊松方程第二类边界条件(Neumann boundary)，如果把这些公式贴上来，估计很多人还没看到算法是怎么实现的，就已经看不下去了。因此我将直接给出离散形式实现方法，算法流程。


开始这个算法前，我需要先讲解一个数学问题：


一、散度计算


现在假设一幅图像为3*3的单通道灰度图像：






我们假设每一点的像素值为V，V(1)表示像素点1的值，那么我们可以定义像素点5的散度的计算公式为：


div(5)=[V(2)+V(4)+V(6)+V(8)]-4*V(5)


说白了就是通过拉普拉斯卷积核，进行卷积，就可以求解散度了。






拉普拉斯卷积核


当然正规的过程应该是先求解像素点5的梯度值，然后在对梯度求导，这样就能得到散度，不过得到的结果其实就是上面的计算公式。


二、泊松重建


OK，现在如果我给定一张图像，那么是不是可以利用拉普拉卷积核，求解每个点的散度(这里需要先说一下，后面用于泊松方程求解散度的时候，应该先求梯度，然后再对梯度求导得到散度，不要直接用卷积核，不然会犯我之前的一个错误，我之前就是直接用卷积核求解散度，导致边界的地方出现了过渡不自然的现象)。


现在反过来，如果我给定每个像素点的散度，我要你求解每个像素点的值，要怎么求取。这便是泊松方程的灵魂了。为了更好的理解重建过程，我现在假设图像的大小是4*4的16个像素点图片，如下：






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，而今重新回顾，才发现原来自己的思路没有错，就是一个参数搞错了，计算散度的时候没写对。我们知道，对于一个像素点的散度求解，其实就是拉普拉斯算子滤波的结果：






拉普拉斯算子


因为泊松重建，其实就是求解方程组：


Ax=b


算法的整个过程在于求解系数稀疏矩阵A、及b。只要A、b求出来了，那么我们就可以求解方程组得到x，而x就是我们得到的融合结果的像素颜色值。


因此当时我想当然的以为，b的求解直接用拉普拉斯卷积核对源图像的兴趣区域(ROI)，进行卷积就可以得到散度b的值，因为对于给定的一幅图像散度其实就是通过拉普拉斯卷积核进行卷积，得到的结果，就是每个像素点的散度。这种思路本没有错，然而这样会出现边界过渡不自然的现象。


而正确的思路应该是：求解ROI的梯度场Isrc，及背景图像不被修改的像素区域的梯度场Idst。然后通过Isrc+Idst得到整幅待重建图像的梯度场，最后才根据梯度场求解散度。所以千万不要偷懒，不要一步求解散度，要先把待重建图像的梯度场求好，再进行求解散度。


OK，再啰嗦一遍算法的流程，看一下下面的图片，


1、问题描述：


现在假设我们有图像g，，如下图所示：






待克隆图像区域(ROI)


还有一张背景图片S:






背景图片S


现在我们希望把图片g融合粘贴到s中，且实现自然融合的效果：








2、算法实现：


步骤1、计算图像g的梯度场。通过差分的方法，可以求得图像g的梯度场v：






ROI的梯度场


梯度场的求取知道怎么求吧？如果连这都不会，那真的需要把图像最基本的东西好好看一看，说的简单一点就是卷积，我们平时边缘检测的时候，就有用到过计算梯度的模长。在这里我贴一下opencv的泊松融合这一步的代码：




[cpp] view plain copy 在CODE上查看代码片派生到我的代码片
computeGradientX(patch,patchGradientX);//计算ROI区域转换复制到destination一样大小的patch图片x方向梯度  
computeGradientY(patch,patchGradientY);//计算y方向梯度  
上面的patch变量，你可以简单的把它理解为图像g，然后计算梯度的函数为：


[cpp] view plain copy 在CODE上查看代码片派生到我的代码片
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、计算背景图片的梯度场：




[cpp] view plain copy 在CODE上查看代码片派生到我的代码片
computeGradientX(destination,destinationGradientX);//计算背景图像的x方向梯度  
computeGradientY(destination,destinationGradientY);//计算背景图像y方向的梯度  
变量destination为背景图像。这样就得到了背景图片的梯度场(destinationGradientX,destinationGradientY)，如下图，下图的梯度场我是随便画一画的。






背景图片的梯度场


步骤3、计算融合图像的梯度场。计算完了以后，我们就直接把ROI的梯度场覆盖到S的梯度场上：




[cpp] view plain copy 在CODE上查看代码片派生到我的代码片
   Mat laplacianX = Mat(destination.size(),CV_32FC3);  
   Mat laplacianY = Mat(destination.size(),CV_32FC3);  
  
//因为前面已经对destinationGradientX做了固定区域的mask，patchGradientX做了修改区域的mask  
   laplacianX = destinationGradientX + patchGradientX;//求解整张图片新的梯度场  
   laplacianY = destinationGradientY + patchGradientY;  


上面的代码需要注意的是destinationGradientX、destinationGradientY已经被做了mask操作，所以才能直接相加。具体的mask操作如下：


[cpp] view plain copy 在CODE上查看代码片派生到我的代码片
arrayProduct(destinationGradientX,binaryMaskFloatInverted, destinationGradientX);  
arrayProduct(destinationGradientY,binaryMaskFloatInverted, destinationGradientY);  


[cpp] view plain copy 在CODE上查看代码片派生到我的代码片
//矩阵点乘，将lhs与rhs点乘得到result，因为有三个通道，估计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);//拆分成3个通道的矩阵  
    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。






待重建图像的梯度场


总之你只要把背景图片的Ω区域的梯度场直接替换为g的梯度场v就可以了，因此如果你前面想简化计算，其实背景图片Ω区域的梯度场是不需要计算的，因为这一块迟早会被g的梯度场替换掉，你只需要要计算背景图片不被覆盖的区域的梯度场就可以了。这一步就是得到待重建图像的梯度场。


步骤4、求解融合图像的散度。通过步骤3，我们可以得到每个像素点的梯度值，也就是待重建图像的梯度场，因此接着我们需要对梯度求偏导，从而获得散度。


[cpp] view plain copy 在CODE上查看代码片派生到我的代码片
computeLaplacianX(laplacianX,laplacianX);//求解梯度的散度 也就是拉普拉坐标  
computeLaplacianY(laplacianY,laplacianY);  
其相关调用函数：


[cpp] view plain copy 在CODE上查看代码片派生到我的代码片
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);  
}  


其实：
[cpp] view plain copy 在CODE上查看代码片派生到我的代码片
computeLaplacianX(laplacianX,laplacianX);//求解梯度的散度 也就是拉普拉坐标  
computeLaplacianY(laplacianY,laplacianY);  
这两句代码就是对梯度(laplacianX,laplacianY)在x和y方向上求偏导。因此最后散度的计算为：
[cpp] view plain copy 在CODE上查看代码片派生到我的代码片
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的构造为：






然后最后列出Ax=b的结果为：












这样分别求解三个通道的方程，我们就可以获得每个点的像素R,G,B值了。再啰嗦一遍上面系数矩阵A的特点，图像最外围一圈的边界的对角线元素之为1，因为这些点是约束方程，其它的非边界点就直接根据拉普拉斯的卷积核就可以了。到了这里我觉得我应经讲的没法再详细了，就这样吧。opencv的源码如果你看不懂，建议看一下这个：http://eric-yuan.me/poisson-blending-2/   


最后贴一下这个算法的神器融合效果：






这篇博文只是讲了《Poisson Image Editing》第一个功能，普通无缝融合功能。后面将继续讲解其它神器功能的实现，敬请期待。   


作者：hjimce     联系qq：1393852684   更多资源请关注我的博客：http://blog.csdn.net/hjimce                  原创文章，转载请保留本行信息，如有错误，欢迎指正。


最后把opencv的完整版普通融合的代码贴在这里，供大家学习：




[cpp] view plain copy 在CODE上查看代码片派生到我的代码片
//无缝融合，这个函数其实没什么功能，只是为了减少、方便计算，先把要融合的区域裁剪下来，也就是对_src、_dst、_mask 进行裁剪到最小  
//输入：_src前景图像  _dst背景图像  _mask前景图像的mask，  
//p是用于对应用的，p点指的是在dst中，待融合区域中心点(dst的待融合区域的大小是根据包含熊的矩形的大小确定的)  
//输出：_blend融合结果图片  
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);//背景图像的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(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);  
  
}  


[cpp] view plain copy 在CODE上查看代码片派生到我的代码片
//克隆融合外部接口函数  
//输入：destination背景图片的整张图片  binaryMask为destination待修改的像素的mask    
//patch是由src图片的ROI区域复制过来的图像，其大小与destination相同，只有patch只有binaryMask区域存的是src的ROI图片  
//输出：cloned融合结果整张图片  
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;  
    //计算destination在x，y方向的梯度，获得结果为：destinationGradientX destinationGradientY  
    //计算patch在x，y方向的梯度，获得结果为： patchGradientX patchGradientY  
    //同时对binaryMask进行边界腐蚀，去除毛刺，让边界变得光滑一点，同时把binaryMask归一化为0~1得binaryMaskFloat  
    //其实这样归一化后binaryMaskFloat只有数值0,1       1代表即将被修改的像素点  
    computeDerivatives(destination,patch,binaryMask);  
   
    //因为patch是一个src包含ROI的最小矩形块图片  
    //patchGradientX与binaryMaskFloat相乘，这样patchGradientX就只剩下有用的区域了  
    arrayProduct(patchGradientX,binaryMaskFloat, patchGradientX);  
    arrayProduct(patchGradientY,binaryMaskFloat, patchGradientY);  
  
  
    evaluate(destination,binaryMask,cloned);  
}  


[cpp] view plain copy 在CODE上查看代码片派生到我的代码片
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);  
}  
//输入：img为背景图像、mod_diff为散度，mod_diff大小不包含img最外围的像素点  
//也就是说矩阵mod_diff的大小为(w-2)*(h-2),并且mod_diff最外围值(散度)为0  
//输出：result  
//这个函数其实功能是快速求解泊松方程的一种方法，就是针对AX=B，由于泊松方程系数矩阵的特殊性  
//这个方程的过程我们不需要深入理解  
void Cloning::solve(const Mat &img, Mat& mod_diff, Mat &result)  
{  
      
    const int w = img.cols;  
    const int h = img.rows;  
    //到了这里其实mod_diff的宽为w-2  高为h-2   
      
    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;  
  
     //first col  
    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);  
  
        //first row  
        resLinePtr[0] = imgLinePtr[0];  
  
        for(int i = 1 ; i < w-1 ; ++i)  
        {  
            //saturate cast is not used here, because it behaves differently from the previous implementation  
            //most notable, saturate_cast rounds before truncating, here it's the opposite.  
            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);  
        }  
  
        //last row  
        resLinePtr[w-1] = imgLinePtr[w-1];  
    }  
  
    //last col  
    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];  
}  
//泊松方程求解 输入散度(laplacianX+laplacianY)，及边界点像素即可重建求解   
//输入：img为背景图片，laplacianX+laplacianY 为散度  
//输出：result重建结果  
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();  
    //边界修正，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);  
}  
  
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);  
  
    //init of the filters used in the dst  
    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));  
}  
//binaryMask为图像destination待修改的区域的mask   
void Cloning::computeDerivatives(const Mat& destination, const Mat &patch, const Mat &binaryMask)  
{  
    initVariables(destination,binaryMask);//相关变量初始化，没用的东西  
  
    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);  
}  
  
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);  
}  
//矩阵点乘，将lhs与rhs点乘得到result，因为有三个通道，估计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);//拆分成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 Cloning::poisson(const Mat &destination)  
{  
    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);  
  
    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]);  
    }  
}  
//输入：I背景整张图片  wmask背景中待修改的区域的mask  
//输出：cloned  
void Cloning::evaluate(const Mat &I, const Mat &wmask, const Mat &cloned)  
{  
    bitwise_not(wmask,wmask);//矩阵元素取反操作，这样背景图片保持不变的像素对应的mask值为1  
  
    wmask.convertTo(binaryMaskFloatInverted,CV_32FC1,1.0/255.0);  
    //上面已经对patchGradientX做了mask操作 ，这边也对destinationGradientX做mask操作  
    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》