3

有人对如何使用 SSE/AVX(内在函数或程序集)实现Lanczos 图像重采样(放大和缩小)算法有任何提示吗?

我查看了一些 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 并且根本不转换为实数(错误可能很大)?

4

1 回答 1

3

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

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

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

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


至于将图像数据处理为浮点数与整数,那么如果您可以摆脱 16 位定点,那会更快(除非它需要更多指令)。请参阅我对另一个关于使用浮点数与定点数的图像处理问题的回答。


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

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

如果您只需要计算sin(x*PI),您可以创建一个自定义sin函数来执行此操作,从而省去预乘 PI 的麻烦。由于理想的实现会sin 根据输入的范围选择要使用的算法,因此您可能无法获得精确到尾数最后一位的矢量化实现。幸运的是,你可能不需要那个,所以用谷歌搜索一个 SSE sin(x) 实现。


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

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

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与 SSE 版本相比,AVX VEX 编码版本中的谓词选择范围更广。

于 2015-12-11T00:35:36.510 回答