Улучшение метода для пространственно-известного медианного фильтра для точечных облаков в Python

Question

У меня есть данные облачных облаков от бортового LiDAR. Это шумно, поэтому я хочу запустить медианный фильтр, который собирает точки в пределах N метров каждой точки, находит медианное значение высоты и возвращает среднюю медианную величину в качестве скорректированного значения высоты.

Это примерно аналогично гридингу данных, а также медианности высот в каждом бункере. Или scipy.signal.medfilt.

Но - Я хочу сохранить местоположение (x, y) каждой точки. Также я не уверен, что medfilt сохраняет требуемую пространственную информацию.

У меня есть метод, но он включает несколько циклов. Дорого, когда миллионы пунктов идут в

Обновлено - для каждой итерации первого цикла для оптимальной операции пересечения выбран небольшой патч точек. Первая версия искала все входные точки для пересечения на каждой итерации. Теперь, только небольшой участок, в то время превращается в стройную геометрии и используется для пересечения:

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] > point.x - radius+0.5) &\ 
         (pointcloud[:,0] < point.x + radius+0.5) &\ 
         (pointcloud[:,1] > point.y - radius+0.5) &\ 
         (pointcloud[:,1] < 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 

Это работает:

• глотания список/массив XYZ указывает
• в первый цикл проходит по списку и для каждой точки:
  • выбирает патч облака точек только больше, чем указано окрестности
  • использует стройным, чтобы поместить буфер 3 метра вокруг точки
  • находит пересечение буфера и облако всей точки
  • извлекает множество точек от этой операции в другом цикл
• поиск медианы и добавление его в список новых Z значений
• возвращая список новых Z значений

для 10^4 точек, я получаю результат в 11 секунд. За 10^5 пунктов 3 минуты, и большинство моих наборов данных попадают в 2- 5 * 10^6 баллов. В облаке 2 * 10^6 очков он бежал всю ночь.

Что я хочу - это более быстрый/более эффективный метод!

Я занимаюсь python-pcl, который быстро применяется для фильтрации облаков точек, но я не знаю, как вернуть индексы точек, которые передают/сбой pcl-python-фильтров. Мне нужны эти индексы, потому что каждая точка имеет другие атрибуты, которые должны оставаться прикрепленными к ней.

Если кто-нибудь может предложить более эффективный метод, пожалуйста, сделайте это - я был бы очень признателен за вашу помощь. Если он не может идти быстрее и этот код полезен, не стесняйтесь его использовать.

Спасибо!

Answer 1

После хорошего совета, я попытался это:

#import numpy and scikit-learn's neighbours modulw 
import numpy as np 
from sklearn import neighbors as nb 

#make a little time ticker 
from datetime import datetime 
startTime = datetime.now() 

# generate a KDTree object. This takes ~95% of the 
# processing time 
tree = nb.KDTree(xyzi[:,0:3], leaf_size=60) 

# how long did tree generation take 
print(datetime.now() - startTime) 

#initialise a list 
new_z = [] 

#for each point, collect neighbours within radius r 
nhoods = tree.query_radius(xyzi[:,0:3], r=3) 

# iterate through the list of neighbourhoods, 
# find the median height, and add it to the output list 
for point in nhoods: 
    new_z.append(np.median(xyzi[point,2])) 

# how long did it take? 
print(datetime.now() - startTime) 

Эта версия приняло 33 минут ~ чуть меньше двух миллионов точек. Приемлемо, но все же может быть лучше.

Может ли генерация KDTree ускоряться с использованием метода% jit?

Есть ли лучший способ, чем перебирать все окрестности, чтобы найти окрестности? здесь, NHOOD является массив из-массивов - Я думал, что-то вроде:

median = np.median(nhoods[:][:,2]) 

... но это не так.

Спасибо!