SSE圖像算法優化系列二十五:二值圖像的Euclidean distance map（EDM)特徵圖計算及其優化。

時間 2019-11-06

標籤 sse 圖像算法優化系列二十五 euclidean distance map edm 特徵計算及其简体版

原文原文鏈接

　　Euclidean distance map（EDM)這個概念可能聽過的人也不多，其主要是用在二值圖像中，做爲一個頗有效的中間處理手段存在。通常的處理都是將灰度圖處理成二值圖或者一個二值圖處理成另一個二值圖，而EDM算法確是由一幅二值圖生成一幅灰度圖。其核心定義以下：html

　　The definition is simple enough: each point in the foreground is assigned a brightness value equal to its straight line (hence 「Euclidean」) distance from the nearest point in the background. java

　　或者用一句更簡潔的話說就是：The value of each pixel is the distance to the nearest background pixel (for background pixels, the EDM is 0)。算法

　　EDM由不少用處，好比因爲像素都是分佈在矩形格上，形成一些形態學上的操做存在一些方向偏移，以及在距離計算上的一些誤差均可以使用EDM技術克服，使用EDM來實現的腐蝕、膨脹或者開和閉操做相比普通的逐像素（模板操做）的結果來的更加各向同性化，以下圖所示（很明顯，基於EDM處理的圖結果更加光滑，沒有任何的方向性，而傳統的逐像素的算法則有點棱角分明）：ide

　　　　　　　　　　　　圖一：使用普通的逐像素的處理（從作至右依次是：二值原圖、閉操做、開操做，找到的邊界）函數

　　　　　　　　　　　　圖二：使用基於EDM的處理（從作至右依次是：二值原圖、閉操做、開操做，找到的邊界）post

　　直接從圖像的前景點搜索和其最接近的背景點的距離是一個極其低效和耗時的操做。一些研究者經過只取幾個方向距離也實現了一些不一樣的EDM效果，好比只有45或者90度的角度的距離，可是這些距離的效果和傳統的一些處理方法一樣存在這沒有足夠的各向同性的性質，所以，也不具備很大的意義，以下圖所示：測試

　　爲了快速的實現EDM值的計算，很多做者提出了一些快速算法，其中比較經典的是Danielsson在1980提出的算法，其經過兩次遍歷圖像給出了和原始計算同樣的效果，具體步驟以下所示：優化

　　1. Assign the brightness value of 0 to each pixel in the background and a large positive value (greater than the maximum feature width) to each pixel in a feature.url

　　2. Proceeding from left to right and top to bottom, assign each pixel within a feature a brightness value one greater than the smallest value of any of its neighbors.idea

　　3. Repeat step 2, proceeding from right to left and bottom to top.

　　在相關的資料中尋找參考的代碼，能夠找到以下一些資料：

The EDM algorithm is similar to the 8SSEDT in F. Leymarie, M. D. Levine, in: CVGIP Image Understanding, vol. 55 (1992), pp 84-94

https://linkinghub.elsevier.com/retrieve/pii/104996609290008Q

