trjhtr

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. 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.