Seamless cloning泊松克隆


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


作者：hjimce


本篇博文主要講解2004年Siggraph的經(jīng)典paper：《Poisson Image Editing》，在圖像融合領(lǐng)域，融合效果最牛逼的paper。講這個(gè)算法，我沒打算講太多理論的公式，理論的東西，對(duì)于大部分?jǐn)?shù)學(xué)比較差的人來(lái)說(shuō)看了就頭暈。什么散度、拉普拉斯算子、梯度場(chǎng)、泊松方程、泊松方程第一類邊界條件(Dirichlet boundary)、泊松方程第二類邊界條件(Neumann boundary)，如果把這些公式貼上來(lái)，估計(jì)很多人還沒看到算法是怎么實(shí)現(xiàn)的，就已經(jīng)看不下去了。因此我將直接給出離散形式實(shí)現(xiàn)方法，算法流程。


開始這個(gè)算法前，我需要先講解一個(gè)數(shù)學(xué)問(wèn)題：


一、散度計(jì)算


現(xiàn)在假設(shè)一幅圖像為3*3的單通道灰度圖像：






我們假設(shè)每一點(diǎn)的像素值為V，V(1)表示像素點(diǎn)1的值，那么我們可以定義像素點(diǎn)5的散度的計(jì)算公式為：


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


說(shuō)白了就是通過(guò)拉普拉斯卷積核，進(jìn)行卷積，就可以求解散度了。






拉普拉斯卷積核


當(dāng)然正規(guī)的過(guò)程應(yīng)該是先求解像素點(diǎn)5的梯度值，然后在對(duì)梯度求導(dǎo)，這樣就能得到散度，不過(guò)得到的結(jié)果其實(shí)就是上面的計(jì)算公式。


二、泊松重建


OK，現(xiàn)在如果我給定一張圖像，那么是不是可以利用拉普拉卷積核，求解每個(gè)點(diǎn)的散度(這里需要先說(shuō)一下，后面用于泊松方程求解散度的時(shí)候，應(yīng)該先求梯度，然后再對(duì)梯度求導(dǎo)得到散度，不要直接用卷積核，不然會(huì)犯我之前的一個(gè)錯(cuò)誤，我之前就是直接用卷積核求解散度，導(dǎo)致邊界的地方出現(xiàn)了過(guò)渡不自然的現(xiàn)象)。


現(xiàn)在反過(guò)來(lái)，如果我給定每個(gè)像素點(diǎn)的散度，我要你求解每個(gè)像素點(diǎn)的值，要怎么求取。這便是泊松方程的靈魂了。為了更好的理解重建過(guò)程，我現(xiàn)在假設(shè)圖像的大小是4*4的16個(gè)像素點(diǎn)圖片，如下：






ok，假設(shè)我給你像素點(diǎn)6、7、10、11的散度值div(6)、div(7)、div(10)、div(11)，那么我們是不是可以列出如下4個(gè)方程：


[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)


這個(gè)時(shí)候，如果我們只有四個(gè)方程，可是里面有16個(gè)像素點(diǎn)，也就是說(shuō)有16個(gè)未知數(shù)。因此單單靠上面的4個(gè)方程，就想把所有的像素值求解出來(lái)是不可能的，這樣方程有無(wú)數(shù)多個(gè)解。因此我們需要添加約束方程，這個(gè)便是泊松重建方程的約束條件了。假設(shè)我們添加邊界約束條件，也就是說(shuō)如果我已經(jīng)知道了上面那副圖像最外圍一圈的每個(gè)像素點(diǎn)的值u，這樣我們就可以得到12個(gè)約束方程。即：


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個(gè)方程，外加給定的散度4個(gè)方程，這樣我們有16個(gè)方程。這樣就可以求解方程組了，這樣就能實(shí)現(xiàn)通過(guò)散度+邊界約束條件，實(shí)現(xiàn)圖像重建。這個(gè)便是泊松方程的主要過(guò)程。


