Lanczos SSE/AVX实现

Question

有没有人知道如何使用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)。使用浮子而不是双倍对质量有什么影响？是否可以一直使用u_int8而根本不转换为实数(错误可能很大)？

Answer 1

如果您还不了解asm或SSE/AVX，那么一次学习可能更容易。与直接使用asm相比，用C/C++本质编写向量算法将给您一个更可移植的实现。(编译32 vs.64位，windows或其他任何东西，而不需要2或4个不同的asm版本(或者asm中的#ifdef-等效宏)。

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

调试向量代码已经足够困难了(因为有更多的状态需要跟踪，而且您必须在头脑中通过洗牌来跟踪事情)。

首先，我将从C中找到一些可能向量化的部分，如果编译器还没有实现自动向量化，然后在C中使用内禀函数，然后，一旦这一方法生效，我就可以使用编译器输出，并在编译器没有生成最优代码的地方手动调优它。(见http://agner.org/optimize/)

至于像浮点和int一样处理图像数据，好吧，如果你能摆脱16位定点，那会更快(除非它需要更多的指令)。请参阅我对另一个图像处理问题的回答关于使用浮点与定点的关系。

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

无论如何，您需要找到一个内联的sin实现，或者保存一些寄存器并进行函数调用。根据ABI (windows或其他一切)，一些或所有xmm寄存器可以被函数关闭。使用特定的sin实现会让您知道它需要哪些寄存器，并且只会泄漏这些寄存器。(因为您是用asm编程的，所以您可以使函数相互了解，而不是仅仅遵循ABI。)

如果你只需要sin(x*PI)，你可以做一个定制的sin函数，这样就省去了预先乘以PI的麻烦。由于sin 根据输入范围选择要使用的算法。是一个理想的实现，所以您可能无法获得精确到尾数最后一点的向量化实现。幸运的是，您可能不需要这样做，所以google上有一个sin(x)实现。

SIMD向量中的条件是通过比较来处理的，这些元素的向量要么是全零的，要么是全一的。你可以然后和或或那些到其他向量。对于添加标识值为0的内容，它工作得很好。(x + 0 = x，因此在将向量添加到累加器之前，可以从向量中筛选出一些元素)。如果您需要基于0/ -1的向量在两个源元素之间进行选择，那么可以将它们放在一起，或者使用blendvps (变量混合填充标量，而不是编译时常量混合)来更快地完成相同的工作。

如果你想避免在一开始就计算一个缓慢的零除法，而不是通常的只做所有的计算，然后掩蔽/混合，那么这个想法就会被分解一些。由于您希望在1上显示x == 0.0时的结果，所以最好的选择是在计算任何sin(x*PI)/(x*PI)之前将x的零元素设置为FLT_MIN * 16或其他什么。这样，你就可以避免除以零，除法的结果接近于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编码版本中有比SSE版本更广泛的谓词选择。