最近邻搜索问题的体素，K-D树，八叉树方法

Nearest Neighbour

1、暴力最近邻搜索

暴力最近邻搜索（Brute-force Nearest Neighbour Search）是最简单直观的最近邻计算方法，无需辅助数据结构。如果只搜索单个最近邻，则算法称为暴力单点最近邻搜索；如果搜索$k$点最近邻，则算法称为暴力$k$点近邻搜索，暴力近邻搜索准确率和召回率都为$100\%$，但计算时间过长。

应用暴力$k$点最近邻搜索方法查找点$\mathbf{x}$在点集 Y = { y 1 , ⋯   , y n } Y=\{\mathbf{y}_{1},\cdots,\mathbf{y}_{n}\} Y={y1​,⋯,yn​}中的$k$点最近邻的步骤：

---

1. 计算$\mathbf{x}$到点集$Y$中所有点的距离
2. 对第1步的结果进行排序，选择距离最短的$k$个点

---

2、栅格与体素方法

当查找某点在点集中的单个最近邻时，在二维情况下，可将原始空间划分为正方形栅格（Grid）；在三维情况下，可将原始空间划分为立方体体素（Voxel），并仅在查询点附近栅格或体素中查找。这种方法能很大限度缩小查找的范围。下面以说明三维体素方法所涉及的函数和具体步骤，二维栅格方法相关内容可以类推得到。

通过引入哈希函数，将相同体素中的点映射为相同哈希值，便于快速确定点集中位于查询点附近或体素中的点。哈希函数$\mathrm{hash}:{\mathbb{R}}^3\to \mathbb{R}$定义为

$$\mathrm{hash}(\mathbf{x})=((x_1\ast n_1)\mathrm{xor}(x_2\ast n_2)\mathrm{xor}(x_3\ast n_3))\mathrm{mod}N$$
其中$n_1,n_2,n_3$是三个大质数，$N$是大整数，可取$n_1=73856093,n_2=471943,n_3=83492791,N=10000000$

哈希函数需要与将位置坐标映射到体素序号的函数$\mathbf{Voxel}:{\mathbb{R}}^3\to {\mathbb{R}}^3$搭配使用，其定义为
$$\mathbf{Voxel}(\mathbf{x})=[\mathrm{round}(\frac{x_1}{r});\mathrm{round}(\frac{x_2}{r});\mathrm{round}(\frac{x_3}{r}){]}^T$$
其中$\mathrm{round}:\mathbb{R}\to \mathbb{Z}$为四舍五入函数，$r$表示分辨率，即体素的边在原始空间中的欧式长度

应用体素方法三维查找点$\mathbf{x}$在三维点集 Y = { y 1 , ⋯   , y n } Y=\{\mathbf{y}_{1},\cdots,\mathbf{y}_{n}\} Y={y1​,⋯,yn​}中的$k$点最近邻的步骤：

---

1. 将点集$Y$中每个点$\mathbf{y}$映射为$\mathbf{Voxel}(\mathbf{y})$，再计算哈希值$\mathrm{hash}(\mathbf{Voxel}(\mathbf{y}))$，以此建立哈希表
2. 记$V=\{[0,0,0{]}^T,[1,0,0{]}^T,[-1,0,0{]}^T,[0,1,0{]}^T,[0,-1,0{]}^T,[0,0,1{]}^T,[0,0,-1{]}^T\}$
3. 令$\mathbf{x}$的待查找最近邻集合$NN(\mathbf{x})=\varnothing$
4. 循环选取$V$中的元素$\mathbf{v}$，在每次循环中，记选取的元素为$\mathbf{v}$，计算哈希值$\mathrm{hash}(\mathbf{Voxel}(\mathbf{x})+\mathbf{v})$，通过查询哈希表将点云$Y$中具有相同哈希值的点加入$NN(\mathbf{x})$
5. 应用暴力单点最近邻搜索方法，查找点$\mathbf{x}$在$NN(\mathbf{x})$中的$k$点最近邻，并将该最近邻作为点$\mathbf{x}$在$Y$中的$k$点最近邻

