Warp算法OpenCL异构开发总结_网站优化分享

一、OpenCL简介

OpenCL最初由苹果公司(Apple)提出(其他合作公司有AMD，IBM，Qualcomm(高通)，Intel和NVIDIA)，之后交由非盈利组织Khronos维护。最初的1.0版标准，由Khronos在2008年发布。OpenCL 1.0定义了主机端的接口，以及使用C语言作为OpenCL内核书写的语言，内核就是在不同的异构设备上并行处理数据的单位。之后的几年，发布了OpenCL 1.1和OpenCL 1.2，新标准为OpenCL增加了很多特性，比如：提高与OpenGL的互动性，补充了很多图像格式，同步事件，设备划分等特性。2013年11月，Khronos组织正式发布了OpenCL 2.0标准。为OpenCL添加了更多的新特性，比如：共享虚拟内存、内核嵌套并行和通用地址空间。这些更加高级的功能会让并行开发变得越来越简单，并且提高了OpenCL应用执行的效率。

开源编程标准设计者也要面对很多的挑战，为了形成一套通用的编程标准，要对一些要求进行一定的取舍。Khronos在这方面做得很不错，其设计的API都能很好的兼容不同的架构，并且能让硬件发挥其最大的性能。只要正确的遵循编程标准，那么一套程序几乎不用做什么修改，就可以从一个硬件平台，移植到另一个硬件平台上。供应商和设备分离的编程模型给OpenCL带来了极佳的可移植性，使其能充分发挥不同平台的加速能力。

执行在OpenCL设备上的代码，与执行在CPU上的不同，其使用OpenCL C进行书写。OpenCl C遵循更加严格的C99标准，在此基础上进行了适当的扩展，使其能在各种异构设备上以数据并行的方式执行。新标准中OpenCL C编程实现了C11标准中的原子操作(其子集)和同步操作。因为OpenCL API本身是C API，那么第三方就将其绑定到很多语言上，比如：Java，C++，Python和.NET。除此之外，很多主流库(线性代数和机器视觉)都集成了OpenCL，为的就是在异构平台上获得实质性的性能提升。

二、OpenCL 标准

OpenCL标准分为四部分，每一部分都用“模型”来定义。

平台模型：指定一个host处理器，用于任务的调度。以及一个或多个device处理器，用于执行OpenCL任务(OpenCL C Kernel)。这里将硬件抽象成了对应的设备(host或device)。

执行模型：定义了OpenCL在host上运行的环境应该如何配置，以及host如何指定设备执行某项工作。这里就包括host运行的环境，host-device交互的机制，以及配置内核时使用到的并发模型。并发模型定义了如何将算法分解成OpenCL工作项和工作组。

内核编程模型：定义了并发模型如何映射到实际物理硬件。

内存模型：定义了内存对象的类型，并且抽象了内存层次，这样内核就不用了解其使用内存的实际架构。其也包括内存排序的要求，并且选择性支持host和device的共享虚拟内存。

通常情况下，OpenCL实现的执行平台包括一个x86 CPU主处理器，和一个GPU设备作为加速器。主处理器会将内核放置在GPU上运行，并且发出指令让GPU按照某个特定的并行方式进行执行。内核使用到的内存数据都由编程者依据层级内存模型分配或开辟。运行时和驱动层会将抽象的内存区域映射到物理内存层面。最后，由GPU开辟硬件线程来对内核进行执行，并且将每个线程映射到对应的硬件单元上。

三、OpenCL基本概念

主机（Host）：CPU及其内存组成的计算系统。

平台（Platform）：主机和OpenCL管理框架下的若干个设备构成了一个平台，所有GPU操作都限定于选择的Platform上运行。OpenCL编程的第一步就是选择并初始化一个平台。

上下文（Context）：定义了整个OpenCL的运行环境，包括Kernel、Device、内存管理和指令队列等。

核函数（Kernel）：是在设备程序上执行运算的入口函数，在主机上调用。

设备（Device）： GPU及其显存组成的计算系统。

指令队列（Command-Queue）：一些需要在设备上执行的OpenCL指令的集合。

SIMT（Single Instruction Multi Thread）：单指令多线程，GPU并行运算的主要方式，很多个多线程同时执行相同的运算指令，当然可能每个线程的数据有所不同，但执行的操作一致。

工作项（Work-item）：跟CUDA中的线程（Threads）是同一个概念，N多个工作项（线程）执行同样的核函数，每个Work-item都有一个唯一固定的ID号，一般通过这个ID号来区分需要处理的数据。

工作组（Work-group）：跟CUDA中的线程块（Block）是同一个概念，N多个工作项组成一个工作组，Work-group内的这些Work-item之间可以通信和协作。

ND-Range ：跟CUDA中的网格是同一个概念，定义了Work-group的组织形式。

Warp算法OpenCL异构开发总结,请添加图片描述,第1张

小结：

OpenCL的整个计算空间称为NDRange空间，它根据计算单元数量分为对应数量的work-group,而一个work-group在一个计算单元上执行，有多少个计算单元，就会有多少个work-group可以并行计算，同时每个work-group中又包含了多个work-item，其数量根据一个work-group可以同时运行多少个work-item决定。

四、基本用法

OpenCL在的编程上，通常是以最细粒度的并行，如对于图像来说，最细的粒度就是像素点，OpenCL可以对每一个像素点同时进行同一个核函数的处理。OpenCL的一般性体现在，接口的通用性和底层内核直接映射物理资源。