OK，不管圖像多大，如果我們已經(jīng)知道圖片最外一圈的像素值(約束條件)，以及其它像素點(diǎn)的散度值，我們就能把這個(gè)方程給列出來(lái)，構(gòu)建泊松方程，重建圖像。如果到這里你都看懂了，那么我覺得其實(shí)已經(jīng)可以開始寫圖像融合的算了，是不是覺得算法很簡(jiǎn)單。說(shuō)白了就是要求解一個(gè)方程組。


因此泊松融合，說(shuō)的再簡(jiǎn)單一點(diǎn)，就是構(gòu)建方程組：


Ax=b


然后通過(guò)求解這個(gè)方程組得到每個(gè)像素點(diǎn)的值。而算法的整個(gè)過(guò)程可以說(shuō)是怎么構(gòu)建方程組的b值 ，而系數(shù)矩陣其實(shí)是一個(gè)系數(shù)矩陣，矩陣的每一行有五個(gè)非零元素，對(duì)應(yīng)于拉普拉斯算法的卷積核。


三、泊松圖像融合


泊松融合可以說(shuō)是目前融合效果上等的算法，泊松融合對(duì)應(yīng)的文獻(xiàn)為《Poisson Image Editing》，這篇文獻(xiàn)叫基于泊松方程的圖像編輯，沒有叫融合，是因?yàn)樗纳衿婀δ懿粌H僅用于簡(jiǎn)單的融合，還有一大堆的神器功能，當(dāng)年我看到這篇文獻(xiàn)的時(shí)候，感覺相當(dāng)神奇，這個(gè)算法唯一的缺點(diǎn)是求解泊松方程需要一定的時(shí)間，速度比較慢。我之前自己看這paper把這篇文獻(xiàn)的代碼寫過(guò)一遍，然而當(dāng)時(shí)結(jié)果會(huì)出現(xiàn)偏色現(xiàn)象，所以一直以為自己沒有真正看懂這篇paper，而今重新回顧，才發(fā)現(xiàn)原來(lái)自己的思路沒有錯(cuò)，就是一個(gè)參數(shù)搞錯(cuò)了，計(jì)算散度的時(shí)候沒寫對(duì)。我們知道，對(duì)于一個(gè)像素點(diǎn)的散度求解，其實(shí)就是拉普拉斯算子濾波的結(jié)果：






拉普拉斯算子


因?yàn)椴此芍亟?#xff0c;其實(shí)就是求解方程組：


Ax=b


算法的整個(gè)過(guò)程在于求解系數(shù)稀疏矩陣A、及b。只要A、b求出來(lái)了，那么我們就可以求解方程組得到x，而x就是我們得到的融合結(jié)果的像素顏色值。


因此當(dāng)時(shí)我想當(dāng)然的以為，b的求解直接用拉普拉斯卷積核對(duì)源圖像的興趣區(qū)域(ROI)，進(jìn)行卷積就可以得到散度b的值，因?yàn)閷?duì)于給定的一幅圖像散度其實(shí)就是通過(guò)拉普拉斯卷積核進(jìn)行卷積，得到的結(jié)果，就是每個(gè)像素點(diǎn)的散度。這種思路本沒有錯(cuò)，然而這樣會(huì)出現(xiàn)邊界過(guò)渡不自然的現(xiàn)象。


而正確的思路應(yīng)該是：求解ROI的梯度場(chǎng)Isrc，及背景圖像不被修改的像素區(qū)域的梯度場(chǎng)Idst。然后通過(guò)Isrc+Idst得到整幅待重建圖像的梯度場(chǎng)，最后才根據(jù)梯度場(chǎng)求解散度。所以千萬(wàn)不要偷懶，不要一步求解散度，要先把待重建圖像的梯度場(chǎng)求好，再進(jìn)行求解散度。


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


1、問(wèn)題描述：


現(xiàn)在假設(shè)我們有圖像g，，如下圖所示：






