Image convolution: direct space versus Fourier space

In image processing, the convolution with a small matrix is used for blurring, sharpening, embossing, edge-detection, and more. This is accomplished by means of convolution between a kernel and an image. http://en.wikipedia.org/wiki/Kernel_%28image_processing%29

In this exercise we will focus on the gaussian blur and compare the execution time of two appoaches: direct space vs Fourier space for different width of gaussians: 1, 2, 4, 8, and 16pixels.

The raw image used is the famous 512x512 pixel image from Lena, and available as part of scipy.



In [1]:

    
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib



In [2]:

    
from scipy.misc import lena
img = lena()
imshow(img, cmap="gray")









    Out[2]:





<matplotlib.image.AxesImage at 0x7f665144bc10>

Properties of the gaussian for bluring an image



In [3]:

    
from IPython.display import Image
Image(filename="files/gaussian.png")









    Out[3]:

Very smooth bluring
Separable in x/y (hence much faster)
The Fourier transform is a gaussian of width proportional to 1/sigma



In [4]:

    
from scipy.signal import gaussian
plot(gaussian(21,1))









    Out[4]:





[<matplotlib.lines.Line2D at 0x7f6641a575d0>]

Direct space approach

Choice of the window size ?



In [5]:

    
(gaussian(21,1)>1e-4).sum()









    Out[5]:





9

It is important to have an odd window size, we choose w = 8 sigma+1.

Scipy.signal has a special function for direct space convolution: sepfir2d https://docs.scipy.org/doc/scipy-0.7.x/reference/generated/scipy.signal.sepfir2d.html



In [6]:

    
from scipy.signal import sepfir2d
def direct(img, sigma):
    gaus = gaussian(8*sigma+1, sigma)
    return sepfir2d(img, gaus, gaus)



In [7]:

    
imshow(direct(img,5),cmap="gray")









    Out[7]:





<matplotlib.image.AxesImage at 0x7f66419efdd0>



In [8]:

    
%timeit direct(img,2)









    



100 loops, best of 3: 14.1 ms per loop

Fourier space convolution

The convolution in direct space is equivalent to a multiplication in Fourier space. http://en.wikipedia.org/wiki/Convolution_theorem



In [9]:

    
from IPython.display import Image
Image(filename="files/convolution.png")









    Out[9]:

This is especially appealing for gaussian blurring as the FT of a gaussian is a gaussian. The total cost is reduced to one FT and its inverse.



In [10]:

    
from scipy.fftpack import fft2, ifft2, fftshift
def ft(img, sigma):
    size = len(img)
    gaus = gaussian(size, size / (2*pi*sigma))
    sgau = fftshift(gaus)
    return ifft2(fft2(img) * outer(sgau,sgau)).real



In [11]:

    
imshow(ft(img,5),cmap="gray")









    Out[11]:





<matplotlib.image.AxesImage at 0x7f664191e390>



In [12]:

    
%timeit ft(img,2)









    



10 loops, best of 3: 30 ms per loop

Benchmark:



In [13]:

    
import time
sigmas = []
d = []
f = []
for i in range(5):
    sigma = 1<<i
    t0 = time.time()
    direct(img, sigma)
    t1 = time.time()
    ft(img, sigma)
    t2 = time.time()
    
    print("sigma=%s \t Direct: %.3fs \t Fourier: %.3fs"%(sigma,t1-t0,t2-t1))
    sigmas.append(sigma)
    d.append(t1-t0)
    f.append(t2-t1)
plot(sigmas,d,label="Direct")
plot(sigmas,f,label="Fourier")
legend()









    



sigma=1 	 Direct: 0.007s 	 Fourier: 0.035s
sigma=2 	 Direct: 0.012s 	 Fourier: 0.035s
sigma=4 	 Direct: 0.027s 	 Fourier: 0.037s
sigma=8 	 Direct: 0.053s 	 Fourier: 0.036s
sigma=16 	 Direct: 0.102s 	 Fourier: 0.032s






    Out[13]:





<matplotlib.legend.Legend at 0x7f6641853910>

Nota

The gaussian blur filter exists also as "ready to use" function as part of scipy.ndimage.filters.gaussian_filter. We can compare how it behaves compared to our implementations:



In [14]:

    
from scipy.ndimage import filters



In [15]:

    
%timeit filters.gaussian_filter(img, 2)









    



100 loops, best of 3: 6.9 ms per loop



In [16]:

    
%timeit filters.gaussian_filter(img, 4)









    



100 loops, best of 3: 10.3 ms per loop



In [17]:

    
import time
from scipy.ndimage import filters
sigmas = []
d = []
f = []
n = []
for i in range(5):
    sigma = 1<<i
    t0 = time.time()
    direct(img, sigma)
    t1 = time.time()
    ft(img, sigma)
    t2 = time.time()
    filters.gaussian_filter(img, sigma)
    t3=time.time()
    print("sigma=%s \t Direct: %.3fs \t Fourier: %.3fs, ndimage: %.3fs"%(sigma,t1-t0,t2-t1, t3-t2))
    sigmas.append(sigma)
    d.append(t1-t0)
    f.append(t2-t1)
    n.append(t3-t2)
plot(sigmas,d,label="Direct")
plot(sigmas,f,label="Fourier")
plot(sigmas,n,label="Ndimage")
legend()









    



sigma=1 	 Direct: 0.006s 	 Fourier: 0.035s, ndimage: 0.005s
sigma=2 	 Direct: 0.010s 	 Fourier: 0.027s, ndimage: 0.007s
sigma=4 	 Direct: 0.019s 	 Fourier: 0.027s, ndimage: 0.011s
sigma=8 	 Direct: 0.049s 	 Fourier: 0.038s, ndimage: 0.022s
sigma=16 	 Direct: 0.106s 	 Fourier: 0.043s, ndimage: 0.046s






    Out[17]:





<matplotlib.legend.Legend at 0x7f6641795f90>



In [17]: