Improving performance iterating in 2d numpy array

Question

I have two 2d numpy array (images). First one defined by image, is storing the sum of a movement at the pixel (i,j)

Second one define by nbCameras, is storing the number of cameras who can see a movement at this pixel (i,j)

I want to create a third image imgFinal which only store the value of the pixel (i,j) and it's neighbours (3 x 3) mask, if the number of cameras who can see the pixel (i,j) is greater than 1.

For now I'm using two for loops which is not the best way. I'd like to increase the speed of the computation but I didn't find the best way to do it yet. Also I'm a bit blocked as the fact I want to converse the neighbours of the pixel (i, j)

I also tried to use bumpy.vectorize but i can keep the neighbours of my pixel in this case.

What would be the best way to increase the speed of this function?

Thanks for your help!

maskWidth = 3
dstCenterMask = int( (maskWidth - 1) / 2)

imgFinal = np.zeros((image.shape),dtype = np.float32) 

 for j in range(dstCenterMask,image.shape[0] - dstCenterMask):
    for i in range(dstCenterMask,image.shape[1] - dstCenterMask):
        if nbCameras[j,i] > 1
           imgFinal[j - dstCenterMask : j + dstCenterMask + 1, i - dstCenterMask : i +     dstCenterMask + 1] =
           image[j - dstCenterMask : j + dstCenterMask + 1, i - dstCenterMask : i + dstCenterMask + 1]

Answer 1

Option 1: Try to rewrite the code in a vectorized way. You could convolve with a 3x3 mask like this:

import numpy as np
from scipy.signal import convolve2d

image = np.random.random((100,100))
nbCameras = np.abs(np.random.normal(size=(100,100)).round())

maskWidth = 3
mask = np.ones((maskWidth, maskWidth))

visibilityMask = (nbCameras>1).astype(np.float)
visibilityMask = convolve2d(visibilityMask, mask, mode="same").astype(np.bool)

imgFinal = image.copy()
imgFinal[~visibilityMask] *= 0

import matplotlib.pyplot as plt
for i, (im, title) in enumerate([(image, "image"), 
                                 (nbCameras, "nbCameras"), 
                                 (visibilityMask, "visibilityMask"), 
                                 (imgFinal, "imgFinal")]):
    plt.subplot(2,2,i+1)
    plt.title(title)
    plt.imshow(im, cmap=plt.cm.gray)

plt.show()

This will result in this plot:

Option 2: Use Numba. This uses an advanced just-in-time optimization technique and is specifically useful for speeding up loops.

Answer 2

This doesn't handle cameras on the edge of the array, but neither does your code:

import numpy as np
from numpy.lib.stride_tricks import as_strided

rows, cols, mask_width = 10, 10, 3
mask_radius = mask_width // 2

image = np.random.rand(rows, cols)
nb_cameras = np.random.randint(3 ,size=(rows, cols))

image_view = as_strided(image, shape=image.shape + (mask_width, mask_width),
                        strides=image.strides*2)
img_final = np.zeros_like(image)
img_final_view = as_strided(img_final,
                            shape=img_final.shape + (mask_width, mask_width),
                            strides=img_final.strides*2)
copy_mask = nb_cameras[mask_radius:-mask_radius,
                       mask_radius:-mask_radius] > 1

img_final_view[copy_mask] = image_view[copy_mask]

After running the above code:

>>> nb_cameras
array([[0, 2, 1, 0, 2, 0, 1, 2, 1, 0],
       [0, 1, 1, 1, 1, 2, 1, 1, 2, 1],
       [1, 2, 2, 2, 1, 2, 1, 0, 2, 0],
       [0, 2, 2, 0, 1, 2, 1, 0, 1, 0],
       [1, 2, 0, 1, 2, 0, 1, 0, 0, 2],
       [2, 0, 1, 1, 1, 1, 1, 1, 0, 1],
       [1, 0, 2, 2, 0, 1, 1, 1, 0, 0],
       [0, 0, 1, 0, 1, 0, 1, 0, 2, 2],
       [0, 1, 0, 1, 1, 2, 1, 1, 2, 2],
       [2, 2, 0, 1, 0, 0, 1, 2, 1, 0]])
