# Remove pixel patch in image which is stored as array

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?

## 1 Answer

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.

• Awesome! Now I try to understand everything. Jun 7, 2018 at 13:57
• I am still amazed by your answer. Did not think that it would go so fast using Python. There is one thing I noted using your approach. Given a odd patch size, for example ps=3 and coordinates x0 and y0 your approach deletes 9 pixels from top left to bottom right (that is the square from [x0,y0] to [x0+ps,y0+ps]). But how can I delete the surrounding pixels instead, that is [x0-1, y0-1] to [x0+1, y0+1]? I tried just doing x0=x0-1 and y0=y0-1which results quite often in images where no pixels got deleted. What would be the right approach to do that? Jun 7, 2018 at 21:44
• Use np.maximum(x0 - 1, 0) instead of x0 - 1. Jun 7, 2018 at 22:04
• I wanted to replace the patches with random numbers instead of zeros. So I tried images[i, x:x1, y:y1] = np.random.rand(x1-x,y1-y) which does not work. Any suggestions? Jun 13, 2018 at 6:57