學習了一篇shared sample技術實作圖像摳圖效果的文章，并就文章中一些細節進行了解釋，提供了大家一些關于該技術配套的資源和代碼，并給出該算法的部分測試效果。
一、序言
陸陸續續的如果累計起來，我估計至少有二十來位左右的朋友加我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：
公式（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，并提供了三個示例程式，給有興趣的朋友測試下效果： https://files.cnblogs.com/Imageshop/SharedMatting.rar
*********************************作者： laviewpbt 時間： 2014.2.15 聯系QQ: 33184777 轉載請保留本行資訊************************