首页 关于
树枝想去撕裂天空 / 却只戳了几个微小的窟窿 / 它透出天外的光亮 / 人们把它叫做月亮和星星
目录

7.5 多面体中的点
(Point in Polyhedron)

判断一个点是否位于多面体内部有许多应用，包括碰撞检测。例如，确定一个移动点(如工具尖端)是否穿透了环境中的物体。与二维情况一样，如果多面体是凸的，这个问题就很容易解决。 事实上，第四章求解凸包的代码在 AddOne(Code4.15)的第一步就解决了这个问题。 对于非凸情况，同样可以采用二维的两种解决方案: 绕数计算的推广和射线交叉计数。

立体角 (Solid Angles)

绕卷数 (winding number)的概念在三维空间和二维空间中都适用。这依赖于有符号立体角(signed solid angle)的概念，它衡量的是以某个点为顶点的圆锥所占据球面面积的比例。 立体角以球面度(steradians)为单位。全球面的立体角为 \(4\pi\)。顶点为 \(q\) 底面为 \(T\) 的四面体的立体角，是指当 \(q\) 位于单位球 \(S\) 的球心， 并将入射于 \(q\) 的面(如有必要)延伸以切割 \(S\) 时，落在四面体内的单位球 \(S\) 的表面积。该角度的符号取决于 \(T\) 的方向。 如果将 \(q\) 与多面体 \(P\) 的每个面形成的立体角求和，当 \(q \in P\) 时，结果为 \(4\pi\)。当 \(q \notin P\)，结果为零。 这引出了一个优雅的算法，判定点在多面体内。但遗憾的是，它与其二维对应算法一样，存在同样的实际缺陷。容易受到数值误差的影响，且速度较慢。 下文将介绍的射线交叉法(ray-crossing)与立体角方法的一种实现(Carvalho & Cavalcanti 1995)之间的计时比较，后者的速度慢二十五倍。然而，他们的代码要短得多。

射线交叉 (Ray Crossing)

三维射线交叉算法背后的逻辑与二维版本完全相同。当且仅当从 \(q\) 射向无穷远的射线穿过 \(P\) 的边界奇数次时，\(q\) 才在 \(P\) 的内部。 射向无穷远的射线可以有效地用一条足够长的线段 \(qr\) 来模拟，只要 \(r\) 足够远，确保它在 \(P\) 的外部即可。 正如我们在 7.3 节中已经解决了线段与三角形的相交问题，统计射线交叉次数似乎很容易。 这种方法的难点在于，如何在各种可能出现的退化情况下，准确的统计交叉次数。\(qr\) 可能落在 \(P\) 的一个面上，可能击中一个顶点，可能与一条边共线重叠，可能与某条边横向相交，等等。 似乎可以做出适当的统计，但我不知道有任何这方面的尝试。我将此作为一个开放性问题(练习 7.5.2[1])。

我们基于两点观察展开讨论。首先，对于 \(P\) 的每个面 \(f\)，很容易统计一条非退化的射线穿过 \(f\) 的次数。射线要么完全错过(miss) \(f\)，要么穿过其内部的一点。 其次，在这个意义上，大多数射线都是非退化的，所以随机射线很可能是非退化的。我们的计划是生成一条随机射线并检查是否有退化情况。 如果没有，则通过射线穿过 \(f\) 的次数即可回答问题。如果发现退化，则丢弃该射线再随机生成一条射线。 退化可以用 7.3 节开发的 SegTriInt(Code7.11)来检测。 由此可得到算法 7.1 中所示的伪代码。

现在我们讨论随机射线的生成。设 \(D\) 为围绕 \(P\) 的最小坐标轴对齐(smallest coordinate-aligned)包围盒 \(B\) 的对角线长度。 \(D\) 可由 \(P\) 顶点的最大和最小坐标轻松计算得出。令 \(R = \lceil D \rceil + 1\)。如果查询点 \(q\) 在 \(B\) 之外，那么它就在 \(P\) 之外。 对于 \(B\) 内的任何查询点，从 \(q\) 出发长度为 \(R\) 的射线必然到达 \(B\) 之外，因为 \(D\) 是 \(B\) 内任意两点之间的最大距离，因此也就到达了 \(P\) 之外。 我们将使用 \(R\) 的值来保证查询射线/线段 \(qr\) 严格到达 \(P\) 之外。

生成长度为 \(R\) 的随机射线可以看作是在半径为 \(R\) 的球面上生成随机点。这是一个已被充分研究的问题，我们在此不作探讨。 本书附带的 sphere.c(用于生成图 4.15 中的 10,000 个点)实现了一种生成此类点的方法。我们将在代码中采用该方法，不再赘述其细节。