　　最直接的就是ImageJ了，在其開源的目錄：ImageJ\source\ij\plugin\filter下有Edm.java文件分享EDM算法的實現，不過那個代碼寫得比較噁心，我後面也不知道在什麼地方找到了一個比較清潔的代碼（也許是我本身改清潔的），核心思想是和這裏的一致的：

        public static int ONE = 100;            // 這個值和一個像素的距離對應 
        public static int SQRT2 = 141;          // 這個值和Sqr(2)個像素的距離對應
        public static int SQRT5 = 224;          // 這個值和Sqr(5)個像素的距離對應

        public static void CalaEuclideanDistanceMap(FastBitmap Bmp) { int X, Y, Index, Value,Max; int Width, Height, Xmax, Ymax ; byte* Pointer; Width = Bmp.Width; Height = Bmp.Height; Xmax = Width - 2; Ymax = Height - 2; // 1. Assign the brightness value of 0 to each pixel in the background and a large positive // value (greater than the maximum feature width) to each pixel in a feature.

            int[] ImageData = new int[Width * Height]; for (Y = 0, Index = 0; Y < Height; Y++) { Pointer = Bmp.Pointer + Y * Bmp.Stride; for (X = 0; X < Width; X++, Index++)  ImageData[Index] = Pointer[X]<<16; } // 2. Proceeding from left to right and top to bottom, assign each pixel within a feature a // brightness value one greater than the smallest value of any of its neighbors.

            for (Y = 0,Index=0; Y < Height; Y++) { Pointer = Bmp.Pointer + Y * Bmp.Stride; for (X = 0; X < Width; X++,Index++) { if (Pointer[X] > 0) { if ((X <= 1) || (X >= Xmax) || (Y <= 1) || (Y >= Ymax)) SetEdgeValue(Index, Width, ImageData, X, Y, Xmax, Ymax); else SetValue(Index, Width, ImageData); } } } // 3. Repeat step 2, proceeding from right to left and bottom to top.
            for (Y = Height - 1,Index=Width*Height-1; Y >= 0; Y--) { Pointer = Bmp.Pointer + Y * Bmp.Stride; for (X = Width - 1; X >= 0; X--,Index--) { if (Pointer[X] > 0) if ((X <= 1) || (X >= Xmax) || (Y <= 1) || (Y >= Ymax)) SetEdgeValue(Index, Width, ImageData, X, Y, Xmax, Ymax); else SetValue(Index, Width, ImageData); } } // Find the max value of the data 
            for (Y = 0,Max=0, Index = 0; Y < Height; Y++) for (X = 0; X < Width; X++, Index++) if (Max < ImageData[Index]) Max = ImageData[Index]; for (Y = 0, Index = 0; Y < Height; Y++) { Pointer = Bmp.Pointer + Y * Bmp.Stride; for (X = 0; X < Width; X++) { Value = (ImageData[Index]) * 255 / Max; Pointer[X] = (byte)Value; Index++; } } } private static void SetValue(int Offset, int Width, int[] ImageData) { int Value; int Index1 = Offset - Width - Width - 2; int Index2 = Index1 + Width; int Index3 = Index2 + Width; int Index4 = Index3 + Width; int Index5 = Index4 + Width; int Min = int.MaxValue; // * // * x * // *
            Value = ImageData[Index2 + 2] + ONE; if (Value < Min) Min = Value; Value = ImageData[Index3 + 1] + ONE; if (Value < Min) Min = Value; Value = ImageData[Index3 + 3] + ONE; if (Value < Min) Min = Value; Value = ImageData[Index4 + 2] + ONE; if (Value < Min) Min = Value; // * * // x // * *
            Value = ImageData[Index2 + 1] + SQRT2; if (Value < Min) Min = Value; Value = ImageData[Index2 + 3] + SQRT2; if (Value < Min) Min = Value; Value = ImageData[Index4 + 1] + SQRT2; if (Value < Min) Min = Value; Value = ImageData[Index4 + 3] + SQRT2; if (Value < Min) Min = Value; // * * // * * // x // * * // * * 
            Value = ImageData[Index1 + 1] + SQRT5; if (Value < Min) Min = Value; Value = ImageData[Index1 + 3] + SQRT5; if (Value < Min) Min = Value; Value = ImageData[Index2 + 4] + SQRT5; if (Value < Min) Min = Value; Value = ImageData[Index4 + 4] + SQRT5; if (Value < Min) Min = Value; Value = ImageData[Index5 + 3] + SQRT5; if (Value < Min) Min = Value; Value = ImageData[Index5 + 1] + SQRT5; if (Value < Min) Min = Value; Value = ImageData[Index4] + SQRT5; if (Value < Min) Min = Value; Value = ImageData[Index2] + SQRT5; if (Value < Min) Min = Value; ImageData[Offset] = Min; } private static void SetEdgeValue(int Offset, int Width, int[] ImageData, int X, int Y, int xmax, int ymax) { }

