A new point containment test algorithm based on preprocessing and determining triangles

摘要

在本文中，我们重新讨论了多面体问题中的点。在回顾了之前的工作之后，我们进一步深入了解了这个问题。然后我们声称，对于给定的测试点和三维多面体，可以找到一个确定三角形，这足以确定该点是在多面体内部还是外部。

这项工作可以被认为是霍恩工作的延伸和实现，霍恩的工作启发我们提出了一个定理来获得确定三角形。在这个定理的基础上，提出、实现和测试算法。结果表明，尽管我们的代码具有与常用的基于八叉树的射线交叉方法相同的渐近时间效率，但在实践中，它通常快几倍，有时快十倍以上，而其他成本（如预处理时间和内存需求）保持不变。

本文提出的想法简单而普遍。因此，它们自然扩展到多材料模型，即通过内部边界细分为较小区域的多面体

介绍

在本文中，我们在大规模问题的背景下讨论了点在多面体测试，其中点包含测试必须针对具有大量三角形面的多面体重复调用。此外，我们还将以更一般的方式考虑这个问题，在多材料模型的情况下，内部边界将多面体划分为许多较小的区域，我们需要知道测试点位于哪个区域。

多面体点测试（以及其2D模拟，多边形点测试）是计算机图形学中的基本操作之一，有许多好的作品解决了这个问题[1-16]。最著名的算法是光线交叉法。不幸的是，它遇到了奇点问题[1]。因此，一些研究人员建议随机重复投射光线[2]。为了提高效率，光线交叉通常与层次结构（如八叉树）相结合。

最近，Wang等人[3,4]提出了一种不同的加速策略通过将三角形分类为层来提高光线交叉测试。他们的方法需要找到直线和三角形之间的交点，可能涉及较大舍入误差的操作。另一种众所周知的方法是基于梯形/BSP结构的方法[5,6]。当处理非常复杂的模型时，无论是在时间上还是在内存上，它都有非常高的预处理成本的缺点。Feito等人在一系列论文中开发了另一种算法[5,7,8]。他们的方法是基于累积由多面体的三角形和测试点形成的四面体的符号。它非常简单可靠；然而，由于它需要整个多面体的几何信息，因此对于大规模问题来说效率很低。

有一些方法使用测试点周围的测试多面体的局部信息来执行操作。其中包括Horn等人的以下定理[9]；下面，T表示多面体上所有点集合的补码：
定理：多面体的放大定向定理。

1.如果测试点p∈ T最接近多面体H的面f，如果f的向外法线与源于f点并终止于p的向量之间的内积为负，则p位于H的内部。

2.如果测试点p∈ T最接近多面体H的边e，如果e是凹的，则p位于H的内部。

3.如果测试点p∈ T最接近多面体H的顶点v，如果两条边er和es在v处相交并形成具有空半空间的判别平面L，则p位于H的内部，如果er是凹的。
上述定理是Nordbeck等人工作[10]的扩展，实际上是多面体点算法的理论基础。奇怪的是，显然没有人试图实施它。

闭合平滑曲面的法线长期以来一直用于确定点是在此类曲面内部还是外部。对于表示为三角形网格的多面体，也很容易使用这种方法。Bærentzen等人[11]证明了角度加权伪法线可用于确定点是否位于多面体内部。具体地说，他们表明，角度加权伪法线和矢量之间的内积的符号，从曲面上的最近点到任意点，唯一地决定了该给定点是否在多面体内部。实现Bærentzen等人和Horn等人的想法的关键步骤是有效地找到给定测试点的多面体上的最近点。分层的AABB/OBB树结构用于加快Bærentzen等人的搜索速度。不幸的是，他们代码的测试结果比基于相同结构的光线交叉慢得多。

在本文中，我们遵循同样的想法，尝试使用局部几何结构来定位点。然而，由于引入了一些新的概念和思想，我们的算法通常比基于八叉树或AABB/OBB树的光线交叉性能要好得多。我们还发现，对于多面体测试中的点，八叉树比AABB/OBB更好。

