View source for PIL Cuda

{{GPU610/DPS915 Index | 20143}}
= PIL Cuda=
== Team Members == 
# [mailto:gcastrolondono@senecacollege.ca?subject=gpu610 Gabriel Castro]
== Pillow and Guassian Blur in CUDA ==

=== Pillow ===

[https://python-pillow.github.io/ Pillow] is an imaging processing is a a library written in python with much of the processing done in C.
Pillow has many filters that can be applied to an image to change or distort it.
For this project I chose to optimize the Gaussian Blur filter,
because after some profiling I determined that it is the single filter with the highest run time.

==== sample ====

 #!/usr/bin/env python
 from PIL import Image, ImageFilter
 im = Image.open("test.png")
 filtered = im.filter(ImageFilter.GaussianBlur(5))
 filtered.save("out.png")

===The Gaussian Blur===

The Gaussian Blur filter is an algorithm applied to images to make them blurry and/or reducing noise.
The algorithm works on every pixel and colour channel (RGBA) in the image individually.
For any given pixel (x,y) take the weighted average of all the pixels in a given radius (r).

Given an image of size `x * y` with `c` channels and a blur radius of `r` pixels.

There are `x * y * c * (r * r + 1)` floating point operations.

This means a change in any dimension will result in an exponential increase.
 ex. a 800x600 RGB image with a 5px blur radius will have 600*800*3*(5*5+1) = 37,440,000 floating point ops

====optimization 1====

We can actually reduce the number of radius based operations by process sing the image in 2 steps.

1. Blur every pixel but only along the `x` axis, and store in a temporary buffer.

2. Blur every pixel in the temporary buffer along the `y` axis, store in result buffer. 



===CUDA with Python setup tools===

Since Pillow is a python library that uses [http://cython.org/ cython] built via setup tools we need to find a way to inject the nvcc compiler for CUDA. I was able to get everything working by modifying a script I found online, see the final version [https://github.com/GabrielCastro/Pillow/blob/master/setup_cuda.py here]

and adding the this code to the existing build script


 _LIB_IMAGING_CUDA = ["UnsharpMaskCuda"]
 have_cuda = False
 CUDA = None
 try:
     import setup_cuda
     CUDA = setup_cuda.CUDA
     have_cuda = True
 except EnvironmentError as e:
     print('CUDA not found')
     print(e)
 if have_cuda:
     defs.append(("HAVE_CUDA", 1))
     setup_cuda.customize_compiler_for_nvcc(self.compiler)
     for f in _LIB_IMAGING_CUDA:
         files.append('libImaging/' + f + '.cu')
     libs.append('cudart')
     libs.append('stdc++')
     _add_directory(self.compiler.include_dirs, 'libImaging/')


 NOTE.
 This requires $CUDAHOME to be set to /usr/local/cuda or appropriate
  and PATH to include $CUDAHOME/bin
  and LD_LIBRARY_PATH to include $CUDAHOME/lib64
  * tested on Ubuntu 14.04 with CUDA 6.5


===Pillow C Api===

The Pillow library has a C api that can for the purposes of this post be described with the following code.

 struct _Imaing {
     int xsize;
     int ysize;
     // y-sized array of pointers to x-sized arrays of (UINT8) pixel data
     UINT8** image8;
 }
 typedef _Imaging* Imaging;
 Imaging gblur(Imaging im, Imaging imOut, float floatRadius);

 This is rather simplified, as the `Imaging` struct contains many more properties and we're taking advantage of the fact that the library separates an incoming image into separate `Imaging` structs for each colour channel.

=== Gaussian Blur in CUDA===
While all the code is available on github https://github.com/GabrielCastro/Pillow this is a simplified psudo code version of the blur process

 __global__ static void blurRows(const px8_t* __restrict__ in, float* __restrict__ buff,
     const size_t xSize, const size_t ySize,
     const float* __restrict__ mask,  const size_t radius) {
 
     const size_t y = blockIdx.y * blockDim.y + threadIdx.y;
     const size_t x = blockIdx.x * blockDim.x + threadIdx.x;
     if (y >= ySize || x >= xSize) {
         return;
     }
     const size_t buffIdx = y * xSize + x;
 
     float sum = 0;
 
     for (size_t p = 0; p < radius; ++p) {
         float maskVal = mask[p];
         int offset = (int)(-((float)radius / 2.0) + (float)p + 0.5);
         int xOff = x + offset;
         if (xOff < 0) {
             offset = -x;
         } else if (xOff >= xSize) {
             offset = xSize - x - 1;
         }
         size_t pxIndex = buffIdx + offset;
         sum += in[pxIndex] * maskVal;
     }
 
     buff[buffIdx] = sum;
 }


 __global__ static void blurCols(const float* __restrict__ buff, px8_t* __restrict__ out,
     const size_t xSize, const size_t ySize,
     const float* __restrict__ mask,  const size_t radius) {
 
     const size_t y = blockIdx.y * blockDim.y + threadIdx.y;
     const size_t x = blockIdx.x * blockDim.x + threadIdx.x;
     if (y >= ySize || x >= xSize) {
         return;
     }
 
     size_t outIdx = y * xSize + x;
     float sum = 0;
 
     for (size_t p = 0; p < radius; ++p) {
         float maskVal = mask[p];
         int offset = (int)(-((float)radius / 2.0) + (float)p + 0.5);
         int lOff = y + offset;
         if (lOff < 0) {
             offset = -y;
         } else if (lOff >= ySize) {
             offset = ySize - y - 1;
         }
 
         size_t buffIdx = outIdx + offset * xSize;
         sum += buff[buffIdx] * maskVal;
     }
     out[outIdx] = (px8_t)CLIP(sum);
 }
 


 gblur(Imaging in, Imagin out, int radius) {
     UINT8* d_img =  // allocate in->xsize * in->ysize on device
     float* d_buff = // allocate a tmp buffer on device of size in->ysize * in->xsize
     float* d_mask = // allocate and create a guassian blur mask on device
 
     // copy the image into the device
     size_t rowSize = in->xsize * sizeof(UINT8);
     for (int i = 0; i < in->ysize; ++i) {
         cudaMemcpyAsync(d_img + i * rowSize, in->image8[i], rowSize, cudaMemcpyHostToDevice, 0);
     }
     size_t xBlocks = (in->xsize + ntpb - 1) / ntpb;
     size_t yBlocks = (in->ysize + ntpb - 1) / ntpb;
 
     dim3 grid(xBlocks, yBlocks);
     dim3 block(ntpb, ntpb);
 
     blurRows<<<grid,block>>>(d_img, d_buff, in->xsize, in->ysize, d_mask, radius);
     blurCols<<<grid,block>>>(d_buff, d_img, in->xsize, in->ysize, d_mask, radius);
 
     for (int i = 0; i < in->ysize; ++i) {
         cudaMemcpyAsync(out->image8[i], d_img + i * rowSize, rowSize, cudaMemcpyDeviceToHost, 0);
     }
     cudaDeviceSynchronize();
     cudaFree(d_img);
     cudaFree(d_buff);
     cudaFree((void*) d_mask);
 }