一个 c++版本的 yolo 部署 （四）-- 基于cuda的预处理

重新理解预处理任务

在上一篇中，已经粗略的实现了一个基于opencv的cpu版本预处理，但是里面很多细节都是一带而过的，因为opencv本身封装了许多细节，而这些细节是我们在cuda程序中必须自己完成的。
我们待处理的原图可能是任意尺寸的，比如1280x640，但是网络的输入是固定的，比如640x640，这就是我们要对图像进行缩放的原因。一般来说我们训练和推理用的图都会保持宽高比，那么显然，我们会对原图宽高各缩1/2，变为640x320，然后居中显示：

这里一共有两步：

1. 将图片高宽各缩0.5
2. 将缩小后的图片左上角移到（0，160）

仿射变换 Affine Transform

只要学过初中数学，缩放+平移的计算就一目了然：
$$\{\begin{array}{l}{x}^{\prime }=\mathbf{x}\ast 0.5+0\\ y^{\prime }=\mathbf{y}\ast 0.5+160\end{array}$$
仔细观察这组公式，会发现粗体部分（缩放）是线性变换，而斜体部分（平移）是 非线性变换。
线性代数告诉我们，旋转，镜像，缩放，都是可以表示为对应矩阵的。比如长宽各缩0.5的缩放变换可以用矩阵运算表示：
$$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right]=\left[\begin{array}{cc}0.5& 0\\ 0& 0.5\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]$$
其他变换矩阵可以自己搜一下，这些都是计算机图形学基础，比如在3D游戏中，一个人物移动，就是大量的三角形进行各类变换。

但是刚刚提到位移是非线性变换，所以为了能同时进行线性变换+平移操作，我们需要了解一个新概念：齐次项，通过补充齐次项，变换矩阵变为：
$$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{ccc}0.5& 0& 0\\ 0& 0.5& 160\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$$
这个新的变换就叫仿射变换 Affine Transform

双线性插值缩放算法

假设我要将任意图缩放为640x640的目标尺寸，实际上我要完成的任务到底是什么？一个非常经典的想法是，将目标尺寸看成是640x640个网格，我真正要做的是计算出每个网格的颜色。当每个网格都填上正确的颜色后，自然一张图片就缩放完成了。

于是问题就变成了如何计算一个网格的颜色。如果我能计算出这个格子相对于原图的位置，那么就能根据其周围4个点的颜色，以及距离权重，加权计算出这个点的实际颜色

权重在高和宽上都是线性变化的，所以这个经典算法就是双线性插值算法。

理解了这个算法的本质后，我们会发现，这个算法其实就是进行640x640次简单的计算，这是一个天然适合用并行计算加速的任务。并且，在计算完每个格子的颜色后，归一化和通道转置都可以顺便完成，前者就是对计算出来的颜色除255，后者就是将顺序存储的3个通道颜色换个存储位置。

现在唯一的问题是，我上面只介绍了从原图坐标到目标图坐标的转换，如何计算目标图坐标到原图坐标的转换呢？

逆矩阵

从目标图到原图的转换当然也可以自己尝试推导，只是坐标偏移（齐次项）稍微需要费点脑细胞。但是线性代数对这类问题（不只是缩放，对其他各类线性变化都适用）有更标准的解决方案：逆矩阵。

若存在一个变换矩阵A，使得
$$\left[P^{\prime }\right]=\left[A\right]\left[P\right]$$
那其逆矩阵（若存在）$A^{-1}$将满足以下性质：
$$\left[A^{-1}\right]\left[P^{\prime }\right]=\left[A^{-1}\right]\left[A\right]\left[P\right]=\left[I\right]\left[P\right]=\left[P\right]$$
其中$I$是单位矩阵，一个矩阵和其逆矩阵相乘得到的就是单位矩阵，任何矩阵和单位矩阵相乘都是其本身。

计算逆矩阵的方法，线性代数中也有专门的研究，比如高斯消元法。

更多的内容并不是本文想讲述的重点，读者可以自行学习相关内容。opencv中有专门的函数获取一个矩阵的逆矩阵，说这么多也只是为了解释以下这一行代码：

	// 获取逆变换矩阵，即目标图坐标到原图的变换矩阵A'
    cv::invertAffineTransform(s2d_mat, d2s_mat);

cuda代码

