孤独明镜

菩提本无树,明镜亦非台,本来无一物,何处惹尘埃。

嗨,我是Dragon,一名全栈开发者。现居上海,就职于一家游戏开发公司。做过Linux/QT、WEB(PHP/Python)、游戏开发,目前主攻图形/图像方向,对数学、文学、哲学、摄影等非常的喜爱。


积分图快速运算及CUDA的实现

我是使用按行求前缀和后再按列前缀和来生成积分图的。
1)先按行求每行的前缀和,例如求第x行的前缀和。
x行中y位置的值为x行中<=y的点的值累加,按行生成一个新的图。
2)对上面新生成的图按列做上面的操作。x行y列的点的值是<=x并在y列的新图中所有点的累加。此步生成的图即是积分图。
更合适的算法可以在这里找到《积分图像的快速GPU 计算》 ,这里我使用CUDA实现我上面所说的算法。

__global__ void
integrogram(int *A, int *B, int w)
{
    extern __shared__ int sdata[];
    int i = threadIdx.x;
    int j = threadIdx.y;

    int k = 0;

    //此处是按行求前缀和
    for(k = 0; k <= j; k++)
    {
        B[i * w + j] += A[i * w + k];
    }

    sdata[i * w + j] = B[i * w + j];

    __syncthreads();  //此处等待末完成的线程,所有线程执行到这里后才向下执行

    //此处是按列求前缀和
    for(k = 0; k < i; k++)
    {
        B[i * w + j] += sdata[k * w + j];
    }
}

主机端程序如下:

int N = 5;
int size = N * N;

int h_A[size];
int h_B[size];

int i, j;

memset(&h_A, 0, sizeof(int) * size);
memset(&h_B, 0, sizeof(int) * size);

printf("print host data A before deal:\n");
for (i = 0; i < N; ++i)
{
    for (j = 0; j < N; ++j)
    {
        h_A[i * N + j] = 1+(int)(10.0*rand()/(RAND_MAX+1.0));
        printf("%i ", h_A[i * N + j]);
    }
    printf("\n");
}

int *d_A = NULL;
err = cudaMalloc((void **)&d_A, sizeof(int) * size);
if (err != cudaSuccess)
{
    fprintf(stderr, "Failed to allocate device A (error code %s)!\n", cudaGetErrorString(err));
    exit(EXIT_FAILURE);
}

int *d_B = NULL;
err = cudaMalloc((void **)&d_B, sizeof(int) * size);
if (err != cudaSuccess)
{
    fprintf(stderr, "Failed to allocate device B (error code %s)!\n", cudaGetErrorString(err));
    exit(EXIT_FAILURE);
}

err = cudaMemcpy(d_A, &h_A, sizeof(int) * size, cudaMemcpyHostToDevice);
if (err != cudaSuccess)
{
    fprintf(stderr, "Failed to copy A from host to device (error code %s)!\n", cudaGetErrorString(err));
    exit(EXIT_FAILURE);
}

err = cudaMemcpy(d_B, &h_B, sizeof(int) * size, cudaMemcpyHostToDevice);
if (err != cudaSuccess)
{
    fprintf(stderr, "Failed to copy B from host to device (error code %s)!\n", cudaGetErrorString(err));
    exit(EXIT_FAILURE);
}

dim3 dimBlock(N, N);
integrogram<<<1, dimBlock, sizeof(int) * size>>>(d_A, d_B, N);

执行结果如下图:

转载请注明地址:孤独明镜