点云光线投射(Point Cloud Ray Casting)方案简述与BVH树构树思路
需求
最近尝试做一个点云上色的东西，发现目前互联网上对于点云光线碰撞的内容似乎寥寥无几，遂写下该文。本文将以光线碰撞切入题目，简要阐述各个点云树形的优缺点，最后逐步阐述BVH树的构建。
给点云上色，其本质为光线碰撞(Ray casting)的一种应用，具体的说：给定一个光线（即光源坐标+方向向量），求哪些点云最先与光线碰撞，然后对碰撞的点云进行操作。在本文中，该操作为给点云上色。显然如果检测光线和每个点云是否发生了碰撞，那么它的搜索平均时间复杂度为O(n)，那对于给定的上亿个点云搜索的效率是难以接受的，所以我们需要找到一种树的存储方式，将搜索平均时间复杂度降到O(logN)。
在开始阐述算法前，我们需要明确点云分辨率的概念，它是指单个点云对应现实中扫描的精度，由于雷达位姿估计不一致等原因，我们的扫描所得到的精度可能只有1cm或者2cm，我们期望在点云的扫描精度范围内完成上色操作。在我们的需求中，如果只检测点是否在光线所在的直线上那么我们大概率搜索不到点云，所以我们需要以点云为中心，建立一个边与坐标轴平行且边的长度与点云分辨率相同的正方体，我们称之为体素(Voxels)。那么此时我们的检测光线是否与点云发生碰撞的需求就变成了检测光线是否与体素(也就是一个边与坐标轴平行的正方体)发生了碰撞，又可能是由于对某一目标多次不同角度扫描的缘故，某一立方厘米的空间内存在很多个点，则此时我们对该空间上的点进行批量上色。
回到点云碰撞需求上，为了加速点云碰撞，我们期望先检测光线是否与一个大的包围盒是否发生了碰撞，如果存在碰撞再去搜索里面的点云，否则则直接跳过该包围里里面的点。成组地接受或者拒绝碰撞检测是必要的，同时当包围盒小到一定程度(如达到最小包围盒限制或者达到了包围盒点云个数的下线)时直接给整个包围盒的点云上色或者继续搜索内部点云的包围盒，顺着这个思路，我们需要思考该如何选择构建一个"大包围盒"的树呢？
方案
目前存储点云的树形主要有以空间划分的Oct-tree、KD-tree和以对象为划分的BVH-tree（层次包围盒）、Ball-tree(层次包围球)，其中空间划分的应用方法在计算机图形学光线投射的领域中较为主流。在这里我们先逐一介绍下每一种树形的构树内容。
首先所有树在初始化时都需要搜索点云在X、Y、Z轴方向上的最小值和最大值，然后将其作为对角线坐标构建一个最大的包围盒(这种包围盒也被称为AABB盒，包围盒的边与坐标轴平行)。关于AABB盒的碰撞检测和优化问题，Another View on the Classic Ray-AABB Intersection Algorithm for BVH Traversal 讲述得较为详细。
- Oct-tree，八叉树的原理很简单，在构树时对每一个大的包围盒于沿着中点“砍两刀”并生成了八个子包围盒，直到砍到某一个深度或者体素为止，它作为最为广泛使用的点云存储树形，它的优点在于数据存储方便，增减点云的操作只需要搜索到对应的大包围盒然后进行修改即可，可以很好地描述移动的物体(或者点云)，而缺点则是由于它是沿空间均分点云，而点云数据并非均匀地分布在空间中，所以构建出来的树并非平衡树，最差的搜索效率可能达到O(n)。八叉树的搜索流程也很简单，先对最大的包围盒进行碰撞检测，如果检测到了碰撞，则对其八个子节点逐一碰撞，并选取与光线最先碰撞的节点继续展开树碰撞，它的搜索流程是一个深度优先的搜索流程。特别的，点云库PCL和Ray casting最相关的函数getIntersectedVoxelIndices()使用的即为八叉树的方案。
- Kd-tree，KD树相比八叉树具有更强的空间适应性，它在构树时只切分一次，沿哪个轴切分、在轴的哪个位置切分则需要根据实际情况考量，如果需要更快的搜索效率则每次沿着方差最大的轴按中位数或者平均数切分，如果需要更好的删减效率则每次按X、Y、Z的顺序按中位数或者平均数切分，这里多一句话，按照中位数、平均数或者等数量的划分方法有时候也会生成一个很差的树形，如果选择表面积启发法贪婪地搜索最优切分位置也可以生成性能更加优越的树。KD树的优点在于搜索效率高，是解决K-最近邻(K-NearestNeighbor)最朴素的方法。KD树的搜索流程与八叉树相近，也是深度优先的搜索方式。
- BVH树，包围盒层次(Bounding Volume Hierarchies)的空间适应性比KD树更强，它在构树时的切分方式KD树一样，只切分一次，但是在切分完成后会重新针对切分的点云生成两个包围盒。KD树在碰撞检测时是按照轴划分生成的包围盒，而BVH树在碰撞检测时生成的包围盒更小，可以更快地拒绝父节点相交但两个子节点未相交的情况。
- Ball-Tree，将BVH树中的包围盒修改为包围球，对于点云而言也许用球来表述包围效果更好(见后记)。
从点云增加、删减等维护性操作来说，维护性排序为: Oct-tree > KD-tree > BVH
从搜索效率来说，BVH > KD-tree > Oct-tree
从建树效率来说，Oct-tree > KD-tree(中点划分) > BVH(中点划分)
在这里我们为什么说搜索效率是BVH>KD>OCT呢？可以想象，一群点云构成的包围盒中，BVH的包围盒最为紧凑，KD树次之，OCT最差，而生成的包围盒如果不够紧凑就很容易出现光线与父节点相交而子节点不相交的情况，在该种情况下我们可能需要遍历更多次包围盒，显然这样是得不偿失的行为。同时，越紧凑的包围盒意味着增减点云需要修改的东西越多，所以维护性会变差。
构树
好了，再次回到我们的需求上，我们现在开始明确下光线相交的开销来源，一是树形的选择，对于上亿点云我们是不考虑实时性(因为如果实时构树得用GPU计算了)，即我们期望的是光线碰撞的速度最快，而不优先考虑树的维护性和建树效率，则选择BVH方案，通过提前计算好包围盒信息和较少无效检测的次数是三种方案中的最优选择；二是相交检测时的开销，根据局部性原理我们不希望通过节点信息离散的存储在内存中，所以我们需要先实现一个内存池；我们也不希望通过递归搜索树，所以我们还需要线性化树。
我们先设计一个简单的内存池Memory.h，内存池的作用是向系统申请一大块连续的内存，然后自行管理内存信息，在这里我们参照了PBRT的内存池设计思路。其中，MemoryArena类，类的作用是申请一定颗粒度的内存(由Blocksize指定)并由指针currentBlock指向，整数currentBlockPos标识分配位置，当内存池分待分配的内存不足时，则将该段信息放入列表usedBlocks中，然后为currentBlock再申请一定颗粒度的内存，并将currentBlockPos设置为0。MemoryArena的作用是在当申请的内存大小不明确时，按照由Blocksize指定的颗粒度申请一大块内存，其主要用于递归构树时候的内存分配，使用了内存池后构树时就不用因为频繁向系统申请内存了。我们还需要直接向内存申请一大块连续内存用于连续存储线性化后的节点等信息，函数AllocAligned使用了std::aligned_alloc，其中实参 alignment=64 是指代CPU架构的缓存行大小为64(详情见CPP REFERENCE)。
其次我们需要完成基础设施的各种类，如点类、光线类(Baseshape.h)、包围盒类(Bounds.h)并提供如计算表面积、判断相交、计算相交等方法，其中包围盒的相交测试IntetsectP和IntetsectPB函数的形参 scale 默认值为 1. + 2*gamma(3)，作用为增大检测阈值避免因为浮点计算误差而错过相交检测。
最后我们选择Glog库作为日志库输出测试信息(Glog真的好用),并使用valgrind检测内存消耗与内存泄露情况(虽然后面的调试信息没用到这个)；使用chrono计时构树以及碰撞所需的时间(据说chrono的计时精度为纳秒级，也就是\(10^{-3}\)毫秒，在后文中，光线测试碰撞的时间几乎都为纳秒级别，我也不懂那个时间到底准不准....)。
下面我们将实现一个BVH树
BVH树的构树分为三个阶段：首先，初始化空间，为所有点云添加下标信息存储于pointInfo中，便于后续排序。其次使用选择的SplitMethod算法构建树，结果是构建得到的节点都有指向两个孩子的指针、节点的包围盒信息、划分轴、叶子存储的第一个点云的下标偏移量以及点云个数(若该节点不是叶子则两个数值均为0)。最后，将树转化为线性树一遍后续更快地进行光线碰撞。这里说一下为什么叶子对应的点云信息要选择存储下标而不用vector<point3f *>存储呢？因为存储vector<point3f *>的内存占用太大了，并且会使得单个节点的内存占用空间超过32字节，为了对齐内存而扩充到64字节这种消耗是难以接受的，而存储下标的方式，每个叶子用于存储点云所用的内存为4字节的下标偏移+4字节的点云个数=8个字节，无论叶子存储多少点云数量都是这个稳定的内存消耗，更省内存而且是顺序存储，读取速度也更快，性价比满满好吗！
我们首先确定建立BVH树时所需要的节点信息
struct BVHBuildNode
{
Bounds3<float> bounds;//包围盒，使用对角线的两个坐标pMin和pMax描述
BVHBuildNode *children[2];//子节点
int splitAxis, firstPrimOffset, nPrimitives;//分割轴，如果该节点为非叶子，则nPrimitives、firstPrimOffset等于0。
};
/*如果firstPrimOffset == 245，nPrimitives == 30，则代表这个叶子存放了30个点云，点云对应的下标存储在orderdata[245],orderdata[246]...orderdata[274]中*/
BVH在建树完成后只需要维护以下三个变量
LinearBVHNode *nodes; //根节点
std::vector<Point3<float>> &_pointcloud;//所有点云数据
std::vector<uint> orderdata;//存放叶子的点云下标
构树第一阶段由构造函数前半部分负责，这里不过多介绍。从构树第二阶段开始，递归地调用recursiveBuild建立树，递归流程中，首先是建立当前点云所在的包围盒，然后判断该包围盒数据和点云数据是否符合叶子标准，如果符合则返回节点，如果不符合则继续按包围盒最大方向轴按特定方法划分，划分后继续递归迭代。
- 在分割轴的选择上，本文选取的是包围盒的最大值方向作为分割方向，传统KD树是选择方差最大的方向作为分割轴，然而在本文中经测试沿方差最大方向建树耗时更长(对于九千万个点云数据则构树需要多10s的时间)但树的搜索速度相差无几，其中原因见后记。
- 在分割轴位置的选择上，本文提供了三种选择方式：包围盒中值划分、等数量划分、SAH划分。在构树时间上，包围盒中值划分<SAH划分<等数量划分; 在性能上，等数量划分 < 中值划分 < SAH划分。也就是说等数量划分最差，中值划分在构树和性能上适中，SAH在性能上最好。
构树的第三阶段是将指针存放的树线性化，我们先构造线性化树的节点结构体。
struct LinearBVHNode
{
Bounds3<float> bounds;
union
{
int primitivesOffset; // 如果是叶子则表示该叶子第一个点云在orderdata的下标位置
int secondChildOffset; // 如果是节点则指代右子节点位置
};
uint16_t nPrimitives; // 如果不是0则代表该节点为叶子，否则为节点
uint8_t axis; // 分割轴
uint8_t pad[1]; //标识是由最小包围盒围成的点云盒还是由最小点云数目包围成的点云盒
};
我们期望依照深度优先的顺序压平树，也就是说每个节点的第一颗子节点位于该节点位置的右边的第一个位置，而第二个子节点的位置则由secondChildOffset描述。这种结构的树在搜索时也有更快的搜索速度(万能的局部性原理)。
那么此时写压平树的递归函数的思路则显而易见了，我们在递归构树同时统计总叶子节点数量total，在完成递归构树后，我们直接调用nodes = AllocAligned
uint flattenBVHTree(BVHBuildNode *node, uint *offset)
{
LinearBVHNode *linearNode = &nodes[*offset];
linearNode->bounds = node->bounds;
int myOffset = (*offset)++;
if (node->nPrimitives > 0)
{
linearNode->primitivesOffset = node->firstPrimOffset;
linearNode->nPrimitives = node->nPrimitives;
}
else
{
linearNode->axis = node->splitAxis;
linearNode->nPrimitives = 0;
flattenBVHTree(node->children[0], offset);
linearNode->secondChildOffset = flattenBVHTree(node->children[1], offset);
}
return myOffset;
}
终于到了最后，我们回到最开始的需求，如何实现光线碰撞，其实这本质上就是一个树的深度优先搜索问题，线性化后的树遍历非常简单，只需要少量维护树位置的数据即可，我们需要使用一个栈来描述下一个扫描的位置。由于在检测光线是否与包围盒发生碰撞时总是需要光线方向向量的倒数，我们提前计算好该数据invDir并作为参数传入bounds.IntersectPD函数中，bounds->IntersectPD函数的作用是返回光线和包围盒的最短相交长度，如果返回值为0则表示并不相交。我们通过判断currentNodeIndex指代的节点的两个子节点是否和包围盒相交来决定栈的弹入。如果两个子节点相交则按距离近远的顺序将子节点弹入栈中，如果只有某一个节点相交则只将该节点弹入栈中，如果都不相交则不对栈做操作，最后从取出栈顶作为currentNodeIndex的内容，循环至栈为空或者达到最大深度位置。
/*
该函数用于计算光线和点云相交的叶子，并将叶子以下标的形式放入vector中，depth = 1 表示当检测到第一个碰撞时立刻结束搜索，即为本文光线最先碰撞的需求,scale表示相交检测，默认值不修改。
*/
void IntersectP(const Ray &ray, std::vector<uint> &ret_nodes, uint depth = 1, double scale = 1 + 2 * gamma(3))
{
if (!nodes)
return;
Point3d invDir(1.f / ray.d.x, 1.f / ray.d.y, 1.f / ray.d.z);
int dirIsNeg[3] = {invDir.x < 0, invDir.y < 0, invDir.z < 0};
uint nodesToVisit[1024] = {0};
uint toVisitOffset = 1, currentNodeIndex = 0;
uint leftnode_index = 1, rightnode_index = nodes->secondChildOffset;
uint currentdepth = 0;
while (toVisitOffset != 0)
{
LinearBVHNode *node = &nodes[currentNodeIndex];
if (node->nPrimitives > 0)
{
ret_nodes.push_back(currentNodeIndex);
currentdepth++;
if (currentdepth == depth)
{
break;
}
}
else
{
leftnode_index = currentNodeIndex + 1;
rightnode_index = node->secondChildOffset;
const LinearBVHNode *leftnode = &nodes[leftnode_index];
const LinearBVHNode *rightnode = &nodes[rightnode_index];
double leftInsect_ = leftnode->bounds.IntersectPD(ray, invDir, dirIsNeg, scale);
double rightInsect_ = rightnode->bounds.IntersectPD(ray, invDir, dirIsNeg, scale);
bool leftInsect = leftInsect_ > 0 ? 1 : 0;
bool rightInsect = rightInsect_ > 0 ? 1 : 0;
if (leftInsect && !rightInsect)
{
nodesToVisit[++toVisitOffset] = leftnode_index;
}
else if (!leftInsect && rightInsect)
{
nodesToVisit[++toVisitOffset] = rightnode_index;
}
else if (leftInsect && rightInsect)
{
if (leftInsect_ < rightInsect_)
{
nodesToVisit[++toVisitOffset] = rightnode_index;
nodesToVisit[++toVisitOffset] = leftnode_index;
}
else
{
nodesToVisit[++toVisitOffset] = leftnode_index;
nodesToVisit[++toVisitOffset] = rightnode_index;
}
}
}
currentNodeIndex = nodesToVisit[toVisitOffset--];
}
}
测试
一、斯坦福兔子
首先当然是对著名的斯坦福兔子进行实验啦！
BVH构树传参点云分辨率为0.1，最小包围盒边长为0.2，分割方法为中点分割，内存池颗粒度为128Mb。输入点云为95947个点云的兔子点云数据。如图完成的构树如下：
此时，调试所得信息如下，构树时间为9毫秒，BVH树运行时维护53937个节点，并需要消耗1.65MB的内存，而构树过程中最高消耗了128.00MB的内存(构树时选择内存的颗粒度太大了，其实没必要嗷)。在光线碰撞问题上，搜索光线方向上所有的点云共计碰撞了59次包围盒，如果当碰撞到最近距离时立刻返回则耗时0.003毫秒。
15:55:21.236143 155514 main.cpp:101] Point sizes : 35947
15:55:21.245301 155514 BVH_1.h:96] BVH created with 53937 nodes for 35947 points (1.65 MB), arena allocated 128.00 MB
15:55:21.245569 155514 main.cpp:107] The time to build tree :9.00269 ms
15:55:21.245608 155514 main.cpp:123] Search block: 59
15:55:21.245620 155514 main.cpp:124] The time to search tree :0.003366 ms
二、车库大规模点云
BVH构树传参点云分辨率为0.01，最小包围盒边长为0.02，内存池颗粒度为128Mb。输入点云为93226744个点云某一车库扫描所得的点云数据。
下图为全部点云信息
下图描述了光线与哪些包围盒碰撞。
本文提供的BVH树分割方法共有3种：中点分割、等量分割和SAH分割，以下为各种分割方法的调试信息。其中IntersectPB搜索了光线沿途所有与之碰撞的点云所用的时间，而IntersectP则搜索了最先与光线碰撞的点云所用的时间，碰撞的调试信息时间单位均为ms。由于搜索量级过小，难以从时间上看出显著的差异，比较不同分割方法的搜索性能可以比较 搜索的包围盒个数Search Block或者比较建立的节点个数即可。
分割方法为等数量划分
10:59:13.234164 360541 BVH_1.h:110] BVH created with 46255157 nodes for 93226744 points (1411.60 MB), arena allocated 3072.00 MB
10:59:13.370687 360541 t1.cpp:101] The time to build tree :37974.2
10:59:13.468467 360541 t1.cpp:116] Block found: 558
10:59:13.468501 360541 t1.cpp:117] The time of IntersectPB :0.044483
10:59:15.738780 360541 t1.cpp:130] The time of IntersectP :0.005764
分割方法为中点划分
10:57:24.640923 360107 BVH_1.h:110] BVH created with 32964309 nodes for 93226744 points (1005.99 MB), arena allocated 2048.00 MB
10:57:24.743727 360107 t1.cpp:101] The time to build tree :19096.3
10:57:24.743813 360107 t1.cpp:116] Block found: 301
10:57:24.743822 360107 t1.cpp:117] The time of IntersectPB :0.026437
10:57:26.011411 360107 t1.cpp:130] The time of IntersectP :0.007619
分割方法为SAH
10:54:04.558768 359100 BVH_1.h:110] BVH created with 32887887 nodes for 93226744 points (1003.66 MB), arena allocated 2048.00 MB
10:54:04.669757 359100 t1.cpp:101] The time to build tree :42497.5
10:54:04.669829 359100 t1.cpp:116] Block found: 260
10:54:04.669837 359100 t1.cpp:117] The time of IntersectPB :0.018552
10:54:05.431996 359100 t1.cpp:130] The time of IntersectP :0.004501
可以看到换了不同的分割方法，构树时间却有很大的差异(废话，改的参数都是构树时才会用到的参数)，从总碰撞的包围盒来看，按等数量划分总是最差的选择，仔细一想原因也很简单，如果点云分布不均匀，那按等数量划分就是逮着点云密集的地方中间砍一刀，而我们期望的其实是在点云不密集的地方进行切分，也就是使用SAH算法进行划分才更加合理。
结论
基于BVH树构建的点云树具有比oct、kd树更加优越的光线碰撞性能，但是其代价为单个节点占用了较多的内存(32字节)，在划分轴和分割位置的选择上，默认的中点划分、包围盒轴最大方向即可满足大部分情况。笔记本平台i5-8265u 16g-RAM，针对近亿点云的数据，构树时间需要23秒左右，而碰撞时间为纳秒级别，由于此时的计时精度已经达到chrono库的最高精度，所以实际结果可能有一定误差。在后续实验中，本文比较了单条光线搜索沿途所有点云的情况，最终结果表明，相比于中点划分法，基于SAH的光线碰撞可以实现接近13%的碰撞速度提升(粗略的估计，从碰撞盒计算得来的不严谨的数据)。
后记
本文代码的GitHub地址为： https://github.com/EuNll0/PointCloud_Tree/tree/main 本文在写的过程中也包含了一些自己的想法，然而我也不总是能猜对，所以如有纰漏，请私信指正。
更好的树
前几天睡觉的时候就在想，用包围盒描述节点的话，那我一个节点的内存占用为包围盒信息(2*12=24)+4+3 = 31个字节。为了内存对齐则该节点的上限为32个字节，除去必要的信息，可以把1个字节描述的分割坐标轴位置换成用2个比特存储，换句话说BVH节点只有6+8=14个比特存储额外信息，这并不方便添加如树深度等信息，那有没有更加节省内存的节点描述方式呢？当然有，那就是用球树描述节点信息，一个球只需要由球心和半径描述，也就是说节点的内存占用为16+4+3=23个字节，则距离32个字节还有很大一部分空间可以用于存储额外信息。同时，判断包围盒是否与光线相交需要三次计算判断才能真正确认是否相交并得到最小相交长度，而判断球是否与光线相交只需要二次计算判断即可得到结果，而且不用特别处理光线方向的分量为0的情况(在本文代码中没有特别处理该情况，也就是会出现除以0的情况)，同时KNN计算也更加快速。
总结一下，如果针对点云数据，相比于BVH树，使用球树有三个优点
- 更多空余的内存
- 更快的相交判断(但因为球会重叠，也可能要更多次相交判断)
- 更快的KNN计算
当然，这些只是理论上的猜想，既然球树这么多优点，那为什么我不写球树呢？一方面是目前PCL库做法都是用大包围盒包围点云(算是PCL误导了我(bushi))，另一方面则是因为我写完BVH树才意识到球树的性能更加优越，所以还得再改改代码(其实主要是懒)，仔细想想构建一个球树会遇到哪些问题呢？联立球(简单问题)、离群点处理(这可以用PCA搞定)....看起来这些问题也不难。
对了，如果使用KD树，那一个节点可以用联合体存储,8个字节足以存储节点信息，相比BVH节点也能节省很多内存。据PBRT所述，在三角形碰撞实验中，相比16字节的节点，8字节的节点可以在碰撞检测上提速近20%(因为CPU的缓存行可以存放更多节点)，不过相对的，在碰撞检测的时候需要生成包围盒并且可能会多一些无效碰撞，这是一种取舍，读取得更快但碰撞时需要计算额外的数据。那我为什么不选择该方案呢？因为使用位操作比较麻烦(bushi)，是因为我还需要存储在叶子中存储一些其他的信息(x)，BVH树的思路更清晰更好写代码(√)。在参考文献中有一篇论文也描述了如何更快地构建KD树，大致思路就是先对XYZ三个方向进行排序，然后再构树，当然这个思路也可以应用在BVH树中。
更好的分割方法
- 分隔位置：
我们使用中点分割又或者使用等数量分割生成的树不一定是一颗最好的树，为什么这么说呢？因为如果点云在包围框中间部分特别密集，那么按中点划分产生的两个节点在碰撞或者搜索时性能会很差(明明他们近在咫尺，只是刚好站在了中线两边，为什么要把他们分开！)我们期望的分割位置其实是按聚类来分割，也就是聚类后在稀疏部分划分生成两个子节点，这就乐了，我们构树的目的之一就是为了聚类的目的，结果要求聚类了才能构树，这合理吗？解决这个问题的方法有两个，一是多次构树多次聚类多次优化，这似乎在OBB形状的树中中可以用到，但是做光线碰撞就很明显纯纯不可行，因为数据量太大的时候时间消耗太大了，这是得不偿失的行为。这时候第二种方法就很有意思了，那就是动态划分的SAH算法。它根据包围框表面积与点云个数构建一个代价函数，对包围盒表面积和包围盒中的点云个数加权求和，然后取代价最小的位置作为分割位置，点云的SAH算法贼简单，因为包围盒重叠的问题完全可以不用考虑(数量级原因)，这里再提一句思路，对分割轴上的点云建立桶，然后再去执行SAH算法。 - 分割轴：
在KD树中，我们选取的分割轴为方差最大的轴，这在思路上是合理的，因为在方差反映了数据离散程度，同时也在一定程度上反映了数量级，例如在地下车库中，我们的Z轴范围可能就在0-5米间，而在X、Y轴上却很大，此时X、Y方向上的方差也很大，此时也符合我们的希望——尽量先在X、Y轴上划分从而使得单个包围盒能够包围足够多的数据，但是计算方差则需要遍历两次点云数据，构树时间会边长，而我用最大包围盒边长的方向近似也能达到很好的构树效果，综合了性能和准确性，我们在程序代码中我并没有选择方差最大方向作为分割方向。而从聚类的角度来说，其实我们更期望车库的天花板直接聚成一类，假设天花板与O-X-Y平行，也就是说我们希望一开始可能就是沿着Z轴切割，那又该如何调参确定选取的分割轴呢？那就是在构树的时候对X、Y、Z三个方向上都计算一遍最佳的分割位置代价然后选择最好的轴和分隔位置，这样也能构造好一个期望的树类型。
为什么在计算机图形学中说KD树的光线相交比BVH树更快
在计算机图形学中，每个图元(如三角形)在BVH树中只能分配到一个叶子，而在KD树中同一个图元可以分配到不同的叶子(因为如果某些图元的包围框横跨了分割轴则将该图元分配到两个叶子中)。在光线碰撞的需求中，我们期望的是当碰撞到一个图元时就立即返回，如果光线在分割轴附近相交，那么在BVH树中就可能出现要碰撞的图元不在最近的包围盒中，而KD树因为将一个图元分配到两个子树中就可以通过相交快速得到结果了。如下图所示，蓝色为光线，KD树的光线碰撞会在检索左边包围框时就搜索到对应三角形(该三角形在左右盒子均可以被搜索到)，而BVH树搜索左盒子时未检测到碰撞，从而需要搜索整个右盒子，这里就存在多余的开销了。而在点云数据中，如果我们为了纯粹地光线加速也可以这样实现KD树，但此时KD树的扩展性能将会下降(例如最近邻搜索能力)。
为什么要选择AABB盒：
我们上述用到的包围盒均为AABB盒，虽然它做光线碰撞的代码很简单但是却丢失了旋转的特性，也许也可以使用OBB盒(也就是盒子的边不一定与轴平行)可以更好地描述点云对象，那么该如何检测碰撞呢？我个人的思路是使用投影变换的方法，也就是在光源处建立一个以光线方向为Z轴，然后在光线的法平面上再任意选取两个正交的向量最终组合成正交矩阵，然后用这个矩阵左乘包围盒的顶点(也就是对包围盒作投影变换)，最后判断光心在不在投影形成的四边形内即可。仔细一想，用OBB盒做光线碰撞也是吃饱了撑的人才会做的事，用OBB盒做聚类或者建模好像可行，以后试试。
参考文献
PBRT
How to create awesome accelerators: The Surface Area Heuristic
Another View on the Classic Ray-AABB Intersection Algorithm for BVH Traversal
Building a Balanced k-d Tree in O(kn log n) Time
RTFC(雾) STFW(逃)