Remove pixel patch in image which is stored as array

Question

I have an array I which stores N images of size P (number of pixels). Every image is of size P = q*q.

Now I want to delete patches of size ps around a selected index IDX (set all values to zero).

My approach was to reshape every single image using reshape(q,q) and delete the pixels around IDX. I also have to check if the index is not outside the image.

Here is an example:

BEFORE:

AFTER:

My code is a real bottleneck and I would like to know if there is a way to improve the performance of my approach.

import numpy as np import matplotlib.pyplot as plt import time def myplot(I): imgs = 5 for i in range(imgs**2): plt.subplot(imgs,imgs,(i+1)) plt.imshow(I[i].reshape(q,q), cmap="viridis", interpolation="none") plt.axis("off") plt.show() N = 10000 q = 28 P = q*q I = np.ones((N,P)) myplot(I) ps = 5 IDX = np.random.randint(0,P,(N,1)) x0, y0 = np.unravel_index(IDX,(q,q)) t0 = time.time() # HOW TO IMPROVE THIS PART ? # for i in range(N): img = I[i].reshape(q,q) for x in range(ps): for y in range(ps): if (x0[i]+x < q) and (y0[i]+y < q): img[x0[i]+x,y0[i]+y] = 0.0 I[i] = img.reshape(1,q*q) print(time.time()-t0) myplot(I) 

I call this code (without the plotting procedure) about one million times from another code. Every call takes about 1 second on my system. This makes the code so far quite useless.

Any advice?

Answer 1

1. On my computer it takes 1.745 seconds to run the code in the post.

2. There's no need for the array of random indexes to be two-dimensional:

   IDX = np.random.randint(0,P,(N,1)) 

   In fact this is harmful for performance, because it means that x0[i] is an array of length 1 (not a scalar) and so img[x0[i]+x,y0[i]+y] requires "fancy indexing" which is slower than normal indexing.

   It would be simpler to make the array of indexes one-dimensional:

   IDX = np.random.randint(P, size=N) 

   This reduces the runtime to about 0.459 seconds (26.3% of the original).

3. There is no need to reassign I[i] at the end of the loop. When you call the reshape method on a NumPy array, what you get is a view onto the original array (not a copy) if possible. (And it is possible in this case.) So updating the view also updates the original.

   This reduces the runtime to about 0.449 seconds (25.8%).

4. Instead of looping over range(N) and then looking up I[i] and x0[i] and y0[i], use zip to loop over all the arrays simultaneously:

   for img, xx, yy in zip(I, x0, y0): img = img.reshape(q,q) for x in range(ps): for y in range(ps): if xx + x < q and yy + y < q: img[xy + x, yy + y] = 0.0 

   This reduces the runtime to about 0.358 seconds (20.5%).

5. Instead of looping over all the pixels in the patch and updating each pixel individually, use slices to update the whole region in one step:

   for image, x, y in zip(I, x0, y0): image.reshape(q, q)[x:x + ps, y:y + ps] = 0.0 

   This works because NumPy (and Python generally) ensures that the bounds of a slice do not go beyond the end of the array. See the slicing documentation:

   The slice of \$s\$ from \$i\$ to \$j\$ is defined as the sequence of items with index \$k\$ such that \$i \le k < j\$. If \$i\$ or \$j\$ is greater than len(s), use len(s).

   This reduces the runtime to about 0.025 seconds (1.4%).

6. We can vectorize the additions x + ps and y + ps:

   for image, x, y, x1, y1 in zip(I, x0, y0, x0 + ps, y0 + ps): image.reshape(q, q)[x:x1, y:y1] = 0.0 

   This reduces the runtime to about 0.021 seconds (1.2%).

7. We could avoid the reshape inside the loop by doing a single reshape of the whole I array:

   images = I.reshape(N, q, q) 

   and then:

   for image, x, y, x1, y1 in zip(images, x0, y0, x0 + ps, y0 + ps): image[x:x1, y:y1] = 0.0 

   This reduces the runtime to about 0.018 seconds (1.0%).

8. We can halve the number of indexing operations by indexing the images array just once on each loop iteration:

   for i, x, y, x1, y1 in zip(range(N), x0, y0, x0 + ps, y0 + ps): images[i, x:x1, y:y1] = 0.0 

   This reduces the runtime to about 0.011 seconds (0.6%).

That's about 150 times speedup overall, so calling this a million times will still take about 3 hours on my computer. There may be other improvements to be had if only we could see more of your code, but you'll need to make a new post for that.