2

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]
Romanzo Criminale
  • 1,289
  • 2
  • 14
  • 21

3 Answers3

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)
M4rtini
  • 13,186
  • 4
  • 35
  • 42
  • Thanks for your answer. But I don't think in my case I want to use a dilation as I want keep the same pixel values as my image. And doing a dilation would increase the brightness of my image – Romanzo Criminale Jan 14 '14 at 01:32
  • @RomanzoCriminale The dilation is not changing the image at all, it is just used to create a mask from nbCameras with all pixels, and neighbors, where nbCameras > 1. This mask is used to get the values from image. – M4rtini Jan 14 '14 at 05:06
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: Result

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

Thorsten Kranz
  • 12,492
  • 2
  • 39
  • 56
1

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]])
Jaime
  • 65,696
  • 17
  • 124
  • 159
  • Thanks Jaime. I didn't know about strides and it seem pretty interested to replace for loop. But I don't really understand why your are using shape=image.shape + (mask_width, mask_width) because this gives me a shape of (10,10,3,3). Also why are you using image.strides*2? shouldn't it be image.strides to go through every pixels? – Romanzo Criminale Jan 14 '14 at 01:37
  • That's the nice thing: it creates a view of your data with shape (10,10,3,3) where the (3,3) array in a[i,j] is a[i:i+3,j:j+3] – Jaime Jan 14 '14 at 03:03
  • Ok I get it. But what about this line : copy_mask = nb_cameras[mask_radius:-mask_radius, mask_radius:-mask_radius] > 1 Here in what I understand copy_mask is fixed from nb_cameras but shouldn't it move along nb_cameras? – Romanzo Criminale Jan 14 '14 at 03:27