---

3、K-D树

K-D树（K-Dimensional Tree）是借鉴了”对排序后的容器进行查找可以大幅节省时间“这一思想，

K-D树中规定：

• 每个非叶子节点都有左右两个分枝
• 叶子节点表示原始空间中的点
• 非叶子节点存储一个分割轴和一个分割阈值，来表达在节点中包含的点不完全一致的情况下，如何分割左分枝和右分枝，例如$y_i=thre$存储为轴$i$和阈值$thre$，表示左分枝取节点中在第$j$维上小于$thre$的元素，右分枝取节点中在第$j$维上大于等于$thre$的元素

根据点集 Y = { y 1 , ⋯   , y n } Y=\{\mathbf{y}_{1},\cdots,\mathbf{y}_{n}\} Y={y1​,⋯,yn​}构建K-D树的过程可以看成执行建树函数$\mathrm{BuildKDTree}(Y,\mathrm{nullptr})$，该函数采用递归方式定义，最终可通过$\mathrm{root}$访问整个K-D树，其定义为：

---

$\mathrm{BuildKDTree}(S,\ast \mathrm{root})$

1. 记当前节点$\mathrm{n}$为根节点$\mathrm{root}$

2. 若$|S|=1$

   • 记$\mathrm{n}$为叶子节点，退出

3. 若$|S|>1$

   • 记$\mathrm{n}$为非叶子节点，计算$S$中的点在各轴上的方差，选择方差最大的轴$i$作为分割轴，取平均数$m_i=\frac{1}{|S|}{\sum }_{\mathbf{y}\in S}{y}_i$作为分割阈值

   • 遍历$\mathbf{y}\in S$，若$y_i<m_i$，则插入$\mathrm{n}$的左子节点；若$y_i⩾m_i$，则插入$\mathrm{n}$的右子节点；若$S$中的多个点完全一致，则只保留一个点，并将节点$\mathrm{n}$记为叶子节点

   • 记全部插入到$\mathrm{n}$左子节点上的点组成的点集为$S_L$，递归调用$\mathrm{BuildKDTree}(S_L)$，构建K-D树；记全部插入到$\mathrm{n}$右子节点上的点组成的点集为$S_R$，递归调用$\mathrm{BuildKDTree}(S_R)$，构建K-D树

---

通过包含$6$个点的点集$\{[2,3{]}^T,[5,4{]}^T,[9,6{]}^T,[4,7{]}^T,[8,1{]}^T,[7,2{]}^T\}$构建的K-D树如下

从已构建完成的K-D树$\mathrm{T}$中查找点$\mathbf{x}$的$k$个最近邻的过程可看成首先定义$D_k=\varnothing$，$d_{max}=+\infty$，再执行$k$点最近邻搜索树函数$SearchKDTree\_k(\mathrm{T},\mathbf{x},k,D_k,d_{max})$，最终$D_k$成为$\mathbf{x}$的$k$点最近邻集合，该函数采用递归方式定义，其定义为：

---

$SearchKDTree\_k(\mathrm{T},\mathbf{x},k,\&D_k,\&d_{max})$

1. 初始化当前节点$\mathrm{n}$为$\mathrm{T}$的根节点

