论文阅读：VoxelNet（3D-detection）+代码复现

Brief

这是一篇发布在CVPR18年的点云检测方面的文章，这一篇文章算是第一篇只用点云去做3D检测的文章，类似于pointnet的点云分割中的地位。
这里是 paper
这里是 code
个人觉得仍然可以提高的部分：

• anchor的Overlap问题
• 文章的针对性比较强，anchor固定检测某一类的object
• 使用7个参数表达Bbox，实际上可以用6个。
• 引入二维热力图机制。
• anchor机制有点傻，一趟实验下来是70400个anchor对应2个gt
• 个人觉得目前的gt的方式和二维上都存在的问题是：并不是很紧致的，也就是有一些背景，三维的bbox会好一些，但为什么一定要用box来包裹而不是采用其他的结构。

1 Abstruct

• 本文直接处理点云输入，是end-to-end的网络；和以往的多视角投影的不同。
• 采用VFE（voxel feature encoding）网络对每个体素内的点进行编码再接上3DRPN检测。不使用手工特征处理，这好比CNN加入R-CNN一样。（感叹一下。作者写文章真的强）
• 该方法在KITTI car detection取得了sota的效果，超过了同时使用mino和lidar数据的方法。

2 Introduction

 大致文章思路都是这样：3D检测应用和其重要性（自动驾驶，housekeeping robots），lidar数据特性，目前的研究方案：

• 将lidar投影到某个视角作为输入的方法。
• 采用3D体素，但是对体素进行手工特征提取的方法。
• Pointnet和pointnet++的介绍和无法处理大量的点云问题，这也是本文要解决的主要问题之一。

RPN网络：

• 这里的RPN也就是2D中faster-RCNN的那个RPN，看了他的引用。
• 这种RPN是对密集型结构有效，这和原始的点云差距有点大，本文close the gap（作者再一次提到本文工作）

本文网络处理过程：

待补充

1.1 Related work

 

• 手工特征：在3D形状细节上很丰富的时候是能得到比较好的效果，但是对于变化后的场景和复杂的场景则不能。最大的问题是不能自己学习会存在bottleneck
• 2DBBOX推广到3D的做法：精准度会受到深度估计的限制。
• 一些方法利用雷达输入对每个voxel进行表示。一般是对voxel进行encode，再处理。
• Image 结合3D数据

1.2 贡献

• 一个直接处理点云，并且不采用人工特征的点云网络voxelnet
• 提出了一种voxel算法的快速实现方法，该想法得益于点云的稀疏特征结构和在voxel上的并行计算。
• sota效果

2 VoxelNet

2.1 architucture

• Feature Learning Network
• Convolutional Middle Layers
• RPN

2.1.1 Feature Learning Network

• 体素分区：对于输入点云所在的空间，我们采用$D,H,W$表示；定义每一个小的voxel的大小为$V_D,V_H,V_W$，通过除法就能得到网格的size。这里作者假设了$D,H,W$是体素网格的整数倍。（这会不会导致roi align和ROI pooling的问题）。这里可以加入先验知识进行分区吗，就比如说一些normal的一致性。
• grouping：根据点所在的voxel进行分组。由于数据的不均匀性，导致分组出来的组别的个数不尽一致。比如图上的一样，有的有三个，有的一个也没有。（这一步具体作用是什么作者没有表达清楚）
• 随机采样：一是为了减小计算负担，而是为了消除因为密度不一致带来的判别问题。采样方案只针对那些voxel数量比T个多的voxel，这种策略两点好处：省计算资源和（降低voxel之间的不平衡性，这可以降低采样偏差，在训练时增加更多的变化性。）
• Stacked Voxel Feature Encoding：
  这是比较重要的一点，作者表示最大的创新就在于VFE层级结构，具体用下面这个图来表达：

对每一个voxel来说，对里边的所有点而言均是如下操作，这里只描述对点$p_i(x_i,y_i,z_i,r_i)$，第四维的$r_i$表达的是反射率（1）计算所有点的均值,记做$(v_x,v_y,v_z)$（2）对于点$p_i$我们采用对中心的偏移量来增大表示其坐标$\widehat{p_i}=[x_i,y_i,z_i,r_i,x_i-v_x,y_i-v_y,z_i-v_z]$，（这样表达有什么好处，实则并没有信息增加,但是却把全局特征$x_i，y_i,z_i$和局部残差结合到了一起）(3)对每一个的$\widehat{p_i}$通过FCN转化到特征空间，特征空间表示可以将点特征整合来描述voxel所表达的表面特性。FCN层由一个线性层，一个BN，一个RELU层组成。这样就得到了point-wise-feature。得到的特征记做$f_i$（4）再通过element-wise maxpooling层,结合周围的所有的点得到locally-aggregated feature:$\tilde{f}$。（5）最后一步是将$\tilde{f}$和$f_i$结合为$f_{out}$。那么对于一个含有t个点的voxel,最终的特征我们采用集合表达方式为： V o u t = { f i o u t } i … t \mathbf{V}_{\mathrm{out}}=\left\{\mathbf{f}_{i}^{o u t}\right\}_{i \ldots t} Vout​={fiout​}i…t​

这一结构比较大的作用在于（1）既保留点的特征又结合了局部结构特征，使得其可以学习shape特征。(2)不同voxel之间结合可以得到voxel-wise特征。类比于point-wise特征。

• 稀疏张量表示:这是第一部分（特征提取）的最后一个层次结构。上面的结构通过作用于每一个非空体素，每一个都唯一关联于一个体素。得到的voxel-wise-feature可以用大小为$C\times D^{\prime }\times H^{\prime }\times W^{\prime }$表示，可以看最上面整体结构图。尽管可能会有100k个点，但是可能90%都是空的；将非空体素特征表示为稀疏张量大大降低了反向传播期间的内存使用和计算成本，并且是我们有效实现中的关键步骤。（能否使用voxel坐标表示其对应的一维而不是提取出来的那一种）

2.1.2 Convolutional Middle Layers

采用3D卷积的方式，采用$ConvMD(C_{in},C_{out},k,s,p)$的结构来表达，含义分别是$channel-in,channel-out,keisize(k\ast k\ast k),stride-size,padding$。基础结构有3D卷积，BN和Relu。

2.1.3 RPN

本文对RPN网络做一些关键性的修改。网络的整体结构如下。网路做如下描述：

• 网路的输入是通过 Convolutional Middle Layers出来的特征图。
• 存在3个全卷积的block

如图，每一个block都存在两个cov2D层，通过第一个层的stride=2来使得feature map的$w，l和h$变成 { w / 2 , l / 2 , h / 2 } \{w/2,l/2,h/2\} {w/2,l/2,h/2}；在每一个CON2D后面也都会有对应的BN和RELU层。

• concate

讲每一层的输出都上采样到一个统一的size（第一个block不用上采样）

• 映射为两个学习的目标：1. probability score map。2.regression map

2.2 loss Function

注意到上面的输出是对于socores的通道是2，分别对应中$positive$和$negative$的分数，而后的regression map存在14维，那是因为对于每一个回归的Bbox都用7维来表示：$\left(x_c^g,y_c^g,z_c^g,l^g,w^g,h^g,{\theta }^g\right)$,也即是中心位置+长宽高+yaw轴。这里另外两个旋转轴没有加入是默认为0，原因是地面水平。（假设加入进去是不是可以更好的定位物体，感觉是个不错的方向）。同理，我们假设预测的anchor位置为：$\left(x_c^a,y_c^a,z_c^a,l^a,w^a,h^a,{\theta }^a\right)$。因此我们可以定义如下的残差：
$$\begin{array}{l}\Delta x=\frac{x_c^g-x_c^a}{d^a},\Delta y=\frac{y_c^g-y_c^a}{d^a},\Delta z=\frac{z_c^g-z_c^a}{h^a}\\ \Delta l=\log \left(\frac{l^g}{l^a}\right),\Delta w=\log \left(\frac{w^g}{w^a}\right),\Delta h=\log \left(\frac{h^g}{h^a}\right)\\ \Delta \theta ={\theta }^g-{\theta }^a\end{array}$$
其中$d^a$是anchor框底部的对角线长度，表示为：
$$d^a=\sqrt{{\left(l^a\right)}^2+{\left(w^a\right)}^2}$$

文章说采用$d^a$来对$\Delta x和\Delta y$进行统一的归一化，为啥要采用$d^a$呢？（参考其引用的文献）

于是作者定义了损失函数：
$$\begin{array}{rl}L& =\alpha \frac{1}{N_{\mathrm{p}\mathrm{o}\mathrm{s}}}\sum _i{L}_{\mathrm{c}\mathrm{l}\mathrm{s}}\left(p_i^{\mathrm{p}\mathrm{o}\mathrm{s}},1\right)+\beta \frac{1}{N_{\mathrm{n}\mathrm{e}\mathrm{g}}}\sum _j{L}_{\mathrm{c}\mathrm{l}\mathrm{s}}\left(p_j^{\mathrm{n}\mathrm{e}\mathrm{g}},0\right)\\ & +\frac{1}{N_{\mathrm{p}\mathrm{o}\mathrm{s}}}\sum _i{L}_{\mathrm{r}\mathrm{e}\mathrm{g}}\left({\mathbf{u}}_i,{\mathbf{u}}_i^{\ast }\right)\end{array}$$
损失函数中的前面两项是归一化判别损失，其中的$p_i^{pos}$和$p_j^{neg}$分别表示$softmax$层对$a_i^{pos}$和$a_j^{neg}$d的分数。采用的是交叉熵表示，$\alpha$和$\beta$为正定平衡系数。最后一项是回归损失，采用的是SmoothL1损失。

