# Remove outliers from a point cloud

This function accepts a cloud of points, and returns those points that are within delta distance of the average (mean) position.

def points_average(points,delta):
""" this function will check, for every point in points
what are the points that are near the point (below a distance delta)
it will then average every such points, creating a new list of points.
This will continue until we arrive at a fixed amount of points """

L = math.inf
n = 0

while L > len(points):

L = len(points)

points2 = []

for i in xrange(len(points)):

pt = points[i]
d = (pt-points[:,0])**2.+(pt-points[:,1])**2.
pts = points[d<delta**2.]

x = np.average(pts[:,0])
y = np.average(pts[:,1])

points2 += [[x,y]]

points2 = np.array(points2)
points = np.unique(points2,axis=0)
print len(points)

return points


The function can be tested with:

import numpy as np

delta = 1
points = np.random.uniform(0,100,(500,2))
new_pts = simplify_by_avg_weighted(points,delta)


Please review the code for any improvements. In particular, I would like to replace the for loop with a numpy vectorized version - any advice towards that would be welcome.

• I've made a substantial edit to your description of the problem solved by this program, to bring it into line with site standards - please check whether I've misunderstood anything, and if necessary edit to correct it. – Toby Speight Dec 21 '17 at 14:42
• Take a look on an KDTree approach. For example stackoverflow.com/questions/35459306/… – max9111 Jan 2 '18 at 14:42
• I guess that is the function you are searching for scipy.spatial.KDTree.query_ball_tree docs.scipy.org/doc/scipy-0.14.0/reference/generated/… – max9111 Jan 2 '18 at 14:46
• @max9111 thank you for the suggestion, I made it work – adrienlucca.wordpress.com Jan 10 '18 at 8:37

Simplify_by_avg_weighted might be named points_near_centroid (or _center if you prefer).

Please initialize with L = math.inf, to clarify your intent.

Rather than a variable length list for [x,y,1], please store it as a 3-tuple (x, y, 1), or more naturally as a pair ((x, y), 1), where 1 is a count rather than a spatial coordinate.

    for i in xrange(len(points)):


This would more naturally be expressed as for pt in points:.

Now we come to the heart of the matter. It's a little hard to evaluate the code in the loop, as it seems to differ from your English language description, and it's hard to declare a "violates spec" bug when there's no docstring offering a spec.

When you assign x, y, they might better have names like x_centroid. And since z is not a spatial coordinate, it would more appropriately be named count. It's fine that you don't take square root, but d should have a name like dsq for distance squared.

To "return those points that are within delta of the mean position", does the interpreter even need to loop? It should suffice to compute centroid of the whole cloud (numpy does the looping), then select on distance (more numpy looping).

This is the result I got thanks to the KDTree approch suggested by max9111. It is indeed much faster as I do not compute the distances but only specify a cutoff distance:

def simplify_by_avg(points, cutoff):

points = np.unique(points, axis = 0)
points = np.hstack((points, np.ones((points.shape, 1), dtype=points.dtype)))

L = np.inf
n = 0
print len(points),"points, n=",n

while L > len(points):
n += 1
L = len(points)

tree = spatial.cKDTree(points[:,0:2])
groups = tree.query_ball_point(points[:,0:2], delta)

new_pts = []

for i in xrange(len(groups)):

ind = groups[i]
pts = points[ind]
x = np.average(pts[:,0])
y = np.average(pts[:,1])
z = np.sum(pts[:,2])

new_pts += [[x,y,z]]

points = np.unique(new_pts,axis=0)

print len(points),"points, n=",n

return points