Image convolution at specific points

Question

Is there a way in scipy (or other similar library) to get the convolution of an image with a given kernel only at some desired points?

I'm looking for something like:

ndimage.convolve(image, kernel, mask=mask)

Where mask contains True (or 1) whenever the kernel needs to be applied, False (or 0) otherwise.

EDIT: Example python code that does what I'm trying to do (but not faster than a whole image convolution using scipy):

def kernel_responses(im, kernel, mask=None, flatten=True):
    if mask is None:
        mask = np.ones(im.shape[:2], dtype=np.bool)

    ks = kernel.shape[0]//2

    data = np.pad(im, ks, mode='reflect')
    y, x = np.where(mask)

    responses = np.empty(y.shape[0], float)

    for k, (i, j) in enumerate(zip(y, x)):
        responses[k] = (data[i:i+ks*2+1, j:j+ks*2+1] * kernel).sum()

    if flatten:
        return responses

    result = np.zeros(im.shape[:2], dtype=float)
    result[y, x] = responses
    return result

The above code does the job with a wrap boundary conditions, but the inner loop is in python, and thus, slow. I was wondering if there is something faster already implemented in scipy/opencv/skimage.

Answer 1

I know that I'm responding to my own answer, I hope the code bellow carries further improvements, or it might be useful for other users.

The code bellow is a cython/python implementation:

PYTHON:

def py_convolve(im, kernel, points):
    ks = kernel.shape[0]//2
    data = np.pad(im, ks, mode='constant', constant_values=0)
    return cy_convolve(data, kernel, points)

CYTHON:

import numpy as np
cimport cython

@cython.boundscheck(False)
def cy_convolve(unsigned char[:, ::1] im, double[:, ::1] kernel, Py_ssize_t[:, ::1] points):
    cdef Py_ssize_t i, j, y, x, n, ks = kernel.shape[0]
    cdef Py_ssize_t npoints = points.shape[0]
    cdef double[::1] responses = np.zeros(npoints, dtype='f8')

    for n in range(npoints):
        y = points[n, 0]
        x = points[n, 1]
        for i in range(ks):
            for j in range(ks):
                responses[n] += im[y+i, x+j] * kernel[i, j]

     return np.asarray(responses)

Comparision with other methods

The following tables shows evaluation of 4 methods:

1. My python method in the question
2. The method from @Vighnesh Birodkar
3. Complete image convolution with scipy
4. My python/cython implementation in this post

Each rows, in order, correspond to those methods for 3 different images (coins, camera and lena from skimage.data respectively) and each of the columns corresponds to a different ammount of points to calculate the kernel responses (is in percentages as meaning "calculate response in x% of the points of the image").

For calculating the kernel response in less than 50% of the points, my implementation is faster than convolving the whole image, but is not faster otherwise..

EDIT: kernel windows for the tests are 5x5 uniform windows (np.ones((5,5))).

['303x384']    1%     2%     5%    10%     20%     50%
1            4.97   9.58  24.32  48.28  100.39  245.77
2            7.60  15.09  37.42  75.17  150.09  375.60
3            3.05   2.99   3.04   2.88    2.96    2.98
4            0.17   0.22   0.38   0.60    1.10    2.49

['512x512']     1%     2%     5%     10%     20%     50%
1            10.68  21.87  55.47  109.16  223.58  543.73
2            17.90  34.59  86.02  171.20  345.46  858.24
3             6.52   6.53   6.74    6.63    6.43    6.60
4             0.31   0.43   0.78    1.34    2.73    6.82

['512x512']     1%     2%     5%     10%     20%     50%
1            13.21  21.45  54.98  110.80  217.11  554.96
2            19.55  34.78  87.09  172.33  344.58  893.02
3             6.87   6.82   7.00    6.60    6.64    7.71
4             0.35   0.47   0.87    1.57    2.47    6.07

NOTE: times are in ms.

Answer 2

I don't know of any function that does exactly what you're asking. If instead of providing a mask of points to be convolved you provided a list of points ex. [(7, 7), (100, 100)] then it might be as simple as getting the appropriate image patch (say the same size as your provided kernel), convolve the image patch and kernel, and insert back into the original image.

Here's a coded example, hopefully it's close enough for you to modify lightly:

[EDIT: I noticed a couple errors I had in my padding and patch arithmetic. Previously, you could not convolve with a point right on the boarder (say (0, 0)), I doubled the padding, fixed some arithmetic, and now all is well.]

import cv2
import numpy as np
from scipy import ndimage
from matplotlib import pyplot as plt

def image_convolve_mask(image, list_points, kernel):
# list_points ex. [(7, 7), (100, 100)]
# assuming kernels of dims 2n+1 x 2n+1
rows, cols = image.shape
k_rows, k_cols = kernel.shape
r_pad = int(k_rows/2)
c_pad = int(k_cols/2)
# zero-pad the image in case desired point is close to border
padded_image = np.zeros((rows + 2*k_rows, cols + 2*k_cols))
# set the original image in the center
padded_image[k_rows: rows + k_rows, k_cols: cols + k_cols] = image
# should you prefer to use np.pad:
# padded_image = np.pad(image, (k_rows, k_cols), 'constant', constant_values=(0, 0))

for p in list_points:
    # extract pertinent patch from image
    # arbitrarily choosing the patch as same size as the kernel; change as needed
    patch = padded_image[p[0] + k_rows - r_pad: p[0] + 2*k_rows - r_pad, p[1] + k_cols - c_pad: p[1] + 2*k_cols - c_pad]

    # here use whatever function for convolution; I prefer cv2filter2D()
    # commented out is another option
    # conv = ndimage.convolve(patch, kernel, mode='constant', cval=0.0)
    conv = cv2.filter2D(patch, -1, kernel)
    # set the convolved patch back in to the image
    padded_image[p[0] + k_rows - r_pad: p[0] + 2*k_rows - r_pad, p[1] + k_cols - c_pad: p[1] + 2*k_cols - c_pad] = conv

return padded_image[k_rows: rows + k_rows, k_cols: cols + k_cols]

Now to try it out on an image:

penguins = cv2.imread('penguins.png', 0)
kernel = np.ones((5,5),np.float32)/25
# kernel = np.array([[-1, -1, -1], [-1, 8, -1], [-1, -1, -1]], np.float32)
conv_image = image_convolve_mask(penguins, [(7, 7), (36, 192), (48, 207)], kernel)
plt.imshow(conv_image, cmap = 'gray', interpolation = 'bicubic')
plt.xticks([]), plt.yticks([])
plt.show()

I applied a 5x5 box smoother and can't really see any change around pixel (7, 7), but I chose the other two points to be the tips of the two left-most penguin's beaks. So you can see the smoothed patches.

Here is a lena512 image with 21 convolution points (time:0.006177 sec).

[EDIT 2: An example of using a mask to generate a list of row, col tuples to feed in to the function.]

mask = np.eye(512)
k = np.ones((25, 25), np.float32)/625
list_mask = zip(np.where(mask==1)[0], np.where(mask==1)[1])
tic = time.time()
conv_image = image_convolve_mask(lena, list_mask, k)
print 'time: ', time.time()-tic # 0.08136 sec