在开始编写整个算法的代码之前，我们最后要说明一点。算法 7.1 中的 for 循环对 \(P\) 的每个面都调用 SegTriInt。 射线不仅会错过 \(P\) 的大部分面，而且错过的距离也相当远。这种情况需要一个快速的未命中测试(miss-test)，而无需执行函数 SegTriInt所有的计算。 我们将在实现中包含一个非常简单的包围盒(bounding-box)测试，具体如下: 当读取每个面时(由 ReadFaces 读取，这是一个未显示的简单例程)， 计算一个包围盒并将其存储为两个最小和最大角点 tPointi Box[PMAX][2]。在进行针与面 \(f\) 的完整相交测试之前， 我们首先调用 BoxTest(f, q, r)(Code7.14)来检查查询射线 \(qr\) 是否完全位于包围 \(f\) 的盒子的六个面之一的同侧。 如果因此保证不相交，则返回'0'，否则返回'?'。我的测试表明，这种简单的剔除操作，处理了超过一半的相交检查，对于剩下那一半来说，这点轻微的额外开销是非常值得的。

        char BoxTest ( int n, tPointi a, tPointi b )
        {
            int i;      /* Coordinate index */
            int w;
        
            for ( i=0; i < DIM; i++ ) {
                w = Box[ n ][0][i];     /* min: lower left */
                if ( (a[i] < w) && (b[i] < w) ) return '0';
                w = Box[ n ][1][i];     /* max: upper right */
                if ( (a[i] > w) && (b[i] > w) ) return '0';
            }
            return '?';
        }

现在我们介绍用于测试查询点 \(q\) 是否在多面体内的函数 InPolyhedron。我们将返回值设计成如下形式:

返回值 \(\{V, E, F\}\) 继承自 SegTriInt。因为我们已经确保 \(r\) 严格位于 \(P\) 之外，如果 \(qr\) 有一个端点位于 \(P\) 的顶点、边或面上， 那么它必然是 \(q\) 端点。因此，我们可以立即区分 \(q\) 在边界上的情况。\(\{i, o\}\) 则通过交叉计数器(crossings counter)的奇偶性来区分。

        char InPolyhedron( int F, tPointi q, tPointi bmin, tPointi bmax, int radius )
        {
            tPointi r;      /* Ray endpoint. */
            tPointd p;      /* Intersection point; not used. */
            int f, k = 0, crossings = 0;
            char code = '?';
        
            /* If query point is outside bounding box, finished. */
            if ( !InBox( q, bmin, bmax ) )
                return 'o';
        
        LOOP:
            while( k++ < F ) {
                crossings = 0;
        
                RandomRay( r, radius );
                AddVec( q, r, r );
        
                for ( f = 0; f < F; f++ ) { /* Begin check each face */
                    /* Intersect ray with face f and increment crossings: see BELOW. */
                } /* End check each face */
        
                /* No degeneracies encountered: ray is generic, so finished. */
                break;
            } /* End while loop */
        
            /* q strictly interior to polyhedron iff an odd number of crossings. */
            if( ( crossings % 2 ) == 1 )
                return 'i';
            else return 'o';
        }

函数 InPolyhedron 的整体结构如 Code7.15 所示。 首先，对于 \(P\) 的包围盒之外的查询点，它会在执行 InBox( q, bmin, bmax ) 测试后会立即返回。 简单函数 InBox 未在此处显示。然后在一个近乎无限的循环里，向 \(q\) 添加一条随机射线得到 \(r\)。 纯粹是出于编程实践的考虑，重复次数设置了上限。接下来，使用 for 循环对 \(qr\) 与 \(P\) 的每个面进行测试。for 循环的代码单独显示在 Code7.16 中。 如果 for 循环运行完成，那么我们可以确定该射线是通用的(generic)，穿越次数的奇偶性将决定结果。

        for ( f = 0; f < F; f++ ) {     /* Begin check each face */
            if ( BoxTest( f, q, r ) == '0' )
                code = '0';
            else code = SegTriInt( Faces[f], q, r, p );
        
            /* If ray is degenerate, then goto outer while to generate another. */
            if ( code == 'p' || code == 'v' || code == 'e' ) {
                printf("Degenerate ray\n");
                goto LOOP;
            }
        
            /* If ray hits face at interior point, increment crossings. */
            else if ( code == 'f' ) {
                crossings++;
                printf( "crossings = %d\n", crossings );
            }
        
            /* If query endpoint q sits on a V/E/F, return that code. */
            else if ( code == 'V' || code == 'E' || code == 'F' )
                return( code );
        
            /* If ray misses triangle, do nothing. */
            else if ( code == '0' )
                ;
        
            else
                fprintf( stderr, "Error, exit(EXIT_FAILURE)\n" ),
                exit (1);
        
        /* End check each face */
        }