2.3 Efficient Implementation

• VFE可以并行的处理voxel
  具体的方法图下图所示；

（1）第一步是初始化一个tensor，储存为$K\times T\times 7$的feature buffer。K表示一个有存在的非空的体素的最大值，T表示的每个voxel中点的最大值。7就是每一个点编码的特征维度。（2）第二步是对于点云中的任意一个点去查找它所在的voxel是否初始化，如果已经初始化了那么就检查是否大于T个点，没有的话就插入进去，否则不要这个点。假如这个voxel没有被初始化，那么久需要初始化并把这个点加入进去。（3）建立了input buffer之后呢，采用GPU的并行计算对点级别和voxel级别的VFE计算，作者提到，在VFE中的concate之后，把没有点的特征再次设置为0。最后，使用存储的坐标缓冲区，我们将计算的稀疏体素结构重新组织到密集体素网格

值得理解的是上面的图的储存示意图，Input buffer的正对面应该是$T\times K$，经过VFE提取后，加大了特征维数，然后再采用一种稀疏的表达方式。

3 Training Details

3.1 Network Details

3.1.1 Car Detection

输入：

• 输入的z,y,x范围是[−3, 1] × [−40, 40] × [0, 70.4]，多余的点就不要了。
• 每一个 voxel_size为：$v_D$ = 0.4, $v_H$ = 0.2, $v_W$ = 0.2。因此计算的个数为$10(D^{\prime })\times 400(H^{\prime })\times 352(W^{\prime })$。
  设置：
• T=35,最多的体素格子数。
• VFE层的：$VF{E}_1(7,32)andVF{E}_2(32,128)$。最后的$VF{E}_2$输出的特征维度为128.因此我们的稀疏点云tensor为：128 × 10 × 400 × 352。
• 进一步为了整合voxel-wise features。3个convolution middle layers按顺序如下：Conv3D(128, 64, 3,(2,1,1), (1,1,1)), Conv3D(64, 64, 3, (1,1,1), (0,1,1))和Conv3D(64, 64, 3, (2,1,1), (1,1,1))；最后会得到一个4D的tensor,大小为64 × 2 × 400 × 352。在送入RPN之前先整合一下，变成:128 × 400 × 352;
• anchor_size:文章仅仅只用了一种size：$l^a=3.9,w^a=1.6,h^a=1.56m,中心在z_c^a=-1.0，具有0度和90度两个旋转$（anchor在哪个feature上产生？这里的anchor我从代码中看到，是对于每一个体素也就是一个anchor，就是类比二维的feature map的每一个像素吧，就是128 × 400 × 352上，即对应中$C\times H\times W$）。

假设预测出的anchor在鸟瞰图中具有和gt最大的IOU或则它与gt的IOU值大于0.6就为positive.如果和所有的gt都小于0.45，就是negtative，处于0.45~0.6之间的不做处理。

• 损失函数中的α = 1.5；β = 1。

3.2 Data Augmentation

代码复现

这里下载的是pytorch版本的代码，新手学习可能比较快一点。后面就是记录一步一步怎么从一点不懂到把代码运行起来。

README.md

阅读程序第一步都是先找到这个文件，看看环境和文件架构等。voxelnet的pytorch版本的代码的需求如下：

# Dependencies
- `python3.5+`
- `Pytorch` (tested on 0.4.1)
- `TensorBoardX` (tested on 1.4)
- `OpenCV`
- `Pillow` (for add_image in TensorBoardX)
- `Boost` (for compiling evaluation code)
1
2
3
4
5
6
7

我们都采用 pip install *** -i https://pypi.tuna.tsinghua.edu.cn/simple加上清华源就快很多了。

运行一些pre文件

按照readme.md文件的指示，第一步先执行如下程序，编译box_overlaps的Cython模块：

python utils/setup.py build_ext --inplace
1

接着编译nms模块：

cd eval/KITTI
g++ -I path/to/boost/include -o evaluate_object_3d_offline evaluate_object_3d_offline.cpp
1
2

接着给eval文件一个可执行权限:

 cd eval/KITTI
chmod +x launch_test.sh
1
2

数据集准备

从KITTI网页一共需要下载4个zip文件：

• Velodyne point clouds (29 GB)
• Training labels of object data set (5 MB) for input labels of VoxelNet.
• Camera calibration matrices of object data set (16 MB) for visualization of predictions
• Left color images of object data set (12 GB) for visualization of predictions

文件下载下来后，源代码中如下图表示，我们是放在KITTI下面，不是MD_KITTI，KITTI后面的结构和MD_KITTI也不一样。我们此时要组织的结构就是后面这个KITTI的。这里的MD_KITTI后续程序会自行安排。这里先不用管。