　　我刪除了邊緣部分的代碼（不但願懶人直接使用，直接拷貝過去就能用），代碼也很簡單，第一步是將原始圖像數據放大，放大的尺度文章中講的必定要大於最大的特徵的尺度，其實就是圖像中前景和背景最遠的那個距離值，實際上這個值不可能超過圖像的對角線長度，通常咱們隨便放大一個倍數就能夠了，上述代碼放大了65536倍，徹底足夠了，實際還能夠縮小的，對於背景色，原始值爲0，放大後仍是0。

　　那麼第二步，是從作到右，從上到下進行處理，處理時的原則先計算半徑爲1時周邊的4個點，爲了不浮點計算，咱們把距離的定義放大了100倍，好比距離爲1，就變爲100，距離爲根號2，就變爲141。這四個點若是他們的值有一個黑色，則確定就是這四個點中最小的值了，也是靠近改點的最小的距離了，若是四個點都是白色，咱們在把半徑擴展到根號2，查看這個時候四個點的顏色狀況，注意或者時候求最小值時須要加上141了，接着在看看根號5時的狀況了。

　　爲何只須要求這三個距離的最小值，我還一時沒有想清楚，可是這裏有一點就是，這個處理過程對每一個點是有前後依賴順序的，咱們看到在SetValue函數中當前點ImageData的更新會影響到下一個點的情況的，這也是有道理的，由於當前點的值處理完後就表示這個點距離最近的背景的距離，那麼下一個點計算時其和當前的距離在加上當前點和背景的距離就反應了下一個點和背景的距離。而這種先後依賴的關係對於咱們的後續的SSE形成了必定的影響。

　　咱們看下這段代碼的顯示結果：

　　　　　　　　　　原圖　　　　　　　　　　　　　　　　　　　　二值化　　　　　　　　　　　　　　　　　　　　EDM圖

　　EDM圖中越是白色的地方說明此處距離背景越遠，而途中的紅色圓圈則表示一個小小的孤立黑點也會對EDM圖產生不利的影響，好比，咱們若是把二值化後的圖進行簡單的去除白色和黑色孤點處理後，在進行EDM處理，則效果會平滑不少，以下所示：

　　　　　　　　二值化　　　　　　　　　　　　　　　　　　　　去除孤點　　　　　　　　　　　　　　　　　　EDM圖

　　此時的EDM圖中清晰的能夠看到5根明顯的分界線，這對於後續的一些分割等等識別工做頗有裨益。

　　爲測試速度，咱們選擇了一幅3000*2000的二值圖進行測試，由於這個算法的時間是和圖像內容有必定關係的（黑色的部分不須要處理），因此咱們的二值圖中有一半的顏色是白色，大概處理時間在160ms左右（上述C#的代碼)，並且未作太多的代碼優化，應該說時間仍是很快的。

　　時間沒有極限，我仍是但願能嘗試進一步加速，因而我仍是構思了實現了改算法的SSE優化。

　　先從簡單的搞起，第一，我前面說過放大的係數沒有必要到65536倍，只要這裏的值大於了對角線大小就OK了，好比咱們放大64倍，那白色的值就變爲255 * 64 = 16320，足夠日常使用了，這個時候有個好處就是咱們能夠不用int類型數據來記錄中間的距離了，而能夠直接使用unsigned short類型。此時最後一部分的量化咱們可飛快的用SSE實現：

