Euclidean distance map（EDM)這個概念可能聽過的人也很少，其主要是用在二值影象中，作為一個很有效的中間處理手段存在。一般的處理都是將灰度圖處理成二值圖或者一個二值圖處理成另外一個二值圖，而EDM演算法確是由一幅二值圖生成一幅灰度圖。其核心定義如下：

　　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.

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

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

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

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

　　直接從影象的前景點搜尋和其最接近的背景點的距離是一個極其低效和耗時的操作。一些研究者通過只取幾個方向距離也實現了一些不同的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.

　　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.

　　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