数据预处理

下载下来的数据，按照原文的说法是不能都直接用的，是需要进行裁剪的。源码中裁剪的方式大致是如下：

第一步是选择具有反射率大于0的点，也就是第四维r值大于0.然后将其置为1，这样就可以使得点云数据直接和相机投影矩阵相乘。方便后续通过image的坐标来除去多余的点云。（2）通过image的坐标参数，只把在相机前面的点（z>0）的点投影到二维上面。投影的方式是直接乘上一个从calib中得到的参数。（3）通过Image的大小来投影的2维图片来对点云进行裁剪，裁剪的同时加入了很多新的信息，包括了color（这个通过投影图片处理得到），归一化后的投影图片等。。（具体的还需要再研究）

下面为插入内容

对上面的内容进行更新一下：（1）我们下载的点云投影到相机平面的数据是calib_velo_to_cam.txt，表示的是点云到相机的定位文件。在KITTI中还有文件calib_cam_to_cam.txt(相机到相机的标定)。（2）相机和点云的坐标定义：相机（x：右，y：下，z：前） 。点云（x：前，y：左，z：上），也就证实了上文中投影在相机前面的点时采用z>0为判定条件。（3）计算点云到图像的投影矩阵，如下展开说：

参考这篇博客：

1核心思想

计算点云到图像的投影矩阵需要三个参数，分别是P_rect(相机内参矩阵)和R_rect(参考相机0到相机xx图像平面的旋转矩阵)以及Tr_velo_to_cam(点云到相机的[R T]外参矩阵)。

2计算投影矩阵（matlab）参考博客

% 计算点云到图像平面的投影矩阵
R_cam_to_rect = eye(4);
R_cam_to_rect(1:3,1:3) = calib.R_rect{1};  % 参考相机0到相机xx图像平面的旋转矩阵
P_velo_to_img = calib.P_rect{cam+1}*R_cam_to_rect*Tr_velo_to_cam; % 内外参数 
1
2
3
4

3点云投影到图像

投影矩阵乘以点云坐标。在此之前需要把点云填充到四维的齐次坐标，也就是加1。原博客是把前三维作为输入因此需要转化为齐次坐标。投影矩阵是4* 4的。
T*p2_in’

插入内容结束

预处理的下面就只执行该文件，执行前，先修改该文件下的路径信息为你自己的路径信息。
执行：

python crop.py
1

这样就会输出裁剪后的点云了，并且是覆盖在之前的点云文件中。如下：

接下来就是执行生成MD_KITTI文件夹了，在这之前需要去下载一个分数据的 文件；并且把路径放置如下面图的59行和60行的意思。当然要换成你自己的路径。分割的数据在刚刚crop的training 的数据上进行的。在执行之前，需要把文件夹preproc下的split.py文件下的文件路径改成你自己的路径，如下：

然后再运行split.py就可以了

python split.py
1

MD_KITTI文件就此生成了。

讲一下分数据的方法。下载分割数据txt文件中随机选了很多数据名称，分割就按照这个名字分为训练数据和训练时的测试数据。

train

作者在这里说，这个代码只能单个GPU跑。在train之前先去config.py中改一下对应的文件路径和gpu的使用。下面就可以跑程序了，对于car的损失函数的值，在原文中的设置如下：

python train.py --alpha 1.5 --beta 1
1

自己对代码的理解

自己是个新手，以前很少用python，也不知道理解的是否正确，有问题的话希望可以只纠正指出。主要把对代码的初步理解分为下面几个部分：

• 数据在模型中的流动
• 3D数据2D数据的转换
• anchors和gt
• nms操作

数据在模型中的流动。

数据的数据在网络中流动，按照每一个层的方式进行记录

1 Feature learning layer

 网络的输入是$batch\times N\times 4$的数据格式，即是batch（代码中设置为1，即每次值输入一个Bin文件）,一个bin文件中的点云数量N、每个点云的维度$(x,y,z,r)$。
我们知道该层的处理是首先分区，输入点云经过裁剪后的$D\times W\times H=4\times 80\times 70.4$，实际上的范围可以在config.py文件中看到：$X轴：(0,70.4)；Y轴：(-4.,40)；Z轴：(-3,1)$；然后经过partion，代码中把每一个voxel的最大的点云量设置为$T=35$；所以我们假设一共有$k$个非空的体素，经过这样处理后我们将点云表示为$K\times T\times 4$。
接下来就是对全局信息和局部信息的综合。通过求每一个局部voxel的点云的均值，再通过残差的方式加入到每个点的信息中，即会增加三个局部信息的维度：$K\times T\times 7$。
下面就是输入到本文最重要的VFE结构了。这里的数据变换如下：

