2. 程式人生 > >SSE影象演算法優化系列二十六:和時間賽跑之優化高斯金字塔建立的計算過程。

SSE影象演算法優化系列二十六:和時間賽跑之優化高斯金字塔建立的計算過程。

阿新 • • 發佈：2019-01-14

　　影象金字塔技術在很多層面上都有著廣泛的應用，很多開源的工具也都有對他們的建立寫了專門的函式，比如IPP，比如OpenCV等等，這方面的理論文章特別多，我不需要贅述，但是我發現大部多分開源的程式碼的實現都不是嚴格意義上的金字塔，而是做了一定的變通，這種變通常常為了快捷的實現類似的效果，雖然這種變通不太會影響金字塔的效果，但是我這裡希望從嚴格意義上對該演算法進行優化，比如簡要貼一下下面的某個高斯金字塔的程式碼：

    public static Mat[] build(Mat img, int level) {
        Mat[] gaussPyr = new Mat[level];
        Mat mask 
 
= filterMask(img);
        Mat tmp = new Mat();
        Imgproc.filter2D(img, tmp, -1, mask);
        gaussPyr[0] = tmp.clone();
        Mat tmpImg = img.clone();
        for (int i = 1; i < level; i++) {
            // resize image
            Imgproc.resize(tmpImg, tmpImg, new Size(), 0.5, 0.5, Imgproc.INTER_LINEAR);
            
 
// blur image
            tmp = new Mat();
            Imgproc.filter2D(tmpImg, tmp, -1, mask);
            gaussPyr[i] = tmp.clone();
        }
        return gaussPyr;
    }
    
    private static Mat filterMask(Mat img) {
        double[] h = { 1.0 / 16.0, 4.0 / 16.0, 6.0 / 16.0, 4.0 / 16.0, 1.0 / 16.0 };
        Mat mask 
 
= new Mat(h.length, h.length, img.type());
        for (int i = 0; i < h.length; i++) {
            for (int j = 0; j < h.length; j++) {
                mask.put(i, j, h[i] * h[j]);
            }
        }
        return mask;
    }

　　上面的過程是對原圖線進行雙線性的下采樣插值，然後再對取樣後的圖進行一定的高斯模糊。我們說這樣做未嘗不可，而真正的高斯金字塔的建立過程是：對上一級的資料進行高斯模糊（5x5)得到結果T，然後刪除T的奇數行和奇數列資料作為下一級的結果（以0為下標起點的行列）。在此過程中使用到的高斯模糊的權重矩陣的形式如下所示：  

        

　　應該說，上面的過程用偽程式碼表示應該是:

Imgproc.filter2D(tmpImg, tmp, -1, mask);
Imgproc.resize(tmp, tmp, new Size(), 0.5, 0.5, Imgproc.INTER_NEARESTNEIGHBOR);
　即先高斯模糊，然後在使用最近領插值縮小一半，和上面程式碼中的先雙性縮小一半，再進行高斯模糊還是有區別的。

  　從程式設計優化的角度考慮，我們沒有必要完整的對上一級進行高斯模糊，而只需要進行偶數行和偶數列的計算，然後賦值到下一層資料中，而進一步考慮上述矩陣的特殊性，可以通過減少重複計算以及合併相同係數項等手段來優化。對於邊緣（2行2列）的畫素，把他單獨提取出來，減少中間資料的邊緣判斷可進一步提高速度。

　　再提出第一版C語言程式碼之前，我們來要看下各層金字塔的大小問題，如果上一層的大小位偶數，則下一層直接就除以2，但是如果是奇數，則除2時需要向上取整，比如寬某層的寬度尺寸為101，則之後的各層依次為101->51->26->13->7->4->2->1。

　　再看看前面權重矩陣的式子最右邊那個乘法，那表示這個權重矩陣是行列可分離的，我們可以先計算行的加權，然後再利用這個加權值計算列的加權，也可以先計算列然後再計算行，這樣原本每個畫素處的25個乘法和多24次加法就可以減少為10次乘法和9次加法。對於沿著寬度方向的更新計算，我們還可以充分利用列方向的重疊累計資訊，減少計算量，一個可行的簡單的程式碼如下所示：

