计算不规则点密度的有效方法

Question

我正在尝试生成地图覆盖图，这将有助于识别热点，即地图上具有高密度数据点的区域。我尝试过的方法都不够快，不能满足我的需要。 注意：我忘了提到算法应该在低和高缩放场景（或低和高数据点密度）下都能很好地工作。

我查看了numpy、pyplot和scipy库，找到的最接近的是numpy.histogram2d，如下图所示，histogram2d的输出相当粗糙。（每张图片都包括覆盖热图的点，以便更好地理解）

我的第二个尝试是遍历所有数据点，然后计算热点值作为距离的函数。这产生了一个更好看的图像，但它太慢，无法在我的应用程序中使用。因为它是O（n），所以它可以在100个点上工作，但是当我使用我的30000个点的实际数据集时，它会崩溃。

我最后的尝试是将数据存储在KDTree中，并使用最近的5个点来计算热点值。这个算法是O（1），对于大数据集来说要快得多。它仍然不够快，生成一个256x256位图大约需要20秒，我希望在大约1秒的时间内完成。

编辑

6502提供的boxsum平滑解决方案在所有缩放级别都能很好地工作，并且比我原来的方法快得多。

Luke和Neil G提出的高斯滤波解是最快的。

您可以看到下面的四种方法，总共使用1000个数据点，在3倍缩放时，大约有60个点可见。

生成我最初3次尝试的完整代码，6502提供的boxsum平滑解决方案和Luke建议的高斯滤波器（改进以更好地处理边缘并允许放大）如下：

import matplotlib
import numpy as np
from matplotlib.mlab import griddata
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import math
from scipy.spatial import KDTree
import time
import scipy.ndimage as ndi


def grid_density_kdtree(xl, yl, xi, yi, dfactor):
    zz = np.empty([len(xi),len(yi)], dtype=np.uint8)
    zipped = zip(xl, yl)
    kdtree = KDTree(zipped)
    for xci in range(0, len(xi)):
        xc = xi[xci]
        for yci in range(0, len(yi)):
            yc = yi[yci]
            density = 0.
            retvalset = kdtree.query((xc,yc), k=5)
            for dist in retvalset[0]:
                density = density + math.exp(-dfactor * pow(dist, 2)) / 5
            zz[yci][xci] = min(density, 1.0) * 255
    return zz

def grid_density(xl, yl, xi, yi):
    ximin, ximax = min(xi), max(xi)
    yimin, yimax = min(yi), max(yi)
    xxi,yyi = np.meshgrid(xi,yi)
    #zz = np.empty_like(xxi)
    zz = np.empty([len(xi),len(yi)])
    for xci in range(0, len(xi)):
        xc = xi[xci]
        for yci in range(0, len(yi)):
            yc = yi[yci]
            density = 0.
            for i in range(0,len(xl)):
                xd = math.fabs(xl[i] - xc)
                yd = math.fabs(yl[i] - yc)
                if xd < 1 and yd < 1:
                    dist = math.sqrt(math.pow(xd, 2) + math.pow(yd, 2))
                    density = density + math.exp(-5.0 * pow(dist, 2))
            zz[yci][xci] = density
    return zz

def boxsum(img, w, h, r):
    st = [0] * (w+1) * (h+1)
    for x in xrange(w):
        st[x+1] = st[x] + img[x]
    for y in xrange(h):
        st[(y+1)*(w+1)] = st[y*(w+1)] + img[y*w]
        for x in xrange(w):
            st[(y+1)*(w+1)+(x+1)] = st[(y+1)*(w+1)+x] + st[y*(w+1)+(x+1)] - st[y*(w+1)+x] + img[y*w+x]
    for y in xrange(h):
        y0 = max(0, y - r)
        y1 = min(h, y + r + 1)
        for x in xrange(w):
            x0 = max(0, x - r)
            x1 = min(w, x + r + 1)
            img[y*w+x] = st[y0*(w+1)+x0] + st[y1*(w+1)+x1] - st[y1*(w+1)+x0] - st[y0*(w+1)+x1]

