MIT-BEVFusion系列九--CUDA-BEVFusion部署5 update、矩阵逆运算、预计算

该系列文章与qwe、Dorothea一同创作，喜欢的话不妨点个赞。

接下来，是bevfusion中，另一个非常重要的部分- - -update。

create_core 是配置参数，反序列化engine，为输出分配内存。

update就是在正式进行forward前，做好其余工作。包括以下几点：

• 为矩阵赋值，矩阵计算、矩阵的逆运算
• 预计算（认为视锥点与360 * 360 * 1的体素格子之间的映射关系，在模型推理前就能计算出来。）

更新 bevfusion::Core 中的变换矩阵数据并进行预计算

• 构建 core 时，很多数据在GPU开辟好了内存，但是没赋值。调用 update 更新矩阵参数，进行拷贝。

更新图像变换矩阵 camera_depth_->update

上图的228行，主要内容见下图。update内部主要就是一个cuda的异步拷贝操作。将cpu上的增广矩阵、lidar2image矩阵数据，拷贝到GPU上。

将 图像增强矩阵 和 lidar2image 矩阵拷贝到在 device 上分配好的内存中。

更新真实世界中几何空间(体素网格)的变换矩阵 camera_geometry_->update

上上图229行，具体内容如下。
270行，求6个视角的相机的内参矩阵的逆矩阵。
271行，求6个视角的相机的图像增广矩阵的逆矩阵。
275-279行，把6个视角相机的矩阵拷贝到GPU
280行，为GPU上，名字是keep_count_的用于计数的变量赋值

求矩阵的逆的方法实现。

• 计算 逆相机内参矩阵 和 逆图像增强矩阵，并将这两个矩阵和 camera2lidar 拷贝到在 device 上分配好的内

• 解释
  

  • 为了方便观看，上图 89 行处，增加下方代码调试。重新编译（bash tools/run.sh）

  for (int i=0; i<4; i++){
    for (int j=0; j<4; j++){
      printf("%f ", m[i*4+j]);
    }
    printf("\n");
  }

• 调试结果(取一个矩阵)
  
  

(0)1266.417236 (1)0.000000     (2)816.267029 (3)0.000000 
(4)0.000000    (5)1266.417236  (6)491.507080 (7)0.000000 
(8)0.000000    (9)0.000000     (10)1.000000  (11)0.000000 
(12)0.000000   (13)0.000000    (14)0.000000  (15)1.000000

• 逆矩阵计算解释
• • 第一步，求行列式
    • 拉普拉斯展开
      对于一个 $4\times 4$ 矩阵:
      M = ( m [ 0 ] m [ 1 ] m [ 2 ] m [ 3 ] m [ 4 ] m [ 5 ] m [ 6 ] m [ 7 ] m [ 8 ] m [ 9 ] m [ 10 ] m [ 11 ] m [ 12 ] m [ 13 ] m [ 14 ] m [ 15 ] ) M=\left(\begin{array}{cccc} m[0] & m[1] & m[2] & m[3] \\ m[4] & m[5] & m[6] & m[7] \\ m[8] & m[9] & m[10] & m[11] \\ m[12] & m[13] & m[14] & m[15] \end{array}\right) M= ​m[0]m[4]m[8]m[12]​m[1]m[5]m[9]m[13]​m[2]m[6]m[10]m[14]​m[3]m[7]m[11]m[15]​ ​
    • 行列式的拉普拉斯展开（按第一列展开）是:$$\begin{array}{rl}& \det (M)=m[0]\times (-1{)}^{1+1}\times \det \left(M_{11}\right)+m[1]\times (-1{)}^{1+2}\times \det \left(M_{12}\right)+m[2]\times (-1{)}^{1+3}\times \det \left(M_{13}\right)+\\ & m[3]\times (-1{)}^{1+4}\times \det \left(M_{14}\right)\end{array}$$
    • 这里, $M_{11},M_{12},M_{13},M_{14}$ 是 $3\times 3$ 的子矩阵, 分别得到通过从原始矩阵 $\mathrm{M}$ 中去掉第一列和第 $\mathrm{i}$ 行 $(\mathrm{i}=1,2,3,4)$ 来形成。
      例如, $M_{11}$ 是去掉第一行和第一列的子矩阵:
      M 11 = ( m [ 5 ] m [ 6 ] m [ 7 ] m [ 9 ] m [ 10 ] m [ 11 ] m [ 13 ] m [ 14 ] m [ 15 ] ) M_{11}=\left(\begin{array}{ccc} m[5] & m[6] & m[7] \\ m[9] & m[10] & m[11] \\ m[13] & m[14] & m[15] \end{array}\right) M11​= ​m[5]m[9]m[13]​m[6]m[10]m[14]​m[7]m[11]m[15]​ ​
    • 其中内参矩阵是3 * 3 的形状，变成4 * 4矩阵。一定存在 m[15]=1，m[12]=m[13]=m[14]=m[3]=m[7]=m[11]=0，
    • 所以行列式$det=m[0]\ast (m[5]\ast m[10]-m[9]\ast m[6])-m[1]\ast (m[4]\ast m[10]-m[6]\ast m[8])+m[2]\ast (m[4]\ast m[9]-m[5]\ast m[8]);$
