Improving a method for a spatially aware median filter for point clouds in Python
<p>I have point cloud data from airborne LiDAR. It is noisy, so I want to run a median filter which collects points within N metres of each point, finds the median elevation value, and returns the neighbourhood median as an adjusted elevation value.</p>

<p>This is roughly analogous to gridding the data, and taking the median of elevations within each bin. Or scipy.signal.medfilt.</p>

<p>But - I want to preserve the location (x,y) of each point. Also I'm not sure that medfilt preserves the spatial information required.</p>

<p>I have a method, but it involves multiple for loops. Expensive when millions of points go in</p>

<p><strong>Updated</strong> - for each iteration of the first loop, a small patch of points is selected for the shapely intersection operation. The first version searched all input points for an intersection at every iteration. Now, only a small patch at a time is converted to a shapely geometry and used for the intersection:</p>

<pre><code>import numpy as np
from shapely import geometry

def spatialmedian(pointcloud,radius):
"""
Using shapely geometires, replace every point in a cloud with the
median value of points within 'radius' units of the point
'pointcloud' must have no more than 3 dimensions (x,y,z)
"""
    new_z = []
    i = 0

    for point in pointcloud:
        #pick a point and make it a shapely Point
        point = geometry.Point(pointcloud[i,:])

        #select a patch around the point and make it a shapely
        # MultiPoint
        patch = geometry.MultiPoint(list(pointcloud[\
                       (pointcloud[:,0] &gt; point.x - radius+0.5) &amp;\
                       (pointcloud[:,0] &lt; point.x + radius+0.5) &amp;\
                       (pointcloud[:,1] &gt; point.y - radius+0.5) &amp;\
                       (pointcloud[:,1] &lt; point.y + radius+0.5)\
                       ]))      

        #buffer the Point by radius
        pbuff = point.buffer(radius)

        #use the intersection method to find points in our
        # patch that lie inside the Point buffer
        isect = pbuff.intersection(patch)
        #print(isect.geom_type)

        #initialise another list
        plist = []

        #for every intersection set,
        # unpack it into a list and collect the median
        # Z value.
        if isect.geom_type == 'MultiPoint':
            #print('point has neightbours')
            for p in isect:
                plist.append(p.z)

            new_z.append(np.median(plist)) 
        else:
            # if the intersection set isn't MultiPoint,
            # it is an isolated point, whose median Z value 
            # is it's own.
            #print('isolated point')

            #append it to the big list
            new_z.append(isect.z) 
    #iterate i  
    i += 1        

    #print(i)    
    #return a list of new median filtered Z coordinates
    return new_z
</code></pre>

<p>This works by:</p>

<ul>
<li>ingesting a list/array of XYZ points</li>
<li>the first for loop goes through the list and for every point:

<ul>
<li>picks out a patch of the point cloud just bigger than the neighbourhood specified</li>
<li>uses shapely to place a 3 metre buffer around the point</li>
<li>finds the intersection of the buffer and the whole point cloud</li>
<li>extracts the set of points from that operation in another for loop</li>
</ul></li>
<li>finding the median and appending it to a list of new Z values</li>
<li>returning the list of new Z values</li>
</ul>

<p>For 10^4 points, I get a result in 11 seconds. For 10^5 points 3 minutes, and most of my datasets run into 2- 5 * 10^6 points. On a 2 * 10^6 point cloud it's been running overnight.</p>

<p>What I want is a faster/more efficient method!</p>

<p>I've been tinkering with python-pcl, which is fast for filtering point clouds, but I don't know how to return indices of points which pass/fail pcl-python filters. I need those indices because each point has other attributes which must remain attached to it.</p>

<p>If anyone can suggest a more efficient method, please do so - I would highly appreciate your help. If it can't go faster and this code is helpful, feel free to use it.</p>

<p>Thanks!</p>