我们的主要思想是在测试点p的附近找到一个可见三角形，然后使用这个可见三角形来找到一个确定三角形，这个确定三角形又被用来确定p的位置。可见三角形和确定三角形的完整定义将在第2节中给出，确定三角形是足以确定测试点位置的三角形。

在描述我们的方法之前，我们给出了一些术语和符号，并定义了两个基本操作。

1.1术语、符号和两个基本操作

多面体中有三个组成元素：顶点、边和三角形（即面）。在开发算法时，我们主要关注三角形。

假设我们有一个多面体Q，它由一组n个顶点V={vi:i=1，2，…，n}和一组m个三角形T={ti:i=1，2，，…，m}组成。假设这些三角形的方向是，任何三角形t的正邻域都在Q内，t的负邻域都位于Q外。要定义三角形的正边，我们遵循右手规则：将手指卷曲，以遵循三角形中三个顶点的顺序，拇指指向的方向是正边。

多边形定向边的2D图解见图1。对于任意测试点p，多面体测试中的点是确定其在Q内还是在Q外。对于多材料多面体，每个三角形都属于两个区域。我们通过在三角形的数据结构中记录两个区域编号，即正区域编号和负区域编号，将三角形的每一侧的区域关联起来。多材料多面体Q的多面体测试点是确定测试点p位于哪个区域。我们的算法需要定义以下两个基本几何运算。

• 定向操作
  这是点p和三角形之间的运算。对于给定的定向三角形，确定p所在三角形的边，即正边或负边。

方向运算通常通过获得4×4矩阵的行列式的符号来完成。这又可以通过精确计算获得，其成本仅略高于近似浮点计算（约1.2倍）；参考Shewchuk提出的方法[17]

• 最近点分类操作
  对于三角形t，此操作将t上最接近p的点分类。有三种可能的分类：在t的一个顶点上（“在顶点上”）；在t的一条边上（“在边上”）；或在t的内部（“面上”）。它们如图2所示。由于p相对于t的关系是所需的，因此不必计算t与p的最近点的位置；相反，在定向操作中，使用行列式的符号来执行操作，这意味着可以避免舍入误差，并且可以调用精确的计算来获得结果。
  

本文的其余部分由五个部分组成。在第2节中，我们给出了确定三角形和可见三角形的定义，然后给出并证明了可见三角形定理，该定理通过演示如何在给定可见三角形的情况下获得确定三角形，为我们的工作奠定了理论基础。我们将给出三种方法来获得点p的可见三角形；第一个在第3节中讨论。我们的算法和实现主要基于八叉树，因此，关于八叉树的一些术语和符号也在第3部分中给出。第二种和第三种获得可见三角形的方法构成第4节；它们分别用最近的三角形和所谓的相对最近的三角形替换可见的三角形；相对最近三角形的定义将在第4节中给出。在第5节中，我们报告了我们的实现代码的一些实验结果。在第6节中，我们讨论了代码中分层树结构的选择，并给出了一些结论。

确定三角形：用于定位p的三角形

Feito等人[4,5]提出了一种使用方向操作定位点的方法。然而，他们的方法必须对每个三角形执行此操作，然后通过累积所有结果来做出决定。只根据一个三角形做出正确的决策应该更有效。

确定三角形（DT）：如果一个三角形足以使用方向操作确定p相对于多面体Q的位置，则称为确定三角形。

可见三角形（VT）：如果我们可以画一条从p到三角形t的线段，而不与Q相交，那么t就是一个可见三角形。可视三角形还可以标记为顶点可见、边可见或面可见，具体取决于线段在t上的结束位置；也就是说，t是顶点可见的，如果线段碰到它的一个顶点，则边可见，如果线段击中它的一条边，则面可见。

