Méthode efficace de calcul de la densité de points irrégulièrement espacés

Question

J'essaie de générer des images de superposition de cartes qui permettraient d'identifier les points chauds, c'est-à-dire les zones de la carte qui présentent une forte densité de points de données. Aucune des approches que j'ai essayées n'est assez rapide pour mes besoins. Note : J'ai oublié de mentionner que l'algorithme devrait fonctionner correctement dans des scénarios de zoom faible ou élevé (ou de densité de points de données faible ou élevée).

J'ai cherché dans les bibliothèques numpy, pyplot et scipy, et le plus proche que j'ai pu trouver était numpy.histogram2d. Comme vous pouvez le voir dans l'image ci-dessous, la sortie de l'histogramme2d est plutôt grossière. (Chaque image comprend des points superposés à la carte thermique pour une meilleure compréhension)

Ma deuxième tentative a consisté à itérer sur tous les points de données, puis à calculer la valeur du point chaud en fonction de la distance. Cela a produit une image de meilleure qualité, mais c'est trop lent pour être utilisé dans mon application. Puisque c'est O(n), ça fonctionne bien avec 100 points, mais ça explose quand j'utilise mon jeu de données réel de 30 000 points.

Ma dernière tentative a été de stocker les données dans un KDTree, et d'utiliser les 5 points les plus proches pour calculer la valeur du point chaud. Cet algorithme est O(1), donc beaucoup plus rapide avec de grands ensembles de données. Ce n'est pas encore assez rapide, il faut environ 20 secondes pour générer un bitmap de 256x256, et j'aimerais que cela se fasse en 1 seconde environ.

Modifier

La solution de lissage boxsum fournie par 6502 fonctionne bien à tous les niveaux de zoom et est beaucoup plus rapide que mes méthodes originales.

La solution du filtre gaussien suggérée par Luke et Neil G est la plus rapide.

Vous pouvez voir les quatre approches ci-dessous, en utilisant 1000 points de données au total, au zoom 3x il y a environ 60 points visibles.

Le code complet qui génère mes 3 essais originaux, la solution de lissage boxsum fournie par 6502 et le filtre gaussien suggéré par Luke (amélioré pour mieux gérer les bords et permettre le zoom) est ici :

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

Cette approche va dans le sens de certaines réponses précédentes : incrémenter un pixel pour chaque point, puis lisser l'image avec un filtre gaussien. Une image de 256x256 s'exécute en environ 350 ms sur mon ordinateur portable de 6 ans.

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

Une mise en œuvre très simple qui pourrait être réalisée (avec C) en temps réel et qui ne prend que quelques fractions de seconde en python pur consiste à calculer le résultat dans l'espace écran.

L'algorithme est

1. Allouer la matrice finale (par exemple 256x256) avec tous les zéros.
2. Pour chaque point du jeu de données, incrémentez la cellule correspondante.
3. Remplace chaque cellule de la matrice par la somme des valeurs de la matrice dans une boîte NxN centrée sur la cellule. Répétez cette étape plusieurs fois.
4. Résultat et sortie de l'échelle

Le calcul de la somme des boîtes peut être rendu très rapide et indépendant de N en utilisant une table de somme. Chaque calcul ne nécessite que deux balayages de la matrice... la complexité totale est O(S + W H P) où S est le nombre de points ; W, H sont la largeur et la hauteur de la sortie et P est le nombre de passages de lissage.

Voici le code d'une implémentation purement python (également très peu optimisée) ; avec 30000 points et une image en niveaux de gris de 256x256 en sortie, le calcul est de 0.5sec incluant la mise à l'échelle linéaire à 0..255 et la sauvegarde d'un fichier .pgm (N = 5, 4 passes).

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)

Modifier

Bien sûr, la définition même de la densité dans votre cas dépend d'un rayon de résolution, ou la densité est-elle simplement +inf lorsque vous atteignez un point et zéro lorsque vous ne le faites pas ?

Voici une animation réalisée à l'aide du programme ci-dessus, avec seulement quelques modifications cosmétiques :

1. utilisé sqrt(average of squared values) au lieu de sum pour le passage de la moyenne
2. codage couleur des résultats
3. étirer le résultat pour toujours utiliser l'échelle de couleurs complète
4. points noirs anti-calculés dessinés où les points de données sont
5. fait une animation en incrémentant le rayon de 2 à 40

Le temps de calcul total des 39 images de l'animation suivante avec cette version cosmétique est de 5,4 secondes avec PyPy et de 26 secondes avec Python standard.

Answer 3

Histogrammes

