c++ - 使用SSE最快缩小8位灰度图像

Question

我有一个将8位图像缩小两倍的功能。我以前有optimised the rgb32 case with SSE。现在，我想对gray8案例做同样的事情。

在核心处，有一个函数接受两行像素数据，其工作方式如下:

/** 
 * Calculates the average of two rows of gray8 pixels by averaging four pixels.
 */
void average2Rows(const uint8_t* row1, const uint8_t* row2, uint8_t* dst, int size)
{
    for (int i = 0; i < size - 1; i += 2)
        *(dst++) = ((row1[i]+row1[i+1]+row2[i]+row2[i+1])/4)&0xFF;
}

现在，我想出了一种SSE变体，它的速度快了大约三倍，但是它确实涉及很多改组，我认为这样做可能会更好。有人在这里看到可以优化的地方吗？

/* row1: 16 8-bit values A-P
 * row2: 16 8-bit values a-p
 * returns 16 8-bit values (A+B+a+b)/4, (C+D+c+d)/4, ..., (O+P+o+p)/4
 */
__m128i avg16Bytes(const __m128i& row1, const __m128i& row2)
{
    static const __m128i  zero = _mm_setzero_si128(); 

    __m128i ABCDEFGHIJKLMNOP = _mm_avg_epu8(row1_u8, row2);

    __m128i ABCDEFGH  = _mm_unpacklo_epi8(ABCDEFGHIJKLMNOP, zero);
    __m128i IJKLMNOP  = _mm_unpackhi_epi8(ABCDEFGHIJKLMNOP, zero);

    __m128i AIBJCKDL = _mm_unpacklo_epi16( ABCDEFGH, IJKLMNOP );
    __m128i EMFNGOHP = _mm_unpackhi_epi16( ABCDEFGH, IJKLMNOP );

    __m128i AEIMBFJN = _mm_unpacklo_epi16( AIBJCKDL, EMFNGOHP );
    __m128i CGKODHLP = _mm_unpackhi_epi16( AIBJCKDL, EMFNGOHP );

    __m128i ACEGIKMO = _mm_unpacklo_epi16( AEIMBFJN, CGKODHLP );
    __m128i BDFHJLNP = _mm_unpackhi_epi16( AEIMBFJN, CGKODHLP );

    return _mm_avg_epu8(ACEGIKMO, BDFHJLNP);
}

/*
 * Calculates the average of two rows of gray8 pixels by averaging four pixels.
 */
void average2Rows(const uint8_t* src1, const uint8_t* src2, uint8_t* dst, int size)
{
    for(int i = 0;i<size-31; i+=32)
    {
        __m128i tl = _mm_loadu_si128((__m128i const*)(src1+i));
        __m128i tr = _mm_loadu_si128((__m128i const*)(src1+i+16));
        __m128i bl = _mm_loadu_si128((__m128i const*)(src2+i));
        __m128i br = _mm_loadu_si128((__m128i const*)(src2+i+16)))

        __m128i l_avg = avg16Bytes(tl, bl);
        __m128i r_avg = avg16Bytes(tr, br);

        _mm_storeu_si128((__m128i *)(dst+(i/2)), _mm_packus_epi16(l_avg, r_avg));
    }
}

笔记:

我意识到我的函数有轻微的(四舍五入)舍入错误，但是我愿意接受这一点。

为了清楚起见，我假设大小是32的倍数。

编辑:现在有一个 github repository实现此问题的答案。用户 Peter Cordes提供了最快的解决方案。有关详细信息，请参见下面的他的文章:

__m128i avg16Bytes(const __m128i& row1, const __m128i& row2)
{
    // Average the first 16 values of src1 and src2:
    __m128i avg = _mm_avg_epu8(row1, row2);

    // Unpack and horizontal add:
    avg = _mm_maddubs_epi16(avg, _mm_set1_epi8(1));

    // Divide by 2:
    return  _mm_srli_epi16(avg, 1);
}

通过计算 (a+b)/2 + (c+d)/2而不是 (a+b+c+d)/4，它可以作为我的原始实现，因此它具有相同的一乘四舍五入误差。

