2017-04-21 66 views
1

我有兩個數據文件,每個文件都包含的3維點一個大數目(文件A存儲約50000點,文件B存儲大約50萬個點)。我的目標是爲文件A中的每個點(a)找出文件B中與(a)具有最小距離的點(b)。我在兩個列表點存儲這樣的:查找最近的三維點

列表A 節點

(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], 
             ... ] 

列表B 數據

(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], 
            ...] 

我的第一種方法是簡單地通過第一循環和第二列表嵌套循環,並計算這樣每點之間的距離:

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(); 

我認爲,Python有運行循環的量太大了(〜250億),所以我不得不來固定我的代碼。我想先計算與列表理解所有的距離,但代碼仍然太慢:

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]; 

說實話,在這一點上,我不知道這兩種方法的比較好,反正沒有一個兩種可能性似乎可行。我甚至不確定是否可以編寫一個代碼,在可接受的時間內計算所有距離。是否有另一種方法來解決我的問題,而不計算所有的距離?

編輯:我忘了提,我上的Python 2.5.1運行,我不可以安裝或添加任何新的圖書館...

回答

0

我花了一點心思,但在最後我想我找到了一個解決方案。 你的問題不在你寫的代碼中,而是在它實現的算法中。 有一種叫做Dijkstra算法的算法,這裏是它的要點:https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm
現在你需要做的是以一種聰明的方式使用這個算法:

創建一個節點S(代表源代碼)。 現在將S中的邊緣鏈接到B組中的所有節點。 後你做了,你應該從每一個點b在B到每個點在A.鏈接邊緣 你應該設置從源頭上鍊接到0成本和其他的距離2點之間(僅適用於3D) 。 現在,如果我們將使用Dijkstra算法,我們將得到的輸出將是從S到圖中每個點的傳播成本(我們只關心到A組中點的距離)。 如此以來,成本是0到每個點b在B和S僅在B相連接的點,以便在路上以A任意點必須包括在B(節點實際上正好一個自對點之間的最短距離爲單線)。

我不知道這是否會繫好你的代碼,但據我所知,這是一種在不計算所有距離不存在解決這個問題,這個算法是最好的時間複雜度的指望。

0

萬一有人在溶液interrested:

我發現了一種通過不計算所有的距離,加快了整個過程:

我創建了一個三維列表,較一格給定3D空間,按給定步長(例如(最大 - 最小)/ 1,000)分爲X,Y和Z兩種。然後我遍歷每個3D點將其放入我的網格。之後,我再次遍歷集合A的點,看看在同一個立方體中是否有來自B的點,如果不是,我會增加搜索半徑,因此該過程正在相鄰的26個立方體中尋找點。半徑正在增加,直到找到至少一個點。結果列表相對較小,可以在短時間內訂購,並找到最近的點。

處理時間縮短到幾分鐘,工作正常。

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(); 

函數來獲得一個點的位置:

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; 

gridcount變量提供了一個固定的過程中,用小gridcount點排序進入的過程中,機會網格非常快,但搜索循環中的鄰居列表變得更大,並且需要更多時間來完成這部分過程。有一個大的網格計數在開始時需要更多的時間,但循環運行速度更快。

我現在唯一的問題是,有些情況下,進程找到鄰居,但還有其他點,但尚未找到,但更接近點(見picture)。到目前爲止,我已經解決了這個問題,通過增加搜索半徑的另一個時間已經neighbours。儘管數量非常少(約100,000人中有92人),但我仍然得到更接近但不在鄰居名單中的分數。在找到鄰居之後,我可以通過增加半徑來解決這個問題,但是這個解決方案看起來不太聰明。也許你們有一個想法...

這是該過程的第一個工作草案,我認爲它可以改進它甚至更多,只是爲了讓你知道它是如何工作的...

0

看看這個普通的3D數據結構:

https://github.com/m4nh/skimap_ros

它有一個非常快的RadiusSearch功能只是隨時可以使用。這種解決方案(類似於八叉樹,但速度更快)避免了您首先創建常規網格(您不必固定每個軸的最大/最小尺寸),並節省大量內存