一、前言
轉置操作在很多算法上都有着廣泛的應用，在數學上矩陣轉置更有着特殊的意義。而在圖像處理上，如果說圖像數據本身的轉置，除了顯示外，本身並無特殊含義，但是在某些情況下，確能有效的提高算法效率，比如很多行列可分離的算法，在很多情況下，行和列方向的算法邏輯隨相同，但是由於多方面原因（比如Cache miss， 優化水平等）行列處理時間還是由很大的差異的，這個時候如果轉置的耗時和處理時間相比所占比例甚小，則可以考慮在進行耗時處理前先轉置數據，然后調用不耗時的方向的算法，處理完后再次進行轉置。因此，一個高效的圖像轉置算法的設計時非常有必要的。
二、目前的狀況
沒怎么搜集這方面的資料，不過在百度上看到的優化的帖子也有幾篇：
1、利用SSE優化圖像轉置 這篇文章講到了SSE優化轉置操作，講的很簡單，我只是稍微看了下他的代碼，他似乎處理的不是普通的8位圖像，而是16位的，反正我是沒有看懂，並且他的提供比較的C代碼本身寫法就完全沒有考慮到C語言自身的優化，因此最后提出SSE代碼比C快5倍說服力就大為打折扣，不過他這里可以值得學習的地方就是這個轉置支持In-Place操作，就是Src和Dest可以相同。
2、圖像轉置的Neon優化代碼 Neon的代碼，沒看懂，不過后面說10倍左右的提速，其實也要看原始的C代碼是怎么寫的了，不過原文也明確的說，只支持RGBA 32位的圖，顯然作者也避而不談灰度或者24位的圖，當然這於手機端似乎沒有24位的概念有關。
3、CUDA學習筆記一：CUDA+OpenCV的圖像轉置，采用Shared Memory進行CUDA程序優化 這個文章是說GPU的優化，不過最后給出的GPU時間和CPU相比真的很慘。
也就是說國內網絡上的優化文章其實都還是停留在皮毛階段，也或者是真正的具有優化意義的代碼都還雪藏在某個公司或者某個人的硬盤里，特別是針對灰度和24位圖像的轉置優化在PC上有更多的使用場景。
三、我的貢獻
普通的C代碼的轉置很簡單，也曾嘗試過各種優化方案，但是最后都無啥特別大的改進，因此考慮使用SSE的方案。
（1）、由奢入儉難啊，我們先挑最簡單來實現，我說的最簡單的是32位圖像。
32位圖像由B/G/R/A 4個分量組成，轉置時我們需要把他們看成一個整體，以4*4大小的轉置威力，如下所示：
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
