平方差和的 SSE 优化答案

【问题标题】：SSE optimization of sum of squared differences平方差和的 SSE 优化
【发布时间】：2017-07-31 18:51:00
【问题描述】：

我最近发现我的程序大部分时间都花在了以下简单的函数上：

void SumOfSquaredDifference(
    const uint8_t * a, size_t aStride, const uint8_t * b, size_t bStride, 
    size_t width, size_t height, uint64_t * sum)
{
    *sum = 0;
    for(size_t row = 0; row < height; ++row)
    {
        int rowSum = 0;
        for(size_t col = 0; col < width; ++col)
        {
            int d = a[col] - b[col];
            rowSum += d*d;
        }
        *sum += rowSum;
        a += aStride;
        b += bStride;
    }
}

此函数求两个 8 位灰度图像的平方差之和。我认为使用 SSE 可以提高其性能，但我在这方面没有经验。有人可以帮帮我吗？

【问题讨论】：

在开始手动优化之前，您需要检查您的编译器是否已经生成矢量化代码。确保您已启用优化和 SIMD，并检查生成的 SSE 指令代码，否则您可能会浪费时间做傻事。（顺便说一句，请说明您使用的 CPU 硬件、操作系统和编译器。）
你的外循环似乎没用。这是一个错误吗？
@A.S.H：我知道这看起来很奇怪，但他在每次行迭代时都会碰到 a、b 指针。
Clang 和 GCC 自动矢量化这个，但不是很好（使用 vpmulld，这是低吞吐量）
我希望编译器对使用 SSE 犹豫不决。 height、width、aStride 和 bStride 都是运行时参数。这意味着对于 SSE（16 字节），每行的开头和结尾都不能适当对齐

标签： c++ image-processing optimization sse simd

【解决方案1】：

当然，您可以改进您的代码。这是使用 SSE2 优化函数的示例：

const __m128i Z = _mm_setzero_si128();
const size_t A = sizeof(__m128i);

inline __m128i SquaredDifference(__m128i a, __m128i b)
{
    const __m128i aLo = _mm_unpacklo_epi8(a, Z);
    const __m128i bLo = _mm_unpacklo_epi8(b, Z);
    const __m128i dLo = _mm_sub_epi16(aLo, bLo);

    const __m128i aHi = _mm_unpackhi_epi8(a, Z);
    const __m128i bHi = _mm_unpackhi_epi8(b, Z);
    const __m128i dHi = _mm_sub_epi16(aHi, bHi);

    return _mm_add_epi32(_mm_madd_epi16(dLo, dLo), _mm_madd_epi16(dHi, dHi));
}

inline __m128i HorizontalSum32(__m128i a)
{
    return _mm_add_epi64(_mm_unpacklo_epi32(a, Z), _mm_unpackhi_epi32(a, Z));
}

inline uint64_t ExtractSum64(__m128i a)
{
    uint64_t  _a[2];
    _mm_storeu_si128((__m128i*)_a, a);
    return _a[0] + _a[1];
}

void SumOfSquaredDifference(
    const uint8_t *a, size_t aStride, const uint8_t *b, size_t bStride, 
    size_t width, size_t height, uint64_t * sum)
{
    assert(width%A == 0 && width < 0x10000);
    __m128i fullSum = Z;
    for(size_t row = 0; row < height; ++row)
    {
        __m128i rowSum = Z;
        for(size_t col = 0; col < width; col += A)
        {
            const __m128i a_ = _mm_loadu_si128((__m128i*)(a + col));
            const __m128i b_ = _mm_loadu_si128((__m128i*)(b + col)); 
            rowSum = _mm_add_epi32(rowSum, SquaredDifference(a_, b_));
        }
        fullSum = _mm_add_epi64(fullSum, HorizontalSum32(rowSum));
        a += aStride;
        b += bStride;
    }
    *sum = ExtractSum64(fullSum);
}

这个例子是一些简化的（如果图像宽度不是 16 的倍数，它就不起作用）。算法的完整版是here。

还有来自SSSE3 版本的一些魔法：

const __m128i K_1FF = _mm_set1_epi16(0x1FF);

inline __m128i SquaredDifference(__m128i a, __m128i b)
{
    const __m128i lo = _mm_maddubs_epi16(_mm_unpacklo_epi8(a, b), K_1FF);
    const __m128i hi = _mm_maddubs_epi16(_mm_unpackhi_epi8(a, b), K_1FF);
    return _mm_add_epi32(_mm_madd_epi16(lo, lo), _mm_madd_epi16(hi, hi));
}

魔法描述（见_mm_maddubs_epi16）：

K_1FF -> {-1, 1, -1, 1, ...};
_mm_unpacklo_epi8(a, b) -> {a0, b0, a1, b1, ...};
_mm_maddubs_epi16(_mm_unpacklo_epi8(a, b), K_1FF) -> {b0 - a0, b1 - a1, ...};

【讨论】：

SSSE3 版本发生了什么？
感谢图书馆链接！太棒了！
它和我想要的很接近。

【解决方案2】：

GCC 有一些开关可以鼓励它对代码进行矢量化处理。例如，-mfma 开关让我在像这样的简单循环中使用双打提高了大约 25% 的速度。我想使用 8 位整数会更好。我更喜欢手写优化，因为您的代码保持可读性。

也就是说，有一些老技巧可以加快你的循环：

不要索引，在每次循环迭代中递增指针。您在外循环中执行此操作，您应该在内循环中执行相同操作。您可以在进入内部循环之前创建一个新指针，因此+=stride 保持有效。
不要分配给循环内的 sum 指针，使用局部变量进行累加并在完成后复制到输出。您使用rowSum，但仅在内部循环中。而是在两个循环中使用该变量。

【讨论】：