题目:Vector instructions. Part II: Vectorization 作者:Oleg Ponomarev 文章地址:https://www.elecard.com/page/article_vector_instructions_part2 内容整理:刘潮磊 本文是对 ELECARD Video Compression Book 第十三章的翻译,本章节介绍了主要的矢量化技术,并举例说明了几种简单的视频信号编码/解码算法。
本章节所有示例都将使用某个图像的像素块作为输入数据。为简单起见,考虑一个像素值范围为
的单色图像,其中
为像素位深度。图像由一维数组表示。在示例中,位深度要么为8位,要么为9 ~ 16位。因此,图像被存储为字节数组或无符号16位整数变量。代码特意简化了:所有示例中的图像块都是正方形的,大小分别为4、8或16像素。每个函数都使用特定大小和像素位深度的块。所有示例中使用的标识符如下:
最简单的图像处理问题之一是复制图像的一部分。下面的例子是一个复制8位图像块的函数,作为c++模板实现。
template <int SIZE>
void copy_mb(const uint8_t * src,
uint8_t * dst,
size_t src_stride,
size_t dst_stride)
{
for(int i = 0; i < SIZE; i++)
{
for(int j = 0; j < SIZE; j++)
{
dst[j] = src[j];
}
src += src_stride;
dst += dst_stride;
}
}
复制是在一个内部循环中完成的,这个内部循环可以被替换,例如,通过调用标准memset函数来替换。但是,我们不调用这个函数,而是使用向量指令,将数据从RAM读入寄存器,然后将结果写回RAM。在SSE2指令集中有三个指令对用于复制数据块:
_mm_cvtsi32_si128
和_mm_cvtsi128_si32
用于4字节块_mm_loadl_epi64
和_mm_store_epi64
用于8字节块,_mm_loadu_si128
和_mm_storeu_si128
用于16字节块。
在后一种情况下,也可以使用_mm_load_si128
和_mm_store_si128
,前提是读写地址与16字节边界对齐——即16字节的倍数。
由此,我们得到了针对不同块大小的复制函数的三种实现。
#include <emmintrin.h>
void copy_mb_4(const uint8_t * src,
uint8_t * dst,
size_t src_stride,
size_t dst_stride)
{
__m128i x0;
for(int i = 0; i < 4; i++)
{
x0 = _mm_cvtsi32_si128(*(int32_t*) src);
*(int32_t*) dst = _mm_cvtsi128_si32(x0);
src += src_stride;
dst += dst_stride;
}
}// copy_mb_4
void copy_mb_8(const uint8_t * src,
uint8_t * dst,
size_t src_stride,
size_t dst_stride)
{
__m128i x0;
for(int i = 0; i < 8; i++)
{
x0 = _mm_loadl_epi64((__m128i*)src);
_mm_storel_epi64((_m128i*)dst, x0);
src += src_stride;
dst += dst_stride;
}
}// copy_mb_8
void copy_mb_16(const uint8_t * src,
uint8_t * dst,
size_t src_stride,
size_t dst_stride)
{
__m128i x0;
for(int i = 0; i < 16; i++)
{
x0 = _mm_loadu_si128((__m128i*)src);
_mm_storeu_si128((_m128i*)dst, x0);
src += src_stride;
dst += dst_stride;
}
}// copy_mb_16
也可以去掉外部循环,如下面的示例3所示。这样的“循环展开”是一种常用的技术,在迭代次数较少且事先已知且循环体相对较小的情况下非常有用。
#include <emmintrin.h>
void copy_mb_4(const uint8_t * src,
uint8_t * dst,
size_t src_stride,
size_t dst_stride)
{
__m128i x0, x1, x2, x3;
x0 = _mm_cvtsi32_si128(*(int32_t*)(src + 0 * src_stride));
x1 = _mm_cvtsi32_si128(*(int32_t*)(src + 1 * src_stride));
x2 = _mm_cvtsi32_si128(*(int32_t*)(src + 2 * src_stride));
x3 = _mm_cvtsi32_si128(*(int32_t*)(src + 3 * src_stride));
*(int32_t*)(dst + 0 * dst_stride) = _mm_cvtsi128_si32(x0);
*(int32_t*)(dst + 1 * dst_stride) = _mm_cvtsi128_si32(x1);
*(int32_t*)(dst + 2 * dst_stride) = _mm_cvtsi128_si32(x2);
*(int32_t*)(dst + 3 * dst_stride) = _mm_cvtsi128_si32(x3);
}
复制9位到16位块的函数的实现与8位的情况类似。所需要的只是一条复制大量字节的指令(例如,_mm_loadl_epi64
而不是_mm_cvtsi32_si128)
,或者重复几次相同的指令(最适合16x16像素块)。
考虑从两个块中添加像素。这是视频解码的一个组成部分,其中一个像素块(在本例中为dst)通过插值当前或以前帧的像素来计算,另一个块(src)通过对解码系数应用反向离散余弦变换来计算。示例省略了 src 和 dst 块的计算,只显示了最后的补偿阶段,其中这些块的像素值被加在一起,结果被写入 dst 块。
在最简单的情况下(相同的像素块变量类型,没有整数溢出),补偿实现如下。
template <int SIZE>
void compensate_mb(const uint16_t * src,
uint16_t * dst,
size_t src_stride,
size_t dst_stride)
{
for(int i = 0; i < SIZE; i++)
{
for(int j = 0; j < SIZE; j++)
{
dst[j] = dst[j] + src[j];
}
src += src_stride;
dst += dst_stride;
}
}
我们已经在上面讨论了块复制指令。现在,我们将用一个附加指令(在本例中是_mm_add_epi16
)来补充它们。得到的向量化代码如下所示。(本例中使用了一个8x8的块)。
void compensate_8(const uint16_t * src
uint16_t * dst,
size_t src_stride,
size_t dst_stride)
{
__m128i x0, x1;
for(int i = 0; i < 8; i++)
{
x0 = _mm_loadu_si128((__m128i*)src); // 8 pixels
x1 = _mm_loadu_si128((__m128i*)dst);
x0 = _mm_add_epi16( x0, x1);
_mm_storeu_si128((_m128i*)dst, x0);
src += src_stride;
dst += dst_stride;
}
}
在实际的视频解码中,反向离散余弦变换(src)输出的块大小将超过插值块(dst)的大小。此外,src块元素可以取负值。更实际的补偿实现如下所示。
template <int SIZE>
void compensate_mb(const int16_t * src,
uint8_t * dst,
size_t src_stride,
size_t dst_stride)
{
for(int i = 0; i < SIZE; i++)
{
for(int j = 0; j < SIZE; j++)
{
int tmp = dst[j] + src[j];
if(tmp > MAX(uint8_t))
dst[j] = MAX(uint8_t);
else if (tmp < 0)
dst[j] = 0;
else
dst[j] = tmp;
}
src += src_stride;
dst += dst_stride;
}
}
在这里,tmp变量的值不应超过dst像素值的允许范围(在本例中为0…255),因为它受到上限或下限的限制。
与示例5相比,向量化实现变得更加复杂,因为dst值必须从8位转换为有符号或无符号的16位值。还需要实现反向转换为带有上限和下限的无符号8位值。将无符号整数的大小加倍的最简单方法是使用_mm_unpacklo_epiX
和_mm_unpackhi_epiX
指令,其中X = 8、16、32或64。例如,一个由8位无符号整数组成的向量可以转换为两个16位整数的向量,如下所示:
0 = _mm_setzero_si128();
X1 = x0;
X0 = _mm_unpacklo_epi8(X0, 0);
X1 = _mm_unpackhi_epi8(X1, 0);
较大的数据项以类似的方式进行转换。
要将16位值转换为无符号8位值,可以使用_mm_packus_epi16
指令,该指令将两个矢量寄存器的内容打包为一个。对于任何值在0之外的16位项…255范围内,它也将值截断到该范围。因此,8х8块的矢量化实现将类似于例7。
void compensate_8(const int16_t * src,
uint8_t * dst,
size_t src_stride,
size_t dst_stride)
{
__m128i x0, x1, zero;
zero = _mm_setzero_si128();
for(int i = 0; i < 8; i++)
{
x0 = _mm_loadu_si128((__m128i*)src); // 8 pixels
x1 = _mm_loadl_epi64((__m128i*)dst); // 8 bit !
x1 = _mm_unpacklo_epi8(x1, zero); // from 8 to 16 bit
x0 = _mm_add_epi16( x0, x1);
x0 = _mm_packus_epi16(x0, x0); // back to 8 bit
_mm_storel_epi64((_m128i*)dst, x0);
src += src_stride;
dst += dst_stride;
}
}
现在考虑这样一个情况,dst的像素位深度在9到16位之间,并且使用反向离散余弦变换计算的src块由32位值组成。将32位值转换为
范围并不那么简单,因为对于任意数据大小,没有类似于_mm_packus_epi16的指令。转换分两步完成。首先,使用_mm_packs_epi32
指令将两个32位无符号整数向量打包成一个16位有符号整数向量(在-32,768..32,767范围内)。接下来,使用_mm_min_epi16
和_mm_max_epi16
指令通过比较两个寄存器的相应元素并分别取最小值和最大值来截断超出范围的值。(注意,如果使用来自SSE4.1集的_mm_packus_epi32
指令将32位值转换为无符号16位值,那么进行一次比较就足够了。)一个4x4块的完整代码类似于示例8(假设32位加法没有溢出)。
void compensate_4(const int32_t * src,
uint16_t *dst,
size_t src_stride,
size_t dst_stride,
int bitdepth)
{
__m128i x0, x1, zero, max_val;
zero = _mm_setzero_si128();
max_val = _mm_set1_epi16((1 << bitdepth) — 1);
for(int i = 0; i < 4; i++)
{
x0 = _mm_loadu_si128((__m128i*)src); // 4 x 32
x1 = _mm_loadl_epi64((__m128i*)dst); // 4 x 16
x1 = _mm_unpacklo_epi16(x1, zero); // from 16 bit to 32 bit
x0 = _mm_add_epi32(x0, x1);
x0 = _mm_packs_epi32( x0, x0 ); // from 32 bit to 16 bit
/* if x0[k] < max_val, then x0[k]. else max_val */
x0 = _mm_min_epi16(x0, max_val);
x0 = _mm_max_epi16(x0, zero);
_mm_storel_epi64((_m128i*)dst, x0);
src += src_stride;
dst += dst_stride;
}
}
该函数使用_mm_set1_epi16
将向量寄存器的16位元素设置为相同的值。实际上,_mm_set1_epi16
是一个伪指令(就像其他类似的指令一样),具有依赖于编译器的实现。当需要高性能时,应该避免使用此类伪指令。
两幅图像之间的差异程度是使用性能指标来确定的,性能指标是包含两幅图像像素值的表达式。性能指标可以用来确定最佳的编码方法。
考虑两个性能指标的计算,绝对差和(SAD)和均方误差(MSE)。对于两个大小相同的图像,这些指标由以下公式定义:
其中
和
分别为图像的高度和宽度。
在“纯”C/С++中计算这些表达式的代码非常明显,因此这里没有显示。相反,让我们看看如何使用向量指令来实现这段代码。
有一条特殊的SSE2指令_mm_sad_epu8
,用于计算两个8位数据块之间的绝对差之和。它分别为寄存器的最小有效部分和最高有效部分计算SAD,并以相同的方式存储总和。由于两个8位值的绝对差不能超过255,所以每次和不能超过2040。
对于16x16像素的块,SAD的计算方法如下:
#include <emmintrin.h>
#include <stdint.h>
int32_t sad_16_8bit(const uint8_t* src0,
const uint8_t* src1,
size_t src0_stride,
size_t src1_stride)
{
__m128i x0, x1, sum;
sum = _mm_setzero_si128();
for(int i = 0; i < 16; i++)
{
x0 = _mm_loadu_si128((__m128i*)src0);
x1 = _mm_loadu_si128((__m128i*)src1);
x0 = _mm_sad_epu8(x0, x1);
sum = _mm_add_epi32(sum, x0); // sum for lower and upper halves
src0 += src0_stride;
src1 += src1_stride;
}
x0 = _mm_shuffle_epi32(sum, _MM_SHUFFLE(1,0,3,2));
sum = _mm_add_epi32(sum, x0);
int32_t s = _mm_cvtsi128_si32(sum); // result
return s;
}
由于_mm_sad_epu8
对块的左右两部分分别计算两个单独的和,所以sum寄存器也会累积两个和。因此,最终结果等于sum寄存器中所有32位元素的和(更准确地说,是第一个和第三个元素)。因此,我们使用_mm_shuffle_epi32
交换sum寄存器的最低有效位和最高有效位,并将结果存储在x0寄存器中。通过将x0和sum寄存器加在一起,我们得到了后者的最低有效位32位的最终结果。
对于较小的块,计算是类似的,除了求和寄存器的最有效的一半是零,并且不需要额外的加法。这是因为_mm_cvtsi32_si128
和_mm_loadl_epi64
分别用零填充寄存器的最高位96和64位。
目前,对于大于8位的数据大小,没有类似于_mm_sad_epu8
的指令。对于9到15位的输入数据,计算绝对差的和很简单:
X0 = _mm_sub_epi16(X0, x1);// x0 - x1
X0 = _mm_abs_epi16(X0);/ / SSE3
这种方法要求CPU支持SSE3指令集,并且不适合输入刚好16位的数据(在实践中很少遇到)。绝对差值可以计算如下:
x2 = x0;
x0 = _mm_subs_epu16(x0, x1); // x0 — x1 or 0
x1 = _mm_subs_epu16(x1, x2); // x1 — x2 or 0
x0 = _mm_xor_si128(x0, x1); // | x0 - x1 | or 0
如果减号大于减号,则_mm_subs_epu16
指令(无符号饱和减法)的结果为零。在这种情况下,_mm_subs_epu16(x0,x1)
和_mm_subs_epu16(x1,x2)
产生x0和x1中所有元素对的绝对差值或零。剩下的就是使用OR或XOR来组合这些元素。通过将16位指令替换为8位指令(即_mm_subs_epu8
),可以类似地计算8位数据的绝对差异。
下面是为每像素16位的8x8像素图像块计算SAD的示例。
__m128i x0, x1, x2, sum, zero;
zero = _mm_setzero_si128();
sum = zero;
for(int i = 0; i < 8; i++)
{
x0 = _mm_loadu_si128((__m128i*)src0);
x1 = _mm_loadu_si128((__m128i*)src1);
/* | x0 — x1 | */
x2 = x0;
x0 = _mm_subs_epu16(x0, x1);
x1 = _mm_subs_epu16(x1, x2);
x0 = _mm_xor_si128(x0, x1);
x1 = x0;
x0 = _mm_unpacklo_epi16(x0, zero); // 16 bit to 32 bit
x1 = _mm_unpackhi_epi16(x1, zero);
sum = _mm_add_epi32(sum, x0);
sum = _mm_add_epi32(sum, x1);
src0 += src0_stride;
src1 += src1_stride;
}
/* sum is a0,a1,a2,a3 */
x0 = _mm_shuffle_epi32(sum, _MM_SHUFFLE(2,3,0,1)); // x0 is a1,a0,a3,a2
sum = _mm_add_epi32(sum, x0);
x0 = _mm_shuffle_epi32(sum, _MM_SHUFFLE(1,0,3,2));
sum = _mm_add_epi32(sum, x0);
int32_t s = _mm_cvtsi128_si32(sum); // result
这里,像素值之间的绝对差异计算如前面所示。然而,它们不能以这种(16位)形式求和,因为可能会溢出。因此,绝对值被转换为32位,然后求和。循环完成后,剩下的唯一工作就是计算sum寄存器中所有32位元素的和。作为上述示例中使用的替代方法,可以使用_mm_hadd_epi32
指令完成此操作,该指令对相邻的32位寄存器元素求和。这需要CPU支持SSSE3。
zero = _mm_setzero_si128();
sum = _mm_hadd_epi32(sum, zero); // sum is a0+a1,a2+a3,0,0
sum = _mm_hadd_epi32(sum, zero); // sum is a0+a1+a2+a3,0,0,0
严格来说,下面所有例子中计算的数量都不是MSE,而是:
由于涉及到平方,因此必须考虑到8位的潜在溢出,对于更大的数据大小更是如此。这意味着需要进行额外的数据大小转换。
下面显示了一个为8位数据块计算SED的示例。
__m128i x0, x1, x2, sum, zero;
zero = _mm_setzero_si128();
sum = zero;
for(int i = 0; i < 8; i++)
{
x0 = _mm_loadl_epi64((__m128i*)src0);
x1 = _mm_loadl_epi64((__m128i*)src1);
x0 = _mm_unpacklo_epi8(x0, zero); // 8 to 16 bit
x1 = _mm_unpacklo_epi8(x1, zero);
x0 = _mm_sub_epi16(x0, x1); // x0 — x1
x0 = _mm_madd_epi16(x0, x0); // (x0 - x1)^2
sum = _mm_add_epi32(sum, x0);
src0 += src0_stride;
src1 += src1_stride;
}
// sum of sum elements
在本例中,首先将8位数据转换为16位,然后计算像素值差。获得平方差的最方便的方法是_mm_madd_epi16
指令,该指令将16位数据直接转换为32位,并执行一些所需的加法。当循环完成时,只需要将sum寄存器中所有元素的值相加,如例10所示。
SED可以以类似的方式计算大小不超过12位的数据和16x16像素的块。对于较大的数据,在循环中求和时需要将32位转换为64位。此外,由于可能溢出,_mm_madd_epi16
不能用于16位数据。应该使用_mm_mullo_epi16
和_mm_mulhi_epu16
指令。这些将在下面的示例12中使用。此外,计算两个像素值之间的绝对差值而不是简单差值,以避免不必要的数据大小转换。
__m128i x0, x1, x2, sum, zero;
zero = _mm_setzero_si128();
sum = zero;
for(int i = 0; i < 8; i++)
{
x0 = _mm_loadu_si128((__m128i*)src0);
x1 = _mm_loadu_si128((__m128i*)src1);
/* | x0 — x1 | */
x2 = x0;
x0 = _mm_subs_epu16(x0, x1);
x1 = _mm_subs_epu16(x1, x2);
x0 = _mm_xor_si128(x0, x1);
/* x0^2 */
x1 = x0;
x0 = _mm_mullo_epi16( x0, x0 );
x1 = _mm_mulhi_epu16( x1, x1 );
x2 = x0;
x0 = _mm_unpacklo_epi16( x0, x1 ); // x0[i]^2, i = 0..3
x2 = _mm_unpackhi_epi16( x2, x1 ); // x0[i]^2, i = 4..7
x0 = _mm_add_epi32( x0, x2 );
x2 = x0;
x0 = _mm_unpacklo_epi32(x0, zero); // from 32 to 64 bit
x2 = _mm_unpackhi_epi32(x0, zero);
sum = _mm_add_epi64(sum, x0);
sum = _mm_add_epi64(sum, x2);
src0 += src0_stride;
src1 += src1_stride;
}
// sum of sum elements
x0 = _mm_shuffle_epi32(sum, _MM_SHUFFLE(1,0,3,2));
sum = _mm_add_epi64(sum, x0);
uint64_t result;
_mm_storel_epi64((__m128i*)&result, sum);
在前面的示例中,图像块的大小固定为4、8或16像素,但是矢量化也可以直接用于任意块大小。让我们回顾一下例子1和6。在这两种方法中,所有的数据处理和复制都是在内部循环中完成的,而外部循环只是将指针移动到下一个像素行。这里,内循环只包含一次迭代,因此省略了循环操作符。然而,对于任意宽度的块,需要多次迭代来覆盖宽度和相同数量的指针移动。在一般情况下,内部循环看起来与示例13相同,其中为8位图像计算SAD。
#include <emmintrin.h>
#include <stdint.h>
#include <stdlib.h>
uint64_t sad_8bit(const uint8_t* src0,
const uint8_t* src1,
size_t width,
size_t height,
size_t src0_stride,
size_t src1_stride)
{
size_t width16 = width — (width % 16); // width16 == 16*x
__m128i x0, x1, sum;
sum = _mm_setzero_si128();
uint64_t sum_tail = 0;
for(int i = 0; i < height; i++)
{
for(int j = 0; j < width16; j += 16)
{
x0 = _mm_loadu_si128((__m128i*)(src0 + j));
x1 = _mm_loadu_si128((__m128i*)(src1 + j));
x0 = _mm_sad_epu8(x0, x1);
sum = _mm_add_epi64(sum, x0);
}
for(int j = width16; j < width; j ++)
{
sum_tail += abs(src0[j] — src1[j]);
}
src0 += src0_stride;
src1 += src1_stride;
}
x0 = _mm_shuffle_epi32(sum, _MM_SHUFFLE(1,0,3,2));
sum = _mm_add_epi64(sum, x0);
uint64_t sum_total;
_mm_storel_epi64((__m128i*)&sum_total, sum);
sum_total += sum_tail;
return sum_total;
}
行长度通常不会是4、8或16字节的倍数。因此,产生的“尾巴”将被单独处理。没有指令从RAM中加载任意数量的字节。最简单的解决方法是不向量化这段代码,如例9和13所示。在下面的示例中,“尾部”永远不会超过15字节,并且对于足够大的图像块的性能损失很小。然而,如果要实现的算法需要大量的CPU时间,则希望使用未向量化的代码处理尽可能少的数据。在这种情况下,可以使用例14所示的技术。
#include <emmintrin.h>
#include <stdint.h>
#include <stdlib.h>
uint64_t sad_8bit(const uint8_t* src0,
const uint8_t* src1,
size_t width,
size_t height,
size_t src0_stride,
size_t src1_stride)
{
size_t width_r = width % 16;
size_t width16 = width — width_r; // width16 == 16*x
size_t width8 = width_r — (width_r % 8); // 8 or 0
width_r -= width8;
size_t width4 = width_r — (width_r % 4); // 4 or 0
width_r -= width4; // 0, 1, 2, or 3
__m128i x0, x1, sum;
sum = _mm_setzero_si128();
uint64_t sum_tail = 0;
for(int i = 0; i < height; i++)
{
for(int j = 0; j < width16; j += 16)
{
/* SAD calculation */
}
if( width8)
{
x0 = _mm_loadl_epi64((__m128i*)(src0 + width16));
x1 = _mm_loadl_epi64((__m128i*)(src1 + width16));
x0 = _mm_sad_epu8(x0, x1);
sum = _mm_add_epi64(sum, x0);
}
if( width4)
{
x0 = _mm_cvtsi32_si128(*(int32_t*)(src0 + width16 + width8));
x1 = _mm_cvtsi32_si128(*(int32_t*)(src1 + width16 + width8));
x0 = _mm_sad_epu8(x0, x1);
sum = _mm_add_epi64(sum, x0);
}
for(int j = width - width_r; j < width; j ++)
{
sum_tail += abs(src0[j] — src1[j]);
}
src0 += src0_stride;
src1 += src1_stride;
}
/**/
}
在本例中,未经向量化处理的像素数永远不会超过3个。