Halide编程

免责申明：本文所用的材料来自于网络，仅供个人学习和交流之用。如果用作其他用途，那么请自行承担法律负责。

Halide是一种特定领域的编程语言（DSL），专门用于加速图像处理，当前已支持x86 SSE，ARM v7 NEON，CUDA，OpenCL等。Halide的主要特点在于分开算法描述和调度描述。一个算法的Halide算法描述不依赖于硬件，其调度描述使用并行化、向量化（SIMD指令）和处理本地化（block）等函数（为便于表述，本文使用了原语），通过这些原语的某种特定组合，使得算法描述可在选定的硬件计算平台上取得最佳调度，或者说高效地执行程序里的计算任务。

Halide程序的执行流程如下：

1. C++代码的编写使用Halide提供的类和函数，如Halide :: Func、Halide :: Expr和Halide :: Var等。
2. 用C++编译器生成目标文件。
3. 链接到Halide的lib文件（libHalide.a），执行用Halide :: Func、Halide :: Expr、Halide :: Var等编写的代码。
4. 首次执行Halide代码时，将动态启动Halide编译器，生成二进制代码，执行图像处理。现有的二进制代码用于第二和后续计算。
   Halide编译器是用C++编写的程序，以lib文件（静态库）和头文件形式提供。当执行主程序时，Halide编译器会动态启动，另一种方法是静态地执行Halide代码（预先编译Halide代码）。

以放大/缩小图像尺寸的过程为例，在解释算法的内容后，本文给出参考代码。
放大/缩小算法的步骤如下所述：

1. 选定任一个放大处理后的像素A；
2. 推算像素A在源图像中的位置，记为B。通常，该位置坐标由两个浮点数（单位距离是一个像素）表示。
3. 位置B上的像素值将关联于其周围像素的像素值，如最近邻法、双线性法和双三次法等。本文的算法实现以双线性插值为例，其中双线性的计算区域由四个子区域组成，记为C0、C1、C2和C3，其中像素B和左上角像素之间的水平和垂直距离分别表示为di和dj。这种计算方式又称为面积法，其中像素B的值通过下式计算。
   V(B) = C0V(P0)+ C1V(P1)+ C2V(P2)+ C3V(P3)
   其中V(x)表示某个像素上的数值，通常灰度图像用一个标量值，彩色图像有三个通道值。注意等式C0+C1+C2+C3=1成立。

出于比较的考虑，先用C++里的指针等实现双线性插值算法，而不用Halide实现该算法。

#include "normal.hpp"
#include <iostream>
#include <opencv2/opencv.hpp>

inline const uint8_t* get_pixel(const uint8_t* data, int row, int col, size_t step, size_t elem_size)
{
   
    return data + (row * step + col * elem_size);
}

inline uint8_t interpolate(float c0, float c1, float c2, float c3, const uint8_t* p0, const uint8_t* p1, const uint8_t* p2, const uint8_t* p3, int i)
{
   
    return c0 * p0[i] + c1 * p1[i] + c2 * p2[i] + c3 * p3[i];
}

void resize_with_raw_access(const cv::Mat& src_image, cv::Mat& dst_image)
{
   
    //使用OpenCV 读取和保存图像。

    // 返回输入源图像的参数
    const int src_cols = src_image.cols;
    const int src_rows = src_image.rows;
    const size_t src_step = src_image.step;
    const size_t src_elem_size = src_image.elemSize();
    const uint8_t* src_data = src_image.data;

    // 返回计算结果的目标图像的参数
    const int dst_cols = dst_image.cols;
    const int dst_rows = dst_image.rows;
    const size_t dst_step = dst_image.step;
    const size_t dst_elem_size = dst_image.elemSize();
    uint8_t* dst_data = dst_image.data;
    // 返回图像拉伸和缩小的比例参数
    const float sc = static_cast<float>(src_cols) / dst_cols;
    const float sr = static_cast<float>(src_rows) / dst_rows;
    // 由目标图像的左上角像素开始图像伸缩操作
    // 遍历y方向，纵向遍历图像的各行
    for (auto j = 0, jd = 0; j < dst_rows; ++j, jd += dst_step)
    {
   
        // 返回该像素在源图像里的所在行
        const float fj = j * sr;
        const int cj0 = static_cast<int>(fj); // 整数部分。
        const float dj = fj - cj0;	//小数部分
        // 整数值部分+1，限定像素在源图像里的位置不超高源图像的行数。
        const int cj1 = cj0 + 1 >= src_rows ? cj0 : cj0 + 1;
        // 返回该像素在目标图像里的位置。
        uint8_t* dst_p = dst_data + jd;
        // 遍历x方向
        for (auto i = 0, id = 0; i < dst_cols; ++i, id += dst_elem_size)
        {
   
            // 返回每个像素在源图像里的位置
            const float fi = i * sc;
            const int ci0 = static_cast<int>(fi);//整数部分
            const float di = fi - ci0; // dx 小数部分
            // 限定像素在源图像的读取位置
            const int ci1 = ci0 + 1 >= src_cols ? ci0 : ci0 + 1;
            // 计算四个子区域的面积
            const