如果找到了VT，我们建议对VT和DT之间的关系进行以下观察：如果存在从p到三角形t的线段，而不与Q相交，那么t本身或与t共享边或顶点的三角形将是DT。受Horn等人定理[9]的启发，我们将从VT获得DT的具体规则表示为以下定理：

可见三角定理。假设t是VT，即，假设存在连接p到t的无障碍线段，则可以使用以下三条规则找到DT：

（1） 面部规则：如果t是面部可见的，那么t本身就是DT。

（2） 边规则：如果t的边ab可见，则入射到ab的三角形之一，即与三角形pab成最小角度的三角形，将是DT。

（3） 顶点规则：假设t的一个顶点v是可见的，线段vp和所有入射到v的三角形所跨越的角度都已计算出来，vp和入射到v上的所有边所跨越的角也已计算出来。然后有两个子规则：

a、 如果三角形vxy跨越的角度vp小于所有其他三角形和入射到v的边的角度，则vxy将是DT。

b、 如果与vp形成角度的边va小于所有三角形和其他入射到v的边形成的角度，则va将可见，然后可以使用边缘标尺
证据face规则是自我证明的，并且构成了其他规则的证明基础。对于另一种情况，我们需要找到一个三角形，该三角形具有从p到其内部的无障碍线段。

边缘规则的证明：我们知道存在一个连接p和ab的无障碍线段px，x位于ab的内部。从ab周围的所有三角形中，我们选择一个与三角形pab成最小角度的线段，用t0表示（图3（a）中最右边的三角形），将点x向t0的内部扰动无穷小，表示新位置为x′，如图3

（b）所示。现在我们将讨论px′没有被任何三角形阻挡。

有两类三角形：包围并将ab作为其边之一的三角形和不包围的三角形。考虑第二类中的三角形；如果它阻塞px′，它也会阻塞px，因为扰动是无穷小的，这是不可能的。现在考虑第一类中的一个三角形；如果它阻塞px′，它将跨越pab小于t0的角度。这也是不可能的，因为我们将t0定义为围绕ab的三角形，它跨越pab的最小角度。因此，px′是无障碍的，t0是DT。

顶点规则的证明：我们知道三角形t的顶点v之间存在一个无阻碍的线段pv。构造一个以v为顶点，vp为轴的无限小圆锥体，并逐渐扩展圆锥体（如图3（b）所示），直到它接触到（a）三角形t0的内部或（b）三角形t2的边ab。对于情况（a），与t0相切的圆锥及其交点是t0内部的线段；则px是无障碍的，t0是DT。这是顶点规则第（a）部分中的情况。对于情况（b），t0在边ab处与圆锥体相交。拾取点x i，ab.px的内部无障碍，因此ab是可见边。可以使用边缘规则在ab周围找到DT。这是顶点规则部分（b）中的情况。这证明了可见三角定理。请注意，可见三角形定理对于单材料和多材料多面体情况都是正确的。可见三角形定理指出，即使只有一个顶点可见，VT也总是会产生DT。现在的问题是如何找到VT。我们将给出三种获得VT的方法，都是基于使用八叉树。

在展示我们的算法之前，我们必须首先定义一些关于八叉树的术语和符号。

与八叉树相关的术语和符号，以及找到可见三角形的理论方法

树是有效地在Q的三角形中找到VT的自然方法。在许多层次结构中，我们选择了一个与PM树相似的八叉树[18]。我们将在第6节稍后讨论我们的选择。

八叉树的根对应于长方体区域，即包含整个多面体的边界框。构建树的规则非常简单。顶点一个接一个地插入树中，其中每个叶节点最多包含k个顶点，其中k是预先选择的常数；否则它被分割。构建树后，三角形被插入其中，因此每个叶节点将存储与其重叠的所有顶点和三角形。图4是二维情况的示例；2D情况下的分层结构是四叉树，其中存储有多边形。一个单元实际上是八叉树的一个节点，八叉树也表示一个长方体区域，它是由边界框表示的整个域的一部分。通常一个单元指的是叶节点，但当它在伪代码中作为变量出现时，它有时可能指的是非叶节点。特别感兴趣的单元是包含测试点p的单元；该单元在下文中将被表示为单元。通常cellp指的是叶节点，但在算法2的伪代码中使用时，它也可能指的是非叶节点（可能滥用了约定）。在八叉树的帮助下，多面体点测试的框架如下：