向用户 Paul R表示祝贺，以实现比我的速度快两倍的解决方案，但准确的是:

__m128i avg16Bytes(const __m128i& row1, const __m128i& row2)
{
    // Unpack and horizontal add:
    __m128i row1 = _mm_maddubs_epi16(row1_u8, _mm_set1_epi8(1));
    __m128i row2 = _mm_maddubs_epi16(row2_u8, _mm_set1_epi8(1));

    // vertical add:
    __m128i avg = _mm_add_epi16(row1_avg, row2_avg);              

    // divide by 4:
    return _mm_srli_epi16(avg, 2);                     
}

Answer 1

如果您愿意接受两次使用pavgb的双舍入处理，则可以通过首先使用pavgb进行垂直平均来比Paul R的答案更快，将需要解压缩为16位元素的数据量减少了一半。 (并将一半的负载折叠到pavgb的内存操作数中，从而减少了某些CPU的前端瓶颈。)

对于水平平均，最好的选择可能仍然是pmaddubsw和set1(1)并移位1，然后打包。

// SSSE3 version
// I used `__restrict__` to give the compiler more flexibility in unrolling
void average2Rows_doubleround(const uint8_t* __restrict__ src1, const uint8_t*__restrict__ src2,
                              uint8_t*__restrict__ dst, size_t size)
{
    const __m128i vk1 = _mm_set1_epi8(1);
    size_t dstsize = size/2;
    for (size_t i = 0; i < dstsize - 15; i += 16)
    {
        __m128i v0 = _mm_load_si128((const __m128i *)&src1[i*2]);
        __m128i v1 = _mm_load_si128((const __m128i *)&src1[i*2 + 16]);
        __m128i v2 = _mm_load_si128((const __m128i *)&src2[i*2]);
        __m128i v3 = _mm_load_si128((const __m128i *)&src2[i*2 + 16]);
        __m128i left  = _mm_avg_epu8(v0, v2);
        __m128i right = _mm_avg_epu8(v1, v3);

        __m128i w0 = _mm_maddubs_epi16(left, vk1);        // unpack and horizontal add
        __m128i w1 = _mm_maddubs_epi16(right, vk1);
        w0 = _mm_srli_epi16(w0, 1);                     // divide by 2
        w1 = _mm_srli_epi16(w1, 1);
        w0 = _mm_packus_epi16(w0, w1);                  // pack

        _mm_storeu_si128((__m128i *)&dst[i], w0);
    }
}

另一个选项是 _mm_srli_epi16(v, 8)，用于将每个水平对中的奇数元素与偶数元素对齐。但是，由于没有水平截断的水平包装，因此在包装之前，您必须先对两半进行 _mm_and_si128(v, _mm_set1_epi16(0x00FF))编码。事实证明，这比使用SSSE3 pmaddubsw慢，尤其是在没有AVX的情况下，AVX需要额外的MOVDQA指令来复制寄存器。

void average2Rows_doubleround_SSE2(const uint8_t* __restrict__ src1, const uint8_t* __restrict__ src2, uint8_t* __restrict__ dst, size_t size)
{
    size /= 2;
    for (size_t i = 0; i < size - 15; i += 16)
    {
        __m128i v0 = _mm_load_si128((__m128i *)&src1[i*2]);
        __m128i v1 = _mm_load_si128((__m128i *)&src1[i*2 + 16]);
        __m128i v2 = _mm_load_si128((__m128i *)&src2[i*2]);
        __m128i v3 = _mm_load_si128((__m128i *)&src2[i*2 + 16]);
        __m128i left  = _mm_avg_epu8(v0, v2);
        __m128i right = _mm_avg_epu8(v1, v3);

        __m128i l_odd  = _mm_srli_epi16(left, 8);   // line up horizontal pairs
        __m128i r_odd  = _mm_srli_epi16(right, 8);

        __m128i l_avg = _mm_avg_epu8(left, l_odd);  // leaves garbage in the high halves
        __m128i r_avg = _mm_avg_epu8(right, r_odd);

        l_avg = _mm_and_si128(l_avg, _mm_set1_epi16(0x00FF));
        r_avg = _mm_and_si128(r_avg, _mm_set1_epi16(0x00FF));
        __m128i avg   = _mm_packus_epi16(l_avg, r_avg);          // pack
        _mm_storeu_si128((__m128i *)&dst[i], avg);
    }
}