void IM_GetMinMaxValue(unsigned short *Src, int Length, unsigned short &Min, unsigned short &Max)
{
    int BlockSize = 8, Block = Length / BlockSize;
    __m128i MinValue = _mm_set1_epi16(65535);
    __m128i MaxValue = _mm_set1_epi16(0);
    for (int Y = 0; Y < Block * BlockSize; Y += BlockSize)
    {
        __m128i SrcV = _mm_loadu_si128((__m128i *)(Src + Y));
        MinValue = _mm_min_epu16(MinValue, SrcV);
        MaxValue = _mm_max_epu16(MaxValue, SrcV);
    }
    Min = _mm_extract_epi16(_mm_minpos_epu16(MinValue), 0);
    Max = 65535 - _mm_extract_epi16(_mm_minpos_epu16(_mm_subs_epu16(_mm_set1_epi16(65535), MaxValue)), 0);        //    通過測試結果正確
    for (int Y = Block * BlockSize; Y < Length; Y++)
    {
        if (Min > Src[Y])    Min = Src[Y];
        if (Max < Src[Y])    Max = Src[Y];
    }
}
int IM_Normalize(unsigned short *Src, unsigned char*Dest, int Width, int Height, bool BmpFormat = false)
{
    if ((Src == NULL) || (Dest == NULL))                    return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                        return IM_STATUS_INVALIDPARAMETER;

    int Stride = BmpFormat == false ? Width : WIDTHBYTES(Width);
    unsigned short Min = 0, Max = 0;
    IM_GetMinMaxValue(Src, Width * Height, Min, Max);
    if (Max == Min)
    {
        memset(Dest, 128, Height * Stride);
    }
    else
    {
        float Inv = 255.0f / (Max - Min);          
        int BlockSize = 8, Block = Width / BlockSize;
        __m128i Zero = _mm_setzero_si128();
        for (int Y = 0; Y < Height; Y++)
        {
            unsigned short *LinePS = Src + Y * Width;
            unsigned char *LinePD = Dest + Y * Stride;
            for (int X = 0; X < Block * BlockSize; X += BlockSize)        //    正規化
            {
                __m128i SrcV = _mm_loadu_si128((__m128i *)(LinePS + X));
                __m128i Diff = _mm_subs_epu16(SrcV, _mm_set1_epi16(Min));
                __m128i Value1 = _mm_cvtps_epi32(_mm_mul_ps(_mm_cvtepi32_ps(_mm_cvtepu16_epi32(Diff)), _mm_set1_ps(Inv)));
                __m128i Value2 = _mm_cvtps_epi32(_mm_mul_ps(_mm_cvtepi32_ps(_mm_cvtepu16_epi32(_mm_srli_si128(Diff, 8))), _mm_set1_ps(Inv)));
                _mm_storel_epi64((__m128i *)(LinePD + X), _mm_packus_epi16(_mm_packus_epi32(Value1, Value2), Zero));
            }
            for (int X = Block * BlockSize; X < Width; X++)
            {
                LinePD[X] = (int)((LinePS[X] - Min) * Inv + 0.5f);
            }
        }
    }
    return IM_STATUS_OK;
}

　　由於在原始的代碼中，最小值必然爲0，因此原始的C#代碼沒有計算最小值，而本代碼爲了通用性，也計算了最小值，其過程也至關簡單，就是_mm_max_epu16和_mm_min_epu16的調用，而因爲數據是unsigned short的，也能夠快捷的使用_mm_minpos_epu16（其餘數據類型都沒有這個函數哦）得到8個ushort類型的最小值。後續的量化到unsigned char過程的也基本就是普通函數的調用，不存在其餘技巧，主要是要注意_mm_cvtepu16_epi32這種數據類型的轉換函數，要用的恰當。

　　注意到前面講的ImageData的更新是有先後依賴的，這就很是不利於我一次性計算多個像素的值，從而相似於SSE圖像算法優化系列九：靈活運用SIMD指令16倍提高Sobel邊緣檢測的速度（4000*3000的24位圖像時間由480ms下降到30ms）這篇文章的優化方式就沒法直接復現。可是咱們仔細觀察，發如今SetValue函數中，除了在半徑爲1的領域求最小值的計算中涉及到了本行值，其餘的在一行值的更新過程當中是不相互干涉的，那麼咱們能夠把這些單獨提取出來計算做爲中間值存儲，而後在用普通的C代碼和領域爲1的左右兩個位置的量進行比較。這樣同樣能起到不小的加速做用。

