如何在数据数组中找到距离给定点最近的点？
例如，假设我有一个带有x、y和z坐标的3D点数组A，以及一个点(x_p，y_p，z_p)。如何找到A中距离(x_p，y_p，z_p)最近的点？
据我所知，最慢的方法是使用线性搜索。是否有更好的解决方案？
可以添加任何辅助数据结构。
9个回答
31
你可以使用Octree来组织你的要点。然后，你只需要搜索一个小的子集。
Octree是一种相当简单的数据结构，你可以自己实现（这将是一个有价值的学习经验），或者你可以找到一些有用的库来帮助你入门。
Octree是一种相当简单的数据结构，你可以自己实现（这将是一个有价值的学习经验），或者你可以找到一些有用的库来帮助你入门。
- dkamins
3
9这里提出的算法仅在需要为大量点重复搜索最近邻时才有效。如果我们只需要一个点的信息，线性搜索更加高效。 - efficiencyIsBliss
2进一步解释我的评论，构建树本身（KD Tree或OC Tree）将比线性更糟。我不确定OC树如何，但是KD树需要O（NlogN）。因此，对于单个查询，线性搜索更好。 - efficiencyIsBliss
1@efficiencyIsBliss 但是你知道为了单个点的kNN计算而计算数十万次乘法的代价吗？这可以通过Octree中的几个距离计算来减少，而建立Octree的成本远远低于线性搜索中距离计算的众多乘法开销，即使只针对一个点，因为现在点云很容易就能有超过100k个点。 - Gab是好人
21
- moinudin
6
我建议使用KD树，这也适用于最近邻搜索。
- Vedang Joshi
3
我需要在实时环境中进行许多最近邻居搜索，因此我想出了一种更简单和更快的算法。将所有点复制到d个列表中，其中d是空间的维度，在您的情况下为3。按其维度对这三个列表进行排序。这需要d(nlog(n))时间。至于数据结构，就是维护这些每个点在每个维度上的适当排序列表。关键是根据定义，某个方向的距离必须小于或等于欧几里得距离。因此，如果一个方向的距离大于我们当前已知的最近点的距离，则该点不能更接近，更重要的是，该方向上的所有点都不能更接近。一旦对于我们拥有的2 * d个方向成立，我们就可以找到最近的点。对于每个特定元素，我们可以二分查找排好序的列表，以找到所需点可能在两个不同维度中的最近位置。数学上，我们知道，如果+x-x+y-y（其他维度易于添加）方向上的距离超过到点的最小已知欧几里得距离，则该点必须超过该距离，并且由于它是已排序的数组，因此根据定义，当我们在该方向上超过该距离时，我们知道我们可以中止该方向，因为在该方向上不可能有更好的答案。但是，当我们沿着这四个方向扩展时，我们可以减少m的值，因为它等于我们找到的最接近点的欧几里得距离。因此，我们只需对每个轴进行排序，并按照该轴排序即可。
然后查询列表：
- 对于每个列表，我们进行二分查找（dlog(n)）。 - 找到当前的最小距离m（初始可以是无限大）。 - 对于每个列表，我们向正负方向移动。 - 对于我们拥有的2 * d个方向，
- 当我们发现更接近的点时，我们遍历列表并降低m的值。
- 当一个方向被证明是数学上没有意义的时，我们停止搜索该方向。 - 当没有方向剩余时，我们已经找到了最近的点。
我们需要在每个方向上的列表中找到正在搜索的点。我们进行二分搜索以保持时间复杂度为log(n)。然后我们有我们当前的最佳距离（可能是无穷大），然后我们在每个可用方向中移动。随着我们发现新点，我们更新到目前为止我们找到的最近点。关键是，一旦仅在一方向上的距离大于我们当前已知的最近点，我们就会停止搜索。
如果我们已知一个最近距离为13的点，那么当沿着某一方向的距离超过我们所知的最近距离时，我们可以立即停止在+x、-x、+y、-y方向上的检查。因为如果它比我们当前的m值更远，那么所有剩余的+x值都可以被数学证明为更远。随着我们得到越来越好的最近点，我们需要搜索的空间量会越来越小。
如果我们在某个方向上没有找到点，那么该方向就结束了。如果沿着线的某一维度到一个点的距离本身大于m，那么该方向也结束了。
当所有已经被证明只有更远的点的方向都完成时，解决方案是m。
由于我们逐渐减小m，每个维度所需的距离整体上迅速下降，尽管像所有算法一样，在更高的维度上下降得不太快。但是，如果仅在一个维度中的距离大于到目前为止最佳距离，那么在该方向上的其余点必然不能更好。
在时间复杂度方面，它似乎与更好的算法相当。但是，在数据结构的简单性方面，这个算法显然是胜出的。有很多其他属性使得这个算法成为一个严肃的竞争者。当你更新东西时，你可以用非常好的性能重新排序列表，因为你往往在排序已经排序或近乎排序的列表。你正在迭代数组。实际性能方面，大多数数据结构都很差。通常由于缓存和内存布局的原因，我们应该对此视而不见，但这却很重要。与您当前相关数据旁边的数据访问速度要快得多。如果我们已经知道要查找的点在列表中的位置，我们甚至可以更快地解决它（因为我们不必用二分搜索找到它）。还有其他允许重复利用上一次迭代信息的技巧。而额外的维度基本上是免费的（保存值不会更快地收敛，但这是因为球体中比半径相同的圆中具有更多随机分布的点）。
更新：关于评论中的k-points问题。您会注意到很少有变化。唯一相关的是，如果发现点v小于当前m（findingDistanceMaxSq），则将该点添加到堆中，并将m的值设置为查找位置和第k个元素之间的欧几里得距离。算法的常规版本可以看作是k = 1的情况。当v被发现更接近时，我们搜索我们想要的1个元素并更新m等于唯一的（k = 1）元素。
请记住，我只在距离平方形式下进行距离比较，因为我只需要知道它是否更远，并且不浪费时钟周期在平方根函数上。
我知道有一个完美的数据结构用于存储k个元素在大小限制堆中。显然，数组插入对此不是最优解。但是，除了太多的java依赖api外，那个特定类别没有一个。但是，鉴于您的k很可能不那么大，因此您实际上不会注意到这一点。但是，它确实使点插入的时间复杂度存储在k时间中。还有诸如缓存元素与查找点之间距离之类的事情。
最后，也可能是最紧迫的是，我将用于测试代码的项目正在转换中，因此我还没有测试过这一点。但是，它肯定展示了您如何做到这一点：您存储到目前为止最好的k个结果，并使m等于到第k个最接近点的距离。——其他所有内容保持不变。
示例源代码。
然后查询列表：
- 对于每个列表，我们进行二分查找（dlog(n)）。 - 找到当前的最小距离m（初始可以是无限大）。 - 对于每个列表，我们向正负方向移动。 - 对于我们拥有的2 * d个方向，
- 当我们发现更接近的点时，我们遍历列表并降低m的值。
- 当一个方向被证明是数学上没有意义的时，我们停止搜索该方向。 - 当没有方向剩余时，我们已经找到了最近的点。
我们需要在每个方向上的列表中找到正在搜索的点。我们进行二分搜索以保持时间复杂度为log(n)。然后我们有我们当前的最佳距离（可能是无穷大），然后我们在每个可用方向中移动。随着我们发现新点，我们更新到目前为止我们找到的最近点。关键是，一旦仅在一方向上的距离大于我们当前已知的最近点，我们就会停止搜索。
如果我们已知一个最近距离为13的点，那么当沿着某一方向的距离超过我们所知的最近距离时，我们可以立即停止在+x、-x、+y、-y方向上的检查。因为如果它比我们当前的m值更远，那么所有剩余的+x值都可以被数学证明为更远。随着我们得到越来越好的最近点，我们需要搜索的空间量会越来越小。
如果我们在某个方向上没有找到点，那么该方向就结束了。如果沿着线的某一维度到一个点的距离本身大于m，那么该方向也结束了。
当所有已经被证明只有更远的点的方向都完成时，解决方案是m。
由于我们逐渐减小m，每个维度所需的距离整体上迅速下降，尽管像所有算法一样，在更高的维度上下降得不太快。但是，如果仅在一个维度中的距离大于到目前为止最佳距离，那么在该方向上的其余点必然不能更好。
在时间复杂度方面，它似乎与更好的算法相当。但是，在数据结构的简单性方面，这个算法显然是胜出的。有很多其他属性使得这个算法成为一个严肃的竞争者。当你更新东西时，你可以用非常好的性能重新排序列表，因为你往往在排序已经排序或近乎排序的列表。你正在迭代数组。实际性能方面，大多数数据结构都很差。通常由于缓存和内存布局的原因，我们应该对此视而不见，但这却很重要。与您当前相关数据旁边的数据访问速度要快得多。如果我们已经知道要查找的点在列表中的位置，我们甚至可以更快地解决它（因为我们不必用二分搜索找到它）。还有其他允许重复利用上一次迭代信息的技巧。而额外的维度基本上是免费的（保存值不会更快地收敛，但这是因为球体中比半径相同的圆中具有更多随机分布的点）。
public class EuclideanNeighborSearch2D {
public static final int INVALID = -1;
static final Comparator<Point> xsort = new Comparator<Point>() {
@Override
public int compare(Point o1, Point o2) {
return Double.compare(o1.x, o2.x);
}
};
static final Comparator<Point> ysort = new Comparator<Point>() {
@Override
public int compare(Point o1, Point o2) {
return Double.compare(o1.y, o2.y);
}
};
ArrayList<Point> xaxis = new ArrayList<>();
ArrayList<Point> yaxis = new ArrayList<>();
boolean dirtySortX = false;
boolean dirtySortY = false;
public Point findNearest(float x, float y, float minDistance, float maxDistance) {
Point find = new Point(x,y);
sortXAxisList();
sortYAxisList();
double findingDistanceMaxSq = maxDistance * maxDistance;
double findingDistanceMinSq = minDistance * minDistance;
Point findingIndex = null;
int posx = Collections.binarySearch(xaxis, find, xsort);
int posy = Collections.binarySearch(yaxis, find, ysort);
if (posx < 0) posx = ~posx;
if (posy < 0) posy = ~posy;
int mask = 0b1111;
Point v;
double vx, vy;
int o;
int itr = 0;
while (mask != 0) {
if ((mask & (1 << (itr & 3))) == 0) {
itr++;
continue; //if that direction is no longer used.
}
switch (itr & 3) {
default:
case 0: //+x
o = posx + (itr++ >> 2);
if (o >= xaxis.size()) {
mask &= 0b1110;
continue;
}
v = xaxis.get(o);
vx = x - v.x;
vy = y - v.y;
vx *= vx;
vy *= vy;
if (vx > findingDistanceMaxSq) {
mask &= 0b1110;
continue;
}
break;
case 1: //+y
o = posy + (itr++ >> 2);
if (o >= yaxis.size()) {
mask &= 0b1101;
continue;
}
v = yaxis.get(o);
vx = x - v.x;
vy = y - v.y;
vx *= vx;
vy *= vy;
if (vy > findingDistanceMaxSq) {
mask &= 0b1101;
continue;
}
break;
case 2: //-x
o = posx + ~(itr++ >> 2);
if (o < 0) {
mask &= 0b1011;
continue;
}
v = xaxis.get(o);
vx = x - v.x;
vy = y - v.y;
vx *= vx;
vy *= vy;
if (vx > findingDistanceMaxSq) {
mask &= 0b1011;
continue;
}
break;
case 3: //-y
o = posy + ~(itr++ >> 2);
if (o < 0) {
mask = mask & 0b0111;
continue;
}
v = yaxis.get(o);
vx = x - v.x;
vy = y - v.y;
vx *= vx;
vy *= vy;
if (vy > findingDistanceMaxSq) {
mask = mask & 0b0111;
continue;
}
break;
}
double d = vx + vy;
if (d <= findingDistanceMinSq) continue;
if (d < findingDistanceMaxSq) {
findingDistanceMaxSq = d;
findingIndex = v;
}
}
return findingIndex;
}
private void sortXAxisList() {
if (!dirtySortX) return;
Collections.sort(xaxis, xsort);
dirtySortX = false;
}
private void sortYAxisList() {
if (!dirtySortY) return;
Collections.sort(yaxis,ysort);
dirtySortY = false;
}
/**
* Called if something should have invalidated the points for some reason.
* Such as being moved outside of this class or otherwise updated.
*/
public void update() {
dirtySortX = true;
dirtySortY = true;
}
/**
* Called to add a point to the sorted list without needing to resort the list.
* @param p Point to add.
*/
public final void add(Point p) {
sortXAxisList();
sortYAxisList();
int posx = Collections.binarySearch(xaxis, p, xsort);
int posy = Collections.binarySearch(yaxis, p, ysort);
if (posx < 0) posx = ~posx;
if (posy < 0) posy = ~posy;
xaxis.add(posx, p);
yaxis.add(posy, p);
}
/**
* Called to remove a point to the sorted list without needing to resort the list.
* @param p Point to add.
*/
public final void remove(Point p) {
sortXAxisList();
sortYAxisList();
int posx = Collections.binarySearch(xaxis, p, xsort);
int posy = Collections.binarySearch(yaxis, p, ysort);
if (posx < 0) posx = ~posx;
if (posy < 0) posy = ~posy;
xaxis.remove(posx);
yaxis.remove(posy);
}
}
更新：关于评论中的k-points问题。您会注意到很少有变化。唯一相关的是，如果发现点v小于当前m（findingDistanceMaxSq），则将该点添加到堆中，并将m的值设置为查找位置和第k个元素之间的欧几里得距离。算法的常规版本可以看作是k = 1的情况。当v被发现更接近时，我们搜索我们想要的1个元素并更新m等于唯一的（k = 1）元素。
请记住，我只在距离平方形式下进行距离比较，因为我只需要知道它是否更远，并且不浪费时钟周期在平方根函数上。
我知道有一个完美的数据结构用于存储k个元素在大小限制堆中。显然，数组插入对此不是最优解。但是，除了太多的java依赖api外，那个特定类别没有一个。但是，鉴于您的k很可能不那么大，因此您实际上不会注意到这一点。但是，它确实使点插入的时间复杂度存储在k时间中。还有诸如缓存元素与查找点之间距离之类的事情。
最后，也可能是最紧迫的是，我将用于测试代码的项目正在转换中，因此我还没有测试过这一点。但是，它肯定展示了您如何做到这一点：您存储到目前为止最好的k个结果，并使m等于到第k个最接近点的距离。——其他所有内容保持不变。
示例源代码。
public static double distanceSq(double x0, double y0, double x1, double y1) {
double dx = x1 - x0;
double dy = y1 - y0;
dx *= dx;
dy *= dy;
return dx + dy;
}
public Collection<Point> findNearest(int k, final float x, final float y, float minDistance, float maxDistance) {
sortXAxisList();
sortYAxisList();
double findingDistanceMaxSq = maxDistance * maxDistance;
double findingDistanceMinSq = minDistance * minDistance;
ArrayList<Point> kpointsShouldBeHeap = new ArrayList<>(k);
Comparator<Point> euclideanCompare = new Comparator<Point>() {
@Override
public int compare(Point o1, Point o2) {
return Double.compare(distanceSq(x, y, o1.x, o1.y), distanceSq(x, y, o2.x, o2.y));
}
};
Point find = new Point(x, y);
int posx = Collections.binarySearch(xaxis, find, xsort);
int posy = Collections.binarySearch(yaxis, find, ysort);
if (posx < 0) posx = ~posx;
if (posy < 0) posy = ~posy;
int mask = 0b1111;
Point v;
double vx, vy;
int o;
int itr = 0;
while (mask != 0) {
if ((mask & (1 << (itr & 3))) == 0) {
itr++;
continue; //if that direction is no longer used.
}
switch (itr & 3) {
default:
case 0: //+x
o = posx + (itr++ >> 2);
if (o >= xaxis.size()) {
mask &= 0b1110;
continue;
}
v = xaxis.get(o);
vx = x - v.x;
vy = y - v.y;
vx *= vx;
vy *= vy;
if (vx > findingDistanceMaxSq) {
mask &= 0b1110;
continue;
}
break;
case 1: //+y
o = posy + (itr++ >> 2);
if (o >= yaxis.size()) {
mask &= 0b1101;
continue;
}
v = yaxis.get(o);
vx = x - v.x;
vy = y - v.y;
vx *= vx;
vy *= vy;
if (vy > findingDistanceMaxSq) {
mask &= 0b1101;
continue;
}
break;
case 2: //-x
o = posx + ~(itr++ >> 2);
if (o < 0) {
mask &= 0b1011;
continue;
}
v = xaxis.get(o);
vx = x - v.x;
vy = y - v.y;
vx *= vx;
vy *= vy;
if (vx > findingDistanceMaxSq) {
mask &= 0b1011;
continue;
}
break;
case 3: //-y
o = posy + ~(itr++ >> 2);
if (o < 0) {
mask = mask & 0b0111;
continue;
}
v = yaxis.get(o);
vx = x - v.x;
vy = y - v.y;
vx *= vx;
vy *= vy;
if (vy > findingDistanceMaxSq) {
mask = mask & 0b0111;
continue;
}
break;
}
double d = vx + vy;
if (d <= findingDistanceMinSq) continue;
if (d < findingDistanceMaxSq) {
int insert = Collections.binarySearch(kpointsShouldBeHeap, v, euclideanCompare);
if (insert < 0) insert = ~insert;
kpointsShouldBeHeap.add(insert, v);
if (k < kpointsShouldBeHeap.size()) {
Point kthPoint = kpointsShouldBeHeap.get(k);
findingDistanceMaxSq = distanceSq(x, y, kthPoint.x, kthPoint.y);
}
}
}
//if (kpointsShouldBeHeap.size() > k) {
// kpointsShouldBeHeap.subList(0,k);
//}
return kpointsShouldBeHeap;
}
- Tatarize
9
1嗯，这是个有趣的想法。显然我们可以将k个项存储在堆（或优先队列）中，只需改变m的定义，使其成为堆中k个点的最远距离，而不是当前找到的最佳点的距离。因此，同样的技巧适用。我们只需保留到目前为止找到的最佳点的堆，m是我们找到的最差最佳点的距离。对于这个问题，堆将给我们带来最好的结果，因为我们需要不断地剔除k个项中的最差点。 - Tatarize
1算法的操作将是单次遍历，比您建议的knlog(n)排除方法要快得多。我们所承受的打击仅仅是因为这些项目的k值会使我们的m值收敛得更慢，因此我们最终考虑了我们本来不会考虑的点，这在根本上是正确的。这很像我们在更多维度中考虑这个算法时发生的情况。算法仍然完全相同，但截断那个方向的技巧不会那么快地生效。 - Tatarize
1@MyStackRunnethOver 看，这是概念证明代码。是的，这只是一个微不足道的更改。相同的核心算法，但我们跟踪前k个元素并保持m等于k个元素中最远的元素。相同的技巧适用。如果任何单个方向上的距离超过我们需要打败当前最佳元素列表的距离，我们就知道该点不能成为前k个元素列表中的元素，并且该方向上没有任何点可以成为前k个元素。 - Tatarize
1很大一部分的技巧就是尽可能地数学证明你所能证明的点，因为你已经有的点是无法更接近的。这也是为什么该算法最终会像所有其他NNS算法一样检查3D中的大多数点，即使是不太优雅的空间划分算法也是如此。许多其他相关算法也可以使用不同的度量或轻微修改来实现。 - Tatarize
1当我开始修改生产代码时，我会再次进行编辑。还有许多其他的速度技巧（尽管这比我用于我的目的的所有其他NNS都要快得多），目前还没有实现。例如，实际上情况是，当同一轴中的两个方向都被证明无果时，我们可以停止。如果在+x和-x、+y和-y或+z和-z中没有更多要检查的数字，则没有更多要检查的数字了。那就是整个数线。--而且由于每个列表都有自己的相同点的副本，如果我们可以标记该点为已检查，并且不重新处理。 - Tatarize
显示剩余4条评论
2
我会使用KD树来实现这个功能，假设点是随机分布的或者你有一种方法来保持树的平衡，时间复杂度为O(log(n))。
http://en.wikipedia.org/wiki/Kd-tree
KD树非常适合这种空间查询，甚至可以让你检索到距离查询点最近的k个邻居。
- Tom
1
据我理解，四叉树是用于二维的，但你可以计算类似于三维的东西。这将加快搜索速度，但如果实时计算索引，则需要更多时间。我建议先计算索引，然后存储它。在每次查找时，您需要找出所有外部四分之一，然后逐步查找命中点...就像剥橙子一样。随着四分之一变小，速度会大大增加。每件事都有一个权衡。
- CrazyDart
2
顺便说一句，如果你在同一个四边形中有大量的点，通常会在一个四边形中再做一个四边形...并继续嵌套到合理的分辨率。对于3D来说，这可能会花费很多...2D通常不会太糟糕。 - CrazyDart
13D结构被称为八叉树。 - moinudin
1
除非它们以适当的数据结构组织起来，否则唯一的方法就是线性搜索。
- ruslik
-1
从搜索的角度考虑，“最快”的方法是使用voxels。通过使用1:1点-体素映射，访问时间是恒定的且非常快速，只需将坐标移动到将点原点居中于体素原点（如果需要），然后向下舍入位置并使用该值访问体素数组。对于某些情况，这是一个不错的选择。正如我之前所解释的那样，当很难获得1:1映射时（点太多，体素分辨率太低，自由空间太大）八叉树更好。
- Santiago Pacheco
1
如果您的输入数据有一簇非常接近的点，那么您需要一个非常细密的网格。对我来说，这似乎不是一个明智的方法。“对于某些情况，这是一个好选择”太模糊了，根本没有说什么。您有任何链接到使用这种体素地图的论文或文章吗？您自己是否曾经使用过这种方法？ - Andreas Haferburg
-2
看这个.. 你也可以参考CLRS计算几何章节.. http://www.cs.ucsb.edu/~suri/cs235/ClosestPair.pdf
- Anwit
1
找到当前点最近的点是一个不同于找出数据集中最接近彼此的两个点的问题。 - Tatarize
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接