image-processing - SIMD 像素对比度 : sum of differences between a pixel and its neighbors (uint16_t color components, 浮点总和)？

Question

使用 SIMD/汇编器减去 2 个 uint16 的绝对值(最大差值)并将结果添加 (+=) 到 float 的最佳方法是什么？

类似于这个 C'ish 示例

c0 += fabs((float)a0 - (float)b0);  // C is Float accumulator, a+b pixels

其中 a 和 b 是无符号 16 位字，c 是 float 。只有 1 个字 -> 浮点转换而不是 3 个。

该应用程序正在处理尽可能多的完整 RGB 像素上的原始 16 位无符号整型图像数据。

可能在 Skylake Xeon E3-1275 v5 上使用 AVX2/SSE4.2？

---

5分钟评论限制？无法保存或重新编辑？？？

Are you sure you need float? Uint16 can't accumulate more than 1 subtraction. I want to do a neighborhood contrast calc so I need to sum at least 8 differences. There are (2D+1)^2-1 neighbors in a neighborhood with D depth. I also want to be able to square the difference which is where a uint32 can be too small. I think the floats look smoother too.

这里有一些关于已经有效的内容以及我想如何改进它的更多背景信息。

澄清一下，我当前的 C 代码计算固定主像素与 8 个或更多邻居像素之间的每 channel 差异。它有一个5层深度的嵌套循环结构: 图像中每个像素的 Y 行和 X 列(3600 万) channel 、R.G 和 B 为 Loop3 循环 4 和 5 用于邻域的行和列。

对于每个 HOME 像素 清除 R、G 和 B 累加器 对于每个邻居，
将abs(home_red - nabr_red) 添加到red_float_accumulator 绿色和蓝色相同 将累积值复制到主内存

我的下一步是将 channel 移至第 5 级，并使用 SIMD 同时执行所有 3 个减法:R、G 和 B。由于每个 MMX 寄存器可使用 48 位/像素和 128 位，因此可以一次完成 2 个操作，而不是仅完成 1 个操作。

使用 Skylake Xeon 上的 AVX2 中的 512 位寄存器，可以完成 10 个。我正在寻找一种好的策略来平衡复杂性与性能，并了解有关这些向量运算的更多信息。

我需要每个“home”像素的 R、G 和 B 累加器。然后将 RGB 移动到与 uint16/ channel RAW、RGB 文件具有相同 XY 分辨率的“ float 图像”中。对每个像素进行相同的对比度计算。

Answer 1

来自Justin W's answer关于无符号整数绝对值的问题:做两次饱和减法。一个结果为零，另一个结果为绝对值。用逻辑或将它们组合起来。在将 16b 整数解包为 32b 整数或 float 之前，这样做是最便宜的。

我们肯定想在从字解包到双字之前进行减法，因此只有一个值需要解包。进行单个非饱和减法然后稍后对其进行排序(例如，在转换为 FP 后将符号位屏蔽为零)是行不通的。存在范围问题:有符号的 int16_t 无法保存两个 uint16_t 值之间差异的完整范围。 UINT16_MAX 将环绕并看起来像 -1 -> 1 的绝对值之差。此外，它还有一个缺点，即需要符号扩展解包。

与 AVX2 一样，解包到不同的向量宽度是最令人头痛的问题，因为某些指令的 channel 内行为和其他指令的跨 channel 行为。

vpunpcklwd 在每个 128b channel 内解包，与 vpackusdw 匹配。 vpmovzxwd ymm, xmm 没有单指令逆，因此我将使用 punpck 指令，假设您可以以这种方式排列 float 。 (对于 PMOVZX，您无法直接从上半部分开始。您必须 vextracti128/vpmovzx。)

#include <immintrin.h>

// untested
__m256 add_pixdiff(__m256 c[2], __m256i a, __m256i b)
{
    __m256i ab = _mm256_subs_epu16(a, b);  // 0        or abs(a-b)
    __m256i ba = _mm256_subs_epu16(b, a);  // abs(a-b) or 0
    __m256i abs_ab_diffs = _mm256_or_si256(ab, ba);

    __m256i lo_uints = _mm256_unpacklo_epi16(abs_ab_diffs, _mm256_setzero_si256());
    __m256i hi_uints = _mm256_unpackhi_epi16(abs_ab_diffs, _mm256_setzero_si256());

    __m256 lo_floats = _mm256_cvtepi32_ps(lo_uints);
    __m256 hi_floats = _mm256_cvtepi32_ps(hi_uints);

    // use load and store intrinsics if the data might not be aligned.
    c[0] = _mm256_add_ps(c[0], lo_floats);
    c[1] = _mm256_add_ps(c[1], hi_floats);
    return c[0];
 }

它compiles exactly like you'd expect ，在上帝 bolt 上。为了在循环中使用，最好手动内联它。您不希望愚蠢的编译器实际上使用内存中的数组来对两个浮点向量进行按引用调用。我只是将其包装在一个函数中以保持简单并隐藏加载/存储。

请注意，一个 uint16 输入向量会生成两个浮点结果向量。我们可以一次对 128b 个整数向量进行操作，但一次对 256b 个整数向量进行操作意味着我们会通过 3 个 insns(不计算负载)+ 2 个解包得到 2 个结果，而不是 6 个 insns + 2 个 pmovzx。这里有相当多的并行性:两个减法可以并行发生，并且有两个解包和转换的依赖链。 (不过，Skylake 只有一个 shuffle 端口，因此您无法一次获得两个解包。这是一条 channel 内指令 with a latency of 1c ，而 vpmovzx ymm, xmm 的延迟为 3，就像其他车道交叉指令一样。如果您需要将 float 保持与整数相同的顺序，而不是最后重新打包回相同的顺序，那么并行性就很重要。)

Vector FP add 的延迟在 Skylake 上增加到 4(之前的英特尔设计为 3)，但每个时钟的吞吐量增加到 2。 (它在 FMA 单元上运行它们。没错，skylake 的 FMA 延迟低至 4c)。

---

我假设您实际上并不想要 SAD(所有差异的一个累加器)，因为您以标量形式编写了 c0，而不是 c。

如果您确实想要一个SAD(绝对差之和)，那就容易得多，并且您应该在整数域中累积。使用相同的 sub 方法和无符号饱和技巧，将它们组合在一起，解压为 32 位 int，然后添加到向量累加器。最后进行水平求和。