待克隆圖像區(qū)域(ROI)


還有一張背景圖片S:






背景圖片S


現(xiàn)在我們希望把圖片g融合粘貼到s中，且實(shí)現(xiàn)自然融合的效果：








2、算法實(shí)現(xiàn)：


步驟1、計(jì)算圖像g的梯度場(chǎng)。通過(guò)差分的方法，可以求得圖像g的梯度場(chǎng)v：






ROI的梯度場(chǎng)


梯度場(chǎng)的求取知道怎么求吧？如果連這都不會(huì)，那真的需要把圖像最基本的東西好好看一看，說(shuō)的簡(jiǎn)單一點(diǎn)就是卷積，我們平時(shí)邊緣檢測(cè)的時(shí)候，就有用到過(guò)計(jì)算梯度的模長(zhǎng)。在這里我貼一下opencv的泊松融合這一步的代碼：




[cpp] view plain copy 在CODE上查看代碼片派生到我的代碼片
computeGradientX(patch,patchGradientX);//計(jì)算ROI區(qū)域轉(zhuǎn)換復(fù)制到destination一樣大小的patch圖片x方向梯度 ?
computeGradientY(patch,patchGradientY);//計(jì)算y方向梯度 ?
上面的patch變量，你可以簡(jiǎn)單的把它理解為圖像g，然后計(jì)算梯度的函數(shù)為：


[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); ?
? ? } ?
} ?


這樣我們就可以計(jì)算出g的梯度場(chǎng)V(patchGradientX,patchGradientY)。
步驟2、計(jì)算背景圖片的梯度場(chǎng)：




[cpp] view plain copy 在CODE上查看代碼片派生到我的代碼片
computeGradientX(destination,destinationGradientX);//計(jì)算背景圖像的x方向梯度 ?
computeGradientY(destination,destinationGradientY);//計(jì)算背景圖像y方向的梯度 ?
變量destination為背景圖像。這樣就得到了背景圖片的梯度場(chǎng)(destinationGradientX,destinationGradientY)，如下圖，下圖的梯度場(chǎng)我是隨便畫一畫的。






背景圖片的梯度場(chǎng)


步驟3、計(jì)算融合圖像的梯度場(chǎng)。計(jì)算完了以后，我們就直接把ROI的梯度場(chǎng)覆蓋到S的梯度場(chǎng)上：




[cpp] view plain copy 在CODE上查看代碼片派生到我的代碼片
? ?Mat laplacianX = Mat(destination.size(),CV_32FC3); ?
? ?Mat laplacianY = Mat(destination.size(),CV_32FC3); ?
??
//因?yàn)榍懊嬉呀?jīng)對(duì)destinationGradientX做了固定區(qū)域的mask，patchGradientX做了修改區(qū)域的mask ?
? ?laplacianX = destinationGradientX + patchGradientX;//求解整張圖片新的梯度場(chǎng) ?
? ?laplacianY = destinationGradientY + patchGradientY; ?


上面的代碼需要注意的是destinationGradientX、destinationGradientY已經(jīng)被做了mask操作，所以才能直接相加。具體的mask操作如下：


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


[cpp] view plain copy 在CODE上查看代碼片派生到我的代碼片
//矩陣點(diǎn)乘，將lhs與rhs點(diǎn)乘得到result，因?yàn)橛腥齻€(gè)通道，估計(jì)mat不能實(shí)現(xiàn)三通道的矩陣的一次性點(diǎn)乘，所以才有這個(gè)函數(shù) ?
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個(gè)通道的矩陣 ?
? ? split(result,result_channels); ?
??
? ? for(int chan = 0 ; chan < 3 ; ++chan)//三個(gè)矩陣進(jìn)行分別相乘 ?
? ? ? ? multiply(lhs_channels[chan],rhs,result_channels[chan]); ?
??
? ? merge(result_channels,result);//合成為一個(gè) ?
} ?


