一、序言
陸陸續續的如果累計起來，我估計至少有二十來位左右的朋友加我QQ，向我咨詢有關摳圖方面的算法，可惜的是，我對這方面之前一直是沒有研究過的。除了利用和Photoshop中的魔棒一樣的技術或者Photoshop中的選區菜單中的色彩范圍類似的算法（這兩個我有何PS至少90%一致的代碼）是實現簡單的摳圖外，現在一些state of art 方面的算法我都不了解。因此，也浪費了不少的將知識轉換為資產的機會。年30那天，偶然的一個機會，有位朋友推薦我看了一篇關於摳圖的文章，並有配套的實現代碼，於是我就決定從這篇文章開始我的摳圖算法研究之旅。
這篇文章就是Shared Sampling for Real-Time Alpha Matting，關於這篇文章的一些信息，可以在這個網站里找到很多：http://www.inf.ufrgs.br/~eslgastal/SharedMatting/，配套的一個代碼在CSDN中可以下載,具體見：http://download.csdn.net/detail/jlwyc/4676516
這篇文章的標題很具有吸引力，發表日期為2010，也算是比較新的。在大家繼續看下去之前，我要提醒的是，這里的Real - Time有比較多的限制：主要是（1）必須依賴於強勁的GPU；（2）應用的摳圖場合的背景應該比較簡單。
不管如何，因為有配套的實現代碼，作為起步的研究來說，該文還是算不錯的。
從目前流行的摳圖技術來看，這篇文章的思路算是比較落伍的一種。
二、技術細節
好了，不管那么多，我先貼些論文中的公式及一些說明將文章的主體細路描述一下。
簡單的說，摳圖問題就是要解決如下的一個超級病態的方程：
式中：Cp是我們觀察到的圖像的顏色，FP、BP、αp均是未知量，可分別稱之為前景、背景及透明度。
要解決這樣的一個病態的方程，就必須給其增加一些附加的約束，通常，這種約束可以是和待分割圖像同等大小的TriMap或者是用戶收工划定的scribbles的形式存在，如下兩圖所示(如未特別說明，一般白色部分表示前景，黑色表示背景，灰色表示待識別的部分）:
TriMap scribbles
這樣的約束條件使得我們知道了那一部分是明確屬於前景（αp=1），而那一部分是屬於背景（αp=0），那么下面的主要任務就是搞定那些未知區域的αp值 。
按照論文的說法，在2010年前后解決matting問題的主要方法是基於sampling, pixel affinities 或者兩者的結合，特別是后兩種是主流的方式。但是這兩種都需要求解一個大型的線性系統，這個系統的大小和未知點的個數成正比（我簡單看了下closed form那篇摳圖文檔的代碼，就用到了一個龐大的稀疏矩陣），因此對於1MB左右大小的圖，求解時間在幾秒到幾分鍾不等。這篇論文提出的算法應該說是基於sampling技術的，他充分利用了相鄰像素之間的相似性，並利用了算法內在的並行性，結合GPU編程，實現摳圖的實時展示。
總的來說，論文提出的算法可以分成4個步驟：
第一步：Expansion，針對用戶的輸入，對已知區域（前景或背景）進行小規模的擴展；
第二步：Sample and Gather，對剩余的未知區域內的每個點按一定的規則取樣，並選擇出最佳的一對前景和背景取樣點；
第三步：Refinement，在一定的領域范圍內，對未知區域內的每個點的最佳配對重新進行組合。
第四步：Local Smoothing，對得到的前景和背景對以及透明度值進行局部平滑，以減少噪音。
2.1 Expansion
這一步，按照我的經驗，可以不做，他唯一的作用就是減少未知點的個數，可能在一定程度上減小后期的計算量，原理也很簡單，就是對一個未知點，在其一定的鄰域半徑內（文中推薦值10，
並且是圓形半徑），如果有已知的背景點或前景點，則計算其顏色和這些已知點顏色的距離，然后把這個未知點歸屬於和其顏色距離小於某個值並且最靠近該點的對象（前景或背景）。
在CSDN提供的參考代碼中，這一部分的編碼其實寫的還是很有特色的，他的循環方式不同於我們普通的鄰域編碼，他是從像素點逐漸向外部循環開來，有點類似左圖的這種循環方式（實際上還是有點區別的，實際是上下兩行一起處理，在左右兩列處理，然后再向外層擴散），這種處理方式的明顯好處就是，只要找到某個點顏色距離小於設定的值，就可以停止循環了，因為這個點肯定是第一個符合顏色距離條件又同時符合物理距離最小的要求的。
這一步做不做，最最終的結果又一定的影響，但是他不具有質的影響。
2.2 Sample and Gather
總的來說，這一步是算法的核心部分，也是對結果影響最大的，他的步驟說起來其實也很簡單，我們先看下圖。
在這個圖中，P和q點都是未知區域，我們需要通過一定的原則在已知區域為其取得一定的樣本對，論文中提出的提取方法是：
設定一個參數Kg,其意義為一個點最多可能取樣的前景點和背景點的個數，也就意味着最多的取樣對為Kg*Kg組，通常這個值可以取為4或者更多，論文建議取4就可以了，越大則程序越耗時。
這樣對於每個未知點，從該點出發，引出Kg條路徑，每個路徑之間成360/Kg的夾角，記錄下每條路徑經過的路線中首次遇到的前景或背景點，直到超出圖像的邊緣。
為了算法的穩定性，每3*3的矩形區域內（4*4或者5*5也沒說不可以的），路徑的起始角度周期性的改變，這樣相鄰像素的Kg條路徑經過的區域就有着較大的不同能得到更為有效的結果集。
由上圖可以看到，在不少情況下，未知點的前景和背景取樣數並不能達到Kg個，甚至極端情況下，找不到任何一個取樣點，這樣該點就無法進行透明度的計算了，這就要靠后面的過程了。
前景取樣點數量分布 背景取樣點數量分布 前景+背景取樣點數量分布
上圖繪制了前面列舉的TriMap圖中未知區域每個部位的取樣點數量分布情況，顏色越靠近白色，表明取樣點的數量越大，從圖中可以明顯看出，處於圖像角落的一些未知點取樣情況並不是特別理想，但基本上未出現沒有取到樣的情況，那我們在來看看scribbles那張圖的結果。
前景取樣點數量分布 背景取樣點數量分布 前景+背景取樣點數量分布
特別是前景取樣分布的結果似乎不太令人滿意，有些部分取樣數為0了，這個問題下面還會談到。
在完成取樣計算后，我們就需要找出這些取樣點中那些是最佳的組合，這個時候就涉及到一般優化時常談到的目標函數了，在這篇論文中，對目標函數用了四個小函數的乘積來計算，分別如下：
1： 
其中 
為了全面，我們將上式中αp的計算公式列出： 
公式（2）的道理很為明顯，用一對F/B算出的α值如果很合理的話，那么用α結合F/B重新計算出的顏色應該和原始顏色的差距很小。公式（3）在表明在一定的領域內，由於像素一般不會有突變，差值的平均值也應該很小。
為方便理解，我貼出計算α的部分代碼：
/// <summary> /// 通過當前點、前景點以及背景點的顏色值計算對應的Alpha值，對應論文的公式（12）。 /// </summary> /// <param name="BC、GC、RC">當前點的BGR顏色分量值。</param> /// <param name="BF、GF、RF">前景點的BGR顏色分量值。</param> /// <param name="BF、GF、RF">背景點的BGR顏色分量值。</param> /// <remarks>Alpha會出現不在[0,1]區間的情況，因此需要抑制。</remarks> double CalcAlpha(int BC, int GC, int RC, int BF, int GF, int RF, int BB, int GB, int RB) { double Alpha =(double) ((BC - BB) * (BF - BB) + (GC - GB) * (GF - GB) + (RC - RB) * (RF - RB)) / ((BF - BB) * (BF - BB) + (GF - GB) * (GF - GB) + (RF - RB) * (RF - RB) + 0.0000001); // 這里0.0000001換成Eps在LocalSmooth階段似乎就不對了，有反常的噪點產生 if (Alpha > 1) Alpha = 1; else if (Alpha < 0) Alpha = 0; return Alpha; }
2： 作者考慮在未知點到取樣的前景和背景點之間的直線路徑上，應該盡量要少有像素的突變，比如如果這條路徑需要經過圖像的邊緣區域，則應該設計一個函數使得該函數的返回值較大，於是作者使用了下面的公式：
上式即沿着路徑對像素顏色進行積分，離散化后也就是一些累加，CSDN的提供的代碼在這個函數的處理過程中是有錯誤的，因為他最后一個判斷條件使得循環只會進行一次，有興趣的朋友可以自己去改改。
按照公式（4）的意義，一個未知點屬於前景的可能性可由下式表示：
、
而一個好的組合也應該最小化下式：
3、未知點和前景點之間的物理距離，一個好的組合中的前景點應該要盡量靠近未知點；
4、未知點和背景點之間的物理距離，一個好的組合中的背景點也應該要盡量靠近未知點；
將這四個條件組合起來，最終得到如下的目標函數：
各子項的指數數據可詳見論文本身。
按照這個要求，對前面進行取樣得到數據進行處理，並記錄下使上式最小的那一對組合，就初步確定了最佳的取樣點。
其實，這個時候我們也就可以初步獲得處理后的α值了，比如對於我們前面所說的Trimap，其原始圖像及經過sample和gather處理后的結果如下圖：
從處理結果看，已經可以粗略的得到處理的效果了。
2.3、Refinement
初步的gather處理后，正如前文所說，得到的結果還不夠細膩，並且有些未知點由於采樣的過程未收集到有效的前景和背景數據，造成該點無法進行處理，因此，在Refinement階段需要進一步解決這個問題。
論文提出，首先，在一定的鄰域內，比如半徑為5的領域內，首先統計出公式（2）對應的MP值最小的3個點相關顏色數據，並對這些數據進行加權平均，得到數據對：
然后按照下面這些公式計算新的前景、背景、透明度及可信度的計算。
可信度的計算是為下一步的局部平滑做准備的，他反應了我們在這一步確定的取樣點是否合理程度的一個度量，經由此步驟，我們可得到的透明度和合成圖如下所示：
可見在這一步得到的結果對於上圖來說已經相當完美了。
2.4 Local Smoothing
這一步說實在的我沒有花太多的精力去看，他的實現過程大概有點類似於高斯模糊，但里面多了很多其他方面的處理，一個很好的事情就是在CSDN提供的代碼中對這部分每個公式的實現都是正確的，也是完整的，因此，有興趣的朋友需要自己多看下論文和對應的代碼了。
三、算法的效果
按照論文提供的相關資料集我自己搜集的一些圖及配套的Trimap測試了該算法的一些結果，現貼出如下所示：
原圖 Trimap 合成后的效果圖
可見，對於這些Trimap圖，在很多情況下是能獲得較為滿意的效果的。
我還找了一些簡單的圖，使用scribble的方式進行處理，效果如下所示：
原圖 操作界面 結果圖
我選的這些都是背景比較簡單的圖，因此還能獲得較為理想的效果，如果是比較復雜的圖，使用scribble是基本上獲取不到很理想的效果的，除非人工仔細的划分邊界。
四、編程實現
在編程實現方面，CSDN提供的那個代碼基本的意思已經達到了，並且里面的函數意義也非常清晰，只不過他使用的opencv的庫的，相信專心去研究摳圖的人，把他改成其他語言也不是個難題，比如里面用到的vector在C#中就可以用list代替，那些Opencv的結構體也可以在C#中重新定義。
不過那個代碼占用的內存非常厲害，這主要是由於VECTOR等數據類型決定的，實際上這里完全可以用數組來搞定的。
我貼一部分代碼大家看看：
/// <summary> /// 對每個未知區域的像素按設定的循環角度搜索有效的前景和背景取樣點 /// </summary> /// <param name="X">未知區域的X坐標。</param> /// <param name="Y">未知區域的Y坐標。</param> /// <param name="F">用於保存前景點的內存區域。</param> /// <param name="B">用於保存背景點的內存區域。</param> /// <param name="CountF">最終獲取的前景點的數量。</param> /// <param name="CountB">最終獲取的背景點的數量。</param> /// <remarks>對於有些點，是有可能獲取不到有效的點的。</remarks> void Sample(int X, int Y, Point *F, Point *B, int &CountF, int &CountB) { int Z, XX, YY, Alpha; bool F1, F2; double InitAngle, IncAngle, Angle; double Dx, Dy, Step, Value, XD, YD, StepX, StepY; IncAngle = 360 / KG; // 每次掃描增加的角度 InitAngle = (Y % 5 * 5 + X % 25) * IncAngle / 25; // 起始角度，范圍在[0,IncAngle]之間，按照3*3的方式輪流替換，有利於提供結果的穩定性，如果KG取值較大，也可以使用4*4或者5*5 CountF = 0; CountB = 0; // 起步時需要記為0，注意參數中的引用（&） for (Z = 0; Z < KG; Z++) { F1 = false; F2 = false; // 開始尋找，暫時未找到任何前景點和背景點 Angle = (InitAngle + Z * IncAngle) / 180.0f * 3.1415926f; // 每次搜索的角度 Dx = cos(Angle); Dy = sin(Angle); Step = min(1.0f / (abs(Dx) + 1e-10), 1.0f / (abs(Dy) + 1e-10)); XD = X + 0.5; YD = Y + 0.5; // +0.5使用於int四舍五入取整，相當於其他語言的round函數 StepX = Step * Dx ; StepY = Step * Dy ; // StepX和StepY中必然有一個為1，另外一個小於1，這樣保證每次搜索的像素點都會不同，一個好的建議是加大這個循環步長 // 這樣有兩個好處。1：加快查找速度；2：讓搜索點可以進入已知點的內部，從而避免了只從已知點的邊緣取樣。但如果已知點的區域狹長，可能丟失有效取樣點。 for (; ;) { XX = int(XD); if (XX < 0 || XX >= Width) break; // 已經超出了圖像的邊界了，結束循環 YY = int(YD); if (YY < 0 || YY >= Height) break; XD += StepX; YD += StepY; Alpha = Mask[YY * MaskStride + XX]; // 得到路徑中該點的特征值（前景/背景/未知) if (F1 == false && Alpha == 0) // 如果沿這條路徑尚未找到背景點，並且改點具有背景點的特征，則記錄下改點到背景點序列中 { B[CountB].X = XX; // 背景點的X坐標 B[CountB].Y = YY; // 背景點的Y坐標 CountB++; // 背景點數量增加1 F1 = true; // 在此路徑已經找到了背景點，不用再找背景點了。 } else if (F2 == false && Alpha == 255) // 如果沿這條路徑尚未找到前景點，並且改點具有前景點的特征，則記錄下改點到前景點序列中 { F[CountF].X = XX; // 前景點的X坐標 F[CountF].Y = YY; // 前景點的X坐標 CountF++; // 前景點數量增加1 F2 = true; // 在此路徑已經找到了前景點，不用再找前景點了。 } else if (F1 == true && F2 == true) // 如果前景點和背景點都已經找到了，則結束循環 { break; } } } }
通過以上的Sample代碼，我們就可以避免使用Vector之類的數據結構了，速度和內存占用都會得到改進。
然后在Gather階段的代碼改成如下的方式。
void Gathering() { int X, Y, K, L, Index, Speed; int CountF, CountB, Fx, Fy, Bx, By; double Pfp, Min, Dpf, Gp; bool Flag; Point *F = (Point *) malloc(KG * sizeof(Point)); // 采用這種方式占用的內存小很多 Point *B = (Point *) malloc(KG * sizeof(Point)); for (Y = 0; Y < Height; Y++) { Index = Y * MaskStride; for (X = 0; X < Width; X++) { tuple[Index].Flag = -1; // 先都設置為無效的點 if (Mask[Index] != 0 && Mask[Index] != 255) // 只處理未知點 { Sample(X, Y, F, B, CountF, CountB); // 對當前點進行前景和背景取樣 Pfp = CalcPFP(X, Y, F, B, CountF, CountB); // 計算公式（5），因為公式（5）只於前景和背景取樣點的整體有關，無需放到下面的循環內部 Min = 1e100; Flag = false; for (K = 0; K < CountF; K++) // 對於每一個前景點 { Dpf = CalcDp(X, Y, F[K].X, F[K].Y, true); // 計算前景點到中心點的歐式距離 for (L = 0; L < CountB; L++) // 對於每一個背景點 { Gp = CalcGp(X, Y, F[K].X, F[K].Y, B[L].X, B[L].Y, Pfp, Dpf); // 按照公式（7）計算目標函數 if (Gp < Min) { Min = Gp; Fx = F[K].X; Fy = F[K].Y; // 記錄下目標函數為最小值處的前景和背景點的坐標 Bx = B[L].X; By = B[L].Y; Flag = true; } } } if (Flag == true) // 說明找到了最好的組合了，如果找不到,則原因可能是：（1）Sample過程為找到任何有效的前景和背景點；（2）Sample過程只找到前景點或只找到背景點 { // 某個點找不到也不用怕，可能在下面的Refine過程中（其算法為領域處理）得以彌補。 Speed = Fy * Stride + Fx * 3; // 記錄下最佳前景點的顏色值 tuple[Index].BF = ImageData[Speed]; tuple[Index].GF = ImageData[Speed + 1]; tuple[Index].RF = ImageData[Speed + 2]; Speed = By * Stride + Bx * 3; tuple[Index].BB = ImageData[Speed]; // 記錄下最佳背景點的顏色值 tuple[Index].GB = ImageData[Speed + 1]; tuple[Index].RB = ImageData[Speed + 2]; tuple[Index].SigmaF = CaclSigma(Fx, Fy); // 計算前景點周邊像素的均方差值 tuple[Index].SigmaB = CaclSigma(Bx, By); // 計算背景點周邊像素的均方差值 tuple[Index].Flag = 1; // 這個像素是個有效的處理點 } } Index++; } } free(F); free(B); }
由於CSDN提供了代碼，其他的代碼我這里就不提供了。
五、算法的缺陷
以下內容純屬個人意見，請各位閱讀的朋友根據自己的知識去判斷。
1、Sample過程存在潛在的問題：論文的圖2闡述了對某點進行取樣的操作過程，這個過程在第一次遇到前景或背景點時就把該點視為前景或背景的一個取樣點。這樣的話所有的取樣點的數據都只可能取在Trimap或者scripple圖的前景和背景區域的邊緣處。如果邊緣處的像素和內部的有很大的區別，顯然會出現部分取樣點取到的數據很不合理，從而導致最終的透明度信息錯誤。一種可行的方案就是在Sample的過程中修改每次取樣的步長，比如，下面的方式：
StepX = Step * Dx * 10 ; StepY = Step * Dy * 10 ;
即每次沿X或者Y方向移動10個像素，這樣就會使得取樣點會進入已知區域內部最多10個像素左右，同時也會加速取樣的速度。
不過這樣做也隱形的帶來一個新的問題，即有可能會丟失一些取樣點，當已知點的寬度小於10像素時出現這種情況。
2、還是會存在某些點無法獲取有效的取樣點的，即使有了refine過程。
3、算法的復雜度還是很高，雖然整體的並行性很好，如果沒有GPU參與，純CPU實現起來一樣很吃力。
六、效果驗證
我用C寫了個DLL，並提供了三個示例程序，給有興趣的朋友測試下效果： http://files.cnblogs.com/Imageshop/SharedMatting.rar
*********************************作者： laviewpbt 時間： 2014.2.15 聯系QQ: 33184777 轉載請保留本行信息************************