title: Improving a method for a spatially aware median filter for point clouds in Python
body: <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>