int IM_DownSample8U_C1(unsigned char *Src, unsigned char *Dest, int SrcW, int SrcH, int StrideS, int DstW, int DstH, int StrideD, int Channel)
{
    if ((Channel != 1) && (Channel != 3) && (Channel != 4))        return IM_STATUS_NOTSUPPORTED;
    if (DstW != ((SrcW + 1) >> 1) || DstH != ((SrcH + 1) >> 1)) return IM_STATUS_INVALIDPARAMETER;            //    尺寸匹配
    int Status = IM_STATUS_OK;
    if (Channel == 1)
    {
        int Sum1, Sum2, Sum3, Sum4, Sum5;
        for (int Y = 1; Y < DstH - 1; Y++)                        //    不處理邊緣部分資料
        {
            unsigned char *LinePD = Dest + Y * StrideD ;
            unsigned char *P1 = Src + (Y * 2 - 2) * StrideS;
            unsigned char *P2 = P1 + StrideS;
            unsigned char *P3 = P2 + StrideS;
            unsigned char *P4 = P3 + StrideS;
            unsigned char *P5 = P4 + StrideS;
            Sum3 = P1[0] + ((P2[0] + P4[0]) << 2) + P3[0] * 6 + P5[0];
            Sum4 = P1[1] + ((P2[1] + P4[1]) << 2) + P3[1] * 6 + P5[1];
            Sum5 = P1[2] + ((P2[2] + P4[2]) << 2) + P3[2] * 6 + P5[2];
            for (int X = 1; X < DstW - 1; X++)
            {
                Sum1 = Sum3;    Sum2 = Sum4;    Sum3 = Sum5;
                Sum4 = P1[3] + ((P2[3] + P4[3]) << 2) + P3[3] * 6 + P5[3];
                Sum5 = P1[4] + ((P2[4] + P4[4]) << 2) + P3[4] * 6 + P5[4];
                LinePD[X] = (Sum1 + ((Sum2 + Sum4) << 2) + Sum3 * 6 + Sum5 + 128) >> 8;        //    注意四捨五入
                P1 += 2;    P2 += 2;    P3 += 2;    P4 += 2;    P5 += 2;
            }
        }
    }
    else if (Channel == 3)
    {

    }
    else if (Channel == 4)
    {

    }

    for (int X = 0; X < DstW; X++)
    {
        //IM_DownSamplePerPixel8U(Src, Dest, SrcW, SrcH, StrideS, DstW, DstH, StrideD, Channel, X, 0);                    //    第一行及最後一行
        //IM_DownSamplePerPixel8U(Src, Dest, SrcW, SrcH, StrideS, DstW, DstH, StrideD, Channel, X, DstH - 1);
    }
    for (int Y = 1; Y < DstH - 1; Y++)
    {
        //IM_DownSamplePerPixel8U(Src, Dest, SrcW, SrcH, StrideS, DstW, DstH, StrideD, Channel, 0, Y);                    //    第一列及最後一列
        //IM_DownSamplePerPixel8U(Src, Dest, SrcW, SrcH, StrideS, DstW, DstH, StrideD, Channel, DstW - 1, Y);
    }
　　return IM_STATUS_OK;
}

 　　注意到，在單通道的示例裡，我們用了5箇中間變數，分別記錄了某個位置5列畫素的累加和，在移動到下一個目標畫素時，由於我們是隔行隔列取樣的，因此移動到下一個畫素時只有3個位置時重疊的，也就是說只有3個分量可以重複利用，另外兩個必須重新計算。

 　　上面的程式碼是先計算列方向，然後在計算行方向的，基本上不需要額外的記憶體，我們再來試下先計算行方向的累積值，再處理列方向的方案：

int IM_DownSample8U_C2(unsigned char *Src, unsigned char *Dest, int SrcW, int SrcH, int StrideS, int DstW, int DstH, int StrideD, int Channel)
{
    if ((Channel != 1) && (Channel != 3) && (Channel != 4))        return IM_STATUS_NOTSUPPORTED;
    if (DstW != ((SrcW + 1) >> 1) || DstH != ((SrcH + 1) >> 1)) return IM_STATUS_INVALIDPARAMETER;            //    尺寸匹配
    int Status = IM_STATUS_OK;
    unsigned short *Sum0 = (unsigned short *)malloc(DstW * Channel * sizeof(unsigned short));
    unsigned short *Sum1 = (unsigned short *)malloc(DstW * Channel * sizeof(unsigned short));
    unsigned short *Sum2 = (unsigned short *)malloc(DstW * Channel * sizeof(unsigned short));
    unsigned short *Sum3 = (unsigned short *)malloc(DstW * Channel * sizeof(unsigned short));
    unsigned short *Sum4 = (unsigned short *)malloc(DstW * Channel * sizeof(unsigned short));
    if ((Sum0 == NULL) || (Sum1 == NULL) || (Sum2 == NULL) || (Sum3 == NULL) || (Sum4 == NULL))
    {
        Status = IM_STATUS_OUTOFMEMORY;
        goto FreeMemroy;
    }
    
    for (int Y = 1; Y < DstH - 1; Y++)                        //    不處理邊緣部分資料
    {
        unsigned char *LinePD = Dest + Y * StrideD;
        if (Y == 1)
        {
            IM_DownSampleLine8U_C(Src + 0 * StrideS, Sum0, DstW, Channel);
            IM_DownSampleLine8U_C(Src + 1 * StrideS, Sum1, DstW, Channel);
            IM_DownSampleLine8U_C(Src + 2 * StrideS, Sum2, DstW, Channel);
            IM_DownSampleLine8U_C(Src + 3 * StrideS, Sum3, DstW, Channel);
            IM_DownSampleLine8U_C(Src + 4 * StrideS, Sum4, DstW, Channel);
        }
        else
        {
            unsigned short *Temp1 = Sum0, *Temp2 = Sum1;
            Sum0 = Sum2;    Sum1 = Sum3;    Sum2 = Sum4;
            Sum3 = Temp1;    Sum4 = Temp2;
            IM_DownSampleLine8U_C(Src + (Y * 2 + 1) * StrideS, Sum3, DstW, Channel);
            IM_DownSampleLine8U_C(Src + (Y * 2 + 2) * StrideS, Sum4, DstW, Channel);
        }
        for (int X = Channel; X < (DstW - 1) * Channel; X++)
        {
            LinePD[X] = (Sum0[X] + ((Sum1[X] + Sum3[X]) << 2) + Sum2[X] * 6 + Sum4[X] + 127) >> 8;
        }
    }

    for (int X = 0; X < DstW; X++)
    {
        //IM_DownSamplePerPixel8U(Src, Dest, SrcW, SrcH, StrideS, DstW, DstH, StrideD, Channel, X, 0);                    //    第一行及最後一行
        //IM_DownSamplePerPixel8U(Src, Dest, SrcW, SrcH, StrideS, DstW, DstH, StrideD, Channel, X, DstH - 1);
    }

    for (int Y = 1; Y < DstH - 1; Y++)
    {
        //IM_DownSamplePerPixel8U(Src, Dest, SrcW, SrcH, StrideS, DstW, DstH, StrideD, Channel, 0, Y);                    //    第一列及最後一列
        //IM_DownSamplePerPixel8U(Src, Dest, SrcW, SrcH, StrideS, DstW, DstH, StrideD, Channel, DstW - 1, Y);
    }

FreeMemroy:
    if (Sum0 != NULL)        free(Sum0);
    if (Sum1 != NULL)        free(Sum1);
    if (Sum2 != NULL)        free(Sum2);
    if (Sum3 != NULL)        free(Sum3);
    if (Sum4 != NULL)        free(Sum4);
    return Status;
}

