I have the set of all 8 possible signs that non-zero 3D (euclidean). I want to efficiently categorise 3D points pts into octants based on their signs, as part of an ongoing project on octahedral spaces, as per this and many other questions.

The entire project is written in numpy, and it does not make sense to use another tool for this part of the project. The purpose of the partitioning is to handle pipelines when dealing with octant specific transforms, as well as projections.

I want to use those to partition/bin a numpy array points of x,y,z (euclidean) coordinates such that subsequent array belongs to the octant governed by the rule that the first sign will filter x, second y, third z, and a -1 means <0, while a +1 means >= 0.

I can write it all out...

import numpy as np

def octate(pts):
    results = {}
    results[(-1, -1, -1)] = pts[np.where((pts[:, 0] < 0) & (pts[:, 1] < 0) & (pts[:, 2] < 0))]
    results[(-1, -1, +1)] = pts[np.where((pts[:, 0] < 0) & (pts[:, 1] < 0) & (pts[:, 2] >= 0))]
    results[(-1, +1, -1)] = pts[np.where((pts[:, 0] < 0) & (pts[:, 1] >= 0) & (pts[:, 2] < 0))]
    results[(-1, +1, +1)] = pts[np.where((pts[:, 0] < 0) & (pts[:, 1] >= 0) & (pts[:, 2] >= 0))]
    results[(+1, -1, -1)] = pts[np.where((pts[:, 0] >= 0) & (pts[:, 1] < 0) & (pts[:, 2] < 0))]
    results[(+1, -1, +1)] = pts[np.where((pts[:, 0] >= 0) & (pts[:, 1] < 0) & (pts[:, 2] >= 0))]
    results[(+1, +1, -1)] = pts[np.where((pts[:, 0] >= 0) & (pts[:, 1] >= 0) & (pts[:, 2] < 0))]
    results[(+1, +1, +1)] = pts[np.where((pts[:, 0] >= 0) & (pts[:, 1] >= 0) & (pts[:, 2] >= 0))]
    return results


if __name__ == '__main__':
    rxy = np.random.uniform(-4, 4, [1000, 3])
    roc = octate(rxy)

But it seems pretty clunky and redundant.

Likewise, the above uses 8 operations, but it's got to be possible to subdivide based on the <0 / >=0 dichotomy across 3 axes? I would have thought that 3 operations should be feasible.

A manual method involves code redundancy, adds unnecessary fragility to the code, and it lacks scalability to higher dimensions. I am looking for something that is concise, efficient, and readable. I have had a look at binning, partition, and arg_partition but they don't seem to be a good fit, though I might be mistaken.

I have the set of all 8 possible signs that non-zero 3D (euclidean). I want to efficiently categorise 3D points pts into octants based on their signs, as part of an ongoing project on octahedral spaces, as per this and many other questions.

The entire project is written in numpy, and it does not make sense to use another tool for this part of the project. The purpose of the partitioning is to handle pipelines when dealing with octant specific transforms, as well as projections.

I want to use those to partition/bin a numpy array points of x,y,z (euclidean) coordinates such that subsequent array belongs to the octant governed by the rule that the first sign will filter x, second y, third z, and a -1 means <0, while a +1 means >= 0.

I can write it all out...

import numpy as np

def octate(pts):
    results = {}
    results[(-1, -1, -1)] = pts[np.where((pts[:, 0] < 0) & (pts[:, 1] < 0) & (pts[:, 2] < 0))]
    results[(-1, -1, +1)] = pts[np.where((pts[:, 0] < 0) & (pts[:, 1] < 0) & (pts[:, 2] >= 0))]
    results[(-1, +1, -1)] = pts[np.where((pts[:, 0] < 0) & (pts[:, 1] >= 0) & (pts[:, 2] < 0))]
    results[(-1, +1, +1)] = pts[np.where((pts[:, 0] < 0) & (pts[:, 1] >= 0) & (pts[:, 2] >= 0))]
    results[(+1, -1, -1)] = pts[np.where((pts[:, 0] >= 0) & (pts[:, 1] < 0) & (pts[:, 2] < 0))]
    results[(+1, -1, +1)] = pts[np.where((pts[:, 0] >= 0) & (pts[:, 1] < 0) & (pts[:, 2] >= 0))]
    results[(+1, +1, -1)] = pts[np.where((pts[:, 0] >= 0) & (pts[:, 1] >= 0) & (pts[:, 2] < 0))]
    results[(+1, +1, +1)] = pts[np.where((pts[:, 0] >= 0) & (pts[:, 1] >= 0) & (pts[:, 2] >= 0))]
    return results


