OpenCV使用通用内部函数对代码进行矢量化(35)

    返回：OpenCV系列文章目录（持续更新中......）

上一篇：如何使用 OpenCV parallel_for_并行化代码(34)

下一篇：OpenCV 4.9使用通用内部函数对代码进行矢量化(36)

​

• 目标

  本教程的目标是提供使用通用内部函数功能矢量化 C++ 代码以实现更快运行时的指南。我们将简要介绍 SIMD 内部函数以及如何使用宽寄存器，然后是关于使用宽寄存器的基本操作的教程。

  理论

  在本节中，我们将简要介绍一些概念，以更好地帮助理解该功能。

  内在因素

  内部函数是由编译器单独处理的函数。这些功能通常经过优化，以最有效的方式执行，因此运行速度比正常实现更快。但是，由于这些函数依赖于编译器，因此很难编写可移植应用程序。

  SIMD的

  SIMD 代表 单指令，多数据。SIMD 内部函数允许处理器对计算进行矢量化。数据存储在所谓的寄存器中。寄存器可以是 128 位、256 位或 512 位宽。每个寄存器存储相同数据类型的多个值。寄存器的大小和每个值的大小决定了总共存储的值数。

  根据您的 CPU 支持的指令集，您可以使用不同的寄存器。要了解更多信息，请看这里

  通用内在函数

  OpenCV 的通用内部函数为 SIMD 矢量化方法提供了抽象，并允许用户使用内部函数，而无需编写特定于系统的代码。

  OpenCV Universal Intrinsics 支持以下指令集：

• 支持各种类型的 128 位寄存器，适用于各种架构，包括
  • x86（SSE/SSE2/SSE4.2），
  • 手臂（霓虹灯），
  • PowerPC（VSX），
  • MIPS（MSA）。
• x86（AVX2） 支持 256 位寄存器和
• x86（AVX512） 支持 512 位寄存器
• 寄存器结构
• 加载和存储
• 数学运算
• 减少和遮罩

寄存器结构

通用内部函数集将每个寄存器实现为基于特定 SIMD 寄存器的结构。所有nlanes类型都包含枚举，该枚举提供该类型可以保存的确切值数。这样就无需在实现过程中对值的数量进行硬编码。

注意

每个寄存器结构都位于命名空间cv下。

• 可变大小的寄存器：这些结构没有固定的大小，它们的确切位长度是在编译过程中根据可用的 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 位浮点数（双精度）

  •  v_uint8 a;                            // a is a register supporting uint8(char) data
      int n = a.nlanes;                     // n holds 32

     可用的数据类型和大小：

    类型	大小（以位为单位）
    uint	8, 16, 32, 64
    int	8, 16, 32, 64
    float	32, 64

• 恒定大小的寄存器：这些结构具有固定的位大小，并保存恒定数量的值。我们需要知道系统支持哪些SIMD指令集，并选择兼容的寄存器。仅当需要精确的位长度时才使用它们。

  每个结构都遵循以下约定：

  v_[type of value][size of each value in bits]x[number of values]

  假设我们要存储

  • 128 位寄存器中的 32 位（以位为单位）有符号整数。由于寄存器大小已经已知，我们可以找出寄存器中的数据点数 （128/32 = 4）：

  v_int32x8 reg1                       // holds 8 32-bit signed integers.

512 位寄存器中的 64 位浮点数：

  v_float64x8 reg2                     // reg2.nlanes = 8

载和存储操作

现在我们知道了寄存器的工作原理，让我们看看用于用值填充这些寄存器的函数。

• 加载：加载函数允许您将值加载到寄存器中。
  • 构造函数 - 在声明寄存器结构时，我们可以提供寄存器将从中获取连续值的内存地址，也可以将值显式提供为多个参数（显式多个参数仅适用于常量大小的寄存器）：

  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 hold

  v_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 函数 - 我们可以使用 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]

注意

load 函数假定数据未对齐。如果您的数据是对齐的，则可以使用vx_load_aligned()函数。