• • 第二步，m中元素为0的代数余子式 乘以 逆行列式det

图像增广矩阵我也截图放在了下方

视锥中的点转换到lidar坐标系下 (预计算)

283行compute_geometry_kernel核函数

• 核心思想：不同相机计算视锥点的方式是一样的。

• 入参：

  • numel_frustum_：视锥点数量

  • frustum_ ： 视锥点的坐标，类型float3*

  • reinterpret_cast(camera2lidar_) 矩阵：将 float*类型数据，转换成float4_

  • reinterpret_cast<const float4*>(camera_intrinsics_inverse_) 相机内参逆矩阵

  • reinterpret_cast<const float4*>(img_aug_matrix_inverse_) 图像增光逆矩阵

  • param_.bx： nvtype::Float3类型。begin x
    

  • param_.dx：nvtype::Float3类型。delta x
    

  • param_.nx: nvtype::Int3类型
    

  • keep_count_ ： unsigned int* 类型指针。 默认值是0

  • param_.geometry_dim : nvtype::Int3类型
    

• 结果：

  • 结论：
    • ranks_与geometry_计算方式化简后一致，ranks_考虑了边界，超出边界的索引为0，而geometry_没有考虑。ranks_的计算公式更为清晰。
  • ranks_：在核函数里赋值。
    • 数量：1993728=1 * 6 * 118 * 32 * 88
    • _pid做为偏移量：int _pid = icamerea * numel_frustum + tid;表示6个相机的视锥点。
    • 偏移量对应的值：在边界内的视锥点，经过边界筛选
    • 值为0的点表示超出边界的。值不为0的是在边界的点
    • 数据存在，/datav/NVIDIA-AI-IOT/Lidar_AI_Solution/CUDA-BEVFusion/Mydebug/tmp_rank_before_sort.txt
  • geometry_: int32_t类型指针，存的6个相机视锥点在360*360大小的体素索引。所有都包含
    • 数量：1993728
    • 已开辟内存：numel_geometry_ * sizeof(int32_t)
    • 数据存在：/datav/NVIDIA-AI-IOT/Lidar_AI_Solution/CUDA-BEVFusion/Mydebug/tmp_geometry_.txt
  • keep_count_:统计落在体素范围内点的数量。

• 开辟线程：至少numel_frustum_个，即至少88 * 32 * 118=332288。512个线程为一个block。一共numel_frustum_ / 512向上取整个 。

• cuda_linear_launch宏值得看一下。

• 上图137行，计算全局线程索引，并获得当前线程（3288118=332288个线程中的一个）对应的视锥点的位置 (w, h, d)，是一个float3类型的数据。

  • 注意：这里没有相机个数的6，是因为6个相机视锥点计算相同。不同的是内外参矩阵。

• 上图141行，进行逆图像增强。

• 142、143行，2维到3维转换。

• 144行，逆相机内参转换

• 146行，camera2lidar转换这样视锥点位置就转换到lidar坐标系下了。

• 151～153行，将lidar坐标系下的点转换到体素网格中，先将lidar坐标系下的点转换到体素坐标系中。

图中虚线为360360的格子范围，超出这个范围的点不需要。*

geometry_out需要后续理解。