如 Code7.16 所示，for 循环首先使用函数 BoxTest，希望快速判定射线是否错过了 \(f\)。 否则，将调用函数 SegTriInt 并使用其返回码进行后续决策。 只有当返回 'f' 时，\(qr\) 与面在其相对内部相交，穿越计数器才会递增。 \(\{p, v, e\}\) 都表示退化情况: 射线位于 \(f\) 所在的平面内，或者穿过 \(f\) 的顶点或边。 我们不需要对其进一步区分，因为即使射线完全错过 \(f\)，我们仍然可以安全地将其作为退化情况，并等待一条更好的射线。 对于所有这些退化情况，for 循环都会被放弃，控制权返回到无限 while 循环以获取另一条随机射线。如前所述，\(\{V, E, F\}\) 允许立即退出。

例子：立方体。 图 7.9 展示了一个简单的例子。这里 \(q = (5, 5, 5)\) 位于一个由 12 个三角形面组成的 \(10 \times 10 \times 10\) 立方体的中心。 当 \(D = \sqrt{300}\) 时，射线半径为 \(R = 19\)。在某次特定试验中，调用函数 RandomRay 得到 \(r = (23, 6, 11)\)， 该点远在 \(P\) 之外。对 12 个面分别进行的测试，BoxTest 判定了 8 个面， 函数 SegTriInt 判定了4 个面，其中只有一个返回 '1'。因此恰好有一次射线穿越，判定 \(q\) 在内部。
例子：非凸多面体。 我们在如图 7.10 所示的多面体 \(P\) 进行了一个更严格的测试。该多面体具有 \(V = 400\) 个顶点和 \(F = 796\) 个三角形面。 通过在 \(P\) 的包围盒子内生成随机查询点来测试代码的性能。在生成的 1,000,000 条随机射线中，有 8,121 条(0.8%)是退化的，导致 while 循环重试。 在 110 次尝试(0.01%)下，循环再次失败。而在这一百万次试验中，仅有一次实例是循环在找到一条通用射线之前生成了三条随机射线。 虽然这个多面体很难说是"典型的(typical)"(不管那意味着什么)，但我预计性能不会比这 99% 的"命中率"差多少。

7.5.1. 分析 (Analysis)

算法 7.1 的期望运行时间为 \(O(\rho n)\)，其中 \(\rho\) 是 while 循环找到通用射线(generic ray)所需的期望迭代次数。 虽然我们刚刚看到，对于一个组合稠密的样本多面体中 \(\rho \approx 1.01\)，但如果能证明一个理论界限会更令人放心。 我没有进行精确的分析(练习 7.5.2[2])，但会提供一个论证，来说明对于任意 \(\epsilon > 0\)，都可以满足 \(\rho = 1 + \epsilon\)。

我们从两个简化的观察开始。首先，分析那些端点的整数坐标落在包围立方体 \(C\) 表面上的随机射线，比分析落在包围球面上的射线更容易。 这并不影响一般性，因为我们可以修改实现方式以采用这种不太美观的策略，或者选择一个足够大的球体，以包含指向所有立方体表面点的射线。 设立方体的每条边长为 \(L\)。

其次，我们只需要关注点 \(q\) 和多面体 \(P\) 的一条边 \(e\) 之间的退化情况(degeneracy)。如果 \(q\) 位于某个面的平面内， 那么存在射线 \(qr\) 与该面的边产生 \(q-e\) 退化。如果从 \(q\) 发出的射线穿过一个顶点，它也会穿过每条(闭)关联边。因此，让我们只需关注多面体的一条边 \(e\)。

如果 \(e\) 离 \(q\) 足够近，那么它将"投影(projects)"成一条完全横切包围立方体 \(C\) 某个面的线段，如图 7.11 所示。 在最坏的情况下，the segment renders \(L\) integer points on that fae of \(C\) unusable as ray tips，因为它们在某种意义上会导致射线退化。 如果 \(P\) 有 \(E\) 条边，那么 \(C\) 的一个面上最多有 \(EL\) 个点会变得不可用。实际上，\(P\) 的边在立方体面上产生了一个线排列(arrangement of lines)。 只有错过该排列的射线才是通用的。但这样的射线有很多，因为立方体的这个面包含 \(L^2\) 个整数点。因此，命中退化情况的概率至多为 \(EL/L^2 = E/L\)。 因为 \(E\) 是一个常数(图 7.10 中为 1,194)，选择足够大的 \(L\) 就能保证达到任意期望的 \(\epsilon = E/L\)。

然而，实际应用时就又是另外一回事。\(L\) 越大，就约容易出现溢出的问题。在实践中，我发现最好选择尽可能小的 \(R\)，在此对应于 \(L/2\)，仅比盒子对角线 \(D\) 大 1。

7.5.2. 练习 (Exercises)




Copyright @ 高乙超. All Rights Reserved. 京ICP备16033081号-1