使用 SSE2（作为浮点数）缩放字节像素值（y=ax+b）？

Question

我要计算y = ax + b，其中x和y是像素值[即，值范围为0~255的字节]，而a和b是浮点数

由于我需要对图像中的每个像素应用这个公式，此外，a 和 b 对于不同的像素是不同的。 C++中的直接计算速度很慢，所以我有兴趣了解c++中的sse2指令..

经过搜索，发现float与sse2的乘法和加法与_mm_mul_ps和_mm_add_ps一样。但首先我需要将字节中的 x 转换为浮点数（4 字节）。

问题是，在我从字节数据源 (_mm_load_si128) 加载数据后，如何将数据从字节转换为浮点数？

Answer 1

a 和 b 每个像素都不一样？这将使向量化变得困难，除非有模式或者你可以生成它们

有什么方法可以有效地生成向量中的a 和b，无论是定点还是浮点？否则，插入 4 个 FP 值或 8 个 16 位整数可能比仅使用标量操作更糟糕。

---

固定点

如果 a 和 b 可以完全重复使用，或者首先使用定点生成，这可能是定点数学的一个很好的用例。 （即表示值 * 2^scale 的整数）。 SSE/AVX 没有 8b*8b->16b 乘法；最小的元素是单词，因此您必须将字节解压缩为单词，但不能一直到 32 位。这意味着每条指令可以处理两倍的数据。

如果b 和a 的变化不够频繁，或者您可以轻松地生成具有交替a*2^4 和b*2^1 字节的向量，则有一个_mm_maddubs_epi16 指令可能很有用。显然它确实是 handy for bilinear interpolation，但如果我们可以准备 a 和 b 向量，它仍然可以为我们完成工作，只需最少的洗牌。

float a, b;
const int logascale = 4, logbscale=1;
const int ascale = 1<<logascale;  // fixed point scale for a: 2^4
const int bscale = 1<<logbscale;  // fixed point scale for b: 2^1

const __m128i brescale = _mm_set1_epi8(1<<(logascale-logbscale));  // re-scale b to match a in the 16bit temporary result

for (i=0 ; i<n; i+=16) {
    //__m128i avec = get_scaled_a(i);  
    //__m128i bvec = get_scaled_b(i);
    //__m128i ab_lo = _mm_unpacklo_epi8(avec, bvec);
    //__m128i ab_hi = _mm_unpackhi_epi8(avec, bvec);

    __m128i abvec = _mm_set1_epi16( ((int8_t)(bscale*b) << 8) | (int8_t)(ascale*a) );  // integer promotion rules might do sign-extension in the wrong place here, so check this if you actually write it this way.

    __m128i block = _mm_load_si128(&buf[i]);  // call this { v[0] .. v[15] }

    __m128i lo = _mm_unpacklo_epi8(block, brescale);  // {v[0], 8, v[1], 8, ...}
    __m128i hi = _mm_unpackhi_epi8(block, brescale);  // {v[8], 8, v[9], 8, ...
    lo = _mm_maddubs_epi16(lo, abvec);  // first arg is unsigned bytes, 2nd arg is signed bytes
    hi = _mm_maddubs_epi16(hi, abvec);
    // lo = { v[0]*(2^4*a) + 8*(2^1*b), ... }

    lo = _mm_srli_epi16(lo, logascale);  // truncate from scaled fixed-point to integer
    hi = _mm_srli_epi16(hi, logascale);

    // and re-pack.  Logical, not arithmetic right shift means sign bits can't be set
    block = _mm_packuswb(lo, hi);  
    _mm_store_si128(&buf[i], block);
}
// then a scalar cleanup loop

2^4 是任意选择。它为 a 的整数部分留下 3 个非符号位和 4 个小数位。所以它有效地将a 舍入到最接近的第 16 位，如果它的幅度大于 8 和 15/16 则溢出。 2^6 将提供更多小数位，并允许 a 从 -2 到 +1 和 63/64ths。

由于b 是相加而不是相乘，它的有用范围要大得多，而小数部分的用处要小得多。用 8 位表示，四舍五入到最接近的一半仍然保留一点小数信息，但允许它为 [-64 : 63.5] 而不会溢出。

为了获得更高的精度，16b 定点是一个不错的选择。您可以将a 和b 向上缩放 2^7 或其他东西，以获得 7b 的小数精度，并且仍然允许整数部分为 [-256 .. 255]。这种情况下没有乘加指令，因此您必须单独执行此操作。做乘法的好选择包括：

• _mm_mulhi_epu16：无符号 16b*16b->high16（位 [31:16]）。如果a 不能为负数，则很有用
• _mm_mulhi_epi16：有符号 16b*16b->high16（位 [31:16]）。
• _mm_mulhrs_epi16：带符号的 16b*16b->bits [30:15] 的 32b 临时，四舍五入。为a 选择一个很好的比例因子，这应该会更好。据我了解，SSSE3 正是针对这种用途引入了这条指令。
• _mm_mullo_epi16：有符号 16b*16b->low16（位 [15:0]）。在 low16 结果溢出之前，这仅允许 a 的 8 个有效位，所以我认为您从 _mm_maddubs_epi16 8 位解决方案中获得的所有收益对于 b 来说更精确。

要使用这些，您将获得 a 和 b 值的缩放 16b 向量，然后：