将超出体素范围x∊[0，360)，y∊[0，360)，z∊[0, 1) 的点通过掩码的形式进行筛选，超出范围的使用0表示，在范围内的通过体素索引表示。keep_count通过原子操作来进行加1，记录未超出范围的点的个数。可以将点的坐标打印出来printf("coords.x :%d, coords.y: %d\n", coords.x , coords.y);，发现的确有大于360的情况出现:

关联视锥中点的索引与真实世界体素的索引

• 上图330行，一个arange的kernel版本，结果存在indices_中，从0-numel_geometry(1993728)
  • 初始化indices_为6个视锥中点的索引，一个arange的核函数版
  • numel_geometry_ 6*32*88*118 indices_ 6*32*88*118 个 int
  • 初始化indices用于表示6个相机的视锥点的索引。

ranks_、indices_排序

• 排序前indices
  

• 排序后indices
  

  

• 排序前rank_结果
  

• 排序后rank_结果
  

void stable_sort_by_key(const thrust::detail::execution_policy_base<DerivedPolicy> &exec, RandomAccessIterator1 keys_first, RandomAccessIterator1 keys_last, RandomAccessIterator2 values_first, StrictWeakOrdering comp)

• thrust::stable_sort_by_key是Thrust库中的一个函数，用于在GPU中执行稳定的键值对排序。对提供的键序列进行排序，同时重新排列相应的值，以保持键-值对关系。当两个元素的键相等时，它们在排序后的序列中的相对顺序保持不变。这是一个并行化的排序算法，可以显著提高大数据集的排序速度。

• 体素索引ranks_(键序列)进行升序排序，indices_(值)也会根据排序进行相应的顺序变化，这样依然可以通过视锥索引indices来找到升序后的体素索引, 可以按照下面的公式来理解 这个排序api的作用: ranks[sorted_indices[i]] = sorted_ranks[i], 这个式子中的i 来自下面这个式子(见下文):i = start_interval + remain_ranks(159963)。ranks_的意义是视锥点到体素点的映射关系，可能存在多个视锥点映射到同一个体素点的情况。也可能存在体素点超出范围的情况，此时ranks_里的值是0

• gpt解释
  这行代码是使用了 Thrust 库在 CUDA 环境中对数据进行排序的一个示例。Thrust 是一个基于 C++ 的并行编程库，专为 CUDA 设计，提供了类似于标准模板库（STL）的接口来处理 GPU 上的数据。让我们逐部分解释这行代码：

  1. thrust::stable_sort_by_key：

     • 这是 Thrust 库中的一个函数，用于根据“键”对一组“键-值”对进行排序。它保证具有相同键值的元素将保持它们原始的相对顺序，这就是“稳定排序”（stable sort）的含义。

  2. **thrust::cuda::par.on(_stream)**：

  • thrust::cuda::par 是一个执行策略，它告诉 Thrust 函数应该在 CUDA 并行环境下执行。
  • .on(_stream) 是一个方法，它指定了 CUDA 流 _stream 用于执行排序操作。CUDA 流代表了一系列异步执行的操作队列，使用流可以实现在 GPU 上的并行计算和数据传输。
  3. **ranks_****, ****ranks_ + numel_geometry_**：
  • 这些是指向“键”数组的指针。ranks_ 是指向数组起始的指针，ranks_ + numel_geometry_ 是指向数组末尾的指针。这个范围内的元素将根据它们的键值进行排序。
  4. **indices_**：
  • 这是一个与“键”数组相对应的“值”数组。当键数组中的元素根据键值重新排序时，indices_ 数组中的元素也将相应地重新排列。
  5. **thrust::less<int>()**：
  • 这是一个比较函数对象，用于定义排序的标准。在这个例子中，它表示使用小于（less than）的关系对整数键进行排序。

  总的来说，这行代码的作用是在 CUDA 环境下，使用指定的 CUDA 流 _stream 对键数组 ranks_ 进行稳定排序，并以此重新排列对应的值数组 indices_。使用 thrust::less<int>() 意味着排序是按照键的升序进行的。这种排序操作在并行计算中非常常见，尤其是在处理大量数据时。

确认体素间隔点

计算剩余的空间remain_ranks，主要用于后续计算间隔点时作为偏移量。

• 间隔点表示的是唯一的体素。

  • 因为在将视锥中的点转换到真实世界体素网格中，不同的点可能会占用同一个体素，通过间隔点来表示这些唯一的体素。
  • 上面这张图中，假如视锥中的一部分点转换到右边的体素网格中后，分布在了5个体素中，这里的间隔点就是5个。