def grid_density_boxsum(x0, y0, x1, y1, w, h, data):
    kx = (w - 1) / (x1 - x0)
    ky = (h - 1) / (y1 - y0)
    r = 15
    border = r * 2
    imgw = (w + 2 * border)
    imgh = (h + 2 * border)
    img = [0] * (imgw * imgh)
    for x, y in data:
        ix = int((x - x0) * kx) + border
        iy = int((y - y0) * ky) + border
        if 0 <= ix < imgw and 0 <= iy < imgh:
            img[iy * imgw + ix] += 1
    for p in xrange(4):
        boxsum(img, imgw, imgh, r)
    a = np.array(img).reshape(imgh,imgw)
    b = a[border:(border+h),border:(border+w)]
    return b

def grid_density_gaussian_filter(x0, y0, x1, y1, w, h, data):
    kx = (w - 1) / (x1 - x0)
    ky = (h - 1) / (y1 - y0)
    r = 20
    border = r
    imgw = (w + 2 * border)
    imgh = (h + 2 * border)
    img = np.zeros((imgh,imgw))
    for x, y in data:
        ix = int((x - x0) * kx) + border
        iy = int((y - y0) * ky) + border
        if 0 <= ix < imgw and 0 <= iy < imgh:
            img[iy][ix] += 1
    return ndi.gaussian_filter(img, (r,r))  ## gaussian convolution

def generate_graph():    
    n = 1000
    # data points range
    data_ymin = -2.
    data_ymax = 2.
    data_xmin = -2.
    data_xmax = 2.
    # view area range
    view_ymin = -.5
    view_ymax = .5
    view_xmin = -.5
    view_xmax = .5
    # generate data
    xl = np.random.uniform(data_xmin, data_xmax, n)    
    yl = np.random.uniform(data_ymin, data_ymax, n)
    zl = np.random.uniform(0, 1, n)

    # get visible data points
    xlvis = []
    ylvis = []
    for i in range(0,len(xl)):
        if view_xmin < xl[i] < view_xmax and view_ymin < yl[i] < view_ymax:
            xlvis.append(xl[i])
            ylvis.append(yl[i])

    fig = plt.figure()


    # plot histogram
    plt1 = fig.add_subplot(221)
    plt1.set_axis_off()
    t0 = time.clock()
    zd, xe, ye = np.histogram2d(yl, xl, bins=10, range=[[view_ymin, view_ymax],[view_xmin, view_xmax]], normed=True)
    plt.title('numpy.histogram2d - '+str(time.clock()-t0)+"sec")
    plt.imshow(zd, origin='lower', extent=[view_xmin, view_xmax, view_ymin, view_ymax])
    plt.scatter(xlvis, ylvis)


    # plot density calculated with kdtree
    plt2 = fig.add_subplot(222)
    plt2.set_axis_off()
    xi = np.linspace(view_xmin, view_xmax, 256)
    yi = np.linspace(view_ymin, view_ymax, 256)
    t0 = time.clock()
    zd = grid_density_kdtree(xl, yl, xi, yi, 70)
    plt.title('function of 5 nearest using kdtree\n'+str(time.clock()-t0)+"sec")
    cmap=cm.jet
    A = (cmap(zd/256.0)*255).astype(np.uint8)
    #A[:,:,3] = zd  
    plt.imshow(A , origin='lower', extent=[view_xmin, view_xmax, view_ymin, view_ymax])
    plt.scatter(xlvis, ylvis)

    # gaussian filter
    plt3 = fig.add_subplot(223)
    plt3.set_axis_off()
    t0 = time.clock()
    zd = grid_density_gaussian_filter(view_xmin, view_ymin, view_xmax, view_ymax, 256, 256, zip(xl, yl))
    plt.title('ndi.gaussian_filter - '+str(time.clock()-t0)+"sec")
    plt.imshow(zd , origin='lower', extent=[view_xmin, view_xmax, view_ymin, view_ymax])
    plt.scatter(xlvis, ylvis)

    # boxsum smoothing
    plt3 = fig.add_subplot(224)
    plt3.set_axis_off()
    t0 = time.clock()
    zd = grid_density_boxsum(view_xmin, view_ymin, view_xmax, view_ymax, 256, 256, zip(xl, yl))
    plt.title('boxsum smoothing - '+str(time.clock()-t0)+"sec")
    plt.imshow(zd, origin='lower', extent=[view_xmin, view_xmax, view_ymin, view_ymax])
    plt.scatter(xlvis, ylvis)

if __name__=='__main__':
    generate_graph()
    plt.show()

Answer 1

直方图

直方图方法不是最快的，并且不能区分任意小的点间距和2 * sqrt(2) * b（其中b是bin宽度）之间的区别。

