現在這個社會發展的太快，到處都充斥著各種各樣的資源，各種開源的平臺，如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