通过对比三个不同版本的向量加法例子来说明OpenCL的处理过程：

（1）串行C语言实现的版本

使用串行的方式实现向量加法(C语言)，其使用一个循环对每个元素进行计算。每次循环将两个输入数组对应位置的数相加，然后存储在输出数组中。

void vecadd(int *C, int *A, int *B, int N){
  for (int i = 0; i < N; ++i){
    C[i] = A[i] + B[i];
  }
}

（2）多线程C语言实现的版本

对于一个多核设备，我们要不就使用底层粗粒度线程API(比如，Win32或POSIX线程)，要不就使用数据并行方式(比如，OpenMP)。粗粒度多线程需要对任务进行划分(循环次数)。因为循环的迭代次数特别多，并且每次迭代的任务量很少，这时我们就需要增大循环迭代的粒度，这种技术叫做“条带处理”。

void vecadd(int *C, int *A, int *B, int N, int NP, int tid){
  int ept = N / NP; // 每个线程所要处理的元素个数
  for (int i = tid * ept; i < (tid + 1) * ept; ++i){
    C[i] = A[i] + B[i];
  }
}

（3）OpenCL C实现的版本

OpenCL C上的并发执行单元称为工作项(work-item)。每一个工作项都会执行内核函数体。这里就不用手动的去划分任务，这里将每一次循环操作映射到一个工作项中。OpenCL运行时可以创建很多工作项，其个数可以和输入输出数组的长度相对应，并且工作项在运行时，以一种默认合适的方式映射到底层硬件上(CPU或GPU核)。

概念上，这种方式与并行机制中原有的功能性“映射”操作(可参考mapReduce)和OpenMP中对for循环进行数据并行类似。当OpenCL设备开始执行内核，OpenCL C中提供的内置函数可以让工作项知道自己的编号。

下面的代码中，调用get_global_id(0)来获取当前工作项的位置，以及访问到的数据位于数组中的位置。get_global_id()的参数用于获取指定维度上的工作项编号，其中“0”这个参数，可获取当前第一维上工作项的ID信息。

__kernel
void vecadd(__global int *C, __global int *A, __global int *B){
  int tid = get_global_id(0); // OpenCL intrinsic函数
  c[tid] = A[tid] + B[tid];
}

1.clEnqueueNDRangeKernel的用法

OpenCL采用NDRange(Global Dimemsion Index Ranges)来组织所有work-item进行工作的。

cl_int clEnqueueNDRangeKernel ( cl_command_queue command_queue,
                                cl_kernel kernel,
                                cl_uint work_dim,
                                const size_t *global_work_offset,
                                const size_t *global_work_size,
                                const size_t *local_work_size,
                                cl_uint num_events_in_wait_list,
                                const cl_event *event_wait_list,
                                cl_event *event)

参数：

[command_queue] 一个有效的命令队列。内核将在与command_queue关联的设备上排队执行。

[kernel] 一个有效的内核对象。与kernel和command_queue关联的OpenCL上下文必须是相同的。

[work_dim] 用于指定全局工作项和工作组中的工作项的维数。Work_dim必须大于0且小于或等于3。

[global_work_offset] 当前必须是NULL值。在OpenCL的未来版本中，global_work_offset可以用来指定一个work_dim无符号值数组，用来描述用于计算工作项全局ID的偏移量，而不是让全局ID总是从offset(0,0，…0)．

[global_work_size] 指向一个work_dim无符号值数组，该数组描述work_dim维度中将执行内核函数的**全局工作项的数量。**全局工作项的总数被计算为global_work_size[0] … global_work_size[work_dim - 1]。global_work_size中指定的值不能超过内核执行将在其上排队的设备的sizeof(size_t)给出的范围。设备的sizeof(size_t)可以通过使用clGetDeviceInfo的OpenCL设备查询表中的CL_DEVICE_ADDRESS_BITS来确定。例如，如果CL_DEVICE_ADDRESS_BITS = 32，即设备使用32位地址空间，则size_t是32位无符号整数，global_work_size值必须在1范围内。2 ^ 32 - 1。超出此范围的值将返回CL_OUT_OF_RESOURCES错误。

[local_work_size] 指向work_dim无符号值的数组，该数组描述组成工作组的工作项的数量(也称为工作组的大小)，这些工作项将执行由kernel指定的内核。工作组中工作项的总数计算为local_work_size[0] … local_work_size[work_dim - 1]。工作组中工作项的总数必须小于或等于在OpenCL设备查询clGetDeviceInfo表中指定的CL_DEVICE_MAX_WORK_GROUP_SIZE值和在local_work_size[0]，…local_work_size [work_id]。实际情况由于设备和内核的不同，导致该参数在设置时候不同，要根据硬件配置来设置。尝试clGetKernelWorkGroupInfo()，获取local_work_size

[event] 返回标识此特定内核执行实例的事件对象。事件对象是惟一的，以后可以用来标识特定的内核执行实例。如果event为NULL，则不会为这个内核执行实例创建任何事件，因此应用程序将不可能查询或排队等待这个特定的内核执行实例。一般用于GPU 运算耗时统计。

下面举一个NDRange为二维的例子

Warp算法OpenCL异构开发总结,请添加图片描述,第2张

从该图可知：

每一个小方格代表了一个工作项work-item
每个4×4的中方格代表了一个工作组work-group

整个12×12的大方格代表了整个NDRange空间

该图中可以通过以下代码得到相应的数据：

