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》

總結

以上是生活随笔為你收集整理的Seamless cloning泊松克隆的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。