Remove outliers from a point cloud

Question

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[0]-points[:,0])**2.+(pt[1]-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.

Answer 1

Please add a docstring.

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

Answer 2

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[0], 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