　　其中IM_DownSampleLine8U_C函式如下所示：

//    一行畫素的下取樣演算法，權重係數為1 4 6 4 1，並且是隔列取樣
void IM_DownSampleLine8U_C(unsigned char *Src, unsigned short *Dest, int DstLen, int Channel)
{
    //    只處理中間有效範圍內的數，第一個和最後一個不處理
    if (Channel == 1)
    {
        for (int Y = 1, Index = 2; Y < DstLen - 1; Y++, Index += 2)
        {
            Dest[Y] = Src[Index - 2] + ((Src[Index - 1] + Src[Index + 1] ) << 2) + Src[Index] * 6 + Src[Index + 2];
        }
    }
    else if (Channel == 3)            //    一個load語句可以包含5個完整畫素，可以處理1個目標畫素，3個分量佔6個位元組，不好儲存，因此，一次性處理2個目標畫素就好儲存了
    {
        
    }
    else if (Channel == 4)        //    4個通道的一次性只處理一個畫素的，需要訪問源影象20個位元組範圍
    {
        
    }
}

　　注意到上述程式碼在結構上和第一版本其實差不多，不過多了5行臨時記憶體，在更新行權重的時候也是隻需要更新2行的，而無需整體更新。

 　  我們對這兩種方案進行了速度測試，由於本身這個的執行速度就比較塊，因此我們對演算法進行了100次計算，對於第一級為1920*1080大小的灰度圖，下一級的高斯金字塔大小為960*540畫素，演算法C1測試的結果為267ms，演算法C2的執行速度約為256ms，這說明他們本質上不存在大的速度差異（這裡的時間都不包括處理邊緣畫素的，後面還要討論，並且高斯金字塔也是採用unsigned char型別資料來儲存的）。

       在某些場合，我們還需要加快這個過程的速度，因此我考慮使用SSE優化他，考慮以上兩種實現方式，哪一種更有利於SSE的處理呢，由於第一種方式前後的依賴比較強，用SSE做不是不可以，但估計效率不會有提升，需要太多次資料重組了，而第二種方式的由中間資料計算最後的結果很明顯可以使用SSE處理，即下面的這三行程式碼：

for (int X = Channel; X < (DstW - 1) * Channel; X++)
{
     LinePD[X] = (Sum0[X] + ((Sum1[X] + Sum3[X]) << 2) + Sum2[X] * 6 + Sum4[X] + 127) >> 8;
}

　　這裡的Sum是16位的，LinePD是8位的。替換的程式碼如下所示：

int BlockSize = 8, Block = (DstW - 1 - 1) * Channel / BlockSize;
for (int X = Channel; X < Block * BlockSize + Channel; X += BlockSize)
{
    __m128i S0 = _mm_loadu_si128((__m128i *)(Sum0 + X));
    __m128i S1 = _mm_loadu_si128((__m128i *)(Sum1 + X));
    __m128i S2 = _mm_loadu_si128((__m128i *)(Sum2 + X));
    __m128i S3 = _mm_loadu_si128((__m128i *)(Sum3 + X));
    __m128i S4 = _mm_loadu_si128((__m128i *)(Sum4 + X));
    __m128i Sum = _mm_add_epi16(_mm_add_epi16(_mm_add_epi16(_mm_add_epi16(S0, S4), _mm_slli_epi16(_mm_add_epi16(S1, S3), 2)), _mm_mullo_epi16(S2, _mm_set1_epi16(6))), _mm_set1_epi16(127));
    _mm_storel_epi64((__m128i *)(LinePD + X), _mm_packus_epi16(_mm_srli_epi16(Sum, 8), _mm_setzero_si128()));
}
for (int X = Block * BlockSize + Channel; X < (DstW - 1) * Channel; X++)
{
    LinePD[X] = (Sum0[X] + ((Sum1[X] + Sum3[X]) << 2) + Sum2[X] * 6 + Sum4[X] + 127) >> 8;
}