uint dim=get_work_dim（）；//dim=2
size_t global_id_0=get_global_id(0);//从参数global_offset(0,0)第一个参数0开始，个数为global_size（12,12）的第一参数12
size_t global_id_1=get_global_id(1);//从参数global_offset(0,0）第二个参数0开始，个数为global_size(12,12)的第二个参数12
size_t global_size_0=get_global_size(0);//大小为global_size（12,12）的第一个参数12，即全局空间第一个维度的大小
size_t global_size_1=get_global_size(1);//大小为global_size（12,12）的第二个参数12，即全局空间第二个维度的大小
size_t offset_0=get_global_offset(0);//获取global_offset(0,0)的第一个参数0
size_t offset_1=get_global_offset(1);//获取global_offset(0,0)的第二个参数0
size_t local_id_0=get_local_id(0);//获取local_size(4,4)的第一个参数个数4，即局部空间（组内空间）第一个维度的大小
size_t local_id_1=get_local_id(1);//获取local_size(4,4)的第二个参数个数4，即局部空间（组内空间）第二个维度的大小

2.创建并执行一个简单的OpenCL应用大致需要以下几步：

查询平台和设备信息
创建一个上下文
为每个设备创建一个命令队列
创建一个内存对象(数组)用于存储数据
拷贝输入数据到设备端
使用OpenCL C代码创建并编译出一个程序
从编译好的OpenCL程序中提取内核
执行内核
拷贝输出数据到主机端
释放资源

OpenCL向量相加的例子：

// This program implements a vector addition using OpenCL
 
// System includes
#include 
#include 
#include 
// OpenCL includes
#include 
 
// OpenCL kernel to perform an element-wise addition 
const char* programSource =
"__kernel                                            \n"
"void vecadd(__global int *A,                        \n"
"            __global int *B,                        \n"
"            __global int *C)                        \n"
"{    \n"
"     \n"
"   // Get the work-item’s unique ID                 \n"
"   int idx = get_global_id(0);                      \n"
"     \n"
"   // Add the corresponding locations of            \n"
"   // 'A' and 'B', and store the result in 'C'.     \n"
"   C[idx] = A[idx] + B[idx];                        \n"
"}    \n"
;
 
