6

I need to check if a point lies inside a bounding cuboid. The number of cuboids is very large (~4M). The code I come up with is:

import numpy as np

# set the numbers of points and cuboids
n_points = 64
n_cuboid = 4000000

# generate the test data
points = np.random.rand(1, 3, n_points)*512
cuboid_min = np.random.rand(n_cuboid, 3, 1)*512
cuboid_max = cuboid_min + np.random.rand(n_cuboid, 3, 1)*8

# main body: check if the points are inside the cuboids
inside_cuboid = np.all((points > cuboid_min) & (points < cuboid_max), axis=1)
indices = np.nonzero(inside_cuboid)

It takes 8 seconds to run np.all and 3 seconds to run np.nonzero on my computer. Any idea to speed up the code?

zell
  • 9,830
  • 10
  • 62
  • 115
f. c.
  • 1,095
  • 11
  • 26

1 Answers1

3

We can reduce memory congestion for all-reduction with slicing along the smallest axis length of 3 to get inside_cuboid -

out = (points[0,0,:] > cuboid_min[:,0]) & (points[0,0,:] < cuboid_max[:,0]) & \
      (points[0,1,:] > cuboid_min[:,1]) & (points[0,1,:] < cuboid_max[:,1]) & \
      (points[0,2,:] > cuboid_min[:,2]) & (points[0,2,:] < cuboid_max[:,2])

Timings -

In [43]: %timeit np.all((points > cuboid_min) & (points < cuboid_max), axis=1)
2.49 s ± 20 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [51]: %%timeit
    ...: out = (points[0,0,:] > cuboid_min[:,0]) & (points[0,0,:] < cuboid_max[:,0]) & \
    ...:       (points[0,1,:] > cuboid_min[:,1]) & (points[0,1,:] < cuboid_max[:,1]) & \
    ...:       (points[0,2,:] > cuboid_min[:,2]) & (points[0,2,:] < cuboid_max[:,2])
1.95 s ± 10.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Divakar
  • 218,885
  • 19
  • 262
  • 358
  • Thanks. It helps a lot! It saves 60% of the time on my computer. But I need to make it at least 10 times faster. Any other ideas? – f. c. Sep 19 '20 at 12:32
  • @f.c. Given your specific requirement of 10x speedup, a bounty might be fair from you, so that we all get the added motivation and also more attention. – Divakar Sep 22 '20 at 18:29
  • I have turned to another solution to make it faster. Please check out my new post https://stackoverflow.com/questions/64043572/how-to-speed-up-an-n-dimensional-interval-tree-in-python. Will add a bounty to that question if needed. – f. c. Sep 24 '20 at 10:55