Box filter with cuda c
i have little experiance with cuda and trying to write a box filter. I read that a box filter is a filter in which each pixel in the resulting image has a value equal to the average value of its neighboring pixels in the input image. I have foud this document http://www.nvidia.com/content/nvision2008/tech_presentations/Game_Developer_Track/NVISION08-Image_Processing_and_Video_with_CUDA.pdf and changed the code a little bit. Here is my function.
#define TILE_W 16
#define TILE_H 16
#define R 2 // filter radius
#define D (R*2+1) // filter diameter
#define S (D*D) // filter size
#define BLOCK_W (TILE_W+(2*R))
#define BLOCK_H (TILE_H+(2*R))
__global__ void d_filter(unsigned char *g_idata, unsigned char *g_odata, unsigned int width, unsigned int height)
{
__shared__ unsigned char smem[BLOCK_W*BLOCK_H];
int x = blockIdx.x*TILE_W + threadIdx.x - R;
int y = blockIdx.y*TILE_H + threadIdx.y - R;
// clamp to edge of image
x = max(0, x);
x = min(x, width-1);
y = max(y, 0);
y = min(y, height-1);
unsigned int index = y*width + x;
unsigned int bindex = threadIdx.y*blockDim.y+threadIdx.x;
// each thread copies its pixel of the block to shared memory
smem[bindex] = g_idata[index];
__syncthreads();
// only threads inside the apron will write results
if ((threadIdx.x >= R) && (threadIdx.x < (BLOCK_W-R)) && (threadIdx.y >= R) && (threadIdx.y < (BLOCK_H-R))) {
float sum = 0;
for(int dy=-R; dy<=R; dy++) {
for(int dx=-R; dx<=R; dx++) {
float i = smem[bindex + (dy*blockDim.x) + dx];
sum += i;
}
}
g_odata[index] = sum / S;
}
}
EDIT: Here is a newer version that works. The Problem was in the kernel launch.
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#include <assert.h>
#define PGMHeaderSize 0x40
inline bool loadPPM(const char *file, unsigned char **data, unsigned int *w, unsigned int *h, unsigned int *channels)
{
FILE *fp = NULL;
fp = fopen(file, "rb");
if (!fp) {
fprintf(stderr, "__LoadPPM() : unable to open filen" );
return false;
}
// check header
char header[PGMHeaderSize];
if (fgets(header, PGMHeaderSize, fp) == NULL)
{
fprintf(stderr,"__LoadPPM() : reading PGM header returned NULLn" );
return false;
}
if (strncmp(header, "P5", 2) == 0)
{
*channels = 1;
}
else if (strncmp(header, "P6", 2) == 0)
{
*channels = 3;
}
else
{
fprintf(stderr,"__LoadPPM() : File is not a PPM or PGM imagen" );
*channels = 0;
return false;
}
// parse header, read maxval, width and height
unsigned int width = 0;
unsigned int height = 0;
unsigned int maxval = 0;
unsigned int i = 0;
while (i < 3)
{
if (fgets(header, PGMHeaderSize, fp) == NULL)
{
fprintf(stderr,"__LoadPPM() : reading PGM header returned NULLn" );
return false;
}
if (header[0] == '#')
{
continue;
}
if (i == 0)
{
i += sscanf(header, "%u %u %u", &width, &height, &maxval);
}
else if (i == 1)
{
i += sscanf(header, "%u %u", &height, &maxval);
}
else if (i == 2)
{
i += sscanf(header, "%u", &maxval);
}
}
// check if given handle for the data is initialized
if (NULL != *data)
{
if (*w != width || *h != height)
{
fprintf(stderr, "__LoadPPM() : Invalid image dimensions.n" );
}
}
else
{
*data = (unsigned char *) malloc(sizeof(unsigned char) * width * height * *channels);
if (!data) {
fprintf(stderr, "Unable to allocate hostmemoryn");
return false;
}
*w = width;
*h = height;
}
// read and close file
if (fread(*data, sizeof(unsigned char), width * height * *channels, fp) == 0)
{
fprintf(stderr, "__LoadPPM() : read data returned error.n" );
fclose(fp);
return false;
}
fclose(fp);
return true;
}
inline bool savePPM(const char *file, unsigned char *data, unsigned int w, unsigned int h, unsigned int channels)
{
assert(NULL != data);
assert(w > 0);
assert(h > 0);
std::fstream fh(file, std::fstream::out | std::fstream::binary);
if (fh.bad())
{
fprintf(stderr, "__savePPM() : Opening file failed.n" );
return false;
}
if (channels == 1)
{
fh << "P5n";
}
else if (channels == 3)
{
fh << "P6n";
}
else
{
fprintf(stderr, "__savePPM() : Invalid number of channels.n" );
return false;
}
fh << w << "n" << h << "n" << 0xff << std::endl;
for (unsigned int i = 0; (i < (w*h*channels)) && fh.good(); ++i)
{
fh << data[i];
}
fh.flush();
if (fh.bad())
{
fprintf(stderr,"__savePPM() : Writing data failed.n" );
return false;
}
fh.close();
return true;
}
#define TILE_W 16
#define TILE_H 16
#define Rx 2 // filter radius in x direction
#define Ry 2 // filter radius in y direction
#define FILTER_W (Rx*2+1) // filter diameter in x direction
#define FILTER_H (Ry*2+1) // filter diameter in y direction
#define S (FILTER_W*FILTER_H) // filter size
#define BLOCK_W (TILE_W+(2*Rx))
#define BLOCK_H (TILE_H+(2*Ry))
__global__ void box_filter(const unsigned char *in, unsigned char *out, const unsigned int w, const unsigned int h){
//Indexes
const int x = blockIdx.x * TILE_W + threadIdx.x - Rx; // x image index
const int y = blockIdx.y * TILE_H + threadIdx.y - Ry; // y image index
const int d = y * w + x; // data index
//shared mem
__shared__ float shMem[BLOCK_W][BLOCK_H];
if(x<0 || y<0 || x>=w || y>=h) { // Threads which are not in the picture just write 0 to the shared mem
shMem[threadIdx.x][threadIdx.y] = 0;
return;
}
shMem[threadIdx.x][threadIdx.y] = in[d];
__syncthreads();
// box filter (only for threads inside the tile)
if ((threadIdx.x >= Rx) && (threadIdx.x < (BLOCK_W-Rx)) && (threadIdx.y >= Ry) && (threadIdx.y < (BLOCK_H-Ry))) {
float sum = 0;
for(int dx=-Rx; dx<=Rx; dx++) {
for(int dy=-Ry; dy<=Ry; dy++) {
sum += shMem[threadIdx.x+dx][threadIdx.y+dy];
}
}
out[d] = sum / S;
}
}
#define checkCudaErrors(err) __checkCudaErrors (err, __FILE__, __LINE__)
inline void __checkCudaErrors(cudaError err, const char *file, const int line)
{
if (cudaSuccess != err)
{
fprintf(stderr, "%s(%i) : CUDA Runtime API error %d: %s.n",
file, line, (int)err, cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
}
int main(){
unsigned char *data=NULL, *d_idata=NULL, *d_odata=NULL;
unsigned int w,h,channels;
if(! loadPPM("../../data/lena_bw.pgm", &data, &w, &h, &channels)){
fprintf(stderr, "Failed to open Filen");
exit(EXIT_FAILURE);
}
printf("Loaded file with w:%d h:%d channels:%d n",w,h,channels);
unsigned int numElements = w*h*channels;
size_t datasize = numElements * sizeof(unsigned char);
// Allocate the Device Memory
printf("Allocate Devicememory for datan");
checkCudaErrors(cudaMalloc((void **)&d_idata, datasize));
checkCudaErrors(cudaMalloc((void **)&d_odata, datasize));
// Copy to device
printf("Copy idata from the host memory to the CUDA devicen");
checkCudaErrors(cudaMemcpy(d_idata, data, datasize, cudaMemcpyHostToDevice));
// Launch Kernel
int GRID_W = w/TILE_W +1;
int GRID_H = h/TILE_H +1;
dim3 threadsPerBlock(BLOCK_W, BLOCK_H);
dim3 blocksPerGrid(GRID_W,GRID_H);
printf("CUDA kernel launch with [%d %d] blocks of [%d %d] threadsn", blocksPerGrid.x, blocksPerGrid.y, threadsPerBlock.x, threadsPerBlock.y);
box_filter<<<blocksPerGrid, threadsPerBlock>>>(d_idata, d_odata, w,h);
checkCudaErrors(cudaGetLastError());
// Copy data from device to host
printf("Copy odata from the CUDA device to the host memoryn");
checkCudaErrors(cudaMemcpy(data, d_odata, datasize, cudaMemcpyDeviceToHost));
// Free Device memory
printf("Free Device memoryn");
checkCudaErrors(cudaFree(d_idata));
checkCudaErrors(cudaFree(d_odata));
// Save Picture
printf("Save Picturen");
bool saved = false;
if (channels==1)
saved = savePPM("output.pgm", data, w, h, channels);
else if (channels==3)
saved = savePPM("output.ppm", data, w, h, channels);
else fprintf(stderr, "ERROR: Unable to save file - wrong channel!n");
// Free Host memory
printf("Free Host memoryn");
free(data);
if (!saved){
fprintf(stderr, "Failed to save Filen");
exit(EXIT_FAILURE);
}
printf("Donen");
}
Something is wrong with the filter function. The loadPPM and savePPM (part of the cuda samples) are working with an other kernel function, but with this filterfunction I get an black image.
So the question is: What did I wrong?
Some other comprehension questions: Here https://www.nvidia.com/docs/IO/116711/sc11-cuda-c-basics.pdf I read that threads can only communicate within a block (shared memory, syncthreads, ..). So in my function the image is split into rectangular blocks and the picture on page 9 of the Image Processing slides is about one block? What about the pixels at the edge of a block? Are they unchanged?
Thanks for your answers.
One problem in your code is that your kernel is expecting a 2D grid and 2D threadblocks:
int x = blockIdx.x*TILE_W + threadIdx.x - R;
int y = blockIdx.y*TILE_H + threadIdx.y - R;
^^^^^^^^^^ ^^^^^^^^^^^
2D grid 2D threadblock
but you are launching a kernel with a 1D grid and threadblock definition:
int threadsPerBlock = 256; // creates 1D threadblock
int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; //1D grid
....
d_filter<<<blocksPerGrid, threadsPerBlock>>>(d_idata, d_odata, w,h);
So when you launch that kernel, threadIdx.y will always be zero as will blockIdx.y
When I make a modified version of your code that does not depend on PPM image load/store (so, using synthetic data), and make the changes necessary to launch a 2D grid and threadblock, to be consistent with your kernel, the code seems to run correctly for me and produce output that seems like it might be filtered output, instead of zeros:
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#define TILE_W 16
#define TILE_H 16
#define R 2 // filter radius
#define D (R*2+1) // filter diameter
#define S (D*D) // filter size
#define BLOCK_W (TILE_W+(2*R))
#define BLOCK_H (TILE_H+(2*R))
__global__ void d_filter(unsigned char *g_idata, unsigned char *g_odata, unsigned int width, unsigned int height)
{
__shared__ unsigned char smem[BLOCK_W*BLOCK_H];
int x = blockIdx.x*TILE_W + threadIdx.x - R;
int y = blockIdx.y*TILE_H + threadIdx.y - R;
// clamp to edge of image
x = max(0, x);
x = min(x, width-1);
y = max(y, 0);
y = min(y, height-1);
unsigned int index = y*width + x;
unsigned int bindex = threadIdx.y*blockDim.y+threadIdx.x;
// each thread copies its pixel of the block to shared memory
smem[bindex] = g_idata[index];
__syncthreads();
// only threads inside the apron will write results
if ((threadIdx.x >= R) && (threadIdx.x < (BLOCK_W-R)) && (threadIdx.y >= R) && (threadIdx.y < (BLOCK_H-R))) {
float sum = 0;
for(int dy=-R; dy<=R; dy++) {
for(int dx=-R; dx<=R; dx++) {
float i = smem[bindex + (dy*blockDim.x) + dx];
sum += i;
}
}
g_odata[index] = sum / S;
}
}
const unsigned int imgw = 512;
const unsigned int imgh = 256;
void loadImg(unsigned char **data, unsigned int *w, unsigned int *h, unsigned int *ch){
*w = imgw;
*h = imgh;
*ch = 1;
*data = (unsigned char *)malloc(imgw*imgh*sizeof(unsigned char));
for (int i = 0; i < imgw*imgh; i++) (*data)[i] = i%8;
}
int main(){
unsigned char *data=NULL, *d_idata=NULL, *d_odata;
unsigned int w,h,channels;
loadImg(&data, &w, &h, &channels);
printf("Loaded file with w:%d h:%d channels:%d n",w,h,channels);
unsigned int numElements = w*h*channels;
size_t datasize = numElements * sizeof(unsigned char);
cudaError_t err = cudaSuccess;
// Allocate the Device Memory
printf("Allocate Devicememory for datan");
err = cudaMalloc((void **)&d_idata, datasize);
if ( err != cudaSuccess)
{
fprintf(stderr, "Failed to allocate device memory for idata (error code %s)!n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
err = cudaMalloc((void **)&d_odata, datasize);
if ( err != cudaSuccess)
{
fprintf(stderr, "Failed to allocate device memory for odata (error code %s)!n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
// Copy to device
printf("Copy idata from the host memory to the CUDA devicen");
err =cudaMemcpy(d_idata, data, datasize, cudaMemcpyHostToDevice);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to copy idata from host to device (error code %s)!n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
// Launch Kernel
dim3 threadsPerBlock(BLOCK_W, BLOCK_H);
dim3 blocksPerGrid((w+threadsPerBlock.x-1)/threadsPerBlock.x, (h+threadsPerBlock.y-1)/threadsPerBlock.y);
printf("CUDA kernel launch with %d,%d blocks of %d,%d threadsn", blocksPerGrid.x, blocksPerGrid.y, threadsPerBlock.x, threadsPerBlock.y);
d_filter<<<blocksPerGrid, threadsPerBlock>>>(d_idata, d_odata, w,h);
err=cudaGetLastError();
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to launch kernel (error code %s)!n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
// Copy data from device to host
printf("Copy odata from the CUDA device to the host memoryn");
err=cudaMemcpy(data, d_odata, datasize, cudaMemcpyDeviceToHost);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to copy odata from device to host (error code %s)!n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
// Free Device memory
printf("Free Device memoryn");
err=cudaFree(d_idata);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to free device idata (error code %s)!n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
err=cudaFree(d_odata);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to free device odata (error code %s)!n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
printf("results:n");
for (int i = 0; i < 16; i++){
for (int j = 0; j < 16; j++) printf("%d ", data[i*w+j]);
printf("n");}
// Free Host memory
printf("Free Host memoryn");
free(data);
printf("nDonen");
}
When I run the above code with cuda-memcheck , I get this:
C:ProgramDataNVIDIA CorporationCUDA Samplesv5.0binwin32Debug>cuda-memcheck test
========= CUDA-MEMCHECK
Loaded file with w:512 h:256 channels:1
Allocate Devicememory for data
Copy idata from the host memory to the CUDA device
CUDA kernel launch with 26,13 blocks of 20,20 threads
Copy odata from the CUDA device to the host memory
Free Device memory
results:
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
0 1 2 3 4 5 4 3 3 2 2 3 4 5 4 3
Free Host memory
Done
========= ERROR SUMMARY: 0 errors
C:ProgramDataNVIDIA CorporationCUDA Samplesv5.0binwin32Debug>
上一篇: GPU blob边框连接组件标签
下一篇: 带有cuda的盒式过滤器c