-1

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.

Ricardo Lordelo
  • 1
  • 1

2 Answers2

0

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

points coloured by quadrant

flawr
  • 10,814
  • 3
  • 41
  • 71
-1

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.

Ricardo Lordelo
  • 1
  • 1
  • 1
    If you believe one version is faster than another, why don't you just time it? Note that you use the same boolean indexing array two times, (once bare, once negated) so you immediately get some overhead, furthermore you have to think abou the cost of creating so many new arrays. But finally, seeing variable names like these should immediately make you rethink your code, as there *must* be better structures that you can use. – flawr Sep 26 '22 at 22:01