Lanczos SSE/AVX 实现

Question

有人对如何使用 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，我应该如何 return 1？在这些索引上，我在 CMPPD 之后有 11...11 (true)。

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

Answer 1

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

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

调试矢量代码已经够难的了（因为要跟踪的状态很多，而且您必须在头脑中通过随机播放来跟进）。

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

就将图像数据处理为 float 与 int 相比，如果您可以使用 16 位定点，那会更快（除非它需要更多指令）。请参阅关于使用浮点数与定点数。

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

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

如果您只需要 calculate sin(x*PI), you could make a custom sin function that does that, saving the trouble of pre-multiplying by PI. Since an ideal implementation of sin chooses what algorithm to use based on the range of the input，您可能无法获得精确到尾数最后一位的矢量化实现。幸运的是，你可能不需要那个，所以 google 建立一个 SSE sin(x) 实现。

SIMD 向量中的条件通过比较来处理，从而使元素向量为全零或全一。然后您可以将它们与其他向量进行 AND 或 OR 运算。它适用于添加标识值为 0 的地方。（x + 0 = x，因此您可以在将向量添加到累加器之前从向量中过滤掉一些元素）。如果您需要 select 基于向量 0 / -1 的两个源元素之间，您可以 AND/OR 东西在一起，或者使用 blendvps 更快地完成相同的工作（可变混合压缩标量，而不是编译时常量混合）。

如果您想首先避免计算缓慢的除以零，而不是通常只对所有内容进行计算然后 masking/blending，那么这个想法有点不合理。由于您希望 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 版本中有更多的谓词选择。

Lanczos SSE/AVX 实现

Lanczos SSE/AVX implementation

assembly

sse

image-processing

avx

lanczos