　　這麼個簡單的改動，速度大概到了180ms。

　  有點麻煩的是IM_DownSampleLine8U_C這個函式的優化，其核心程式碼如下：

  for (int Y = 1, Index = 2; Y < DstLen - 1; Y++, Index += 2)
  {
 　　Dest[Y] = Src[Index - 2] + ((Src[Index - 1] + Src[Index + 1] ) << 2) + Src[Index] * 6 + Src[Index + 2];
  }

　　最簡單的SSE處理方式是載入5次不同位置的Src值，然後將資料轉換為16位，再進行加法和乘法計算。

int ValidLen = DstLen - 2;                                //    -2是因為第一和最後一個點的部分取樣值是不在有效範圍內的
int BlockSize = 8, Block = ValidLen / BlockSize;
for (int Y = 1; Y < Block * BlockSize + 1; Y += BlockSize)
{
    int Index = (Y - 1) * 2;
    __m128i SrcV0 = _mm_cvtepu8_epi16(_mm_loadl_epi64((__m128i *)(Src + Index)));
    __m128i SrcV1 = _mm_cvtepu8_epi16(_mm_loadl_epi64((__m128i *)(Src + Index + 1)));
    __m128i SrcV2 = _mm_cvtepu8_epi16(_mm_loadl_epi64((__m128i *)(Src + Index + 2)));
    __m128i SrcV3 = _mm_cvtepu8_epi16(_mm_loadl_epi64((__m128i *)(Src + Index + 3)));
    __m128i SrcV4 = _mm_cvtepu8_epi16(_mm_loadl_epi64((__m128i *)(Src + Index + 4)));
    _mm_storeu_si128((__m128i *)(Dest + Y), _mm_add_epi16(_mm_add_epi16(SrcV0, SrcV1), _mm_add_epi16(_mm_slli_epi16(_mm_add_epi16(SrcV1, SrcV3), 2), _mm_mullo_epi16(SrcV2, _mm_set1_epi16(6)))));
}
for (int Y = Block * BlockSize + 1; Y < DstLen - 1; Y++)
{
    Dest[Y] = Src[Y * 2 - 2] + Src[Y * 2 - 1] * 4 + Src[Y * 2] * 6 + Src[Y * 2 + 1] * 4 + Src[Y * 2 + 2];
}

　　一次性處理8個畫素，需要多次載入記憶體，原以為速度會有問題，結果一測，速度居然飆升到40ms，單次只需要0.4ms了。真的很高興。

　　但是和普通的C比較一下，似乎結果不對啊，仔細分析，原來是因為這個對Src取樣計算時每次是隔一個點取一個樣的，而上述程式碼側連續取樣的，那怎麼辦呢？

　　也很簡單，我們使用_mm_loadu_si128一次性載入16個位元組，然後每隔一個畫素就置0，這樣就相當於把剩下的8個畫素的值直接變為了16位資料了，一舉兩得。如下所示：

        int ValidLen = DstLen - 2;                                //    -2是因為第一和最後一個點的部分取樣值是不在有效範圍內的
        int BlockSize = 8, Block = ValidLen / BlockSize;
        __m128i Mask = _mm_setr_epi8(255, 0, 255, 0, 255, 0, 255, 0, 255, 0, 255, 0, 255, 0, 255, 0);
        for (int Y = 1; Y < Block * BlockSize + 1; Y += BlockSize)
        {
            int Index = (Y - 1) * 2;
            __m128i SrcV0 = _mm_and_si128(_mm_loadu_si128((__m128i *)(Src + Index + 0)), Mask);
            __m128i SrcV1 = _mm_and_si128(_mm_loadu_si128((__m128i *)(Src + Index + 1)), Mask);
            __m128i SrcV2 = _mm_and_si128(_mm_loadu_si128((__m128i *)(Src + Index + 2)), Mask);
            __m128i SrcV3 = _mm_and_si128(_mm_loadu_si128((__m128i *)(Src + Index + 3)), Mask);
            __m128i SrcV4 = _mm_and_si128(_mm_loadu_si128((__m128i *)(Src + Index + 4)), Mask);
            _mm_storeu_si128((__m128i *)(Dest + Y), _mm_add_epi16(_mm_add_epi16(SrcV0, SrcV4), _mm_add_epi16(_mm_slli_epi16(_mm_add_epi16(SrcV1, SrcV3), 2), _mm_mullo_epi16(SrcV2, _mm_set1_epi16(6)))));
        }
        for (int Y = Block * BlockSize + 1; Y < DstLen - 1; Y++)
        {
            Dest[Y] = Src[Y * 2 - 2] + Src[Y * 2 - 1] * 4 + Src[Y * 2] * 6 + Src[Y * 2 + 1] * 4 + Src[Y * 2 + 2];
        }

　　我們上面用的and運算來將有關位置0，當然還可以使用 shuffle指令來得到同樣的結果，速度大概也就稍微慢一點，大概再45ms。

