Efficiently computing standard deviation of all points in a given area for every point in a dense point cloud

Question

I have a structure from motion (SfM) 3D dense point cloud of a landscape that has 200,000,000+ points (and this is one of the smaller point clouds) and am trying to compute the standard deviation of every point given a user-specified search radius so that I can feed these values into a machine learning model alongside other metrics and information. I'm using a cKDTree and the query_ball_point function from scipy to perform the neighborhood search. I think the real slowdown is coming from having to iterate over every point in the point cloud but don't know of a more efficient way to accomplish this.

Here is the code that I currently have to compute the standard deviation of every point in the dense point cloud and return a numpy array of standard deviation values. Is there a more efficient way to iterate over all points in the dense point cloud?

import math
from scipy.spatial import cKDTree

def calc_3d_sd(coords, rad=0.5):
    # build the KDTree
    tree = cKDTree(coords, leafsize=5)
    # intiialize an empty numpy array of the same length as coords
    sd = np.zeros(len(coords))
    # iterate over every point in the dense point cloud
    # (I think this is where the real slowdown is happening)
    for count,elem in enumerate(coords):
        # perform spatial query on point
        result = tree.query_ball_point(elem, r=rad)
        # if at least one other point is returned, then continue
        if len(result) > 0:
            # intialize 'sums' var to track sum
            sums = 0
            # compute standard deviation of X, Y, and Z separately
            # then combine these through the sqrt of the sum of squares to get 3D standard deviation
            for x in np.std(coords[result],axis=0):
                sums += x**2
            sd[count] = math.sqrt(sums)
        # otherwise, no points were found in the search radius, return zero
        else:
            sd[count] = 0
    # return the numpy array of computed standard deviation values all points in the cloud
    return sd

# create 'coords' variable containing all X, Y, and Z coordinates from the dense cloud
coords = np.stack([testcloud.x, testcloud.y, testcloud.z]).transpose()

# iterate over all points in the dense cloud and return the standard deviation of all other points within a 0.5m radius
sd3d = calc_3d_sd(coords, rad=0.5)

In testing a couple of different search radius for the spatial query, it is not surprising to find that the computation time increases substantially with larger radius queries. This doesn't address the primary question though, which is still searching for a more efficient approach to iterating over every point in the dense point cloud (or some equivalent function).

user_1609514 · Answer

Passing all the points as an array to query_ball_point achieves about a 100% speed up I think:

def calc_3d_sd2(coords, rad=0.5):
    # build the KDTree
    tree = cKDTree(coords, leafsize=5)
    # intiialize an empty numpy array of the same length as coords
    sd = np.zeros(len(coords))
    # iterate over every point in the dense point cloud
    # (I think this is where the real slowdown is happening)
    # perform spatial query on all points
    results = tree.query_ball_point(coords, r=rad)
    for count, result in enumerate(results):
        # if at least one other point is returned, then continue
        if len(result) > 0:
            # intialize 'sums' var to track sum
            sums = 0
            # compute standard deviation of X, Y, and Z separately
            # then combine these through the sqrt of the sum of squares to get
            # 3D standard deviation
            for x in np.std(coords[result],axis=0):
                sums += x**2
            sd[count] = math.sqrt(sums)
        # otherwise, no points were found in the search radius, return zero
        else:
            sd[count] = 0
    # return the numpy array of computed standard deviation values all points
    # in the cloud
    return sd

My test script is only with 10,000 points though:

class testcloud:
    np.random.seed(0)
    n_pts = 10000
    data = 10*np.random.rand(3, n_pts)
    x, y, z = data

# create 'coords' variable containing all X, Y, and Z coordinates from the
# dense cloud
coords = np.stack([testcloud.x, testcloud.y, testcloud.z]).transpose()

sd3d = calc_3d_sd(coords, rad=0.5)
# %timeit: 987 ms ± 21.2 ms per loop 

sd3d2 = calc_3d_sd2(coords, rad=0.5)
# %timeit: 414 ms ± 8.49 ms per loop

assert(np.isclose(sd3d2, sd3d).all())

UPDATE:

Managed to get a slight additional speed improvement by vectorizing the std. dev. calculations. Note that the condition should be len(result) > 1 because result always contains the point at the centre of the ball and you only need to compute the std. deviation if there are two points or more:

    sd = np.zeros((len(results), 3))
    for i, result in enumerate(results):
        if len(result) > 1:
            sd[i] = np.std(coords[result], axis=0)
    sd = np.sqrt(np.sum(sd ** 2, axis=1))

Efficiently computing standard deviation of all points in a given area for every point in a dense point cloud

1 Answers