　　固然要注意一點，因爲這種依賴關係，在跳到下行計算時，更新過的這一行數據必須得以保留，不然又會獲得錯誤的結果。

　　咱們採用相似Sobel優化博文中的優化方式，採用了5行臨時數據緩衝區來解決邊緣處的處理問題，這樣就無需爲邊緣處在單獨寫代碼，同時，咱們沒有必要先對擴大後的數據進行計算，而時在5行數據更新的同時更新此部分數據，這樣在整個過程當中能夠少一部分拷貝工做。

//  1. Assign the brightness value of 0 to each pixel in the background and a large positive 
//  value (greater than the maximum feature width) to each pixel in a feature.
__forceinline void IM_GetExpandAndZoomInData(unsigned char *Src, unsigned short *Dest, int Width, int Height)
{
    const int BlockSize = 8, Block = Width / BlockSize;
    for (int X = 0; X < Block * BlockSize; X += BlockSize)
    {
        _mm_storeu_si128((__m128i *)(Dest + X + 2), _mm_slli_epi16(_mm_cvtepu8_epi16(_mm_loadl_epi64((__m128i *)(Src + X))), 6));
    }
    for (int X = Block * BlockSize; X < Width; X++)        Dest[X + 2] = Src[X] << 6;
    Dest[0] = Dest[2];       　　                 Dest[1] = Dest[Width > 1 ? 3 : 2];                        //    鏡像數據
    Dest[Width + 2] = Dest[Width + 1];        Dest[Width + 3] = Dest[Width > 1 ? Width : 2];
}

　　從左到右，從右到左方向的更新。

