已知若干個點的坐標（0<x<L, 0<y<L), 如何高效地找出所有相對距離小於r的點對？

02-03

r，L為定值，且r&<0.1*L。因為點的數目大（幾千個點），所以不希望用逐個檢查某點與所有其他點的距離的辦法。跪求大神提供可以高效解決這問題的思路。

剛剛學了divide and conqure演算法Closest Pair: A Divide-and-Conquer Approach，用來找點集中距離最小的兩個點，改一改應該可以用來做這個，只需要把最小距離d固定為1就可以了吧。問題詳細的wiki描述Nearest neighbor search

把點分為兩組，所有小於一的數目等於第一組中小於一的+第二組中小於一的+第一二組之間小於一的，如此遞歸下去是 $O(nlog n)$ 複雜度。這個演算法裡面最tricky的是如何計算一二組間距離小於一的點對，這個很多大學的演算法課件都有，我剛剛學就不獻醜了。

不過點數不是太多所以log n的優勢不知道能不能體現出來。

=======================

Wiki上給出了一個貌似是 $O(n)$ 的演算法，對所有點做rounding，保留整數位，相當於畫了1x1的網格，然後做hash，就是網格標號，所有rounding之後在一個網格內的計算一下實際距離，相鄰網格所有點計算一下實際距離。

幾千個點，純粹枚舉求距離O(n^2)，單線程每秒可以處理數百次。

排序後稍微優化一下，平均情況變得好很多，再快一到兩個數量級。

沒有任何特異性，可以丟給GPU算，每秒可以處理幾十萬輪。

可以將L * L的區域劃分成r * r的網格，每個點只用檢查與其所處網格相鄰網格內的點與自己的距離。假定r &<&< L且點期望是均勻分布的，那麼可以用O(L^2/r^2)的空間換取O(L^2/r^2)的加速。

不過不確定Matlab實現起來方不方便……

kdtree

瀉藥，不怎麼會matlab。不過，你只要堅持只用向量矩陣運算就不會慢。

講真的，就算是最慢的寫法（兩重循環），配上matlab那麼垃圾的runtime，也不會多慢的，有你糾結這個問題的時間，程序早就跑完了……

演算法上：

1. 先對數據排序

2. 利用數據有序的性質減少比較次數

此外：

我猜測這步計算只是一個大演算法中的一部分，一種可考慮的思路是在演算法之前的步驟中預處理一些可以簡化本部分計算的信息。

ps. 幾千個點其實很少的，感覺可以直接暴力去求的。

二維問題的話可以先生成三角形網格，然後用最短距離的方法做。我記得有文獻說這個的。這類問題好像叫nearest neighborhood problem。在計算力學領域蠻常見的。

如果小心地不使用循環的話，感覺幾千個點不會很慢呀？

比如，設X, Y是存儲坐標的向量。

[X1, X2]=meshgrid(X, X);

[Y1, Y2]=meshgrid(Y, Y);

R=sqrt((X1-X2).^2+(Y1-Y2).^2);

R(R&R(R&>=0)=0;

R=-R;

現在矩陣R所有等於1的元素的行列指標對就是所求點對。

R是對稱矩陣，所以這個演算法多了一倍運算量，但寫著方便...

如果這樣還很慢，可以先分別卡x y坐標，這樣在避免計算乘法的情況下先排除掉大部分可能性，然後再用上面的方法進行挑選。...其實還是很笨的辦法。

上面的代碼我沒跑，也許有錯。

Kdtree或者Cell list

C++庫就是flann，nanoflann，point cloud library

看到了很多很複雜的方法，都太裝逼了，說點實際的，可以寫成點算的。

假如距離為x，下面我個人的一個想法：

Step One

先做外包矩形判斷，超出距離x直接就排除，小於x的進入Step Two。

Step Two

再進行計算歐式幾何距離判斷，或者用大地線距離計算公式Geodetic distance的計算公式計算。

注意

由於距離是兩個點相對的，循環其實可以減半，如

(1,2,3,4) (2,3,4) (3,4)

在性能上可以進一步優化一下。

四分樹，逐點在周圍的2r*2r正方形內枚舉（其實幾千個點你有功夫發這個問題寫個暴力都跑完了

如果單純是一維的，那很簡單，用二分法，用你要比較的值和整體量的中間值進行比較，如果是二維的話，思路也是一樣的，關鍵看你怎麼設計。這樣的話，一般來說，效率應該會快很多

1.背景網格。把所有點組織到網格里，搜索時只需要考慮鄰近網格內的點。

如果點的位置逐步更新，那麼每步都要更新網格。網格尺寸與L的關係需要注意，可以優化。

2.如果只是千量級，另一種方式是對每個點預先建立neighbor list。

list包含兩個距離r1(=L)和r2，r2表示接下來一段時間內可能進入r1，這樣沒必要每步更新所有list。

首先，幾千個點的話你O(n^2)枚舉任意兩點O(1)判定的話還是很快的呀，如果數據範圍再大一些的話，kd-tree大法吼啊。

如果是python的話有一個自帶的叫map的function，主要define一個小function，小於你的設定值return 1，else return 0，然後數有多少個0就ok了。不過matlab就不知道了。

幾千個點的意思就是小於一萬嗎？10000*10000的矩陣還是能開出來的，matlab用矩陣算就不會太慢。用筆記本試了一下5000個點1秒，10000個點5秒，不知道你想要多快。

L=100;

r=5;

n=10000;

p=rand(n,2)*L;

d=(ones(n,1)*p(:,1)-p(:,1)*ones(1,n)).^2+(ones(n,1)*p(:,2)-p(:,2)*ones(1,n)).^2;

idx1=find(d&pair1=[floor((idx1-1)/n)+1 mod(idx1-1,n)+1];