int main() {
    // This code executes on the OpenCL host
    
    // Host data
    int *A = NULL;  // Input array
    int *B = NULL;  // Input array
    int *C = NULL;  // Output array
    
    // Elements in each array
    const int elements = 2048;   
    
    // Compute the size of the data 
    size_t datasize = sizeof(int)*elements;
 
    // Allocate space for input/output data
    A = (int*)malloc(datasize);
    B = (int*)malloc(datasize);
    C = (int*)malloc(datasize);
 
    // Initialize the input data
    int i;
    for(i = 0; i < elements; i++) {
        A[i] = i;
        B[i] = i;
    }
 
    // Use this to check the output of each API call
    cl_int status;  
     
    // Retrieve the number of platforms
    /* 1. 查询平台和设备 */
    cl_uint numPlatforms = 0;
    status = clGetPlatformIDs(0, NULL, &numPlatforms); // 检索平台的数量
    std::cout << "numPlatforms: " << numPlatforms << std::endl;
    // Allocate enough space for each platform
    cl_platform_id *platforms = NULL;
    platforms = (cl_platform_id*)malloc(
        numPlatforms*sizeof(cl_platform_id)); // 为每个平台对象分配足够的空间
    std::cout << "platforms: " << platforms << std::endl;
 
    // Fill in the platforms
    status = clGetPlatformIDs(numPlatforms, platforms, NULL); // 将具体的平台对象填充其中
 
    // Retrieve the number of devices
    cl_uint numDevices = 0;
    cl_device_id* device_id;
    status = clGetDeviceIDs(platforms[0], CL_DEVICE_TYPE_ALL, 0, 
        NULL, &numDevices); // 检索设备的数量
    std::cout << "numDevices: " << numDevices << std::endl;
    // Allocate enough space for each device
    cl_device_id *devices;
    devices = (cl_device_id*)malloc(
        numDevices*sizeof(cl_device_id)); // 为每个设备对象分配足够的空间
 
    // Fill in the devices 
    status = clGetDeviceIDs(platforms[0], CL_DEVICE_TYPE_ALL,
        numDevices, devices, NULL); // 将具体的设备对象填充其中
  /* 获取当前计算单元的信息 */
    size_t maxWorkItemPerGroup;
    cl_uint maxComputeUnits = 0;
    cl_ulong maxGlobalMemSize = 0;
    cl_ulong maxConstantBufferSize = 0;
    cl_ulong maxLocalMemSize = 0;
    /// print parallel compute units(CU)
    clGetDeviceInfo(*devices, CL_DEVICE_MAX_COMPUTE_UNITS, sizeof(maxComputeUnits),
    &maxComputeUnits, NULL);
    printf("Parallel compute units: %u\n", maxComputeUnits);
    /// maxWorkItemPerGroup
    clGetDeviceInfo(*devices, CL_DEVICE_MAX_WORK_GROUP_SIZE,
    sizeof(maxWorkItemPerGroup), &maxWorkItemPerGroup, NULL);
    printf("maxWorkItemPerGroup: %zd\n", maxWorkItemPerGroup);
    /// print maxGlobalMemSize
    clGetDeviceInfo(*devices, CL_DEVICE_GLOBAL_MEM_SIZE, sizeof(maxGlobalMemSize),
    &maxGlobalMemSize, NULL);
    printf("maxGlobalMemSize: %lu(MB)\n", maxGlobalMemSize / 1024 / 1024);
    /// print maxConstantBufferSize
    clGetDeviceInfo(*devices, CL_DEVICE_MAX_CONSTANT_BUFFER_SIZE,
    sizeof(maxConstantBufferSize), &maxConstantBufferSize, NULL);
    printf("maxConstantBufferSize: %lu(KB)\n", maxConstantBufferSize / 1024);
    /// print maxLocalMemSize
    clGetDeviceInfo(*devices, CL_DEVICE_LOCAL_MEM_SIZE, sizeof(maxLocalMemSize),
    &maxLocalMemSize, NULL);
    printf("maxLocalMemSize: %lu(KB)\n", maxLocalMemSize / 1024);
    // Create a context and associate it with the devices
    /* 2. 创建一个上下文 */
    cl_context context;
    context = clCreateContext(NULL, numDevices, devices, NULL, 
        NULL, &status); // 创建的上下文包含所有找到的设备
 
    // Create a command queue and associate it with the device
    /* 3. 为每个设备创建一个命令队列 */
    cl_command_queue cmdQueue;
    cmdQueue = clCreateCommandQueue(context, devices[0], 0, 
        &status); // 为第一个发现的设备创建命令队列
 
    // Create a buffer object that will contain the data 
    /* 4. 创建设备数组用于存储数据 */
    /* 向量加法的三个向量，2个输入数组和1个输出数组 */
    // from the host array A
    cl_mem bufA;
    bufA = clCreateBuffer(context, CL_MEM_READ_ONLY, datasize,                       
       NULL, &status);
 
    // Create a buffer object that will contain the data 
    // from the host array B
    cl_mem bufB;
    bufB = clCreateBuffer(context, CL_MEM_READ_ONLY, datasize,                        
        NULL, &status);
 
    // Create a buffer object that will hold the output data
    cl_mem bufC;
    bufC = clCreateBuffer(context, CL_MEM_WRITE_ONLY, datasize,
        NULL, &status); 
    /* 5. 拷贝输入数据到设备端 */
    /* 将输入数据填充到数组中 */
    // Write input array A to the device buffer bufferA
    status = clEnqueueWriteBuffer(cmdQueue, bufA, CL_FALSE, 
        0, datasize, A, 0, NULL, NULL);
    
    // Write input array B to the device buffer bufferB
    status = clEnqueueWriteBuffer(cmdQueue, bufB, CL_FALSE, 
        0, datasize, B, 0, NULL, NULL);
 
    // Create a program with source code
    /* 6. 使用OpenCL C代码创建并编译出一个程序 */
    cl_program program = clCreateProgramWithSource(context, 1, 
        (const char**)&programSource, NULL, &status); // 使用源码创建程序
 
    // Build (compile) the program for the device
    status = clBuildProgram(program, numDevices, devices, 
        NULL, NULL, NULL); // 为设备构建(编译)程序
 
    // Create the vector addition kernel
    /* 7. 从编译好的OpenCL程序中提取内核 */
    cl_kernel kernel;
    kernel = clCreateKernel(program, "vecadd", &status); // 创建向量相加内核
 
    // Associate the input and output buffers with the kernel
    /* 8. 执行内核 */
    status = clSetKernelArg(kernel, 0, sizeof(cl_mem), &bufA); // 设置内核参数
    status = clSetKernelArg(kernel, 1, sizeof(cl_mem), &bufB);
    status = clSetKernelArg(kernel, 2, sizeof(cl_mem), &bufC);
 
    // Define an index space (global work size) of work 
    // items for execution. A workgroup size (local work size) 
    // is not required, but can be used.
    size_t globalWorkSize[1];   
 
    // There are 'elements' work-items 
    globalWorkSize[0] = elements; // 定义工作项的空间维度和空间大小
 
    // Execute the kernel for execution
    status = clEnqueueNDRangeKernel(cmdQueue, kernel, 1, NULL, 
        globalWorkSize, NULL, 0, NULL, NULL); // 通过执行API执行内核
 
    // Read the device output buffer to the host output array
    /* 9. 拷贝输出数据到主机端 */
    clEnqueueReadBuffer(cmdQueue, bufC, CL_TRUE, 0, 
        datasize, C, 0, NULL, NULL); // 将输出数组拷贝到主机端内存中
 
    // Verify the output
    int result = 1;
    for(i = 0; i < elements; i++) {
        if(C[i] != i+i) {
            result = 0;
            break;
        }
    }
    if(result) {
        std::cout << "result: " << *C << std::endl;
        printf("Output is correct\n");
    } else {
        printf("Output is incorrect\n");
    }
 
    // Free OpenCL resources
    /* 10. 释放资源 */
    clReleaseKernel(kernel);
    clReleaseProgram(program);
    clReleaseCommandQueue(cmdQueue);
    clReleaseMemObject(bufA);
    clReleaseMemObject(bufB);
    clReleaseMemObject(bufC);
    clReleaseContext(context);
 
    // Free host resources
    free(A);
    free(B);
    free(C);
    free(platforms);
    free(devices);
 
    return 0;
}

3.用OpenCL来操作映射内存中的数据

用GPU进行加速运行运算时，通常首先将数据copy（clEnqueueWriteBuffer）到GPU缓存对象，运算结束后，再将数据copy（clEnqueueReadBuffer）到内存；OpenCL提供了内存映射机制，无需读写操作，将设备上的内存映射到主机上，便可以通过指针方式直接修改主机上的内存对象，内存映射的运行性能远远高于普通的读写函数。

