使用通用内联函数向量化你的代码
本教程的目标是提供一份指南,介绍如何使用通用内联函数(Universal Intrinsics)特性来向量化你的 C++ 代码,以获得更快的运行速度。我们将简要了解 SIMD 内联函数(SIMD intrinsics)以及如何使用宽寄存器(register),随后通过教程学习使用宽寄存器进行基本操作。
在本节中,我们将简要了解几个概念,以帮助更好地理解其功能。
内联函数(Intrinsics)
Section titled “内联函数(Intrinsics)”内联函数是由编译器单独处理的函数。这些函数通常经过优化,能够以尽可能高效的方式执行,因此比普通实现运行得更快。然而,由于这些函数依赖于编译器,这使得编写可移植的应用程序变得困难。
SIMD 代表单指令多数据(Single Instruction, Multiple Data)。SIMD 内联函数允许处理器对计算进行向量化。数据存储在所谓的寄存器中。一个寄存器可以是 128 位、256 位或 512 位宽。每个寄存器存储相同数据类型的多个值。寄存器的宽度和每个值的大小共同决定了总共能存储的值的数量。
根据你的 CPU 支持的指令集,你可能能够使用不同的寄存器。要了解更多,请参见此处。
VLA 代表向量长度无关(Vector Length Agnostic)。这是一种寄存器宽度由硬件在运行时决定、而不是在编译时固定的机制。这使得单个二进制文件能够在同一架构(例如 RVV 或 SVE)内的不同 CPU 上扩展其性能。
通用内联函数(Universal Intrinsics)
Section titled “通用内联函数(Universal Intrinsics)”OpenCV 的通用内联函数为 SIMD 和 VLA 向量化方法提供了抽象,让用户无需编写系统特定的代码即可使用内联函数。
现在我们将介绍可用的结构和函数:
- 寄存器结构
- 加载与存储
- 数学运算
- 归约与掩码
通用内联函数集将每个寄存器实现为一个基于特定 SIMD 寄存器的结构。所有类型都包含 nlanes 枚举,它给出了该类型可以容纳的值的确切数量。这消除了在实现过程中硬编码值数量的需要。
注: 每个寄存器结构都在
cv命名空间下。
寄存器有两种类型:
-
可变大小寄存器(Variable sized registers):这些结构没有固定大小,其确切位长会在编译期间根据可用的 SIMD 能力推导得出。因此,
nlanes枚举的值在编译时确定。每个结构遵循以下命名约定:
v_[type of value][size of each value in bits]例如,
v_uint8持有 8 位无符号整数,v_float32持有 32 位浮点值。我们像在 C++ 中声明任何对象一样声明一个寄存器。根据可用的 SIMD 指令集,一个特定的寄存器将容纳不同数量的值。例如:如果你的计算机最多支持 256 位寄存器:
-
v_uint8将容纳 32 个 8 位无符号整数 -
v_float64将容纳 4 个 64 位浮点数(double)v_uint8 a; // a is a register supporting uint8(char) dataint n = a.nlanes; // n holds 32可用的数据类型和大小:
类型 位宽 uint 8, 16, 32, 64 int 8, 16, 32, 64 float 32, 64
-
-
固定大小寄存器(Constant sized registers):这些结构具有固定的位宽,并容纳固定数量的值。我们需要知道系统支持什么 SIMD 指令集,并选择兼容的寄存器。仅在需要确切位长时才使用它们。
每个结构遵循以下命名约定:
v_[type of value][size of each value in bits]x[number of values]假设我们想存储:
-
一个 128 位寄存器中的 32 位(以位为单位的大小)有符号整数。由于寄存器大小已知,我们可以求出寄存器中的数据点数量(128/32 = 4):
v_int32x4 reg1 // holds 4 32-bit signed integers. -
512 位寄存器中的 64 位浮点数:
v_float64x8 reg2 // reg2.nlanes = 8
-
加载与存储操作
Section titled “加载与存储操作”现在我们知道了寄存器的工作方式,让我们看看用于向这些寄存器填充值的函数。
-
加载(Load):加载函数允许你把值加载到一个寄存器中。
-
构造函数——在声明寄存器结构时,我们既可以提供一个内存地址(寄存器将从那里读取连续的值),也可以显式地以多个参数的形式提供值(显式多参数形式仅适用于固定大小寄存器):
float ptr[32] = {1, 2, 3 ..., 32}; // ptr is a pointer to a contiguous memory block of 32 floats// Variable Sized Registers //int x = v_float32().nlanes; // set x as the number of values the register can holdv_float32 reg1(ptr); // reg1 stores first x values according to the maximum register size available.v_float32 reg2(ptr + x); // reg stores the next x values// Constant Sized Registers //v_float32x4 reg1(ptr); // reg1 stores the first 4 floats (1, 2, 3, 4)v_float32x4 reg2(ptr + 4); // reg2 stores the next 4 floats (5, 6, 7, 8)// Or we can explicitly write down the values.v_float32x4(1, 2, 3, 4); -
加载函数——我们可以使用 load 方法并提供数据的内存地址:
float ptr[32] = {1, 2, 3, ..., 32};v_float32 reg_var;reg_var = vx_load(ptr); // loads values from ptr[0] upto ptr[reg_var.nlanes - 1]v_float32x4 reg_128;reg_128 = v_load(ptr); // loads values from ptr[0] upto ptr[3]v_float32x8 reg_256;reg_256 = v256_load(ptr); // loads values from ptr[0] upto ptr[7]v_float32x16 reg_512;reg_512 = v512_load(ptr); // loads values from ptr[0] upto ptr[15]注: 加载函数假定数据是未对齐的。如果你的数据是对齐的,可以使用
vx_load_aligned()函数。
-
-
存储(Store):存储函数允许你把寄存器中的值存储到特定的内存位置。
-
要将值从寄存器存储到内存位置,可以使用 v_store() 函数:
float ptr[4];v_store(ptr, reg); // store the first 128 bits(interpreted as 4x32-bit floats) of reg into ptr.
-
注: 确保 ptr 与寄存器具有相同的类型。你也可以在执行操作之前将寄存器转换为正确的类型。简单地把指针强制转换为特定类型会导致对数据的错误解释。
二元和一元运算符
Section titled “二元和一元运算符”通用内联函数集提供了逐元素的二元和一元运算。
注: 自 OpenCV 4.11 起,通用内联函数中的 C++ 运算符重载(例如
+、*)已被弃用,取而代之的是显式的包装函数(例如v_add、v_mul),以确保与 VLA 架构的兼容性。另见:https://github.com/opencv/opencv/issues/27267
-
算术运算:我们可以对两个寄存器进行逐元素的加、减、乘、除。寄存器必须具有相同的宽度并持有相同的类型。例如,将两个寄存器相加、相乘:
v_float32 a, b; // {a1, ..., an}, {b1, ..., bn}v_float32 c = v_add(a, b); // {a1 + b1, ..., an + bn}v_flaot32 d = v_mul(a, b); // {a1 * b1, ..., an * bn} -
按位逻辑与移位:我们可以对寄存器中每个元素的位进行左移或右移。我们还可以在两个寄存器之间逐元素地应用按位与、或、异或和非运算符:
v_int32 as; // {a1, ..., an}v_int32 al = v_shl(as, 2); // {a1 << 2, ..., an << 2}v_int32 bl = v_shr(as, 2); // {a1 >> 2, ..., an >> 2}v_int32 a, b;v_int32 a_and_b = v_and(a, b); // {a1 & b1, ..., an & bn} -
比较运算符:我们可以使用
v_lt(<)、v_gt(>)、v_le(<=)、v_ge(>=)、v_eq(==)和v_ne(!=)来比较两个寄存器之间的值。由于每个寄存器包含多个值,这些操作不会得到单个 bool 值。相反,对于真值,所有位都被置为一(8 位为 0xff,16 位为 0xffff 等),而假值返回的位全部为零。// let us consider the following code is run in a 128-bit registerv_uint8 a; // a = {0, 1, 2, ..., 13, 14, 15}v_uint8 b; // b = {15, 14, 13, ..., 2, 1, 0}v_uint8 c = v_lt(a, b); // c = {255, 255, 255, ..., 0, 0, 0}/*let us look at the first 4 values in binarya = |00000000|00000001|00000010|00000011|b = |00001111|00001110|00001101|00001100|c = |11111111|11111111|11111111|11111111|If we store the values of c and print them as integers, we will get 255 for true values and 0 for false values.*/在一台支持 256 位寄存器的计算机中:
v_int32 a; // a = {1, 2, 3, 4, 5, 6, 7, 8}v_int32 b; // b = {8, 7, 6, 5, 4, 3, 2, 1}v_int32 c = v_lt(a, b); // c = {-1, -1, -1, -1, 0, 0, 0, 0}/*The true values are 0xffffffff, which in signed 32-bit integer representation is equal to -1.*/ -
最小/最大运算:我们可以使用 v_min() 和 v_max() 函数,返回包含两个寄存器逐元素最小值或最大值的寄存器:
v_int32 a; // {a1, ..., an}v_int32 b; // {b1, ..., bn}v_int32 mn = v_min(a, b); // {min(a1, b1), ..., min(an, bn)}v_int32 mx = v_max(a, b); // {max(a1, b1), ..., max(an, bn)}
注: 比较和最小/最大运算符不适用于 64 位整数。按位移位和逻辑运算符仅适用于整数值。按位移位仅适用于 16、32 和 64 位寄存器。
-
归约运算:v_reduce_min()、v_reduce_max() 和 v_reduce_sum() 返回一个单一值,表示整个寄存器的最小值、最大值或总和:
v_int32 a; // a = {a1, ..., a4}int mn = v_reduce_min(a); // mn = min(a1, ..., an)int sum = v_reduce_sum(a); // sum = a1 + ... + an -
掩码运算:掩码运算允许我们在宽寄存器中复制条件判断。这些包括:
-
v_check_all()——返回一个 bool 值,如果寄存器中的所有值都小于零则为真。
-
v_check_any()——返回一个 bool 值,如果寄存器中有任何一个值小于零则为真。
-
v_select()——返回一个寄存器,它基于掩码混合两个寄存器。
v_uint8 a; // {a1, .., an}v_uint8 b; // {b1, ..., bn}v_int32x4 mask: // {0xff, 0, 0, 0xff, ..., 0xff, 0}v_uint8 Res = v_select(mask, a, b) // {a1, b2, b3, a4, ..., an-1, bn}/*"Res" will contain the value from "a" if mask is true (all bits set to 1),and value from "b" if mask is false (all bits set to 0)We can use comparison operators to generate mask and v_select to obtain results based on conditionals.It is common to set all values of b to 0. Thus, v_select will give values of "a" or 0 based on the mask.*/
-
在以下部分中,我们将对一个单通道的简单卷积函数进行向量化,并将结果与标量实现进行比较。
注: 并非所有算法都能通过手动向量化得到改进。事实上,在某些情况下,编译器可能会对代码进行自动向量化(autovectorize),从而让标量实现产生更快的结果。
你可以从上一个教程中了解更多关于卷积的知识。我们使用与上一个教程相同的朴素实现,并将其与向量化版本进行比较。
完整的教程代码在此处。
我们将首先实现一维卷积,然后对其进行向量化。二维向量化卷积将跨行执行一维卷积以产生正确的结果。
一维卷积:标量
Section titled “一维卷积:标量”void conv1d(Mat src, Mat &dst, Mat kernel){ int len = src.cols; dst = Mat(1, len, CV_8UC1);
int sz = kernel.cols / 2; copyMakeBorder(src, src, 0, 0, sz, sz, BORDER_REPLICATE);
for (int i = 0; i < len; i++) { double value = 0; for (int k = -sz; k <= sz; k++) value += src.ptr<uchar>(0)[i + k + sz] * kernel.ptr<float>(0)[k + sz];
dst.ptr<uchar>(0)[i] = saturate_cast<uchar>(value); }}-
我们首先设置变量,并在 src 矩阵的两侧添加边界,以处理边缘情况:
int len = src.cols;dst = Mat(1, len, CV_8UC1);int sz = kernel.cols / 2;copyMakeBorder(src, src, 0, 0, sz, sz, BORDER_REPLICATE); -
对于主循环,我们选择一个索引 i,并使用 k 变量把它连同卷积核一起向两侧偏移。我们把值存储在 value 中并将其写入 dst 矩阵:
for (int i = 0; i < len; i++){double value = 0;for (int k = -sz; k <= sz; k++)value += src.ptr<uchar>(0)[i + k + sz] * kernel.ptr<float>(0)[k + sz];dst.ptr<uchar>(0)[i] = saturate_cast<uchar>(value);}
一维卷积:向量
Section titled “一维卷积:向量”现在我们来看一维卷积的向量化版本:
void conv1dsimd(Mat src, Mat kernel, float *ans, int row = 0, int rowk = 0, int len = -1){ if (len == -1) len = src.cols;
Mat src_32, kernel_32;
const int alpha = 1; src.convertTo(src_32, CV_32FC1, alpha);
int ksize = kernel.cols, sz = kernel.cols / 2; copyMakeBorder(src_32, src_32, 0, 0, sz, sz, BORDER_REPLICATE);
int step = VTraits<v_float32x4>::vlanes(); float *sptr = src_32.ptr<float>(row), *kptr = kernel.ptr<float>(rowk); for (int k = 0; k < ksize; k++) { v_float32 kernel_wide = vx_setall_f32(kptr[k]); int i; for (i = 0; i + step < len; i += step) { v_float32 window = vx_load(sptr + i + k); v_float32 sum = v_add(vx_load(ans + i), v_mul(kernel_wide, window)); v_store(ans + i, sum); }
for (; i < len; i++) { *(ans + i) += sptr[i + k]*kptr[k]; } }}-
在我们的例子中,卷积核是 float 类型。由于卷积核的数据类型最大,我们把 src 转换为 float32,形成 src_32。我们也像朴素情况那样添加边界:
Mat src_32, kernel_32;const int alpha = 1;src.convertTo(src_32, CV_32FC1, alpha);int ksize = kernel.cols, sz = kernel.cols / 2;copyMakeBorder(src_32, src_32, 0, 0, sz, sz, BORDER_REPLICATE); -
现在,对于 kernel 中的每一列,我们计算该值与所有长度为
step的 window 向量的标量积。我们把这些值累加到 ans 中已存储的值上:- 我们声明指向 src_32 和 kernel 的指针,并为每个卷积核元素运行一个循环:
int step = VTraits<v_float32x4>::vlanes();float *sptr = src_32.ptr<float>(row), *kptr = kernel.ptr<float>(rowk);for (int k = 0; k < ksize; k++){- 我们用当前卷积核元素加载一个寄存器。一个窗口从 0 移动到 len - step,它与 kernel_wide 的乘积被累加到 ans 中存储的值上。我们把值存回 ans:
v_float32 kernel_wide = vx_setall_f32(kptr[k]);int i;for (i = 0; i + step < len; i += step){v_float32 window = vx_load(sptr + i + k);v_float32 sum = v_add(vx_load(ans + i), v_mul(kernel_wide, window));v_store(ans + i, sum);}- 由于长度可能不能被步长整除,我们直接处理剩余的值。尾部值的数量将始终小于 step,不会显著影响性能。我们把所有值存储到 ans(它是一个 float 指针)中。我们也可以直接把它们存储在
Mat对象中:
for (; i < len; i++){*(ans + i) += sptr[i + k]*kptr[k];}- 下面是一个迭代示例:
For example:kernel: {k1, k2, k3}src: ...|a1|a2|a3|a4|...iter1:for each idx i in (0, len), 'step' idx at a timekernel_wide: |k1|k1|k1|k1|window: |a0|a1|a2|a3|ans: ...| 0| 0| 0| 0|...sum = ans + window * kernel_wide= |a0 * k1|a1 * k1|a2 * k1|a3 * k1|iter2:kernel_wide: |k2|k2|k2|k2|window: |a1|a2|a3|a4|ans: ...|a0 * k1|a1 * k1|a2 * k1|a3 * k1|...sum = ans + window * kernel_wide= |a0 * k1 + a1 * k2|a1 * k1 + a2 * k2|a2 * k1 + a3 * k2|a3 * k1 + a4 * k2|iter3:kernel_wide: |k3|k3|k3|k3|window: |a2|a3|a4|a5|ans: ...|a0 * k1 + a1 * k2|a1 * k1 + a2 * k2|a2 * k1 + a3 * k2|a3 * k1 + a4 * k2|...sum = sum + window * kernel_wide= |a0*k1 + a1*k2 + a2*k3|a1*k1 + a2*k2 + a3*k3|a2*k1 + a3*k2 + a4*k3|a3*k1 + a4*k2 + a5*k3|
注: 函数参数还包括 row、rowk 和 len。当把该函数用作二维卷积的中间步骤时,会用到这些值。
假设我们的卷积核有 ksize 行。为了计算某一特定行的值,我们计算前 ksize/2 行和后 ksize/2 行与对应卷积核行的一维卷积。最终值就是各个一维卷积之和:
void convolute_simd(Mat src, Mat &dst, Mat kernel){ int rows = src.rows, cols = src.cols; int ksize = kernel.rows, sz = ksize / 2; dst = Mat(rows, cols, CV_32FC1);
copyMakeBorder(src, src, sz, sz, 0, 0, BORDER_REPLICATE);
int step = VTraits<v_float32x4>::vlanes();
for (int i = 0; i < rows; i++) { for (int k = 0; k < ksize; k++) { float ans[N] = {0}; conv1dsimd(src, kernel, ans, i + k, k, cols); int j; for (j = 0; j + step < cols; j += step) { v_float32 sum = v_add(vx_load(&dst.ptr<float>(i)[j]), vx_load(&ans[j])); v_store(&dst.ptr<float>(i)[j], sum); }
for (; j < cols; j++) dst.ptr<float>(i)[j] += ans[j]; } }
const int alpha = 1; dst.convertTo(dst, CV_8UC1, alpha);}-
我们首先初始化变量,并在 src 矩阵的上下添加边界。左右两侧由一维卷积函数处理:
int rows = src.rows, cols = src.cols;int ksize = kernel.rows, sz = ksize / 2;dst = Mat(rows, cols, CV_32FC1);copyMakeBorder(src, src, sz, sz, 0, 0, BORDER_REPLICATE);int step = VTraits<v_float32x4>::vlanes(); -
对于每一行,我们计算其上下各行的一维卷积。然后把值累加到 dst 矩阵中:
for (int i = 0; i < rows; i++){for (int k = 0; k < ksize; k++){float ans[N] = {0};conv1dsimd(src, kernel, ans, i + k, k, cols);int j;for (j = 0; j + step < cols; j += step){v_float32 sum = v_add(vx_load(&dst.ptr<float>(i)[j]), vx_load(&ans[j]));v_store(&dst.ptr<float>(i)[j], sum);}for (; j < cols; j++)dst.ptr<float>(i)[j] += ans[j];}} -
最后我们把 dst 矩阵转换为 8 位
unsigned char矩阵:const int alpha = 1;dst.convertTo(dst, CV_8UC1, alpha);
在本教程中,我们使用了一个水平梯度卷积核。两种方法得到相同的输出图像。
运行时间的提升因情况而异,取决于你的 CPU 中可用的 SIMD 能力。