　　實際上，在這裡的由於權重有一些特殊性，比如有2個1，2個4，4還可以用移位實現，如果是一些其他不太有規律的權重，比如 3 7 9 7 3這種，我們實際上還有一種優化方式來處理，因為在SSE裡還有一個_mm_maddubs_epi16這個指令，他可以一次性實現16個位元組數*16個signed char，然後再兩兩相加儲存到8個short型別中去，比如上面的程式碼也可以用下面的方式實現：

        int ValidLen = DstLen - 2;                                //    -2是因為第一和最後一個點的部分取樣值是不在有效範圍內的
        int BlockSize = 6, Block = ValidLen / BlockSize;
        //    Src     S0    S1    S2    S3    S4    S5    S6    S7    S8    S9    S10    S11    S12    S13    S14    S15                //    16個畫素    
        //    Dst                D1        D2        D3        D4        D5        D6                            //    有效位置的只有6個結果

        //    1   4   6   4   1   4   6   4   1   4   6   4   1   4   6   4
        __m128i Cof = _mm_setr_epi8(1, 4, 6, 4, 1, 4, 6, 4, 1, 4, 6, 4, 1, 4, 6, 4);    //    用同一個係數不影響，因為後面反正拋棄了後半部分的累加        //__m128i Cof1 = _mm_setr_epi8(1, 4, 6, 4, 1, 4, 6, 4, 1, 4, 6, 4, 1, 4, 6, 4);
        //    1   4   6   4   1   4   6   4   1   4   6   4   1   4   6   4
        //__m128i Cof2 = _mm_setr_epi8(1, 4, 6, 4, 1, 4, 6, 4, 0, 0, 0, 0, 0, 0, 0, 0);
        for (int Y = 1; Y < Block * BlockSize + 1; Y += BlockSize)
        {
            //    S0    S1    S2    S3    S4    S5    S6    S7    S8    S9    S10    S11    S12    S13    S14    S15
            __m128i SrcV = _mm_loadu_si128((__m128i *)(Src + (Y - 1) * 2));
            //    S0    S1    S2    S3    S2    S3    S4    S5  S4    S5    S6    S7    S6    S7    S8    S9
            __m128i Src1 = _mm_shuffle_epi8(SrcV, _mm_setr_epi8(0, 1, 2, 3, 2, 3, 4, 5, 4, 5, 6, 7, 6, 7, 8, 9));
            //    S8    S9    S10    S11    S10    S11    S12    S13 S4    S6    S8    S10    S12    S14    0    0
            __m128i Src2 = _mm_shuffle_epi8(SrcV, _mm_setr_epi8(8, 9, 10, 11, 10, 11, 12, 13, 4, 6, 8, 10, 12, 14, -1, -1));
            //    S0 + S1 * 4        S2 * 6 + S3 * 4        S2 + S3 * 4        S4 * 6 + S5 * 4        S4 + S5 * 4        S6 * 6 + S7 * 4  S6 + S7 * 4        S8 * 6 + S9 * 4
            __m128i Dst1 = _mm_maddubs_epi16(Src1, Cof);    //     _mm_maddubs_epi16(Src1, Cof1);
            //    S8 + S9 * 4        S10 * 6 + S11 * 4    S10 + S11 * 4    S12 * 6 + S13 * 4        0        0        0        0        0        0        0        0
            __m128i Dst2 = _mm_maddubs_epi16(Src2, Cof);    //     _mm_maddubs_epi16(Src1, Cof2);
            //    S0 + S1 * 4    + S2 * 6 + S3 * 4        S2 + S3 * 4    + S4 * 6 + S5 * 4        S4 + S5 * 4    + S6 * 6 + S7 * 4    S6 + S7 * 4    + S8 * 6 + S9 * 4        S8 + S9 * 4    + S10 * 6 + S11 * 4        S10 + S11 * 4 + S12 * 6 + S13 * 4
            __m128i Dst12 = _mm_hadd_epi16(Dst1, Dst2);
            //    S0 + S1 * 4    + S2 * 6 + S3 * 4 + S4        S2 + S3 * 4    + S4 * 6 + S5 * 4 + S6        S4 + S5 * 4    + S6 * 6 + S7 * 4+ S8    S6 + S7 * 4    + S8 * 6 + S9 * 4 + S10        S8 + S9 * 4    + S10 * 6 + S11 * 4    + S12    S10 + S11 * 4 + S12 * 6 + S13 * 4 + S14
            __m128i Dst = _mm_add_epi16(Dst12, _mm_unpackhi_epi8(Src2, _mm_setzero_si128()));

            _mm_storeu_epi96((__m128i *)(Dest + Y), Dst);
        }
        
        for (int Y = Block * BlockSize + 1; Y < DstLen - 1; Y++)
        {
            Dest[Y] = Src[Y * 2 - 2] + Src[Y * 2 - 1] * 4 + Src[Y * 2] * 6 + Src[Y * 2 + 1] * 4 + Src[Y * 2 + 2];
        }

　　在本例中，上述程式碼執行100次要50幾毫秒，比前面的慢了，但是這裡的組合確實是蠻有味道的，各種資料的靈活應用也是值得參考的。我反而更欣賞這段程式碼。

 　  以上談及的均是單通道的演算法，如果是BGR 3個通道或者BGRA 4個通道的影象資料，情況就會複雜一些，但是同樣的道理，可以使用shuffle來調整位置，然後使用類似的方式處理。