对于AVX512BW，有 _mm_cvtepi16_epi8，但是IACA表示它在Skylake-AVX512上是2 oups，并且只需要输入1并产生半角输出。根据IACA，内存目标形式是4个未融合域的总和(与reg，reg +单独存储相同)。我必须使用 _mm_mask_cvtepi16_storeu_epi8(&dst\[i+0\], -1, l_avg); 来获取它，因为gcc和clang无法将单独的 _mm_store折叠到 vpmovwb的存储目标中。 (没有非屏蔽存储内在函数，因为编译器应该像为典型ALU指令将 _mm_load折叠到内存操作数中一样为您做到这一点)。

仅在缩小到1/4或1/8( cvtepi64_epi8)而不是一半时，它才有用。也许对避免再次洗牌以处理 _mm512_packus_epi16的车内行为很有用。使用AVX2，在 _mm256_packus_epi16上的 [D C] [B A]之后，您拥有 [D B | C A]，可以使用AVX2 _mm256_permute4x64_epi64 (__m256i a, const int imm8)对其进行修复，以按64位块进行随机播放。但是对于AVX512，您需要对 vpermq进行 vector 随机播放控制。 packus +修正改组可能仍然是一个更好的选择。

完成此操作后，循环中将没有太多 vector 指令，而 可以让编译器使asm 更加紧密，从而有很多收获。不幸的是，对于编译器来说，完成循环很困难。
(这也有助于Paul R的解决方案，因为他从问题中复制了对编译器不利的循环结构。)

使用循环计数器的方式应使gcc / clang可以更好地优化，并使用避免每次循环都重新进行符号扩展的类型。

在当前循环中，gcc / clang实际上对 i/2进行算术右移，而不是增加16(而不是32)并使用缩放索引寻址模式进行加载。看来他们没有意识到 i总是偶数。

