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

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

这种方法与以前的一些答案一致:为每个点增加一个像素,然后用高斯滤波器平滑图像.在我6岁的笔记本电脑上,256x256图像在大约350ms内运行.

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 2

一个非常简单的实现可以实时完成(使用C)并且在纯python中只需要几分之一秒就可以计算屏幕空间中的结果.

算法是

1. 用全零分配最终矩阵(例如256x256)
2. 对于数据集中的每个点,递增相应的单元格
3. 将单元格中的每个单元格替换为以单元格为中心的NxN框中矩阵值的总和.重复此步骤几次.
4. 缩放结果和输出

通过使用和表,可以非常快速地并且在N上独立地计算箱和.每次计算只需要对矩阵进行两次扫描......总复杂度为O(S + W H P),其中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而在你没有时是零？

以下是使用上述程序构建的动画,只有一些修饰:

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

使用此化妆品版本的以下动画的39帧的总计算时间使用PyPy为5.4秒,使用标准Python为26秒.