• 用零（或pmovzx byte->word）解压你的字节，得到仍在[0..255]范围内的有符号字
• 将单词左移 7 位。
• 乘以你的a 16b 个单词的向量，取每个 16*16->32 结果的上半部分。 （例如 mul
• 如果您想要a 和b 的不同比例，请在此处右移，以获得a 的更多小数精度
• 添加b。
• 右移将最后的截断从固定点返回到 [0..255]。

有了定点刻度的好选择，这应该能够处理更广泛的a和b，以及比8位定点更高的小数精度。

如果在将字节解包为单词后不左移字节，a 必须是全范围的，才能在结果的高 16 位中设置 8 位。这意味着您可以支持的a 范围非常有限，而不会在乘法期间将您的临时值截断为小于 8 位。即使_mm_mulhrs_epi16 也不会留下太多空间，因为它从第 30 位开始。

---

将字节扩展为浮点数

如果您不能有效地为每个像素生成定点 a 和 b 值，最好将像素转换为浮点数。这需要更多的拆包/重新打包，因此延迟和吞吐量更差。值得研究生成具有固定点的 a 和 b。

要让packed-float 工作，您仍然必须为4 个相邻像素有效地构建一个a 值的向量。

这是pmovzx (SSE4.1) 的一个很好的用例，因为它可以直接从 8b 元素到 32b。其他选项是具有多个步骤的 SSE2 punpck[l/h]bw/punpck[l/h]wd，或模拟 pmovzx 的 SSSE3 pshufb。 （你可以做一个 16B 的加载，并用 4 种不同的方式将它解压成四个 32b 整数的向量。）

char *buf;

// const __m128i zero = _mm_setzero_si128();
for (i=0 ; i<n; i+=16) {
    __m128 a = get_a(i);
    __m128 b = get_b(i);

    // IDK why there isn't an intrinsic for using `pmovzx` as a load, because it takes a m32 or m64 operand, not m128.  (unlike punpck*)
    __m128i unsigned_dwords = _mm_cvtepu8_epi32((__m128i)(buf+i));  // load 4B at once.
    __m128 floats = _mm_cvtepi32_ps(unsigned_dwords);

    floats = _mm_fmadd_ps(floats, a, b);  // with FMA available, this might as well be 256b vectors, even with the inconvenience of the different lane-crossing semantics of pmovzx vs. punpck
    // or without FMA, do this with _mm_mul_ps and _mm_add_ps
    unsigned_dwords = _mm_cvtps_epi32(floats);

    // repeat 3 more times for buf+4, buf+8, and buf+12, then:

    __m128i packed01 = _mm_packss_epi32(dwords0, dwords1);  // SSE2
    __m128i packed23 = _mm_packss_epi32(dwords2, dwords3);
    // packuswb wants SIGNED input, so do signed saturation on the first step

    // saturate into [0..255] range
    __m12i8 packedbytes=_mm_packus_epi16(packed01, packed23);  // SSE2
    _mm_store_si128(buf+i, packedbytes);  // or storeu if buf isn't aligned.
}

// cleanup code to handle the odd up-to-15 leftover bytes, if n%16 != 0

此答案的先前版本来自带有 packusdw/packuswb 的 float->uint8 向量，并且有一整节介绍了没有 SSE4.1 的解决方法。如果您只是stay in the signed integer domain until the last pack，则不需要在无符号包之后屏蔽符号位。我认为这就是 SSE2 仅包含从 dword 到 word 的签名包，但包含从 word 到 byte 的有符号和无符号包的原因。 packuswd 仅在您的最终目标是 uint16_t 时才有用，而不是进一步打包。

---

最后一个没有 SSE4.1 的 CPU 是 Intel Conroe/merom（第一代 Core2，从 2007 年底之前）和 AMD 之前的巴塞罗那（2007 年底之前）。如果这些 CPU 可以接受工作但速度慢，只需为 AVX2 编写一个版本，为 SSE4.1 编写一个版本。或者 SSSE3（使用 4x pshufb 来模拟寄存器的四个 32b 元素的 pmovzxbd）pshufb 在 Conroe 上很慢，所以如果您关心没有 SSE4.1 的 CPU，请编写特定版本。实际上，Conroe/merom 也有慢 xmm punpcklbw 等等（q->dq 除外）。 4x 慢 pshufb 仍然应该比 6x 慢解包要好。由于拆包和重新包装的缓慢洗牌，矢量化在沃尔夫代尔之前的胜利要少得多。定点版本，解包/重新打包少得多，将有更大的优势。

在我意识到它需要多少额外指令之前，请查看编辑历史以了解使用 punpck 的未完成尝试。删除它是因为这个答案已经很长了，另一个代码块会令人困惑。

Answer 2

我猜你正在寻找 __m128 _mm_cvtpi8_ps(__m64 a ) 复合内在函数。

这是一个最小的例子：

#include <xmmintrin.h>
#include <stdio.h>

int main() {
  unsigned char  a[4] __attribute__((aligned(32)))= {1,2,3,4};
  float b[4] __attribute__((aligned(32)));
  _mm_store_ps(b, _mm_cvtpi8_ps(*(__m64*)a));
  printf("%f %f, %f, %f\n", b[0], b[1], b[2], b[3]);
  return 0;
}