SSE图像算法优化系列四：图像转置的SSE优化（支持8位、24位、32位），提速4-6倍

一、前言

      转置操作在很多算法上都有着广泛的应用，在数学上矩阵转置更有着特殊的意义。而在图像处理上，如果说图像数据本身的转置，除了显示外，本身并无特殊含义，但是在某些情况下，确能有效的提高算法效率，比如很多行列可分离的算法，在很多情况下，行和列方向的算法逻辑随相同，但是由于多方面原因（比如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）   

     可见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的优化速度。

     比较工程：   https://files.cnblogs.com/files/Imageshop/Transpose.rar