原始数据为$[B,K,T,7]$；在生成大小为$[B,K,T,128]$的pointwise_feature经过了两个VFE层，在原文中说只用学习一半的参数，对于我这种新手来说就很难理解，看了代码才明白，实际上在每一层都有一个最大池化的cat操作。VEF层的定义如下：

        self.dense = nn.Sequential(nn.Linear(self.in_channels, self.out_channels/2), nn.ReLU())#
        self.batch_norm = nn.BatchNorm1d(self.units)
1
2

仅仅看上面的我们知道输出的通道值仅仅为目标输出值的一半，在后续的forward()函数中，有这样一句操作：

concatenated = torch.cat([pointwise, repeated], dim = 2)#整合了最大值结构
1

也就是如此操作的：

假设输入的通道数为7，目标通道数为32，VFE先学习$7\times 16$的矩阵变换，此时的数据shape变为$[K,T,16]$（这里数据为什么少了一维，后面会说到）;然后对第三个维度$dim=1$的temp，也就是每一个特征维都找最大值。最后组合一下这两个，生成shape为$[K,T,32]$的输出

然后VFE_1，VFE_2的定义：

self.vfe1 = VFELayer(7, 32)
self.vfe2 = VFELayer(32, 128)
1
2

经过了VFE_2层后shape变为$[B,K,T，128]$每一个点都是含有128维的。这里需要说明的是在这一份pythorch的实现代码中是没有带Batch维度的，而是采用的for循环。原因是源代码在tf下是可以在每个batch的点数目不一致的时候使用，而pytorch不行。

在FeatureNet的forward函数中传入的feature大小是$[B,K,T,7]$的原始数据，经过了torch.cat(dim=0)后减少一个维度feature变为$[K,T,7]$(本来B=1);然后feature再取最后一个维度的最大值vmax$[K,T,1]（为了产生mask，具体作用后面研究）$，在经过了刚刚的两个VFE后feature变为为$[K,T，128]$；voxelwise在feature第二维度（dim=1）上取最大值也就是voxelwise的shape为$[K,1,128]$。在pytorch中也就是$[K,128]$

下面这句是这一层的最后一步，首先前面的coordinate（$shape[K,4]$，这个表示的voxel的坐标，需要进一步研究），voxelwise的shape为$[K,128]$，这两个参数也正正好好的表达了每一个voxel的坐标和特征维度。（具体这么用四维来表达的一个voxel的坐标，需要进一步研究）；最后将特征和坐标表达为稀疏表示。

outputs = torch.sparse.FloatTensor(coordinate.t(), voxelwise, torch.Size([batch_size, cfg.INPUT_DEPTH, cfg.INPUT_HEIGHT, cfg.INPUT_WIDTH, 128]))
1

最后得到的outpu的shape为$[batc{h}_{s}ize,INPUT-DEPTH,INPUT-HEIGHT,INPUT-WIDTH,128]$,大小也就是对应的体素总个数的shape为：$[1,10,400,352,128]$,其中最后一维的含义为特征维度。上一句代码的含义就是说：每一个voxel的坐标的地方都有一个128维度的特征向量。需要指出的是原文说了会把剩下的空voxel的特征补上0.

总结一下这一层，这里做的操作就是对每一个voxel的特征进行提取，但是为了不操作空的呢，就只对非空的进行提取；同时VFE层采用了局部和全局坐标结合的方式特取特征；同时又为了解决有的voxel的点不足T个，所有最后采用最大池化表示了一个voxel的特征。

2 Middle layer

这一层在文中仅仅是一个概述，实际上呢，也是唯一一个用到3D卷积的地方了。也十分简单。
上一层的输出为：$[1,10,400,352,128]$，代码为：

        self.middle_layer = nn.Sequential(ConvMD(3, 128, 64, 3, (2, 1, 1,), (1, 1, 1)),
                                          ConvMD(3, 64, 64, 3, (1, 1, 1), (0, 1, 1)),
                                          ConvMD(3, 64, 64, 3, (2, 1, 1), (1, 1, 1)))
1
2
3

这里的第一个参数3表示的3D卷积的含义，原文有进行封装处理；所以这一层的shape变化为：
$[1,10,400,352,128]首先交换一下维度的次序为[1,128,10,400,352]经过第一个3D卷积后为[1,64,5,400,352]$，这里说明一下交换后的通道表示的含义为$(B,D,H,W,C)->(B,C,D,H,W)$深度这一维度是有在减少。（2）经过第二个卷积后变成$[1,64,4,400,352]$，因为其padding在深度上为0；（3）经过第三个3D卷积后为$[1,64,2,400,352]$。

总结一下，这一层用到了3D卷积，没有其他的了。

3 RPN

