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

時間 2019-11-06

標籤 sse 圖像算法優化系列二十六時間賽跑高斯金字塔創建計算過程简体版

原文原文鏈接

　　圖像金字塔技術在不少層面上都有着普遍的應用，不少開源的工具也都有對他們的創建寫了專門的函數，好比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爲下標起點的行列）。在此過程當中使用到的高斯模糊的權重矩陣的形式以下所示：編程

　　應該說，上面的過程用僞代碼表示應該是:ide

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個份量能夠重複利用，另外兩個必須從新計算。優化

　　上面的代碼是先計算列方向，而後在計算行方向的，基本上不須要額外的內存，咱們再來試下先計算行方向的累積值，再處理列方向的方案：ui

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函數以下所示：spa

// 一行像素的下取樣算法，權重係數爲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行的，而無需總體更新。code

　咱們對這兩種方案進行了速度測試，因爲自己這個的執行速度就比較塊，所以咱們對算法進行了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