　　我們來談談浮點版本的高斯金字塔，這個再很多情況下也有需求，畢竟有很多演算法是再浮點上進行處理的，那浮點版本的普通C的程式碼其實和C語言是差不多的，只需要將有關資料型別改為浮點就可以了，那對於其核心的DownSampleLine函式，也是我們優化的關鍵和難點，由於SSE一次性只能載入4個浮點數，如果還是和剛才處理位元組資料那樣，隔一個數取一個數，那麼利用SSE一次性只能處理2個畫素，而我們通過下面的美好的優化方式，一次性就能處理4個畫素了，而且程式碼也很優美，我很是喜歡。

void IM_DownSampleLine32F(float *Src, float *Dest, int DstLen, int Channel)
{
    //    只處理中間有效範圍內的數，第一個和最後一個不處理
    if (Channel == 1)
    {
        int ValidLen = DstLen - 2;                                                //    -2是因為第一和最後一個點的部分取樣值是不在有效範圍內的
        int BlockSize = 4, Block = ValidLen / BlockSize;
        //    Src    S0    S1    S2    S3    S4    S5    S6    S7    S8    S9    S10    S11   　　//    12個數據    
        //    Dst                D1           D2          D3         D4                            　　//    有效位置的只有4個結果

        //    1   4   6   4  
        __m128 Cof = _mm_setr_ps(1.0f, 4.0f, 6.0f, 4.0f);      
        for (int Y = 1; Y < Block * BlockSize + 1; Y += BlockSize)
        {
            //    S0    S1    S2    S3
            __m128 SrcV0 = _mm_loadu_ps(Src + (Y - 1) * 2);
            //    S4    S5    S6    S7
            __m128 SrcV1 = _mm_loadu_ps(Src + (Y - 1) * 2 + 4);
            //    S8    S9    S10    S11
            __m128 SrcV2 = _mm_loadu_ps(Src + (Y - 1) * 2 + 8);            //    下一次載入時的SrcV0和本次的SrcV2時相同的，測試過用變數賦值，結果區別不大
            //    S0    S1 * 4  S2 * 6  S3 * 4
            __m128 Sum0 = _mm_mul_ps(SrcV0, Cof);
            //    S2  S3 * 4  S4 * 6  S5 * 4
            __m128 Sum1 = _mm_mul_ps(_mm_shuffle_ps(SrcV0, SrcV1, _MM_SHUFFLE(1, 0, 3, 2)), Cof);
            //    S4  S5 * 4  S6 * 6  S7 * 4
            __m128 Sum2 = _mm_mul_ps(SrcV1, Cof);
            //    S6  S7 * 4  S8 * 6  S9 * 4
            __m128 Sum3 = _mm_mul_ps(_mm_shuffle_ps(SrcV1, SrcV2, _MM_SHUFFLE(1, 0, 3, 2)), Cof);
            //    S0 + S1 * 4 + S2 * 6 + S3 * 4            S2 + S3 * 4 + S4 * 6 + S5 * 4        S4 + S5 * 4 + S6 * 6 + S7 * 4        S6 + S7 * 4 + S8 * 6 + S9 * 4
            __m128 Dst = _mm_hadd_ps(_mm_hadd_ps(Sum0, Sum1), _mm_hadd_ps(Sum2, Sum3));
            //    S0 + S1 * 4 + S2 * 6 + S3 * 4 + S4            S2 + S3 * 4 + S4 * 6 + S5 * 4 + S6        S4 + S5 * 4 + S6 * 6 + S7 * 4 + S8        S6 + S7 * 4 + S8 * 6 + S9 * 4 + S10
            Dst = _mm_add_ps(Dst, _mm_shuffle_ps(SrcV1, SrcV2, _MM_SHUFFLE(2, 0, 2, 0)));
            _mm_storeu_ps(Dest + Y, Dst);
        }
        for (int Y = Block * BlockSize + 1; Y < DstLen - 1; Y++)
        {
            Dest[Y] = Src[Y * 2 - 2] + Src[Y * 2 - 1] * 4 + Src[Y * 2] * 6 + Src[Y * 2 + 1] * 4 + Src[Y * 2 + 2];
        }
    }
}

　　具體的每句程式碼的意思根據我上面的註釋應該很容易能弄懂的，我就不加解釋了。

　　最後，我們來關注下邊緣的處理，邊緣部分由於取樣時會超出影象邊界，因此，需要做判斷，一種合理的方式是採用映象資料，此時可以保證權重一定是256，我做了一個簡單的函式：