Step1:调用函数clEnqueueMapBuffer或函数clEnqueueMapImage,将内存映射命令入队;

Step2:使用诸如memcpy之类的函数，对内存中的数据进行传输操作。

Step3:调用clEnqueueUnmapMemObject函数解映射。

clCreateBuffer(cl_context,   //上下文
              cl_mem_flags,    //内存对象的性质，见下表
              size_t,          //内存对象数据块大小
              void *,          //host_ptr主机数据内存地址（可以为空）
              cl_int *)        // err错误码

cl_mem_flags	含义
CL_MEM_READ_WRITE（默认）	指明的是在kernel中对该缓冲区的访问权限：读写。限制设备的读写权限的，而不是主机端读写权限。
CL_MEM_WRITE_ONLY	指明的是在kernel中对该缓冲区的访问权限：只写。限制设备的读写权限的，而不是主机端读写权限。
CL_MEM_READ_ONLY	指明的是在kernel中对该缓冲区的访问权限：只读。限制设备的读写权限的，而不是主机端读写权限。
CL_MEM_COPY_HOST_PTR	1.host_ptr指针不能为空； 2.由于是直接拷贝，Gpu操作buffer时不会影响host_ptr 3.CL_MEM_USE_HOST_PTR与CL_MEM_COPY_HOST_PTR不能够同时使用
CL_MEM_USE_HOST_PTR	1.host_ptr指针不能为空； 2.在device中开辟空间用host指针内存区域内的数据来初始化,device执行的时候直接使用device上的内存。 3.在Mali GPU上,将数据拷贝到另一块区域,map的时候再拷贝回来,浪费时间,不建议使用该选项分配内存对象。 4.但是对于有单独显存的设备,此种分配内存对象的方式还有意义的。 5.Gpu操作buffer时会影响host_ptr 6.CL_MEM_USE_HOST_PTR与CL_MEM_ALLOC_HOST_PTR不能够同时使用 7.CL_MEM_USE_HOST_PTR与CL_MEM_COPY_HOST_PTR不能够同时使用
CL_MEM_ALLOC_HOST_PTR	1.零拷贝内存 2.从host指针可以访问的内存区域申请内存 3.在Mali系列的GPU上,这将有利于系统性能提升；不适合单独显存的设备。 4.独显存的设备CL_MEM_ALLOC_HOST_PTR和CL_MEM_COPY_HOST_PTR一起使用，比使用CL_MEM_USE_HOST_PTR要速度快。

后三个参数中，如果主机端的数据直接传输后不需要再读取结果，就使用CL_MEM_COPY_HOST_PTR，如果数据需要再读取回来，则可以使用CL_MEM_ALLOC_HOST_PTR和CL_MEM_COPY_HOST_PTR一起或者CL_MEM_USE_HOST_PTR

五、透视变换WarpPerspective算法原理

所谓的透视变换，就是三维世界中的物体经过一个变换矩阵转化为一张二维图像，这个过程就是几何投影的过程，这个矩阵就是一个透视变换矩阵，包括尺度、旋转、剪切、平移的变化。比如从不同角度去拍摄一个正方形，那么他所成的图像会得到不同的平行四边形，该成像的投影过程就是透视变换。

对于仿射变换来说，其投影向量(透视变换)等于0。可以认为是透视变换的特例。

透视变换主要用于处理三维场景中的透视投影问题，而仿射变换则更适用于处理二维场景中的平移、旋转和缩放等问题。

透视变换是将图片投影到一个新的视平面，也称作投影映射。通用的变换公式为：

[ x ′ , y ′ , w ′ ] = [ u , v , w ] [ h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ] \begin{bmatrix}x',y',w'\end{bmatrix}=\begin{bmatrix}u,v,w\end{bmatrix}\begin{bmatrix}h_{11}&h_{12}&h_{13}\\h_{21}&h_{22}&h_{23}\\h_{31}&h_{32}&h_{33}\end{bmatrix} [x′,y′,w′]=[u,v,w] h11h21h31h12h22h32h13h23h33

(u,v)是原始图片坐标，对应得到变换后的图片坐标(x,y),其中 x = x ′ / w ′ , y = y ′ / w ′ x=x'/w',y=y'/w' x=x′/w′,y=y′/w′。

变换矩阵 [ h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ] \begin{bmatrix}h_{11}&h_{12}&h_{13}\\h_{21}&h_{22}&h_{23}\\h_{31}&h_{32}&h_{33}\end{bmatrix} h11h21h31h12h22h32h13h23h33 可以拆成4部分， [ h 11 h 12 h 21 h 22 ] \begin{bmatrix}h_{11}&h_{12}\\h_{21}&h_{22}\end{bmatrix} [h11h21h12h22]表示线性变换，比如scaling，shearing和ratotion。 [ h 31 h 32 ] \begin{bmatrix}h_{31}&h_{32}\end{bmatrix} [h31h32]用于平移， [ h 13 h 23 ] \begin{bmatrix}h_{13}\\h_{23}\end{bmatrix} [h13h23]产生透视变换。所以可以理解成仿射等是透视变换的特殊形式。经过透视变换之后的图片通常不是平行四边形（除非映射视平面和原来平面平行的情况)。

重写之前的变换公式可以得到：

