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: enter image description here

AFTER: enter image description here

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.imshow(I[i].reshape(q,q), cmap="viridis", interpolation="none")

N = 10000
q = 28
P = q*q
I = np.ones((N,P))

ps = 5
IDX = np.random.randint(0,P,(N,1))
x0, y0 = np.unravel_index(IDX,(q,q))

t0 = time.time()

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)


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

  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.

  • \$\begingroup\$ Awesome! Now I try to understand everything. \$\endgroup\$
    – Gilfoyle
    Jun 7, 2018 at 13:57
  • \$\begingroup\$ 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? \$\endgroup\$
    – Gilfoyle
    Jun 7, 2018 at 21:44
  • \$\begingroup\$ Use np.maximum(x0 - 1, 0) instead of x0 - 1. \$\endgroup\$ Jun 7, 2018 at 22:04
  • \$\begingroup\$ 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? \$\endgroup\$
    – Gilfoyle
    Jun 13, 2018 at 6:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.