(full code + asm on Matt Godbolt's compiler explorer):

.LBB1_2:     ## clang's inner loop for int i, dst[i/2] version
    movdqu  xmm1, xmmword ptr [rdi + rcx]
    movdqu  xmm2, xmmword ptr [rdi + rcx + 16]
    movdqu  xmm3, xmmword ptr [rsi + rcx]
    movdqu  xmm4, xmmword ptr [rsi + rcx + 16]
    pavgb   xmm3, xmm1
    pavgb   xmm4, xmm2
    pmaddubsw       xmm3, xmm0
    pmaddubsw       xmm4, xmm0
    psrlw   xmm3, 1
    psrlw   xmm4, 1
    packuswb        xmm3, xmm4

    mov     eax, ecx         # This whole block is wasted instructions!!!
    shr     eax, 31
    add     eax, ecx
    sar     eax              # eax = ecx/2, with correct rounding even for negative `i`
    cdqe                     # sign-extend EAX into RAX

    movdqu  xmmword ptr [rdx + rax], xmm3
    add     rcx, 32          # i += 32
    cmp     rcx, r8
    jl      .LBB1_2          # }while(i < size-31)

gcc7.1并不是很糟糕(只是 mov / sar / movsx)，但是gcc5.x和6.x分别为src1和src2以及商店的计数器/索引单独增加了指针。 (完全脑残行为，尤其是因为他们仍然使用 -march=sandybridge这样做。索引的 movdqu存储和未索引的 movdqu加载为您提供了最大的循环开销。)

无论如何，在循环中使用 dstsize并乘以 i而不是将其除以给出更好的结果。 gcc和clang的不同版本将其可靠地编译为单个循环计数器，它们与缩放索引寻址模式一起用于负载。您得到的代码如下:

    movdqa  xmm1, xmmword ptr [rdi + 2*rax]
    movdqa  xmm2, xmmword ptr [rdi + 2*rax + 16]
    pavgb   xmm1, xmmword ptr [rsi + 2*rax]
    pavgb   xmm2, xmmword ptr [rsi + 2*rax + 16]   # saving instructions with aligned loads, see below
    ...
    movdqu  xmmword ptr [rdx + rax], xmm1
    add     rax, 16
    cmp     rax, rcx
    jb      .LBB0_2

我使用 size_t i来匹配 size_t的大小，以确保gcc不会浪费任何指令将其符号扩展或零扩展到指针的宽度。 (不过，零扩展通常是免费的，因此 unsigned size和 unsigned i可能还行，并保存了几个REX前缀。)

您仍然可以摆脱 cmp，但将索引的计数提高到0，这比我所做的要快得多。我不确定让编译器变得不愚蠢并且如果您算上零就忽略 cmp指令会多么容易。不过，从对象末尾开始索引是没有问题的。 src1+=size;。但是，如果您要使用unaligned-cleanup循环，它会使事情复杂化。

在我的Skylake i7-6700k上(最大睿频速度为4.4GHz，但请查看时钟周期计数而不是时间)。使用g++ 7.1时，对于1024字节的100M次重复，此间隔为2.7秒，而间隔为3.3秒。

 Performance counter stats for './grayscale-dowscale-by-2.inline.gcc-skylake-noavx' (2 runs):

   2731.607950      task-clock (msec)         #    1.000 CPUs utilized            ( +-  0.40% )
             2      context-switches          #    0.001 K/sec                    ( +- 20.00% )
             0      cpu-migrations            #    0.000 K/sec                  
            88      page-faults:u             #    0.032 K/sec                    ( +-  0.57% )
11,917,723,707      cycles                    #    4.363 GHz                      ( +-  0.07% )
42,006,654,015      instructions              #    3.52  insn per cycle           ( +-  0.00% )
41,908,837,143      uops_issued_any           # 15342.186 M/sec                   ( +-  0.00% )
49,409,631,052      uops_executed_thread      # 18088.112 M/sec                   ( +-  0.00% )
 3,301,193,901      branches                  # 1208.517 M/sec                    ( +-  0.00% )
   100,013,629      branch-misses             #    3.03% of all branches          ( +-  0.01% )

   2.731715466 seconds time elapsed                                          ( +-  0.40% )

与相同的向量化，但是 与int i和dst[i/2]一起创建了更高的循环开销(更多标量指令):

 Performance counter stats for './grayscale-dowscale-by-2.loopoverhead-aligned-inline.gcc-skylake-noavx' (2 runs):

   3314.335833      task-clock (msec)         #    1.000 CPUs utilized            ( +-  0.02% )
             4      context-switches          #    0.001 K/sec                    ( +- 14.29% )
             0      cpu-migrations            #    0.000 K/sec                  
            88      page-faults:u             #    0.026 K/sec                    ( +-  0.57% )
14,531,925,552      cycles                    #    4.385 GHz                      ( +-  0.06% )
51,607,478,414      instructions              #    3.55  insn per cycle           ( +-  0.00% )
51,109,303,460      uops_issued_any           # 15420.677 M/sec                   ( +-  0.00% )
55,810,234,508      uops_executed_thread      # 16839.040 M/sec                   ( +-  0.00% )
 3,301,344,602      branches                  #  996.080 M/sec                    ( +-  0.00% )
   100,025,451      branch-misses             #    3.03% of all branches          ( +-  0.00% )

   3.314418952 seconds time elapsed                                          ( +-  0.02% )

与Paul R的版本(针对较低的循环开销进行了优化):准确但较慢

Performance counter stats for './grayscale-dowscale-by-2.paulr-inline.gcc-skylake-noavx' (2 runs):

   3751.990587      task-clock (msec)         #    1.000 CPUs utilized            ( +-  0.03% )
             3      context-switches          #    0.001 K/sec                  
             0      cpu-migrations            #    0.000 K/sec                  
            88      page-faults:u             #    0.024 K/sec                    ( +-  0.56% )
16,323,525,446      cycles                    #    4.351 GHz                      ( +-  0.04% )
58,008,101,634      instructions              #    3.55  insn per cycle           ( +-  0.00% )
57,610,721,806      uops_issued_any           # 15354.709 M/sec                   ( +-  0.00% )
55,505,321,456      uops_executed_thread      # 14793.566 M/sec                   ( +-  0.00% )
 3,301,456,435      branches                  #  879.921 M/sec                    ( +-  0.00% )
   100,001,954      branch-misses             #    3.03% of all branches          ( +-  0.02% )

   3.752086635 seconds time elapsed                                          ( +-  0.03% )

与 对比Paul R的原始版本带有额外的循环开销:

Performance counter stats for './grayscale-dowscale-by-2.loopoverhead-paulr-inline.gcc-skylake-noavx' (2 runs):

   4154.300887      task-clock (msec)         #    1.000 CPUs utilized            ( +-  0.01% )
             3      context-switches          #    0.001 K/sec                  
             0      cpu-migrations            #    0.000 K/sec                  
            90      page-faults:u             #    0.022 K/sec                    ( +-  1.68% )
18,174,791,383      cycles                    #    4.375 GHz                      ( +-  0.03% )
67,608,724,157      instructions              #    3.72  insn per cycle           ( +-  0.00% )
66,937,292,129      uops_issued_any           # 16112.769 M/sec                   ( +-  0.00% )
61,875,610,759      uops_executed_thread      # 14894.350 M/sec                   ( +-  0.00% )
 3,301,571,922      branches                  #  794.736 M/sec                    ( +-  0.00% )
   100,029,270      branch-misses             #    3.03% of all branches          ( +-  0.00% )

   4.154441330 seconds time elapsed                                          ( +-  0.01% )

请注意，分支未命中次数与重复次数大致相同:每次循环结束时，内部循环都会发生错误的预测。展开以使循环迭代次数保持在约22以下将使该模式足够短，以使Skylake的分支预测器在大多数时间正确预测未采取的情况。分支错误预测是我们每个周期不能在管道中获得约4.0 ups的唯一原因，因此避免分支未命中会将IPC从3.5提高到超过4.0(cmp / jcc宏融合将2条指令放入一个uop中)。

即使您在二级缓存带宽(而不是前端)上遇到瓶颈，这些分支丢失也可能会受到伤害。不过，我没有进行测试:我的测试只是在Paul R的测试工具的函数调用周围包裹了一个 for()循环，因此L1D缓存中的所有内容都很热门。内部循环的32次迭代在这里接近最坏的情况:对于频繁的错误预测而言足够低，但又不那么低，以至于分支预测可以选择并避免这种模式。

我的版本每次迭代应运行3个周期，仅在前端，英特尔Sandybridge及更高版本上成为瓶颈。 (Nehalem每个时钟将在一个负载上成为瓶颈。)

请参阅 http://agner.org/optimize/以及 Can x86's MOV really be "free"? Why can't I reproduce this at all?，以获取有关融合域与未融合域uops和perf计数器的更多信息。

更新:clang至少在大小为编译时常量时为您解压缩 ...奇怪的是，它甚至解压缩了 dst[i/2]函数的非内联版本(具有未知的 size)，但未展开下部循环-开销版本。

使用 clang++-4.0 -O3 -march=skylake -mno-avx，我的版本(编译器将其展开为2)运行于:9.61G周期，适用于100M iters(2.2s)。 (发出了35.6G的uops(融合域)，执行了45.0G的uops(融合域)，接近零的分支未命中。)可能不再是前端的瓶颈，但是AVX仍然会受到伤害。

Paul R(也展开了2个)以12.29G的周期运行100M次迭代(2.8s)。已发出48.4G uops(融合域)，已执行51.4G uops(未融合域)。对于4.08 IPC，50.1G指令可能仍然在前端出现瓶颈(因为在销毁寄存器之前，它需要一对 movdqa指令来复制寄存器)。即使没有AVX2来实现更宽的整数 vector ，AVX也将有助于无损 vector 指令。

通过仔细的编码，您应该能够很好地解决运行时可变大小的问题。

使用对齐的指针和对齐的加载，因此编译器可以将pavgb与内存操作数一起使用，而不是使用单独的未对齐加载指令。这意味着更少的指令和更少的前端操作，这是此循环的瓶颈。

这对Paul的版本无济于事，因为 pmaddubsw的第二个操作数只能来自内存，这就是被视为有符号字节的那个。如果我们使用 _mm_maddubs_epi16(_mm_set1_epi8(1), v0);，则16位乘法结果将被符号扩展而不是零扩展。因此 1+255会显示为0，而不是256。

折叠负载需要与SSE对齐，而不需要与AVX对齐。但是， on Intel Haswell/Skylake, indexed addressing modes can only stay micro-fused with instructions which read-modify-write their destination register. vpavgb xmm0, xmm0, [rsi+rax*2]在发布到内核的乱序部分之前，在Haswell / Skylake上未分层为2 oups，但是 pavgb xmm1, [rsi+rax*2]可以一直保持微融合状态，因此它作为单个uop发出。前端问题瓶颈是除Ryzen(即非Atom / Silvermont)以外的主流x86 CPU上每个时钟4个融合域uops。将一半的负载折叠到内存操作数中，可以帮助除Sandybridge / Ivybridge之外的所有Intel CPU和所有AMD CPU减轻负载。

当内联到使用 alignas(32)的测试函数中时，即使您使用 _mm_loadu内部函数，gcc和clang也会折叠负载。他们知道数据是一致的，并且可以利用。

奇怪的事实:在启用了AVX代码源( -march=native)的情况下编译128b向量化的代码实际上会减慢Haswell / Skylake的速度，因为即使这四个负载都是 vpavgb的内存操作数，它也会使所有4个负载以单独的uops发出没有AVX会避免的任何 movdqa寄存器复制指令。 (由于3操作数指令不会破坏其输入之一，因此，即使对于仍仅使用128b vector 的手动矢量化代码，AVX通常还是会领先于此。)在这种情况下， 13,53G cycles ( +- 0.05% )或 3094.195773 ms ( +- 0.20% )是从 11.92G循环开始的约2.7秒。 uops_issued = 48.508G，高于 41,908。指令计数和uops_exected计数基本相同。

OTOH，实际的256b AVX2版本运行速度略低于两倍。采取一些措施来减少前端瓶颈肯定会有所帮助。根据@Mysticial的测试，AVX512版本在Skylake-AVX512 Xeon上的运行速度几乎是Skyx-AVX512 Xeon的4倍，但可能会限制ALU吞吐量，因为当RS中有512b联运处理器等待执行时，SKX会关闭执行端口1。 (这说明了 pavgb zmm has 1 per clock throughput while pavgb ymm is 2 per clock.的原因。)

要使两个输入行对齐， 将您的图像数据以行距为16的倍数的格式存储，即使实际图像尺寸是奇数。您的存储跨度不必与实际图像尺寸匹配。

如果您只能对齐源或目标(例如，因为您正在缩小从源图像中奇数列开始的区域的比例)，则可能仍应对齐源指针。

如果无法同时调整目标位置，则英特尔的优化手册建议将目标位置而不是源位置对齐，但是将负载加载为商店的4倍可能会改变余额。

要处理起点/终点的未对齐，请从起点和终点进行可能重叠的未对齐像素 vector 。商店与其他商店重叠是很好的，并且由于dst与src是分开的，因此您可以重做部分重叠的 vector 。

在Paul的测试 main()中，我只是在每个数组的前面添加了 alignas(32)。

AVX2:

由于使用了 compile one version with -march=native ，因此可以使用 #ifdef __AVX2__在编译时轻松检测AVX2。没有简单的方法对128b和256b手动矢量化使用完全相同的代码。所有内部函数都有不同的名称，因此即使没有其他差异，您通常也需要复制所有内容。

(有一些用于内部函数的C++包装库，这些库使用运算符重载和函数重载来让您编写在不同宽度的 vector 上使用相同逻辑的模板化版本。例如，Agner Fog的VCL是好的，但是除非您的软件是开放式的，源，您不能使用它，因为它已获得GPL许可，并且您想分发二进制文件。)

要利用二进制发行版中的AVX2，必须执行运行时检测/调度。在这种情况下，您希望分派(dispatch)到循环行的函数版本，因此不必在行循环内分派(dispatch)开销。或者只是让该版本使用SSSE3。