x = x ′ w ′ = h 11 u + h 21 v + h 31 h 13 u + h 23 v + h 33 y = y ′ w ′ = h 12 u + h 22 v + h 32 h 13 u + h 23 v + h 33 \begin{aligned}x&=\frac{x'}{w'}=\frac{h_{11}u+h_{21}v+h_{31}}{h_{13}u+h_{23}v+h_{33}}\\y&=\frac{y'}{w'}=\frac{h_{12}u+h_{22}v+h_{32}}{h_{13}u+h_{23}v+h_{33}}\end{aligned} xy=w′x′=h13u+h23v+h33h11u+h21v+h31=w′y′=h13u+h23v+h33h12u+h22v+h32

在OpenCV提供了关于透视变换的函数warpPerspective。

其实OpenCV中很多图像变换的映射关系都是反直觉的，直觉告诉我们，该函数的输入是原图的像素坐标，通过映射表或矩阵运算，输出的是目标图像的像素坐标。其实不然，OpenCV是先取一个目标图像的坐标，然后根据映射关系定位到原图中去，再从原图中得到该坐标的像素值。然而通过映射关系得到的坐标通常不是一个整数，即并不是原图的一个像素坐标，所以还需根据原图中该坐标周围的像素值用到插值算法计算出该坐标应有的像素值。也就是：

det ⁡ ( x , y ) = src ⁡ ( f x ( x , y ) , f y ( x , y ) ) \det(x,y)=\operatorname{src}(f_x(x,y),f_y(x,y)) det(x,y)=src(fx(x,y),fy(x,y))

对于warpPerspective这个函数，我们已知它表示的映射关系了：

d s t ( x , y ) = s r c ( h 11 x + h 12 y + h 13 h 31 x + h 32 y + h 33 , h 21 x + h 22 y + h 23 h 31 x + h 32 y + h 33 ) \mathrm{dst}(x,y)=\mathrm{src}\left(\frac{h_{11}x+h_{12}y+h_{13}}{h_{31}x+h_{32}y+h_{33}},\frac{h_{21}x+h_{22}y+h_{23}}{h_{31}x+h_{32}y+h_{33}}\right) dst(x,y)=src(h31x+h32y+h33h11x+h12y+h13,h31x+h32y+h33h21x+h22y+h23)

特别提醒，公式里的x，y是目标图像的坐标！那么我们设原图坐标是 x o ， y o x_o，y_o xo，yo，则有：

x o = h 11 x + h 12 y + h 13 h 31 x + h 32 y + h 33 x_o=\frac{h_{11}x+h_{12}y+h_{13}}{h_{31}x+h_{32}y+h_{33}} xo=h31x+h32y+h33h11x+h12y+h13

y o = h 21 x + h 22 y + h 23 h 31 x + h 32 y + h 33 y_o=\frac{h_{21}x+h_{22}y+h_{23}}{h_{31}x+h_{32}y+h_{33}} yo=h31x+h32y+h33h21x+h22y+h23

当我们想计算原图中某个特定的点在目标图像中的位置时，就需要解这个二元一次方程组了，解得：

x = ( h 22 − h 32 y o ) ( h 33 x o − h 13 ) − ( h 12 − h 32 x o ) ( h 33 y o − h 23 ) ( h 22 − h 32 y o ) ( h 11 − h 31 x o ) − ( h 12 − h 32 x o ) ( h 21 − h 31 y o ) x=\frac{(h_{22}-h_{32}y_o)(h_{33}x_o-h_{13})-(h_{12}-h_{32}x_o)(h_{33}y_o-h_{23})}{(h_{22}-h_{32}y_o)(h_{11}-h_{31}x_o)-(h_{12}-h_{32}x_o)(h_{21}-h_{31}y_o)} x=(h22−h32yo)(h11−h31xo)−(h12−h32xo)(h21−h31yo)(h22−h32yo)(h33xo−h13)−(h12−h32xo)(h33yo−h23)

y = ( h 21 − h 31 y o ) ( h 33 x o − h 13 ) − ( h 11 − h 31 x o ) ( h 33 y o − h 23 ) ( h 21 − h 31 y o ) ( h 12 − h 32 x o ) − ( h 11 − h 31 x o ) ( h 22 − h 32 y o ) y=\frac{(h_{21}-h_{31}y_o)(h_{33}x_o-h_{13})-(h_{11}-h_{31}x_o)(h_{33}y_o-h_{23})}{(h_{21}-h_{31}y_o)(h_{12}-h_{32}x_o)-(h_{11}-h_{31}x_o)(h_{22}-h_{32}y_o)} y=(h21−h31yo)(h12−h32xo)−(h11−h31xo)(h22−h32yo)(h21−h31yo)(h33xo−h13)−(h11−h31xo)(h33yo−h23)

其中所用到的插值算法为双线性插值，其算法原理如下：

一般的图像插值算法有最近邻插值、线性插值、双线性插值、高阶插值等算法，本文针对的对象是图像，所用到的是双线性插值法。

在介绍双线性插值之前，我们先介绍一下什么是单线性插值。

最近邻插值只需要在原图像矩阵找到离它最近的一个像素点，然后进行赋值即可。而单线性插值需要找到周围最近的两个点进行插值运算。

Warp算法OpenCL异构开发总结,请添加图片描述,第4张

如上图所示：已知点 ( x 0 , y 0 ) (x_0,y_0) (x0,y0)， ( x , y ) (x,y) (x,y) ， ( x 1 , y 1 ) (x_1,y_1) (x1,y1)要得到区间 [ x 1 , x 2 ] [x_1,x_2] [x1,x2]中某一位置 x x x在直线上的 y y y值. 根据上图，我们可以得到：

y − y 0 x − x 0 = y 1 − y 0 x 1 − x 0 \frac{y-y_0}{x-x_0}=\frac{y_1-y_0}{x_1-x_0} x−x0y−y0=x1−x0y1−y0

进而得到：

y = x 1 − x x 1 − x 0 y 0 + x − x 0 x 1 − x 0 y 1 y=\frac{x_1-x}{x_1-x_0}y_0+\frac{x-x_0}{x_1-x_0}y_1 y=x1−x0x1−xy0+x1−x0x−x0y1

从上式中我们可以知道，分子可以看成 x x x分别到 x 1 x_1 x1和 x 0 x_0 x0的距离. 把整个分式分别看成 y 0 y_0 y0与 y 1 y_1 y1的加权系数. (离哪个点近，权重就高)

我们要根据原图像中2个点的像素值求未知点的像素值. 因此我们可以这样做，将上式中的 y 0 y_0 y0与 y 1 y_1 y1分别代表原图像中的像素值，上面的公式可写成如下形式：

f ( p ) = x 1 − x x 1 − x 0 f ( p 0 ) + x − x 0 x 1 − x 0 f ( p 1 ) f(p)=\frac{x_1-x}{x_1-x_0}f(p_0)+\frac{x-x_0}{x_1-x_0}f(p_1) f(p)=x1−x0x1−xf(p0)+x1−x0x−x0f(p1)

其本质就是在 x x x方向进行了一次线性插值，从而得到目标图像的像素值。

下面我们开始介绍双线性插值算法，双线性插值中的"双"，指的是对x方向与y方向分别进行插值，其实共进行了三次单线性插值，所以这个"双"应该理解成两个方向的插值，而非插值的次数.

双线性插值的原理就是基于单线性插值拓展的，其原理如下图所示：

Warp算法OpenCL异构开发总结,请添加图片描述,第5张

通过周围四个像素点进行插值得到目标点，假设已知四个像素值 f ( P 11 ) , f ( P 12 ) , f ( P 21 ) , f ( P 22 ) f(P_{11}),f(P_{12}),f(P_{21}),f(P_{22}) f(P11),f(P12),f(P21),f(P22)其求解步骤如下：

第一步：利用单线性插值在 x x x方向求 f ( Q 1 ) f(Q_1) f(Q1)：

f ( Q 1 ) = x 1 − x x 1 − x 0 f ( P 11 ) + x − x 0 x 1 − x 0 f ( P 21 ) f(Q_1)=\frac{x_1-x}{x_1-x_0}f(P_{11})+\frac{x-x_0}{x_1-x_0}f(P_{21}) f(Q1)=x1−x0x1−xf(P11)+x1−x0x−x0f(P21)

第二步：利用单线性插值在 x x x方向求 f ( Q 2 ) f(Q_2) f(Q2)：

f ( Q 2 ) = x 1 − x x 1 − x 0 f ( P 12 ) + x − x 0 x 1 − x 0 f ( P 22 ) f(Q_2)=\frac{x_1-x}{x_1-x_0}f(P_{12})+\frac{x-x_0}{x_1-x_0}f(P_{22}) f(Q2)=x1−x0x1−xf(P12)+x1−x0x−x0f(P22)

第三步：利用单线性插值在 y y y方向求 f ( R ) f(R) f(R)：

f ( R ) = y 1 − y y 1 − y 0 f ( Q 1 ) + y − y 0 y 1 − y 0 f ( Q 2 ) f(R)=\frac{y_1-y}{y_1-y_0}f(Q_1)+\frac{y-y_0}{y_1-y_0}f(Q_2) f(R)=y1−y0y1−yf(Q1)+y1−y0y−y0f(Q2)

联合以上三式，可求得：

f ( R ) = y 1 − y y 1 − y 0 [ x 1 − x x 1 − x 0 f ( P 11 ) + x − x 0 x 1 − x 0 f ( P 21 ) ] + y − y 0 y 1 − y 0 [ x 1 − x x 1 − x 0 f ( P 12 ) + x − x 0 x 1 − x 0 f ( P 22 ) ] f(R)=\frac{y_1-y}{y_1-y_0}[\frac{x_1-x}{x_1-x_0}f(P_{11})+\frac{x-x_0}{x_1-x_0}f(P_{21})]+\frac{y-y_0}{y_1-y_0}[\frac{x_1-x}{x_1-x_0}f(P_{12})+\frac{x-x_0}{x_1-x_0}f(P_{22})] f(R)=y1−y0y1−y[x1−x0x1−xf(P11)+x1−x0x−x0f(P21)]+y1−y0y−y0[x1−x0x1−xf(P12)+x1−x0x−x0f(P22)]

在图像中，其实不难发现，在计算中以下关系：