3.1.多面体算法中点的伪码

前两步是预处理。这里构建的树包含三种细胞（叶节点）：黑色、白色和灰色。与任何三角形重叠的单元格标记为灰色；这些重叠的三角形将在单元格中有一个记录。如果一个单元不包含任何三角形的记录，它将被称为空，根据它是在Q的内部还是外部，它将进一步标记为黑色或白色。如果Q是一个多材质多面体，一个黑色单元将与它所属的区域编号相关联。由于八叉树，位于黑色或白色单元中的点的位置将立即可用。如果cellp为灰色，我们需要获取VT，从中生成DT，然后根据Orientation Operation报告p的位置。

现在，我们给出以下简单算法作为获得VT的第一种方法。

3.2.获得VT的理论算法

使用构造的八叉树，可以在本地高效地执行搜索。保证cellp中至少有一个VT，我们寻求获得这样的VT。算法2，在cellp内获得VT：

由算法2找到的三角形t保证是VT。上述算法可以直接从可见三角形定理导出。三角形t是面可见的、边可见的或顶点可见的，这取决于光线击中t的位置，即在其内部、边或顶点。这种找到VT的方法非常有效，因为它只涉及单元格中的三角形。不幸的是，它需要计算光线和三角形之间的交点，这种操作很容易舍入错误。我们称之为“理论”方法，因为我们还没有尝试实现它。

由于p的最近三角形显然是一个VT，因此在实际实现中，我们求助于Horn等人[9]和Bærentzen等人[11]提出的最近三角形方法。假设我们使用最近的三角形t作为VT；t将通过最近点分类操作标记为面可见、边可见或顶点可见，然后我们可以使用可见

三角定理。我们不仅必须识别普通的最近三角形（下文称为全局最近三角形），而且还必须识别另一种称为相对最近三角形（RCT）的最近三角形。下一节将解释RCT在算法中的含义和作用。

4.全局和相对最近的三角形，以及它们作为可见三角形的作用

全局最近三角形（GCT）：Q的所有三角形中最接近p的三角形。

相对最近三角形（RCT）：设x为三角形t上与p最近的点。如果t在与一组相邻单元重叠的所有三角形中与p最接近，由p到x的线段相交，则t是相对最近的三角形。

GCT保证是VT。我们将在本节稍后介绍，也可以将RCT用作VT；RCT具有比GCT更快获得的优点。然而，由于RCT的概念来源于GCT的概念，我们首先给出了获得GCT的算法。

4.1.GCT算法

在实践中，我们通过应用最初来自[18]的想法改进了算法3：当必须搜索多个子树时，表示更靠近p的区域的子树会在距离更远的子树之前搜索。

使用GCT，我们可以使用可见三角形定理从相邻三角形中找到DT。此外，由于GCT是一种特殊的VT，可见三角定理的顶点规则可以简化如下：
GCT的顶点规则：假设GCT t的顶点v可见。如果边va与vp形成的角度小于入射到v的所有其他边的角度，则va将可见，并且在va周围有一个DT，可以根据边规则找到。

算法3和基于它的多面体点测试已经非常有效。然而，我们将证明，获得RCT的效率更高，这个RCT可以承担与GCT相同的角色。

4.2.RCT算法

我们已经给出了上述RCT的定义；现在我们将对其进行说明，并给出找到它的算法。为了便于演示，我们将其显示在2D案例中。GCT和RCT的2D类似物分别是GCE（全局最近边缘）和RCE（相对最近边缘）。例如，考虑确定图5（a）中所示的边ab是否为RCE的问题。首先，画一条从p到x的线段，ab到p的最近点。请注意，px穿过四个单元格，这些单元格中的其他四条边比ab离p更远。因此，ab是RCE。请注意，GCE保证是RCE，但RCE不一定是GCE。在图5（a）中，cd是GCE，也是RCE；ab是RCE但不是GCE。