• 通过threads来分配计算后续计算间隔点需要的线程数。

  • 这里减去了1是因为后面计算时从第二个间隔点开始计算，第一个间隔点索引会设置为0。（可以参考python学习的时候的boardmix）
  • 初始化间隔点数量interval_starts_size_为0，用于储存间隔点的数量。将interval_starts_中的第一个间隔点的索引初始化为0。

interval_starts_kernel核函数

• cuda_linear_launch宏命令启动核函数

• 入参：

  • threads: unsigned int，即下图num，线程数。含义：成功落在体素范围内的点的数量 减1
  • remain_ranks：unsigned int，即下图remain，含义：落在体素范围外的点的个数。
  • numel_geometry_：unsigned int，即下图total，含义：63288*118
  • ranks_：即下图ranks，含义：点对应的体素索引，超出边界的点对应ranks_值为0.此时已升序排序。
    • 已在gpu开辟内存 numel_geometry_ * sizeof(int32_t)。numel_geometry_63288*118
  • indices_:指针，即下图indices，含义：arange生成的规则的所有点的索引（ranks_排序前），ranks_排序后已经跟随ranks_变化而发生变化。
    • 已在gpu开辟内存 numel_geometry_ * sizeof(int32_t)
  • interval_starts_ + 1指针：即下图interval_starts。interval_starts_指针默认偏移1传入，数据从这里开始记录。含义：间隔点索引（落在体素范围内的点的间隔点索引）。
    • 已在gpu开辟内存 numel_geometry_ * sizeof(int32_t)
  • interval_starts_size_：指针，即下图interval_starts_size。含义：间隔点数量。
    • 已在gpu开辟内存 sizeof(int32_t)大小

• 结果：

  • **interval_starts_ + 1 ，**开始存入间隔点索引，不包含第一个间隔点（83）的索引。
  • **interval_starts_size_ ，**统计有多少个间隔点，不包含0-83的那个间隔点。

• 调试：

  • 建议如图所示，在排序后加入调试打印代码
    

thrust::stable_sort_by_key(thrust::cuda::par.on(_stream), ranks_, ranks_ + numel_geometry_, indices_, thrust::less<int>());       // 将ranks_进行升序排序，indices记录变化后的值的序列。通过视锥中点的索引可以ranks_中找到转换到体素中的索引，这里将体素索引进行排序，将视锥中点的索引打乱，未被占用的体素中没有索引，默认值为0，被排序到前端。
    // test part start
    int32_t* tmp_h_ranks = new int32_t[numel_geometry_];
    checkRuntime(cudaMemcpy(tmp_h_ranks, ranks_, numel_geometry_ * sizeof(int32_t), cudaMemcpyDeviceToHost));
    FILE* f = fopen("saved_data_from_prog/sorted_ranks.txt", "w");
    for (int i = 0; i < numel_geometry_; ++i) fprintf(f, "%d\n", tmp_h_ranks[i]);
    delete [] tmp_h_ranks;
    fclose(f);

• 打印信息保存在/datav/NVIDIA-AI-IOT/Lidar_AI_Solution/CUDA-BEVFusion/Mydebug/tmp_interval_starts_before_sort.txt
  - 结果第一个数是0 ，后续存入
  
  
  

• 上图77行，查看当前体素索引与前一个体素索引是否相同，如果不相同，interval_starts_size通过原子操作加1。

  • 例如ranks=[0,0,83,86]，第一次ranks[i] != ranks[i - 1]条件成立，是对比86与83。此时，0没有计算

• 并在interval_starts中添加间隔点索引，这样通过存储的相邻的两个间隔点索引就可以计算出重叠了多少个视锥点。后续更新计数器counter_host_为间隔点数量，并对流进行同步，保证核函数执行完毕。
  

• 上图346行，更新间隔点数量n_intervals，这里的加1是添加了第一个间隔点（0与83这个间隔）。

interval_starts排序

• 上图355行，对间隔点索引interval_starts_进行升序排序。
  • 这里只对前间隔点个数量的间隔点索引进行排序，即在interval_starts_kernel核函数中进行赋值的部分。
  • 结果：前96801个点进行了排序, 一共有1993728个数据，其余都是0。
    

