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

Question

我有一组 (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

这个函数运行得不错，能得到我想要的结果（见下面的图片）。不过，生成一张大约100万像素的图需要大约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 的结果），但实际上生成一张1MP的距离图需要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)，我是不是漏掉了什么？

Answer 1

你可以用一种叫做层次聚类的方法来把这些点分组，这种方法会生成一个树状图。对于每一棵树，你还可以计算出一个叫做Voronoi图的东西，以及一个判断点是否在多边形内的功能。

Answer 2

也许 Voronoi图 和 福琼算法 能帮助你找到想要的解决方案。

大概的思路是：

• 先计算出你那些M个点的Voronoi图，时间复杂度是O(M ln(M))
• 然后对于每一个N个像素：
  • 找出对应的区域
  • 计算这个区域到点的距离

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

希望这对你有帮助！

Answer 3

这个问题听起来即使是用最简单的暴力破解方法配合显卡，也能有不错的提升。所以我试了一下，结果提升确实很明显。

我使用了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倍。如果你需要更快的速度，还有一些更高级的算法可以用显卡来做最近邻搜索。