if __name__ == '__main__':
    rxy = np.random.uniform(-4, 4, [1000, 3])
    roc = octate(rxy)

But it seems pretty clunky and redundant.

Likewise, the above uses 8 operations, but it's got to be possible to subdivide based on the <0 / >=0 dichotomy across 3 axes? I would have thought that 3 operations should be feasible.

A manual method involves code redundancy, adds unnecessary fragility to the code, and it lacks scalability to higher dimensions. I am looking for something that is concise, efficient, and readable. I have had a look at binning, partition, and arg_partition but they don't seem to be a good fit, though I might be mistaken.

• python
• numpy
• binning

Share Improve this question Follow edited Mar 17 at 19:09 Konchog asked Mar 17 at 18:54 KonchogKonchog 2,22211 gold badge2222 silver badges2828 bronze badges

Add a comment  | 

3 Answers 3

Sorted by: Reset to default 3

You could get an index for your octant from binary 0-7:

i = int(p[0] > 0) + 2 * int( p[1] > 0 ) + 4 * int( p[2] > 0 )

For example:

import numpy as np

def octate(points):
    results = [[] for _ in range(8)]
    for p in points:
        i = int(p[0] > 0) + 2 * int( p[1] > 0 ) + 4 * int( p[2] > 0 )
        results[i].append( p )
    return results

rxy = np.random.uniform(-4, 4, [20, 3])
roc = octate(rxy)
for quadrant in range( 8 ):
    print( quadrant, ': ', end='' )
    for p in roc[quadrant]: print( (3*'{:4.1f} ').format( p[0],p[1],p[2] ), ' | ', end='' )
    print()

Output:

0 : -2.2 -0.5 -1.8   | -2.5 -2.7 -2.9   | -0.5 -0.2 -1.4   | 
1 :  3.8 -2.0 -1.4   |  2.0 -2.6 -0.4   |  3.6 -0.8 -2.3   |  3.4 -2.7 -0.1   |  2.8 -1.4 -3.5   |  2.5 -2.1 -1.2   | 
2 : -0.2  2.4 -1.4   | -0.2  0.4 -2.2   | 
3 :  2.8  3.0 -2.3   |  3.5  0.7 -0.4   |  0.3  0.7 -3.9   | 
4 : -1.3 -1.3  1.9   | -0.7 -1.0  2.3   | 
5 :  2.1 -0.9  1.3   | 
6 : -2.2  1.3  1.9   | -2.2  3.6  4.0   | 
7 :  2.1  1.3  3.6   | 

This does more splits, but can handle an arbitrary number of dimensions:

def octate(pts):
    results = {tuple():pts}
    for d in range(0,pts.shape[1]):
        old_results = results
        results = {}
        for k,v in old_results.items():
            results[tuple(list(k)+[-1])] = v[np.where((v[:,d] <  0))]
            results[tuple(list(k)+[+1])] = v[np.where((v[:,d] >= 0))]
    return results

Note that it *does* more splits, but they are over smaller lists with each successive value of `d`, so effectively each list gets scanned twice for each dimension.

You can speed up the code from the previous answer lastchance a little. The variable choices can simply contain the names of the octants ("NEA",...).

import numpy as np
points=np.random.randn(6, 3)
v=np.array([1,2,4])
choices=np.array([[-1,-1,-1],[1,-1,-1],[-1,1,-1],[1,1,-1],[-1,-1,1],[1,-1,1],[-1,1,1],[1,1,1]])
result=choices[(points=>0)@v]