• 总结interval_starts的作用 ：

  1. 记录rank中每个起始点，相对于第一个不是0的位置的偏移量。理解需要参考下图。

     • 以rank中的83为例子。他是第一个不是0的位置。83又是第一个起始点。所以偏移量为0

       • 86是第二个起始点，距离83偏移量为1
       • 87是第三个起始点，距离83偏移量为7
       • 排序后rank_结果
         

       

  2. 这个**偏移量**** 加上 reamain_ranks的数量（159963）就是这个数据在排序后的rank中距离首地址的****偏移量 **，

     1. 以interval_starts_图片为例，第一个数据是0， 0+159963 = 159963。对应ranks导出数据 159963索引位置的数据，即83。 （图片是文本编辑器，行号从1开始，忽略它）

核函数collect_starts_kernel

• 入参
  • n_intervals_: unsigned int，即下图num，为96801，线程数。含义：总的间隔数（包含筛选掉的点和留下来的点的间隔，即0与83的间隔）。
  • remain_ranks：unsigned int，即下图remain，含义：落在体素范围外的点的个数159663。
  • numel_geometry_：unsigned int，即下图numel_geometry，含义：63288*118=1993728
  • indices_:指针，即下图indices，含义：arange生成的规则的所有点的索引（ranks_排序前），ranks_排序后已经跟随ranks_变化而发生变化。
    • 已在gpu开辟内存 numel_geometry_ * sizeof(int32_t)
  • **interval_starts_ **指针：即下图interval_starts，此时已是完全体。含义：间隔点索引（包含全部点，即包含0与83的那个间隔）。
    • 已在gpu开辟内存 numel_geometry_ * sizeof(int32_t)
  • geometry_：指针，即下图geometry。含义：存的6个相机视锥点在360*360大小的体素索引。
    • 已在gpu开辟内存 sizeof(int32_t)大小，已在1.8.3.1中存入参数。
  • intervals_：指针，即下图intervals，核函数结果存这里
    • 已在gpu开辟内存 numel_geometry_ * sizeof(int3)
• 结果：intervals中存入了数据。

• intervals三个参数
  val.x表示当前间隔点对应ranks中体素的索引。
  val.y表示下一个间隔点对应ranks中体素的索引，对最后一个间隔点计算时，remain和重叠点个数相加。
  val.z表示geometry的索引。bevpool输出特征索引(在360*360*1的上特征图上的索引)，这里需要通过视锥索引来获取，在compute_geometry_kernel函数中定义。
  最后将这三个值通过三元结构体来存储，保存在intervals_中，索引为间隔点索引。

• 理解上图92行，需要知道indices与interval_start、geometry的关系

  • 以，排序前ranks为例，索引为0，值为63905，对应排序前indices为0。

    • ranks排序后，如果想
    • 因此需要用排序后 indices值为0这个条件去找。发现排序后 indices值为0
      

  • 第一步：这里解释了interval_starts[i] + remain得到这个起始点在排序后的ranks中相对首地址的偏移量。

    • 以i=1为例，interval_starts[i]=1，interval_starts[i] + remain=159964。在排序后的ranks中,ranks[159964]=86
      

• 第二步：indices存的是ranks排序前的索引。indices[interval_starts[i] + remain]即indices[159964]=1658655。得到的是ranks排序前的索引1658655
  

第三步：拿到排序前的索引，就能得到对应视锥上的点的索引。注意：geometry是没排序的。
geometry[indices[interval_starts[i] + remain]]= 86

• 为什么用indices再geometry，而不直接用ranks，结果不都是86吗？？？
  - 因为geometry 保留了所有的视锥点 转换到 体素上的索引。ranks则经过边界筛选
  

文章来源地址https://www.toymoban.com/news/detail-838413.html

• intervals后续的作用
  • 假如从intervals中取第一个间隔点数据，例如a=intervals[0]，a是间隔点起始点索引（已经加了remain了），b=indices[a](这里的indices是已经排序过了的)就能找到这个间隔点没排序前，在geometry中的索引b。通过geometry[b]就能找到对应的点， 也就是intervals.z的值。

到了这里，关于MIT-BEVFusion系列九--CUDA-BEVFusion部署5 update、矩阵逆运算、预计算的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！