上面函數(shù)中binaryMaskFloatInverted是一個(gè)mask，即Ω區(qū)域的值為0，非Ω區(qū)域的值為1。






待重建圖像的梯度場(chǎng)


總之你只要把背景圖片的Ω區(qū)域的梯度場(chǎng)直接替換為g的梯度場(chǎng)v就可以了，因此如果你前面想簡(jiǎn)化計(jì)算，其實(shí)背景圖片Ω區(qū)域的梯度場(chǎng)是不需要計(jì)算的，因?yàn)檫@一塊遲早會(huì)被g的梯度場(chǎng)替換掉，你只需要要計(jì)算背景圖片不被覆蓋的區(qū)域的梯度場(chǎng)就可以了。這一步就是得到待重建圖像的梯度場(chǎng)。


步驟4、求解融合圖像的散度。通過(guò)步驟3，我們可以得到每個(gè)像素點(diǎn)的梯度值，也就是待重建圖像的梯度場(chǎng)，因此接著我們需要對(duì)梯度求偏導(dǎo)，從而獲得散度。


[cpp] view plain copy 在CODE上查看代碼片派生到我的代碼片
computeLaplacianX(laplacianX,laplacianX);//求解梯度的散度 也就是拉普拉坐標(biāo) ?
computeLaplacianY(laplacianY,laplacianY); ?
其相關(guān)調(diào)用函數(shù)：


[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); ?
} ?


其實(shí)：
[cpp] view plain copy 在CODE上查看代碼片派生到我的代碼片
computeLaplacianX(laplacianX,laplacianX);//求解梯度的散度 也就是拉普拉坐標(biāo) ?
computeLaplacianY(laplacianY,laplacianY); ?
這兩句代碼就是對(duì)梯度(laplacianX,laplacianY)在x和y方向上求偏導(dǎo)。因此最后散度的計(jì)算為：
[cpp] view plain copy 在CODE上查看代碼片派生到我的代碼片
lap = laplacianX + laplacianY;//散度 ?
步驟5、求解系數(shù)矩陣。OK，第4步我們已經(jīng)把散度計(jì)算完畢，回顧一下前面的泊松重建方程，Ax=b，b便是散度，因此接著我們需要只要構(gòu)建系數(shù)矩陣，還有約束方程就ok了，這一步因?yàn)閛pencv的源碼是用了泊松方程的快速求解的方法，它沒有直接按我們的一般理解去求A，然后x=A-1*b。因?yàn)椴此煞匠逃锌焖俚那蠼夥椒?#xff0c;如果直接用求解A，然后求A得逆矩陣，那計(jì)算真不是一般的大。假如待重建圖像的大小是1000*1000的，那么系數(shù)矩陣的大小就是(1000*1000)X(1000*1000)的方陣。雖然A最后是稀疏矩陣，但是這么龐大的矩陣，搞起來(lái)也要崩潰啊，其實(shí)也不是很慢，差不多也就幾十秒鐘的時(shí)間，計(jì)算機(jī)的計(jì)算速度感覺還是挺快的。這一步我貼一下其它的代碼，因?yàn)閛pencv沒有直接構(gòu)建A矩陣，它是用了泊松方程的快速求解算法進(jìn)行求解的，求解算法里面有正弦、余弦函數(shù)，因此我猜它是用FFT方法求解泊松方程的，具體我沒有細(xì)看，也沒有必要細(xì)看，因?yàn)檫@個(gè)求解方程不該是我們關(guān)注的重點(diǎn)，我們需要關(guān)注的是怎么構(gòu)建這個(gè)方程。所以我還是得講一下普通的解法，系數(shù)矩陣A到底是個(gè)什么玩意。


其實(shí)矩陣A，我前面已經(jīng)提到過(guò)了。矩陣A的對(duì)角線的元素為-4，然后每行有對(duì)應(yīng)的其它4個(gè)非零元素，其值為1，因?yàn)槲覀兝绽咕矸e核的時(shí)候，就是這樣搞的。還有一點(diǎn)我們圖像邊界像素點(diǎn)的值應(yīng)該為1。為了簡(jiǎn)單理解，我現(xiàn)在回到博文最開始的部分：




