opencv-cuda开发（5）当工作量远大于线程量时的跨步处理，拓展到图像处理领域

参考下一个动图

接下来这部分是引述的博文： 

 有时，工作量比网格大，或者出于某种原因，一个网格中的线程数量可能会小于实际工作量的大小。请思考一下包含 1000 个元素的数组和包含 250 个线程的网格（此处使用极小的规模以便于说明）。此网格中的每个线程将需使用 4 次。如要实现此操作，一种常用方法便是在核函数中使用跨网格循环。

在跨网格循环中，每个线程将在网格内使用 threadIdx + blockIdx*blockDim 计算自身唯一的索引，并对数组内该索引的元素执行相应运算，然后用网格中的线程数加上自身索引值，并重复此操作，直至超出数组范围。

例如，对于包含 500 个元素的数组 a 和包含 250 个线程的网格，网格中索引为 20 的线程将执行如下操作：

• 对 a[20] 执行相应运算；
• 将线程索引增加 250，使网格的大小达到 270
• 对a[270] 执行相应运算；
• 将线程索引增加 250，使网格的大小达到 520
• 由于 520 现已超出数组范围，因此线程将停止工作。

CUDA 提供一个记录了网格中线程块数的变量：gridDim.x。然后可以利用它来计算网格中的总线程数，即网格中的线程块数乘以每个线程块中的线程数：gridDim.x * blockDim.x。现在来看看以下核函数中网格跨度循环的示例：

#include <stdio.h>

// 初始化数组a
void init(int *a, int N) {
  int i;
  for (i = 0; i < N; ++i) {
    a[i] = i;
  }
}

__global__ void doubleElements(int *a, int N) {

  // 使用grid-stride循环，这样每个线程可以处理数组中的多个元素。
  int idx = blockIdx.x * blockDim.x + threadIdx.x;
  int stride = gridDim.x * blockDim.x; // grid 的一个跨步

  for (int i = idx; i < N; i += stride) {
    a[i] *= 2;
  }
}

// 检查数组内所有元素的值是否均为复数
bool checkElementsAreDoubled(int *a, int N) {
  int i;
  for (i = 0; i < N; ++i) {
    if (a[i] != i*2) return false;
  }
  return true;
}

int main() {
  int N = 10000;
  int *a;
  size_t size = N * sizeof(int);
  cudaMallocManaged(&a, size);

  init(a, N); // 初始化数组a

  size_t threads_per_block = 256; // 每个block的thread数量
  size_t number_of_blocks = 32; // block数量

  doubleElements<<<number_of_blocks, threads_per_block>>>(a, N);
  cudaDeviceSynchronize();

  bool areDoubled = checkElementsAreDoubled(a, N); 
// 检查数组内所有元素的值是否均为复数
  printf("All elements were doubled? %s\n", areDoubled ? "TRUE" : "FALSE");

  cudaFree(a);
}

引述结束，以上是对一维数组的跨步处理，接下来将思路拓展到二维

例子还是将两张100*100的图相减

线程数量足够的时候是这样写的核函数，一个线程对应一个图片像素，需要100*100个线程

    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;

    image1_dev[y * width + x] = abs(iamge1_dev[y * width + x] - image2_dev[y * width + x]); 

在线程足够的情况下，线程块0的线程0仅仅需要处理（0，0）的像素点

当把线程数量限制到100个的时候，加上跨度处理的思路

那线程块0的线程0需要处理的像素点变成了（0，0）（0，1）（0，2）.......（0，99）这100个像素点

那么跨度就可以看出来了，是在y轴上的有1的跨度，那么去内存里面怎么去遍历这100个像素点呢，这个时候就需要用到 idx = ny * iy + ix这个式子，即线程块0的线程0需要处理的像素点有0 99 99*1 一直到9999。这个思路太死板了。

从遍历内存的角度上来理解更为容易，这里有10000个像素点需要遍历，不关心它曾经是二维还是三维的，现在都是一维，那么通过任何线程数量完成对它的遍历即算结束，100个线程需要100次，15个线程就需要667次，99个线程就需要102次

所以说stride是根据当前设计的线程总量来进行计算的，根据行来展开，需要关心的只有X轴向 即stride = gridDim.x * blockDim.x; 

比如32*64的情况下，stride=32*64=2048

那代码就可以这样写

    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    int idx = y * width + x;
    int stride = gridDim.x * blockDim.x; // grid 的一个跨步

    
    for (int i = idx; i < N; i += stride)   //这里的N是总像素点个数width*height
    {
        image1_dev[i] = abs(iamge1_dev[i] - image2_dev[i]); 
    }

在每次线程里都去跳着完成多次求差运算，进而可以循环覆盖整个数据存储区域，最终遍历完所有图片数据

这里最好加一个内存指针总长的判定，在处理完成图片数据后停止计算，不这样判定会存在一个浪费计算量的问题，比如1000个图像数据，给分配了32个线程，那么遍历到最后一组的线程的时候，计算只有前面8次是有效的，空跑了24次。

提示：最好不要在核函数里写循环操作，这个效率会低不少，除非线程数量不够用的情况下

参考博客

CUDA C/C++ 教程一：加速应用程序_c++ cuda_白水baishui的博客-CSDN博客