唠唠叨叨说了这么多，我们知道了预处理实际需要做哪些事情：
我们需要进行640x640次计算任务，每个任务：

• 计算自己负责的点映射到原图的坐标，并找到相邻的4个点
• 根据距离加权计算自己的rgb颜色值
• 对计算结果进行归一化
• 将结果存储到合适的位置

于是有了这一段cuda代码(有简化，完整内容可到这里查看）

__global__ void preprocess_kernel(uint8_t* img_buffer_device, float* input_buff_device, int src_w, int src_h, int dst_w, int dst_h, float* d2s_matrix, int edge){
    int position = blockDim.x * blockIdx.x + threadIdx.x;
    if(position >= edge) return;    // 因为cuda并行计算的特性，可能有部分线程是多余的，不需要实际参与计算
    int dx = position % dst_w;
    int dy = position / dst_w;
    uint8_t default_value = 128;


	float sx = d2s_matrix[0] * dx + d2s_matrix[1] * dy + d2s_matrix[2];
	float sy = d2s_matrix[3] * dx + d2s_matrix[4] * dy + d2s_matrix[5];

    
    int channel = 3;    // cv::Mat应该都是3通道的
    float c0,c1,c2;
	bool done = false;
    // 有效范围为实际范围外延一个像素。对于不在有效范围内的点，直接填充默认值（128）
    if(sx <= -1  || sx >= src_w || sy <= -1 || sy >= src_h){
        c0 = c1 = c2 = default_value;
		done = true;
	} 
	if(!done) {
        int y_low = sy;
        int x_low = sx;
        int y_high = y_low + 1;
        int x_high = x_low + 1;

        // 每个目标点的实际颜色是周围四个点的加权值之和
        float x_h_weight = sx - x_low;
        float y_h_weight = sy - y_low;
        float x_l_weight = 1 - x_h_weight;
        float y_l_weight = 1 - y_h_weight;

        uint8_t default_point[3] = {default_value, default_value, default_value};
        uint8_t* left_top = default_point;
        uint8_t* right_top = default_point;
        uint8_t* left_bottom = default_point;
        uint8_t* right_bottom = default_point;
        if(y_low >= 0){
            if(x_low >= 0){
                left_top = img_buffer_device + (y_low * src_w + x_low) * channel;
            }
            if (x_high < src_w){
                right_top = img_buffer_device + (y_low * src_w + x_high) * channel;
            }
        }
        if(y_high < src_h){
            if(x_low >= 0){
                left_bottom = img_buffer_device + (y_high * src_w + x_low) * channel;
            }
            if(x_high < src_w){
                right_bottom = img_buffer_device + (y_high * src_w + x_high) * channel;
            }
        }
        // cv::Mat的通道顺序是BGR，而yolo的输入是RGB，所以需要交换
        c2 = left_top[0] * x_l_weight * y_l_weight + right_top[0] * x_h_weight * y_l_weight + left_bottom[0] * x_l_weight * y_h_weight + right_bottom[0] * x_h_weight * y_h_weight;
        c1 = left_top[1] * x_l_weight * y_l_weight + right_top[1] * x_h_weight * y_l_weight + left_bottom[1] * x_l_weight * y_h_weight + right_bottom[1] * x_h_weight * y_h_weight;
        c0 = left_top[2] * x_l_weight * y_l_weight + right_top[2] * x_h_weight * y_l_weight + left_bottom[2] * x_l_weight * y_h_weight + right_bottom[2] * x_h_weight * y_h_weight;
    }
    input_buff_device[position] = c0 / 255.0;
    input_buff_device[position + edge] = c1 / 255.0;
    input_buff_device[position + edge * 2] = c2 / 255.0;
}

void preprocess_by_cuda(unsigned char* img, float* device_input, int dst_h, int dst_w, float* d2s_matrix, cudaStream_t& stream){
    int src_w = img.cols;
    int src_h = img.rows;
    // 分配cuda线程
    int jobs = dst_h * dst_w;  // 每个像素给线程，负责完成该像素的生成、rgb拆解、归一化操作
    int threads = 256;
    int blocks = ceil(jobs / (float)threads);

    preprocess_kernel<<<blocks, threads, 0, stream>>>(img, device_input, src_w, src_h, dst_w, dst_h, d2s_matrix, jobs);

    CUDA_CHECK(cudaStreamSynchronize(stream));
}

接下来会写cuda后处理
to be continue…