如果一幅圖像，除了邊界像素點(diǎn)之外，上面3*3圖像的邊界像素點(diǎn)為1、2、3、4、6、7、8、9。其它像素點(diǎn)的散度(上圖中的像素5)我都已經(jīng)知道了。那么我就可以列出泊松方程：
[V(2)+V(4)+V(6)+V(8)]-4*V(5)=div(5)


然后如果在把一幅圖像的邊界像素點(diǎn)的像素值告訴你，那么你就可以求解泊松方程了，假設(shè)約束點(diǎn)的值為u。以上面3*3的圖像為例，最后系數(shù)矩陣A的構(gòu)造為：






然后最后列出Ax=b的結(jié)果為：












這樣分別求解三個(gè)通道的方程，我們就可以獲得每個(gè)點(diǎn)的像素R,G,B值了。再啰嗦一遍上面系數(shù)矩陣A的特點(diǎn)，圖像最外圍一圈的邊界的對(duì)角線元素之為1，因?yàn)檫@些點(diǎn)是約束方程，其它的非邊界點(diǎn)就直接根據(jù)拉普拉斯的卷積核就可以了。到了這里我覺得我應(yīng)經(jīng)講的沒法再詳細(xì)了，就這樣吧。opencv的源碼如果你看不懂，建議看一下這個(gè)：http://eric-yuan.me/poisson-blending-2/ ??


最后貼一下這個(gè)算法的神器融合效果：






這篇博文只是講了《Poisson Image Editing》第一個(gè)功能，普通無(wú)縫融合功能。后面將繼續(xù)講解其它神器功能的實(shí)現(xiàn)，敬請(qǐng)期待。 ??


作者：hjimce ? ? 聯(lián)系qq：1393852684 ? 更多資源請(qǐng)關(guān)注我的博客：http://blog.csdn.net/hjimce ? ? ? ? ? ? ? ? ?原創(chuàng)文章，轉(zhuǎn)載請(qǐng)保留本行信息，如有錯(cuò)誤，歡迎指正。


最后把opencv的完整版普通融合的代碼貼在這里，供大家學(xué)習(xí)：




[cpp] view plain copy 在CODE上查看代碼片派生到我的代碼片
//無(wú)縫融合，這個(gè)函數(shù)其實(shí)沒什么功能，只是為了減少、方便計(jì)算，先把要融合的區(qū)域裁剪下來(lái)，也就是對(duì)_src、_dst、_mask 進(jìn)行裁剪到最小 ?
//輸入：_src前景圖像 ?_dst背景圖像 ?_mask前景圖像的mask， ?
//p是用于對(duì)應(yīng)用的，p點(diǎn)指的是在dst中，待融合區(qū)域中心點(diǎn)(dst的待融合區(qū)域的大小是根據(jù)包含熊的矩形的大小確定的) ?
//輸出：_blend融合結(jié)果圖片 ?
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是彩色圖 需要轉(zhuǎn)換成單通道灰度圖 ?
? ? ? ? cvtColor(mask, gray, COLOR_BGR2GRAY ); ?
? ? else ?
? ? ? ? gray = mask; ?
? ? //計(jì)算包含mask的最小矩形，也就是把那只熊包含起來(lái)的最小矩形框，這個(gè)矩形是位于src的，后面還有一個(gè)對(duì)應(yīng)的矩形位于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;//計(jì)算矩形的寬 ?
? ? int leny = maxy - miny;//計(jì)算矩形的高 ?
? ? ??
??
? ? Mat patch = Mat::zeros(Size(leny, lenx), CV_8UC3);//根據(jù)上面的矩形區(qū)域，創(chuàng)建一個(gè)大小相同矩陣 ?
??
? ? int minxd = p.y - lenx/2;//計(jì)算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 興趣區(qū)域的矩形 ?
? ? Rect roi_s(miny,minx,leny,lenx);//src 興趣區(qū)域矩形 ?
??
? ? 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為把前景圖片的矩形區(qū)域圖像 轉(zhuǎn)換到背景圖片矩形中的圖片 ?
??
? ? ??
? ? Cloning obj; ?
? ? obj.normalClone(dest,cd_mask,dst_mask,blend,flags); ?
??
} ?