现在的问题是如何找到RCE。在GCE的算法中，我们从cellp开始。在cellp中，找到最近的边，用e表示。找到e上最接近p的点，用x表示

图5（b）），则e为RCE。如果x不在cellp中，我们必须按照RCE的定义检查它。沿着线段px从p移动到x，检查途中所有单元格的边缘；如果没有发现比e更近的边，如图5（a）所示，那么e是RCE。否则，如图5（c）所示，将有一条较近的边缘作为RCE的新候选。我们用更近的边更新e（参见图5（c）中的边cd），然后再次开始检查过程。该过程将收敛到RCE，因为候选边越来越接近p。

尽管这些想法在2D案例中得到了展示，但在3D案例中却完全相同。如前所述，RCT可以扮演与GCT相同的角色并充当VT，并且找到RCT比GCT更快。

5.实验结果、与基于八叉树/AABB的光线交叉的比较以及成本的理论分析

我们在Windows XP下使用MS Visual Studio Plus编程环境，在C++中实现了该算法。该算法在IBM笔记本电脑T2400、1.83 GHz Intel CPU和512 MB RAM上进行了测试。

为了验证单一材料情况下的代码，我们下载了John Burkardt编码的球形多边形方法[12]。在单材料情况下，我们考虑多材料情况。

对于多材料多面体，我们在每个区域内创建相等数量的测试点，将其存储在一个文件中，然后，由于其位置已知，因此使用该文件来测试我们的代码。已经测试了数十亿个点，到目前为止没有发现错误。

为了比较效率，下面将在本节中讨论八个示例的测试结果。我们首先测试了由AutoCAD创建的三个机械模型（图6至8），并由第一作者开发的网格生成器进行网格划分，以及兔子、马和龙的三个自然主义模型

（图9，来自JA Bærentzen）。为了显示GCT和基于八叉树的射线交叉（RC）的几何依赖性，我们构建了图10所示的盒模型；这个盒子包含五层，每层都是一个长方体，其左边的密度是右边的十倍。图11所示的最后一个例子是一个具有15个区域和分形特征的多材料模型。

该模型是从一家开发硅模拟软件的公司获得的，代表了模拟芯片蚀刻过程的最终结果。

5.1.基于GCT/RCT/八叉树的RC之间的效率比较

为了进行比较，除了基于八叉树的GCT和RCT之外，我们还实现了基于八叉树的RC。对于基于八叉树的RC，当p位于黑色或白色单元格中时，采用相同的想法直接报告结果。我们的实现不处理基于八叉树的RC的奇点，因此它偶尔会给出错误的结果。如果这些奇点得到解决，实现可能会运行得更慢。

在测试模型时，我们随机抽取了分布在模型边界框中的100000个点。表1和表2列出了这些测试的统计数据。从表1所示的实验结果可以得出结论：RCT确实比GCT快；例如，在齿轮箱模型中，RCT比GCT快2.1倍。当几何形状从机械模型、动物模型到芯片模型（具有分形特征）发生变化时，它们之间的差距从两倍增加到几乎五倍。

对于基于八叉树的RC，自然主义动物模型的查询时间最快，因为有巨大的空白空间，因此光线通常只会碰到边界一两次。

例如，Horse模型的查询时间仅为0.25秒。即便如此，基于八叉树的RC仍然比动物模型中的RCT慢两倍。对于其他型号，尤其是机械部件，差距很大，比RCT慢3到16倍。基于RCT/GCT/八叉树的RC的预处理时间和内存需求是相同的，因此仅在表2中列出，我们将RCT与[11]中基于AABB的RC进行了比较。