代码中RPN和Middle层是放在一起的。代码不长，如下：

            self.block1 = nn.Sequential(ConvMD(2, 128, 128, 3, (2, 2), (1, 1)),
                                        ConvMD(2, 128, 128, 3, (1, 1), (1, 1)),
                                        ConvMD(2, 128, 128, 3, (1, 1), (1, 1)),
                                        ConvMD(2, 128, 128, 3, (1, 1), (1, 1)))
            self.deconv1 = Deconv2D(128, 256, 3, (1, 1), (1, 1))                  
            self.block2 = nn.Sequential(ConvMD(2, 128, 128, 3, (2, 2), (1, 1)),
                                    ConvMD(2, 128, 128, 3, (1, 1), (1, 1)),
                                    ConvMD(2, 128, 128, 3, (1, 1), (1, 1)),
                                    ConvMD(2, 128, 128, 3, (1, 1), (1, 1)),
                                    ConvMD(2, 128, 128, 3, (1, 1), (1, 1)),
                                    ConvMD(2, 128, 128, 3, (1, 1), (1, 1)))

        self.deconv2 = Deconv2D(128, 256, 2, (2, 2), (0, 0))

        self.block3 = nn.Sequential(ConvMD(2, 128, 256, 3, (2, 2), (1, 1)),
                                    ConvMD(2, 256, 256, 3, (1, 1), (1, 1)),
                                    ConvMD(2, 256, 256, 3, (1, 1), (1, 1)),
                                    ConvMD(2, 256, 256, 3, (1, 1), (1, 1)),
                                    ConvMD(2, 256, 256, 3, (1, 1), (1, 1)),
                                    ConvMD(2, 256, 256, 3, (1, 1), (1, 1)))

        self.deconv3 = Deconv2D(256, 256, 4, (4, 4), (0, 0))

        self.prob_conv = ConvMD(2, 768, 2, 1, (1, 1), (0, 0), bn = False, activation = False)

        self.reg_conv = ConvMD(2, 768, 14, 1, (1, 1), (0, 0), bn = False, activation = False)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

这里顺便再把原文的RPN贴出来，对比一看就很清楚了

需要注意的是：这份代码貌似第一次卷积的时候多了一层，应该是我上面文章的fliter个数，代码中比原文多一层。
这里没什么好说的，数据在Middle layer最后输出的大小为：
$[batch,64,2,400,352]$,为了可以用二维的卷积，我们把D维和特征维度融合（这里的D为什么是2，是为了和pos,neg对应起来吗？），变成：$[batch,128,400,352]$;
经过block1输出为$[batch,128,200,176]$;
经过block2输出为$[batch,128,100,88]$;
经过block3输出为$[batch,128,50,44]$;
这里的反卷积操作：
第一个block的输出$[batch,128,200,176]$经过deconv1后输出为$[batch,256,200,176]$；
第二个block的输出$[batch,128,100,88]$经过deconv2后输出为$[batch,256,200,176]$；
第三个block的输出$[batch,256,50,44]$经过deconv3后输出为$[batch,256,200,176]$；
也就是上图所示的那样，但是呢。原文图上和代码的deconv1大小不一样，论文应该有误。
组合块$[1,768,200,176]$

接下来这一步就是在对最后得到的组合块$[1,768,200,176]$进行压缩为cls和reg map:
直接暴力压缩到一个通道数为2$[batch,2,200,176]$，另外一个通道数为14$[batch,14,200,176]$。（这里为什么是2和14这个数字呢，我试着解释一下，我们知道最初的voxel个数是$400\times 352$个，我们的anchor是根据这个数目来的，采样的频率是2，也就是在x和y轴上每隔一个单位选择一个为anchor，所以个数是$200\times 176$个，每一个anchor具有七个特征属性（这里考虑0度和90度的旋转的情形），那么每一个anchor就存在14个特征维度，每一个角度分配一组特征）。

anchors和gt

1生成anchors

def cal_anchors():
    # Output:
    # Anchors: (w, l, 2, 7) x y z h w l r
    x = np.linspace(cfg.X_MIN, cfg.X_MAX, cfg.FEATURE_WIDTH)
    y = np.linspace(cfg.Y_MIN, cfg.Y_MAX, cfg.FEATURE_HEIGHT)
    cx, cy = np.meshgrid(x, y)
    # All are (w, l, 2)
    cx = np.tile(cx[..., np.newaxis], 2)
    cy = np.tile(cy[..., np.newaxis], 2)
    cz = np.ones_like(cx) * cfg.ANCHOR_Z
    w = np.ones_like(cx) * cfg.ANCHOR_W
    l = np.ones_like(cx) * cfg.ANCHOR_L
    h = np.ones_like(cx) * cfg.ANCHOR_H
    r = np.ones_like(cx)
    r[..., 0] = 0  # 0
    r[..., 1] = 90 / 180 * np.pi  # 90

    # 7 * (w, l, 2) -> (w, l, 2, 7)
    anchors = np.stack([cx, cy, cz, h, w, l, r], axis = -1)
    return anchors
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