• 存储：存储函数允许您将寄存器中的值存储到特定的内存位置。
• 要将寄存器中的值存储到内存位置，可以使用 v_store（） 函数:

• •   float ptr[4];
      v_store(ptr, reg); // store the first 128 bits(interpreted as 4x32-bit floats) of reg into ptr.

    注意

    确保 ptr 与寄存器具有相同的类型。您还可以在执行操作之前将寄存器转换为正确的类型。简单地将指针类型转换为特定类型将导致对数据的错误解释。

二元运算符和一元运算符

通用内部函数集提供元素二元和一元运算。

• 算术：我们可以按元素对两个寄存器进行加、减、乘和除。寄存器必须具有相同的宽度并保存相同的类型。要将两个寄存器相乘，例如：

  v_float32 a, b;                          // {a1, ..., an}, {b1, ..., bn}
  v_float32 c;
  c = a + b                                // {a1 + b1, ..., an + bn}
  c = a * b;                               // {a1 * b1, ..., an * bn}

按位逻辑和移位：我们可以左移或右移寄存器中每个元素的位。我们还可以在两个寄存器之间按元素应用按位 &、|、^ 和 ~ 运算符：

  v_int32 as;                              // {a1, ..., an}
  v_int32 al = as << 2;                    // {a1 << 2, ..., an << 2}
  v_int32 bl = as >> 2;                    // {a1 >> 2, ..., an >> 2}

  v_int32 a, b;
  v_int32 a_and_b = a & b;                 // {a1 & b1, ..., an & bn}

比较运算符：我们可以使用 <、>、<=、>=、== 和 ！= 运算符来比较两个寄存器之间的值。由于每个寄存器都包含多个值，因此对于这些操作，我们没有得到一个布尔值。相反，对于真值，所有位都转换为 1（0xff 表示 8 位，0xffff 表示 16 位等），而 false 值返回转换为零的位。

  // let us consider the following code is run in a 128-bit register
  v_uint8 a;                               // a = {0, 1, 2, ..., 15}
  v_uint8 b;                               // b = {15, 14, 13, ..., 0}

  v_uint8 c = a < b;

  /*
      let us look at the first 4 values in binary

      a = |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.
  */
  ---
  // In a computer supporting 256-bit registers
  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 = (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 位寄存器

减少和遮罩

• Reduce 运算：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（） - 返回一个布尔值，如果寄存器中的所有值都小于零，则为 true。
• v_check_any（） - 返回一个布尔值，如果寄存器中的任何值小于零，则为 true。
• 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.
  */

示范

在下一节中，我们将对单通道的简单卷积函数进行矢量化，并将结果与标量实现进行比较。

注意

并非所有算法都可以通过手动矢量化进行改进。事实上，在某些情况下，编译器可能会自动矢量化代码，从而为标量实现生成更快的结果。

您可以从上一教程中了解有关卷积的更多信息。我们使用与上一教程相同的朴素实现，并将其与矢量化版本进行比较。

完整的教程代码在这里。

矢量化卷积

我们将首先实现一维卷积，然后对其进行矢量化。二维矢量化卷积将对行执行一维卷积，以产生正确的结果。

一维卷积：标量

  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.
  */

我们首先设置变量并在 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);
 }

一维卷积：向量

现在，我们将研究一维卷积的矢量化版本。

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];
 }
 }
}

在我们的例子中，内核是一个浮点数。由于内核的数据类型最大，我们将 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);

现在，对于内核中的每一列，我们计算该值的标量乘积，其中包含所有长度为 的窗口向量。我们将这些值添加到已存储step = VTraits<v_float32x4>::vlanes();的值

 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];
 }
 }

我们声明一个指向src_32和内核的指针，并为每个内核元素运行一个循环

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 - 步长，其带有 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);
 }

由于长度可能无法按步长整除，因此我们直接处理剩余值。尾值的数量将始终小于步进，并且不会显着影响性能。我们将所有值存储到 ans，这是一个浮点指针。我们也可以直接将它们存储在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 time
      kernel_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 功能。

 参考文献：

1、《Vectorizing your code using Universal Intrinsics》