The purpose of this code is to demonstrate how Cython can be used to reclaim much of the performance lost when using naive Python code. To demonstrate this, we write a naive image processing loop in Python. Then, we write the same thing using Cython. By running each cell separately, we see that the naive Python code takes about 30 seconds to run, whereas the Cython code runs almost instantly.
#
# this line of code is neccesary when writing Cython code
#
%load_ext cythonmagic
import numpy
def python_convolve(f, g):
if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1:
raise ValueError("Only odd dimensions on filter supported")
vmax = f.shape[0]
wmax = f.shape[1]
smax = g.shape[0]
tmax = g.shape[1]
smid = smax // 2
tmid = tmax // 2
xmax = vmax + 2*smid
ymax = wmax + 2*tmid
h = numpy.zeros([xmax, ymax], dtype=f.dtype)
for x in range(xmax):
for y in range(ymax):
s_from = max(smid - x, -smid)
s_to = min((xmax - x) - smid, smid + 1)
t_from = max(tmid - y, -tmid)
t_to = min((ymax - y) - tmid, tmid + 1)
value = 0
for s in range(s_from, s_to):
for t in range(t_from, t_to):
v = x - smid + s
w = y - tmid + t
value += g[smid - s, tmid - t] * f[v, w]
h[x, y] = value
return h
%%cython
#
# the line above is neccesary when writing Cython code (and must come before any other code or comments)
#
import numpy
#
# begin cython boilerplate code
#
cimport numpy
FLOAT32_DTYPE = numpy.float32
ctypedef numpy.float32_t FLOAT32_DTYPE_t
cdef inline int _max(int a, int b): return a if a >= b else b
cdef inline int _min(int a, int b): return a if a <= b else b
#
# end cython boilerplate code
#
#
# note the explicit type declarations in the function signature
#
def cython_convolve(numpy.ndarray[FLOAT32_DTYPE_t, ndim=2] f, numpy.ndarray[FLOAT32_DTYPE_t, ndim=2] g):
if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1:
raise ValueError("Only odd dimensions on filter supported")
assert f.dtype == FLOAT32_DTYPE and g.dtype == FLOAT32_DTYPE
#
# note the explicit type declarations
#
cdef int vmax = f.shape[0]
cdef int wmax = f.shape[1]
cdef int smax = g.shape[0]
cdef int tmax = g.shape[1]
cdef int smid = smax // 2
cdef int tmid = tmax // 2
cdef int xmax = vmax + 2*smid
cdef int ymax = wmax + 2*tmid
cdef numpy.ndarray[FLOAT32_DTYPE_t, ndim=2] h = numpy.zeros([xmax, ymax], dtype=FLOAT32_DTYPE)
cdef int x, y, s, t, v, w, s_from, s_to, t_from, t_to
cdef FLOAT32_DTYPE_t value
for x in range(xmax):
for y in range(ymax):
s_from = _max(smid - x, -smid)
s_to = _min((xmax - x) - smid, smid + 1)
t_from = _max(tmid - y, -tmid)
t_to = _min((ymax - y) - tmid, tmid + 1)
value = 0
for s in range(s_from, s_to):
for t in range(t_from, t_to):
v = x - smid + s
w = y - tmid + t
value += g[smid - s, tmid - t] * f[v, w]
h[x, y] = value
return h
import scipy
gaussian_blur_kernel_width = numpy.int32(9)
gaussian_blur_kernel_half_width = numpy.int32(4)
gaussian_blur_sigma = numpy.float32(2)
y, x = \
scipy.mgrid[-gaussian_blur_kernel_half_width:gaussian_blur_kernel_half_width+1,-gaussian_blur_kernel_half_width:gaussian_blur_kernel_half_width+1]
gaussian_blur_kernel_not_normalized = numpy.exp( ( - ( x**2 + y**2 ) ) / ( 2 * gaussian_blur_sigma**2 ) )
normalization_constant = numpy.float32(1) / numpy.sum(gaussian_blur_kernel_not_normalized)
gaussian_blur_kernel = (normalization_constant * gaussian_blur_kernel_not_normalized).astype(numpy.float32)
imshow(gaussian_blur_kernel, cmap="gray", interpolation="nearest");
title("gaussian_blur_kernel");
colorbar();
import skimage
import skimage.data
img = numpy.array(skimage.data.coins(), dtype=numpy.float32)
imshow(img, cmap="gray", interpolation="nearest");
title("img");
colorbar();
python_result = python_convolve(img, gaussian_blur_kernel)
imshow(python_result, cmap="gray", interpolation="nearest");
cython_result = cython_convolve(img, gaussian_blur_kernel)
imshow(cython_result, cmap="gray", interpolation="nearest");
diff = abs( cython_result.astype(float32) - python_result.astype(float32) );
imshow(diff, cmap="gray", interpolation="nearest");
colorbar();