<jupyter_start><jupyter_text>See <NAME>'s blog entries ([Part3](https://www.vincent-lunot.com/post/an-introduction-to-cuda-in-python-part-3/), [Part5](https://www.vincent-lunot.com/post/an-introduction-to-cuda-in-python-part-5/)) on image convolution for more detail.<jupyter_code>from numba import cuda
import numpy as np
import skimage.data
from skimage.color import rgb2gray
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Reference image<jupyter_code>full_image = rgb2gray(skimage.data.coffee()).astype(np.float32) / 255
plt.figure()
plt.imshow(full_image, cmap='gray')
plt.title("Full size image:")
image = full_image[150:350, 200:400].copy() # We don't want a view but an array and therefore use copy()
plt.figure()
plt.imshow(image, cmap='gray')
plt.title("Part of the image we use:")
plt.show()<jupyter_output><empty_output><jupyter_text>Convolution Kernel<jupyter_code>@cuda.jit
def convolve(result, mask, image):
    # expects a 2D grid and 2D blocks,
    # a mask with odd numbers of rows and columns, (-1-) 
    # a grayscale image
    
    # (-2-) 2D coordinates of the current thread:
    i, j = cuda.grid(2) 
    
    # (-3-) if the thread coordinates are outside of the image, we ignore the thread:
    image_rows, image_cols = image.shape
    if (i >= image_rows) or (j >= image_cols): 
        return
    
    # To compute the result at coordinates (i, j), we need to use delta_rows rows of the image 
    # before and after the i_th row, 
    # as well as delta_cols columns of the image before and after the j_th column:
    delta_rows = mask.shape[0] // 2 
    delta_cols = mask.shape[1] // 2
    
    # The result at coordinates (i, j) is equal to 
    # sum_{k, l} mask[k, l] * image[i - k + delta_rows, j - l + delta_cols]
    # with k and l going through the whole mask array:
    s = 0
    for k in range(mask.shape[0]):
        for l in range(mask.shape[1]):
            i_k = i - k + delta_rows
            j_l = j - l + delta_cols
            # (-4-) Check if (i_k, j_k) coordinates are inside the image: 
            if (i_k >= 0) and (i_k < image_rows) and (j_l >= 0) and (j_l < image_cols):  
                s += mask[k, l] * image[i_k, j_l]
    result[i, j] = s<jupyter_output><empty_output><jupyter_text>Run the kernel<jupyter_code># We preallocate the result array:
result = np.empty_like(image)

# We choose a random mask:
mask = np.random.rand(13, 13).astype(np.float32) 
mask /= mask.sum()  # We normalize the mask
print('Mask shape:', mask.shape)
print('Mask first (3, 3) elements:\n', mask[:3, :3])

# We use blocks of 32x32 pixels:
blockdim = (32, 32)
print('Blocks dimensions:', blockdim)

# We compute grid dimensions big enough to cover the whole image:
griddim = (image.shape[0] // blockdim[0] + 1, image.shape[1] // blockdim[1] + 1)
print('Grid dimensions:', griddim)

# We apply our convolution to our image:
convolve[griddim, blockdim](result, mask, image)

# We plot the result:
plt.figure()
plt.imshow(image, cmap='gray')
plt.title("Before convolution:")
plt.figure()
plt.imshow(result, cmap='gray')
plt.title("After convolution:")
plt.show()<jupyter_output><empty_output><jupyter_text>Benchmark / compareOur custom numba kernel vs one in `scipy`  Error<jupyter_code>from scipy.ndimage.filters import convolve as scipy_convolve

scipy_result = scipy_convolve(image, mask, mode='constant', cval=0.0, origin=0)
print('Maximum relative error:', np.max(np.abs(result - scipy_result) / np.abs(scipy_result)))<jupyter_output><empty_output><jupyter_text>Timing<jupyter_code>%timeit convolve[griddim, blockdim](result, mask, image)

scipy_result = np.empty_like(image)
%timeit scipy_convolve(image, mask, output=scipy_result, mode='constant', cval=0.0, origin=0)
d_image = cuda.to_device(image)
d_mask = cuda.to_device(mask)
d_result = cuda.to_device(result)
%timeit convolve[griddim, blockdim](d_result, d_mask, d_image)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13172588832487309, 'fraction_numerical': 0.017766497461928935, 'mean_word_length': 4.226790450928382, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 30, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '373268', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 1493, 'n_tokens_neox': 1408, 'n_words': 386, 'path': '4_numba_image_convolution.ipynb', 'repo_name': 'geoffwoollard/pydata-amsterdam2019-numba', 'size': '3947'}