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

我正在尝试生成地图叠加图像，以帮助识别热点，即地图上具有高密度数据点的区域。 我所尝试的方法都不足以满足我的需求。 注意：我忘记提及algorithm在低和高缩放情况（或低和高数据点密度）下都能正常工作。

我通过numpy，pyplot和scipy库来查看，最接近我能find的是numpy.histogram2d。 正如你在下面的图片中看到的，histogram2d的输出是比较粗糙的。 （每个图像包含覆盖热图的点以便更好地理解）

我的第二个尝试是迭代所有的数据点，然后计算热点值作为距离的函数。 这产生了一个更好看的形象，但它在我的应用程序中使用太慢。 由于它是O（n），所以它可以正常工作100点，但是当我使用我的30000点的实际数据集时会发生。

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

编辑

由6502提供的boxsum平滑解决scheme在所有缩放级别都能正常运行，比我的原始方法快得多。

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

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

完整的代码，生成我原来的3次尝试，6502提供boxsum平滑解决scheme和卢克build议的高斯filter（改进处理边缘更好，允许放大）在这里：

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岁的笔记本电脑上运行一个256×256的图像大约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 2

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

algorithm是

1. 分配最终matrix（例如256×256）全零
2. 对于数据集中的每个点增加相应的单元格
3. 将matrix中的每个单元格replace为以单元格为中心的NxN框中的matrix值的总和。 重复这个步骤几次。
4. 缩放结果和输出

盒子和的计算可以通过使用和表非常快速和独立于N. 每个计算只需要两次matrix扫描…总的复杂度是O（S + W H P）其中S是点的数量; W，H是输出的宽度和高度，P是平滑通道的数量。

下面是一个纯python实现的代码（也非常优化）。 与30000点和一个256×256输出灰度图像计算是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时，密度刚刚达到一个点，当你不达到时，密度为零。

以下是上述程序构build的animation，仅进行了一些修改：

1. 使用sqrt(average of squared values)而不是平均值
2. 对结果进行颜色编码
3. 将结果拉伸至始终使用全色标
4. 绘制数据点所在的抗锯齿黑点
5. 通过将半径从2增加到40来制作animation

以下化妆版本的39帧animation的总计算时间是PyPy 5.4秒和标准Python 26秒。

Answer 3

直方图

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

即使你分别构造x b和y bins（O（N）），你仍然必须执行一些ab卷积（每个bins的数目），对于任何密集系统来说，接近N ^ 2，甚至更大一个稀疏的系统（在一个稀疏系统中，ab >> N ^ 2）。

看上面的代码，你似乎在grid_density()有一个循环，它运行在x的bin数目循环内的y中的bin数量，这就是为什么你得到O（N ^ 2）的性能（尽pipe如果你已经订购了N，你应该在不同数量的元素上绘制，那么你将不得不在每个周期运行更less的代码 ）。

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

接触检测

天真的接触检测algorithm在O（N ^ 2）的RAM或CPU时间内出现，但有一个algorithm，正确或错误地归因于在伦敦圣玛丽学院的Munjiza，它运行在线性时间和内存。

如果你愿意的话，你可以从他的书中阅读并实施它。

事实上，我自己写了这个代码

我用2D编写了一个python包装的C实现，这个实现并没有准备好用于生产（它仍然是单线程的，等等），但是它将在你的数据集允许的最接近O（N）的地方运行。 你设置了“元素大小”，它作为一个bin大小（代码将调用另一点b的所有内容的交互，有时在b和2 * sqrt(2) * b ）之间的所有内容，给它一个数组（本地python列表）的对象与x和y属性和我的C模块将callback到您select的Python函数运行匹配的元素对交互function。 它被devise用于运行接触力DEM模拟，但是它也可以很好地处理这个问题。

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

你需要使用Cython和ac编译器才能使它工作，而且它只是在* nix环境下进行testing和工作，如果你在windows上，你需要使用mingw c编译器来运行Cython 。

一旦Cython的安装，构build/安装pynet应该是运行setup.py的情况。

你感兴趣的函数是pynet.d2.run_contact_detection(py_elements, py_interaction_function, py_simulation_parameters) （如果你想要抛出更less的错误，你应该检查出同样级别的Element和SimulationParameters类， archive-root/pynet/d2/__init__.py来查看类实现，它们是有用的构造函数的简单数据持有者。）

（我会更新这个答案，当公共代码已经准备好更普遍的发布…）

Answer 4

你的解决scheme是可以的，但是一个明显的问题是，尽pipe在它们中间有一个点，

相反，我会在每个点上以n维高斯为中心，并计算每个要显示的点的总和。 为了在常见情况下将其减less到线性时间，请使用query_ball_point仅考虑几个标准偏差内的点。

如果你发现他KDTree真的很慢，为什么不调用query_ball_point每五个像素一个稍大的阈值？ 评估一些太多的高斯人并没有太多的伤害。

Answer 5

你可以用一个高斯形状的内核对原始图像进行2D分离卷积（scipy.ndimage.convolve1d）。 在图像大小为M×M，滤波器大小为P的情况下，使用可分离滤波的复杂度为O（PM ^ 2）。 “大哦”的复杂性无疑更大，但你可以利用numpy的高效数组操作，这应该大大加快你的计算。

Answer 6

只需要注意， histogram2d函数应该可以正常工作。 你玩过不同的垃圾桶大小吗？ 你的初始histogram2d图似乎只是使用默认的bin大小…但没有理由期望默认大小给你你想要的表示。 话虽如此，其他许多解决scheme也令人印象深刻。