int IM_DownSamplePerPixel8U(unsigned char *Src, unsigned char *Dest, int SrcW, int SrcH, int StrideS, int DstW, int DstH, int StrideD, int Channel, int X, int Y)
{
    if ((Channel != 1) && (Channel != 3) && (Channel != 4))        return IM_STATUS_NOTSUPPORTED;
    if (X < 0 || X >= DstW || Y < 0 || Y >= DstH)                return IM_STATUS_INVALIDPARAMETER;

    int Sum, SumB, SumG, SumR, SumA;
    const int Kernel[25] = {
        1, 4, 6, 4, 1,
        4, 16, 24, 16, 4,
        6, 24, 36, 24, 6,
        4, 16, 24, 16, 4,
        1, 4, 6, 4, 1 };                                                //    高斯卷積核

    Sum = SumB = SumG = SumR = SumA = 128;                                        //     128是用於四捨五入的
    for (int J = -2; J <= 2; J++)
    {
        int YY = IM_GetMirrorPos_Ex(SrcH, Y * 2 + J);                            //    YY = Clamp(Y * 2 + J, 0, SrcH - 1);        在上一級中的Y座標要乘以2, 使用Clamp速度會稍微慢一點
        for (int I = -2; I <= 2; I++)
        {
            int Weight = Kernel[(J + 2) * 5 + (I + 2)];                        //    嚴格的按照卷積公式進行計算，沒必要用中間變數去優化他了
            int XX = IM_GetMirrorPos_Ex(SrcW, X * 2 + I);                        //    XX = Clamp(X * 2 + I, 0, SrcW - 1); 用EX這個版本的保證邊緣的部分不會有重複的畫素，不然取樣就不對了
            unsigned char *Sample = Src + YY * StrideS + XX * Channel;
            if (Channel == 1)
            {
                Sum += Sample[0] * Weight;
            }
            else if (Channel == 3)
            {

            }
            else if (Channel == 4)
            {

            }
        }
    }
    int Index = Y * StrideD + X * Channel;

    if (Channel == 1)
    {
        Dest[Index] = Sum >> 8;
    }
    else if (Channel == 3)
    {

    }
    else if (Channel == 4)
    {
    }
    return IM_STATUS_OK;
}

　　為了程式完整，我們最後再加上週邊畫素的處理，然而我們發現一個嚴重的問題，再沒有處理四周的函式中，我們執行100次SSE的耗時大約是45ms,一旦加入邊緣畫素的處理，這個耗時我們發現75ms，而普通C語言版本里由原來的260ms變為290ms，我們可能感受不到大的區別，但SSE的優化後，邊緣部分居然佔用了40%的耗時，因此，此時邊緣特殊畫素的處理就成了核心的事情了。

　　一種可行的優化方式就是類似於我前面做的Sobel邊緣檢測時方式，先對資料進行擴充套件，然後對擴充套件後的資料進行處理，此時邊緣部分的處理已經被包括到SSE裡去了，我尚未實踐此方案的可行性和速度效果，相信應該不成問題。

　　附本文相關工程程式碼供參考：https://files.cnblogs.com/files/Imageshop/GaussPyramid.rar

 

 

 

相關推薦

SSE影象演算法優化系列二十六:和時間賽跑之優化高斯金字塔建立的計算過程。

　　影象金字塔技術在很多層面上都有著廣泛的應用，很多開源的工具也都有對他們的建立寫了專門的函式，比如IPP，比如OpenCV等等，這方面的理論文章特別多，我不需要贅述，但是我發現大部多分開源的程式碼的實現都不是嚴格意義上的金字塔，而是做了一定的變通，這種變通常常為了快捷的實現類似的效果，雖然這種變通不太會影響

SSE影象演算法優化系列二十四: 基於形態學的影象後期抗鋸齒演算法--MLAA優化研究。

       偶爾看到這樣的一個演算法，覺得還是蠻有意思的，花了將近10天多的時間研究了下相關程式碼。        以下為百度的結果：MLAA全稱Morphological Antialiasing，意為形態抗鋸齒是AMD推出

SSE影象演算法優化系列二十五:二值影象的Euclidean distance map（EDM)特徵圖計算及其優化。 SSE影象演算法優化系列九：靈活運用SIMD指令16倍提升Sobel邊緣檢測的速度（4000*3000的24點陣圖像時間由480ms降低到30ms）

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

SSE圖像算法優化系列二十四: 基於形態學的圖像後期抗鋸齒算法--MLAA優化研究。

遊戲 targe component 數據類型 ice computing walk 處理方式 soft 偶爾看到這樣的一個算法，覺得還是蠻有意思的，花了將近10天多的時間研究了下相關代碼。 以下為百度的結果：MLAA全稱Morphologic

演算法題（二十六）利用堆排序解決找出最小的k個值問題

題目描述 輸入n個整數，找出其中最小的K個數。例如輸入4,5,1,6,2,7,3,8這8個數字，則最小的4個數字是1,2,3,4,。 分析 問題很簡單，升序排序後直接輸出前k個，不過要考慮時間複雜度的問題。可以用堆排序，構建有k個值的大頂堆，然後用堆頭部與其他值比較，堆

二十六、閒談Tomcat效能優化

Tomcat在各位JavaWeb從業者常常就是預設的開發環境，但是Tomcat的預設配置作為生產環境，尤其是記憶體和執行緒的配置，預設都很低，容易成為效能瓶頸. 幸好Tomcat還有很多的提升空間.下文介紹一下Tomcat優化,可以分為記憶體,執行緒,IO. 一:Tomcat記憶體優

Expo大作戰(二十六)--expo sdk api之Video和WebBrowser

ast download 公眾 根據 pen 正在 但我 AI volume 簡要：本系列文章講會對expo進行全面的介紹，本人從2017年6月份接觸expo以來，對expo的研究斷斷續續，一路走來將近10個月，廢話不多說，接下來你看到內容，講全部來與官網