即使您分别构造x个箱和y个箱（O（N）），您仍然需要执行一些ab卷积（每种方法的箱数），对于任何稠密系统，该卷积接近N^2，对于稀疏系统，该卷积甚至更大（在稀疏系统中，ab>；>；N^2）

看上面的代码，您似乎在grid_density()中有一个循环，它在x中的一个循环中运行超过y中的容器数量，这就是为什么您获得O（N^2）性能（尽管如果您已经订购了N，您应该在不同数量的元素上绘制N，那么您只需要每个循环run less code）。

如果你想要一个实际的距离函数，那么你需要开始研究接触检测算法。

接触检测

朴素的接触检测算法在0（N^2）的RAM或CPU时间内出现，但是有一个算法，正确或错误地归于伦敦圣玛丽学院的Munjiza，它运行在线性时间和RAM中。

如果愿意的话，您可以从his book中阅读并自己实现它。

实际上，这段代码是我自己写的

我已经在2D中编写了一个python包装的C实现，它还没有真正准备好投入生产（它仍然是单线程的，等等），但是它将在数据集允许的接近O（N）的地方运行。设置“元素大小”，它充当一个bin大小（代码将调用另一个点的b内的所有内容的交互，有时调用b和2 * sqrt(2) * b之间的交互），给它一个带有x和y属性的对象数组（原生python列表），我的C模块将回调到您选择的python函数，为匹配的元素对运行交互函数。它是为运行接触力DEM模拟而设计的，但它也能很好地解决这个问题。

由于我还没有发布它，因为库中的其他部分还没有准备好，我将不得不给你一个我当前源的zip，但接触检测部分是可靠的。代码是LGPL'd

你需要Cython和一个c编译器才能使它工作，而且它只在*nix环境下测试和工作过，如果你在windows上，你需要the mingw c compiler for Cython to work at all。

安装Cython之后，构建/安装pynet应该是运行setup.py的情况。

您感兴趣的函数是pynet.d2.run_contact_detection(py_elements, py_interaction_function, py_simulation_parameters)（如果希望类元素和模拟参数抛出更少的错误，您应该在同一级别签出它们——请查看archive-root/pynet/d2/__init__.py中的文件以查看类实现，它们是带有有用构造函数的简单数据持有者。）

（当代码准备好进行更一般的发布时，我将用公共mercurial repo更新这个答案…）

Answer 2

这种方法与前面的一些答案大致相同：为每个点增加一个像素，然后用高斯滤波器平滑图像。在我6岁的笔记本电脑上，一张256x256的图像大约需要350毫秒。

import numpy as np
import scipy.ndimage as ndi

data = np.random.rand(30000,2)           ## create random dataset
inds = (data * 255).astype('uint')       ## convert to indices

img = np.zeros((256,256))                ## blank image
for i in xrange(data.shape[0]):          ## draw pixels
    img[inds[i,0], inds[i,1]] += 1

img = ndi.gaussian_filter(img, (10,10))

Answer 3

在纯python中，一个非常简单的实现（用C）可以实时完成，只需几秒钟的时间，就是在屏幕空间计算结果。

算法是

1. 将最后一个矩阵（如256x256）分配为所有零
2. 对于数据集中的每个点，递增相应的单元格
3. 用以单元格为中心的NxN框中的矩阵值之和替换矩阵中的每个单元格。重复此步骤几次。
4. 缩放结果和输出

使用求和表可以使盒和的计算非常快速且独立于N。每次计算只需要对矩阵进行两次扫描。。。总复杂度为O（S+WHP），其中S是点的数目；W，H是输出的宽度和高度，P是平滑过程的数目。

下面是纯python实现的代码（也非常未优化）；有30000个点和256x256个输出灰度图像，计算时间为0.5秒，包括线性缩放到0..255和保存.pgm文件（N=5，4个过程）。

def boxsum(img, w, h, r):
    st = [0] * (w+1) * (h+1)
    for x in xrange(w):
        st[x+1] = st[x] + img[x]
    for y in xrange(h):
        st[(y+1)*(w+1)] = st[y*(w+1)] + img[y*w]
        for x in xrange(w):
            st[(y+1)*(w+1)+(x+1)] = st[(y+1)*(w+1)+x] + st[y*(w+1)+(x+1)] - st[y*(w+1)+x] + img[y*w+x]
    for y in xrange(h):
        y0 = max(0, y - r)
        y1 = min(h, y + r + 1)
        for x in xrange(w):
            x0 = max(0, x - r)
            x1 = min(w, x + r + 1)
            img[y*w+x] = st[y0*(w+1)+x0] + st[y1*(w+1)+x1] - st[y1*(w+1)+x0] - st[y0*(w+1)+x1]

