2

我有一个点云 C，其中每个点都有一个关联的值。假设这些点位于二维空间中，因此每个点都可以用三元组 (x, y, v) 表示。

我想找到局部最大值的点子集。也就是说，对于某个半径 R，我想找到 C 中点 S 的子集，使得对于 S 中的任何点 Pi（值为 vi），在 Pi 的 R 距离内没有 C 中的点 Pj，其值 vj 为大于六。

我知道如何在 O(N^2) 时间内做到这一点，但这似乎很浪费。有没有一种有效的方法来做到这一点？

旁注：

  • 这个问题的根源是我试图在稀疏矩阵中找到局部最大值，所以在我的例子中，x，y 是有序整数 indeces - 如果这简化了问题，请告诉我！
  • 如果解决方案仅适用于曼哈顿距离或其他什么，我将非常高兴。
  • 我在 python 中，所以如果有某种很好的矢量化 numpy 方法可以做到这一点，那就太好了。
4

3 回答 3

5

遵循 Yves 的建议，这是一个使用 scipy 的KDTree的答案：

from scipy.spatial.kdtree import KDTree
import numpy as np

def locally_extreme_points(coords, data, neighbourhood, lookfor = 'max', p_norm = 2.):
    '''
    Find local maxima of points in a pointcloud.  Ties result in both points passing through the filter.

    Not to be used for high-dimensional data.  It will be slow.

    coords: A shape (n_points, n_dims) array of point locations
    data: A shape (n_points, ) vector of point values
    neighbourhood: The (scalar) size of the neighbourhood in which to search.
    lookfor: Either 'max', or 'min', depending on whether you want local maxima or minima
    p_norm: The p-norm to use for measuring distance (e.g. 1=Manhattan, 2=Euclidian)

    returns
        filtered_coords: The coordinates of locally extreme points
        filtered_data: The values of these points
    '''
    assert coords.shape[0] == data.shape[0], 'You must have one coordinate per data point'
    extreme_fcn = {'min': np.min, 'max': np.max}[lookfor]
    kdtree = KDTree(coords)
    neighbours = kdtree.query_ball_tree(kdtree, r=neighbourhood, p = p_norm)
    i_am_extreme = [data[i]==extreme_fcn(data[n]) for i, n in enumerate(neighbours)]
    extrema, = np.nonzero(i_am_extreme)  # This line just saves time on indexing
    return coords[extrema], data[extrema]
于 2014-11-24T23:43:11.320 回答
3

使用 2D-tree（kD-tree的 2D 实例）。经过 N.Log(N) 时间预处理后，它将允许您在大约 Log(N) + K 时间（平均找到 K 个邻居）内围绕所有点执行固定半径近邻搜索，总共 N。 Log(N)+ KN 它将完美地适应曼哈顿距离。

于 2014-11-22T20:18:04.947 回答
0

我找到了这个解决方案，但它可能是 O(N^2)：

import numpy as np

# generate test data
n = 10
foo = np.random.rand(n,n)

# fixed test data for visual back-checking
# foo = np.array([[ 0.12439309,  0.88878825,  0.21675684,  0.21422532,  0.7016789 ],
#                 [ 0.14486462,  0.40642871,  0.4898418 ,  0.41611303,  0.12764404],
#                 [ 0.41853585,  0.22216484,  0.36113181,  0.5708699 ,  0.3874901 ],
#                 [ 0.24314391,  0.22488507,  0.22054467,  0.25387521,  0.46272496],
#                 [ 0.99097341,  0.76083447,  0.37941783,  0.932519  ,  0.9668254 ]])

# list to collect local maxima
local_maxima = []

# distance in x / y to define region of interest around current center coordinate
# roi = 1 corresponds to a region of interest of 3x3 (except at borders)
roi = 1

# give pseudo-coordinates
x,y = np.meshgrid(range(foo.shape[0]), range(foo.shape[1]))

for i in range(foo.shape[0]):
    for j in range(foo.shape[1]):
        x0 = x[i,j]
        y0 = y[i,j]
        z0 = foo[i,j]
        # index calculation to avoid out-of-bounds error when taking sub-matrix
        mask_x = abs(x - x0) <= roi
        mask_y = abs(y - y0) <= roi
        mask = mask_x & mask_y
        if np.max(foo[mask]) == z0:
            local_maxima.append((i, j))

print local_maxima

这一切都是关于在矩阵上定义滑动窗口/过滤器。我想到的所有其他解决方案都指向绝对最大值（例如直方图）......

但是我希望我的 ansatz 在某种程度上有用......

编辑：这里的另一个解决方案应该比第一个更快，但仍然是 O(N^2)，它不依赖于直线网格数据：

import numpy as np

# generate test data
# points = np.random.rand(10,3)

points = np.array([[ 0.08198248,  0.25999721,  0.07041999],
                   [ 0.19091977,  0.05404123,  0.25826508],
                   [ 0.8842875 ,  0.90132467,  0.50512316],
                   [ 0.33320528,  0.74069399,  0.36643752],
                   [ 0.27789568,  0.14381512,  0.13405309],
                   [ 0.73586202,  0.4406952 ,  0.52345838],
                   [ 0.76639731,  0.70796547,  0.70692905],
                   [ 0.09164532,  0.53234394,  0.88298593],
                   [ 0.96164975,  0.60700481,  0.22605181],
                   [ 0.53892635,  0.95173308,  0.22371167]])

# list to collect local maxima
local_maxima = []

# distance in x / y to define region of interest around current center coordinate
radius = 0.25

for i in range(points.shape[0]):
        # radial mask with radius radius, could be beautified via numpy.linalg
        mask = np.sqrt((points[:,0] - points[i,0])**2 + (points[:,1] - points[i,1])**2) <= radius
        # if current z value equals z_max in current region of interest, append to result list
        if points[i,2] == np.max(points[mask], axis = 0)[2]:
            local_maxima.append(tuple(points[i]))

结果：

local_maxima = [
 (0.19091976999999999, 0.054041230000000003, 0.25826507999999998), 
 (0.33320527999999999, 0.74069399000000002, 0.36643752000000002), 
 (0.73586202000000001, 0.44069520000000001, 0.52345838), 
 (0.76639731, 0.70796546999999999, 0.70692904999999995), 
 (0.091645320000000002, 0.53234393999999996, 0.88298593000000003), 
 (0.53892635, 0.95173308000000001, 0.22371167)
]
于 2014-11-20T22:47:41.893 回答