本文轉自：http://www.cnblogs.com/Imageshop/p/6193433.html
現在這個社會發展的太快，到處都充斥着各種各樣的資源，各種開源的平台，如github,codeproject,pudn等等，加上一些大型的官方的開源軟體，基本上能找到各個類型的代碼。很多初創業的老闆可能都曾經說過基本上我的程式員不需要自己寫算法，但是他們要學會搜尋，強有力的搜尋能力基本能解決可能會遇到的一切問題，比如前一段時間流行的prisma濾鏡，現在你在github上能找到一大堆了。
可我想表達的并不是這個，上述這個情況确實是正确的。但是，我相信這個世界上90%的開源代碼是垃圾代碼，還有9%的代碼是理想代碼，能夠有1%的是優質代碼就已經是相當不錯了。基本上我們拿到的網絡中的某個參考代碼都還是要經過自己的細心改良和優化才能真正的應用于項目中，而這部分能力并不是每個人都有。
超像素經典的算法SLIC就屬于上述1%的一員，他有論文的介紹原理性的東西，有數學公式的推導，有和其他算法的比較資料，更重要的是他還有和論文完全對應的參考代碼，而且有C++、matlab以及GPU版本，可以說是非常有良心的一篇論文。雖然是優質代碼，但是當你真正的去研究他的代碼時，你就會發現離實際的應用距離還有很遠的路要走：可怕的記憶體占用，大量的浮點計算還是很客觀的時間開銷。你在網絡上搜尋，包括github上，可以找到一些相關的用SLIC做圖像分割的代碼，而百度上所搜SLIC也能找到一些部落格園或者CSDN的部落格的介紹，不過大部分都停留在對源代碼的解釋上。
我這裡不對論文的思想做詳細的解釋，相關的論文和源代碼可以參考：http://ivrl.epfl.ch/supplementary_material/RK_SLICSuperpixels/index.html
為了描述友善，這裡貼下原文對算法流程的描述：
我們來一條一條的優化和整理。
首先，論文提到圖像資訊是取自于CIELAB空間而不是RGB空間，在論文中有相關的描述如下：CIELAB color space is widely considered as perceptually uniform for small color space。我們實際程式設計測試時也發現在相同的參數下，很明顯CIELAB空間的分割效果明顯要比RGB空間的規則很多（如下兩圖所示）。這個可以從LAB空間的球形模型（左圖）得到闡釋吧。但是也不是說RGB空間的結果就完全不行，我們如果适當的調大M參數的值，也能獲得不錯的結果。唯一不同的時，使用LAB空間在時間和空間上會稍微多一點。
接着，我們來看看作者的代碼，首先是RGB2XYZ、RGB2LAB、DoRGBtoLABConversion三個函數，很明顯，這是RGB空間轉換為LAB空間的算法，XYZ空間隻是作為兩者轉換的中間件存在的，作者代碼裡采用了double類型資料來儲存轉換的結果，這是嚴格意義上的轉換，作為論文的算法驗證來說，是無可厚非的，但是拿到工程中使用可以說是一場災難。8倍于原始圖像資料的額外的記憶體占用在很多場景下就已經無法使用了，還有可怕的浮點計算的硬體屏障（不要以為大家都是在PC上使用圖像^^），我們必須想辦法降低這個記憶體占用和複雜度。
在我的博文 顔色空間系列2: RGB和CIELAB顔色空間的轉換及優化算法 中，提出了一種快速算法，可以無任何浮點計算快速的将RGB轉換到和原圖占用記憶體一樣大小的記憶體空間中，而後續的編碼也證明這種轉換的精度損失對于結果的影響是完全在可以接受的範圍内的。
LAB空間結果 RGB空間結果
下一步就是初始化聚類中心，對應的函數在GetKValues_LABXYZ或者GetLABXYSeeds_ForGivenStepSize中，按照使用者提供的超像素的大小或者超像素的個數來均勻的分布取樣的XY位置和LAB的值，作者用vector來儲存這些值，在源代碼的其他很多地方，作者也使用了vector，就我看來，如果不是迫不得已，我基本不使用這個東西。特别作為本例，并沒有使用到vector的啥進階特性，如果直接自己動态的配置設定記憶體，在很多方面都有着特别的靈活性，而且封裝的東西效率很多情況下并不如意。
同樣的道理，從速度方面考慮，再一次，我們也不使用源代碼中的double類型來儲存中心點的坐标，而是使用整形，是否可行，一切皆以最後的結果說話，實際就證明此法可行。在GetKValues_LABXYZ最後，作者為了避免初始中心點落入圖像的邊緣或者噪音點，進行了一次重定位操作，即取初始中心點周邊3*3領域梯度值最小的那個像素為新的中心點，原文的話為：The centers are moved to seed locations corresponding to the lowest gradient position in a 3 × 3 neighborhood. This is done to avoid centering a superpixel on an edge, and to reduce the chance of seeding a superpixel with a noisy pixel 。坦白的說，個人這個其實沒有啥必要，因為一般情況下圖像的邊緣寬度也非一個像素寬，而噪音一般也非一個孤立點，但是無論如何，做一下也好。
于是就涉及到了圖像梯度值的計算了，其實說白了也就是圖像的邊緣算法，對應于源代碼的DetectLabEdges函數，注意這裡是在LAB空間的邊緣檢測了。哇，又是一堆的浮點計算。實際上，計算這個邊緣資訊隻用L通道的資訊就完全足夠了，L就是圖像的明度嗎，也就表面了一個像素的強度資訊。鑒于前面我們使用的整形的資料，這裡我給出一個簡單的單通道的邊緣計算方式，還支援In-Place操作，以及對圖像四周的一圈像素也是進行了計算的。
// 基于Sobel算子的邊緣檢測算法
// 支援In-Place操作，即Dest可以等于Src,當然處理後Src的資料就被改變了
// 四周邊緣像素同樣得到了處理
void GetEdgeGradient(unsigned char *Src, unsigned char *Dest, int Width, int Height)
{
const int Channel = 1;
unsigned char *RowCopy = (unsigned char *)malloc((Width + 2) * 3 * Channel);
unsigned char *SqrValue = (unsigned char *)malloc(65026 * sizeof(unsigned char));
unsigned char *First = RowCopy;
unsigned char *Second = RowCopy + (Width + 2) * Channel;
unsigned char *Third = RowCopy + (Width + 2) * 2 * Channel;
for (int Y = 0; Y < 65026; Y++) SqrValue[Y] = (int)(sqrtf(Y * 1.0) + 0.49999f);
memcpy(Second, Src, Channel);
memcpy(Second + Channel, Src, Width * Channel); // 拷貝資料到中間位置
memcpy(Second + (Width + 1) * Channel, Src + (Width - 1) * Channel, Channel);
memcpy(First, Second, (Width + 2) * Channel); // 第一行和第二行一樣
memcpy(Third, Src + Width * Channel, Channel); // 拷貝第二行資料
memcpy(Third + Channel, Src + Width * Channel, Width * Channel);
memcpy(Third + (Width + 1) * Channel, Src + Width * Channel + (Width - 1) * Channel, Channel);
for (int Y = 0; Y < Height; Y++)
{
unsigned char *LinePD = Dest + Y * Width;
if (Y != 0)
{
unsigned char *Temp = First; First = Second; Second = Third; Third = Temp;
}
if (Y == Height - 1)
{
memcpy(Third, Second, (Width + 2) * Channel);
}
else
{
memcpy(Third, Src + (Y + 1) * Width * Channel, Channel);
memcpy(Third + Channel, Src + (Y + 1) * Width * Channel, Width * Channel); // 由于備份了前面一行的資料，這裡即使Src和Dest相同也是沒有問題的
memcpy(Third + (Width + 1) * Channel, Src + (Y + 1) * Width * Channel + (Width - 1) * Channel, Channel);
}
for (int X = 0; X < Width; X++)
{
int GradientH, GradientV;
if (X == 0)
{
GradientH = First[X + 0] + First[X + 1] + First[X + 2] - (Third[X + 0] + Third[X + 1] + Third[X + 2]);
}
else
{
GradientH = GradientH - First[X - 1] + First[X + 2] + Third[X - 1] - Third[X + 2];
}
GradientV = First[X + 0] + Second[X + 0] * 2 + Third[X + 0] - (First[X + 2] + Second[X + 2] * 2 + Third[X + 2]);
int Value = (GradientH * GradientH + GradientV * GradientV) >> 1;
if (Value > 65025)
LinePD[X] = 255;
else
LinePD[X] = SqrValue[Value];
}
}
free(RowCopy);
free(SqrValue);
} 至于上述的最後的量化操作，也可以考慮使用abs(GradientH) + abs(GradientV)來實作。
下一步就是最核心的計算了，通過聚類的方式疊代計算新的聚類中心，具體過程見源代碼中的PerformSuperpixelSLIC函數，我們這裡提下優化問題。
疊代計算聚類中心是本算法的核心，而疊代的核心就是計算距離，因為SLIC的聚類的資料源是LAB和XY的綜合體，而LAB和XY各自的取值範圍不同，如果直接采用歐式距離計算，則會産生非常混亂的效果，是以，作者提出了如下的計算公式：
通過兩個參數m和S來協調兩種距離的比例配置設定。參數S表示XY空間的最大的可能值，這個通過使用者輸入的超像素大小可以自動計算出來（詳見論文），而參數M為LAB空間的距離可能最大值，論文提出其可取的範圍建議為[1,40]。
讓我們來仔細看看上述公式，似乎運算量比較大， 有三次開平方，有除法還有多次乘法。 首先，由于隻是距離比較，而不是基于距離資料的某種名額計量，是以，D'公式中的開平方計算是無需的。其實，在D’公式的開平方内部，Dc和Ds都有平方計算，是以，Dc和Ds計算式中的開平方計算是無需的。這已經大大的減少了計算量，接着，我們把D'開平方内部資料乘以m*m*s*s可以得到下式：
是以，我們隻需要比較這個式子的值就是比較距離的相對大小了，前面已經說了LAB和XY我們都采用整形資料，聚類的中心也是整形，是以，如果S和M也是整形值，則上述計算的所有計算都是整形，而S值當然可以是整形，M值得範圍[1,40]，當然也可以取整形了, So, 這裡我們優化後就隻有整形計算了，完全避免了浮點的計算。
但是要仔細的分析下上述式子的取值範圍，如果超出了int32所能表達的範圍，那在32位系統上就有所麻煩了，衆所周知，在32位系統上int64資料類型的計算速度是不快的。當然如果是64位系統這個問題是不會存在的。我們經過上述的LAB顔色空間轉換後，LAB各分量的範圍都在[0,255]之間，是以,
這個值最大不會超過255 * 255 * 3，而S*S等于超像素的大小，超像素一般不會大于100*100把，再大也沒啥意義了。同理
表示了超像素中心點和其他的距離，也一般不會大于100*100，而m為使用者輸入值，前面說了m最大取值40，綜合所屬上面的式子的最大值是不會超過int32所能表達的最大範圍的，具體的代碼如下：
//// ******* 疊代更新中心點 ******** //
int PowerStep = Step * Step;
int PowerM = M * M;
int *SumL = (int *)_aligned_malloc(ActualSPAmount * sizeof(int), 16);
int *SumA = (int *)_aligned_malloc(ActualSPAmount * sizeof(int), 16);
int *SumB = (int *)_aligned_malloc(ActualSPAmount * sizeof(int), 16);
int *SumX = (int *)_aligned_malloc(ActualSPAmount * sizeof(int), 16);
int *SumY = (int *)_aligned_malloc(ActualSPAmount * sizeof(int), 16);
int *SumC = (int *)_aligned_malloc(ActualSPAmount * sizeof(int), 16);
//int *SumMem = (int *)_aligned_malloc(ActualSPAmount * 6 * sizeof(int), 16);
//int *SumL = SumMem + ActualSPAmount * 0, *SumA = SumMem + ActualSPAmount * 1, *SumB = SumMem + ActualSPAmount * 2;
//int *SumX = SumMem + ActualSPAmount * 3, *SumY = SumMem + ActualSPAmount * 4, *SumC = SumMem + ActualSPAmount * 5;
unsigned short *TempLabel = (unsigned short *)_aligned_malloc(Width * Height * sizeof(unsigned short), 16);
unsigned int *Distance = (unsigned int *)_aligned_malloc(Width * Height * sizeof(unsigned int), 16);
for (int Iteration = 0; Iteration < Iter; Iteration++)
{
memset(Distance, 255, Width * Height * sizeof(int)); // 把距離設定為最大值
for (int Index = 0; Index < ActualSPAmount; Index++) // 周遊所有的超像素
{
int CenterL = SeedL[Index];
int CenterA = SeedA[Index];
int CenterB = SeedB[Index];
int CenterX = SeedX[Index];
int CenterY = SeedY[Index];
int X1 = CenterX - Step; // 以目前超像素的中心點位中心，向四周擴散Step距離
int Y1 = CenterY - Step;
int X2 = CenterX + Step;
int Y2 = CenterY + Step;
if (X1 < 0) X1 = 0; // 注意不要越界
if (X2 > Width) X2 = Width;
if (Y1 < 0) Y1 = 0;
if (Y2 > Height) Y2 = Height;
for (int Y = Y1; Y < Y2; Y++)
{
int Pos = Y * Width + X1;
int PowerY = (Y - CenterY) * (Y - CenterY);
for (int X = X1; X < X2; X++)
{
int CurrentL = L[Pos];
int CurrentA = A[Pos];
int CurrentB = B[Pos];
int DistLAB = (CurrentL - CenterL) * (CurrentL - CenterL) + (CurrentA - CenterA) * (CurrentA - CenterA) + (CurrentB - CenterB) * (CurrentB - CenterB);
int DistXY = (X - CenterX) * (X - CenterX) + PowerY;
int Dist = PowerStep * DistLAB + DistXY * PowerM; // 整形化
if (Dist < Distance[Pos])
{
Distance[Pos] = Dist; // 記錄最小的距離
TempLabel[Pos] = Index;
}
Pos++;
}
}
}
memset(SumL, 0, ActualSPAmount * sizeof(int)); // 把所有的累加資料清零
memset(SumA, 0, ActualSPAmount * sizeof(int));
memset(SumB, 0, ActualSPAmount * sizeof(int));
memset(SumX, 0, ActualSPAmount * sizeof(int));
memset(SumY, 0, ActualSPAmount * sizeof(int));
memset(SumC, 0, ActualSPAmount * sizeof(int));
//memset(SumMem, 0, ActualSPAmount * 6 * sizeof(int)); // 把所有的累加資料清零
for (int Y = 0; Y < Height; Y++)
{
int Pos = Y * Width;
for (int X = 0; X < Width; X++)
{
int Index = TempLabel[Pos];
SumL[Index] += L[Pos];
SumA[Index] += A[Pos];
SumB[Index] += B[Pos]; // 累加
SumX[Index] += X;
SumY[Index] += Y;
SumC[Index]++;
Pos++;
}
}
for (int Index = 0; Index < ActualSPAmount; Index++)
{
int Amount = SumC[Index];
if (Amount == 0) Amount = 1;
SeedL[Index] = SumL[Index] / Amount;
SeedA[Index] = SumA[Index] / Amount; // 重新計算新的中心點
SeedB[Index] = SumB[Index] / Amount;
SeedX[Index] = SumX[Index] / Amount;
SeedY[Index] = SumY[Index] / Amount;
}
}
_aligned_free(L); // 減少綜合記憶體占用量
_aligned_free(A);
_aligned_free(B);
_aligned_free(SumL);
_aligned_free(SumA);
_aligned_free(SumB);
_aligned_free(SumX);
_aligned_free(SumY);
_aligned_free(SumC);
//_aligned_free(TempLabel);
_aligned_free(Distance);
//_aligned_free(SumMem); 幾個細節描述下：
（1）上述代碼有幾句被我注釋掉了，我的本意是一次性配置設定足夠的記憶體，然後在分給其他的變量，這樣有些代碼寫起來簡潔些比如清零和釋放，但是我測試時發現對我預設的那個測試圖，被注釋掉的代碼會慢的比較明顯（170ms和140ms），但是換一副圖差別就不太明顯了，後面我分析了下對于我那些測試圖，其變量ActualSPAmount恰好等于1024，這個看似和對齊有很大的關系，不知道各位童鞋是怎麼認為的。
（2）上述代碼的中Distance資料我使用了unsigned類型的，這主要是為了後面能友善使用：
memset(Distance, 255, Width * Height * sizeof(int)); // 把距離設定為最大值
這一句代碼，因為如果使用的是signed類型的，則隻能用一個for循環給每個變量指派為int類型的最大值了，memset肯定要比for快一些。
在疊代完成後，大部分工作就已經基本完成了，論文中提出大概10次疊代就足夠了，其實我們實測，經過一次疊代就基本已經比較靠譜了，上述流程中那個計算residual error E的過程是不需要的，耗時有耗力，實際應用中，我覺得取疊代次數為4基本能平衡速度和精度了。
在最後，論文提出了一些後處理的過程，這主要是為了去除前面分割過程中的産生的一些比較小的分割塊，相關的代碼在EnforceLabelConnectivity中，這個算法的核心思想就是利用區域生長法找到圖像中每個超像素的大小，然後把過于小的超像素合并到周邊大的超像素中，這裡有幾個問題其實值得商榷：
（1）過小的超像素合并到周邊哪個超像素中呢，論文是采用的是找到的最後的相鄰的超像素，其實抑或是采用找到的第一個相信的超像素也好，都存在一個問題，這合理嗎？我們看下面3個圖：
局部原圖 不進行EnforceLabelConnectivity前的結果 進行EnforceLabelConnectivity後
大家注意上述的第三張圖，進行了EnforceLabelConnectivity後，紅色圓圈内的結果，背景和那個項鍊已經練成了同一個區域，顯然這個是不合理的，也是不正确的，對于後續的分割或者其他操作都是不利的，而未進行前，這兩個個區域明顯是分開的，是以這個操作到底要如何進行呢？？？？？？？？？？？？？？
（2）很多人可能都沒有注意到進行了EnforceLabelConnectivity後，超像素的個數可能不會有任何減少，甚至還會出現增加的情況，這豈不逆天了。
EnforceLabelConnectivity不就是合并一些小的超像素到周邊去嗎，怎麼可能增加呢，我為了這個問題也是困惑了很久很久，甚至有好幾天走路吃飯都在想這個問題，有一次差點被車撞倒。也我曾懷疑作者的這個函數寫的有問題，也用自己以前一個非常靠譜的區域生長來代替他的代碼結果還是一樣。後來有一天我突然恍然大悟，大家都沒有錯，問題是出在超像素分割個本身過程的，由于聚類過程的特性，并不能保證每一類在XY空間裡都能連續的，比如上面中間那個圖裡，大量白色邊界重疊在一起的區域裡就有很多超像素在空間上已經分離了（其實這個時候就不能叫他為1個超像素了，而是2個或者多個），而在EnforceLabelConnectivity函數裡，是基于區域生長的，也就是基于空間連續性做的判斷，這樣就出現了上述問題。
最後，這代碼值得參考的還有各個标簽之間分界線的繪制的技巧，對應的函數在DrawContoursAroundSegments裡。
通過上述優化，加上其他的一些程式設計技巧（比如及時釋放不在需要的記憶體），在整個函數的執行過程，大概需要3.5倍額外的圖像記憶體，就可以完成整個過程，而作者提供的代碼，少說也要20倍吧，呵呵，速度方面的，1000*1000的圖，産生1600個超像素，疊代10次，大概需要300ms，而源代碼大概需要3S, 10倍的提速，實際上，如果疊代4次，150ms就可以獲得結果。
論文中提出，當m越大時，空間相似性越重要，超像素的結果就越緊湊，而m小時，則超像素則越靠近圖像的邊界，但是形狀越不規則，這個我們也可以從下面的實驗圖檔中看到：
局部原圖 m = 5 m = 2 m = 40
論文中還提出了自動參數的實作，即m值無需使用者輸入，在疊代過程自動更新，我也實作了下，覺得這個意義不大，而且或拖慢計算和速度和增加記憶體消耗。 通常我們取m固定為32左右就可以了。
我做了UI界面，如下圖，優化思路已經描述的相當清晰了，有能力者自可為之。
DEMO下載下傳位址：http://files.cnblogs.com/files/Imageshop/Slic.rar
超像素有很多應用，其中一個方向就是摳圖和圖像分割，下一階段我将抽空研究下這方面的東西。