def saveGraph(w, h, data):
    X = [x for x, y in data]
    Y = [y for x, y in data]
    x0, y0, x1, y1 = min(X), min(Y), max(X), max(Y)
    kx = (w - 1) / (x1 - x0)
    ky = (h - 1) / (y1 - y0)

    img = [0] * (w * h)
    for x, y in data:
        ix = int((x - x0) * kx)
        iy = int((y - y0) * ky)
        img[iy * w + ix] += 1

    for p in xrange(4):
        boxsum(img, w, h, 2)

    mx = max(img)
    k = 255.0 / mx

    out = open("result.pgm", "wb")
    out.write("P5\n%i %i 255\n" % (w, h))
    out.write("".join(map(chr, [int(v*k) for v in img])))
    out.close()

import random

data = [(random.random(), random.random())
        for i in xrange(30000)]

saveGraph(256, 256, data)

编辑

当然，密度的定义取决于一个分辨率半径，或者当你到达一个点的时候密度是+inf，而当你没有到达一个点的时候密度是0？

以下是用上述程序制作的动画，只需稍作改动：

1. 使用sqrt(average of squared values)代替sum进行平均通过
2. 对结果进行颜色编码
3. 将结果拉伸为始终使用全色比例
4. 在数据点所在的位置绘制抗锯齿黑点
5. 通过将半径从2增加到40来制作动画

以下39帧动画的总计算时间在PyPy中为5.4秒，在标准Python中为26秒。

Answer 4

直方图

直方图方法不是最快的，并且不能区分任意小的点间距和2 * sqrt(2) * b（其中b是bin宽度）之间的区别。

即使您分别构造x个箱和y个箱（O（N）），您仍然需要执行一些ab卷积（每种方法的箱数），对于任何稠密系统，该卷积接近N^2，对于稀疏系统，该卷积甚至更大（在稀疏系统中，ab>；>N^2）

看上面的代码，您似乎在grid_density()中有一个循环，它在x中的一个循环中运行超过y中的容器数量，这就是为什么您获得O（N^2）性能（尽管如果您已经订购了N，您应该在不同数量的元素上绘制N，那么您只需要每个循环run less code）。

如果你想要一个实际的距离函数，那么你需要开始研究接触检测算法。

接触检测

朴素的接触检测算法在0（N^2）的RAM或CPU时间内出现，但是有一个算法，正确或错误地归于伦敦圣玛丽学院的Munjiza，它在线性时间和RAM中运行。

如果愿意的话，您可以从his book中阅读并自己实现它。

实际上，这段代码是我自己写的

我已经在2D中编写了一个python包装的C实现，它还没有真正准备好投入生产（它仍然是单线程的，等等），但是它将在数据集允许的接近O（N）的地方运行。设置“元素大小”，它充当一个bin大小（代码将调用另一个点的b内的所有内容的交互，有时调用b和2 * sqrt(2) * b之间的交互），给它一个带有x和y属性的对象数组（原生python列表），我的C模块将回调到您选择的python函数，以便为匹配的元素对运行交互函数。它是为运行接触力DEM模拟而设计的，但它也能很好地解决这个问题。

由于我还没有发布它，因为库的其他部分还没有准备好，我将不得不给你一个我当前源的zip，但接触检测部分是可靠的。代码是LGPL'd

你需要Cython和一个c编译器才能使它工作，而且它只在*nix环境下测试和工作过，如果你在windows上，你需要the mingw c compiler for Cython to work at all。

安装Cython之后，构建/安装pynet应该是运行setup.py的情况。

您感兴趣的函数是pynet.d2.run_contact_detection(py_elements, py_interaction_function, py_simulation_parameters)（如果希望类元素和模拟参数抛出更少的错误，您应该在同一级别签出它们——请查看archive-root/pynet/d2/__init__.py中的文件以查看类实现，它们是带有有用构造函数的简单数据持有者。）

（当代码准备好进行更一般的发布时，我将用一个公共的mercurial repo更新这个答案…）