2. 若$\mathrm{n}$是叶子节点

   • 计算$\mathrm{n}$上的点$\mathbf{y}$与$\mathbf{x}$的距离$d(\mathbf{x},\mathbf{y})$，若$d(\mathbf{x},\mathbf{y})<d_{max}$，则 D k = D k ∪ { y } D_{k}=D_{k}\cup\{\mathbf{y}\} Dk​=Dk​∪{y}，若此时$|D_k|>k$，则删除$D_k$中到查询点距离最大的点，并重新计算$d_{max}=\{\begin{array}{ll}{\max }_{\mathbf{y}\in D_k}d(\mathbf{x},\mathbf{y})& |D_k|=k\\ +\infty & |D_k|<k\end{array}$
   • 退出

3. 若$\mathrm{n}$是非叶子节点

   • 确定$\mathbf{x}$在分割超平面的侧，也即计算$\mathbf{x}$落在的子树 n → C h i l d K D T r e e i , i ∈ { 1 , 2 } \mathrm{n{\to}ChildKDTree}_{i},i\in\{1,2\} n→ChildKDTreei​,i∈{1,2}。同时得到另一侧子树 n → C h i l d K D T r e e j , j = { 1 , 2 } \ i \mathrm{n{\to}ChildKDTree}_{j},j=\{1,2\}\backslash{i} n→ChildKDTreej​,j={1,2}\i

   • 递归调用$SearchKDTree\_k({\mathrm{n}\to \mathrm{ChildKDTree}}_i,\mathbf{x},k,D_k,d_{max})$，在$\mathrm{n}$的该侧子树中搜索$\mathbf{x}$的$k$点最近邻

   • 计算$\mathbf{x}$到分割超平面的垂直距离为$d_{split}$，若$d_{split}<d_{max}$，则递归调用$SearchKDTree\_k({\mathrm{n}\to \mathrm{ChildKDTree}}_j,\mathbf{x},k,D_k,d_{max})$，在$\mathrm{n}$的另一侧子树中搜索$\mathbf{x}$的$k$点最近邻；若$d_{split}⩾d_{max}$，则跳过$\mathrm{n}$的另一侧子树

   • 退出

---

查找点$\mathbf{x}=[6;5{]}^T$在上例点集中的两个最近邻，搜索流程如下图，结果为$\{[5;4{]}^T,[4;7{]}^T\}$

在上述的K-D树搜索过程中，最关键的部分是剪枝，即不继续在另一侧搜索$k$点最近邻，剪枝的条件是可以判定树形结构的另一侧不存在比现有结果更近的最近邻，判定条件如下
$$d_{split}⩾d_{max}$$
在该剪枝条件下，K-D树准确率和召回率方面都可以做到$100\%$，但可能会遇到去很远的分枝查找$k$点最近邻的情况，时间成本高。为解决该问题，可以添加一个比例因子$0<\alpha <1$，可将剪枝判定条件改为
$$d_{split}⩾\alpha d_{max}$$
可见剪枝条件被放宽，无法判定另一侧子树中是否存在比现有结果更近的最近邻。应用该剪枝条件，可使K-D树$k$点最近邻查找速度加快，但不再能保证找到严格的$k$点最近邻，这种做法称为近似最近邻（Approximate Nearest Neighbour，ANN）

4、八叉树

八叉树（Octo Tree，Octree)

根据点集 Y = { y 1 , ⋯   , y n } Y=\{\mathbf{y}_{1},\cdots,\mathbf{y}_{n}\} Y={y1​,⋯,yn​}构建八叉树树的过程可以看成执行建树函数$\mathrm{BuildOcTree}(Y,\mathbf{O},1,\mathrm{nullptr})$，该函数采用递归的方式定义，最终可通过$\mathrm{root}$访问整个八叉树，二者定义为：

---

$\mathrm{BuildOcTree}(S,\mathbf{box},depth,\ast \mathrm{root})$

1. 初始化当前节点$\mathrm{n}$为根节点$\mathrm{root}$

2. 若$depth=1$

   • 计算$Y$中的点在各轴上的最小值${\min }_i={\min }_{\mathbf{y}\in Y}{y}_i,i=1,2,3$和最大值${\max }_i={\max }_{\mathbf{y}\in Y}{y}_i,i=1,2,3$，这些最大值最小值组成表示包围盒的向量$\mathbf{box}=[{\min }_1,{\max }_1,{\min }_2,{\max }_2,{\min }_3,{\max }_3{]}^T$

3. 若$|S|=0$

   • 退出

4. 若$|S|=1$

   • 记$\mathrm{n}$为叶子节点，退出

