本實驗實作了論文《“GrabCut” — Interactive Foreground Extraction using Iterated Graph Cuts》中的圖像分割功能，利用疊代的graph cuts進行互動式的前景提取，是對靜态圖像高效的、互動的前景或背景分割，對圖像編輯具有重大的現實意義。
通過在想要提取的物體外圍拉出一個包圍物體的矩形框，就可以對圖像實作較好的分割，提取出我們像要得到的圖像部分。論文示例如下：
目前景和背景的顔色比較相近的時候，相似的部分可能會難以分割出來，這時我們可以通過前景或背景的畫刷進行互動，通過不同的按鍵告訴程式自己希望圖像該區域所屬的部分，并再次進行疊代，以達到更加理想的效果：
除此之外，為了使分割出來得到的圖像的分割邊緣更加的自然平滑，論文中提出了Border Matting算法，來對得到的圖像邊緣進行平滑，以得到更加理想的圖像效果。論文中的Matting算法對比示例圖像如下：
在本項目中，我實作了論文中的GrabCut算法與Border Matting算法，以助教提供的OpenCV架構代碼為基礎，實作了GrabCut與Border Matting算法的背景代碼，通過與OpenCV中相同的方法進行圖像的分離操作。
- 理論分析
- 圖像分割算法概述
圖像分割就是将圖像分割成若幹個特定的、具有獨特性質的區域并提出感興趣目标的技術和過程，是圖像處理到圖像分析的重要步驟。圖像分割較為常用的特征是圖像的灰階、顔色、紋理、形狀等，這些特征在同一區域呈現出相似性，而在不同區域呈現出明顯的差異性。下面首先介紹一下現有的圖像分割算法。
2.1.1 基于門檻值的圖像分割算法
基于門檻值的分割算法的基本思想是基于圖像的灰階特征來計算一個或者多個灰階門檻值，并将圖像中每個像素的灰階值與門檻值進行比較，最後将比較得到的結果分到合适的類别，是以在本算法中，最重要的部分就是如何找到一個函數來求解最佳的灰階門檻值，來保證圖像的分割效果。
2.1.2 基于邊緣的圖像分割算法
基于邊緣的分割算法涉及到了圖像的邊緣提取，是建立在邊緣灰階值會呈現出階躍型或屋頂型變化這一觀測基礎上的算法。圖像的邊緣指的是圖像中不同區域的邊界線的像素點集合，展現出了圖像灰階、顔色、紋理等特征的突變，通常情況下基于邊緣的圖像分割算法指的是基于灰階值的邊緣檢測，通常展現的特性是灰階值的階躍性或屋頂性。階躍性指的是邊緣兩邊的像素的灰階值出現明顯的差異，而屋頂性指的是灰階值出現上升或下降的轉折處，基于這一特性，可以使用微分算子進行邊緣檢測來确定邊緣，常用的算子有Canny算子、Sobel算子等。
2.1.3 基于區域的圖像分割算法
基于區域的分割算法按照圖像的相似性準則分為不同的區域，一般有種子區域生長法，區域分裂合并法，分水嶺法。
種子區域生長法從一組代表不同圖像區域的種子像素開始，不斷地将目前像素鄰域中符合要求的像素加入到區域像素中，并以該像素為種子繼續疊代，直到找不到新像素為止。
區域分裂合并法的基本思想是首先将圖像任意分為互不相交的區域，再按照相關準則進行分裂或合并，這種方法既适合灰階圖像分割也适用于紋理圖像分割。
分水嶺法是基于拓撲理論的數學形态學的分割方法，其基本思想是将圖像看成地質學上的拓撲打雞毛，圖像上每一點的像素灰階值就是該點的海拔高度，每個局部極小值和其影響的區域為集水盆，集水盆的邊緣為分水嶺。分水嶺算法對邊緣微弱的圖像有着較好的效果，但是圖像中的噪聲會使得分水嶺算法産生過分割的現象。
2.1.4 基于圖論的分割方法
這一類的算法将圖像的分割與圖像的最小割問題相關鍊，首先将圖像映射為帶權的無向圖G=<V, E>，圖中的每個節點對應于圖像中的每一個像素，每條邊連接配接着一對相鄰的像素，邊的權值表示了像素之間灰階、顔色、紋理等内容的非負相似度，而對圖像的一個分割就是對圖像的一次裁剪，被分割的每個區域對應着圖像的一個子圖，分割的最優原則就是讓分割後的子圖在子圖的内部維持相似度最大，而子圖之間的相似度最小，本次項目所實作的GrabCut就是基于圖論的分割方法。
2.1.5 基于能量泛函的分割方法
這一類的方法主要指的是活動輪廓模型以及在其基礎上發展出來的算法，基本思想是使用連續的曲線來表達目标邊緣，并且定義一個能量泛函來使得自變量包括邊緣曲線，是以分割過程就會變成求解能量泛函最小值的過程，一般可以通過求解函數對應的歐拉方程來實作，當能量達到最小的時候就是目标輪廓的位置。
2.GrabCut算法概述
2.2.1 Graph Cut圖割模型
Graph cut圖割模型将圖像的分割問題與圖像的最小割(min cut)問題相結合，将圖像轉化為一個無向圖來表示，無向圖為G=<V,E>，V為無向圖的頂點，E為無向圖的邊。其中無向圖的頂點和邊分為兩類，第一類的每一個頂點對應于圖像中的每一個像素，而第一類的邊為每兩個相鄰像素之間的邊，這種邊叫做n-links；而第二類的頂點有兩個，叫做S（源點source）點和T（彙點Sink）點，而每個第一類頂點都和這兩個終端頂點之間有連接配接，組成了第二類的邊，叫做t-links。
上圖即為圖像對應的無向圖，每個像素在無向圖中對應着一個頂點，另外還會有S點和T點，而在圖像分割中，我們通常使用s代表圖像的前景目标，使用t代表圖像的背景目标。圖像中的每一條邊都有一個非負的權值，每一次圖像的切分就意味着消耗掉所切割邊的所有權值之和，如果一次切割所删除的邊的權值之和最小，則這一次切割就成為最小割，而福特-富克森定理表明，網絡的最大流max flow與最小割min cut相等。是以由Boykov和Kolmogorov發明的max-flow/min-cut算法就可以用來獲得s-t圖的最小割。這個最小割把圖的頂點劃分為兩個不相交的子集S和T，其中s ∈S，t∈ T和S∪T=V 。這兩個子集就對應于圖像的前景像素集和背景像素集，那就相當于完成了圖像分割，是以圖像分割的問題就可以轉化為無向圖中每一條邊的權值的确定問題。
圖像的分割可以看成是像素的标記問題，目标的标記與背景的标記不同，一般來講我們将背景标記為0，而将目标标記為1，而我們想要的分割就是要将背景與目标分割開來，而此時分割的cost也是最小的。我們假設每個像素的label為L={l1,l2,...lp}，其中li為0或者1，代表着背景或者前景，當圖像的分割為L的時候，圖像的能量為：E(L)=aR(L)+B(L)，其中R(L)為區域項，B(L)為邊界項，E(L)就是圖像的能量函數，我們要求得的圖割的最小值來使得能量函數最小。
下面我們分别對區域項和邊界項進行分析和解讀。
區域項：
其中Rp(Lp)為像素p配置設定到标簽lp的懲罰，在GraphCut中，Rp(Lp)的權值通過比較像素p的灰階值和給定目标的和前景的灰階直方圖獲得，我們希望像素p配置設定的标簽是其機率最大的标簽，但是這時我們希望像素的能量最小，是以我們取機率的負對數。當像素p的灰階值屬于前景的機率大于背景的機率時，Rp(1)就會小于Rp(0)，是以像素屬于前景的能量會相對更小，是以像素就會被劃分到前景中。
邊界項：
p和q為鄰域像素，B<p,q>為像素p和q之間不連續的懲罰，如果p和q差異非常大，則B<p,q>就接近于0，反之他們屬于同一個分類的可能性就比較大，B<p,q>就越大，即能量越大。
綜上所述，無向圖中的像素頂點之間的邊的權值由邊界平滑項決定，其權值取決于兩個像素之間的差異，另一類邊為像素定點與S點和T點的邊的權值由區域能量項決定，這時就可以根據mincut對圖進行切割，找到使能量最小的圖割方案。
2.2.2 GrabCut算法分析
與Graph Cut相比，GrabCut使用了BGR三通道的高斯混合模型進行模組化，GrabCut的能量最小化是通過不斷進行分割估計和參數模型學習的互動疊代過程，并且允許不完全的标注。
2.2.2.1 RGB顔色空間與GMM模組化
在GrabCut中，我們使用5個高斯分量的GMM模型來對前景和背景進行模組化，是以存在一個額外的向量k={k1,...kn,...kn}，kn為第n個像素對應于哪個高斯分量，對于每一個像素，可以屬于前景高斯分量或者背景高斯分量。整個圖像的Gibbs能量為：
其格式與GraphCut的能量方程較為類似，其中U為區域能量項，表示了一個像素被分為前景或者背景像素的懲罰，同樣使用像素屬于前景或者背景的機率的負對數作為能量值。将GMM模型應用到公式中後，(7)式就變為了(9)式的格式，參數θ有三個，如(10)式所示，分别為每個高斯分量的權重，均值和協方差。是以确定了這三個值後，我們就可以将RGB顔色值帶入到前景和背景的GMM求得其屬于前景或者背景的機率，即确定了區域能量項。
V為邊界能量項，其計算公式為
邊界能量項的作用與GraphCut的作用相似，展現了鄰域像素m與n之間的不連續懲罰，如果兩個像素的差别大，則它們不太可能同屬于同一區域，則之間的能量較小，反之則可能同屬于同一區域，之間的能量就會較大。在RGB空間中，我們使用歐式距離來衡量兩個像素的相似性，即||zm-zn||。其中β由圖像的對比度決定，如果圖像的對比度較高，則即使兩個像素同屬一個區域，其歐式距離仍然可能比較大，這時候我們就需要一個β來縮小這種差别，反之亦然，這就使得邊界能量項的計算不會在受到圖像對比度的影響。γ為50。
2.2.2.2 疊代能量最小化分割圖像
GrabCut算法使用了疊代最小化，每一次疊代的過程都使得對前景和背景的GMM參數更優，也就讓圖像的分割更優。首先使用者通過手動框選一個初始化的方框，方框内部的像素為圖像的”可能前景”Tu，而方框外部的像素全部被設定為“背景”Tb。我們得到初步的分類後，就可以使用k-means算法将同屬于前景和背景的像素聚類為K類，即K個高斯模型，這樣每個高斯模型就有了像素樣本集，我們就可以根據這些像素求得每一個高斯模型的均值、權值和協方差。
在疊代最小化的過程中，為每個像素配置設定GMM中的高斯分量，找到每一個像素對應的高斯分量中機率最大的，并将其配置設定給相應的高斯分量；對于給定的圖像，我們需要學習其GMM參數，根據每個高斯分量中已經擁有的像素點，估計其均值和協方差和權重；然後我們對圖像進行最大流分割，重複上述操作直到收斂即可。
2.2.2.3 使用者互動操作
在經過了一次框選之後，我們可以獲得一個初步的分割結果，但是對于圖像前景與背景差異不太明顯的圖像，我們往往難以僅僅通過一次框選就得到令人滿意的結果，是以我們需要對圖像進行進一步的互動，手動指出圖像的部分前景部分和部分背景部分，來幫助程式分析得到更加合适的模型。
3. Border Matting算法概述
我們從前一步的GrabCut中得到的圖像前景資訊是硬邊界的，而Border Matting算法通過“邊界消光”機制，使得硬邊界分界周圍的一個條帶上允許透明，以達到讓圖像邊界更加自然的目的。
在Border Matting算法中，我們首先需要找到的是圖像的邊界資訊，因為我們要處理的部分就是圖像的邊界區域，讓圖像的過渡更加自然。我們一般選用圖像的mask來使用Canny算子進行邊緣檢測。擷取到儲存了邊緣資訊的圖像後，我們對圖像進行輪廓參數化，找到同屬于一個連續輪廓上的像素點，并且儲存所有輪廓像素點的資訊，為每個像素點配置設定一個單獨的編号，因為後續需要在每一個輪廓上像素點的基礎上進行擴充，以找到整個需要處理的條帶。在條帶查找的階段，我們一般使用寬度搜尋的方式，擷取到距離每個中心像素點的一定距離以内的像素點，共同組成了圖像前景的邊界條帶。
為了計算出圖像中像素點的alpha值，我們需要确定的參數為sigma與delta，而這兩個參數的擷取我們使用了動态規劃算法，以能量最小化為原則，找到使得能量函數值最小的每個像素點的delta與sigma，最終使用delta與sigma進行計算，就可以得到像素點的alpha值。
能量函數如圖所示：
其中分為D項與V項，D項為資料項，V項為平滑項。其中delta與sigma的參數通過最小化該能量函數來擷取。
V項的計算方法如下：
其中lambda1取值為50，lambda2取值為1000。
V項的作用是使得alpha的值随着t的增加而平滑的變化，以達到圖像的邊緣逐漸模糊的效果。在我們使用動态規劃算法進行最優的參數值計算的時候，我們通常将delta離散化為30個等級，将sigma離散化為10個等級，在計算中我們使用循環将每一個sigma與delta都代入到輪廓像素中進行計算，直到最後找到使得能量函數取值最小的delta和sigma值并儲存。
D項的計算方法如下：
其中N代表了高斯機率密度函數，其中的計算需要使用到均值和協方差，而u和Σ的計算如下：
我們的alpha計算是通過Sigmoid函數進行計算的，其中每個像素使用了其delta和sigma值，擷取到整張圖像的alphaMask後，将其與原圖像進行融合，就得到了border matting後的結果。
3. 代碼細節
3.1 GrabCut算法細節
3.1.1 Mask初始化
首先我們建立一個與圖像大小相同的mask，且mask中每一個像素點的值為GC_BGD，即均設定為背景。接下來我們判斷我們所繪制的矩形的位置，進行放越界限制，并且将矩形内部的像素設定為GC_PR_FGD，即為待定的前景。
3.1.2 初始化GMM模型
我們首先周遊整張圖檔的所有像素，對于可能是背景的像素全部作為背景的樣本像素，而對于所有可能是前景的像素都作為前景的樣本像素。
我們定義矩陣bgdSamplesMat，設定其大小為背景樣本數*3，3為圖像的顔色通道數。我們在這裡使用kmeans聚類的方法，将背景的全部樣本像素點進行聚類，通常類數為5個，輸出為bgdLabels，儲存輸出樣本集中每一個樣本對應的類标簽，對于前景像素采用相同的處理方法，獲得fgdLabels作為标簽。我們在得知每一個像素所屬的高斯模型之後，下一步需要估計高斯模型中每一個高斯模型的參數，調用AddSample函數，将前景和背景中的每一個像素點加入到其對應的bgdLabels中儲存的GMM高斯模型的編号中。所有的像素添加完畢後，我們調用GMMLearning函數，求解每一個高斯模型的權值、均值和方差。權值系數指的是目前樣本像素的個數占全部像素個數的比值，并且計算每一個高斯模型的協方差的逆矩陣與行列式，以便于後續步驟的計算。
3.1.3 求解參數
之後我們需要确定一些後續計算中需要應用到的參數。我們在論文中得到gamma為50，建構無向圖中使用的lambda為9*gamma，而beta需要我們進行計算得到，beta為Gibbs能量項中的平滑項，用于當圖像的對比度較高或較低時對圖像相鄰像素之間的能量數值進行調整。我們需要分析整幅圖像來确定geta的值。首先周遊整個圖像的全部像素，計算每個像素與其相鄰像素的歐式距離，再計算每兩個相鄰像素的歐式距離期望，最後計算完成後按照論文中的公式(5)計算得到beta值。
3.1.4 計算無向圖邊界能量項權重
得到計算所需要的參數後，我們可以開始計算無向圖中的每兩個頂點之間的邊的權重了，這相當于Gibbs能量中的第二項，即平滑項，由于是無向圖，我們需要計算八鄰域，而對于每一個頂點，我們需要建構的是四個方向的邊的權重。首先構造四個矩陣來存儲每個點與其左側、上方、左上方、右上方的權重，然後我們周遊圖像中的全部像素，我們首先求出目前像素點與邊所連接配接的像素點的像素值的差，然後調用論文中公式(4)進行計算，并将結果儲存在對應的矩陣中。
3.1.5 疊代最小化求解
在求得平滑項後，我們便可以開始疊代求解了。首先我們進行疊代最小化的第一步，為每個像素配置設定GMM中所屬的高斯模型，我們周遊所有像素，如果像素屬于前景，則将該像素配置設定給前景GMM中機率最大的高斯模型；如果像素屬于背景，則将該像素配置設定給背景GMM中機率最大的高斯模型，并将結果儲存在與圖像大小相同的矩陣compIdxs中。
在初步将所有像素點配置設定給高斯模型之後，我們調用LearnGMMs來進行疊代最小化求解的第二部，從每個高斯模型的像素樣本集中學習每個高斯模型參數。我們循環每個高斯模型，在每一次循環中，我們周遊每一個圖像像素點，如果找到圖像的所屬高斯模型為目前模型，則将該點的資訊添加到對應的前景或背景高斯模型中。全部循環完成後，我們再針對每一個高斯模型計算其權重、均值和方差。
在高斯模型參數學習完成後，我們需要建立無向圖。在這裡我使用了助教在課程網站上提供的Max-flow/min-cut庫，調用其中的接口完成對無向圖的建立。我們需要通過計算得到的能量項建構圖，圖的頂點為像素點，邊由兩部分組成，第一類邊為每個像素點頂點與S點和T點的邊，其權值由Gibbs能量項的第一項表示，而另一類邊為每個頂點與其鄰域頂點相連接配接的邊，這一類的邊的權值通過Gibbs能量項的第二項平滑項表示。
在計算過程中，我首先周遊圖像的全部像素，對于不确定是前景還是背景的頂點，首先計算第一個能量項，其計算公式為論文的公式(3)，對于已經确定是背景或前景的點，我們可以直接指派：确定為背景的頂點其與S點的權重為0，與T點的權重為lambda；确定為前景的頂點其與S點的權重為lambda，與T點的權重為0，權重計算完畢後我們直接調用函數add_tweights函數，在圖中添加該node與source和和sink點的權值。對于像素頂點，我們要獲得之前計算得到的平滑項能量值，并且調用add_edge函數将目前頂點與相鄰頂點之間的邊的權添加到圖中。
最後我們使用最小割最大流算法進行分割估計。我們周遊整個mask，如果目前像素點是不确定為前景或者背景的像素點，我們調用what_segment函數進行判斷，如果傳回值為1則為背景點，如果傳回值為0即為前景點。
接下來我們根據自己的需要将上述疊代最小化求解過程進行疊代即可，最終就可以得到一個圖像的前景提取結果。
3.2 Border Matting算法細節
3.2.1 圖像的邊緣檢測
首先我們需要對圖像的mask進行邊緣檢測。我們知道mask中的每一個像素值為0、1、2、3中的一個，分别代表着确定背景、确定前景、待定背景、待定前景。我們将mask與1進行與操作，得到新的圖像。接下來我們使用Canny算子對圖像進行邊緣檢測，得到圖像的邊緣資訊。
3.2.2 輪廓資訊參數化
我們使用深度優先的方法來對輪廓進行參數計算，這一步的作用是找到圖像中全部的處于輪廓上的點，并且将其資訊儲存起來，并且存儲目前圖像中所有非連續的輪廓的個數，将每個輪廓上的點進行歸類，即将位于同一輪廓線上的所有像素點歸為同一類，并為每一個輪廓上的像素點配置設定一個獨一無二的編号，為下面配置設定條帶的處理做準備。我們首先周遊輪廓圖像的全部像素點，如果目前像素點沒有被處理過并且是輪廓上的點，我們就使用深度優先搜尋來處理與目前點相鄰的周圍的點，并将目前的點傳入存儲了全部輪廓點的vector：contour中，為後續的處理做準備。在這裡我們處理每個目前輪廓像素點周圍的8個像素點，如果目标點沒有超過圖像邊界、在輪廓線上并且沒有被标記過，就将該點作為新的點，并以之為起點再次進行深度優先的搜尋，直到所有的輪廓點都被處理完畢。
3.2.3 建構條帶資訊
在這一步中，我們需要首先初始化Tu，使用一個無序圖來進行存儲，key值為目标點的坐标資訊，value值為目标點。首先，我們需要将所有輪廓上的點添加到strip中，是以我們周遊contour中的全部點，将所有的點加入到無序圖strip中，并且添加到向量tmpPointOnStrip中，便于後續的寬搜。周遊完成後，我們以目前已經添加到strip中的點為基礎，進行寬度優先的搜尋，尋找符合條件的新的點添加到strip中去。我們周遊目前已經添加到條帶中的點，并且限制目前點到其條帶中心的距離為6，如果超過寬度的限制則不予處理。接下來我們處理目前點相鄰的4個點，分别為上、下、左、右四個方向。如果目标點沒有超出圖像範圍，并且沒有被處理過，則将其到條帶中心的距離絕對值+1，如果目标點是背景點，則距離為負值，設定目标點在條帶中所屬區域和目前點區域相同，并将目标點加入到條帶中。
3.2.4 能量最小化計算
下面需要使用能量最小化方法來求得參數的值，參數即為論文中指出的sigma與delta。首先我們需要使用目标圖像的灰階圖，由目标圖像轉化得到。然後我們周遊在輪廓上全部的像素點，計算目前輪廓上的像素點在40*40的區域上的前背景的均值和方差。根據論文，我們将delta離散化為30個等級，将sigma離散化為10個等級，我們在循環中枚舉每一個delta與sigma，使用三維數組energyFuncVals記錄每個輪廓上的點在每個delta和sigma時的能量值。根據論文，能量函數由兩部分組成，分别為D項與V項。
我們首先調用函數計算D項，根據論文的公式(14)，我們需要計算從一個起始點開始，整個輪廓的資料項之和。首先我們計算起始點的資料項值，然後寬度搜尋以起始點為中心點的區域，周遊目前點相鄰的點，如果目标點屬于中心點的區域，就将該點的資料項值添加到sum中，直到循環結束，傳回D項結果。由于能量方程是累加的結果，是以如果目前輪廓點是第一個，則直接将得到的D值作為能量值添加到energyFuncVals數組中；否則我們需要枚舉index-1時的delta與sigma，如果目前點與上一點屬于同一個輪廓，則可以計算V項。V項的計算由論文中的公式(13)決定。
計算完畢後，我們将新的到的D項和V項累加在index-1時的能量值上，并與目前energyFuncVals[index][d0][s0]相比較，如果目前得到的值更小，就用目前值代替之前的值，并記錄目前的delta與sigma，直到循環結束。
循環結束後，我們需要找到整個圖像總能量的最小值，即能量最小時的delta和sigma值。我們需要周遊所有的deltaLevels和sigmaLevels，找到minEnergy，并記錄目前的delta和sigma。記錄總能量的delta和sigma後，我們利用該值，取得每個區域的delta和sigma，并儲存。
3.2.5 Mask的Alpha值計算
首先我們要周遊輪廓上的所有的點，使用Sigmoid函數(論文中的圖6c)計算alpha值，并将其指派給alpha圖像。接下來我們再次寬度優先搜尋所有條帶上的點，針對目前頂點相鄰的4個頂點，如果目标頂點沒有超出圖像範圍并且沒有被處理過，則同樣使用Sigmoid函數計算該點的alpha值，并添加到alpha圖像中。
在擷取到包含了Alpha的mask後，我們為了讓結果更加明顯，對圖像進行一次高斯平滑，使結果的顯示範圍增加。
3.2.6 圖像效果展示
在這一階段，我們将原圖與儲存了alpha值的mask相融合，主要使用了将每個通道的矩陣與alphaMask相乘，得到的結果merge後就顯示出了border matting的結果。
#include "GMM.h"
GMM::GMM(Mat& _model)
{
if (_model.empty())
{
_model.create(1, modelSize*componentsCount, CV_64FC1);
_model.setTo(Scalar(0));
}
model = _model;
coefs = model.ptr<double>(0);
mean = coefs + componentsCount;
cov = coefs + componentsCount * 4;
for (int i = 0; i < componentsCount; i++)
{
covDeterms.push_back(0);
}
for (int i = 0; i < componentsCount; i++)
{
CalInverseCovAndDeterm(i);
}
}
GMM::~GMM()
{
}
// 第三次被GMMLearning調用時報錯
void GMM::CalInverseCovAndDeterm(int index)
{
if (coefs[index] > 0)
{
double* tmpCovPtr = cov + index * 9;
double dtrm;
Mat m = (Mat_<double>(3, 3) << tmpCovPtr[0], tmpCovPtr[1], tmpCovPtr[2],
tmpCovPtr[3], tmpCovPtr[4], tmpCovPtr[5],
tmpCovPtr[6], tmpCovPtr[7], tmpCovPtr[8]);
dtrm = determinant(m);
covDeterms[index] = dtrm;
Mat inv;
invert(m, inv);
inverseCovMats[index][0][0] = inv.at<double>(0, 0);
inverseCovMats[index][0][1] = inv.at<double>(0, 1);
inverseCovMats[index][0][2] = inv.at<double>(0, 2);
inverseCovMats[index][1][0] = inv.at<double>(1, 0);
inverseCovMats[index][1][1] = inv.at<double>(1, 1);
inverseCovMats[index][1][2] = inv.at<double>(1, 2);
inverseCovMats[index][2][0] = inv.at<double>(2, 0);
inverseCovMats[index][2][1] = inv.at<double>(2, 1);
inverseCovMats[index][2][2] = inv.at<double>(2, 2);
}
}
void GMM::InitGMMLearning()
{
totalSampleNum = 0;
for (int i = 0; i < componentsCount; i++)
{
sampleNum[i] = 0;
for (int j = 0; j < 3; j++)
{
sums[i][j] = 0;
for (int k = 0; k < 3; k++)
{
prods[i][j][k] = 0;
}
}
}
}
// 為前景後者背景GMM的第i個高斯模型增加樣本像素
// 在sum和prod中增加樣本值
void GMM::AddSample(int index, Vec3f color)
{
for (int i = 0; i < 3; i++)
{
sums[index][i] += color[i];
for (int j = 0; j < 3; j++)
{
prods[index][i][j] += color[i] * color[j];
}
}
++sampleNum[index];
++totalSampleNum;
}
// 每個像素有3個通道，是以均值有3個，協方差有9個
void GMM::GMMLearning()
{
double variance = 0.01;
for (int i = 0; i < componentsCount; i++)
{
int num = sampleNum[i];
if (num == 0)
coefs[i] = 0;
else
{
// 權值系數為目前樣本像素個數占全部像素個數的比值
coefs[i] = (double)(num / totalSampleNum);
double* mVal = mean + 3 * i;
double* cVal = cov + 9 * i;
for (int j = 0; j < 3; j++)
{
mVal[j] = sums[i][j] / num;
for (int k = 0; k < 3; k++)
{
cVal[j * 3 + k] = prods[i][j][k] / num - mVal[j] * mVal[k];
}
}
Mat tmp = (Mat_<double>(3, 3) << cVal[0], cVal[1], cVal[2],
cVal[3], cVal[4], cVal[5],
cVal[6], cVal[7], cVal[8]);
double d = determinant(tmp);
// 如果行列式小于等于0，則增加白噪聲，防止退化成沒有逆矩陣的行列式
if (d <= 0)
{
cVal[0] += variance;
cVal[4] += variance;
cVal[8] += variance;
}
CalInverseCovAndDeterm(i);
}
}
}
int GMM::GetComponent(const Vec3d color)
{
int idx = 0;
double max = 0;
for (int i = 0; i < componentsCount; i++)
{
double tmp = GetProbability(i, color);
if (tmp > max)
{
idx = i;
max = tmp;
}
}
return idx;
}
double GMM::GetProbability(const Vec3d color)
{
double result = 0;
for (int i = 0; i < componentsCount; i++)
{
result += coefs[i] * GetProbability(i, color);
}
return result;
}
// 計算目前像素與目前高斯分量的均值的差
// 使用diff轉置*inverseCovMats*diff求得mul
double GMM::GetProbability(int index, const Vec3d color)
{
double result = 0;
if (coefs[index] > 0)
{
if (covDeterms[index] > 0)
{
double* mVal = mean + 3 * index;
Vec3d diff = color;
for (int i = 0; i < 3; i++)
{
diff[i] -= mVal[i];
}
double mul = 0;
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
mul += diff[i] * diff[j] * inverseCovMats[index][j][i];
}
}
// 多元變量的聯合高斯密度函數求解機率
result = 1.0 / sqrt(covDeterms[index]) * exp(-0.5f*mul);
}
else
{
}
}
return result;
}
#pragma once
#include<iostream>
#include<vector>
#include<opencv2\opencv.hpp>
using namespace std;
using namespace cv;
class GMM
{
public:
GMM(Mat& _model);
~GMM();
static const int componentsCount = 5;
void CalInverseCovAndDeterm(int index);
void InitGMMLearning();
void AddSample(int index, Vec3f color);
void GMMLearning();
int GetComponent(const Vec3d color);
double GetProbability(const Vec3d color);
double GetProbability(int index, const Vec3d color);
private:
int modelSize = 13;
Mat model;
double* cov;
double* coefs;
double* mean;
int sampleNum[componentsCount];
int totalSampleNum;
// 儲存了每個component的協方差的逆矩陣
double inverseCovMats[componentsCount][3][3] ;
// 儲存了每個component的協方差的行列式
vector<double> covDeterms;
// 每個component的sum
double sums[componentsCount][3];
double prods[componentsCount][3][3];
};
#include "GrabCut.h"
void GrabCut2D::CheckRectAndMask(Mat &_mask, Rect rect, int mode, Size imgSize)
{
if (mode == GC_INIT_WITH_RECT)
{
InitMaskWithRect(_mask, rect, imgSize);
}
else if (mode == GC_INIT_WITH_MASK)
{
}
else
{
}
}
void GrabCut2D::InitMaskWithRect(Mat & _mask, Rect rect, Size imgSize)
{
if (_mask.empty())
{
_mask.create(imgSize, CV_8UC1);
_mask.setTo(Scalar(GC_BGD));
}
else
{
_mask.setTo(GC_BGD);
}
// 防止出現矩形越界
if (rect.x < 0)
rect.x = 0;
if (rect.y < 0)
rect.y = 0;
if (rect.width > imgSize.width - rect.x)
rect.width = imgSize.width - rect.x;
if (rect.height > imgSize.height - rect.y)
rect.height = imgSize.height - rect.y;
_mask(rect).setTo(Scalar(GC_PR_FGD));
//cout << GC_PR_FGD << endl;
//cout << Scalar(GC_PR_FGD) << endl;
}
GrabCut2D::~GrabCut2D(void)
{
}
// k-means聚類初始化兩個GMM
void GrabCut2D::InitGMMs(Mat& img, Mat& mask, GMM& bgdGMM, GMM& fgdGMM)
{
const int iterCnt = 10;
const int kMeansType = KMEANS_PP_CENTERS;
// 記錄前景和背景的像素樣本集中每個像素對應GMM中的哪一個高斯模型
Mat bgdLabels, fgdLabels;
vector<Vec3f> bgdSamples, fgdSamples;
// 周遊整個圖檔的全部像素
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
// mask中所有可能是背景的像素都作為背景的樣本像素
if (mask.at<uchar>(i, j) == GC_BGD || mask.at<uchar>(i, j) == GC_PR_BGD)
{
bgdSamples.push_back((Vec3d)img.at<Vec3b>(i, j));
}
// mask中所有可能是前景的像素都作為前景的樣本像素
else
{
fgdSamples.push_back((Vec3d)img.at<Vec3b>(i, j));
}
}
}
Mat bgdSamplesMat;
bgdSamplesMat.create((int)bgdSamples.size(), 3, CV_32FC1);
bgdSamplesMat.setTo(bgdSamples[0][0]);
kmeans(bgdSamplesMat, GMM::componentsCount, bgdLabels,
TermCriteria(CV_TERMCRIT_ITER, iterCnt, 0), 0, kMeansType);
Mat fgdSamplesOutput;
fgdSamplesOutput.create((int)fgdSamples.size(), 3, CV_32FC1);
fgdSamplesOutput.setTo(fgdSamples[0][0]);
kmeans(fgdSamplesOutput, GMM::componentsCount, fgdLabels,
TermCriteria(CV_TERMCRIT_ITER, iterCnt, 0), 0, kMeansType);
// 已經确定每個像素所屬的高斯模型，進而估計GMM中每個高斯模型參數
bgdGMM.InitGMMLearning();
fgdGMM.InitGMMLearning();
for (int i = 0; i < bgdSamples.size(); i++)
{
bgdGMM.AddSample(bgdLabels.at<int>(i, 0), bgdSamples[i]);
}
bgdGMM.GMMLearning();
for (int i = 0; i < fgdSamples.size(); i++)
{
fgdGMM.AddSample(fgdLabels.at<int>(i, 0), fgdSamples[i]);
}
fgdGMM.GMMLearning();
//Mat bgdSamplesOutput((int)bgdSamples.size(), 3, CV_32FC1, &bgdSamples[0][0]);
}
double GrabCut2D::CalBeta(const Mat& img)
{
double beta = 0;
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
// 左側
if (j > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i, j - 1);
beta += diff.dot(diff);
}
// 左上方
if (i > 0 && j > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i - 1, j - 1);
beta += diff.dot(diff);
}
// 上方
if (i > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i - 1, j);
beta += diff.dot(diff);
}
// 右上方
if (i > 0 && j < img.cols - 1)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i - 1, j + 1);
beta += diff.dot(diff);
}
}
}
if (beta < 0)
{
beta = 0;
}
else
{
beta = 1.0f / (2 * beta / (4 * img.cols * img.rows - 3 * img.cols - 3 * img.rows + 2));
}
return beta;
}
// 計算圖中每個像素頂點的權值，即與其相連的8個像素的權值
void GrabCut2D::CalAllWeights(const Mat& img, Mat& leftWeight, Mat& upLeftWeight, Mat& upWeight, Mat& upRightWeight,
double beta, double gamma)
{
const double gammaNear = gamma;
const double gammaFar = gamma / sqrt(2.0);
leftWeight.create(img.rows, img.cols, CV_64FC1);
upLeftWeight.create(img.rows, img.cols, CV_64FC1);
upWeight.create(img.rows, img.cols, CV_64FC1);
upRightWeight.create(img.rows, img.cols, CV_64FC1);
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
// 左側權重
if (j > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i, j - 1);
leftWeight.at<double>(i, j) = gammaNear * exp(diff.dot(diff) * (beta * -1));
}
else
leftWeight.at<double>(i, j) = 0;
// 上方權重
if (i > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i - 1, j);
upWeight.at<double>(i, j) = gammaNear * exp(diff.dot(diff) * (beta * -1));
}
else
upWeight.at<double>(i, j) = 0;
// 左上方權重
if (i > 0 && j > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i - 1, j - 1);
upLeftWeight.at<double>(i, j) = gammaFar * exp(diff.dot(diff) * (beta * -1));
}
else
upLeftWeight.at<double>(i, j) = 0;
// 右上方權重
if (j < img.cols - 1 && i > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i - 1, j + 1);
upRightWeight.at<double>(i, j) = gammaFar * exp(diff.dot(diff) * (beta * -1));
}
else
upRightWeight.at<double>(i, j) = 0;
}
}
}
// 疊代最小化的步驟一，給每個像素配置設定GMM中所屬的高斯模型
void GrabCut2D::AssignGMMsComponents(const Mat& img, const Mat& mask, GMM& bgdGMM, GMM& fgdGMM,
Mat& compIdxs)
{
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
if (mask.at<uchar>(i, j) == GC_BGD || mask.at<uchar>(i, j) == GC_PR_BGD)
{
Vec3d v = (Vec3d)(img.at<Vec3b>(i, j));
compIdxs.at<int>(i, j) = bgdGMM.GetComponent(v);
}
else
{
Vec3d v = (Vec3d)(img.at<Vec3b>(i, j));
compIdxs.at<int>(i, j) = fgdGMM.GetComponent(v);
}
}
}
}
// 疊代最小化算法step 2：從每個高斯模型的像素樣本集中學習每個高斯模型的參數
void GrabCut2D::LearnGMMs(const Mat & img, const Mat & mask, const Mat & compIdxs, GMM & bgdGMM, GMM & fgdGMM)
{
bgdGMM.InitGMMLearning();
fgdGMM.InitGMMLearning();
for (int index = 0; index < GMM::componentsCount; index++)
{
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
if (compIdxs.at<int>(i, j) == index)
{
if (mask.at<uchar>(i, j) == GC_BGD || mask.at<uchar>(i, j) == GC_PR_BGD)
{
bgdGMM.AddSample(index, img.at<Vec3b>(i, j));
}
else
{
fgdGMM.AddSample(index, img.at<Vec3b>(i, j));
}
}
}
}
}
bgdGMM.GMMLearning();
fgdGMM.GMMLearning();
}
void GrabCut2D::ConstructGraph(const Mat & img, const Mat & mask, GMM & bgdGMM, GMM & fgdGMM, double lambda, const Mat & leftWeight, const Mat & upLeftWeight, const Mat & upWeight, const Mat upRightWeight, Graph<double, double, double>& graph)
{
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
int nodeIndex = graph.add_node();
// 計算每個頂點與source和sink連接配接的權值，即第一個能量項
double sourceEdge, sinkEdge;
// 對于不确定的點，計算第一個能量項
if (mask.at<uchar>(i, j) == GC_PR_BGD || mask.at<uchar>(i, j) == GC_PR_FGD)
{
sourceEdge = -log(bgdGMM.GetProbability(img.at<Vec3b>(i, j)));
sinkEdge = -log(fgdGMM.GetProbability(img.at<Vec3b>(i, j)));
}
// 對于确定的背景點直接指派
else if (mask.at<uchar>(i, j) == GC_BGD)
{
sourceEdge = 0;
sinkEdge = lambda;
}
else
{
sourceEdge = lambda;
sinkEdge = 0;
}
// 設定該node與source和sink點的權值
graph.add_tweights(nodeIndex, sourceEdge, sinkEdge);
// 左側頂點的edge
if (j > 0)
{
double cap = leftWeight.at<double>(i, j);
graph.add_edge(nodeIndex, nodeIndex - 1, cap, cap);
}
// 上方頂點的edge
if (i > 0)
{
double cap = upWeight.at<double>(i, j);
graph.add_edge(nodeIndex, nodeIndex - img.cols, cap, cap);
}
// 左上方頂點的edge
if (j > 0 && i > 0)
{
double cap = upLeftWeight.at<double>(i, j);
graph.add_edge(nodeIndex, nodeIndex - img.cols - 1, cap, cap);
}
// 右上方頂點的edge
if (i > 0 && j < img.cols - 1)
{
double cap = upRightWeight.at<double>(i, j);
graph.add_edge(nodeIndex, nodeIndex - img.cols + 1, cap, cap);
}
}
}
}
// 使用最小割最大流算法進行分割估計
void GrabCut2D::EstimateSegmentation(Graph<double, double, double>& graph, Mat & mask)
{
graph.maxflow();
for (int i = 0; i < mask.rows; i++)
{
for (int j = 0; j < mask.cols; j++)
{
if (mask.at<uchar>(i, j) == GC_PR_BGD || mask.at<uchar>(i, j) == GC_PR_FGD)
{
if (graph.what_segment(i * mask.cols + j))
{
// 當傳回1時為SINK
mask.at<uchar>(i, j) = GC_PR_BGD;
}
else
{
// 當傳回0時為SOURCE
mask.at<uchar>(i, j) = GC_PR_FGD;
}
}
}
}
}
void GrabCut2D::GrabCut( cv::InputArray _img, cv::InputOutputArray _mask, cv::Rect rect, cv::InputOutputArray _bgdModel,cv::InputOutputArray _fgdModel, int iterCount, int mode )
{
std::cout<<"Execute GrabCut Function: Please finish the code here!"<<std::endl;
GMM bgdGmm(_bgdModel.getMatRef());
GMM fgdGmm(_fgdModel.getMatRef());
CheckRectAndMask(_mask.getMatRef(), rect, mode, _img.size());
InitGMMs(_img.getMat(), _mask.getMatRef(), bgdGmm, fgdGmm);
const double gamma = 50;
const double lambda = 9 * gamma;
const double beta = CalBeta(_img.getMat());
Mat leftWeight, upLeftWeight, upWeight, upRightWeight;
CalAllWeights(_img.getMat(), leftWeight, upLeftWeight, upWeight, upRightWeight, beta, gamma);
Mat compIdxs(_img.getMat().size(), CV_32SC1);
for (int i = 0; i < iterCount; i++)
{
int vertexCnt = _img.getMat().rows * _img.getMat().cols;
int edgeCnt = 2 * (4 * vertexCnt - 3 * (_img.getMat().cols + _img.getMat().rows) + 2);
Graph<double, double, double> graph(vertexCnt, edgeCnt);
AssignGMMsComponents(_img.getMat(), _mask.getMat(), bgdGmm, fgdGmm, compIdxs);
LearnGMMs(_img.getMat(), _mask.getMat(), compIdxs, bgdGmm, fgdGmm);
ConstructGraph(_img.getMat(), _mask.getMatRef(), bgdGmm, fgdGmm, lambda, leftWeight, upLeftWeight,
upWeight, upRightWeight, graph);
EstimateSegmentation(graph, _mask.getMatRef());
// cout << _mask.getMatRef() << endl;
}
//一.參數解釋：
//輸入：
//cv::InputArray _img, :輸入的color圖像(類型-cv:Mat)
//cv::Rect rect :在圖像上畫的矩形框（類型-cv:Rect)
//int iterCount : :每次分割的疊代次數（類型-int)
//中間變量
//cv::InputOutputArray _bgdModel ： 背景模型（推薦GMM)（類型-13*n（元件個數）個double類型的自定義資料結構，可以為cv:Mat，或者Vector/List/數組等）
//cv::InputOutputArray _fgdModel : 前景模型（推薦GMM) （類型-13*n（元件個數）個double類型的自定義資料結構，可以為cv:Mat，或者Vector/List/數組等）
//輸出:
//cv::InputOutputArray _mask : 輸出的分割結果 (類型： cv::Mat)
//二. 僞代碼流程：
//1.Load Input Image: 加載輸入顔色圖像;
//2.Init Mask: 用矩形框初始化Mask的Label值（确定背景：0， 确定前景：1，可能背景：2，可能前景：3）,矩形框以外設定為确定背景，矩形框以内設定為可能前景;
//3.Init GMM: 定義并初始化GMM(其他模型完成分割也可得到基本分數，GMM完成會加分）
//4.Sample Points:前背景顔色采樣并進行聚類（建議用kmeans，其他聚類方法也可)
//5.Learn GMM(根據聚類的樣本更新每個GMM元件中的均值、協方差等參數）
//4.Construct Graph（計算t-weight(資料項）和n-weight（平滑項））
//7.Estimate Segmentation(調用maxFlow庫進行分割)
//8.Save Result輸入結果（将結果mask輸出，将mask中前景區域對應的彩色圖像儲存和顯示在互動界面中）
}
#pragma once
#include<opencv2\opencv.hpp>
#include <iostream>
#include "GMM.h"
#include "graph.h"
#include "block.h"
using namespace cv;
using namespace std;
enum
{
GC_WITH_RECT = 0,
GC_WITH_MASK = 1,
GC_CUT = 2
};
class GrabCut2D
{
public:
void GrabCut( cv::InputArray _img, cv::InputOutputArray _mask, cv::Rect rect,
cv::InputOutputArray _bgdModel,cv::InputOutputArray _fgdModel,
int iterCount, int mode );
void CheckRectAndMask(Mat &_mask, Rect rect, int mode, Size imgSize);
void InitMaskWithRect(Mat &_mask, Rect rect, Size imgSize);
void InitGMMs(Mat& img, Mat& mask, GMM& bgdGMM, GMM& fgdGMM);
double CalBeta(const Mat& img);
void CalAllWeights(const Mat& img, Mat& leftWeight, Mat& upLeftWeight, Mat& upWeight, Mat& upRightWeight,
double beta, double gamma);
void AssignGMMsComponents(const Mat& img, const Mat& mask, GMM& bgdGMM, GMM& fgdGMM,
Mat& compIdxs);
void LearnGMMs(const Mat& img, const Mat& mask, const Mat& compIdxs, GMM& bgdGMM, GMM& fgdGMM);
void ConstructGraph(const Mat& img, const Mat& mask, GMM& bgdGMM, GMM& fgdGMM,
double lambda, const Mat& leftWeight, const Mat& upLeftWeight, const Mat& upWeight,
const Mat upRightWeight, Graph<double, double, double>& graph);
void EstimateSegmentation(Graph<double, double, double>& graph, Mat& mask);
~GrabCut2D(void);
};
#pragma once
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "GrabCut.h"
#include <iostream>
#include "BorderMatting.h"
using namespace std;
using namespace cv;
const Scalar BLUE = Scalar(255,0,0); // Background
const Scalar GREEN = Scalar(0,255,0);//Foreground
const Scalar LIGHTBLUE = Scalar(255,255,160);//ProbBackground
const Scalar PINK = Scalar(230,130,255); //ProbBackground
const Scalar RED = Scalar(0,0,255);//color of Rectangle
const int BGD_KEY = CV_EVENT_FLAG_CTRLKEY;// When press "CTRL" key,the value of flags return.
const int FGD_KEY = CV_EVENT_FLAG_SHIFTKEY;// When press "SHIFT" key, the value of flags return.
//Copy the value of comMask to binMask
static void getBinMask( const Mat& comMask, Mat& binMask )
{
if( comMask.empty() || comMask.type()!=CV_8UC1 )
CV_Error( CV_StsBadArg, "comMask is empty or has incorrect type (not CV_8UC1)" );
if( binMask.empty() || binMask.rows!=comMask.rows || binMask.cols!=comMask.cols )
binMask.create( comMask.size(), CV_8UC1 );
binMask = comMask & 1;
}
class GCApplication
{
public:
enum{ NOT_SET = 0, IN_PROCESS = 1, SET = 2 };
static const int radius = 2;
static const int thickness = -1;
void reset();
void setImageAndWinName( const Mat& _image, const string& _winName );
void showImage() const;
void mouseClick( int event, int x, int y, int flags, void* param );
int nextIter();
int getIterCount() const { return iterCount; }
void borderMatting();
private:
void setRectInMask();
void setLblsInMask( int flags, Point p, bool isPr );
const string* winName;
const Mat* image;
Mat mask, alphaMask;
Mat bgdModel, fgdModel;
uchar rectState, lblsState, prLblsState;
bool isInitialized;
Rect rect;
vector<Point> fgdPxls, bgdPxls, prFgdPxls, prBgdPxls;
int iterCount;
GrabCut2D gc;
BorderMatting bm;
};
#include "GCApplication.h"
//Set value for the class
void GCApplication::reset()
{
if( !mask.empty() )
mask.setTo(Scalar::all(GC_BGD));
bgdPxls.clear(); fgdPxls.clear();
prBgdPxls.clear(); prFgdPxls.clear();
isInitialized = false;
rectState = NOT_SET;
lblsState = NOT_SET;
prLblsState = NOT_SET;
iterCount = 0;
}
//Set image and window name
void GCApplication::setImageAndWinName( const Mat& _image, const string& _winName )
{
if( _image.empty() || _winName.empty() )
return;
image = &_image;
winName = &_winName;
mask.create( image->size(), CV_8UC1);
reset();
}
//Show the result image
void GCApplication::showImage() const
{
if( image->empty() || winName->empty() )
return;
Mat res;
Mat binMask;
if( !isInitialized )
image->copyTo( res );
else
{
getBinMask( mask, binMask );
image->copyTo( res, binMask ); //show the GrabCuted image
}
vector<Point>::const_iterator it;
//Using four different colors show the point which have been selected
for( it = bgdPxls.begin(); it != bgdPxls.end(); ++it )
circle( res, *it, radius, BLUE, thickness );
for( it = fgdPxls.begin(); it != fgdPxls.end(); ++it )
circle( res, *it, radius, GREEN, thickness );
for( it = prBgdPxls.begin(); it != prBgdPxls.end(); ++it )
circle( res, *it, radius, LIGHTBLUE, thickness );
for( it = prFgdPxls.begin(); it != prFgdPxls.end(); ++it )
circle( res, *it, radius, PINK, thickness );
imwrite("result.jpg", res);
//Draw the rectangle
if( rectState == IN_PROCESS || rectState == SET )
rectangle( res, Point( rect.x, rect.y ), Point(rect.x + rect.width, rect.y + rect.height ), RED, 2);
imshow( *winName, res );
}
//Using rect initialize the pixel
void GCApplication::setRectInMask()
{
assert( !mask.empty() );
mask.setTo( GC_BGD ); //GC_BGD == 0
rect.x = max(0, rect.x);
rect.y = max(0, rect.y);
rect.width = min(rect.width, image->cols-rect.x);
rect.height = min(rect.height, image->rows-rect.y);
(mask(rect)).setTo( Scalar(GC_PR_FGD) ); //GC_PR_FGD == 3
}
void GCApplication::setLblsInMask( int flags, Point p, bool isPr )
{
vector<Point> *bpxls, *fpxls;
uchar bvalue, fvalue;
if( !isPr ) //Points which are sure being FGD or BGD
{
bpxls = &bgdPxls;
fpxls = &fgdPxls;
bvalue = GC_BGD; //0
fvalue = GC_FGD; //1
}
else //Probably FGD or Probably BGD
{
bpxls = &prBgdPxls;
fpxls = &prFgdPxls;
bvalue = GC_PR_BGD; //2
fvalue = GC_PR_FGD; //3
}
if( flags & BGD_KEY )
{
bpxls->push_back(p);
circle( mask, p, radius, bvalue, thickness ); //Set point value = 2
}
if( flags & FGD_KEY )
{
fpxls->push_back(p);
circle( mask, p, radius, fvalue, thickness ); //Set point value = 3
}
}
//Mouse Click Function: flags work with CV_EVENT_FLAG
void GCApplication::mouseClick( int event, int x, int y, int flags, void* )
{
switch( event )
{
case CV_EVENT_LBUTTONDOWN: // Set rect or GC_BGD(GC_FGD) labels
{
bool isb = (flags & BGD_KEY) != 0,
isf = (flags & FGD_KEY) != 0;
if( rectState == NOT_SET && !isb && !isf )//Only LEFT_KEY pressed
{
rectState = IN_PROCESS; //Be drawing the rectangle
rect = Rect( x, y, 1, 1 );
}
if ( (isb || isf) && rectState == SET ) //Set the BGD/FGD(labels).after press the "ALT" key or "SHIFT" key,and have finish drawing the rectangle
lblsState = IN_PROCESS;
}
break;
case CV_EVENT_RBUTTONDOWN: // Set GC_PR_BGD(GC_PR_FGD) labels
{
bool isb = (flags & BGD_KEY) != 0,
isf = (flags & FGD_KEY) != 0;
if ( (isb || isf) && rectState == SET ) //Set the probably FGD/BGD labels
prLblsState = IN_PROCESS;
}
break;
case CV_EVENT_LBUTTONUP:
if( rectState == IN_PROCESS )
{
rect = Rect( Point(rect.x, rect.y), Point(x,y) ); //After draw the rectangle
rectState = SET;
setRectInMask();
assert( bgdPxls.empty() && fgdPxls.empty() && prBgdPxls.empty() && prFgdPxls.empty() );
showImage();
}
if( lblsState == IN_PROCESS )
{
setLblsInMask(flags, Point(x,y), false); // Draw the FGD points
lblsState = SET;
showImage();
}
break;
case CV_EVENT_RBUTTONUP:
if( prLblsState == IN_PROCESS )
{
setLblsInMask(flags, Point(x,y), true); //Draw the BGD points
prLblsState = SET;
showImage();
}
break;
case CV_EVENT_MOUSEMOVE:
if( rectState == IN_PROCESS )
{
rect = Rect( Point(rect.x, rect.y), Point(x,y) );
assert( bgdPxls.empty() && fgdPxls.empty() && prBgdPxls.empty() && prFgdPxls.empty() );
showImage(); //Continue showing image
}
else if( lblsState == IN_PROCESS )
{
setLblsInMask(flags, Point(x,y), false);
showImage();
}
else if( prLblsState == IN_PROCESS )
{
setLblsInMask(flags, Point(x,y), true);
showImage();
}
break;
}
}
void GCApplication::borderMatting() {
bm.BorderMattingInterface(*image, mask, alphaMask);
}
//Execute GrabCut algorithm，and return the iter count.
int GCApplication::nextIter()
{
if( isInitialized )
gc.GrabCut(*image, mask, rect, bgdModel, fgdModel,1,GC_CUT);
else
{
if( rectState != SET )
return iterCount;
if( lblsState == SET || prLblsState == SET )
gc.GrabCut( *image, mask, rect, bgdModel, fgdModel, 1, GC_WITH_MASK );
else
gc.GrabCut(*image, mask, rect, bgdModel, fgdModel,1,GC_WITH_RECT);
isInitialized = true;
}
iterCount++;
bgdPxls.clear(); fgdPxls.clear();
prBgdPxls.clear(); prFgdPxls.clear();
return iterCount;
}
#include <iostream>
#include "GCApplication.h"
static void help()
{
std::cout << "\nThis program demonstrates GrabCut segmentation -- select an object in a region\n"
"and then grabcut will attempt to segment it out.\n"
"Call:\n"
"./grabcut <image_name>\n"
"\nSelect a rectangular area around the object you want to segment\n" <<
"\nHot keys: \n"
"\tESC - quit the program\n"
"\tr - restore the original image\n"
"\tn - next iteration\n"
"\n"
"\tleft mouse button - set rectangle\n"
"\n"
"\tCTRL+left mouse button - set GC_BGD pixels\n"
"\tSHIFT+left mouse button - set CG_FGD pixels\n"
"\n"
"\tCTRL+right mouse button - set GC_PR_BGD pixels\n"
"\tSHIFT+right mouse button - set CG_PR_FGD pixels\n" << endl;
}
GCApplication gcapp;
static void on_mouse( int event, int x, int y, int flags, void* param )
{
gcapp.mouseClick( event, x, y, flags, param );
}
int main()
{
string filename = "img.jpg";
Mat image = imread( filename, 1 );
//resize(image, image, Size(image.cols / 4, image.rows / 4));
if( image.empty() )
{
cout << "\n , couldn't read image filename " << filename << endl;
return 1;
}
help();
const string winName = "image";
cvNamedWindow( winName.c_str(), CV_WINDOW_AUTOSIZE );
cvSetMouseCallback( winName.c_str(), on_mouse, 0 );
gcapp.setImageAndWinName( image, winName );
gcapp.showImage();
for(;;)
{
int c = cvWaitKey(0);
switch( (char) c )
{
case '\x1b':
cout << "Exiting ..." << endl;
goto exit_main;
case 'r':
cout << endl;
gcapp.reset();
gcapp.showImage();
break;
case 'b':
gcapp.borderMatting();
cout << "done" << endl;
break;
case 'n':
int iterCount = gcapp.getIterCount();
cout << "<" << iterCount << "... ";
int newIterCount = gcapp.nextIter();
if( newIterCount > iterCount )
{
gcapp.showImage();
cout << iterCount << ">" << endl;
}
else
cout << "rect must be determined>" << endl;
break;
}
}
exit_main:
cvDestroyWindow( winName.c_str() );
return 0;
return 0;
}
/* graph.h */
/*
This software library implements the maxflow algorithm
described in
"An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision."
Yuri Boykov and Vladimir Kolmogorov.
In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI),
September 2004
This algorithm was developed by Yuri Boykov and Vladimir Kolmogorov
at Siemens Corporate Research. To make it available for public use,
it was later reimplemented by Vladimir Kolmogorov based on open publications.
If you use this software for research purposes, you should cite
the aforementioned paper in any resulting publication.
----------------------------------------------------------------------
REUSING TREES:
Starting with version 3.0, there is a also an option of reusing search
trees from one maxflow computation to the next, as described in
"Efficiently Solving Dynamic Markov Random Fields Using Graph Cuts."
Pushmeet Kohli and Philip H.S. Torr
International Conference on Computer Vision (ICCV), 2005
If you use this option, you should cite
the aforementioned paper in any resulting publication.
*/
/*
For description, license, example usage see README.TXT.
*/
#ifndef __GRAPH_H__
#define __GRAPH_H__
#include <string.h>
#include "block.h"
#include <assert.h>
// NOTE: in UNIX you need to use -DNDEBUG preprocessor option to supress assert's!!!
// captype: type of edge capacities (excluding t-links)
// tcaptype: type of t-links (edges between nodes and terminals)
// flowtype: type of total flow
//
// Current instantiations are in instances.inc
template <typename captype, typename tcaptype, typename flowtype> class Graph
{
public:
typedef enum
{
SOURCE = 0,
SINK = 1
} termtype; // terminals
typedef int node_id;
/
// BASIC INTERFACE FUNCTIONS //
// (should be enough for most applications) //
/
// Constructor.
// The first argument gives an estimate of the maximum number of nodes that can be added
// to the graph, and the second argument is an estimate of the maximum number of edges.
// The last (optional) argument is the pointer to the function which will be called
// if an error occurs; an error message is passed to this function.
// If this argument is omitted, exit(1) will be called.
//
// IMPORTANT: It is possible to add more nodes to the graph than node_num_max
// (and node_num_max can be zero). However, if the count is exceeded, then
// the internal memory is reallocated (increased by 50%) which is expensive.
// Also, temporarily the amount of allocated memory would be more than twice than needed.
// Similarly for edges.
// If you wish to avoid this overhead, you can download version 2.2, where nodes and edges are stored in blocks.
Graph(int node_num_max, int edge_num_max, void (*err_function)(char *) = NULL);
// Destructor
~Graph();
// Adds node(s) to the graph. By default, one node is added (num=1); then first call returns 0, second call returns 1, and so on.
// If num>1, then several nodes are added, and node_id of the first one is returned.
// IMPORTANT: see note about the constructor
node_id add_node(int num = 1);
// Adds a bidirectional edge between 'i' and 'j' with the weights 'cap' and 'rev_cap'.
// IMPORTANT: see note about the constructor
void add_edge(node_id i, node_id j, captype cap, captype rev_cap);
// Adds new edges 'SOURCE->i' and 'i->SINK' with corresponding weights.
// Can be called multiple times for each node.
// Weights can be negative.
// NOTE: the number of such edges is not counted in edge_num_max.
// No internal memory is allocated by this call.
void add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink);
// Computes the maxflow. Can be called several times.
// FOR DESCRIPTION OF reuse_trees, SEE mark_node().
// FOR DESCRIPTION OF changed_list, SEE remove_from_changed_list().
flowtype maxflow(bool reuse_trees = false, Block<node_id>* changed_list = NULL);
// After the maxflow is computed, this function returns to which
// segment the node 'i' belongs (Graph<captype,tcaptype,flowtype>::SOURCE or Graph<captype,tcaptype,flowtype>::SINK).
//
// Occasionally there may be several minimum cuts. If a node can be assigned
// to both the source and the sink, then default_segm is returned.
termtype what_segment(node_id i, termtype default_segm = SOURCE);
//
// ADVANCED INTERFACE FUNCTIONS //
// (provide access to the graph) //
//
private:
struct node;
struct arc;
public:
// 1. Reallocating graph. //
// Removes all nodes and edges.
// After that functions add_node() and add_edge() must be called again.
//
// Advantage compared to deleting Graph and allocating it again:
// no calls to delete/new (which could be quite slow).
//
// If the graph structure stays the same, then an alternative
// is to go through all nodes/edges and set new residual capacities
// (see functions below).
void reset();
// 2. Functions for getting pointers to arcs and for reading graph structure. //
// NOTE: adding new arcs may invalidate these pointers (if reallocation //
// happens). So it's best not to add arcs while reading graph structure. //
// The following two functions return arcs in the same order that they
// were added to the graph. NOTE: for each call add_edge(i,j,cap,cap_rev)
// the first arc returned will be i->j, and the second j->i.
// If there are no more arcs, then the function can still be called, but
// the returned arc_id is undetermined.
typedef arc* arc_id;
arc_id get_first_arc();
arc_id get_next_arc(arc_id a);
// other functions for reading graph structure
int get_node_num() { return node_num; }
int get_arc_num() { return (int)(arc_last - arcs); }
void get_arc_ends(arc_id a, node_id& i, node_id& j); // returns i,j to that a = i->j
///
// 3. Functions for reading residual capacities. //
///
// returns residual capacity of SOURCE->i minus residual capacity of i->SINK
tcaptype get_trcap(node_id i);
// returns residual capacity of arc a
captype get_rcap(arc* a);
/
// 4. Functions for setting residual capacities. //
// NOTE: If these functions are used, the value of the flow //
// returned by maxflow() will not be valid! //
/
void set_trcap(node_id i, tcaptype trcap);
void set_rcap(arc* a, captype rcap);
// 5. Functions related to reusing trees & list of changed nodes. //
// If flag reuse_trees is true while calling maxflow(), then search trees
// are reused from previous maxflow computation.
// In this case before calling maxflow() the user must
// specify which parts of the graph have changed by calling mark_node():
// add_tweights(i),set_trcap(i) => call mark_node(i)
// add_edge(i,j),set_rcap(a) => call mark_node(i); mark_node(j)
//
// This option makes sense only if a small part of the graph is changed.
// The initialization procedure goes only through marked nodes then.
//
// mark_node(i) can either be called before or after graph modification.
// Can be called more than once per node, but calls after the first one
// do not have any effect.
//
// NOTE:
// - This option cannot be used in the first call to maxflow().
// - It is not necessary to call mark_node() if the change is ``not essential'',
// i.e. sign(trcap) is preserved for a node and zero/nonzero status is preserved for an arc.
// - To check that you marked all necessary nodes, you can call maxflow(false) after calling maxflow(true).
// If everything is correct, the two calls must return the same value of flow. (Useful for debugging).
void mark_node(node_id i);
// If changed_list is not NULL while calling maxflow(), then the algorithm
// keeps a list of nodes which could potentially have changed their segmentation label.
// Nodes which are not in the list are guaranteed to keep their old segmentation label (SOURCE or SINK).
// Example usage:
//
// typedef Graph<int,int,int> G;
// G* g = new Graph(nodeNum, edgeNum);
// Block<G::node_id>* changed_list = new Block<G::node_id>(128);
//
// ... // add nodes and edges
//
// g->maxflow(); // first call should be without arguments
// for (int iter=0; iter<10; iter++)
// {
// ... // change graph, call mark_node() accordingly
//
// g->maxflow(true, changed_list);
// G::node_id* ptr;
// for (ptr=changed_list->ScanFirst(); ptr; ptr=changed_list->ScanNext())
// {
// G::node_id i = *ptr; assert(i>=0 && i<nodeNum);
// g->remove_from_changed_list(i);
// // do something with node i...
// if (g->what_segment(i) == G::SOURCE) { ... }
// }
// changed_list->Reset();
// }
// delete changed_list;
//
// NOTE:
// - If changed_list option is used, then reuse_trees must be used as well.
// - In the example above, the user may omit calls g->remove_from_changed_list(i) and changed_list->Reset() in a given iteration.
// Then during the next call to maxflow(true, &changed_list) new nodes will be added to changed_list.
// - If the next call to maxflow() does not use option reuse_trees, then calling remove_from_changed_list()
// is not necessary. ("changed_list->Reset()" or "delete changed_list" should still be called, though).
void remove_from_changed_list(node_id i)
{
assert(i>=0 && i<node_num && nodes[i].is_in_changed_list);
nodes[i].is_in_changed_list = 0;
}
/
/
/
private:
// internal variables and functions
struct node
{
arc *first; // first outcoming arc
arc *parent; // node's parent
node *next; // pointer to the next active node
// (or to itself if it is the last node in the list)
int TS; // timestamp showing when DIST was computed
int DIST; // distance to the terminal
int is_sink : 1; // flag showing whether the node is in the source or in the sink tree (if parent!=NULL)
int is_marked : 1; // set by mark_node()
int is_in_changed_list : 1; // set by maxflow if
tcaptype tr_cap; // if tr_cap > 0 then tr_cap is residual capacity of the arc SOURCE->node
// otherwise -tr_cap is residual capacity of the arc node->SINK
};
struct arc
{
node *head; // node the arc points to
arc *next; // next arc with the same originating node
arc *sister; // reverse arc
captype r_cap; // residual capacity
};
struct nodeptr
{
node *ptr;
nodeptr *next;
};
static const int NODEPTR_BLOCK_SIZE = 128;
node *nodes, *node_last, *node_max; // node_last = nodes+node_num, node_max = nodes+node_num_max;
arc *arcs, *arc_last, *arc_max; // arc_last = arcs+2*edge_num, arc_max = arcs+2*edge_num_max;
int node_num;
DBlock<nodeptr> *nodeptr_block;
void (*error_function)(char *); // this function is called if a error occurs,
// with a corresponding error message
// (or exit(1) is called if it's NULL)
flowtype flow; // total flow
// reusing trees & list of changed pixels
int maxflow_iteration; // counter
Block<node_id> *changed_list;
/
node *queue_first[2], *queue_last[2]; // list of active nodes
nodeptr *orphan_first, *orphan_last; // list of pointers to orphans
int TIME; // monotonically increasing global counter
/
void reallocate_nodes(int num); // num is the number of new nodes
void reallocate_arcs();
// functions for processing active list
void set_active(node *i);
node *next_active();
// functions for processing orphans list
void set_orphan_front(node* i); // add to the beginning of the list
void set_orphan_rear(node* i); // add to the end of the list
void add_to_changed_list(node* i);
void maxflow_init(); // called if reuse_trees == false
void maxflow_reuse_trees_init(); // called if reuse_trees == true
void augment(arc *middle_arc);
void process_source_orphan(node *i);
void process_sink_orphan(node *i);
void test_consistency(node* current_node=NULL); // debug function
};
///
// Implementation - inline functions //
///
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::node_id Graph<captype,tcaptype,flowtype>::add_node(int num)
{
assert(num > 0);
if (node_last + num > node_max) reallocate_nodes(num);
if (num == 1)
{
node_last -> first = NULL;
node_last -> tr_cap = 0;
node_last -> is_marked = 0;
node_last -> is_in_changed_list = 0;
node_last ++;
return node_num ++;
}
else
{
memset(node_last, 0, num*sizeof(node));
node_id i = node_num;
node_num += num;
node_last += num;
return i;
}
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink)
{
assert(i >= 0 && i < node_num);
tcaptype delta = nodes[i].tr_cap;
if (delta > 0) cap_source += delta;
else cap_sink -= delta;
flow += (cap_source < cap_sink) ? cap_source : cap_sink;
nodes[i].tr_cap = cap_source - cap_sink;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::add_edge(node_id _i, node_id _j, captype cap, captype rev_cap)
{
assert(_i >= 0 && _i < node_num);
assert(_j >= 0 && _j < node_num);
assert(_i != _j);
assert(cap >= 0);
assert(rev_cap >= 0);
if (arc_last == arc_max) reallocate_arcs();
arc *a = arc_last ++;
arc *a_rev = arc_last ++;
node* i = nodes + _i;
node* j = nodes + _j;
a -> sister = a_rev;
a_rev -> sister = a;
a -> next = i -> first;
i -> first = a;
a_rev -> next = j -> first;
j -> first = a_rev;
a -> head = j;
a_rev -> head = i;
a -> r_cap = cap;
a_rev -> r_cap = rev_cap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::arc* Graph<captype,tcaptype,flowtype>::get_first_arc()
{
return arcs;
}
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::arc* Graph<captype,tcaptype,flowtype>::get_next_arc(arc* a)
{
return a + 1;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::get_arc_ends(arc* a, node_id& i, node_id& j)
{
assert(a >= arcs && a < arc_last);
i = (node_id) (a->sister->head - nodes);
j = (node_id) (a->head - nodes);
}
template <typename captype, typename tcaptype, typename flowtype>
inline tcaptype Graph<captype,tcaptype,flowtype>::get_trcap(node_id i)
{
assert(i>=0 && i<node_num);
return nodes[i].tr_cap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline captype Graph<captype,tcaptype,flowtype>::get_rcap(arc* a)
{
assert(a >= arcs && a < arc_last);
return a->r_cap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_trcap(node_id i, tcaptype trcap)
{
assert(i>=0 && i<node_num);
nodes[i].tr_cap = trcap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_rcap(arc* a, captype rcap)
{
assert(a >= arcs && a < arc_last);
a->r_cap = rcap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::termtype Graph<captype,tcaptype,flowtype>::what_segment(node_id i, termtype default_segm)
{
if (nodes[i].parent)
{
return (nodes[i].is_sink) ? SINK : SOURCE;
}
else
{
return default_segm;
}
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::mark_node(node_id _i)
{
node* i = nodes + _i;
if (!i->next)
{
/* it's not in the list yet */
if (queue_last[1]) queue_last[1] -> next = i;
else queue_first[1] = i;
queue_last[1] = i;
i -> next = i;
}
i->is_marked = 1;
}
#endif
/* graph.cpp */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "graph.h"
template <typename captype, typename tcaptype, typename flowtype>
Graph<captype, tcaptype, flowtype>::Graph(int node_num_max, int edge_num_max, void (*err_function)(char *))
: node_num(0),
nodeptr_block(NULL),
error_function(err_function)
{
if (node_num_max < 16) node_num_max = 16;
if (edge_num_max < 16) edge_num_max = 16;
nodes = (node*) malloc(node_num_max*sizeof(node));
arcs = (arc*) malloc(2*edge_num_max*sizeof(arc));
if (!nodes || !arcs) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
node_last = nodes;
node_max = nodes + node_num_max;
arc_last = arcs;
arc_max = arcs + 2*edge_num_max;
maxflow_iteration = 0;
flow = 0;
}
template <typename captype, typename tcaptype, typename flowtype>
Graph<captype,tcaptype,flowtype>::~Graph()
{
if (nodeptr_block)
{
delete nodeptr_block;
nodeptr_block = NULL;
}
free(nodes);
free(arcs);
}
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::reset()
{
node_last = nodes;
arc_last = arcs;
node_num = 0;
if (nodeptr_block)
{
delete nodeptr_block;
nodeptr_block = NULL;
}
maxflow_iteration = 0;
flow = 0;
}
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::reallocate_nodes(int num)
{
int node_num_max = (int)(node_max - nodes);
node* nodes_old = nodes;
node_num_max += node_num_max / 2;
if (node_num_max < node_num + num) node_num_max = node_num + num;
nodes = (node*) realloc(nodes_old, node_num_max*sizeof(node));
if (!nodes) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
node_last = nodes + node_num;
node_max = nodes + node_num_max;
if (nodes != nodes_old)
{
arc* a;
for (a=arcs; a<arc_last; a++)
{
a->head = (node*) ((char*)a->head + (((char*) nodes) - ((char*) nodes_old)));
}
}
}
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::reallocate_arcs()
{
int arc_num_max = (int)(arc_max - arcs);
int arc_num = (int)(arc_last - arcs);
arc* arcs_old = arcs;
arc_num_max += arc_num_max / 2; if (arc_num_max & 1) arc_num_max ++;
arcs = (arc*) realloc(arcs_old, arc_num_max*sizeof(arc));
if (!arcs) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
arc_last = arcs + arc_num;
arc_max = arcs + arc_num_max;
if (arcs != arcs_old)
{
node* i;
arc* a;
for (i=nodes; i<node_last; i++)
{
if (i->first) i->first = (arc*) ((char*)i->first + (((char*) arcs) - ((char*) arcs_old)));
}
for (a=arcs; a<arc_last; a++)
{
if (a->next) a->next = (arc*) ((char*)a->next + (((char*) arcs) - ((char*) arcs_old)));
a->sister = (arc*) ((char*)a->sister + (((char*) arcs) - ((char*) arcs_old)));
}
}
}
#include "instances.inc"
/* maxflow.cpp */
#include <stdio.h>
#include "graph.h"
/*
special constants for node->parent
*/
#define TERMINAL ( (arc *) 1 ) /* to terminal */
#define ORPHAN ( (arc *) 2 ) /* orphan */
#define INFINITE_D ((int)(((unsigned)-1)/2)) /* infinite distance to the terminal */
/***********************************************************************/
/*
Functions for processing active list.
i->next points to the next node in the list
(or to i, if i is the last node in the list).
If i->next is NULL iff i is not in the list.
There are two queues. Active nodes are added
to the end of the second queue and read from
the front of the first queue. If the first queue
is empty, it is replaced by the second queue
(and the second queue becomes empty).
*/
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_active(node *i)
{
if (!i->next)
{
/* it's not in the list yet */
if (queue_last[1]) queue_last[1] -> next = i;
else queue_first[1] = i;
queue_last[1] = i;
i -> next = i;
}
}
/*
Returns the next active node.
If it is connected to the sink, it stays in the list,
otherwise it is removed from the list
*/
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::node* Graph<captype,tcaptype,flowtype>::next_active()
{
node *i;
while ( 1 )
{
if (!(i=queue_first[0]))
{
queue_first[0] = i = queue_first[1];
queue_last[0] = queue_last[1];
queue_first[1] = NULL;
queue_last[1] = NULL;
if (!i) return NULL;
}
/* remove it from the active list */
if (i->next == i) queue_first[0] = queue_last[0] = NULL;
else queue_first[0] = i -> next;
i -> next = NULL;
/* a node in the list is active iff it has a parent */
if (i->parent) return i;
}
}
/***********************************************************************/
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_orphan_front(node *i)
{
nodeptr *np;
i -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = i;
np -> next = orphan_first;
orphan_first = np;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_orphan_rear(node *i)
{
nodeptr *np;
i -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = i;
if (orphan_last) orphan_last -> next = np;
else orphan_first = np;
orphan_last = np;
np -> next = NULL;
}
/***********************************************************************/
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::add_to_changed_list(node *i)
{
if (changed_list && !i->is_in_changed_list)
{
node_id* ptr = changed_list->New();
*ptr = (node_id)(i - nodes);
i->is_in_changed_list = true;
}
}
/***********************************************************************/
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::maxflow_init()
{
node *i;
queue_first[0] = queue_last[0] = NULL;
queue_first[1] = queue_last[1] = NULL;
orphan_first = NULL;
TIME = 0;
for (i=nodes; i<node_last; i++)
{
i -> next = NULL;
i -> is_marked = 0;
i -> is_in_changed_list = 0;
i -> TS = TIME;
if (i->tr_cap > 0)
{
/* i is connected to the source */
i -> is_sink = 0;
i -> parent = TERMINAL;
set_active(i);
i -> DIST = 1;
}
else if (i->tr_cap < 0)
{
/* i is connected to the sink */
i -> is_sink = 1;
i -> parent = TERMINAL;
set_active(i);
i -> DIST = 1;
}
else
{
i -> parent = NULL;
}
}
}
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::maxflow_reuse_trees_init()
{
node* i;
node* j;
node* queue = queue_first[1];
arc* a;
nodeptr* np;
queue_first[0] = queue_last[0] = NULL;
queue_first[1] = queue_last[1] = NULL;
orphan_first = orphan_last = NULL;
TIME ++;
while ((i=queue))
{
queue = i->next;
if (queue == i) queue = NULL;
i->next = NULL;
i->is_marked = 0;
set_active(i);
if (i->tr_cap == 0)
{
if (i->parent) set_orphan_rear(i);
continue;
}
if (i->tr_cap > 0)
{
if (!i->parent || i->is_sink)
{
i->is_sink = 0;
for (a=i->first; a; a=a->next)
{
j = a->head;
if (!j->is_marked)
{
if (j->parent == a->sister) set_orphan_rear(j);
if (j->parent && j->is_sink && a->r_cap > 0) set_active(j);
}
}
add_to_changed_list(i);
}
}
else
{
if (!i->parent || !i->is_sink)
{
i->is_sink = 1;
for (a=i->first; a; a=a->next)
{
j = a->head;
if (!j->is_marked)
{
if (j->parent == a->sister) set_orphan_rear(j);
if (j->parent && !j->is_sink && a->sister->r_cap > 0) set_active(j);
}
}
add_to_changed_list(i);
}
}
i->parent = TERMINAL;
i -> TS = TIME;
i -> DIST = 1;
}
//test_consistency();
/* adoption */
while ((np=orphan_first))
{
orphan_first = np -> next;
i = np -> ptr;
nodeptr_block -> Delete(np);
if (!orphan_first) orphan_last = NULL;
if (i->is_sink) process_sink_orphan(i);
else process_source_orphan(i);
}
/* adoption end */
//test_consistency();
}
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::augment(arc *middle_arc)
{
node *i;
arc *a;
tcaptype bottleneck;
/* 1. Finding bottleneck capacity */
/* 1a - the source tree */
bottleneck = middle_arc -> r_cap;
for (i=middle_arc->sister->head; ; i=a->head)
{
a = i -> parent;
if (a == TERMINAL) break;
if (bottleneck > a->sister->r_cap) bottleneck = a -> sister -> r_cap;
}
if (bottleneck > i->tr_cap) bottleneck = i -> tr_cap;
/* 1b - the sink tree */
for (i=middle_arc->head; ; i=a->head)
{
a = i -> parent;
if (a == TERMINAL) break;
if (bottleneck > a->r_cap) bottleneck = a -> r_cap;
}
if (bottleneck > - i->tr_cap) bottleneck = - i -> tr_cap;
/* 2. Augmenting */
/* 2a - the source tree */
middle_arc -> sister -> r_cap += bottleneck;
middle_arc -> r_cap -= bottleneck;
for (i=middle_arc->sister->head; ; i=a->head)
{
a = i -> parent;
if (a == TERMINAL) break;
a -> r_cap += bottleneck;
a -> sister -> r_cap -= bottleneck;
if (!a->sister->r_cap)
{
set_orphan_front(i); // add i to the beginning of the adoption list
}
}
i -> tr_cap -= bottleneck;
if (!i->tr_cap)
{
set_orphan_front(i); // add i to the beginning of the adoption list
}
/* 2b - the sink tree */
for (i=middle_arc->head; ; i=a->head)
{
a = i -> parent;
if (a == TERMINAL) break;
a -> sister -> r_cap += bottleneck;
a -> r_cap -= bottleneck;
if (!a->r_cap)
{
set_orphan_front(i); // add i to the beginning of the adoption list
}
}
i -> tr_cap += bottleneck;
if (!i->tr_cap)
{
set_orphan_front(i); // add i to the beginning of the adoption list
}
flow += bottleneck;
}
/***********************************************************************/
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::process_source_orphan(node *i)
{
node *j;
arc *a0, *a0_min = NULL, *a;
int d, d_min = INFINITE_D;
/* trying to find a new parent */
for (a0=i->first; a0; a0=a0->next)
if (a0->sister->r_cap)
{
j = a0 -> head;
if (!j->is_sink && (a=j->parent))
{
/* checking the origin of j */
d = 0;
while ( 1 )
{
if (j->TS == TIME)
{
d += j -> DIST;
break;
}
a = j -> parent;
d ++;
if (a==TERMINAL)
{
j -> TS = TIME;
j -> DIST = 1;
break;
}
if (a==ORPHAN) { d = INFINITE_D; break; }
j = a -> head;
}
if (d<INFINITE_D) /* j originates from the source - done */
{
if (d<d_min)
{
a0_min = a0;
d_min = d;
}
/* set marks along the path */
for (j=a0->head; j->TS!=TIME; j=j->parent->head)
{
j -> TS = TIME;
j -> DIST = d --;
}
}
}
}
if (i->parent = a0_min)
{
i -> TS = TIME;
i -> DIST = d_min + 1;
}
else
{
/* no parent is found */
add_to_changed_list(i);
/* process neighbors */
for (a0=i->first; a0; a0=a0->next)
{
j = a0 -> head;
if (!j->is_sink && (a=j->parent))
{
if (a0->sister->r_cap) set_active(j);
if (a!=TERMINAL && a!=ORPHAN && a->head==i)
{
set_orphan_rear(j); // add j to the end of the adoption list
}
}
}
}
}
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::process_sink_orphan(node *i)
{
node *j;
arc *a0, *a0_min = NULL, *a;
int d, d_min = INFINITE_D;
/* trying to find a new parent */
for (a0=i->first; a0; a0=a0->next)
if (a0->r_cap)
{
j = a0 -> head;
if (j->is_sink && (a=j->parent))
{
/* checking the origin of j */
d = 0;
while ( 1 )
{
if (j->TS == TIME)
{
d += j -> DIST;
break;
}
a = j -> parent;
d ++;
if (a==TERMINAL)
{
j -> TS = TIME;
j -> DIST = 1;
break;
}
if (a==ORPHAN) { d = INFINITE_D; break; }
j = a -> head;
}
if (d<INFINITE_D) /* j originates from the sink - done */
{
if (d<d_min)
{
a0_min = a0;
d_min = d;
}
/* set marks along the path */
for (j=a0->head; j->TS!=TIME; j=j->parent->head)
{
j -> TS = TIME;
j -> DIST = d --;
}
}
}
}
if (i->parent = a0_min)
{
i -> TS = TIME;
i -> DIST = d_min + 1;
}
else
{
/* no parent is found */
add_to_changed_list(i);
/* process neighbors */
for (a0=i->first; a0; a0=a0->next)
{
j = a0 -> head;
if (j->is_sink && (a=j->parent))
{
if (a0->r_cap) set_active(j);
if (a!=TERMINAL && a!=ORPHAN && a->head==i)
{
set_orphan_rear(j); // add j to the end of the adoption list
}
}
}
}
}
/***********************************************************************/
template <typename captype, typename tcaptype, typename flowtype>
flowtype Graph<captype,tcaptype,flowtype>::maxflow(bool reuse_trees, Block<node_id>* _changed_list)
{
node *i, *j, *current_node = NULL;
arc *a;
nodeptr *np, *np_next;
if (!nodeptr_block)
{
nodeptr_block = new DBlock<nodeptr>(NODEPTR_BLOCK_SIZE, error_function);
}
changed_list = _changed_list;
if (maxflow_iteration == 0 && reuse_trees) { if (error_function) (*error_function)("reuse_trees cannot be used in the first call to maxflow()!"); exit(1); }
if (changed_list && !reuse_trees) { if (error_function) (*error_function)("changed_list cannot be used without reuse_trees!"); exit(1); }
if (reuse_trees) maxflow_reuse_trees_init();
else maxflow_init();
// main loop
while ( 1 )
{
// test_consistency(current_node);
if ((i=current_node))
{
i -> next = NULL; /* remove active flag */
if (!i->parent) i = NULL;
}
if (!i)
{
if (!(i = next_active())) break;
}
/* growth */
if (!i->is_sink)
{
/* grow source tree */
for (a=i->first; a; a=a->next)
if (a->r_cap)
{
j = a -> head;
if (!j->parent)
{
j -> is_sink = 0;
j -> parent = a -> sister;
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
set_active(j);
add_to_changed_list(j);
}
else if (j->is_sink) break;
else if (j->TS <= i->TS &&
j->DIST > i->DIST)
{
/* heuristic - trying to make the distance from j to the source shorter */
j -> parent = a -> sister;
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
}
}
}
else
{
/* grow sink tree */
for (a=i->first; a; a=a->next)
if (a->sister->r_cap)
{
j = a -> head;
if (!j->parent)
{
j -> is_sink = 1;
j -> parent = a -> sister;
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
set_active(j);
add_to_changed_list(j);
}
else if (!j->is_sink) { a = a -> sister; break; }
else if (j->TS <= i->TS &&
j->DIST > i->DIST)
{
/* heuristic - trying to make the distance from j to the sink shorter */
j -> parent = a -> sister;
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
}
}
}
TIME ++;
if (a)
{
i -> next = i; /* set active flag */
current_node = i;
/* augmentation */
augment(a);
/* augmentation end */
/* adoption */
while ((np=orphan_first))
{
np_next = np -> next;
np -> next = NULL;
while ((np=orphan_first))
{
orphan_first = np -> next;
i = np -> ptr;
nodeptr_block -> Delete(np);
if (!orphan_first) orphan_last = NULL;
if (i->is_sink) process_sink_orphan(i);
else process_source_orphan(i);
}
orphan_first = np_next;
}
/* adoption end */
}
else current_node = NULL;
}
// test_consistency();
if (!reuse_trees || (maxflow_iteration % 64) == 0)
{
delete nodeptr_block;
nodeptr_block = NULL;
}
maxflow_iteration ++;
return flow;
}
/***********************************************************************/
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::test_consistency(node* current_node)
{
node *i;
arc *a;
int r;
int num1 = 0, num2 = 0;
// test whether all nodes i with i->next!=NULL are indeed in the queue
for (i=nodes; i<node_last; i++)
{
if (i->next || i==current_node) num1 ++;
}
for (r=0; r<3; r++)
{
i = (r == 2) ? current_node : queue_first[r];
if (i)
for ( ; ; i=i->next)
{
num2 ++;
if (i->next == i)
{
if (r<2) assert(i == queue_last[r]);
else assert(i == current_node);
break;
}
}
}
assert(num1 == num2);
for (i=nodes; i<node_last; i++)
{
// test whether all edges in seach trees are non-saturated
if (i->parent == NULL) {}
else if (i->parent == ORPHAN) {}
else if (i->parent == TERMINAL)
{
if (!i->is_sink) assert(i->tr_cap > 0);
else assert(i->tr_cap < 0);
}
else
{
if (!i->is_sink) assert (i->parent->sister->r_cap > 0);
else assert (i->parent->r_cap > 0);
}
// test whether passive nodes in search trees have neighbors in
// a different tree through non-saturated edges
if (i->parent && !i->next)
{
if (!i->is_sink)
{
assert(i->tr_cap >= 0);
for (a=i->first; a; a=a->next)
{
if (a->r_cap > 0) assert(a->head->parent && !a->head->is_sink);
}
}
else
{
assert(i->tr_cap <= 0);
for (a=i->first; a; a=a->next)
{
if (a->sister->r_cap > 0) assert(a->head->parent && a->head->is_sink);
}
}
}
// test marking invariants
if (i->parent && i->parent!=ORPHAN && i->parent!=TERMINAL)
{
assert(i->TS <= i->parent->head->TS);
if (i->TS == i->parent->head->TS) assert(i->DIST > i->parent->head->DIST);
}
}
}
#include "instances.inc"
/* block.h */
/*
Template classes Block and DBlock
Implement adding and deleting items of the same type in blocks.
If there there are many items then using Block or DBlock
is more efficient than using 'new' and 'delete' both in terms
of memory and time since
(1) On some systems there is some minimum amount of memory
that 'new' can allocate (e.g., 64), so if items are
small that a lot of memory is wasted.
(2) 'new' and 'delete' are designed for items of varying size.
If all items has the same size, then an algorithm for
adding and deleting can be made more efficient.
(3) All Block and DBlock functions are inline, so there are
no extra function calls.
Differences between Block and DBlock:
(1) DBlock allows both adding and deleting items,
whereas Block allows only adding items.
(2) Block has an additional operation of scanning
items added so far (in the order in which they were added).
(3) Block allows to allocate several consecutive
items at a time, whereas DBlock can add only a single item.
Note that no constructors or destructors are called for items.
Example usage for items of type 'MyType':
///
#include "block.h"
#define BLOCK_SIZE 1024
typedef struct { int a, b; } MyType;
MyType *ptr, *array[10000];
...
Block<MyType> *block = new Block<MyType>(BLOCK_SIZE);
// adding items
for (int i=0; i<sizeof(array); i++)
{
ptr = block -> New();
ptr -> a = ptr -> b = rand();
}
// reading items
for (ptr=block->ScanFirst(); ptr; ptr=block->ScanNext())
{
printf("%d %d\n", ptr->a, ptr->b);
}
delete block;
...
DBlock<MyType> *dblock = new DBlock<MyType>(BLOCK_SIZE);
// adding items
for (int i=0; i<sizeof(array); i++)
{
array[i] = dblock -> New();
}
// deleting items
for (int i=0; i<sizeof(array); i+=2)
{
dblock -> Delete(array[i]);
}
// adding items
for (int i=0; i<sizeof(array); i++)
{
array[i] = dblock -> New();
}
delete dblock;
///
Note that DBlock deletes items by marking them as
empty (i.e., by adding them to the list of free items),
so that this memory could be used for subsequently
added items. Thus, at each moment the memory allocated
is determined by the maximum number of items allocated
simultaneously at earlier moments. All memory is
deallocated only when the destructor is called.
*/
#ifndef __BLOCK_H__
#define __BLOCK_H__
#include <stdlib.h>
/***********************************************************************/
/***********************************************************************/
/***********************************************************************/
template <class Type> class Block
{
public:
/* Constructor. Arguments are the block size and
(optionally) the pointer to the function which
will be called if allocation failed; the message
passed to this function is "Not enough memory!" */
Block(int size, void (*err_function)(char *) = NULL) { first = last = NULL; block_size = size; error_function = err_function; }
/* Destructor. Deallocates all items added so far */
~Block() { while (first) { block *next = first -> next; delete[] ((char*)first); first = next; } }
/* Allocates 'num' consecutive items; returns pointer
to the first item. 'num' cannot be greater than the
block size since items must fit in one block */
Type *New(int num = 1)
{
Type *t;
if (!last || last->current + num > last->last)
{
if (last && last->next) last = last -> next;
else
{
block *next = (block *) new char [sizeof(block) + (block_size-1)*sizeof(Type)];
if (!next) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
if (last) last -> next = next;
else first = next;
last = next;
last -> current = & ( last -> data[0] );
last -> last = last -> current + block_size;
last -> next = NULL;
}
}
t = last -> current;
last -> current += num;
return t;
}
/* Returns the first item (or NULL, if no items were added) */
Type *ScanFirst()
{
for (scan_current_block=first; scan_current_block; scan_current_block = scan_current_block->next)
{
scan_current_data = & ( scan_current_block -> data[0] );
if (scan_current_data < scan_current_block -> current) return scan_current_data ++;
}
return NULL;
}
/* Returns the next item (or NULL, if all items have been read)
Can be called only if previous ScanFirst() or ScanNext()
call returned not NULL. */
Type *ScanNext()
{
while (scan_current_data >= scan_current_block -> current)
{
scan_current_block = scan_current_block -> next;
if (!scan_current_block) return NULL;
scan_current_data = & ( scan_current_block -> data[0] );
}
return scan_current_data ++;
}
/* Marks all elements as empty */
void Reset()
{
block *b;
if (!first) return;
for (b=first; ; b=b->next)
{
b -> current = & ( b -> data[0] );
if (b == last) break;
}
last = first;
}
/***********************************************************************/
private:
typedef struct block_st
{
Type *current, *last;
struct block_st *next;
Type data[1];
} block;
int block_size;
block *first;
block *last;
block *scan_current_block;
Type *scan_current_data;
void (*error_function)(char *);
};
/***********************************************************************/
/***********************************************************************/
/***********************************************************************/
template <class Type> class DBlock
{
public:
/* Constructor. Arguments are the block size and
(optionally) the pointer to the function which
will be called if allocation failed; the message
passed to this function is "Not enough memory!" */
DBlock(int size, void (*err_function)(char *) = NULL) { first = NULL; first_free = NULL; block_size = size; error_function = err_function; }
/* Destructor. Deallocates all items added so far */
~DBlock() { while (first) { block *next = first -> next; delete[] ((char*)first); first = next; } }
/* Allocates one item */
Type *New()
{
block_item *item;
if (!first_free)
{
block *next = first;
first = (block *) new char [sizeof(block) + (block_size-1)*sizeof(block_item)];
if (!first) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
first_free = & (first -> data[0] );
for (item=first_free; item<first_free+block_size-1; item++)
item -> next_free = item + 1;
item -> next_free = NULL;
first -> next = next;
}
item = first_free;
first_free = item -> next_free;
return (Type *) item;
}
/* Deletes an item allocated previously */
void Delete(Type *t)
{
((block_item *) t) -> next_free = first_free;
first_free = (block_item *) t;
}
/***********************************************************************/
private:
typedef union block_item_st
{
Type t;
block_item_st *next_free;
} block_item;
typedef struct block_st
{
struct block_st *next;
block_item data[1];
} block;
int block_size;
block *first;
block_item *first_free;
void (*error_function)(char *);
};
#endif
/* graph.h */
/*
This software library implements the maxflow algorithm
described in
"An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision."
Yuri Boykov and Vladimir Kolmogorov.
In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI),
September 2004
This algorithm was developed by Yuri Boykov and Vladimir Kolmogorov
at Siemens Corporate Research. To make it available for public use,
it was later reimplemented by Vladimir Kolmogorov based on open publications.
If you use this software for research purposes, you should cite
the aforementioned paper in any resulting publication.
----------------------------------------------------------------------
REUSING TREES:
Starting with version 3.0, there is a also an option of reusing search
trees from one maxflow computation to the next, as described in
"Efficiently Solving Dynamic Markov Random Fields Using Graph Cuts."
Pushmeet Kohli and Philip H.S. Torr
International Conference on Computer Vision (ICCV), 2005
If you use this option, you should cite
the aforementioned paper in any resulting publication.
*/
/*
For description, license, example usage see README.TXT.
*/
#ifndef __GRAPH_H__
#define __GRAPH_H__
#include <string.h>
#include "block.h"
#include <assert.h>
// NOTE: in UNIX you need to use -DNDEBUG preprocessor option to supress assert's!!!
// captype: type of edge capacities (excluding t-links)
// tcaptype: type of t-links (edges between nodes and terminals)
// flowtype: type of total flow
//
// Current instantiations are in instances.inc
template <typename captype, typename tcaptype, typename flowtype> class Graph
{
public:
typedef enum
{
SOURCE = 0,
SINK = 1
} termtype; // terminals
typedef int node_id;
/
// BASIC INTERFACE FUNCTIONS //
// (should be enough for most applications) //
/
// Constructor.
// The first argument gives an estimate of the maximum number of nodes that can be added
// to the graph, and the second argument is an estimate of the maximum number of edges.
// The last (optional) argument is the pointer to the function which will be called
// if an error occurs; an error message is passed to this function.
// If this argument is omitted, exit(1) will be called.
//
// IMPORTANT: It is possible to add more nodes to the graph than node_num_max
// (and node_num_max can be zero). However, if the count is exceeded, then
// the internal memory is reallocated (increased by 50%) which is expensive.
// Also, temporarily the amount of allocated memory would be more than twice than needed.
// Similarly for edges.
// If you wish to avoid this overhead, you can download version 2.2, where nodes and edges are stored in blocks.
Graph(int node_num_max, int edge_num_max, void (*err_function)(char *) = NULL);
// Destructor
~Graph();
// Adds node(s) to the graph. By default, one node is added (num=1); then first call returns 0, second call returns 1, and so on.
// If num>1, then several nodes are added, and node_id of the first one is returned.
// IMPORTANT: see note about the constructor
node_id add_node(int num = 1);
// Adds a bidirectional edge between 'i' and 'j' with the weights 'cap' and 'rev_cap'.
// IMPORTANT: see note about the constructor
void add_edge(node_id i, node_id j, captype cap, captype rev_cap);
// Adds new edges 'SOURCE->i' and 'i->SINK' with corresponding weights.
// Can be called multiple times for each node.
// Weights can be negative.
// NOTE: the number of such edges is not counted in edge_num_max.
// No internal memory is allocated by this call.
void add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink);
// Computes the maxflow. Can be called several times.
// FOR DESCRIPTION OF reuse_trees, SEE mark_node().
// FOR DESCRIPTION OF changed_list, SEE remove_from_changed_list().
flowtype maxflow(bool reuse_trees = false, Block<node_id>* changed_list = NULL);
// After the maxflow is computed, this function returns to which
// segment the node 'i' belongs (Graph<captype,tcaptype,flowtype>::SOURCE or Graph<captype,tcaptype,flowtype>::SINK).
//
// Occasionally there may be several minimum cuts. If a node can be assigned
// to both the source and the sink, then default_segm is returned.
termtype what_segment(node_id i, termtype default_segm = SOURCE);
//
// ADVANCED INTERFACE FUNCTIONS //
// (provide access to the graph) //
//
private:
struct node;
struct arc;
public:
// 1. Reallocating graph. //
// Removes all nodes and edges.
// After that functions add_node() and add_edge() must be called again.
//
// Advantage compared to deleting Graph and allocating it again:
// no calls to delete/new (which could be quite slow).
//
// If the graph structure stays the same, then an alternative
// is to go through all nodes/edges and set new residual capacities
// (see functions below).
void reset();
// 2. Functions for getting pointers to arcs and for reading graph structure. //
// NOTE: adding new arcs may invalidate these pointers (if reallocation //
// happens). So it's best not to add arcs while reading graph structure. //
// The following two functions return arcs in the same order that they
// were added to the graph. NOTE: for each call add_edge(i,j,cap,cap_rev)
// the first arc returned will be i->j, and the second j->i.
// If there are no more arcs, then the function can still be called, but
// the returned arc_id is undetermined.
typedef arc* arc_id;
arc_id get_first_arc();
arc_id get_next_arc(arc_id a);
// other functions for reading graph structure
int get_node_num() { return node_num; }
int get_arc_num() { return (int)(arc_last - arcs); }
void get_arc_ends(arc_id a, node_id& i, node_id& j); // returns i,j to that a = i->j
///
// 3. Functions for reading residual capacities. //
///
// returns residual capacity of SOURCE->i minus residual capacity of i->SINK
tcaptype get_trcap(node_id i);
// returns residual capacity of arc a
captype get_rcap(arc* a);
/
// 4. Functions for setting residual capacities. //
// NOTE: If these functions are used, the value of the flow //
// returned by maxflow() will not be valid! //
/
void set_trcap(node_id i, tcaptype trcap);
void set_rcap(arc* a, captype rcap);
// 5. Functions related to reusing trees & list of changed nodes. //
// If flag reuse_trees is true while calling maxflow(), then search trees
// are reused from previous maxflow computation.
// In this case before calling maxflow() the user must
// specify which parts of the graph have changed by calling mark_node():
// add_tweights(i),set_trcap(i) => call mark_node(i)
// add_edge(i,j),set_rcap(a) => call mark_node(i); mark_node(j)
//
// This option makes sense only if a small part of the graph is changed.
// The initialization procedure goes only through marked nodes then.
//
// mark_node(i) can either be called before or after graph modification.
// Can be called more than once per node, but calls after the first one
// do not have any effect.
//
// NOTE:
// - This option cannot be used in the first call to maxflow().
// - It is not necessary to call mark_node() if the change is ``not essential'',
// i.e. sign(trcap) is preserved for a node and zero/nonzero status is preserved for an arc.
// - To check that you marked all necessary nodes, you can call maxflow(false) after calling maxflow(true).
// If everything is correct, the two calls must return the same value of flow. (Useful for debugging).
void mark_node(node_id i);
// If changed_list is not NULL while calling maxflow(), then the algorithm
// keeps a list of nodes which could potentially have changed their segmentation label.
// Nodes which are not in the list are guaranteed to keep their old segmentation label (SOURCE or SINK).
// Example usage:
//
// typedef Graph<int,int,int> G;
// G* g = new Graph(nodeNum, edgeNum);
// Block<G::node_id>* changed_list = new Block<G::node_id>(128);
//
// ... // add nodes and edges
//
// g->maxflow(); // first call should be without arguments
// for (int iter=0; iter<10; iter++)
// {
// ... // change graph, call mark_node() accordingly
//
// g->maxflow(true, changed_list);
// G::node_id* ptr;
// for (ptr=changed_list->ScanFirst(); ptr; ptr=changed_list->ScanNext())
// {
// G::node_id i = *ptr; assert(i>=0 && i<nodeNum);
// g->remove_from_changed_list(i);
// // do something with node i...
// if (g->what_segment(i) == G::SOURCE) { ... }
// }
// changed_list->Reset();
// }
// delete changed_list;
//
// NOTE:
// - If changed_list option is used, then reuse_trees must be used as well.
// - In the example above, the user may omit calls g->remove_from_changed_list(i) and changed_list->Reset() in a given iteration.
// Then during the next call to maxflow(true, &changed_list) new nodes will be added to changed_list.
// - If the next call to maxflow() does not use option reuse_trees, then calling remove_from_changed_list()
// is not necessary. ("changed_list->Reset()" or "delete changed_list" should still be called, though).
void remove_from_changed_list(node_id i)
{
assert(i>=0 && i<node_num && nodes[i].is_in_changed_list);
nodes[i].is_in_changed_list = 0;
}
/
/
/
private:
// internal variables and functions
struct node
{
arc *first; // first outcoming arc
arc *parent; // node's parent
node *next; // pointer to the next active node
// (or to itself if it is the last node in the list)
int TS; // timestamp showing when DIST was computed
int DIST; // distance to the terminal
int is_sink : 1; // flag showing whether the node is in the source or in the sink tree (if parent!=NULL)
int is_marked : 1; // set by mark_node()
int is_in_changed_list : 1; // set by maxflow if
tcaptype tr_cap; // if tr_cap > 0 then tr_cap is residual capacity of the arc SOURCE->node
// otherwise -tr_cap is residual capacity of the arc node->SINK
};
struct arc
{
node *head; // node the arc points to
arc *next; // next arc with the same originating node
arc *sister; // reverse arc
captype r_cap; // residual capacity
};
struct nodeptr
{
node *ptr;
nodeptr *next;
};
static const int NODEPTR_BLOCK_SIZE = 128;
node *nodes, *node_last, *node_max; // node_last = nodes+node_num, node_max = nodes+node_num_max;
arc *arcs, *arc_last, *arc_max; // arc_last = arcs+2*edge_num, arc_max = arcs+2*edge_num_max;
int node_num;
DBlock<nodeptr> *nodeptr_block;
void (*error_function)(char *); // this function is called if a error occurs,
// with a corresponding error message
// (or exit(1) is called if it's NULL)
flowtype flow; // total flow
// reusing trees & list of changed pixels
int maxflow_iteration; // counter
Block<node_id> *changed_list;
/
node *queue_first[2], *queue_last[2]; // list of active nodes
nodeptr *orphan_first, *orphan_last; // list of pointers to orphans
int TIME; // monotonically increasing global counter
/
void reallocate_nodes(int num); // num is the number of new nodes
void reallocate_arcs();
// functions for processing active list
void set_active(node *i);
node *next_active();
// functions for processing orphans list
void set_orphan_front(node* i); // add to the beginning of the list
void set_orphan_rear(node* i); // add to the end of the list
void add_to_changed_list(node* i);
void maxflow_init(); // called if reuse_trees == false
void maxflow_reuse_trees_init(); // called if reuse_trees == true
void augment(arc *middle_arc);
void process_source_orphan(node *i);
void process_sink_orphan(node *i);
void test_consistency(node* current_node=NULL); // debug function
};
///
// Implementation - inline functions //
///
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::node_id Graph<captype,tcaptype,flowtype>::add_node(int num)
{
assert(num > 0);
if (node_last + num > node_max) reallocate_nodes(num);
if (num == 1)
{
node_last -> first = NULL;
node_last -> tr_cap = 0;
node_last -> is_marked = 0;
node_last -> is_in_changed_list = 0;
node_last ++;
return node_num ++;
}
else
{
memset(node_last, 0, num*sizeof(node));
node_id i = node_num;
node_num += num;
node_last += num;
return i;
}
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink)
{
assert(i >= 0 && i < node_num);
tcaptype delta = nodes[i].tr_cap;
if (delta > 0) cap_source += delta;
else cap_sink -= delta;
flow += (cap_source < cap_sink) ? cap_source : cap_sink;
nodes[i].tr_cap = cap_source - cap_sink;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::add_edge(node_id _i, node_id _j, captype cap, captype rev_cap)
{
assert(_i >= 0 && _i < node_num);
assert(_j >= 0 && _j < node_num);
assert(_i != _j);
assert(cap >= 0);
assert(rev_cap >= 0);
if (arc_last == arc_max) reallocate_arcs();
arc *a = arc_last ++;
arc *a_rev = arc_last ++;
node* i = nodes + _i;
node* j = nodes + _j;
a -> sister = a_rev;
a_rev -> sister = a;
a -> next = i -> first;
i -> first = a;
a_rev -> next = j -> first;
j -> first = a_rev;
a -> head = j;
a_rev -> head = i;
a -> r_cap = cap;
a_rev -> r_cap = rev_cap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::arc* Graph<captype,tcaptype,flowtype>::get_first_arc()
{
return arcs;
}
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::arc* Graph<captype,tcaptype,flowtype>::get_next_arc(arc* a)
{
return a + 1;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::get_arc_ends(arc* a, node_id& i, node_id& j)
{
assert(a >= arcs && a < arc_last);
i = (node_id) (a->sister->head - nodes);
j = (node_id) (a->head - nodes);
}
template <typename captype, typename tcaptype, typename flowtype>
inline tcaptype Graph<captype,tcaptype,flowtype>::get_trcap(node_id i)
{
assert(i>=0 && i<node_num);
return nodes[i].tr_cap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline captype Graph<captype,tcaptype,flowtype>::get_rcap(arc* a)
{
assert(a >= arcs && a < arc_last);
return a->r_cap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_trcap(node_id i, tcaptype trcap)
{
assert(i>=0 && i<node_num);
nodes[i].tr_cap = trcap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_rcap(arc* a, captype rcap)
{
assert(a >= arcs && a < arc_last);
a->r_cap = rcap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::termtype Graph<captype,tcaptype,flowtype>::what_segment(node_id i, termtype default_segm)
{
if (nodes[i].parent)
{
return (nodes[i].is_sink) ? SINK : SOURCE;
}
else
{
return default_segm;
}
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::mark_node(node_id _i)
{
node* i = nodes + _i;
if (!i->next)
{
/* it's not in the list yet */
if (queue_last[1]) queue_last[1] -> next = i;
else queue_first[1] = i;
queue_last[1] = i;
i -> next = i;
}
i->is_marked = 1;
}
#endif
#include "graph.h"
#ifdef _MSC_VER
#pragma warning(disable: 4661)
#endif
// Instantiations: <captype, tcaptype, flowtype>
// IMPORTANT:
// flowtype should be 'larger' than tcaptype
// tcaptype should be 'larger' than captype
template class Graph<int,int,int>;
template class Graph<short,int,int>;
template class Graph<float,float,float>;
template class Graph<double,double,double>;
使用了maxflow庫：http://vision.csd.uwo.ca/code/ 中Max-flow/min-cut