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.