[cpp] view plain copy 在CODE上查看代碼片派生到我的代碼片
//克隆融合外部接口函數(shù) ?
//輸入：destination背景圖片的整張圖片 ?binaryMask為destination待修改的像素的mask ? ?
//patch是由src圖片的ROI區(qū)域復(fù)制過(guò)來(lái)的圖像，其大小與destination相同，只有patch只有binaryMask區(qū)域存的是src的ROI圖片 ?
//輸出：cloned融合結(jié)果整張圖片 ?
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; ?
? ? //計(jì)算destination在x，y方向的梯度，獲得結(jié)果為：destinationGradientX destinationGradientY ?
? ? //計(jì)算patch在x，y方向的梯度，獲得結(jié)果為： patchGradientX patchGradientY ?
? ? //同時(shí)對(duì)binaryMask進(jìn)行邊界腐蝕，去除毛刺，讓邊界變得光滑一點(diǎn)，同時(shí)把binaryMask歸一化為0~1得binaryMaskFloat ?
? ? //其實(shí)這樣歸一化后binaryMaskFloat只有數(shù)值0,1 ? ? ? 1代表即將被修改的像素點(diǎn) ?
? ? computeDerivatives(destination,patch,binaryMask); ?
? ?
? ? //因?yàn)閜atch是一個(gè)src包含ROI的最小矩形塊圖片 ?
? ? //patchGradientX與binaryMaskFloat相乘，這樣patchGradientX就只剩下有用的區(qū)域了 ?
? ? 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最外圍的像素點(diǎn) ?
//也就是說(shuō)矩陣mod_diff的大小為(w-2)*(h-2),并且mod_diff最外圍值(散度)為0 ?
//輸出：result ?
//這個(gè)函數(shù)其實(shí)功能是快速求解泊松方程的一種方法，就是針對(duì)AX=B，由于泊松方程系數(shù)矩陣的特殊性 ?
//這個(gè)方程的過(guò)程我們不需要深入理解 ?
void Cloning::solve(const Mat &img, Mat& mod_diff, Mat &result) ?
{ ?
? ? ??
? ? const int w = img.cols; ?
? ? const int h = img.rows; ?
? ? //到了這里其實(shí)mod_diff的寬為w-2 ?高為h-2 ??
? ? ??
? ? Mat res; ?
? ? dst(mod_diff, res);//這個(gè)函數(shù)？ ?
??
? ? 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)，及邊界點(diǎn)像素即可重建求解 ??
//輸入：img為背景圖片，laplacianX+laplacianY 為散度 ?
//輸出：result重建結(jié)果 ?
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為了方便，直接把圖片最外圍的像素點(diǎn)排除在外，不參與泊松重建 ?
? ? 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待修改的區(qū)域的mask ??
void Cloning::computeDerivatives(const Mat& destination, const Mat &patch, const Mat &binaryMask) ?
{ ?
? ? initVariables(destination,binaryMask);//相關(guān)變量初始化，沒用的東西 ?
??
? ? computeGradientX(destination,destinationGradientX);//計(jì)算背景圖像的x方向梯度 ?
? ? computeGradientY(destination,destinationGradientY);//計(jì)算背景圖像y方向的梯度 ?
??
? ? computeGradientX(patch,patchGradientX);//計(jì)算ROI區(qū)域轉(zhuǎn)換復(fù)制到destination一樣大小的patch圖片x方向梯度 ?
? ? computeGradientY(patch,patchGradientY);//計(jì)算y方向梯度 ?
??
? ? Mat Kernel(Size(3, 3), CV_8UC1); ?
? ? Kernel.setTo(Scalar(1)); ?
? ? erode(binaryMask, binaryMask, Kernel, Point(-1,-1), 3);//對(duì)binaryMask進(jìn)行腐蝕，去掉邊界的毛刺點(diǎn)，讓邊界曲線平滑一點(diǎn) ?
??
? ? 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); ?
} ?
//矩陣點(diǎn)乘，將lhs與rhs點(diǎn)乘得到result，因?yàn)橛腥齻€(gè)通道，估計(jì)mat不能實(shí)現(xiàn)三通道的矩陣的一次性點(diǎn)乘，所以才有這個(gè)函數(shù) ?
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個(gè)通道的矩陣 ?
? ? split(result,result_channels); ?
??
? ? for(int chan = 0 ; chan < 3 ; ++chan)//三個(gè)矩陣進(jìn)行分別相乘 ?
? ? ? ? multiply(lhs_channels[chan],rhs,result_channels[chan]); ?
??
? ? merge(result_channels,result);//合成為一個(gè) ?
} ?
//泊松重建 ?
void Cloning::poisson(const Mat &destination) ?
{ ?
? ? Mat laplacianX = Mat(destination.size(),CV_32FC3); ?
? ? Mat laplacianY = Mat(destination.size(),CV_32FC3); ?
??
? ? //因?yàn)榍懊嬉呀?jīng)對(duì)destinationGradientX做了固定區(qū)域的mask，patchGradientX做了修改區(qū)域的mask ?
? ? laplacianX = destinationGradientX + patchGradientX;//求解整張圖片新的梯度場(chǎng) ?
? ? laplacianY = destinationGradientY + patchGradientY; ?
??
? ? computeLaplacianX(laplacianX,laplacianX);//求解梯度的散度 也就是拉普拉坐標(biāo) ?
? ? 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背景中待修改的區(qū)域的mask ?
//輸出：cloned ?
void Cloning::evaluate(const Mat &I, const Mat &wmask, const Mat &cloned) ?
{ ?
? ? bitwise_not(wmask,wmask);//矩陣元素取反操作，這樣背景圖片保持不變的像素對(duì)應(yīng)的mask值為1 ?
??
? ? wmask.convertTo(binaryMaskFloatInverted,CV_32FC1,1.0/255.0); ?
? ? //上面已經(jīng)對(duì)patchGradientX做了mask操作 ，這邊也對(duì)destinationGradientX做mask操作 ?
? ? arrayProduct(destinationGradientX,binaryMaskFloatInverted, destinationGradientX); ?
? ? arrayProduct(destinationGradientY,binaryMaskFloatInverted, destinationGradientY); ?
??
? ? poisson(I); ?
??
? ? merge(output,cloned); ?
} ?
上面的代碼對(duì)應(yīng)于opencv的這個(gè)算法的第一個(gè)功能“Normal Cloning”，后面還有五大神奇的功能。其具體功能選項(xiàng)如下：


具體使用看文獻(xiàn)《Poisson Image Editing》


* 1- Normal Cloning
* 2- Mixed Cloning ? ?1與2的區(qū)別見文獻(xiàn)圖片6
* 3- Monochrome Transfer ?細(xì)節(jié)風(fēng)格轉(zhuǎn)換文獻(xiàn)中的圖片5，這個(gè)就像paper《Style Transfer for Headshot Portraits》一樣的功能
* 4- Color Change ? 文獻(xiàn)圖片11
* 5- Illumination change 文獻(xiàn)圖片10
* 6- Texture Flattening ?文獻(xiàn)圖片9




參考文獻(xiàn)：


1、Opencv3.0


2、《Poisson Image Editing》

總結(jié)

以上是生活随笔為你收集整理的Seamless cloning泊松克隆的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。