La méthode de l'histogramme n'est pas la plus rapide, et ne peut pas faire la différence entre une séparation arbitrairement petite de points et 2 * sqrt(2) * b (où b est la largeur de la case).

Même si vous construisez les bins x et y séparément (O(N)), vous devez toujours effectuer une convolution ab (nombre de bins dans chaque sens), ce qui est proche de N^2 pour tout système dense, et encore plus grand pour un système clairsemé (bien, ab >> N^2 dans un système clairsemé.)

En regardant le code ci-dessus, il semble que vous ayez une boucle en grid_density() qui passe sur le nombre de cases de y à l'intérieur d'une boucle du nombre de cases de x, ce qui explique pourquoi vous obtenez des performances O(N^2) (bien que si vous êtes déjà de l'ordre de N, ce que vous devriez tracer sur différents nombres d'éléments pour voir, alors vous devrez simplement exécuter moins de code par cycle).

Si vous voulez une fonction de distance réelle, vous devez commencer à examiner les algorithmes de détection des contacts.

Détection des contacts

Les algorithmes naïfs de détection des contacts sont de l'ordre de O(N^2) en temps de RAM ou de CPU, mais il existe un algorithme, attribué à tort ou à raison à Munjiza au St. Mary's College de Londres, qui fonctionne en temps linéaire et en RAM.

vous pouvez vous renseigner sur ce sujet et le mettre en œuvre vous-même à partir de son livre si vous voulez.

J'ai écrit ce code moi-même, en fait

J'ai écrit une implémentation en C enveloppée de python de ceci en 2D, qui n'est pas vraiment prête pour la production (c'est toujours un thread unique, etc.) mais elle fonctionnera aussi près de O(N) que votre ensemble de données le permettra. Vous définissez la "taille de l'élément", qui agit comme une taille de bac (le code appellera des interactions sur tout ce qui se trouve à l'intérieur de b d'un autre point, et parfois entre b et 2 * sqrt(2) * b ), donnez-lui un tableau (liste native python) d'objets avec une propriété x et y et mon module C fera un appel à une fonction python de votre choix pour exécuter une fonction d'interaction pour les paires d'éléments appariées. il est conçu pour exécuter des simulations DEM de force de contact, mais il fonctionnera bien pour ce problème aussi.

Comme je ne l'ai pas encore publié, parce que les autres parties de la bibliothèque ne sont pas encore prêtes, je vais devoir vous donner un zip de ma source actuelle mais la partie détection de contact est solide. Le code est sous LGPL.

Vous aurez besoin de Cython et d'un compilateur c pour le faire fonctionner, et il n'a été testé et ne fonctionne que sous les environnements *nix, si vous êtes sous Windows vous aurez besoin de le compilateur mingw c pour Cython ne fonctionne pas du tout .

Une fois que Cython est installé, construire/installer pynet devrait être un cas d'exécution de setup.py.

La fonction qui vous intéresse est pynet.d2.run_contact_detection(py_elements, py_interaction_function, py_simulation_parameters) (et vous devriez vérifier les classes Element et SimulationParameters au même niveau si vous voulez qu'il jette moins d'erreurs - regardez dans le fichier à l'adresse suivante archive-root/pynet/d2/__init__.py pour voir les implémentations des classes, ce sont des supports de données triviaux avec des constructeurs utiles).

(Je mettrai à jour cette réponse avec un repo mercuriel public lorsque le code sera prêt pour une diffusion plus générale...)

Answer 4

Votre solution est correcte, mais un problème évident est que vous obtenez des régions sombres malgré l'existence d'un point en plein milieu de celles-ci.

Je centrerais plutôt une gaussienne à n dimensions sur chaque point et j'évaluerais la somme sur chaque point que vous voulez afficher. Pour le réduire en temps linéaire dans le cas courant, utilisez query_ball_point pour ne prendre en compte que les points situés dans un intervalle de deux écarts types.

Si vous trouvez que le KDTree est vraiment lent, pourquoi ne pas appeler query_ball_point une fois tous les cinq pixels avec un seuil légèrement plus élevé ? Cela ne fait pas trop de mal d'évaluer un peu trop de gaussiennes.

Answer 5

Vous pouvez le faire avec une convolution séparable 2D (scipy.ndimage.convolve1d) de votre image originale avec un noyau de forme gaussienne. Avec une taille d'image de MxM et une taille de filtre de P, la complexité est O(PM^2) en utilisant un filtrage séparable. La complexité "Big-Oh" est sans doute plus grande, mais vous pouvez profiter des opérations efficaces de numpy sur les tableaux qui devraient accélérer considérablement vos calculs.