SSE图像算法优化系列四:图像转置的SSE优化(支持8位、24位、32位),提速4-6倍

1、前言

      转置操做在不少算法上都有着普遍的应用,在数学上矩阵转置更有着特殊的意义。而在图像处理上,若是说图像数据自己的转置,除了显示外,自己并没有特殊含义,可是在某些状况下,确能有效的提升算法效率,好比不少行列可分离的算法,在不少状况下,行和列方向的算法逻辑随相同,可是因为多方面缘由(好比Cache miss, 优化水平等)行列处理时间仍是由很大的差别的,这个时候若是转置的耗时和处理时间相比所占比例甚小,则能够考虑在进行耗时处理前先转置数据,而后调用不耗时的方向的算法,处理完后再次进行转置。所以,一个高效的图像转置算法的设计时很是有必要的。算法

2、目前的情况

     没怎么搜集这方面的资料,不过在百度上看到的优化的帖子也有几篇:缓存

     一、利用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

 

相关文章
相关标签/搜索