>>> np.round(img_final, 1)
array([[ 0. ,  0. ,  0. ,  0. ,  0.7,  0.5,  0.6,  0.5,  0.6,  0.9],
       [ 0.1,  0.6,  1. ,  0.2,  0.3,  0.6,  0. ,  0.2,  0.9,  0.9],
       [ 0.2,  0.3,  0.3,  0.5,  0.2,  0.3,  0.4,  0.1,  0.7,  0.5],
       [ 0.9,  0.1,  0.7,  0.8,  0.2,  0.9,  0.9,  0.1,  0.3,  0.3],
       [ 0.8,  0.8,  1. ,  0.9,  0.2,  0.5,  1. ,  0. ,  0. ,  0. ],
       [ 0.2,  0.3,  0.5,  0.4,  0.6,  0.2,  0. ,  0. ,  0. ,  0. ],
       [ 0. ,  0.2,  1. ,  0.2,  0.8,  0. ,  0. ,  0.7,  0.9,  0.6],
       [ 0. ,  0.2,  0.9,  0.9,  0.3,  0.4,  0.6,  0.6,  0.3,  0.6],
       [ 0. ,  0. ,  0. ,  0. ,  0.8,  0.8,  0.1,  0.7,  0.4,  0.4],
       [ 0. ,  0. ,  0. ,  0. ,  0. ,  0.5,  0.1,  0.4,  0.3,  0.9]])

---

Another option, to manage the edges, is to use a convolution function from scipy.ndimage:

import scipy.ndimage
mask = scipy.ndimage.convolve(nb_cameras > 1, np.ones((3,3)),
                              mode='constant') != 0
img_final[mask] = image[mask]
>>> np.round(img_final, 1)
array([[ 0.6,  0.8,  0.7,  0.9,  0.7,  0.5,  0.6,  0.5,  0.6,  0.9],
       [ 0.1,  0.6,  1. ,  0.2,  0.3,  0.6,  0. ,  0.2,  0.9,  0.9],
       [ 0.2,  0.3,  0.3,  0.5,  0.2,  0.3,  0.4,  0.1,  0.7,  0.5],
       [ 0.9,  0.1,  0.7,  0.8,  0.2,  0.9,  0.9,  0.1,  0.3,  0.3],
       [ 0.8,  0.8,  1. ,  0.9,  0.2,  0.5,  1. ,  0. ,  0.3,  0.8],
       [ 0.2,  0.3,  0.5,  0.4,  0.6,  0.2,  0. ,  0. ,  0.7,  0.6],
       [ 0.2,  0.2,  1. ,  0.2,  0.8,  0. ,  0. ,  0.7,  0.9,  0.6],
       [ 0. ,  0.2,  0.9,  0.9,  0.3,  0.4,  0.6,  0.6,  0.3,  0.6],
       [ 0.4,  1. ,  0.8,  0. ,  0.8,  0.8,  0.1,  0.7,  0.4,  0.4],
       [ 0.9,  0.5,  0.8,  0. ,  0. ,  0.5,  0.1,  0.4,  0.3,  0.9]])

Answer 3

This got quite elegant using skimage.morphology's binary_dilation function. It will take a binary array, and kinda expand any pixels that are true into a 3x3 grid of true values (or any other size). This should also handle cases at the edges. Which i think your implementation did not.

Using this mask it's quite easy to calculate imgFinal

from skimage.morphology import binary_dilation, square
mask = binary_dilation(nbCameras > 1, square(maskWidth))
imgFinal = np.where(mask, image, 0)

square(3) is just shorthand for np.ones((3,3))

http://scikit-image.org/docs/dev/api/skimage.morphology.html?highlight=dilation#skimage.morphology.dilation

Example use of dilation for better explenation of what it does:

In [27]: a
Out[27]:
array([[ 1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.]])

In [28]: binary_dilation(a, square(3))
Out[28]:
array([[1, 1, 0, 0, 0],
       [1, 1, 0, 0, 0],
       [0, 0, 1, 1, 1],
       [0, 0, 1, 1, 1],
       [0, 0, 1, 1, 1]], dtype=uint8)