{ x 1 − x 0 = 1 y 1 − y 0 = 1 \left.\left\{\begin{array}{l}x_1-x_0=1\\y_1-y_0=1\end{array}\right.\right. {x1−x0=1y1−y0=1

因此，可以简化为：

f ( R ) = ( y 1 − y ) ( x 1 − x ) f ( P 11 ) + ( y 1 − y ) ( x − x 0 ) f ( P 21 ) + ( y − y 0 ) ( x 1 − x ) f ( P 12 ) + ( y − y 0 ) ( x − x 0 ) f ( P 22 ) f(R)=(y_1-y)(x_1-x)f(P_{11})+(y_1-y)(x-x_0)f(P_{21})+(y-y_0)(x_1-x)f(P_{12})+(y-y_0)(x-x_0)f(P_{22}) f(R)=(y1−y)(x1−x)f(P11)+(y1−y)(x−x0)f(P21)+(y−y0)(x1−x)f(P12)+(y−y0)(x−x0)f(P22)

也可以写成权重的形式，如下：

f ( R ) = f ( P 11 ) w 11 + f ( P 12 ) w 12 + f ( P 21 ) w 21 + f ( P 22 ) w 22 f(R)=f(P_{11})w_{11}+f(P_{12})w_{12}+f(P_{21})w_{21}+f(P_{22})w_{22} f(R)=f(P11)w11+f(P12)w12+f(P21)w21+f(P22)w22

观察一下，就会发现每个点的权重和待求点与对角线的距离有关.

六、对RGB图像进行warp操作的OpenCL源码

由于RGB图像是三通道，而OpenCL是对所有像素点并行处理，因此内核参数width=3×图像宽度，global_work_size的第一维也需设为3×图像宽度。具体核函数代码如下：

static const char* kernel_code =
"__kernel void kernel_warp_rgb(__global unsigned char* inputImage, __global unsigned char* outputImage, int width, int height, __global float* H)\n"
"{\n"
    "int x = get_global_id(0);\n"
    "int y = get_global_id(1);\n"
    "int offset = y * width + x;"
    "\n"
    "float w=(H[6] * x + H[7] * y + H[8]);\n"
    "float dstX = (H[0] * x + H[1] * y + H[2]) / w;\n"
    "float dstY = (H[3] * x + H[4] * y + H[5]) / w;\n"
    "\n"
    "int x0 = (int)dstX;\n"
    "int y0 = (int)dstY;\n"
    "\n"
    "float x_diff = dstX - x0;\n"
    "float y_diff = dstY - y0;\n"
    "\n"
    "float fac1=(1-x_diff)*(1-y_diff);\n"
    "float fac2=x_diff*(1-y_diff);\n"
    "float fac3=(1-x_diff)*y_diff;\n"
    "float fac4=x_diff*y_diff;\n"
    "\n"
    "outputImage[y * width + x]=fac1* inputImage[y0 * width + x0] + fac2 * inputImage[y0 * width + x0+1] + fac3 * inputImage[y0+1 * width + x0] + fac4* inputImage[y0+1 * width + x0+1];\n"
    "\n"
"}\n"
"\n"

若内核参数width和global_work_size的第一维不想改动，仍然设为图像宽度，则需要在核函数中利用像素偏移分别对图像的RGB三个通道进行处理。

static const char* kernel_code =
"__kernel void kernel_warp_rgb(__global unsigned char* inputImage, __global unsigned char* outputImage, int width, int height, __global float* H)\n"
"{\n"
    "int x = get_global_id(0);\n"
    "int y = get_global_id(1);\n"
    "int offset = y * width + x;"
    "\n"
    "float w=(H[6] * x + H[7] * y + H[8]);\n"
    "float dstX = (H[0] * x + H[1] * y + H[2]) / w;\n"
    "float dstY = (H[3] * x + H[4] * y + H[5]) / w;\n"
    "\n"
    "int x0 = (int)dstX;\n"
    "int y0 = (int)dstY;\n"
    "\n"
    "float x_diff = dstX - x0;\n"
    "float y_diff = dstY - y0;\n"
    "\n"
    "float fac1=(1-x_diff)*(1-y_diff);\n"
    "float fac2=x_diff*(1-y_diff);\n"
    "float fac3=(1-x_diff)*y_diff;\n"
    "float fac4=x_diff*y_diff;\n"
    "\n"
    "int offset0 = y0 * width + x0;"
    "int offset1 = y0 * width + x0+1;"
    "int offset2 = y0+1 * width + x0;"
    "int offset3 = y0+1 * width + x0+1;"
    // R通道
    "outputImage[3 * offset + 0]=fac1* inputImage[3 * offset0 + 0] + fac2 * inputImage[3 * offset1 + 0] + fac3 * inputImage[3 * offset2 + 0] + fac4* inputImage[3 * offset3 + 0];\n"
    // G通道
    "outputImage[3 * offset + 1]=fac1* inputImage[3 * offset0 + 1] + fac2 * inputImage[3 * offset1 + 1] + fac3 * inputImage[3 * offset2 + 1] + fac4* inputImage[3 * offset3 + 1];\n"
    // B通道
    "outputImage[3 * offset + 2]=fac1* inputImage[3 * offset0 + 2] + fac2 * inputImage[3 * offset1 + 2] + fac3 * inputImage[3 * offset2 + 2] + fac4* inputImage[3 * offset3 + 2];\n"
    "\n"
"}\n"
"\n"

参考资料：

1.《OpenCL 2.0 异构计算 [第三版](Heterogeneous Computing with OpenCL 2.0)》

2.【图像处理】透视变换 Perspective Transformationhttps://blog.csdn.net/xiaowei_cqu/article/details/26471527

3.插值算法https://zhuanlan.zhihu.com/p/266951945

免费加速器看国外网站分销商城系统织梦官网如何成立一个公司流程上海网站建设系统湛江网站

上一篇：分治算法总结（Java）

下一篇：VSCode中写Vue没有代码提示的解决办法