生成的代码如上，我的理解如下:

这里需要配合config.py中的参数一起看，就很好的解释了。np.meshgrid(x, y)生成格网坐标；np.tile(cx[..., np.newaxis], 2)这一句是为了得到两个一样的cx，是不是因为后续的0度和90度各需要一组？暂且这么理解。后续的内容也是每个有两组坐标。这里的z坐标是根据anchor的height来的，使用了先验知识：车都是停在路面上的。所以一开始都只是在x和y轴上进行了采样，z轴一般是固定为车心的高度。同样使用了先验的还有旋转轴，本来可以使yaw,roll和pitch；但是车只能在路面上yaw轴旋转；因此也就少了两个维度。

  anchors = np.stack([cx, cy, cz, h, w, l, r], axis = -1)
1

这里的anchor每一个维度的含义分别是：坐标位置（x,y,z）；文章设置的anchor_size（h,w,l）;可动角度yaw ( r )；这和在VFE中的七个特征维度没有半毛钱的关系。
这里的np.stack()之前各个特征的大小都应该是$200\times 176\times 2$(这里直接把在x和y上的采得到的数据写上了)；但是除r外，每个特征的两组数据都是一样的；经过np.stack()我们的数据会分成两个组，前后两组数据最后一个特征不同，其余的相同；形如：

[x,y,z,,h,w,l,0]和[x,y,z,,h,w,l,pi/2]
1

就此，anchors就算生成成功了，其对应的shape为$[200,176,2,7]$

gt(laels)

每个训练的点云文件都有自己的labels，大小不一，有的是5个，有的有17个。其shape为[bitch,gt_box_num];这里的gt_box_num大小不一，这也是batch在pytorch只能设置为1，其对应的gt_box_num不固定。每一个gt_box都具有15个维度。第一个维度表征是对谁的label（可以是car，也可以是pes等等），剩下的14维是特征。

因此首要的任务是对gtbox的特征转化，在代码中采用label_to_gt_box3d()函数进行操作，具体流程如下：

（1）确定是表征哪一类的gt_box ；对于car的gtbox先转化为carmer_box
（2）转化为carmer_box：14维特征后面七个维度是carme_box的，取出来就可以了
（3）carmer_box到lidar_box:使用函数camera_to_lidar_box；具体是对（x,y,z）的变换是先扩充到4维度再乘以变换矩阵（上文补充资料中讲过）得到在lidar的坐标（x,y,z）；对旋转度r:$r=-r_{camera}-pi/2，这里用角度限制在0到360$;其余三个l,h,w不需要变化。

gt(laels)和anchors之间做loss

既然都得到了3DBBOX，那么就可以按照损失函数求对应的损失了。先给出损失原文函数：
$$\begin{array}{rl}L& =\alpha \frac{1}{N_{\mathrm{p}\mathrm{o}\mathrm{s}}}\sum _i{L}_{\mathrm{c}\mathrm{l}\mathrm{s}}\left(p_i^{\mathrm{p}\mathrm{o}\mathrm{s}},1\right)+\beta \frac{1}{N_{\mathrm{n}\mathrm{e}\mathrm{g}}}\sum _j{L}_{\mathrm{c}\mathrm{l}\mathrm{s}}\left(p_j^{\mathrm{n}\mathrm{e}\mathrm{g}},0\right)\\ & +\frac{1}{N_{\mathrm{p}\mathrm{o}\mathrm{s}}}\sum _i{L}_{\mathrm{r}\mathrm{e}\mathrm{g}}\left({\mathbf{u}}_i,{\mathbf{u}}_i^{\ast }\right)\end{array}$$
代码中是这样求损失的：
总的来说，使用的是函数cal_rpn_target(),这个函数输入的是标签labesls（[bitch,gt_box_num]）和anchors(200,176, 2, 7)，然后把他们都分别转化到鸟视图上，进行RPN操作。这个函数的输出是shape为[batch,200，176,2]的pos_equal_one和[batch,200，176,2]的neg_equal_one以及表示回归的(batch, 200，176, 14)的位置回归。后续会解释为什么是这个shape。那么具体的操作如下：（这里假设已经把labes转化为了gt_box了）
首先是把3维度转化到二维上来：