int IM_ShowEuclideanDistanceMap(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride)
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL))                        return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                            return IM_STATUS_INVALIDPARAMETER;
    if (Channel != 1)                                            return IM_STATUS_INVALIDPARAMETER;
    if (IM_IsBinaryImage(Src, Width, Height, Stride) == false)    return IM_STATUS_INVALIDPARAMETER;
    int Status = IM_STATUS_OK;
    int ExpandW = 2 + Width + 2;

    const int ONE = 100;            //  這個值和一個像素的距離對應     
    const int SQRT2 = 141;          //  這個值和Sqr(2)個像素的距離對應
    const int SQRT5 = 224;          //  這個值和Sqr(5)個像素的距離對應
    const int BlockSize = 8, Block = Width / BlockSize;
    unsigned short *Distance = (unsigned short *)malloc(Width * Height * sizeof(unsigned short));
    unsigned short *Buffer = (unsigned short *)malloc((ExpandW * 5 + Width) * sizeof(unsigned short));        //    5行緩衝用於記錄相鄰5行取樣數據，1行用於記錄中間結果
    if ((Distance == NULL) || (Buffer == NULL))
    {
        Status = IM_STATUS_OUTOFMEMORY;
        goto FreeMemory;
    }

    unsigned short *First = Buffer, *Second = First + ExpandW, *Third = Second + ExpandW;
    unsigned short *Fourth = Third + ExpandW, *Five = Fourth + ExpandW, *Temp = Five + ExpandW;
    
    //  2. Proceeding from left to right and top to bottom, assign each pixel within a feature a 
    //  brightness value one greater than the smallest value of any of its neighbors.

    IM_GetExpandAndZoomInData(Src, First, Width, Height);        //    把他們寫在這個過程當中，能夠減小一次賦值
    IM_GetExpandAndZoomInData(Src, Second, Width, Height);
    IM_GetExpandAndZoomInData(Src, Third, Width, Height);        //    前面三行數據相同
    IM_GetExpandAndZoomInData(Src + IM_ClampI(1, 0, Height - 1) * Stride, Fourth, Width, Height);
    IM_GetExpandAndZoomInData(Src + IM_ClampI(2, 0, Height - 1) * Stride, Five, Width, Height);

    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        if (Y != 0)
        {
            unsigned short *Temp = First; First = Second; Second = Third; Third = Fourth; Fourth = Five; Five = Temp;        //    交換指針
        }
        if (Y == Height - 2)            //    倒數第二行
        {
            memcpy(Five, Fourth, ExpandW * sizeof(unsigned short));
        }
        else if (Y == Height - 1)        //    最後一行，不能用第三行的數據，由於第三行是修改後的
        {
            IM_GetExpandAndZoomInData(Src + (Height - 1) * Stride, Fourth, Width, Height);
            IM_GetExpandAndZoomInData(Src + IM_ClampI(Height - 2, 0, Height - 1) * Stride, Five, Width, Height);
        }
        else
        {
            IM_GetExpandAndZoomInData(Src + (Y + 2) * Stride, Five, Width, Height);
        }
        for (int X = 0; X < Block * BlockSize; X += BlockSize)
        {
            __m128i SrcV = _mm_loadl_epi64((__m128i *)(LinePS + X));
            if (_mm_movemask_epi8(SrcV) != 0)            //    8個字節全是白色，則不須要繼續處理
            {
                __m128i FirstP1 = _mm_loadu_si128((__m128i *)(First + X + 1));
                __m128i FirstP3 = _mm_loadu_si128((__m128i *)(First + X + 3));

                __m128i SecondP0 = _mm_loadu_si128((__m128i *)(Second + X + 0));
                __m128i SecondP1 = _mm_loadu_si128((__m128i *)(Second + X + 1));
                __m128i SecondP2 = _mm_loadu_si128((__m128i *)(Second + X + 2));
                __m128i SecondP3 = _mm_loadu_si128((__m128i *)(Second + X + 3));
                __m128i SecondP4 = _mm_loadu_si128((__m128i *)(Second + X + 4));

                __m128i FourthP0 = _mm_loadu_si128((__m128i *)(Fourth + X + 0));
                __m128i FourthP1 = _mm_loadu_si128((__m128i *)(Fourth + X + 1));
                __m128i FourthP2 = _mm_loadu_si128((__m128i *)(Fourth + X + 2));
                __m128i FourthP3 = _mm_loadu_si128((__m128i *)(Fourth + X + 3));
                __m128i FourthP4 = _mm_loadu_si128((__m128i *)(Fourth + X + 4));

                __m128i FiveP1 = _mm_loadu_si128((__m128i *)(Five + X + 1));
                __m128i FiveP3 = _mm_loadu_si128((__m128i *)(Five + X + 3));
                //      *
                //  *   x   *
                //      *
                __m128i    Min0 = _mm_min_epu16(SecondP2, FourthP2);        //    由於第三行是先後相關的，因此不能在這裏參與計算

                //  *       *   
                //      x   
                //  *       *
                __m128i    Min1 = _mm_min_epu16(_mm_min_epu16(SecondP1, SecondP3), _mm_min_epu16(FourthP1, FourthP3));

                //      *       *   
                //  *               *   
                //          x    
                //  *               *   
                //      *       * 
                __m128i Min2 = _mm_min_epu16( _mm_min_epu16(_mm_min_epu16(FirstP1, FirstP3), _mm_min_epu16(SecondP0, SecondP4)),
                                                _mm_min_epu16(_mm_min_epu16(FourthP0, FourthP4), _mm_min_epu16(FiveP1, FiveP3)));
                
                __m128i Min = _mm_min_epu16(_mm_min_epu16(_mm_adds_epu16(Min0, _mm_set1_epi16(ONE)), _mm_adds_epu16(Min1, _mm_set1_epi16(SQRT2))), _mm_adds_epu16(Min2, _mm_set1_epi16(SQRT5)));

                _mm_storeu_si128((__m128i *)(Temp + X), Min);
            }
            else
            {
                memset(Temp + X, 0, 8 * sizeof(unsigned short));
            }
        }
        for (int X = Block * BlockSize; X < Width; X++)
        {
            if (LinePS[X] == 255)
            {
                unsigned short Min0 = IM_Min(Second[X + 2], Fourth[X + 2]);
                unsigned short Min1 = IM_Min(IM_Min(Second[X + 1], Second[X + 3]), IM_Min(Fourth[X + 1], Fourth[X + 3]));
                unsigned short Min2 = IM_Min(IM_Min(IM_Min(First[X + 1], First[X + 3]), IM_Min(Second[X + 0], Second[X + 4])),
                                             IM_Min(IM_Min(Fourth[X + 0], Fourth[X + 4]), IM_Min(Five[X + 1], Five[X + 3])));
                Temp[X] = IM_Min(IM_Min(Min0 + ONE, Min1 + SQRT2), Min2 + SQRT5);
            }
            else
            {
                Temp[X] = 0;
            }
        }
        for (int X = 0; X < Width; X++)
        {
            Third[X + 2] = IM_Min(IM_Min(Third[X + 1] + ONE, Third[X + 3] + ONE), Temp[X]);
        }
        Third[0] = Third[2];                        Third[1] = Third[Width > 1 ? 3 : 2];                        //    鏡像數據
        Third[Width + 2] = Third[Width + 1];        Third[Width + 3] = Third[Width > 1 ? Width : 2];

        memcpy(Distance + Y * Width, Third + 2, Width * sizeof(unsigned short));                //    複製回去
    }

    //  3. Repeat step 2, proceeding from right to left and bottom to top.