5. 若$|S|>1$

   • 赋值$\{\begin{array}{l}{\min }_1=bo{x}_1,{\max }_1=bo{x}_2\\ {\min }_2=bo{x}_3,{\max }_2=bo{x}_4\\ {\min }_3=bo{x}_5,{\max }_3=bo{x}_6\end{array}$

   • 计算包围盒在各个轴上的中心位置$c_i=\frac{{\min }_i+{\max }_i}{2},i=1,2,3$

   • 计算$8$个子节点的表示包围盒的向量：

     ${\mathbf{box}}_1=[{\min }_1,c_1,{\min }_2,c_2,{\min }_3,c_3{]}^T$，${\mathbf{box}}_2=[c_1,{\max }_1,{\min }_1,c_2,{\min }_3,c_3{]}^T$

     ${\mathbf{box}}_3=[{\min }_1,c_1,c_2,{\max }_2,{\min }_3,c_3{]}^T$，${\mathbf{box}}_4=[c_1,{\max }_1,c_2,{\max }_2,{\min }_3,c_3{]}^T$

     ${\mathbf{box}}_5=[{\min }_1,c_1,{\min }_2,c_2,c_3,{\max }_3{]}^T$，${\mathbf{box}}_6=[c_1,{\max }_1,{\min }_1,c_2,c_3,{\max }_3{]}^T$

     ${\mathbf{box}}_7=[{\min }_1,c_1,c_2,{\max }_2,c_3,{\max }_3{]}^T$，${\mathbf{box}}_8=[c_1,{\max }_1,c_2,{\max }_2,c_3,{\max }_3{]}^T$

   • 遍历$\mathbf{y}\in S$，在循环内部遍历各子节点，根据各子节点的包围盒确定$\mathbf{y}$应当落在的子节点：若$\{\begin{array}{l}{bo{x}_j}_1⩽y_1⩽{bo{x}_j}_2\\ {bo{x}_j}_3⩽y_2⩽{bo{x}_j}_4\\ {bo{x}_j}_5⩽y_3⩽{bo{x}_j}_6\end{array}$，则插入第$j$个子节点

   • 记插入到各子节点的点组成的点集为$S_j,j=1,\cdots ,8$，对$8$个子节点分别递归调用$\mathrm{BuildOcTree}(S_j,{\mathbf{box}}_j,depth+1)$，构建八叉树

---

从已构建完成的八叉树$\mathrm{T}$中查找点$\mathbf{x}$的$k$个最近邻的过程可以看成看成首先定义$D_k=\varnothing$，$d_{max}=+\infty$，再执行$k$点最近邻搜索树函数$SearchOcTree\_k(\mathrm{T},\mathbf{x},k,D_k,d_{max})$，最终$D_k$成为$\mathbf{x}$的$k$点最近邻集合，该函数采用递归方式定义，其定义为：

---

$SearchOcTree\_k(\mathrm{T},\mathbf{x},k,\&D_k,\&d_{max})$

1. 初始化当前节点$\mathrm{n}$为$\mathrm{T}$的根节点，$\mathbf{box}$为表示$\mathrm{n}$的包围盒的向量

2. 若$\mathrm{n}$是叶子节点

   • 计算$\mathrm{n}$上的点$\mathbf{y}$与$\mathbf{x}$的距离$d(\mathbf{x},\mathbf{y})$，若$d(\mathbf{x},\mathbf{y})<d_{max}$，则 D k = D k ∪ { y } D_{k}=D_{k}\cup\{\mathbf{y}\} Dk​=Dk​∪{y}，若此时$|D_k|>k$，则删除$D_k$中到查询点距离最大的点，并重新计算$d_{max}=\{\begin{array}{ll}{\max }_{\mathbf{y}\in D_k}d(\mathbf{x},\mathbf{y})& |D_k|=k\\ +\infty & |D_k|<k\end{array}$
   • 退出

3. 若$\mathrm{n}$是非叶子节点

   • 若$\mathbf{x}$落在$\mathrm{n}$的包围盒的外面

     • 遍历$\mathrm{n}$的每个子节点，递归调用$SearchOcTree\_k(\mathrm{n}\to {\mathrm{ChildOcTree}}_i,\mathbf{x},k,D_k,d_{max})$
     • 退出

   • 若$\mathbf{x}$落在$\mathrm{n}$的包围盒的里面

     • 计算$\mathbf{x}$落在哪个子节点的包围盒中，也即确定$\mathbf{x}$落在的子树 n → C h i l d O c T r e e i , i ∈ { 1 , ⋯   , 8 } \mathrm{n}\mathrm{\to}\mathrm{ChildOcTree}_{i},i\in\{1,\cdots,8\} n→ChildOcTreei​,i∈{1,⋯,8}。同时得到其他子树 n . C h i l d O c T r e e j , j ∈ { 1 , ⋯   , 8 } \ i \mathrm{n}.\mathrm{ChildOcTree}_{j},j\in\{1,\cdots,8\}\backslash{i} n.ChildOcTreej​,j∈{1,⋯,8}\i

     • 递归调用$SearchOcTree\_k({\mathrm{n}\to \mathrm{ChildOcTree}}_i,\mathbf{x},k,D_k,d_{max})$，在$\mathrm{n}$的该子树中搜索$\mathbf{x}$的$k$点最近邻

     • 遍历$\mathrm{n}$的其他子节点，计算$\mathbf{x}$到其他各子节点包围盒各面距离的最大值 d b o x j = max ⁡ { ∣ x 1 − b o x j 1 ∣ , ∣ x 2 − b o x j 2 ∣ , ∣ x 2 − b o x j 3 ∣ , ∣ x 2 − b o x j 4 ∣ , ∣ x 3 − b o x j 5 ∣ , ∣ x 3 − b o x j 6 ∣ } {d_{box}}_{j}=\max\{|x_{1}-{box_{j}}_{1}|,|x_{2}-{box_{j}}_{2}|,|x_{2}-{box_{j}}_{3}|,|x_{2}-{box_{j}}_{4}|,|x_{3}-{box_{j}}_{5}|,|x_{3}-{box_{j}}_{6}|\} dbox​j​=max{∣x1​−boxj​1​∣,∣x2​−boxj​2​∣,∣x2​−boxj​3​∣,∣x2​−boxj​4​∣,∣x3​−boxj​5​∣,∣x3​−boxj​6​∣}，若${d_{box}}_j<d_{max}$，则递归调用$SearchOcTree\_k(\mathrm{n}\to \mathrm{ChildOcTree},\mathbf{x},k,D_k,d_{max})$，在$\mathrm{n}$的子树中搜索$\mathbf{x}$的$k$点最近邻；若${d_{box}}_j⩾d_{max}$，则跳过$\mathrm{n}$的子树${\mathrm{n}\to \mathrm{ChildOcTree}}_j$

     • 退出

---

上述八叉树搜索过程也涉及剪枝，即不继续在对应子树中搜索$k$点最近邻，剪枝判定条件如下
$${d_{box}}_j⩾d_{max}$$
在该剪枝条件下，八叉树准确率和召回率方面都可以做到$100\%$，但可能会遇到去很远的分枝查找$k$点最近邻的情况，时间成本高。为解决该问题，同样可以添加一个比例因子$0<\alpha <1$，可将剪枝判定条件改为
$${d_{box}}_j⩾\alpha d_{max}$$
剪枝条件被放宽，无法判定对应子树中是否存在比现有结果更近的最近邻。应用该剪枝条件，可使八叉树$k$点最近邻查找速度加快，但不再能保证找到严格的$k$点最近邻。

实验：不同单点最近邻搜索方法对比

现有两组三维点云，采用多线程版本的不同方法查找第一组点云中每个点在第二组点云中的最近邻所需的时间，准确率，召回率如下

	暴力搜索	体素方法（$r=0.1$）	K-D树（$\alpha =0.1$）	八叉树（$\alpha =0.1$）
时间（毫秒）	$345$	$1.13$	$1.61$	$4.77525$
准确率	$1$	$0.90$	$0.83$	$0.60$
召回率	$1$	$0.41$	$0.83$	$0.60$