Bærentzen等人基于AABB和OBB实现了两点多面体算法。我们仅将我们的方法与基于AABB的RC进行比较；对这一限制的解释是合理的。Bærentzen等人报告说，OBB比AABB更高效[11]。在我们的机器上安装他们的软件并仔细测试之后，我们相信Bærentzen等人一定混淆了AABB和OBB的结果。在我们手中，当使用AABB树时，他们的代码运行得更快，占用的空间更少，这已经得到了Bærentzen的确认。它们的符号距离计算是确定点位置的另一种方法；由于它比基于AABB的RC慢得多，因此不考虑进行比较。

如表2所示，在预处理时间、内存需求和测试方面，RCT比基于AABB的RC要好得多

黑白丢弃有益于基于八叉树的方法。如果只测试灰色节点中的点，基于八叉树的方法将变慢，而基于AABB的RC将保持稳定；表3列出了RCT的较慢因素。然而，在这种情况下，RCT仍然比基于AABB的RC快。

我们已经下载并测试了最新版本的CGAL，CGAL 3.6。它在这一特定操作中的效率很低：要执行100000点围堵测试，马模型需要5分钟，钻柄需要20分钟。

最好将我们的方法与其他方法进行比较，例如Wang等人的基于层的方法[8]；不幸的是，这些其他方法的源代码尚未公开。

5.2.时间和空间成本的理论分析

基于八叉树的算法有两种最坏的情况。一种情况是，超过k个顶点非常接近，因此树的深度变得非常高。这可以通过设置深度限制并允许某些特殊节点包含超过k个顶点来避免。因此，节点数可以估计为O（n），其中n是顶点数。另一种最糟糕的情况是，节点被每个三角形切割。虽然不清楚实际中会出现这种情况，但在这种情况下，空间成本将为O（nm），查询时间将为0（n），其中m是三角形的数量。然而，对于现代计算机图形学和网格生成社区中流行的网格几何模型，这些场景是极不可能的。

通常，存储在一个节点中的三角形的数量可以被认为是固定的，空间成本为O（n+m）；因此八叉树将大致平衡，深度为O（logn）。那么cellp（一个包含p的节点）的搜索时间将为O（logn）。随后搜索小区内或小区周围的RCT可以被认为是常数，即O（1）；然而，这个估计很难证明。因此，总的估计渐近查询将是O（logn）+O（1）=O（log n）。在实践中（至少在少于一百万个三角形的情况下），获得cellp的O（logn）时间只占整个查询时间的一小部分，因为与后续搜索RCT不同，它不涉及几何计算。

6.讨论和结论

Bærentzen等人的论文[11]的主要目的是计算有符号距离场。AABB树是计算距离场的一个很好的选择，但对于点定位来说是一个糟糕的选择。使用八叉树，可以立即定位黑白单元格中的点。因此，基于八叉树的RC通常比基于AABB的RC运行得更快。

对于基于最近三角形的方法，两个数据结构之间的性能差距甚至更大。这是因为p离边界越远，找到它最近的三角形就越昂贵。使用八叉树，我们可以避免当p在黑色或白色单元格中时计算成本高昂的情况，并且仅当p在灰色单元格中并且因此接近边界时才搜索最近的三角形。这是我们实施GCT时采用的方法，也是它表现相当好的原因。我们使用RCT的想法进一步改进了对VT的搜索，从而进一步提高了整个操作的效率。然而，我们必须承认，我们的方法仍然具有与普通基于八叉树的射线交叉相同的渐近时间效率。

除了实施问题，我们扩展了Nordbeck等人[10]和Horn等人[9]的工作；我们提出了VT、DT和可见三角定理的概念，由于它们的简单性，可以很容易地扩展到多材料多面体。此外，我们相信这种方法也可以很容易地扩展到更高的层面；事实上，我们在流式细胞术的推动下，将这项工作扩展到了N维。最后，我们的代码可以在http://ptinpoly.pbwiki.com/.