Efficiently Compute Voxel Indices for a Point Cloud

user_7805352 0
Viewed 1598

Suppose I have an input point cloud X, represented by an array of dimensions N x 3. Each row in this array corresponds to a single point in XYZ space, between -1 and 1. Now, let k be a parameter which defines the resolution of a voxel grid (e.g. k = 3 corresponds to a voxel grid of dimensions 3 x 3 x 3). I am looking for an efficient way to compute the index for each point's corresponding voxel. This is more or less the current way I'm currently doing it using NumPy (written more expressively for sake of clarity): 

    # Generate some random input point cloud, 100 points
    X = np.random.randn(100, 3)
    # Define the resolution of the grid (e.g. say 3 x 3 x 3)
    k = 3
    # Iterate through points of input point cloud
    for j, point in enumerate(X):
        # Placeholder for voxel membership
        partitions = np.zeros(3,)

        for i, p in enumerate(point):
            for d in range(k):
                # Points are between -1 and 1, so the interval for each dimension is [-1, 1]
                # Length 2, "left"/"lower" end is -1
                if p <= (-1 + (d + 1) * 2 / k):
                    partitions[i] = d

        # Compute the index of the voxel this point belongs to
        # Can think of this as the base 10 representation of the base k number given by (x, y, z) in partitions
        # Voxels are indexed such that (0, 0, 0) --> index 0, (0, 0, 1) --> index 1, (0, 0, 2) --> 
        # index 2, (0, 1, 0) --> index 3, etc.
        p_reversed = np.flip(partitions)
        idx= 0
        for d in range(3):
            idx += (k ** d) * p_reversed[d]

       # Now idx corresponds to the index of the voxel to which point j belongs

This clearly scales poorly with increasing N and increasing k; is there a more efficient implementation?

python-3.x numpy python
2 Answers

You are, in effect, comparing your point values as float against a series of other floats between -1 and 1.

However, what you want to do is compute (once) a function that produces the value. Perform a simple computation instead of iterating.

Ideally, that simple function would be something that numpy could distribute across the columns of your point value.

More ideally, it could be distributed across your entire array, allowing you to use the 2-d array in a single operation or series of operations.

Because you are using a fixed voxel size, and because you are using the same size and same range for all of your dimensions, I think you can do it by simple subtraction and multiplication:

  1. First, subtract the amount required to move the start of your range to zero. In this case, since your range is [-1, 1), you would subtract -1 (or add +1) to the point values to get them to start at 0.

  2. Next, "scale" the point values into the range [0, 1). You can multiply by the inverse of the length of the range (high - low == 1 - -1 == 2), so multiply by 1/2 == 0.5.

  3. At this point, your intermediate values are the fraction of the range where the point occurs. So map them into voxel space by multiplying that fraction by the size of the voxel range (3). And convert the resulting value to an integer either explicitly (via a function) or implicitly (via a datatype). The Zen of Python says that explicit is better than implicit, so that's my preference.

How can you code this?

RANGE_START = -1
RANGE_END = 1
VOXEL_K = 3

# Make up some points in the valid range
CLOUD = np.random.random((100, 3)) * (RANGE_END - RANGE_START) + RANGE_START

# shift points to 0-based:
zero_based_points = CLOUD - RANGE_START

# convert points to [0, 1) fraction of range
fractional_points = zero_based_points / (RANGE_END - RANGE_START)

# project points into voxel space: [0, k)
voxelspace_points = fractional_points * VOXEL_K

# convert voxel space to voxel indices (truncate decimals: 0.1 -> 0)
voxel_indices = voxelspace_points.astype(int)

NOTE: Floating point values can be nan or inf which do not convert well to integers. So you should probably pre-process your points to filter those values somehow (replace them with a sentinel value or eliminate them from the dataset or ...?)

user_40290140

Here's my approach:

N, k = 3, 10
np.random.seed(1)
X = np.random.normal(-1,1, (100,3))

# X:
# [[ 0.62434536 -1.61175641 -1.52817175]
# [-2.07296862 -0.13459237 -3.3015387 ]
# [ 0.74481176 -1.7612069  -0.6809609 ]]

# output:
np.argmax(X[:,:,None] < np.linspace(-1,1,k)[None, None, :], axis=-1)

Output:

array([[8, 0, 0],
       [0, 4, 0],
       [8, 0, 2]])
user_42384080
Related