　　 //
    //　　此處代碼可自行添加
    //
    //
    Status = IM_Normalize(Distance, Dest, Width, Height, true);

FreeMemory:
    if (Distance != NULL)        free(Distance);
    if (Buffer != NULL)            free(Buffer);
    return Status;
}

　　看起來代碼量增長了很多，不過爲了效率都值得啊，通過測試，SSE優化後的速度從以前的160ms提高至60ms左右，加速仍是至關可觀的。

　上面代碼有一處有的BUG，留待有興趣研究的朋友自行查找。

那麼咱們來看看用EDM怎麼實現腐蝕或者膨脹這一類操做。在上述代碼中，最後的Distance數據中其實保存的就是某個點和最近的背景的距離，固然整個距離根據前面的代碼是放大了100倍的，並且注意到距離的最小值必然爲1（100），即在Distance數據中除了0以外的最小值爲100。若是咱們要進行腐蝕操做（即黑色部分向白色擴展），咱們能夠指定一個距離（即所謂的半徑），在Distance中若是數據小於這個距離，則變成了背景，而只有大於整個距離的部分纔會依舊是前景。注意到這裏所謂的距離能夠是小數（除以100後的結果），這就比傳統的半徑只能爲整數的需求進了一步。而膨脹操做，只須要反色後在進行距離變換，而後在和腐蝕同樣的操做便可。

　　咱們來作個測試，以下圖所示，很明顯的EDM版本的腐蝕是各向同性的，原來是個圓，處理後仍是個圓，而普通的算法，則逐漸的在想圓角矩形轉變。

　　　　　　　　　　原圖　　　　　　　　　　　　　　　　　　　半徑爲10的腐蝕（EDM版本）　　　　　　　　半徑爲10的腐蝕（普通版本）

　　固然，普通版本的腐蝕也能夠實現相似EDM那個效果的，好比咱們使用真正圓形的半徑的腐蝕（普通的半徑的意義是都是矩形半徑），可是一個核心的問題是隨着半徑的增長，圓形半徑的計算量會成倍上升，雖然圓形半徑也有快速算法，可是他始終沒法作到矩形矩形半徑那種O(1)算法的。而EDM算法確是和半徑無關的。

　　EDM還有不少其餘的用處，咱們能夠在ImageJ的代碼裏找到其更多的功能，這部分有待於讀者本身去研究。

　 Demo下載地址：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar