SSE圖像算法優化系列九：靈活運用SIMD指令16倍提升Sobel邊緣檢測的速度（4000*3000的24位圖像時間由480ms降低到30ms）。

這半年多時間，基本都在折騰一些基本的優化，有很多都是十幾年前的技術了，從随大流的角度來考慮，研究這些東西在很多人看來是浪費時間了，即不能賺錢，也對工作能力提升無啥幫助。可我覺得人類所謂的幸福，可以分為物質檔次的享受，還有更為複雜的精神上的富有，哪怕這種富有隻是存在于短暫的自我滿足中也是值得的。

閑話少說， SIMD指令集，這個古老的東西，從第一代開始算起，也快有近20年的曆史了，從最開始的MMX技術，到SSE，以及後來的SSE2、SSE3、SSE4、AVX以及11年以後的AVX2，逐漸的成熟和豐富，不過目前考慮通用性方面，AVX的輻射範圍還是有限，大部分在優化時還是考慮使用128位的SSE指令集，我在之前的一系列文章中，也有不少文章涉及到了這個方面的優化了。

今天我們來學習下Sobel算法的優化，首先，我們給出傳統的C++實作的算法代碼：

int IM_Sobel(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) && (Channel != 3))            return IM_STATUS_INVALIDPARAMETER;

    unsigned char *RowCopy = (unsigned char *)malloc((Width + 2) * 3 * Channel);
    if (RowCopy == NULL)    return IM_STATUS_OUTOFMEMORY;

    unsigned char *First = RowCopy;
    unsigned char *Second = RowCopy + (Width + 2) * Channel;
    unsigned char *Third = RowCopy + (Width + 2) * 2 * Channel;

    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 + Stride, Channel);                                                    //    拷貝第二行資料
    memcpy(Third + Channel, Src + Stride, Width * Channel);
    memcpy(Third + (Width + 1) * Channel, Src + Stride + (Width - 1) * Channel, Channel);

    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        unsigned char *LinePD = Dest + Y * Stride;
        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) * Stride, Channel);
            memcpy(Third + Channel, Src + (Y + 1) * Stride, Width * Channel);                            //    由于備份了前面一行的資料，這裡即使Src和Dest相同也是沒有問題的
            memcpy(Third + (Width + 1) * Channel, Src + (Y + 1) * Stride + (Width - 1) * Channel, Channel);
        }
        if (Channel == 1)
        {
            for (int X = 0; X < Width; X++)
            {
                int GX = First[X] - First[X + 2] + (Second[X] - Second[X + 2]) * 2 + Third[X] - Third[X + 2];
                int GY = First[X] + First[X + 2] + (First[X + 1] - Third[X + 1]) * 2 - Third[X] - Third[X + 2];
                LinePD[X] = IM_ClampToByte(sqrtf(GX * GX + GY * GY + 0.0F));
            }
        }
        else
        {
            for (int X = 0; X < Width * 3; X++)
            {
                int GX = First[X] - First[X + 6] + (Second[X] - Second[X + 6]) * 2 + Third[X] - Third[X + 6];
                int GY = First[X] + First[X + 6] + (First[X + 3] - Third[X + 3]) * 2 - Third[X] - Third[X + 6];
                LinePD[X] = IM_ClampToByte(sqrtf(GX * GX + GY * GY + 0.0F));
            }
        }
    }
    free(RowCopy);
    return IM_STATUS_OK;
}

　　代碼很短，但是這段代碼很經典，第一，這個代碼支援In-Place操作，也就是Src和Dest可以是同一塊記憶體，第二，這個代碼本質就支援邊緣。網絡上很多參考代碼都是隻進行中間有效的區域。具體的實作細節我不願意多述，自己看。

　　那麼Sobel的核心就是計算X方向的梯度GX和Y方向的梯度GY，最後有一個耗時的操作是求GX*GX+GY*GY的平方。

上面這段代碼，在不打開編譯器的SSE優化和快速浮點計算的情況，直接使用FPU，對4000*3000的彩色圖約需要480ms，當開啟SSE後，大概時間為220ms ，是以系統編譯器的SSE優化也很厲害，反編譯後可以看到彙編裡這樣的部分：

59AD12F8  movd        xmm0,ecx  
59AD12FC  cvtdq2ps    xmm0,xmm0  
59AD12FF  sqrtss      xmm0,xmm0  
59AD1303  cvttss2si   ecx,xmm0

　　可見他是調用的單浮點數的sqrt優化。

　由于該Sobel的過程最後是把資料用圖像的方式記錄下來，是以，IM_ClampToByte(sqrtf(GX * GX + GY * GY + 0.0F))可以用查表的方式來實作。簡單改成如下的版本，避免了浮點計算。

int IM_SobelTable(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) && (Channel != 3))            return IM_STATUS_INVALIDPARAMETER;

    unsigned char *RowCopy = (unsigned char *)malloc((Width + 2) * 3 * Channel);
    if (RowCopy == NULL)    return IM_STATUS_OUTOFMEMORY;

    unsigned char *First = RowCopy;
    unsigned char *Second = RowCopy + (Width + 2) * Channel;
    unsigned char *Third = RowCopy + (Width + 2) * 2 * Channel;

    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 + Stride, Channel);                                                    //    拷貝第二行資料
    memcpy(Third + Channel, Src + Stride, Width * Channel);
    memcpy(Third + (Width + 1) * Channel, Src + Stride + (Width - 1) * Channel, Channel);

    int BlockSize = 16, Block = (Width * Channel) / BlockSize;

    unsigned char Table[65026];
    for (int Y = 0; Y < 65026; Y++)        Table[Y] = (sqrtf(Y + 0.0f) + 0.5f);

    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        unsigned char *LinePD = Dest + Y * Stride;
        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) * Stride, Channel);
            memcpy(Third + Channel, Src + (Y + 1) * Stride, Width * Channel);                            //    由于備份了前面一行的資料，這裡即使Src和Dest相同也是沒有問題的
            memcpy(Third + (Width + 1) * Channel, Src + (Y + 1) * Stride + (Width - 1) * Channel, Channel);
        }
        if (Channel == 1)
        {
            for (int X = 0; X < Width; X++)
            {
                int GX = First[X] - First[X + 2] + (Second[X] - Second[X + 2]) * 2 + Third[X] - Third[X + 2];
                int GY = First[X] + First[X + 2] + (First[X + 1] - Third[X + 1]) * 2 - Third[X] - Third[X + 2];
                LinePD[X] = Table[IM_Min(GX * GX + GY * GY, 65025)];
            }
        }
        else
        {
            for (int X = 0; X < Width * 3; X++)
            {
                int GX = First[X] - First[X + 6] + (Second[X] - Second[X + 6]) * 2 + Third[X] - Third[X + 6];
                int GY = First[X] + First[X + 6] + (First[X + 3] - Third[X + 3]) * 2 - Third[X] - Third[X + 6];
                LinePD[X] = Table[IM_Min(GX * GX + GY * GY, 65025)];
            }
        }
    }
    free(RowCopy);
    return IM_STATUS_OK;
}

　　對4000*3000的彩色圖約需要180ms，比系統的SSE優化快了40ms，而這個過程完全無浮點計算，是以，可以知道計算GX和GY的耗時在本例中也占據了相當大的部分。

　　這樣的過程最适合于SSE處理了。

　　我們分析之。

　　第一來看一看這兩句：

int GX = First[X] - First[X + 2] + (Second[X] - Second[X + 2]) * 2 + Third[X] - Third[X + 2];
  int GY = First[X] + First[X + 2] + (First[X + 1] - Third[X + 1]) * 2 - Third[X] - Third[X + 2];

裡面涉及到了8個不同的像素，考慮計算的特性和資料的範圍，在内部計算時這個int可以用short代替，也就是要把加載的位元組圖像資料轉換為short類型先，這樣SSE優化方式為用8個SSE變量分别記錄8個連續的像素風量的顔色值，每個顔色值用16位資料表達。

　　這可以使用_mm_unpacklo_epi8配合_mm_loadl_epi64實作：

__m128i FirstP0 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(First + X)), Zero);
    __m128i FirstP1 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(First + X + 3)), Zero);
    __m128i FirstP2 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(First + X + 6)), Zero);

    __m128i SecondP0 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(Second + X)), Zero);
    __m128i SecondP2 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(Second + X + 6)), Zero);

    __m128i ThirdP0 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(Third + X)), Zero);
    __m128i ThirdP1 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(Third + X + 3)), Zero);
    __m128i ThirdP2 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(Third + X + 6)), Zero);

　　接着就是搬積木了，用SSE的指令代替普通的C的函數指令實作GX和GY的計算。

__m128i GX16 = _mm_abs_epi16(_mm_add_epi16(_mm_add_epi16(_mm_sub_epi16(FirstP0, FirstP2), _mm_slli_epi16(_mm_sub_epi16(SecondP0, SecondP2), 1)), _mm_sub_epi16(ThirdP0, ThirdP2)));
    __m128i GY16 = _mm_abs_epi16(_mm_sub_epi16(_mm_add_epi16(_mm_add_epi16(FirstP0, FirstP2), _mm_slli_epi16(_mm_sub_epi16(FirstP1, ThirdP1), 1)), _mm_add_epi16(ThirdP0, ThirdP2)));

　　找個時候的GX16和GY16裡儲存的是8個16位的中間結果，由于SSE隻提供了浮點數的sqrt操作，我們必須将它們轉換為浮點數，那麼這個轉換的第一步就必須是先将它們轉換為int的整形數，這樣，就必須一個拆成2個，即：

__m128i GX32L = _mm_unpacklo_epi16(GX16, Zero);
    __m128i GX32H = _mm_unpackhi_epi16(GX16, Zero);
    __m128i GY32L = _mm_unpacklo_epi16(GY16, Zero);
    __m128i GY32H = _mm_unpackhi_epi16(GY16, Zero);

　　接着又是搬積木了：

__m128i ResultL = _mm_cvtps_epi32(_mm_sqrt_ps(_mm_cvtepi32_ps(_mm_add_epi32(_mm_mullo_epi32(GX32L, GX32L), _mm_mullo_epi32(GY32L, GY32L)))));
    __m128i ResultH = _mm_cvtps_epi32(_mm_sqrt_ps(_mm_cvtepi32_ps(_mm_add_epi32(_mm_mullo_epi32(GX32H, GX32H), _mm_mullo_epi32(GY32H, GY32H)))));

　　整形轉換為浮點數，最後計算完之後又要将浮點數轉換為整形數。

　　最後一步，得到8個int型的結果，這個結果有要轉換為位元組類型的，并且這些資料有可能會超出位元組所能表達的範圍，是以就需要用到SSE自帶的抗飽和向下打包函數了，如下所示：

_mm_storel_epi64((__m128i *)(LinePD + X), _mm_packus_epi16(_mm_packus_epi32(ResultL, ResultH), Zero));

　　Ok，一切搞定了，還有一些細節處理自己慢慢補充吧。

　　運作，哇，隻要37ms了，速度快了N倍，可結果似乎和其他方式實作的不一樣啊，怎麼回事。

　　我也是找了半天都沒有找到問題所在，後來一步一步的測試，最終問題定位在16位轉換為32位整形那裡去了。

　　通常，我們都是對像素的位元組資料進行向上擴充，他們都是正數，是以用unpack之類的配合zero把高8位或高16位的資料填充為0就可以了，但是在本例中，GX16或者GY16很有可能是負數，而負數的最高位是符号位，如果都填充為0，則變為正數了，明顯改變原始的資料了，是以得到了錯誤的結果。

　　那如何解決問題呢，對于本例，很簡單，因為後面隻有一個平方操作，是以，對GX先取絕對值是不會改變計算的結果的，這樣就不會出現負的資料了，修改之後，果然結果正确。

