【IT168 文档】SDK2.3的convolutionSeparable示例，纯代码，零注释，忒血汗。。汗了半小时才o掉，帖出来供大家参考。
离散数据的二维卷积：
其中，Ar、Ac分别是A的行数与列数。应用很多，比如对图像做高斯平滑(去噪)，拿高斯核与输入图像做卷积。
convolutionSeparable之所以”Separable”，是因为它在row、col两个维上分别做了卷积操作。在此先奉上CPU代码，无敌明了<本帖只讲述row方向上的，col上的太类似了，自己看咯>
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代码，顺风顺水了是不是。。
#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");
}
