Lanczos SSE/AVX 实施

Posted 2023-02-16

【中文标题】Lanczos SSE/AVX 实施

【英文标题】：Lanczos SSE/AVX implementation

【发布时间】：2015-12-10 23:15:45

【问题描述】：

有人对如何使用 SSE/AVX（内在函数或汇编）实现Lanczos image resampling（放大和缩小）算法有任何提示吗？

我查看了一些 C 实现，但有很多分支，我真的不知道如何使用 SSE/AVX 巧妙地实现它。

示例 - 归一化的大罪：

// C implementation
if (!x)
  return sin(x*M_PI)/(x*M_PI);
else
  return 1;

// AVX implementation
PXOR ymm0, ymm0
MOVAPD ymm1, [x]     // x - array of double
CMPPD ymm0, ymm1, 0  // if (!x)
// what now?

MOVAPD ymm3, [pi]    // pi - array of double = M_PI - is there better way?
PMULPD ymm1, ymm3    // ymm1 = x*pi
(SINPD ymm2, ymm1)   // found intrinsic _mm256_sin_pd - Intel math library, intrinsic functions are OK with me
DIVPD ymm2, ymm1     // result in ymm2

对于值 x == 0，我应该如何返回 1？在那些索引上，我在 CMPPD 之后有 11...11（真）。

另外，我正在为灰度、8 位图片执行此操作，所以一个像素只有 (0..255)。使用 float 而不是 double 会对质量产生什么影响？是否可以一直使用 u_int8 并且根本不转换为实数（错误可能很大）？

【问题讨论】：

SSE 中唯一的数学指令（除了基本的 add/sub/mul/div）是 sqrt。 Trig / log / exp 都是库函数。请注意，在英特尔的内在指南中，“指令”字段是空白的，这意味着它映射到多个指令。只有英特尔的编译器甚至提供这些复合内在函数。

您是否需要基于每个像素的任意x，即您真的需要为每个像素计算新的 Lanczos 系数吗？ （这个问题的答案取决于您尝试执行什么样的重采样或插值。）

我需要调整图像大小 - 放大和缩小。我是这方面的新手，仍然试图掌握这一切是如何运作的。我有信号处理的基础知识（例如，我知道卷积是如何工作的，等等）。这是参考算法的主要部分（垂直跟随的类似代码）：pastebin.com/XtUspjGW

@PeterCordes 是的，我也很害怕。但这不应该是一个问题。我正在使用 MS Visual Studio 和其他内在函数（用于说明）工作。如果它也适用于他们的数学库，我将不得不尝试。现在当我考虑它时 - 是否有可能找出它映射到什么指令？这意味着我可以用任何东西编译它。但它可能是闭源的，对吧？

您可能想要使用“足够准确”的基于表格的 Lanczos 内核实现。否则性能将很大程度上取决于你的 sin 函数的速度。

【参考方案1】：

如果您还不了解 asm 或 SSE/AVX，一次学习一个可能会更容易。使用 C/C++ 内在函数编写向量算法将为您提供比直接使用 asm 更可移植的实现。 （针对 32 位和 64 位以及 windows 或其他所有版本进行编译，而不需要 2 或 4 个不同的 asm 版本（或 asm 中的 #ifdef-equivalent 宏）。

在编写 C 时查看编译器输出可能会有所帮助，以确保加载/存储按您期望的方式进行，并且编译器不会因为别名/对齐而对臃肿的代码做任何愚蠢的事情（缺乏of) 假设，或存储/生成常量。

调试矢量代码已经够难的了（因为要跟踪的状态要多得多，而且你必须通过随机播放在精神上跟随事物）。

如果编译器还没有自动矢量化，我会先找到 C 中可能矢量化的某些部分，然后在 C 中使用内在函数。然后，一旦工作正常，我可能会获取编译器的输出和手- 在编译器没有生成最佳代码的地方调整它。 （见http://agner.org/optimize/）

---

至于将图像数据处理为浮点数与整数，如果您可以摆脱 16 位定点，那会更快（除非它需要更多指令）。请参阅 my answer on another image-processing question 了解使用浮点数与定点数。

---

SSE 中唯一的数学指令（除了基本的 add/sub/mul/div）是 sqrt。 Trig / log / exp 都是库函数。请注意，在英特尔的内在指南中，“指令”字段是空白的，这意味着它映射到多个指令。只有英特尔的编译器甚至提供这些复合内在函数。

无论如何，您需要找到一个sin 实现来内联，或者保存一些寄存器并进行函数调用。根据 ABI（windows 或其他一切），某些或所有 xmm 寄存器可以被函数破坏。使用特定的sin 实现会让您知道它需要哪些寄存器，并且只会溢出那些。 （由于您是在 asm 中编程，因此您可以使函数彼此了解更多，而不是仅仅遵循 ABI。）

如果您只需要 calculate sin(x*PI)，您可以创建一个自定义的 sin 函数来执行此操作，从而省去预乘 PI 的麻烦。由于sinchooses what algorithm to use based on the range of the input 的理想实现，您可能无法获得精确到尾数最后一位的矢量化实现。幸运的是，你可能不需要那个，所以用谷歌搜索一个 SSE sin(x) 实现。

---

SIMD 向量中的条件是通过比较生成一个全零或全一元素的向量来处理的。然后，您可以将这些与或或与其他向量相结合。它适用于添加标识值为0 的位置。 （x + 0 = x，因此您可以在将向量添加到累加器之前从向量中过滤掉一些元素）。如果您需要基于 0 / -1 的向量在两个源元素之间进行选择，您可以 AND/OR 一起使用，或者使用 blendvps（变量混合打包标量，而不是编译 -时间常数混合）。

如果您想避免首先计算缓慢的除以零，而不是通常只对所有内容进行计算然后进行掩蔽/混合，那么这个想法就会有点失败。由于您希望在x == 0.0 时将结果输出到1，因此最好的选择可能是将x 的零元素设置为FLT_MIN * 16 或其他值，然后再计算sin(x*PI)/(x*PI)。这样，您可以避免除以零，并且除法的结果非常接近 1。如果您需要它精确到 1.0f（并且没有 x 的值使 sin(x*PI) == x*PI使用您的sin 实现），那么您需要混合两次：一次在分子中，一次在分母中。 （将它们都设置为相同的非零值）。

vxorps     xmm15, xmm15, xmm15   ; if you can spare a reg to hold a zero constant



; inside your loop:  xmm0 holds  x3, x2, x1, x0 .
vcmpeqps     xmm1, xmm0, xmm15   ;; mnemonic for vcmpps xmm1, xmm0, xmm15, 0.
;;  Different predicates are an immediate operand, not different opcodes


vblendvps  xmm2, xmm0, [memory_holding_vector_of_float_min], xmm1  ; Or cache it in a reg if you have one to spare
     ; blendv takes elements from the 2nd source operand when the selector (xmm1) has a 1-bit in the MSB (sign bit)

; xmm2 = (x==0.0f) ? FLT_MIN : x
;  xmm1 holds  sin(x3*pi), sin(x2*pi)... 

请注意，cmpps 在 AVX VEX 编码版本中比在 SSE 版本中具有更广泛的谓词选择。