（1）把$200\times 176\times 2$个bbox合成为二维的，即是：$70400\times 7$的大小，而gtbox大小为：$2\times 7$；
（2）计算需要的中间变量：$d_a$、
（3）将anchors的3DBbox转化到2DBbox中，使用函数anchor_standup_box2d()：第一步：只需要用到7维特征中的$(x,y,w,l)$四个维度,因此将其取出；第二步：根据长方形中心点和长宽的关系得到左上和右下的角点坐标；
（4）将gtbox的3Dbox转化到2Dbox中。会使用到$(x,y,w,l,r)$ 5个特征维度，这是因为gtbox的旋转度要对应到anchor的Box当中来，也就是一个gt_box要分成两个对应的feature_box，也就是$(labelnum,5,1)->(labelnum,4,2)$的维度变化。第一步：取出5个特征维度先变成正常的3Dbox$(labelnum,7,1)$(那5个维度以外的数据用0填补就可以了)。第二步：$(labelnum,7)->(labelnum,8,3)$;这一步有超级多的变换，没有深入研究。第三步：$(labelnum,7)->(labelnum,4,2)$;这里的操作就是直接取(labelnum,8, 3)里面的元素就可以了。

此时已经分别得到了anchor的二维box和gt的二维box,下面的步骤就是在二维上对计算的IOU进行一些筛选得到pos的anchor：

（1）使用函数bbox_overlaps得到IOU矩阵，shape为74600*gt_num，第一维度对应到一共有多少个anchor，第二个维度表示gt的个数。
（2）找到每一个gt对应的IOU最大的anchor的编号id_highest
（3）在[0,gt_num]均匀采样出id_highest_gt[0,1,2,...,gt_num-1]
（4）根据id_highest和id_highest_gt创建mask[1,gt_num]；
（5）上面求得了第一种的pos_anchor，文章中说过如果anchor对任意一个gt的IOU大于0.6，那么也算是pos。这一步就是选取那些大于0.6的anchor。得到id_pos, id_pos_gt的组合；再和前面IOU最高的组合成id_pos,id_pos_gt。再除去那些既是最大IOU又大于0.6的anchor。
（6）既然求了pos，也要求neg，对于小于0.45的算是neg。得到id_neg，这里只有neg的编号，因为gt是爸爸，没有neg。
（7）接下来就要在对应的3D坐标中找到数哪些anchor是pos的，采用函数：
index_x, index_y, index_z = np.unravel_index(id_pos, (*feature_map_shape, 2)) pos_equal_one[batch_id, index_x, index_y, index_z] = 1 这里的两句就是把对应Pos的anchor变为1.
（8）下面一步就是根据上面的数据求对应特征的偏差，也就是：
$$\begin{array}{l}\Delta x=\frac{x_c^g-x_c^a}{d^a},\Delta y=\frac{y_c^g-y_c^a}{d^a},\Delta z=\frac{z_c^g-z_c^a}{h^a}\\ \Delta l=\log \left(\frac{l^g}{l^a}\right),\Delta w=\log \left(\frac{w^g}{w^a}\right),\Delta h=\log \left(\frac{h^g}{h^a}\right)\\ \Delta \theta ={\theta }^g-{\theta }^a\end{array}$$
这一堆东西。储存在变量targets中，shape为[1,200,176,14]
（9）对neg的就不要求上面的残差了，只需要知道体素中哪个位置是neg的。
（10）这是后就需要想，我们怎么用得到的得到的pos进行对pos的anchor位置回归。（这个后续再说）

那我们如何根据上面的参数计算loss呢，如下：

        cls_pos_loss = (-pos_equal_one * torch.log(prob_output + small_addon_for_BCE)) / pos_equal_one_sum
        cls_neg_loss = (-neg_equal_one * torch.log(1 - prob_output + small_addon_for_BCE)) / neg_equal_one_sum

        cls_loss = torch.sum(self.alpha * cls_pos_loss + self.beta * cls_neg_loss)
        cls_pos_loss_rec = torch.sum(cls_pos_loss)
        cls_neg_loss_rec = torch.sum(cls_neg_loss)

        reg_loss = smooth_l1(delta_output * pos_equal_one_for_reg, targets * pos_equal_one_for_reg, self.sigma) / pos_equal_one_sum
        reg_loss = torch.sum(reg_loss)

        loss = cls_loss + reg_loss
1
2
3
4
5
6
7
8
9
10
11

（1）这里就引来了RPN的输出的prob_output和delta_output用于计算损失了，这里的prob_output 是在sigmod输出的分数后的值。现在填一下前面说的输出的delta_output怎么去回归每个anchor，在reg_loss中的定义我们看到，传进去的三个参数：delta_output * pos_equal_one_for_reg、targets * pos_equal_one_for_reg、self.sigma。这里的pos_equal_one_for_reg只有1和0两种元素，前两个最重要的参数的含义分别是：输出的预测残差和实际的残差值。这里我们知道targets 是在三维中预测的pos anchor对gt的七个维度的残差值。也就是说回归实际上是对这个残差的拟合。

这里就基本把anchor，gt和Loss讲清楚了。后续有时间再研究一下前面的坑。后续一段时间可能要去搞一下ros了。一周后再做3D检测。