根据相互距离对 2D/3D 点数组进行排序的启发式方法

Posted 2023-02-23

【中文标题】根据相互距离对 2D/3D 点数组进行排序的启发式方法

【英文标题】：Heuristics to sort array of 2D/3D points according their mutual distance

【发布时间】：2016-05-17 08:52:28

【问题描述】：

考虑 2D、3D、(4D...) 空间中的点数组（例如 unstructured mesh 的节点）。最初，数组中点的索引与其在空间中的位置无关。在简单的情况下，假设我已经知道一些最近邻连接图。

我想要一些启发式方法来增加空间中彼此接近的两个点具有相似索引的概率（在数组中会接近）。

我知道精确解决方案非常困难（可能类似于 Travelling salesman problem ），但我不需要精确解决方案，只需要增加概率的东西。

我对解决方案的想法：

一些天真的解决方案就像：

1. for each point "i" compute fitness E_i given by sum of distances in array (i.e. index-wise) from its spatial neighbors (i.e. space-wise)
   E_i = -Sum_k ( abs( index(i)-index(k) ) ) 
   where "k" are spatial nearest neighbors of "i" 
2. for pairs of points (i,j) which have low fitness (E_i,E_j) 
   try to swap them, 
   if fitness improves, accept

但具体实现及其性能优化不是很清楚。

其他不需要预先计算最近邻的解决方案将基于一些Locality-sensitive_hashing

我认为这可能是相当普遍的问题，并且可能存在好的解决方案，我不想重新发明***。

应用：

提高缓存局部性，考虑到内存访问通常是图遍历的瓶颈 它可以加速非结构化网格的插值，更具体地说是搜索靠近样本的节点（例如径向基函数的中心）。

【问题讨论】：

我什至不明白你在“幼稚的解决方案”中想说什么。如果两点接近或不接近，您的计算指标是什么？

一些指标，例如欧几里得。为什么？我使用哪个指标重要吗？最近的邻居也可以有多个定义，但一些自然定义将是距离最小的 N 个点。我不想指定这些细节，因为它会干扰问题的一般性。

gsamaras > 啊哈，混乱的根源是我搞砸了计算适应度的公式（更改为 k 和 j ）。现在我将 i 更正为 E_i = -Sum_k ( abs( index(i)-index(k) ) ) ...希望现在更清楚

你提到的LSH有什么问题？这似乎非常适合您的任务。

我怀疑@Regenschein，检查我的答案。

【参考方案1】：

我想说space filling curves (SPC) 是将空间接近度映射到线性排序的标准解决方案。最常见的是Hilbert-curves 和z-curves (Morton order)。

希尔伯特曲线具有最佳的邻近映射，但计算起来有些昂贵。 Z-ordering 仍然具有良好的邻近映射，但很容易计算。对于 z 排序，交错每个维度的位就足够了。假设整数值，如果您有一个 64 位 3D 点 (x,y,z)，则 z 值为 $x_0,y_0,z_0,x_1,y_1,z_1, ... x_63,y_63,z_63$，即 192位值由每个维度的第一位组成，然后是每个维度的第二位，依此类推。如果您的数组是根据该 z 值排序的，那么在空间中靠近的点通常在数组中也很靠近。

Here 是将 (merge) 值交织成 z 值的示例函数（nBitsPerValue 通常为 32 或 64）：

public static long[] mergeLong(final int nBitsPerValue, long[] src) 
    final int DIM = src.length;
    int intArrayLen = (src.length*nBitsPerValue+63) >>> 6;
    long[] trg = new long[intArrayLen];

    long maskSrc = 1L << (nBitsPerValue-1);
    long maskTrg = 0x8000000000000000L;
    int srcPos = 0;
    int trgPos = 0;
    for (int j = 0; j < nBitsPerValue*DIM; j++) 
        if ((src[srcPos] & maskSrc) != 0) 
            trg[trgPos] |= maskTrg;
         else 
            trg[trgPos] &= ~maskTrg;
        
        maskTrg >>>= 1;
        if (maskTrg == 0) 
            maskTrg = 0x8000000000000000L;
            trgPos++;
        
        if (++srcPos == DIM) 
            srcPos = 0;
            maskSrc >>>= 1;
        
    
    return trg;

您还可以交错浮点值的位（如果使用 IEEE 754 编码，因为它们通常在标准计算机中），但这会导致非欧几里德距离属性。您可能必须先转换负值，请参阅here，第 2.3 节。

编辑 两个回答 cmets 的问题：

1) 我了解如何为常规制作空间填充曲线 矩形网格。但是，如果我随机定位浮动 点，几个点可以映射到一个盒子里。该算法会起作用吗 在那种情况下？

