Octant partition equally a cloud of points in 3D

Question

Is there a better way to do this? Not necessarily prettier, although it would be nice.

P = [N,3] ‘Cloud of points’

P -= np.sum(P, axis=0) / P.shape[0]

Map = [i for i in range(P.shape[0])]

p_0 = Map[P[:,0] <= 0] 
p_1 = Map[P[:,0] > 0]

p_0_0 = p_0[P[p_0,1] <= 0]
p_0_1 = p_0[P[p_0,1] > 0]
p_1_0 = p_1[P[p_1,1] <= 0]
p_1_1 = p_1[P[p_1,1] > 0]

p_0_0_0 = p_0_0[P[p_0_0,2] <= 0]
p_0_0_1 = p_0_0[P[p_0_0,2] > 0]
p_0_1_0 = p_0_1[P[p_0_1,2] <= 0]
p_0_1_1 = p_0_1[P[p_0_1,2] > 0]
p_1_0_0 = p_1_0[P[p_1_0,2] <= 0]
p_1_0_1 = p_1_0[P[p_1_0,2] > 0]
p_1_1_0 = p_1_1[P[p_1_1,2] <= 0]
p_1_1_1 = p_1_1[P[p_1_1,2] > 0]

Or in other words, is there a way to compound conditions like,

Oct_0_0_0 = Map[P[:,0] <= 0 and P[:,1] <= 0 and P[:,2] <= 0]

I’m assuming a loop won’t be better than this… not sure.

Thanks in advance.

Answer 1

Instead of repeatedly slicing and keeping lists of the indices, I'd instead recommend creating a single array that maps the index of the point to the octant it belongs to. I'd argue that this is a more natural way of doing it in numpy. So for instance with

octants = (P>0) @ 2**np.arange(P.shape[2])

the n-th entry of octants is the index (here in the range 0,1,2,...,7) of the octant that P belongs to. This works by checking each coordinate whether it is positive or not. This gives three boolean values per point which we can interpret as the binary expansion of that index. (In fact, the line above works for indexing the 2^d-ants in any number of dimensions d.)

To demonstrate this solution, following snippet makes a point cloud and colours them according to their quadrant:

import numpy as np
P = np.random.rand(10000, 3)*2-1
quadrants = (P > 0) @ 2**np.arange(3)
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(*P.T, c=quadrants, cmap='jet')
plt.show()

If you still need to extract an array of the indices of a specific quadrant, say quadrant 0, 1, 1 this corresponds to finding the corresponding decimal number, which is 0*2^0 + 1*2^1 + 1*2^2 = 6, which you can do with e.g.

p_0_1_1 = np.where(quadrants == np.array([0,1,1]) @ 2**np.arange(3))

Answer 2

All in all, I’m going with an even messier version,

p_x = P[:,0] < 0

p_0 = np.nonzero(p_x)[0]
p_1 = np.nonzero(~p_x)[0]

p_0_y = P[p_0,1] < 0
p_1_y = P[p_1,1] < 0

p_0_0 = p_0[p_0_y]
p_0_1 = p_0[~p_0_y]
p_1_0 = p_1[p_1_y]
p_1_1 = p_1[~p_1_y]

p_0_0_z = P[p_0_0,2] < 0
p_0_1_z = P[p_0_1,2] < 0
p_1_0_z = P[p_1_0,2] < 0
p_1_1_z = P[p_1_1,2] < 0

p_0_0_0 = p_0_0[p_0_0_z]
p_0_0_1 = p_0_0[~p_0_0_z]
p_0_1_0 = p_0_1[p_0_1_z]
p_0_1_1 = p_0_1[~p_0_1_z]
p_1_0_0 = p_1_0[p_1_0_z]
p_1_0_1 = p_1_0[~p_1_0_z]
p_1_1_0 = p_1_1[p_1_1_z]
p_1_1_1 = p_1_1[~p_1_1_z]

I have a religious belief (I mean based on thin air) that a comparison is fundamentally cheaper than an arithmetic operation.