　　還可以繼續優化，我們來看最後的計算GX*GX+GY*GY的過程，我們知道，SSE3提供了一個_mm_madd_epi16指令，其作用為：

r0 := (a0 * b0) + (a1 * b1)
r1 := (a2 * b2) + (a3 * b3)
r2 := (a4 * b4) + (a5 * b5)
r3 := (a6 * b6) + (a7 * b7)

如果我們能把GX和GY的資料拼接成另外兩個資料：

　　　　GXYL = GX0　　GY0　　GX1　　GY1　　GX2　　GY2　　GX3　　GY3

　　　　GXYH = GX4　　GY4　　GX5　　GY5　　GX6　　GY6　　GX7　　GY7

　　那麼調用_mm_madd_epi16(GXYL ，GXYL )和_mm_madd_epi16(GXYH ，GXYH )不就是能得到和之前一樣的結果了嗎，而這個拼接SSE已經有現成的函數了即：

__m128i GXYL = _mm_unpacklo_epi16(GX16, GY16);
__m128i GXYH = _mm_unpackhi_epi16(GX16, GY16);

　　這樣就把原來需要的10個指令變為了4個指令，代碼簡潔了并且速度也更快了。

　　嘗試如此修改，整個的計算過程時間減少到了32ms左右。

　　另外，還有一個可以優化的地方就是借用 _mm_maddubs_epi16 函數實作像素之間的加減乘除和擴充。

　　這個函數的作用如下：

r0 := SATURATE_16((a0 * b0) + (a1 * b1))
r1 := SATURATE_16((a2 * b2) + (a3 * b3))
...
r7 := SATURATE_16((a14 * b14) + (a15 * b15))

　　他的第一個參數是16個無符号的位元組資料，第二個參數是16個有符号的char資料。

　　配合unpack使用類似上面的技術就可以一次性處理16個位元組的像素簡加減了，這樣做整個過程大概能再加速2ms，達到最終的30ms。

　　源代碼位址：https://files.cnblogs.com/files/Imageshop/Sobel.rar （其中的SSE代碼請按照本文的思路自行整理。）

　　https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar，這裡是一個我全部用SSE優化的圖像處理的Demo，有興趣的朋友可以看看。

SSE圖像算法優化系列九：靈活運用SIMD指令16倍提升Sobel邊緣檢測的速度（4000*3000的24位圖像時間由480ms降低到30ms）。

　　歡迎點贊和打賞。

SSE圖像算法優化系列九：靈活運用SIMD指令16倍提升Sobel邊緣檢測的速度（4000*3000的24位圖像時間由480ms降低到30ms）。

繼續閱讀

Spring Boot 二三事：WEB 應用消息推送的那點事

Python+opencv 高斯濾波

Delphi 彙編學習（二）--- 學習 SIMD 的痛苦

單資料多指令(SIMD) MMX SSE AVX概述SIMD資料類型簡介

【論文閱讀】Bi-Directional Cascade Network for Perceptual Edge Detection（2019）摘要1. 介紹2. 模型

【論文閱讀】Dense Extreme Inception Network: Towards a Robust CNN Model for Edge Detection（2020）1. 摘要2. 模型3. 實驗準備步驟4. 實驗結果5. 代碼

深度學習傳統CV算法——邊緣檢測算子總結以及新興邊緣算法

四.使用OpenCv采集圖像及圖像邊緣檢測

對應Intel SSE的android NEON

畢業設計論文-基于matlab的車牌識别系統

Image Processing and Analysis_21_Scale Space

數字圖像處理實驗（2）——中值濾波/均值濾波

如何快速計算圖像梯度、幅值以及梯度方向角 -- 使用SSE指令集

深入淺出各種邊緣檢測算子及其推導

實戰 | OpenCV實作紡織物缺陷檢測->髒污、油漬、線條破損(詳細步驟 + 代碼)

OpenCV-Python 筆記（四）Canny邊緣檢測、圖像金字塔