有几种方法可以使用浮点 (FP) 值。最简单的方法是通过将它们乘以一个大常数来将它们转换为整数值。例如，将所有内容乘以 10^6 以保持 6 位精度。

另一种方法是使用 FP 值的位级表示将其转换为整数。这样做的好处是不会丢失精度，并且您不必确定乘法常数。缺点是欧几里得距离度量不再起作用。

它的工作原理如下：诀窍是浮点值没有无限精度，但仅限于 64 位。因此，它们会自动形成一个网格。与整数值的区别在于浮点值不形成二次网格，而是形成矩形网格，其中矩形随着与 (0,0) 的距离增加而变大。网格大小取决于给定点的可用精度。接近 (0,0)，精度 (=grid_size) 为 10^-28，接近 (1,1)，为 10^-16 见here。这个扭曲的网格仍然有邻近映射，但距离不再是欧几里得了。

这是进行转换的代码（Java，取自 here；在 C++ 中，您可以简单地将 float 转换为 int）：

public static long toSortableLong(double value) 
    long r = Double.doubleToRawLongBits(value);
    return (r >= 0) ? r : r ^ 0x7FFFFFFFFFFFFFFFL;


public static double toDouble(long value) 
    return Double.longBitsToDouble(value >= 0.0 ? value : value ^ 0x7FFFFFFFFFFFFFFFL);

这些转换保留了转换值的顺序，即对于每两个 FP 值，结果整数相对于 ,= 具有相同的顺序。非欧几里得行为是由编码在位串中的指数引起的。如上所述，这也在here，第 2.3 节中进行了讨论，但是代码的优化程度稍差。

2) 是否有一些算法如何对此类空间进行迭代更新 如果我的点在空间中移动，填充曲线？ （即没有重新排序 每次整个数组）

空间填充曲线强加了一个特定的排序，所以对于每一组点只有一个有效的排序。如果一个点被移动，它必须重新插入到由它的 z 值确定的新位置。

好消息是，微小的移动可能意味着一个点可能经常停留在阵列的同一“区域”中。所以如果你真的使用一个固定的数组，你只需要移动它的一小部分。

如果您有很多移动对象并且数组很麻烦，您可能需要查看“移动对象索引”（MX-CIF-quadtree 等）。我个人可以推荐我自己的PH-Tree。它是一种使用 z 曲线进行内部排序的按位基四叉树。它对于更新（和其他操作）非常有效。但是，我通常只推荐它用于较大的数据集，对于小型数据集，简单的四叉树通常就足够了。

【讨论】：

不错的答案。我也发布了一份，附有补充材料。

谢谢，我认为这最适合我的情况。我不确定两个方面：1）我了解如何为规则矩形网格制作空间填充曲线。但是，如果我随机定位浮点数，则可以将多个点映射到一个框中。这种算法在那种情况下会起作用吗？ 2）如果我的点在空间中移动，是否有一些算法如何迭代更新这种空间填充曲线？ （即无需每次都重新排序整个数组）

啊哈，我从 CGAL 补充中看到了这个页面，你的回答很好 doc.cgal.org/latest/Spatial_sorting/…

【参考方案2】：

你试图解决的问题有意义当且仅当，给定一个点 p 和它的 NN q，那么 q 的 NN 是 p 是真的。

这不是微不足道的，因为例如两个点可以代表景观中的位置，所以一个点可以在山上很高，所以从底部到山的成本更高反过来（从山到脚）。因此，请确保您检查这不是您的情况。

---

既然 TilmannZ 已经提出了解决方案，我想强调一下你提到的LSH。我会不选择那个，因为您的积分位于真正的低维空间，它甚至不是 100，那么为什么要使用 LSH？

在这种情况下，我会选择CGAL 的算法，例如2D NNS，甚至是简单的kd-tree。如果速度很关键，但空间不是，那么为什么不选择quadtree（3D 八叉树）？我已经构建了一个，它在 8GB 内存中不会超过 10 个维度。

但是，如果您觉得您的数据将来可能属于更高维度的空间，那么我建议您使用：

LSH 来自 Andoni，真的很酷。 FLANN，它提供了另一种方法。 kd-GeRaF，是我自己开发的。

【讨论】：

谢谢，TilmannZ 的回答比较中肯，不过你的cmets 也很有用。 CGAL 库可能非常有用，但我尽量保持我的代码库小而简单，没有太多依赖。但我可能会研究 CGAL 代码来复制一些算法。

欢迎您@ProkopHapala。是的！嗯，这有点难，因为cgal 一开始有点脾气暴躁！也许自己实现算法会更好，但你说你不想重新发明***，但我看到你在这里提到的权衡。祝你好运！顺便说一句，好问题！