二十六、python學習之前端(九):JQuery資料互動

一.juqery物件: js中的物件,就是python中的字典: 1定義方法: 定義方法1: var obj1 = {name:"張三", "age":18, "address":"三里屯"} 定義方法2: var obj2 = new Object(); ob

Jmeter （二十六）邏輯控制器 之 Module Controller and Include Controller

Module Controller ---模組控制器 　　測試計劃設定“獨立執行沒每個執行緒組” 　　 　　執行緒組2中使用Module Controller執行執行緒組1中的Sampler： 　　 　　緊接著，將執行緒組1disable掉，執行： 　　 　　 　　結果樹中是從執行緒

效能測試二十六：環境部署之Mysql+Redis+Tomcat環境整合

系統中使用了快取+資料庫，通用讀取資料規則1、先從快取讀資料，如果有，直接返回資料；2、如果沒有，去資料庫中讀，然後再插入到快取中，再返回資料 Mysql+Redis+Tomcat環境整合 1、修改PerfTeach/WEB-INF/classes/redis-config.properties2、將re

性能測試二十六：環境部署之Mysql+Redis+Tomcat環境整合

web-inf 技術 提高 tar 從數據 之前 sse 取數 最小 系統中使用了緩存+數據庫，通用讀取數據規則1、先從緩存讀數據，如果有，直接返回數據；2、如果沒有，去數據庫中讀，然後再插入到緩存中，再返回數據 Mysql+Redis+Tomcat環境整合 1、修改Per

Java學習系列(二十一)Java面向物件之註解詳解

一、前言 Java中的註解Annotation運用到很多方面，比如之前講的單元測試中的@Test、Spring、SpringMVC還有其他框架等等。Java本身自帶的註解也有一些，比如：@Override(子類要重寫/覆寫父類的對應方法)、@Deprecated(表

SSE影象演算法優化系列二十三: 基於value-and-criterion structure 系列濾波器（如Kuwahara，MLV，MCV濾波器）的優化。 SSE影象演算法優化系列十四：區域性均方差及區域性平方差演算法的優化 SSE影象演算法優化系列七：基於SSE實現的極速的矩形核腐蝕和膨脹（

       基於value-and-criterion structure方式的實現的濾波器在原理上其實比較簡單，感覺下面論文中得一段話已經描述的比較清晰了，直接貼英文吧，感覺翻譯過來反而失去了原始的韻味了。        T

SSE影象演算法優化系列二十三: 基於value-and-criterion structure 系列濾波器（如Kuwahara，MLV，MCV濾波器）的優化。

int IM_KuwaharaFilter(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Radius) { int Channel = Stride / Width; if

SSE影象演算法優化系列三十：GIMP中的Noise Reduction演算法原理及快速實現。

　　GIMP原始碼連結：https://gitlab.gnome.org/GNOME/gimp/-/archive/master/gimp-master.zip 　　GEGL相關程式碼連結：https://gitlab.gnome.org/GNOME/gegl/-/archive/master/gegl-m

SSE圖像算法優化系列二十三: 基於value-and-criterion structure 系列濾波器（如Kuwahara，MLV，MCV濾波器）的優化。

Once 領域 hat ble 感覺 das 噪音 gre 種類 基於value-and-criterion structure方式的實現的濾波器在原理上其實比較簡單，感覺下面論文中得一段話已經描述的比較清晰了，直接貼英文吧，感覺翻譯過來反而失去了原始的韻味了

演算法題（三十六）：對稱的二叉樹

題目描述 請實現一個函式，用來判斷一顆二叉樹是不是對稱的。注意，如果一個二叉樹同此二叉樹的映象是同樣的，定義其為對稱的。 分析 從上而下遞迴掃描，只要發現左右結點不相同的返回false。 程式碼 public class Symmetry { public static

Django學習筆記（二十六）:閱讀計數和優化計數（防重新整理）

閱讀計數主要就是記錄網頁被瀏覽的次數，防刷的方法就是記錄瀏覽的cookie，使每個cookie只能增加一次閱讀數， 當然，刪除cookie和關閉瀏覽器就可以重新增加閱讀數，更加優化的方法還在學習中，先分享這個我會的 models.py的程式碼 class Books(m

Spring Cloud系列（二十六） Zuul過濾器詳解（Finchley.RC2版本）

Zuul核心過濾器 通過前兩篇對Zuul的介紹，我們對於Zuul的有了一個大概的印象：它包含了對請求的路由和過濾兩個功能，其中路由功能負責將外部請求轉發到具體的微服務例項上，是實現外部訪問統一入口的基礎；而過濾器功能則負責對請求的處理過程進行干預，是實現請求校驗、服務聚

TiDB 原始碼閱讀系列文章（二十一）基於規則的優化 II

在 TiDB 原始碼閱讀系列文章（七）基於規則的優化 一文中，我們介紹了幾種 TiDB 中的邏輯優化規則，包括列剪裁，最大最小消除，投影消除，謂詞下推和構建節點屬性，本篇將繼續介紹更多的優化規則：聚合消除、外連線消除和子查詢優化。 聚合消除 聚合消除會檢查 SQL 查詢中 Group By 語句所使用的列是否