扩展CUDA SDK 2.3 の convolutionSeparable

SDK2.3的convolutionSeparable示例，纯代码，零注释，忒血汗。。汗了半小时才o掉，帖出来供大家参考。

离散数据的二维卷积：

 

其中，Ar、Ac分别是A的行数与列数。应用很多，比如对图像做高斯平滑(去噪)，拿高斯核与输入图像做卷积。

convolutionSeparable之所以”Separable”，是因为它在row、col两个维上分别做了卷积操作。在此先奉上CPU代码，无敌明了<本帖只讲述row方向上的，col上的太类似了，自己看咯>

// Reference row convolution filter

extern "C" void convolutionRowCPU(

    float *h_Dst,

    float *h_Src,

    float *h_Kernel,

    int imageW,

    int imageH,

    int kernelR//-8

){

    for(int y = 0; y < imageH; y++)

        for(int x = 0; x < imageW; x++){

            float sum = 0;

            for(int k = -kernelR; k <= kernelR; k++){

                int d = x + k;//即左右各8个，外加自己本位上的，共17个元素做邻域加权

                if(d >= 0 && d < imageW)

                    sum += h_Src[y * imageW + d] * h_Kernel[kernelR - k];//h_src视为上式中B阵，h_kernel视为A阵（本例中h_kernel为随机给出的核，用户可以自己写高斯核玩）

            }

            h_Dst[y * imageW + x] = sum;

        }

}

好了，下面想想CUDA怎么实现，每个block完成什么样的任务，每个thread又负责完成怎样的任务。下图分别是每个block的共享存储体数组s_Data[4][96]、全局存储器里的输入数组d_Src[3072/2][3072]。线程组织结构是这样的：grid( imagW/(ROWS_RESULT_STEPS*ROWS_BLOCKDIM_X), imagH/ROWS_BLOCKDIM_Y ), block(ROWS_BLOCKDIM_X, ROWS_BLOCKDIM_Y)。其中ROWS_BLOCKDIM_X为16，ROWS_BLOCKDIM_Y为4，表明block内线程组织结构是16*4；ROWS_RESULT_STEPS 为4，表明一个block每次做4轮操作，看图中s_Data的绿色部分，每次一个block即16*4个线程对应计算一个蓝色框框标注的部分，一个线程负责计算一个位置上的数据，做4轮于是就有了4列绿矩嘛。注意，最左边的一列红阵和最右边的一列橙阵特别标出，这是在进行一些“越界处理”，例如计算d_src[][]最左边的元素时，它们的“左边8个元素”(事实上已经不能再左了)，这时红色阵置0，同理理解橙阵。

 

 

 

下面再来看这段kernel代码，顺风顺水了是不是。。

// Row convolution filter

#define   ROWS_BLOCKDIM_X 16

#define   ROWS_BLOCKDIM_Y 4

#define ROWS_RESULT_STEPS 4

#define   ROWS_HALO_STEPS 1

 

__global__ void convolutionRowsKernel(

    float *d_Dst,

    float *d_Src,

    int imageW,

    int imageH,

    int pitch

){

    __shared__ float s_Data[ROWS_BLOCKDIM_Y][(ROWS_RESULT_STEPS + 2 * ROWS_HALO_STEPS) * ROWS_BLOCKDIM_X];

 

    //Offset to the left halo edge

    const int baseX = (blockIdx.x * ROWS_RESULT_STEPS - ROWS_HALO_STEPS) * ROWS_BLOCKDIM_X + threadIdx.x;

    const int baseY = blockIdx.y * ROWS_BLOCKDIM_Y + threadIdx.y;

 

    d_Src += baseY * pitch + baseX;

    d_Dst += baseY * pitch + baseX;

 

    //Main data

    #pragma unroll

    for(int i = ROWS_HALO_STEPS; i < ROWS_HALO_STEPS + ROWS_RESULT_STEPS; i++)//i=1,i<5

        s_Data[threadIdx.y][threadIdx.x + i * ROWS_BLOCKDIM_X] = d_Src[i * ROWS_BLOCKDIM_X];

 

    //Left halo

    for(int i = 0; i < ROWS_HALO_STEPS; i++){//i=0,i<1

        s_Data[threadIdx.y][threadIdx.x + i * ROWS_BLOCKDIM_X] =

            (baseX >= -i * ROWS_BLOCKDIM_X ) ? d_Src[i * ROWS_BLOCKDIM_X] : 0;

    }

 

    //Right halo

    for(int i = ROWS_HALO_STEPS + ROWS_RESULT_STEPS; i < ROWS_HALO_STEPS + ROWS_RESULT_STEPS + ROWS_HALO_STEPS; i++){

        s_Data[threadIdx.y][threadIdx.x + i * ROWS_BLOCKDIM_X] =

            (imageW - baseX > i * ROWS_BLOCKDIM_X) ? d_Src[i * ROWS_BLOCKDIM_X] : 0;

    }

 

    //Compute and store results

    __syncthreads();

    #pragma unroll

    for(int i = ROWS_HALO_STEPS; i < ROWS_HALO_STEPS + ROWS_RESULT_STEPS; i++){

        float sum = 0;

 

        #pragma unroll

        for(int j = -KERNEL_RADIUS; j <= KERNEL_RADIUS; j++)

            sum += c_Kernel[KERNEL_RADIUS - j] * s_Data[threadIdx.y][threadIdx.x + i * ROWS_BLOCKDIM_X + j];

 

        d_Dst[i * ROWS_BLOCKDIM_X] = sum;

    }

}

 

extern "C" void convolutionRowsGPU(

    float *d_Dst,

    float *d_Src,

    int imageW,

    int imageH

){

    assert( ROWS_BLOCKDIM_X * ROWS_HALO_STEPS >= KERNEL_RADIUS );

    assert( imageW % (ROWS_RESULT_STEPS * ROWS_BLOCKDIM_X) == 0 );

    assert( imageH % ROWS_BLOCKDIM_Y == 0 );

 

    dim3 blocks(imageW / (ROWS_RESULT_STEPS * ROWS_BLOCKDIM_X), imageH / ROWS_BLOCKDIM_Y);

    dim3 threads(ROWS_BLOCKDIM_X, ROWS_BLOCKDIM_Y);

 

    convolutionRowsKernel<<<blocks, threads>>>(

        d_Dst, d_Src, imageW, imageH, imageW

    );

    cutilCheckMsg("convolutionRowsKernel() execution failed/n");

}

 

 

另外，偶稀饭这张照片，也要帖出来。。不要拦偶。。

 

 

另外，我对程序进行了一些改进，使其可以处理任意宽*任意高的大小，不一定非要是16或64或什么什么的整数倍，而且对边界做了clamp处理，SDK的程序边界处理不好（越界的填充0，模糊两下边界就白白了），我是将越界的都填充为它最临近的那个像素点的值，这样模糊下来，边界效果几乎不变。具体程序可以参见我上传的这份源码，需要的朋友直接调用就好。