轉置操做在不少算法上都有着普遍的應用,在數學上矩陣轉置更有着特殊的意義。而在圖像處理上,若是說圖像數據自己的轉置,除了顯示外,自己並沒有特殊含義,可是在某些狀況下,確能有效的提升算法效率,好比不少行列可分離的算法,在不少狀況下,行和列方向的算法邏輯隨相同,可是因爲多方面緣由(好比Cache miss, 優化水平等)行列處理時間仍是由很大的差別的,這個時候若是轉置的耗時和處理時間相比所佔比例甚小,則能夠考慮在進行耗時處理前先轉置數據,而後調用不耗時的方向的算法,處理完後再次進行轉置。所以,一個高效的圖像轉置算法的設計時很是有必要的。算法
沒怎麼蒐集這方面的資料,不過在百度上看到的優化的帖子也有幾篇:緩存
一、利用SSE優化圖像轉置 這篇文章講到了SSE優化轉置操做,講的很簡單,我只是稍微看了下他的代碼,他彷佛處理的不是普通的8位圖像,而是16位的,反正我是沒有看懂,而且他的提供比較的C代碼自己寫法就徹底沒有考慮到C語言自身的優化,所以最後提出SSE代碼比C快5倍說服力就大爲打折扣,不過他這裏能夠值得學習的地方就是這個轉置支持In-Place操做,就是Src和Dest能夠相同。網絡
二、圖像轉置的Neon優化代碼 Neon的代碼,沒看懂,不事後面說10倍左右的提速,其實也要看原始的C代碼是怎麼寫的了,不過原文也明確的說,只支持RGBA 32位的圖,顯然做者也避而不談灰度或者24位的圖,固然這於手機端彷佛沒有24位的概念有關。ide
三、CUDA學習筆記一:CUDA+OpenCV的圖像轉置,採用Shared Memory進行CUDA程序優化 這個文章是說GPU的優化,不過最後給出的GPU時間和CPU相比真的很慘。函數
也就是說國內網絡上的優化文章其實都仍是停留在皮毛階段,也或者是真正的具備優化意義的代碼都還雪藏在某個公司或者某我的的硬盤裏,特別是針對灰度和24位圖像的轉置優化在PC上有更多的使用場景。學習
3、個人貢獻測試
普通的C代碼的轉置很簡單,也曾嘗試過各類優化方案,可是最後都無啥特別大的改進,所以考慮使用SSE的方案。優化
(1)、由奢入儉難啊,咱們先挑最簡單來實現,我說的最簡單的是32位圖像。spa
32位圖像由B/G/R/A 4個份量組成,轉置時咱們須要把他們當作一個總體,以4*4大小的轉置威力,以下所示:.net
A0 A1 A2 A3 A0 B0 C0 D0
B0 B1 B2 B3 A1 B1 C1 D1
C0 C1 C2 C3 -------------------> A2 B2 C2 D2
D0 D1 D2 D3 A3 B3 C3 D3
其中每個元素都有4個字節份量組成。
看到這個若是用過SSE的朋友都會想起_MM_TRANSPOSE4_PS這個宏,他已經實現了4*4數據的轉置,可是仔細去看,這個是針對浮點數的一個宏,那好,咱們能夠直接看內部的實現,能夠發現內部主要是_mm_shuffle_ps的靈活應用,很明顯SSE針對長整型也有一整套的shuffle函數,對應的是_mm_shuffle_epi32,可以讓我想不明白的就是_mm_shuffle_ps能夠對兩個__m128數據進行混合shuffle,可是整形的只能處理一個__m128i數據的shuffle,由於沒法把這個宏的算法直接轉移到整形下。
最近對SSE的一些組合和拆解函數有了較爲專一的理解和測試,實現上述功能也沒有耗費我多少時間:
// BFRA的轉置,彷佛作成8*8的並無速度優點 void Transpose4x4_I(int *Src, int *Dest, int WidthS, int WidthD) { __m128i S0 = _mm_loadu_si128((__m128i *)(Src + 0 * WidthS)); // A3 A2 A1 A0 __m128i S1 = _mm_loadu_si128((__m128i *)(Src + 1 * WidthS)); // B3 B2 B1 B0 __m128i S01L = _mm_unpacklo_epi32(S0, S1); // B1 A1 B0 A0 __m128i S01H = _mm_unpackhi_epi32(S0, S1); // B3 A3 B2 A2 __m128i S2 = _mm_loadu_si128((__m128i *)(Src + 2 * WidthS)); // C3 C2 C1 C0 __m128i S3 = _mm_loadu_si128((__m128i *)(Src + 3 * WidthS)); // D3 D2 D1 D0 __m128i S23L = _mm_unpacklo_epi32(S2, S3); // D1 C1 D0 C0 __m128i S23H = _mm_unpackhi_epi32(S2, S3); // D3 C3 D2 C2 _mm_storeu_si128((__m128i *)(Dest + 0 * WidthD), _mm_unpacklo_epi64(S01L, S23L)); // D0 C0 B0 A0 _mm_storeu_si128((__m128i *)(Dest + 1 * WidthD), _mm_unpackhi_epi64(S01L, S23L)); // D1 C1 B1 A1 _mm_storeu_si128((__m128i *)(Dest + 2 * WidthD), _mm_unpacklo_epi64(S01H, S23H)); // D2 C2 B2 A2 _mm_storeu_si128((__m128i *)(Dest + 3 * WidthD), _mm_unpackhi_epi64(S01H, S23H)); // D3 C3 B3 A3 }
上面的代碼彷佛也不須要多作什麼解釋,能看懂我後面註釋的組合順序、能百度MSDN查每一個Intirsic指令的意義就能搞懂代碼的意思,注意SSE指令加載的數據低位在後,高位在前,所以我註釋裏也是這樣表達的。
考慮圖像數據的特性和通用性,是沒法使用16字節對齊的加載和保存的SIMD指令的,可是測試好像結果是這兩個指令對結果的影響差別已經很小的。
以上只是4*4大小的轉置,若是是圖像的轉置,則能夠和利用SSE優化圖像轉置一文提出的方式同樣,把圖像分紅不少個4*4的小塊,而後每一個小塊調用上述模塊。
考慮32位的特殊性,若是用純C語言實現轉置,可使用如下的代碼:
for (int Y = 0; Y < DstH; Y++) { int *LinePS = (int *)(Src + Y * 4); int *LinePD = (int *)(Dest + Y * StrideD); for (int X = 0; X < DstW; X++) { LinePD[X] = LinePS[0]; LinePS += DstH; } }
使用上述SSE的方式,則以下所示:
int BlockH = DstW / 4, BlockV = DstH / 4; for (int Y = 0; Y < BlockV * 4; Y += 4) { unsigned char *LinePS = Src + Y * 4; unsigned char *LinePD = Dest + Y * StrideD; for (int X = 0; X < BlockH * 4; X += 4) { Transpose4x4_I((int *)(LinePS + X * StrideS), (int *)(LinePD + X * 4), SrcW, DstW); } }
// 處理未被SSE覆蓋到的行和列方向的數據。
上述處理未被SSE覆蓋到的行和列方向的數據可由讀者自行完成,這部分的耗時能夠不計。
(2)、 灰度模式的SSE實現
爲何先提灰度,而不是24位是由於24位圖像使用SSE處理始終是個坑,而且是個很難填的坑,咱們把它放在最後。
有了上面的32位的轉置,對灰度模式的轉置基本思路也是定位在各類pack和unpack的組合了,由於SSE支持一次性讀取16個字節的數據,因此最原始的想法也是寫個16*16小塊的灰度轉置函數,可是因爲灰度數據一個像素就是一個字節,這種轉置的組合須要大量的SSE函數才能實現,並且因爲中間須要多個變量保存臨時結果,很難保證XMM寄存器的充分利用,經過一段時間的摸索和實踐,我認爲這不是最佳答案。
最終,我將解決方案鎖定在8*8大小塊的灰度轉置優化中,由於有_mm_loadl_epi64和_mm_storel_epi64兩個SSE函數能夠只加載和保存__m128i數據的低8位,能夠很好的解決保存和加載問題。加上其餘一些組合函數,完美的解決的灰度問題,核心代碼以下所示:
// 灰度數據的8*8轉置 void Transpose8x8_8U(unsigned char *Src, unsigned char *Dest, int StrideS, int StrideD) { __m128i S0 = _mm_loadl_epi64((__m128i *)(Src + 0 * StrideS)); // 0 0 0 0 0 0 0 0 A7 A6 A5 A4 A3 A2 A1 A0 __m128i S1 = _mm_loadl_epi64((__m128i *)(Src + 1 * StrideS)); // 0 0 0 0 0 0 0 0 B7 B6 B5 B4 B3 B2 B1 B0 __m128i S2 = _mm_loadl_epi64((__m128i *)(Src + 2 * StrideS)); // 0 0 0 0 0 0 0 0 C7 C6 C5 C4 C3 C2 C1 C0 __m128i S3 = _mm_loadl_epi64((__m128i *)(Src + 3 * StrideS)); // 0 0 0 0 0 0 0 0 D7 D6 D5 D4 D3 D2 D1 D0 __m128i S01 = _mm_unpacklo_epi8(S0, S1); // B7 A7 B6 A6 B5 A5 B4 A4 B3 A3 B2 A2 B1 A1 B0 A0 __m128i S23 = _mm_unpacklo_epi8(S2, S3); // D7 C7 D6 C6 D5 C5 D4 C4 D3 C3 D2 C2 D1 C1 D0 C0 __m128i S0123L = _mm_unpacklo_epi16(S01, S23); // D3 C3 B3 A3 D2 C2 B2 A2 D1 C1 B1 A1 D0 C0 B0 A0 __m128i S0123H = _mm_unpackhi_epi16(S01, S23); // D7 C7 B7 A7 D6 C6 B6 A6 D5 C5 B5 A5 D4 C4 B4 A4 __m128i S4 = _mm_loadl_epi64((__m128i *)(Src + 4 * StrideS)); // 0 0 0 0 0 0 0 0 A7 A6 A5 A4 A3 A2 A1 A0 __m128i S5 = _mm_loadl_epi64((__m128i *)(Src + 5 * StrideS)); // 0 0 0 0 0 0 0 0 B7 B6 B5 B4 B3 B2 B1 B0 __m128i S6 = _mm_loadl_epi64((__m128i *)(Src + 6 * StrideS)); // 0 0 0 0 0 0 0 0 C7 C6 C5 C4 C3 C2 C1 C0 __m128i S7 = _mm_loadl_epi64((__m128i *)(Src + 7 * StrideS)); // 0 0 0 0 0 0 0 0 D7 D6 D5 D4 D3 D2 D1 D0 __m128i S45 = _mm_unpacklo_epi8(S4, S5); // B7 A7 B6 A6 B5 A5 B4 A4 B3 A3 B2 A2 B1 A1 B0 A0 __m128i S67 = _mm_unpacklo_epi8(S6, S7); // D7 C7 D6 C6 D5 C5 D4 C4 D3 C3 D2 C2 D1 C1 D0 C0 __m128i S4567L = _mm_unpacklo_epi16(S45, S67); // H3 G3 F3 E3 H2 G2 F2 E2 H1 G1 F1 E1 H0 G0 F0 E0 __m128i S4567H = _mm_unpackhi_epi16(S45, S67); // H7 G7 F7 E7 H6 G6 F6 E6 H5 G5 F5 E5 H4 G4 F4 E4 __m128i T0 = _mm_unpacklo_epi32(S0123L, S4567L); // H1 G1 F1 E1 D1 C1 B1 A1 H0 G0 F0 E0 D0 C0 B0 A0 _mm_storel_epi64((__m128i *)(Dest + 0 * StrideD), T0); // H0 G0 F0 E0 D0 C0 B0 A0 _mm_storel_epi64((__m128i *)(Dest + 1 * StrideD), _mm_srli_si128(T0, 8)); // H1 G1 F1 E1 D1 C1 B1 A1 __m128i T1 = _mm_unpackhi_epi32(S0123L, S4567L); // H3 G3 F3 E3 D3 C3 B3 A3 H2 G2 F2 E2 D2 C2 B2 A2 _mm_storel_epi64((__m128i *)(Dest + 2 * StrideD), T1); _mm_storel_epi64((__m128i *)(Dest + 3 * StrideD), _mm_srli_si128(T1, 8)); __m128i T2 = _mm_unpacklo_epi32(S0123H, S4567H); // H5 G5 F5 E5 D5 C5 B5 H4 G4 F4 E4 A5 D4 C4 B4 A4 _mm_storel_epi64((__m128i *)(Dest + 4 * StrideD), T2); _mm_storel_epi64((__m128i *)(Dest + 5 * StrideD), _mm_srli_si128(T2, 8)); __m128i T3 = _mm_unpackhi_epi32(S0123H, S4567H); _mm_storel_epi64((__m128i *)(Dest + 6 * StrideD), T3); _mm_storel_epi64((__m128i *)(Dest + 7 * StrideD), _mm_srli_si128(T3, 8)); }
上述代碼也進行了詳細的註釋,標記了每一步數據是如何變化的,代碼充分利用了8位、16位、32位的pack組合,相信有SSE基礎的人都能看的懂,有的時候本身看着這段代碼都以爲是一種享受。
有幾個問題也在這裏留給你們,一個是保存__m128i數據的高8位有沒有不須要上述移位的方式而更高效的實現方式呢,第二就是咱們不必定拘泥於正方形的轉置,若是使用16*8的轉置效率會不會有變化或者說提高呢。
(3)、 BGR24位的SSE實現
24位咱們在PC上最常遇到的格式(手機上卻是基本不用),是最難以用SSE處理的,一個像素3個字節是的以4爲基本需求的一些SIMD函數難以發揮勇武之地,除了一些和像素成分無關的算法(也就是每一個通道都用相同的算法處理,而且算法和領域無關)外,都很難直接用SIMD處理,不少狀況下必須作一些轉換處理來提升適配性。
對於轉置,因爲一個像素佔用3個字節,若是徹底按照轉置的嚴格意義對24位圖像使用各類unpack來獲得結果,不是說作不到,可是將變得異常複雜,耗時耗力,而且不必定有加速做用,我這裏提出的方案是借用32位的來處理。
咱們也是一次性進行4*4的圖像塊的轉置,首先仍是讀取16個字節的信息,這裏就包括了5個多的24位像素的像素值,咱們只取前4個,並將它們擴展爲4個BGRA的格式,A值填充任何數據均可,而後使用32位的轉置算法,轉置獲得32位的結果,在將結果轉換到4個24位像素信息,因爲這中間只是借用了XMM寄存器或者一級或者二級緩存做爲保存數據的地址,沒有用到普通的中間內存,所以效率也很是之高。
部分代碼以下:
// BGR數據的轉置,藉助了BGRA的中間數據 void Transpose4x4_BGR(unsigned char *Src, unsigned char *Dest, int StrideS, int StrideD) { __m128i MaskBGR2BGRA = _mm_setr_epi8(/* */); // Mask爲-1的地方會自動設置數據爲0 __m128i MaskBGRA2BGR = _mm_setr_epi8(/* */); __m128i S0 = _mm_shuffle_epi8(_mm_loadu_si128((__m128i *)(Src + 0 * StrideS)), MaskBGR2BGRA); __m128i S1 = _mm_shuffle_epi8(_mm_loadu_si128((__m128i *)(Src + 1 * StrideS)), MaskBGR2BGRA); __m128i S01L = _mm_unpacklo_epi32(S0, S1); __m128i S01H = _mm_unpackhi_epi32(S0, S1); __m128i S2 = _mm_shuffle_epi8(_mm_loadu_si128((__m128i *)(Src + 2 * StrideS)), MaskBGR2BGRA); __m128i S3 = _mm_shuffle_epi8(_mm_loadu_si128((__m128i *)(Src + 3 * StrideS)), MaskBGR2BGRA); __m128i S23L = _mm_unpacklo_epi32(S2, S3); __m128i S23H = _mm_unpackhi_epi32(S2, S3); _mm_storeu_si128((__m128i *)(Dest + 0 * StrideD), _mm_shuffle_epi8(_mm_unpacklo_epi64(S01L, S23L), MaskBGRA2BGR)); _mm_storeu_si128((__m128i *)(Dest + 1 * StrideD), _mm_shuffle_epi8(_mm_unpackhi_epi64(S01L, S23L), MaskBGRA2BGR)); _mm_storeu_si128((__m128i *)(Dest + 2 * StrideD), _mm_shuffle_epi8(_mm_unpacklo_epi64(S01H, S23H), MaskBGRA2BGR)); _mm_storeu_si128((__m128i *)(Dest + 3 * StrideD), _mm_shuffle_epi8(_mm_unpackhi_epi64(S01H, S23H), MaskBGRA2BGR)); }
上述代碼中的部分數據被我用/* */給代替了,主要是我不想讓懶人直接使用,能作事的人這個數據確定能搞得定的。
因爲_mm_loadu_si128會一次性加載16個字節的數據,而咱們實際只使用了其前面的12個字節的信息,因此須要考慮程序的嚴謹性,對最後一行圖像分塊時應該注意不要超出圖像能訪問的數據範圍(我想不少人不會明白我這句話的意思的)。
(4)、 循環方式的影響
轉置操做會改變長寬的尺寸,可是必然有DstH = SrcW, DstW = SrcH, 最後的循環也有兩種方式,即按照原圖先行後列,或者按照目的圖先行後列,前一種方式訪問原圖的數據是連續的,可是寫入目的圖的時候是非連續的,後者訪問原圖的數據是非連續的,可是寫入目的圖是的地址是連續的,不管如何寫,都會有一個方向存在較大的Cache miss的可能性,這也是轉置難以提升速度的難點所在,可是通過測試,第二種方式彷佛速度來的仍是快一些,咱們以灰度圖爲例,前一種方式的寫法爲:
for (int Y = 0; Y < SrcH; Y++) { unsigned char *LinePS = Src + Y * StrideS; unsigned char *LinePD = Dest + Y; for (int X = 0; X < SrcW; X++) { LinePD[0] = LinePS[X]; LinePD += StrideD; } }
後一種爲:
for (int Y = 0; Y < DstH; Y++) { unsigned char *LinePS = Src + Y; unsigned char *LinePD = Dest + Y * StrideD; for (int X = 0; X < DstW; X++) { LinePD[X] = LinePS[0]; LinePS += StrideS; } }
在大部分的狀況下,後一種寫法會快不少,對於SSE的優化也是同樣的(因爲轉置的特性,上述兩種方式的SSE對應代碼的塊代碼是同樣的),這一塊的緣由雖然有些想法,但恐怕理解的不到位,這裏也就不闡述了,望有經驗的老司機指點。
(5)、 時間比較
100次重複轉置耗時(單位ms)
圖像大小(W*H)) |
1024*768 |
3000*2000 |
4000*3000 |
|
灰度模式 |
普通C語言 |
92 |
1398 |
7278 |
SSE |
18 |
294 |
1015 |
|
BGR 24位 |
普通C語言 |
145 |
4335 |
9239 |
SSE |
43 |
1067 |
2270 |
|
BGRA 32位 |
普通C語言 |
78 |
4338 |
9797 |
SSE |
51 |
1214 |
2690 |
可見SSE優化後相比普通的C語言仍是至關可觀的,特別是灰度模式的,對於大圖能夠達到6倍左右的提速。
同時由上表也能夠看出,圖像越大,彷佛提速比越大,我分析認爲是當圖像較小時,訪問相鄰行時的Cache miss的可能性要比大圖時爲小,所以SSE優化的提速不是特明顯,而大圖時Cache miss的機率會增長,這個時候SSE一次性處理多個像素的優勢就能充分體現了。
注:做者注意到在部分PC上測試時,SSE的加速沒有如此明顯,特別是對於小圖。
在 CUDA學習筆記一:CUDA+OpenCV的圖像轉置,採用Shared Memory進行CUDA程序優化 一文中提供的Lena灰度測試圖片大小爲512*512的,使用上述算法執行100次只須要6ms,原文提供的時間GPU使用都須要0.7ms,雖然不清楚他的GPU是啥配置的,可是可見本例的優化速度至關可觀的。
總的來講,轉置操做的大部分耗時都是在訪問內存上,這是個很大的瓶頸,使用CPU能優化的空間也是有限的,可是隻要能優化,就應該充分榨取CPU的資源。
核心代碼都已經共享了,由須要的朋友請自行整理成工程。
有興趣的朋友也能夠試試AVX的優化速度。
比較工程: http://files.cnblogs.com/files/Imageshop/Transpose.rar