高效地去交错并将 float 转换为 uint16_t
Deinterleave and convert float to uint16_t efficiently
我需要将浮点的打包图像缓冲区 (YUVA) 解交织到平面缓冲区。我也想将这些 float 转换为 uint16_t，但这真的很慢。我的问题是：如何使用内在函数加快速度？
void deinterleave(char* pixels, int rowBytes, char *bufferY, char *bufferU, char *bufferV, char *bufferA)
{
// Scaling factors (note min. values are actually negative) (limited range)
const float yuva_factors[4][2] = {
{ 0.07306f, 1.09132f }, // Y
{ 0.57143f, 0.57143f }, // U
{ 0.57143f, 0.57143f }, // V
{ 0.00000f, 1.00000f } // A
};
float *frameBuffer = (float*)pixels;
// De-Interleave and convert source buffer / bottom first
for (int r = height - 1, p = 0; r >= 0; r--)
{
for (int c = 0; c < width; c++)
{
// Get beginning of next block
const int pos = r * width * 4 + c * 4;
// VUYA -> YUVA
((uint16_t*)bufferY)[p] = (uint16_t)((frameBuffer[pos + 2] + yuva_factors[0][0]) / (yuva_factors[0][0] + yuva_factors[0][1]) * 65535.0f);
((uint16_t*)bufferU)[p] = (uint16_t)((frameBuffer[pos + 1] + yuva_factors[1][0]) / (yuva_factors[1][0] + yuva_factors[1][1]) * 65535.0f);
((uint16_t*)bufferV)[p] = (uint16_t)((frameBuffer[pos + 0] + yuva_factors[2][0]) / (yuva_factors[2][0] + yuva_factors[2][1]) * 65535.0f);
((uint16_t*)bufferA)[p] = (uint16_t)((frameBuffer[pos + 3] + yuva_factors[3][0]) / (yuva_factors[3][0] + yuva_factors[3][1]) * 65535.0f);
p++;
}
}
}
澄清一下：我从这个 API 函数中得到 "pixels" 缓冲区 ...
// prSuiteError (*GetPixels)(PPixHand inPPixHand, PrPPixBufferAccess inRequestedAccess, char** outPixelAddress);
char *pixels;
ppixSuite->GetPixels(inRenderedFrame, PrPPixBufferAccess_ReadOnly, &pixels);
... 根据所选的像素格式，它可以是从 uint8_t 到 float 的任何值。在这个用例中，它肯定是 floats.
我的简化代码如下所示：
#include <stdint.h>
static const int width = 1920;
static const int height = 1080;
void unpackFloatToUint16(float* pixels, uint16_t *bufferY, uint16_t *bufferU, uint16_t *bufferV, uint16_t *bufferA)
{
for (int r = height - 1; r >= 0; r--)
{
for (int c = 0; c < (int)width * 4; c += 4)
{
const int pos = r * width * 4 + c;
*bufferV++ = (uint16_t)((pixels[pos] + 0.57143f) * 57342.98164f);
*bufferU++ = (uint16_t)((pixels[pos + 1] + 0.57143f) * 57342.98164f);
*bufferY++ = (uint16_t)((pixels[pos + 2] + 0.07306f) * 56283.17216f);
*bufferA++ = (uint16_t)(pixels[pos + 3] * 65535.0f);
}
}
}
有一件事很明显：用乘法代替除法。如果您在 FP 划分上遇到瓶颈，吞吐量应该增加 ~7 到 ~15 倍左右。 （Haswell 的 divss 每 7 个时钟吞吐量一个，但 mulss 是每 0.5 个时钟一个）。
（假设您还没有使用 -ffast-math 来让编译器为您用倒数乘法代替常量除法）。
GCC 和 clang 已经自动向量化你的函数（至少使用编译时常量 height 和 width。问题中的代码没有编译，因为它们没有定义，所以我不知道该假设什么）See it on the Godbolt compiler explorer。在没有 -ffast-math 的情况下，它确实使用 divps 作为除法，但它确实使用 SIMD 进行数学运算（包括转换为 32 位整数），并通过洗牌将第 4 组 Y 值混合在一起对于 64 位商店。我不认为它的工作效率很高，但如果您在 div 吞吐量方面遇到瓶颈，那么它可能比 gcc 好得多。
但是对于 int height = rowBytes >> 1;，clang 不会自动矢量化，而 gcc 仍然可以。
不过，编译器的工作似乎还有改进的余地。
无论如何，假设我们要为 AVX + FMA（例如 Haswell 或 Steamroller / Ryzen）手动矢量化。您也可以制作其他版本，但由于您没有指定任何关于您想要定位的微体系结构（甚至是 x86），我将作为一个有趣的例子来做。
首先，我们可以将(Y + c0) / (c0 + c1) * 65535.0f转化为单个FMA。将加法里面的* (1.0f/(c0+c1)) * 65535.0f分布得到(Y * mul_const + add_const)，可以用单个FMA求值。我们可以一次对一个像素的所有 4 个分量执行此操作，使用 128 位 SIMD FMA 和两个向量常量以匹配内存中 float 布局的顺序保存系数。 （或者使用 256 位 FMA 一次处理所有 2 个像素）。
不幸的是，gcc 和 clang 没有使用 -ffast-math 进行此优化。
在转换为整数之前保存所有改组可能效果最好。它的优点是您只需要两个 FP 常量向量，而不是将每个分量的系数广播到所有元素的单独向量。好吧，我猜你可以在转换为整数之前对 FMA 的结果使用 FP 洗牌。 （例如FMA，然后洗牌，然后将4个Y值的向量转换为整数）。
FP 数学不是严格关联或分配的（因为舍入误差）所以这可能会改变边缘情况的结果。但不一定会使情况变得更糟，只是与您以前的舍入方式不同。如果使用 FMA 将 (Y + const1) * const2 转换为 Y * altconst1 + altconst2 不会丢失任何精度，因为 FMA 在进行加法之前不会舍入内部临时产品。
所以我们知道如何有效地进行数学运算并转换为整数（2 CPU 条指令用于包含 2 个像素的 8 float 向量）。剩下的就是将 Y 与其他 Y 组合在一起，并将 32 位有符号整数压缩为 16 位无符号整数。 （x86 只能在 FP 和有符号整数之间转换，直到 AVX512F 引入直接 FP <-> 无符号（和 SIMD 用于 64 位整数 <-> FP 而不是 64 位模式下的标量）。无论如何，我们不能直接从 float 转换为 16 位无符号整数向量）。
因此，给定一个包含 VUYA 个 32 位整数元素的 128 位向量，我们的第一步可能是缩小到 16 位整数。 x86 有一条指令（SSE4.1 packusdw，内在 _mm_packus_epi32）来打包无符号饱和（因此负输入饱和为 0，而大的正输入饱和为 65535）。大概这就是您想要的，而不是 t运行cating 会使溢出回绕的整数。这需要 2 个 SIMD 向量作为输入，并产生一个输出向量，所以我们会得到 VYUAVYUA（对于 2 个不同的像素）。
即使您不需要饱和行为（例如，如果输入不可能超出范围），packusdw 可能仍然是缩小整数范围的最有效选择。其他洗牌只有一个输入向量，或固定的洗牌模式，它不会丢弃每个 32 位元素的上半部分，因此在 pshufb 之后，结果中只有 64 位有用数据或 punpck 随机播放。
立即从 pack 开始可以很好地将事情减少到只有 2 个向量。我查看了其他洗牌顺序，但如果您从 32 位洗牌开始而不是压缩到 16 位元素，则总是需要更多的总洗牌。 （见下方神马link中的评论）
vpackusdw 的 256 位 AVX2 版本分别在 256 位向量的两个 128 位通道上运行，因此对于 pack(ABCD EFGH, IJKL MNOP) 你得到 ABCD IJKL EFGH MNOP。通常你需要另一次洗牌才能把事情按正确的顺序排列。无论如何，我们需要进一步洗牌，但这仍然很麻烦。不过，我认为您可以在每次循环迭代中处理两倍的数据，而循环中只需再洗几下。
这是我使用 128 位向量得出的结果
源代码 + 编译器输出 on the Godbolt compiler explorer
请注意，它不处理像素数不是 4 的倍数的情况。您可以使用清理循环（加载、缩放和填充饱和度，然后提取四个16 位组件）。或者你可以做一个部分重叠的最后 4 个像素。 （如果像素数确实是 4 的倍数，则不重叠，否则部分重叠存储到 Y、U、V 和 A 数组中。）这很容易，因为它不是就地操作，所以你可以重新读取相同的存储输出后输入。
此外，它假定行步幅与宽度相匹配，因为您在问题中的代码做了同样的事情。因此，宽度是否为 4 像素的倍数并不重要。但是，如果您确实有一个与宽度分开的可变行步幅，则您将不得不担心每行末尾的清理。 （或者使用填充，这样你就不必这样做了）。
#include <stdint.h>
#include <immintrin.h>
static const int height = 1024;
static const int width = 1024;
// helper function for unrolling
static inline __m128i load_and_scale(const float *src)
{ // and convert to 32-bit integer with truncation towards zero.
// Scaling factors (note min. values are actually negative) (limited range)
const float yuvaf[4][2] = {
{ 0.07306f, 1.09132f }, // Y
{ 0.57143f, 0.57143f }, // U
{ 0.57143f, 0.57143f }, // V
{ 0.00000f, 1.00000f } // A
};
// (Y + yuvaf[n][0]) / (yuvaf[n][0] + yuvaf[n][1]) ->
// Y * 1.0f/(yuvaf[n][0] + yuvaf[n][1]) + yuvaf[n][0]/(yuvaf[n][0] + yuvaf[n][1])
// Pixels are in VUYA order in memory, from low to high address
const __m128 scale_mul = _mm_setr_ps(
65535.0f / (yuvaf[2][0] + yuvaf[2][1]), // V
65535.0f / (yuvaf[1][0] + yuvaf[1][1]), // U
65535.0f / (yuvaf[0][0] + yuvaf[0][1]), // Y
65535.0f / (yuvaf[3][0] + yuvaf[3][1]) // A
);
const __m128 scale_add = _mm_setr_ps(
65535.0f * yuvaf[2][0] / (yuvaf[2][0] + yuvaf[2][1]), // V
65535.0f * yuvaf[1][0] / (yuvaf[1][0] + yuvaf[1][1]), // U
65535.0f * yuvaf[0][0] / (yuvaf[0][0] + yuvaf[0][1]), // Y
65535.0f * yuvaf[3][0] / (yuvaf[3][0] + yuvaf[3][1]) // A
);
// prefer having src aligned for performance, but with AVX it won't help the compiler much to know it's aligned.
// So just use an unaligned load intrinsic
__m128 srcv = _mm_loadu_ps(src);
__m128 scaled = _mm_fmadd_ps(srcv, scale_mul, scale_add);
__m128i vuya = _mm_cvttps_epi32(scaled); // truncate toward zero
// for round-to-nearest, use cvtps_epi32 instead
return vuya;
}
void deinterleave_avx_fma(char* __restrict pixels, int rowBytes, char *__restrict bufferY, char *__restrict bufferU, char *__restrict bufferV, char *__restrict bufferA)
{
const float *src = (float*)pixels;
uint16_t *__restrict Y = (uint16_t*)bufferY;
uint16_t *__restrict U = (uint16_t*)bufferU;
uint16_t *__restrict V = (uint16_t*)bufferV;
uint16_t *__restrict A = (uint16_t*)bufferA;
// 4 pixels per loop iteration, loading 4x 16 bytes of floats
// and storing 4x 8 bytes of uint16_t.
for (unsigned pos = 0 ; pos < width*height * 4; pos += 4) {
// pos*4 because each pixel is 4 floats long
__m128i vuya0 = load_and_scale(src+pos*4);
__m128i vuya1 = load_and_scale(src+pos*4 + 1);
__m128i vuya2 = load_and_scale(src+pos*4 + 2);
__m128i vuya3 = load_and_scale(src+pos*4 + 3);
__m128i vuya02 = _mm_packus_epi32(vuya0, vuya2); // vuya0 | vuya2
__m128i vuya13 = _mm_packus_epi32(vuya1, vuya3); // vuya1 | vuya3
__m128i vvuuyyaa01 = _mm_unpacklo_epi16(vuya02, vuya13); // V0V1 U0U1 | Y0Y1 A0A1
__m128i vvuuyyaa23 = _mm_unpackhi_epi16(vuya02, vuya13); // V2V3 U2U3 | Y2Y3 A2A3
__m128i vvvvuuuu = _mm_unpacklo_epi32(vvuuyyaa01, vvuuyyaa23); // v0v1v2v3 | u0u1u2u3
__m128i yyyyaaaa = _mm_unpackhi_epi32(vvuuyyaa01, vvuuyyaa23);
// we have 2 vectors holding our four 64-bit results (or four 16-bit elements each)
// We can most efficiently store with VMOVQ and VMOVHPS, even though MOVHPS is "for" FP vectors
// Further shuffling of another 4 pixels to get 128b vectors wouldn't be a win:
// MOVHPS is a pure store on Intel CPUs, no shuffle uops.
// And we have more shuffles than stores already.
//_mm_storeu_si64(V+pos, vvvvuuuu); // clang doesn't have this (AVX512?) intrinsic
_mm_storel_epi64((__m128i*)(V+pos), vvvvuuuu); // MOVQ
_mm_storeh_pi((__m64*)(U+pos), _mm_castsi128_ps(vvvvuuuu)); // MOVHPS
_mm_storel_epi64((__m128i*)(Y+pos), yyyyaaaa);
_mm_storeh_pi((__m64*)(A+pos), _mm_castsi128_ps(yyyyaaaa));
}
}
希望 shuffle 的变量名 + 注释应该是相当可读的。这是未经测试的；最可能的错误是某些向量作为随机播放参数的错误排序。但是修复它应该只是颠倒 arg 顺序之类的问题，而不需要额外的洗牌来减慢它的速度。
看起来包括包装在内的 6 次洗牌是我能做的最好的了。它基本上包括一个从 vuya x4 到 vvvv uuuu yyyy aaaa 的 4x4 转置，并且 pack 有一个固定的洗牌模式，这对去交错没有帮助，所以我认为我们不能做得比这更好128 位向量。当然，我总是有可能忽略了一些东西。
gcc 和 clang 都编译它略微次优:
Clang 使用 vpextrq 而不是 vmovhps（在 Intel CPUs 上，每次循环迭代总共花费额外的 2 个 shuffle uops）。并且还使用 2 个单独的循环计数器，而不是将同一个计数器缩放 1 或 8，这样会花费 1 个额外的整数 add 指令，没有任何好处。 （要是 gcc 选择这样做而不是使用折叠到 FMA 中的索引加载就好了……愚蠢的编译器。）
gcc，而不是使用 vmovups 负载，处理 FMA3 通过复制向量常量然后使用具有索引寻址模式的内存操作数来破坏其输入。 This doesn't stay micro-fused，所以前端总共有 4 个额外的微指令。
如果编译完美，仅使用一个循环计数器作为 float 源数组索引（按 *8 缩放）和整数目标数组（未缩放），就像 gcc 那样，并且如果 gcc 以 clang 的方式加载，那么整个循环将是 Haswell/Skylake.
上的 24 个融合域 uops
因此它可以每 6 个时钟从前端发出一次迭代，恰好在每个时钟 4 个融合域微指令的限制下。它包含 6 个 shuffle uops，因此它也可以解决 port5 吞吐量瓶颈。 （HSW / SKL 只有 1 个洗牌执行单元）。所以 理论上 ，这可以 运行 每 6 个时钟 4 个像素（16 个浮点数），或者在英特尔 CPUs 上每 1.5 个时钟周期一个像素。在 Ryzen 上可能稍微多一些，虽然 MOVHPS 在那里花费了多个 uops，但管道更宽。 （参见 http://agner.org/optimize/）
每 6 个时钟 4 次加载和 4 次存储很好，并且远离任何瓶颈，但如果您的 src 和 dst 在缓存中不热，则可能存在内存带宽。商店只有一半宽度，一切都是顺序的。来自一个循环的 4 个独立输出流已经足够少，通常不会成为问题。 （如果你有超过 4 个输出流，你可能会考虑裂变循环并只存储一次传递的 U 和 V 值，然后在另一次传递中只存储 Y 和 A 值传递相同的源数据，或者类似的东西。但是就像我说，4 个输出流就可以了，不保证循环裂变。）
这个的 256 位版本最后可能需要的不仅仅是 2 个额外的 vpermq，因为我认为你不能轻易解决你想要相邻的两个 V 值 [= =49=] 卡在同一个 __m256i 向量的高低车道。因此，您可能需要在此过程的早期进行额外的 4 vpermd 或 vperm2i128 洗牌，因为最小粒度的跨车道洗牌是 32 位粒度。这可能会对 Ryzen 造成很大伤害。
在将 4 或 8 个 v 组合在一起后，也许您可以使用 vpblendw 重新排列向量之间的词元素，但 v 是错误的。
，所以我认为这个的 AVX512 版本会更便宜。 AVX512 具有收缩饱和/t运行cate 指令，但这些指令在 Skylake-X 上比 vpackusdw 慢，所以可能不会。
我需要将浮点的打包图像缓冲区 (YUVA) 解交织到平面缓冲区。我也想将这些 float 转换为 uint16_t，但这真的很慢。我的问题是：如何使用内在函数加快速度？
void deinterleave(char* pixels, int rowBytes, char *bufferY, char *bufferU, char *bufferV, char *bufferA)
{
// Scaling factors (note min. values are actually negative) (limited range)
const float yuva_factors[4][2] = {
{ 0.07306f, 1.09132f }, // Y
{ 0.57143f, 0.57143f }, // U
{ 0.57143f, 0.57143f }, // V
{ 0.00000f, 1.00000f } // A
};
float *frameBuffer = (float*)pixels;
// De-Interleave and convert source buffer / bottom first
for (int r = height - 1, p = 0; r >= 0; r--)
{
for (int c = 0; c < width; c++)
{
// Get beginning of next block
const int pos = r * width * 4 + c * 4;
// VUYA -> YUVA
((uint16_t*)bufferY)[p] = (uint16_t)((frameBuffer[pos + 2] + yuva_factors[0][0]) / (yuva_factors[0][0] + yuva_factors[0][1]) * 65535.0f);
((uint16_t*)bufferU)[p] = (uint16_t)((frameBuffer[pos + 1] + yuva_factors[1][0]) / (yuva_factors[1][0] + yuva_factors[1][1]) * 65535.0f);
((uint16_t*)bufferV)[p] = (uint16_t)((frameBuffer[pos + 0] + yuva_factors[2][0]) / (yuva_factors[2][0] + yuva_factors[2][1]) * 65535.0f);
((uint16_t*)bufferA)[p] = (uint16_t)((frameBuffer[pos + 3] + yuva_factors[3][0]) / (yuva_factors[3][0] + yuva_factors[3][1]) * 65535.0f);
p++;
}
}
}
澄清一下：我从这个 API 函数中得到 "pixels" 缓冲区 ...
// prSuiteError (*GetPixels)(PPixHand inPPixHand, PrPPixBufferAccess inRequestedAccess, char** outPixelAddress);
char *pixels;
ppixSuite->GetPixels(inRenderedFrame, PrPPixBufferAccess_ReadOnly, &pixels);
... 根据所选的像素格式，它可以是从 uint8_t 到 float 的任何值。在这个用例中，它肯定是 floats.
我的简化代码如下所示：
#include <stdint.h>
static const int width = 1920;
static const int height = 1080;
void unpackFloatToUint16(float* pixels, uint16_t *bufferY, uint16_t *bufferU, uint16_t *bufferV, uint16_t *bufferA)
{
for (int r = height - 1; r >= 0; r--)
{
for (int c = 0; c < (int)width * 4; c += 4)
{
const int pos = r * width * 4 + c;
*bufferV++ = (uint16_t)((pixels[pos] + 0.57143f) * 57342.98164f);
*bufferU++ = (uint16_t)((pixels[pos + 1] + 0.57143f) * 57342.98164f);
*bufferY++ = (uint16_t)((pixels[pos + 2] + 0.07306f) * 56283.17216f);
*bufferA++ = (uint16_t)(pixels[pos + 3] * 65535.0f);
}
}
}
有一件事很明显：用乘法代替除法。如果您在 FP 划分上遇到瓶颈，吞吐量应该增加 ~7 到 ~15 倍左右。 （Haswell 的 divss 每 7 个时钟吞吐量一个，但 mulss 是每 0.5 个时钟一个）。
（假设您还没有使用 -ffast-math 来让编译器为您用倒数乘法代替常量除法）。
GCC 和 clang 已经自动向量化你的函数（至少使用编译时常量 height 和 width。问题中的代码没有编译，因为它们没有定义，所以我不知道该假设什么）See it on the Godbolt compiler explorer。在没有 -ffast-math 的情况下，它确实使用 divps 作为除法，但它确实使用 SIMD 进行数学运算（包括转换为 32 位整数），并通过洗牌将第 4 组 Y 值混合在一起对于 64 位商店。我不认为它的工作效率很高，但如果您在 div 吞吐量方面遇到瓶颈，那么它可能比 gcc 好得多。
但是对于 int height = rowBytes >> 1;，clang 不会自动矢量化，而 gcc 仍然可以。
不过，编译器的工作似乎还有改进的余地。
无论如何，假设我们要为 AVX + FMA（例如 Haswell 或 Steamroller / Ryzen）手动矢量化。您也可以制作其他版本，但由于您没有指定任何关于您想要定位的微体系结构（甚至是 x86），我将作为一个有趣的例子来做。
首先，我们可以将(Y + c0) / (c0 + c1) * 65535.0f转化为单个FMA。将加法里面的* (1.0f/(c0+c1)) * 65535.0f分布得到(Y * mul_const + add_const)，可以用单个FMA求值。我们可以一次对一个像素的所有 4 个分量执行此操作，使用 128 位 SIMD FMA 和两个向量常量以匹配内存中 float 布局的顺序保存系数。 （或者使用 256 位 FMA 一次处理所有 2 个像素）。
不幸的是，gcc 和 clang 没有使用 -ffast-math 进行此优化。
在转换为整数之前保存所有改组可能效果最好。它的优点是您只需要两个 FP 常量向量，而不是将每个分量的系数广播到所有元素的单独向量。好吧，我猜你可以在转换为整数之前对 FMA 的结果使用 FP 洗牌。 （例如FMA，然后洗牌，然后将4个Y值的向量转换为整数）。
FP 数学不是严格关联或分配的（因为舍入误差）所以这可能会改变边缘情况的结果。但不一定会使情况变得更糟，只是与您以前的舍入方式不同。如果使用 FMA 将 (Y + const1) * const2 转换为 Y * altconst1 + altconst2 不会丢失任何精度，因为 FMA 在进行加法之前不会舍入内部临时产品。
所以我们知道如何有效地进行数学运算并转换为整数（2 CPU 条指令用于包含 2 个像素的 8 float 向量）。剩下的就是将 Y 与其他 Y 组合在一起，并将 32 位有符号整数压缩为 16 位无符号整数。 （x86 只能在 FP 和有符号整数之间转换，直到 AVX512F 引入直接 FP <-> 无符号（和 SIMD 用于 64 位整数 <-> FP 而不是 64 位模式下的标量）。无论如何，我们不能直接从 float 转换为 16 位无符号整数向量）。
因此，给定一个包含 VUYA 个 32 位整数元素的 128 位向量，我们的第一步可能是缩小到 16 位整数。 x86 有一条指令（SSE4.1 packusdw，内在 _mm_packus_epi32）来打包无符号饱和（因此负输入饱和为 0，而大的正输入饱和为 65535）。大概这就是您想要的，而不是 t运行cating 会使溢出回绕的整数。这需要 2 个 SIMD 向量作为输入，并产生一个输出向量，所以我们会得到 VYUAVYUA（对于 2 个不同的像素）。
即使您不需要饱和行为（例如，如果输入不可能超出范围），packusdw 可能仍然是缩小整数范围的最有效选择。其他洗牌只有一个输入向量，或固定的洗牌模式，它不会丢弃每个 32 位元素的上半部分，因此在 pshufb 之后，结果中只有 64 位有用数据或 punpck 随机播放。
立即从 pack 开始可以很好地将事情减少到只有 2 个向量。我查看了其他洗牌顺序，但如果您从 32 位洗牌开始而不是压缩到 16 位元素，则总是需要更多的总洗牌。 （见下方神马link中的评论）
vpackusdw 的 256 位 AVX2 版本分别在 256 位向量的两个 128 位通道上运行，因此对于 pack(ABCD EFGH, IJKL MNOP) 你得到 ABCD IJKL EFGH MNOP。通常你需要另一次洗牌才能把事情按正确的顺序排列。无论如何，我们需要进一步洗牌，但这仍然很麻烦。不过，我认为您可以在每次循环迭代中处理两倍的数据，而循环中只需再洗几下。
这是我使用 128 位向量得出的结果
源代码 + 编译器输出 on the Godbolt compiler explorer
请注意，它不处理像素数不是 4 的倍数的情况。您可以使用清理循环（加载、缩放和填充饱和度，然后提取四个16 位组件）。或者你可以做一个部分重叠的最后 4 个像素。 （如果像素数确实是 4 的倍数，则不重叠，否则部分重叠存储到 Y、U、V 和 A 数组中。）这很容易，因为它不是就地操作，所以你可以重新读取相同的存储输出后输入。
此外，它假定行步幅与宽度相匹配，因为您在问题中的代码做了同样的事情。因此，宽度是否为 4 像素的倍数并不重要。但是，如果您确实有一个与宽度分开的可变行步幅，则您将不得不担心每行末尾的清理。 （或者使用填充，这样你就不必这样做了）。
#include <stdint.h>
#include <immintrin.h>
static const int height = 1024;
static const int width = 1024;
// helper function for unrolling
static inline __m128i load_and_scale(const float *src)
{ // and convert to 32-bit integer with truncation towards zero.
// Scaling factors (note min. values are actually negative) (limited range)
const float yuvaf[4][2] = {
{ 0.07306f, 1.09132f }, // Y
{ 0.57143f, 0.57143f }, // U
{ 0.57143f, 0.57143f }, // V
{ 0.00000f, 1.00000f } // A
};
// (Y + yuvaf[n][0]) / (yuvaf[n][0] + yuvaf[n][1]) ->
// Y * 1.0f/(yuvaf[n][0] + yuvaf[n][1]) + yuvaf[n][0]/(yuvaf[n][0] + yuvaf[n][1])
// Pixels are in VUYA order in memory, from low to high address
const __m128 scale_mul = _mm_setr_ps(
65535.0f / (yuvaf[2][0] + yuvaf[2][1]), // V
65535.0f / (yuvaf[1][0] + yuvaf[1][1]), // U
65535.0f / (yuvaf[0][0] + yuvaf[0][1]), // Y
65535.0f / (yuvaf[3][0] + yuvaf[3][1]) // A
);
const __m128 scale_add = _mm_setr_ps(
65535.0f * yuvaf[2][0] / (yuvaf[2][0] + yuvaf[2][1]), // V
65535.0f * yuvaf[1][0] / (yuvaf[1][0] + yuvaf[1][1]), // U
65535.0f * yuvaf[0][0] / (yuvaf[0][0] + yuvaf[0][1]), // Y
65535.0f * yuvaf[3][0] / (yuvaf[3][0] + yuvaf[3][1]) // A
);
// prefer having src aligned for performance, but with AVX it won't help the compiler much to know it's aligned.
// So just use an unaligned load intrinsic
__m128 srcv = _mm_loadu_ps(src);
__m128 scaled = _mm_fmadd_ps(srcv, scale_mul, scale_add);
__m128i vuya = _mm_cvttps_epi32(scaled); // truncate toward zero
// for round-to-nearest, use cvtps_epi32 instead
return vuya;
}
void deinterleave_avx_fma(char* __restrict pixels, int rowBytes, char *__restrict bufferY, char *__restrict bufferU, char *__restrict bufferV, char *__restrict bufferA)
{
const float *src = (float*)pixels;
uint16_t *__restrict Y = (uint16_t*)bufferY;
uint16_t *__restrict U = (uint16_t*)bufferU;
uint16_t *__restrict V = (uint16_t*)bufferV;
uint16_t *__restrict A = (uint16_t*)bufferA;
// 4 pixels per loop iteration, loading 4x 16 bytes of floats
// and storing 4x 8 bytes of uint16_t.
for (unsigned pos = 0 ; pos < width*height * 4; pos += 4) {
// pos*4 because each pixel is 4 floats long
__m128i vuya0 = load_and_scale(src+pos*4);
__m128i vuya1 = load_and_scale(src+pos*4 + 1);
__m128i vuya2 = load_and_scale(src+pos*4 + 2);
__m128i vuya3 = load_and_scale(src+pos*4 + 3);
__m128i vuya02 = _mm_packus_epi32(vuya0, vuya2); // vuya0 | vuya2
__m128i vuya13 = _mm_packus_epi32(vuya1, vuya3); // vuya1 | vuya3
__m128i vvuuyyaa01 = _mm_unpacklo_epi16(vuya02, vuya13); // V0V1 U0U1 | Y0Y1 A0A1
__m128i vvuuyyaa23 = _mm_unpackhi_epi16(vuya02, vuya13); // V2V3 U2U3 | Y2Y3 A2A3
__m128i vvvvuuuu = _mm_unpacklo_epi32(vvuuyyaa01, vvuuyyaa23); // v0v1v2v3 | u0u1u2u3
__m128i yyyyaaaa = _mm_unpackhi_epi32(vvuuyyaa01, vvuuyyaa23);
// we have 2 vectors holding our four 64-bit results (or four 16-bit elements each)
// We can most efficiently store with VMOVQ and VMOVHPS, even though MOVHPS is "for" FP vectors
// Further shuffling of another 4 pixels to get 128b vectors wouldn't be a win:
// MOVHPS is a pure store on Intel CPUs, no shuffle uops.
// And we have more shuffles than stores already.
//_mm_storeu_si64(V+pos, vvvvuuuu); // clang doesn't have this (AVX512?) intrinsic
_mm_storel_epi64((__m128i*)(V+pos), vvvvuuuu); // MOVQ
_mm_storeh_pi((__m64*)(U+pos), _mm_castsi128_ps(vvvvuuuu)); // MOVHPS
_mm_storel_epi64((__m128i*)(Y+pos), yyyyaaaa);
_mm_storeh_pi((__m64*)(A+pos), _mm_castsi128_ps(yyyyaaaa));
}
}
希望 shuffle 的变量名 + 注释应该是相当可读的。这是未经测试的；最可能的错误是某些向量作为随机播放参数的错误排序。但是修复它应该只是颠倒 arg 顺序之类的问题，而不需要额外的洗牌来减慢它的速度。
看起来包括包装在内的 6 次洗牌是我能做的最好的了。它基本上包括一个从 vuya x4 到 vvvv uuuu yyyy aaaa 的 4x4 转置，并且 pack 有一个固定的洗牌模式，这对去交错没有帮助，所以我认为我们不能做得比这更好128 位向量。当然，我总是有可能忽略了一些东西。
gcc 和 clang 都编译它略微次优:
Clang 使用 vpextrq 而不是 vmovhps（在 Intel CPUs 上，每次循环迭代总共花费额外的 2 个 shuffle uops）。并且还使用 2 个单独的循环计数器，而不是将同一个计数器缩放 1 或 8，这样会花费 1 个额外的整数 add 指令，没有任何好处。 （要是 gcc 选择这样做而不是使用折叠到 FMA 中的索引加载就好了……愚蠢的编译器。）
gcc，而不是使用 vmovups 负载，处理 FMA3 通过复制向量常量然后使用具有索引寻址模式的内存操作数来破坏其输入。 This doesn't stay micro-fused，所以前端总共有 4 个额外的微指令。
如果编译完美，仅使用一个循环计数器作为 float 源数组索引（按 *8 缩放）和整数目标数组（未缩放），就像 gcc 那样，并且如果 gcc 以 clang 的方式加载，那么整个循环将是 Haswell/Skylake.
因此它可以每 6 个时钟从前端发出一次迭代，恰好在每个时钟 4 个融合域微指令的限制下。它包含 6 个 shuffle uops，因此它也可以解决 port5 吞吐量瓶颈。 （HSW / SKL 只有 1 个洗牌执行单元）。所以 理论上 ，这可以 运行 每 6 个时钟 4 个像素（16 个浮点数），或者在英特尔 CPUs 上每 1.5 个时钟周期一个像素。在 Ryzen 上可能稍微多一些，虽然 MOVHPS 在那里花费了多个 uops，但管道更宽。 （参见 http://agner.org/optimize/）
每 6 个时钟 4 次加载和 4 次存储很好，并且远离任何瓶颈，但如果您的 src 和 dst 在缓存中不热，则可能存在内存带宽。商店只有一半宽度，一切都是顺序的。来自一个循环的 4 个独立输出流已经足够少，通常不会成为问题。 （如果你有超过 4 个输出流，你可能会考虑裂变循环并只存储一次传递的 U 和 V 值，然后在另一次传递中只存储 Y 和 A 值传递相同的源数据，或者类似的东西。但是就像我说，4 个输出流就可以了，不保证循环裂变。）
这个的 256 位版本最后可能需要的不仅仅是 2 个额外的 vpermq，因为我认为你不能轻易解决你想要相邻的两个 V 值 [= =49=] 卡在同一个 __m256i 向量的高低车道。因此，您可能需要在此过程的早期进行额外的 4 vpermd 或 vperm2i128 洗牌，因为最小粒度的跨车道洗牌是 32 位粒度。这可能会对 Ryzen 造成很大伤害。
在将 4 或 8 个 v 组合在一起后，也许您可以使用 vpblendw 重新排列向量之间的词元素，但 v 是错误的。
vpackusdw 慢，所以可能不会。