2017-04-21 7 views
1

J'ai deux fichiers de données, chacun d'entre eux contient un grand nombre de points en 3 dimensions (le fichier A stocke environ 50 000 points, le fichier B stocke environ 500 000 points). Mon but est de trouver pour chaque point (a) dans le fichier A le point (b) dans le fichier B qui a la plus petite distance à (a). Je stocke les points dans deux listes comme ceci:Trouver le point 3D le plus proche

Liste A nœuds:

(ID  X  Y   Z) 
[ ['478277', -107.0, 190.5674, 128.1634], 
    ['478279', -107.0, 190.5674, 134.0172], 
    ['478282', -107.0, 190.5674, 131.0903], 
    ['478283', -107.0, 191.9798, 124.6807], 
             ... ] 

Liste B données:

(X  Y  Z  Data) 
[ [-28.102, 173.657, 229.744, 14.318], 
    [-28.265, 175.549, 227.824, 13.648], 
    [-27.695, 175.925, 227.133, 13.142], 
            ...] 

Ma première approche a consisté à parcourons simplement la première et deuxième liste avec une boucle imbriquée et calculer la distance entre tous les points comme ceci:

outfile = open(job[0] + '/' + output, 'wb'); 

dist_min = float(job[5]); 
dist_max = float(job[6]); 

dists = []; 

for node in nodes: 

    shortest_distance = 1000.0; 
    shortest_data = 0.0; 

    for entry in data: 
    dist = math.sqrt((node[1] - entry[0])**2 + (node[2] - entry[1])**2 + (node[3] - entry[2])**2); 
    if (dist_min <= dist <= dist_max) and (dist < shortest_distance): 
     shortest_distance = dist; 
     shortest_data = entry[3]; 

    outfile.write(node[0] + ', ' + str('%10.5f' % shortest_data + '\n')); 

outfile.close(); 

J'ai reconnu que le nombre de boucles que Python doit exécuter est beaucoup trop important (~ 25 000 000 000), j'ai donc dû fixer mon code. J'ai essayé d'abord calculer toutes les distances avec compréhensions de liste mais le code est encore trop lent:

p_x = [row[1] for row in nodes]; 
p_y = [row[2] for row in nodes]; 
p_z = [row[3] for row in nodes]; 

q_x = [row[0] for row in data]; 
q_y = [row[1] for row in data]; 
q_z = [row[2] for row in data]; 

dx = [[(px - qx) for px in p_x] for qx in q_x]; 
dy = [[(py - qy) for py in p_y] for qy in q_y]; 
dz = [[(pz - qz) for pz in p_z] for qz in q_z]; 

dx = [[dxxx * dxxx for dxxx in dxx] for dxx in dx]; 
dy = [[dyyy * dyyy for dyyy in dyy] for dyy in dy]; 
dz = [[dzzz * dzzz for dzzz in dzz] for dzz in dz]; 

D = [[(dx[i][j] + dy[i][j] + dz[i][j]) for j in range(len(dx[0]))] for i in range(len(dx))]; 
D = [[(DDD**(0.5)) for DDD in DD] for DD in D]; 

Pour être honnête, à ce stade, je ne sais pas lequel des deux approches est mieux, de toute façon, aucun des deux possibilités semblent réalisables. Je ne suis même pas sûr s'il est possible d'écrire un code qui calcule toutes les distances dans un temps acceptable. Y a-t-il encore une autre façon de résoudre mon problème sans calculer toutes les distances?

Edit: J'ai oublié de mentionner que je suis en cours d'exécution sur Python 2.5.1 et ne suis pas autorisé à installer ou d'ajouter de nouvelles bibliothèques ...

Répondre

0

Il m'a fallu un peu de réflexion, mais à la fin Je pense que j'ai trouvé une solution pour vous. Votre problème ne réside pas dans le code que vous avez écrit mais dans l'algorithme qu'il implémente. Il existe un algorithme appelé algorithme de Dijkstra et voici l'essentiel: https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm.
Maintenant, ce que vous devez faire est d'utiliser cet algorithme d'une manière intelligente:

créer un nœud S (stand for source). Reliez maintenant les arêtes de S à tous les noeuds du groupe B. Après avoir fait cela, vous devez lier les bords de chaque point b de B à chaque point a de A. Vous devez définir le coût des liens de la source à 0 et l'autre à la distance entre 2 points (uniquement en 3D) . Maintenant, si nous allons utiliser l'algorithme de Dijkstra, la sortie que nous obtiendrons serait le coût pour aller de S à chaque point du graphique (nous ne sommes intéressés que par la distance aux points du groupe A). Donc, puisque le coût est de 0 à chaque point b dans B et S est seulement connecté aux points dans B donc la route vers un point a dans A doit inclure un nœud dans B (en fait exactement un puisque la plus courte distance est un une seule ligne). Je ne suis pas sûr si cela va fixer votre code mais pour autant que je sache, une façon de résoudre ce problème sans calculer toutes les distances n'existe pas et cet algorithme est la meilleure complexité de temps que l'on pourrait espérer.

0

Juste au cas où quelqu'un est interresse dans la solution:

J'ai trouvé un moyen d'accélérer le processus en ne calculant pas toutes les distances:

J'ai créé une liste 3D, ce qui représente une grille dans la espace 3D donné, divisé en X, Y et Z dans une taille de pas donnée (par exemple (Max - Min.)/1000). Puis j'ai itéré sur chaque point 3D pour le mettre dans ma grille. Après cela, j'ai réitéré sur les points de l'ensemble A, en regardant s'il y a des points de B dans le même cube, sinon j'augmenterais le rayon de recherche, donc le processus cherche dans les 26 cubes adjacents des points. Le rayon augmente jusqu'à ce qu'il y ait au moins un point trouvé. La liste qui en résulte est relativement petite et peut être commandée en peu de temps et le point le plus proche est trouvé.

Le temps de traitement est passé à quelques minutes et il fonctionne bien.

p_x = [row[1] for row in nodes]; 
p_y = [row[2] for row in nodes]; 
p_z = [row[3] for row in nodes]; 

q_x = [row[0] for row in data]; 
q_y = [row[1] for row in data]; 
q_z = [row[2] for row in data]; 

min_x = min(p_x + q_x); 
min_y = min(p_y + q_y); 
min_z = min(p_z + q_z); 
max_x = max(p_x + q_x); 
max_y = max(p_y + q_y); 
max_z = max(p_z + q_z); 

max_n = max(max_x, max_y, max_z); 
min_n = min(min_x, min_y, max_z); 

gridcount = 1000; 
step = (max_n - min_n)/gridcount; 

ruler_x = [min_x + (i * step) for i in range(gridcount + 1)]; 
ruler_y = [min_y + (i * step) for i in range(gridcount + 1)]; 
ruler_z = [min_z + (i * step) for i in range(gridcount + 1)]; 

grid = [[[0 for i in range(gridcount)] for j in range(gridcount)] for k in range(gridcount)]; 

for node in nodes: 
    loc_x = self.abatemp_get_cell(node[1], ruler_x); 
    loc_y = self.abatemp_get_cell(node[2], ruler_y); 
    loc_z = self.abatemp_get_cell(node[3], ruler_z); 

    if grid[loc_x][loc_y][loc_z] is 0: 
    grid[loc_x][loc_y][loc_z] = [[node[1], node[2], node[3], node[0]]]; 
    else: 
    grid[loc_x][loc_y][loc_z].append([node[1], node[2], node[3], node[0]]); 

for entry in data: 
    loc_x = self.abatemp_get_cell(entry[0], ruler_x); 
    loc_y = self.abatemp_get_cell(entry[1], ruler_y); 
    loc_z = self.abatemp_get_cell(entry[2], ruler_z); 

    if grid[loc_x][loc_y][loc_z] is 0: 
    grid[loc_x][loc_y][loc_z] = [[entry[0], entry[1], entry[2], entry[3]]]; 
    else: 
    grid[loc_x][loc_y][loc_z].append([entry[0], entry[1], entry[2], entry[3]]); 

out = []; 
outfile = open(job[0] + '/' + output, 'wb'); 

for node in nodes: 
    neighbours = []; 
    radius = -1; 
    loc_nx = self.abatemp_get_cell(node[1], ruler_x); 
    loc_ny = self.abatemp_get_cell(node[2], ruler_y); 
    loc_nz = self.abatemp_get_cell(node[3], ruler_z); 
    reloop = True; 

    while reloop: 
    if neighbours: 
     reloop = False; 
    radius += 1; 
    start_x = 0 if ((loc_nx - radius) < 0) else (loc_nx - radius); 
    start_y = 0 if ((loc_ny - radius) < 0) else (loc_ny - radius); 
    start_z = 0 if ((loc_nz - radius) < 0) else (loc_nz - radius); 
    end_x = (len(ruler_x) - 1) if ((loc_nx + radius + 1) > (len(ruler_x) - 1)) else (loc_nx + radius + 1); 
    end_y = (len(ruler_y) - 1) if ((loc_ny + radius + 1) > (len(ruler_y) - 1)) else (loc_ny + radius + 1); 
    end_z = (len(ruler_z) - 1) if ((loc_nz + radius + 1) > (len(ruler_z) - 1)) else (loc_nz + radius + 1); 

    for i in range(start_x, end_x): 
     for j in range(start_y, end_y): 
     for k in range(start_z, end_z): 
      if not grid[i][j][k] is 0: 
      for grid_entry in grid[i][j][k]: 
       if not isinstance(grid_entry[3], basestring): 
       neighbours.append(grid_entry); 
    dists = []; 
    for n in neighbours: 
    d = math.sqrt((node[1] - n[0])**2 + (node[2] - n[1])**2 + (node[3] - n[2])**2); 
    dists.append([d, n[3]]); 
    dists = sorted(dists); 

    outfile.write(node[0] + ', ' + str(dists[0][-1]) + '\n'); 

outfile.close(); 

Fonction pour obtenir la position d'un point:

def abatemp_get_cell(self, n, ruler): 
    for i in range(len(ruler)): 
     if i >= len(ruler): 
     return False; 
     if ruler[i] <= n <= ruler[i + 1]: 
     return i; 

La variable de gridcount donne l'occasion de fixer le processus, avec une petite gridcount le processus de tri des points en la grille est très rapide, mais les listes de voisins dans la boucle de recherche deviennent plus grandes et plus de temps est nécessaire pour cette partie du processus. Avec un gros gridcount plus de temps est nécessaire au début, mais la boucle tourne plus vite. Le seul problème que j'ai maintenant est le fait qu'il y a des cas où le processus a trouvé des voisins mais il y a d'autres points, qui ne sont pas encore trouvés, mais sont plus proches du point (voir picture). Jusqu'à présent, j'ai résolu ce problème en augmentant le rayon de recherche une autre fois quand il y a déjà des voisins. Et encore, j'ai des points qui sont plus proches mais pas dans la liste des voisins, bien que ce soit une très petite quantité (92 sur 100 000). Je pourrais résoudre ce problème en incrémentant le rayon deux fois après avoir trouvé des voisins, mais cette solution ne semble pas très intelligente. Peut-être que vous avez une idée ...

Ceci est le premier brouillon du processus, je pense qu'il sera possible de l'améliorer encore plus, juste pour vous donner une idée de son fonctionnement ...

0

jeter un oeil à ce générique structure de données 3D:

https://github.com/m4nh/skimap_ros

il dispose d'une fonction RadiusSearch très rapide tout prêt à être utilisé. Cette solution (similaire à Octree mais plus rapide) vous évite d'abord de créer la grille régulière (vous n'avez pas besoin de fixer la taille MAX/MIN le long de chaque axe) et vous économisez beaucoup de mémoire