为一组点生成最近邻距离图

Question

为一组点生成最近邻距离图

3

我有一组（x,y）坐标点。我想要生成一个距离每个点的图像，因此我的简单函数如下：

from scipy.spatial.distance import cdist
from numpy import *

def findDistances(imageSize, points):
    image = zeros(imageSize)
    for x in range(0,imageSize[1]):
        for y in range(0,imageSize[0]):
            image[y,x] = np.min(cdist(array([[x,y]]), points))
    return image

这个函数很好，它提供了我想要的结果（见下图）。对于一个大约1MP的图像，处理时间需要大约100秒，对于只需要一次处理的任务来说已经足够好了，但我认为还有改进的空间。我也尝试过：

image[y,x] = min(linalg.norm(array([[x,y]])- points, axis=1))

我认为这两者的运行时间相当 - 这很有道理，它们可能在内部执行类似的操作，尽管我没有检查源代码以确保。

我研究了Scipy的cKDTree：

from scipy.spatial import cKDTree

def findDistancesTree(imageSize, points):
    tree = cKDTree(points)
    image = zeros(imageSize)
    for x in range(0,imageSize[1]):
        for y in range(0,imageSize[0]):
            image[y,x] = tree.query([x,y])[0]
    return image

调用tree.query大约需要50微秒（来自%timeit），但实际上生成一张100万像素的距离图需要70-80秒。比起受到踢腿的打击，20秒的改进是好的，但我不知道是否可以进一步改进。

%timeit np.min(linalg.norm(array([[random.random()*1000,random.random()*1000]]) - points, axis=1))
%timeit np.min(cdist(array([[random.random()*1000,random.random()*1000]]), points))
%timeit tree.query(array([[random.random()*1000,random.random()*1000]]))

10000 loops, best of 3: 82.8 µs per loop
10000 loops, best of 3: 67.9 µs per loop
10000 loops, best of 3: 52.3 µs per loop

就复杂度而言，暴力搜索应该是类似于 O(NM) 的，其中 N 是图像中的像素数量，M 是要检查的点的数量。我原本期望速度会更快，因为对于每个像素来说，搜索时间应该是类似于 O(N log(M)) 的，每次查找需要 log(M) 的时间 - 我有什么遗漏了吗？

- Josh

我认为基本的kd树最近邻搜索通常不是O(log(M))操作，因为你的点不是随机分布的。 - Fabricator

3个回答

0

也许沃罗诺伊图和福特算法可以帮助你找到想要的解决方案。

具体思路如下：

计算 M 个点的 Voronoi 图，时间复杂度为 O(M ln(M))
对于每个像素 N
- 找到相应的区域
- 计算到该区域点的距离

现在的问题是如何在 O(ln(M)) 的时间内找到与像素对应的区域。

希望这可以帮到你！

- Flabetvibes

我想过了。我尝试使用Scipy的voronoi函数和matplotlib的路径。检查每个多边形是否被击中大约需要10µs。问题是有大约1200个多边形，检查每一个都等同于蛮力算法，我猜肯定有更聪明的方法来决定要检查哪些多边形。实际的图表生成非常快，只需要几毫秒。 - Josh

我猜还有更聪明的方法。在CGAL文档中，"点定位"查询有很好的文档记录。它可能是一个很好的灵感来源。 - Flabetvibes

点定位非常难以理解。最好不要去深究。 - Micromega

困难但高效。值得了解。 - Flabetvibes

0

您可以使用分层聚类对点进行聚类，即树状图。对于每棵树，您可以计算 Voronoi 图和点在多边形内的函数。

- Micromega

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- M4rtini · Accepted Answer

这个问题似乎是一个即使使用GPU进行基本暴力实现也能得到良好改进的问题。所以我尝试了一下。改进效果相当不错。

我使用pyopencl进行测试。

import pyopencl as cl 
import numpy as np 

def findDistances_cl(imageSize, points):
    #Boilerplate opencl code
    ctx = cl.create_some_context()
    queue = cl.CommandQueue(ctx)

    f = open('nn.cl', 'r')
    fstr = "".join(f.readlines())
    program = cl.Program(ctx, fstr).build()

    #Creating buffers for the opencl kernel
    mf = cl.mem_flags
    img = np.empty(imageSize, dtype=np.float32)
    x_buf = cl.Buffer(ctx, mf.READ_ONLY | mf.USE_HOST_PTR, hostbuf=points[:,0].astype(np.float32))
    y_buf = cl.Buffer(ctx, mf.READ_ONLY | mf.USE_HOST_PTR, hostbuf=points[:,1].astype(np.float32))
    n_points = cl.Buffer(ctx, mf.READ_ONLY | mf.USE_HOST_PTR, hostbuf=np.array([len(points)],dtype=np.int))
    img_buf = cl.Buffer(ctx, mf.WRITE_ONLY, img.nbytes)

    #Run the kernel
    exec_evt = program.nn(queue, img.shape, None, img_buf, x_buf, y_buf, n_points)
    exec_evt.wait()
    #read back the result
    cl.enqueue_read_buffer(queue, img_buf, img).wait()

    return img

OpenCL内核（nn.cl）

__kernel void nn(__global float *output, __global constant float *x , __global constant float *y, __global constant int *numPoints)
    {
        int row = get_global_id(0);
        int col = get_global_id(1);

        int numRows = get_global_size(0);
        int numCols = get_global_size(1);

        int gid = col+ row*numCols;

        float minDist = numRows * numCols;

        for(int i = 0; i < *numPoints; i++){
          minDist = min(minDist, sqrt((row - y[i])*(row - y[i]) + (col - x[i])*(col - x[i])));
        }
        output[gid] = minDist;
    }

计时结果。

imageSize = [1000, 1000]
points = np.random.random((1000,2))*imageSize[0]

In [4]: %timeit findDistancesTree(imageSize, points)
1 loops, best of 3: 27.1 s per loop

In [7]: %timeit findDistances_cl(imageSize, points)
10 loops, best of 3: 55.3 ms per loop

大约提高了490倍的速度。如果需要更快的速度，还有更先进的算法可以使用GPU进行最近邻居计算。