4.3 增量算法的实现
(Implementation of Incremental Algorithm)
尽管增量算法在概念上清晰明了，但实现并不简单。在本节中，我们将深入探讨一个完整的实现方案，这也是本书中最复杂的实现。 上一节中省略的细节，本节将结合代码一一补充。对代码不感兴趣的读者， 可直接跳至 4.3.5 节关于"体积溢出 volume overflow"的讨论。
4.3.1 数据结构
找到一种表示多面体表面的最佳方式并不容易，一些文献中提到了多种复杂的方案。我们将在 4.4 节中探讨其中几种思路。 这里，我们将采用一种非常简单的数据结构，不过其适用范围有限。我们假设所处理的多胞体(polytope)表面是三角化的，即每个面都是三角形。这将简化我们的数据结构， 但对非三角化多胞体，例如立方体，这种表示方式多少有些不方便。此外，我们的数据结构不具备某些其他结构所拥有的对称性，也会使部分操作略显繁琐。 尽管存在这些缺点，我仍认为它是最易于理解的。
结构定义(Structure Definitions)。定义三种基本数据类型，顶点(vertices)，边(edges)，面(faces)。 所有顶点保存在一个双向循环链表中，所有边和面也各自保存在双向循环链表中。这些列表的结构与第 1 章中使用的多边形顶点列表相同， 见代码 1.2。列表中元素的顺序并无实际意义，我们应将它们视为集合而非有序列表。 链表中的每个元素是一个固定大小的结构体，其中包括指向其他列表的链接。顶点结构体包含该顶点的整数坐标，但不包含指向与其关联的边或面的指针。 我们并不能直接加入此类指针，因为一个顶点可能与任意数量的边和面相关联。 边结构体包含指向该边两个顶点的指针，以及指向其相邻两个面的指针。这两组指针的顺序是任意的，一些更复杂的数据结构会强制规定顺序。 面结构体包含指向构成该三角形面的顶点的指针，以及指向构成该面的三条边的指针。 请注意，此处我们利用了"所有面均为三角形"这一假设。这三种结构的基本字段如Code 4.1 所示。这些结构还需包含其他杂项字段，稍后将予以讨论。
struct tVertexStructure {
int v[3];
int vnm;
tVertex next, prev;
};
struct tEdgeStructure {
tFace adjface[2];
tVertex endpts[2];
tEdge next, prev;
};
struct tFaceStructure {
tEdge edge[3];
tVertex vertex[3];
tFace next, prev;
};
三种主要结构的类型名称，均以字母 t 开头。顶点的类型名为 tVertexStructure，该名称仅用于声明中。
struct tVertexStructure 类型被赋予别名 tsVertex。
此名称仅在分配内存时作为 sizeof 的参数使用。在整个代码中实际使用的类型是 tVertex，它是一个指向顶点列表中某个元素的指针。
边结构和面结构也有类似的命名方式。这些名称通过位于结构体声明之前的 typedef 语句定义，详见Code 4.2。
typedef struct tVertexStructure tsVertex;
typedef struct tvertexStructure tsVertex;
typedef struct tEdgeStructure tsEdge;
typedef tsEdge *tEdge;
typedef struct tFaceStructure tsFace;
typedef tsFace *tFace;
数据结构实例(Example of Data Structures)。我们将通过一个实例来演示凸包(convex hull)代码。该实例将构造由立方体八个角点组成的点集的凸包。 在最终形成立方体的过程中，会先生成一个包含五个顶点的多胞体，我们用它来说明数据结构的设计。称这个多胞体为 \(P_5\)。
顶点列表包含所有输入点，但并非所有点都会被边和面引用。该立方体边长为 10，位于坐标系的正象限中。所谓的正象限，应该是指 x,y,z 三个坐标都为正数的那个象限。 此处分配给顶点(以及边和面)的索引在代码中不起任何作用，因为所有引用均通过指针进行。
多胞体 \(P_5\) 由 9 条边和 6 个面组成。表 4.5–4.7 中的三个列表完全展现了顶点、边、面的数据结构。标签 \(v, e, f\) 上的下标表示对应记录被创建的先后顺序。 请注意，面列表中不包含 \(f_1\) 或 \(f_4\)。这两个面在图示的数据结构快照之前已被创建并删除。
图 4.12 给出了多胞体 \(P_5\) 的一个视图。其中，面 \(f_2, f_5, f_6\) 是可见的。\(f_0\) 位于 \(xy\)-平面上。\(f_3\) 和 \(f_7\) 在背面并且共面， 它俩共同构成立方体的一个正方形面 \((v_0, v_1, v_5, v_4)\)。
面数据结构需要始终维持的一个重要性质，保证其 vertex 字段中的顶点按逆时针顺序排列。这样根据右手定则，可得到一个指向多胞体外部的法向量。
因此，\(f_6\) 的顶点顺序为 \((v_4, v_2, v_5)\)。字段 edge 也需要保持逆时针顺序。因此，\(f_6\) 的边顺序为 \((e_4, e_6, e_8)\)。
代码经常利用顶点的逆时针顺序，但几乎从未需要使用边的逆时针顺序。尽管如此，为了美观以及未来可能超出本文范围的用途，我们仍通过谨慎的交换操作来维护边的顺序。
头指针(Head Pointers)。我们始终为三个列表(顶点、边、面)各维护一个全局的头指针，就像第 1 章(Code 1.2)那样，
将其初始值设为 NULL。参见Code 4.3。
tVertex vertices = NULL;
tEdge edges = NULL;
tFace faces = NULL;
此前已在Code 1.3中展示过，遍历所有顶点、边或面的循环结构的基本形式。 该循环结构假设列表非空，而在初始多胞体构建完成后，这一条件确实成立。
列表的基本操作(Basic List Processing)。三种数据结构都需要四个列表的基本操作。① 分配新元素内存NEW
② 释放元素所占内存FREE ③ 将新元素添加到列表中ADD ④ 从列表中删除旧元素DELETE。
前三个操作已在第 1–3 章中使用过，参见Code 1.4。
请注意，宏定义NEW适用于所有三种结构，需要将类型 tsVertex, tsEdge, tsFace作为参数传入。
第四个宏定义 DELETE 如Code 4.4 所示，它需要妥善管理指针 head，当它指向的单元格 p 正是被删除的那个时，
需将其更新为 head->next。
#define DELETE( head, p ) if ( head ) { \
if ( head == head->next ) \
head = NULL; \
else if ( p == head ) \
head = head->next; \
p->next->prev = p->prev; \
p->prev->next = p->next; \
FREE( p ); \
}
结构体：完整细节。我们还对基本数据结构扩充了一些标志位(flags)和辅助指针，完整的结构定义见Code 4.5 和 Code 4.6。所有新增字段均已添加注释，将在首次使用时进一步解释。
每种数据结构都对应一个 MakeNull 的程序，用于创建新单元、初始化，并将其添加到相应的列表中。参见代码 4.7。
/* 定义顶点索引 */
#define X 0
#define Y 1
#define Z 2
/* 定义布尔类型 */
typedef enum { FALSE, TRUE } bool;
/* 定义标志位 */
#define ONHULL TRUE /* 点在凸包上 */
#define REMOVED TRUE /* 已被移除 */
#define VISIBLE TRUE /* 可见 */
#define PROCESSED TRUE /* 已处理 */
struct tVertexStructure {
int v[3]; // 三维坐标 (x, y, z)
int vnum; // 顶点编号
tEdge duplicate; // 指向关联边的指针(或 NULL)
bool onhull; // True 当且仅当该点在凸包上
bool mark; // True 当且仅当该点已被处理过
tVertex next, prev; // 双向链表指针
};
struct tEdgeStructure {
tFace adjface[2]; // 相邻的两个面
tVertex endpts[2]; // 两个端点
tFace newface; // 指向关联锥面的指针
bool delete; // True 当且仅当该边应被删除
tEdge next, prev; // 双向链表指针
};
struct tFaceStructure {
tEdge edge[3]; // 构成该面的三条边
tVertex vertex[3]; // 构成该面的三个顶点
bool visible; // True 当且仅当该面对新点可见
tFace next, prev; // 双向链表指针
};
/* 顶点创建与初始化函数 */
tVertex MakeNullVertex(void)
{
tVertex v;
NEW(v, tsVertex);
v->duplicate = NULL;
v->onhull = !ONHULL;
v->mark = !PROCESSED;
ADD(vertices, v);
return v;
}
/* 边创建与初始化函数 */
tEdge MakeNullEdge(void)
{
tEdge e;
NEW(e, tsEdge);
e->adjface[0] = e->adjface[1] = e->newface = NULL;
e->endpts[0] = e->endpts[1] = NULL;
e->delete = !REMOVED;
ADD(edges, e);
return e;
}
/* 面创建与初始化函数 */
tFace MakeNullFace(void)
{
tFace f;
int i;
NEW(f, tsFace);
for (i = 0; i < 3; ++i) {
f->edge[i] = NULL;
f->vertex[i] = NULL;
}
f->visible = !VISIBLE;
ADD(faces, f);
return f;
}
4.3.2 例子: 立方体(Example: Cube)
在本节中，我们将通过一个示例来说明程序的运行过程，其输入为一个立方体的八个顶点。我们将逐段讨论相关代码。
主函数(Main): 整个过程可以分为四个部分(见Code 4.8)：读取数据、创建初始多面体、构建凸包、输出结果。 我们会尽可能按线性顺序讨论它们。代码 4.9 展示了各函数之间的调用关系，并用注释编号标明我们讨论它们的先后顺序。
main(int argc, char *argv[]) {
/* （标志位等此处未显示）*/
ReadVertices();
DoubleTriangle();
ConstructHull();
Print();
}
/* 1 */ ReadVertices()
MakeNullVertex()
/* 2 */ DoubleTriangle()
/* 3 */ Collinear()
/* 4 */ MakeFace()
MakeNullEdge()
MakeNullFace()
VolumeSign()
/* 5 */ ConstructHull()
/* 6 */ AddOne()
/* 7 */ VolumeSign()
/* 8 */ MakeConeFace()
MakeNullEdge()
MakeNullFace()
MakeCcw()
/* 9 */ Cleanup()
/* 10 */ CleanEdges()
/* 11 */ CleanFaces()
/* 12 */ CleanVertices()
/* 13 */ Print()
读取顶点(ReadVertices): 下面是一个立方体的示例输入文件:
0 0 0
0 10 0
10 10 0
10 0 0
0 0 10
0 10 10
10 10 10
10 0 10
上述顶点被标记为 \(v_0, \dots, v_7\)，如之前表 4.5 所示。我们通过函数 ReadVertices (Code4.10) 直接从文件中读取这些顶点坐标， 并通过 MakeNullVertex (Code4.7) 构建成一个列表。稍后，我们将详细解释每个顶点的各个字段的含义。
void ReadVertices( void ) {
tVertex v;
int x, y, z;
int vnum = 0;
while ( scanf ("%d %d %d", &x, &y, &z ) != EOF ) {
v = MakeNullVertex();
v->v[X] = x;
v->v[Y] = y;
v->v[Z] = z;
v->vnum = vnum++;
}
}
双三角形(DoubleTriangle): 下一步，也是第一个实质性的步骤，是创建初始多胞体(polytope)。 按照算法 4.1 的思路，应该从四面体开始。 但我发现从双重覆盖三角形(doubly covered triangle)开始会稍微容易一些，以下简称 d-triangle。这是一种拥有三个顶点和两个面的多面体，这两个面除了顶点的排列顺序相反外，其他完全相同。 尽管根据 4.1 节的定义，这算不上是一个严格意义上的多面体，但它拥有与多面体相同的局部关联结构， 这足以满足代码当前的需求了。
人们可能会认为构建一个 d-triangle 并不困难。但实际上，实现起来既复杂又繁琐，原因有几点。 首先，简单地取顶点列表中的前三个点是不够的，因为这些点可能是共线的。虽然双重覆盖可以容忍这种退化情况，但如果一个面的面积为零，它与后续的点构成的四面体体积也将为零，这是我们无法接受的。 因此，我们必须首先找到三个不共线的点。当然，如果我们假设点处于一般位置，就可以避开这种麻烦，但即便是立方体的顶点也不总是处于一般位置。 其次，需要构建特别的数据结构，来描述一些属性。特别是，必须确保每个面中记录的顶点是逆时针排序的。这似乎也是不可避免的。 第三，这些数据结构有些笨重。我并不认为使用更复杂的数据结构，可以避免这些问题。
d-triangle 的构建分为三个阶段:
- 找到三个不共线的点(\(v_0, v_1, v_2\))。
- 创建并链接两个三角形面 \(f_0\) 和 \(f_1\)。
- 找到第四个点 \(v_3\), 不与 (\(v_0, v_1, v_2\)) 共面。
现在，我们来详细讨论一下 DoubleTriangle 的实现，参见 Code4.11。
void DoubleTriangle( void ) {
tVertex v0, v1, v2, v3, t;
tFace f0, f1 = NULL;
tEdge e0, e1, e2, s;
int vol;
/* 寻找 3 个不共线的点。 */
v0 = vertices;
while (Collinear(v0, v0->next, v0->next->next)) {
if ((v0 = v0->next) == vertices )
printf("DoubleTriangle: All points are collinear!\n"), exit(0);
}
v1 = v0->next; v2 = v1->next;
/* 标记顶点为已处理。 */
v0->mark=PROCESSED; v1->mark=PROCESSED; v2->mark=PROCESSED;
/* 创建两个 twin 面。 */
f0 = MakeFace( v0, v1, v2, f1 );
f1 = MakeFace( v2, v1, v0, f0 );
/* 链接相邻的面字段。 */
f0->edge[0]->adjface[1] = f1;
f0->edge[1]->adjface[1] = f1;
f0->edge[2]->adjface[1] = f1;
f1->edge[0]->adjface[1] = f0;
f1->edge[1]->adjface[1] = f0;
f1->edge[2]->adjface[1] = f0;
/* 寻找第四个不共面的点以形成四面体。 */
v3 = v2->next;
vol = VolumeSign( f0, v3 );
while (!vol) {
if ((v3 = v3->next) == v0)
printf("DoubleTriangle: All points are coplanar!"), exit(0);
vol = VolumeSign( f0, v3 );
}
/* 确保 v3 将是第一个被添加的。 */
vertices = v3;
}
- 三个非共线点(Three noncollinear points): 我们只需要检查顶点列表中所有连续三个点的组合就足够了。因为只要不是所有点都共线，那么在这些组合中至少会有一个是不共线的。
共线的检查方法与第1章中使用的方法相同，但现在因为点是在三维空间中，我们不能仅仅依赖叉积的 \(z\) 坐标。
当且仅当 公式 (1.1) 中叉积的每一个分量都为零时，由这三个点确定的三角形面积为零。
这一逻辑在代码 4.12 中体现。
bool Collinear(tVertex a, tVertex b, tVertex c) { return (c->v[Z] - a->v[Z]) * (b->v[Y] - a->v[Y]) - (b->v[Z] - a->v[Z]) * (c->v[Y] - a->v[Y]) == 0 && (b->v[Z] - a->v[Z]) * (c->v[X] - a->v[X]) - (b->v[X] - a->v[X]) * (c->v[Z] - a->v[Z]) == 0 && (b->v[X] - a->v[X]) * (c->v[Y] - a->v[Y]) - (b->v[Y] - a->v[Y]) * (c->v[X] - a->v[X]) == 0; } - 构建面(Face construction): 每个面都是由一个专门的函数
MakeFace构建的，它以三个顶点的指针和一个面指针fold作为输入，参见Code 4.13，构建一个指向这三个顶点的面。如果面指针fold不为空，我们就会利用它来访问边指针。 这个技巧并不深奥。其目标是按照传入的顺序，用三个顶点指针和三个边指针填充面记录。这些边指针要么是全新构建的(针对第一个三角形)， 要么是从fold复制过来的(针对第二个三角形)，最后还要链接每条边的相邻面(adjface)字段。 注意，为每个面实现初始的正确方向很容易，一个面使用 \((v_0, v_1, v_2)\)，另一个则使用 \((v_2, v_1, v_0)\)。tFace MakeFace(tVertex v0, tVertex v1, tVertex v2, tFace fold) { tFace f; tEdge e0, e1, e2; /* 构建初始三角形的边 */ if( !fold ) { e0 = MakeNullEdge(); e1 = MakeNullEdge(); e2 = MakeNullEdge(); } else { /* 以相反的顺序从 fold 拷贝 */ e0 = fold->edge[2]; e1 = fold->edge[1]; e2 = fold->edge[0]; } e0->endpts[0] = v0; e0->endpts[1] = v1; e1->endpts[0] = v1; e1->endpts[1] = v2; e2->endpts[0] = v2; e2->endpts[1] = v0; /* 构造三角形的面. */ f = MakeNullFace(); f->edge[0] = e0; f->edge[1] = e1; f->edge[2] = e2; f->vertex[0] = v0; f->vertex[1] = v1; f->vertex[2] = v2; /* 建立边与面的邻接关系. */ e0->adjface[0] = e1->adjface[0] = e2->adjface[0] = f; return f; } - 第四个非共面点(Fourth noncoplanar point): 我们将搜索一个点 \(v_3\)，使得四面体 \((v_0, v_1, v_2, v_3)\) 的体积不为零，从而找到一个非共面点。 一旦找到这个点，头指针就会被重新定位到 \(v_3\)，这样它就会成为第一个被添加的点。采用这种策略是为了确保我们在下一步能构建出一个合法的、具有非零体积的多面体。 如果允许它在同一个平面内增长，将会使方向计算变得非常困难。
我们在这个立方体示例上运行 DoubleTriangle 函数时，前三个顶点\(v_0, v_1, v_2\)就是非共线的。实际上，输入数据中没有任何三个点是共线的。
随后构建了面 \(f_0\) 和 \(f_1\)。\(f_1\) 将在后续处理中被删除。接着尝试的第一个 v3 候选点是 \(v_3\)，如表 4.5所示，
它实际上与 \((v_0, v_1, v_2)\) 共面。我们稍后会讨论 VolumeSign。所以将头指针 vertices 设置为 \(v_4\)，
因为它不共面，这就为增量算法添加第一个点做好了准备。
构造凸包(ConstructHull): 现在我们来到了该算法的核心部分。值得注意的是，为了达到这一步，我们需要编写多少"外围peripheral"代码。
函数 ConstructHull，如Code4.14所示，是在初始多胞体构建完成后由主程序调用的，它只是简单地调用函数 AddOne 不断地添加一个点。
有一个小细节值得注意，整个顶点列表的处理过程中，通过 v->mark 字段来避免重复处理点。
我们不能接着函数 DoubleTriangle 结束的地方继续处理顶点列表，因为构成那个 d-triangle的顶点可能分散在列表的各个位置。
void ConstructHull(void) {
tVertex v, vnext;
int vol;
v = vertices;
do {
vnext = v->next;
if (!v->mark) {
v->mark = PROCESSED;
AddOne(v);
CleanUp();
}
v = vnext;
} while (v != vertices);
}
将每个点添加到前一个凸包之后，会调用一个重要的函数 CleanUp。它会删除数据结构中多余的部分，并为下一次迭代做准备。我们将在下文详细讨论这一点。
添加一个点(AddOne): 算法的主要工作就是这个 AddOne 函数，如 Code4.15 所示。该过程将一个点 \(p\) 添加到凸包中，
如果 \(p\) 位于外部，则构建新的面锥(cone of faces)。该过程包含两个步骤:
bool AddOne(tVertex p)
{
tFace f;
tEdge e, temp;
bool vis = FALSE;
/* 标记从 p 点可见的面 */
f = faces;
do {
if (VolumeSign(f, p) < 0) {
f->visible = VISIBLE;
vis = TRUE;
}
f = f->next;
} while (f != faces);
/* 如果 p 点看不到任何面, 那么 p 就在凸包的内部 */
if (!vis) {
p->onhull = !ONHULL;
return FALSE;
}
/* 标记可是区域内的边，方便后续删除。并在每个边上构建一个新的面 */
e = edges;
do {
temp = e->next;
if (e->adjface[0]->visible && e->adjface[1]->visible)
/* e interior: mark for deletion. */
e->delete = REMOVED;
else if (e->adjface[0]->visible || e->adjface[1]->visible)
/* e border: make a new face. */
e->newface = MakeConeFace(e, p);
e = temp;
} while (e != edges);
return TRUE;
}
- 确定之前构建的凸包中哪些面对于 \(p\) 是可见的。当且仅当 \(p\) 严格位于，面 \(f\) 所确定的正半空间内时，面 \(f\) 对 \(p\) 可见。
按照惯例，正侧是由 \(f\) 的逆时针方向确定的。这里的"严格"条件是一个至关重要的细节，如果 \(p\) 只是从侧面照向该面，即点落在面所在的平面上，而非其正半空间内，我们不认为该面是可见的。
可见性条件是通过体积来确定的(我们将在下文中讨论它)，当且仅当由 \(f\) 和 \(p\) 确定的四面体体积为负时，\(f\) 对 \(p\) 可见。
如果点 \(p\) 不能看到任何面，那么 \(p\) 一定位于凸包内部，我们将标记它，以便随后删除。 - 给 \(p\) 添加一个面锥。\(p\) 点可见的多胞体部分，其表面形成了一个连通区域。该区域的内部的点需要删除，而锥体则连接到连通域的边界上。
沿着凸包的每条边依次检查。那些两个相邻面都标记为可见的边，一定位于可见区域的内部。将它们标记下来，以便随后删除，但尚未删除。
那些只有一个相邻面可见的边，则位于可见区域的边界上。这些边正是构成以 \(p\) 为顶点的新三角形面的底边。
构建这个新面的大部分工作由函数
MakeConeFace完成的。
这段代码的棘手之处在于，循环遍历所有边的同时，MakeConeFace还要向列表中添加新边。 所有边都是立即插入到edges列表头之前的。因此，新创建的边会被重新处理。但是两个 if 语句的判定对这些边都会失效， 因为它们的相邻面在创建时，其可见标志(visible flag)被设置为FALSE。
函数 AddOne 被编写为根据凸包是否被修改来返回 TRUE 或 FALSE，但此处显示的代码版本并没有使用这个布尔值。
体积符号(VolumeSign): 在 1.3.8 节提到，顶点为 \((a, b, c, d)\) 的四面体的体积是下列行列式的 \(1/6\)
$$ \begin{equation}\label{f4.6} \begin{vmatrix} a_x & a_y & a_z & 1 \\ b_x & b_y & b_z & 1 \\ c_x & c_y & c_z & 1 \\ d_x & d_y & d_z & 1 \end{vmatrix} \end{equation} $$我们可以直接将该行列式展开得到该体积的代数表达式。函数 VolumeSign
以一种不同于方程(1.12) 的方式来表达其计算过程，如 Code4.16 所示，它们在代数上是等价的，但使用的乘法运算更少。
它的核心思想是，平移四面体使得 \(p\)-corner 位于原点。各个坐标被繁琐地赋值给许多不同的变量，以便更容易且无误地描述体积方程。
int VolumeSign(tFace f, tVertex p) {
double vol;
double ax, ay, az, bx, by, bz, cx, cy, cz;
ax = f->vertex[0]->v[X] - p->v[X];
ay = f->vertex[0]->v[Y] - p->v[Y];
az = f->vertex[0]->v[Z] - p->v[Z];
bx = f->vertex[1]->v[X] - p->v[X];
by = f->vertex[1]->v[Y] - p->v[Y];
bz = f->vertex[1]->v[Z] - p->v[Z];
cx = f->vertex[2]->v[X] - p->v[X];
cy = f->vertex[2]->v[Y] - p->v[Y];
cz = f->vertex[2]->v[Z] - p->v[Z];
vol = ax * (by*cz - bz*cy)
+ ay * (bz*cx - bx*cz)
+ az * (bx*cy - by*cx);
/* The volume should be an integer. */
if (vol > 0.5) return 1;
else if (vol < -0.5) return -1;
else return 0;
}
读者可能注意到代码做了一件奇怪的事。它接受整数坐标作为输入，转换为浮点数进行计算，最后返回集合 \(\{-1, 0, +1\}\) 中的一个整数。我们将在 4.3.5 节中讨论这个技巧。
按照右手定则，当 \(p\) 位于 \(f\) 的负侧时，体积为正。现在将点 \(v_6 = (10, 10, 10)\) 添加到图 4.12 中的多面体 \(P_5\)上。它可以看到面 \(f_6\)， 从"外部"看其顶点的逆时针顺序为 \((v_4, v_2, v_5)\)。\(f_6\) 和 \(v_6\) 的行列式为:
$$ \begin{vmatrix} 0 & 0 & 10 & 1 \\ 10 & 10 & 0 & 1 \\ 0 & 10 & 10 & 1 \\ 10 & 10 & 10 & 1 \end{vmatrix} = -1,000 < 0 $$在函数 AddOne 中，这个负体积被解释为 \(v_6\) 可以看到 \(f_6\)。
AddOne(cube example): 在讨论 AddOne 之前，我们先用立方体示例说明其功能。
顶点列表中的前三个顶点 \(v_0, v_1, v_2\) 被标记为 DoubleTriangle 。如前所述，因为 \(v_3\) 与前三个顶点共面，目前头指针被移动到 \(v_4\)。
然后按 \(v_4, v_5, v_6, v_7, v_3\) 的顺序添加顶点。设 \(P_i\) 为添加顶点 \(v_i\) 后的多面体。随后按 \(P_2, P_4, P_5, P_6, P_7\) 和 \(P_3\) 的顺序生成多面体。它们如图 4.13(a)-(f) 所示。
我们来看看添加 \(v_6\) 时，如何从 \(P_5\) 转换到 \(P_6\) 的。从图 4.13(c) 和图 4.12可以明显看出，\(v_6\) 只能看到面 \(f_6 = (v_4, v_2, v_5)\)。 可见性通过 \(v_6\) 与 \(P_6\) 的所有面组成的四面体的体积来判定，对 \(f_6\) 返回 -1。对面 \(f_0, f_3 f_7\) 返回 +1，对 \(f_2\) 和 \(f_3\) 返回 0。 注意，根据我们的可见性定义，代码不会将共面的 \(f_2\) 和 \(f_5\) 标记为可见。
AddOne 的第二部分在可见区域的内部找不到任何边，因为它仅由 \(f_6\) 组成。并且 \(f_6\) 的每条边，\((e_4, e_6, e_8)\)，都是边界边，
因此以这些边为底构建三个新面 \(f_8, f_9 f_{10}\)。最初，这些面被链接到字段 e->newface 中，从而在添加锥体时可以保持旧凸包数据结构的完整性。
这使得在构建新结构的同时，可以通过 MakeConeFace 对旧结构进行查询处理。只有在整个锥体附加完成后，
才会通过 CleanUp 对数据结构进行清理。
再探共面性Coplanarity Revisited: 现在我们回到共面性的问题上来。如果我们认为从 \(v_6\) 可以看见 \(f_2\)，那么 \(f_2\) 的两条边 \(e_0, e_3\) 就是边界， 而 \(e_4\) 将位于可见区域内部。此时的锥体将基于四条边而不是三条边。因此，我们将共面的面视为不可见的，这样可见区域，以及由此产生的新锥体，将尽可能小。
将 \(vol == 0\) 的面判定为不可见的有两个原因:
- 数据结构的变化被最小化，因为正如刚才所解释的，可见区域被最小化了。
- 落在旧凸包面上的点会被丢弃。
如果我们把零体积的面视为可见的，那么位于面内部的点将"看到"该面，从而导致不必要地将其分割成新的面。
虽然这种可见性处理方法避免了在旧面的内部插入新点，但它并不能排除所有不必要的共面点。如果在构造过程中首先遇到内部点，它将永远不会被删除。 它与我们在3.5 节中介绍的 Graham 二维凸包算法的代码不同，三维凸包代码对于相同输入点的不同排列可能会产生不同的输出。 我们可以通过后处理，删除不必要的共面点，来实现关于排列的不变性(练习 4.3.6[11])。
代码中还有两个主要部分需要解释，它们都用于管理数据结构: MakeConeFace 和 CleanUp。
MakeConeFace: 如 Code4.17 所示，函数 MakeConeFace以一条边 \(e\) 和一个点 \(p\) 作为输入，创建一个由 \(e\) 和 \(p\) 张成的新面，
以及连接 \(p\) 和 \(e\) 的端点之间的两条新边。返回指向该面的指针，并且创建的结构体被正确地链接在一起。
tFace MakeConeFace( tEdge e, tVertex p ) {
tEdge new_edge[2];
tFace new_face;
int i, j;
/* 构建两条新边 */
for ( i=0; i < 2; ++i )
/* 如果边已经存在, 就拷贝到 new_edge. */
if ( !( new_edge[i] = e->endpts[i]->duplicate) ) {
/* 否则, duplicate is NULL, 调用 MakeNullEdge. */
new_edge[i] = MakeNullEdge();
new_edge[i]->endpts[0] = e->endpts[i];
new_edge[i]->endpts[1] = p;
e->endpts[i]->duplicate = new_edge[i];
}
/* 构建新面. */
new_face = MakeNullFace();
new_face->edge[0] = e;
new_face->edge[1] = new_edge[0];
new_face->edge[2] = new_edge[1];
MakeCcw( new_face, e, p );
/* 记录面的邻接顶点 */
for ( i=0; i < 2; ++i )
for ( j=0; j < 2; ++j )
/* Only one NULL link should be set to new_face. */
if ( !new_edge[i]->adjface[j] ) {
new_edge[i]->adjface[j] = new_face;
break;
}
return new_face;
}
这种方式最直接，但有两个细节需要注意。首先，必须避免创建重复的边。因为我们并不构建可见区域的边界结构，所以锥体的面是以任意顺序构建的。 一旦锥体的一个面及其边被创建，随后的面可能与先前创建的面共享两条、一条或没有边。
我们可以通过如下的机制检查它。每创建一条边 \(e_i\)，其一端在 \(p\)，另一端在旧凸包上的顶点 \(v\) 时，就将 \(v\) 的字段 v->duplicate 就指向 \(e_i\)。
对于任何不与已构建的锥体边相关联的顶点，字段 duplicate 为 NULL。注意，每个顶点最多与一条锥体边相关联。
对于可见区域边界上的每条边 \(e\)，总是创建一个新面 \(f\)。但是，只有当 \(e\) 的端点 \(v\) 的 duplicate 字段为 NULL 时，才会为 \(f\) 创建一条新边 \(e\)。 如果该字段不为 NULL，则使用该字段指向的已创建的锥体边来填充 \(f\) 的相应边字段。
MakeConeFace 中的第二个复杂之处是需要将 \(f\) 的 vertex 字段中的数组元素按逆时针顺序排列。
这是由函数 MakeCcw 处理的。基本思想很简单，假设旧凸包的面已经有了正确的朝向，只需使新面与旧方向保持一致即可。
特别是，锥体面 \(f\) 可以继承与形成其底边的旧凸包的边 \(e\) 相邻的可见面相同的方向。这是因为新面隐藏了旧面，并且在某种意义上是它的替代。因此它自然地采用了相同的方向。
正是在这里，我们选择的数据结构的尴尬之处就显现出来了。因为 \(e\) 的方向是任意的，我们必须弄清楚 \(e\) 相对于可见面的方向是什么。
也就是说，可见面的哪个顶点指针 \(i\) 指向 \(e\) 的 base[0]端。然后我们可以从这个索引 \(i\) 开始锚定决策。
虽然在显示的代码中不需要，我们也交换新面 \(f\) 的边以遵循相同的方向。因为 \(e\) 在 MakeConeFace 中被设置为 edge[0]，
当它们与可见面的方向相反时，我们将 edge[1] 与 edge[2] 交换。参见Code4.18。
void MakeCcw( tFace f, tEdge e, tVertex p ) {
tFace fv; /* 与 e 相邻的可见面 */
int i; /* e->endpoint[0] 在 fv 中的索引 */
tEdge s; /* 临时变量，用于交换 */
if ( e->adjface[0]->visible )
fv = e->adjface[0];
else fv = e->adjface[1];
/* 设置 f 的 vertex[0] 和 [1]，使其方向与 fv 的对应顶点一致。 */
for ( i=0; fv->vertex[i] != e->endpts[0]; ++i )
;
/* 使 f 的方向与 fv 相同。 */
if ( fv->vertex[ (i+1) % 3 ] != e->endpts[1] ) {
f->vertex[0] = e->endpts[1];
f->vertex[1] = e->endpts[0];
} else {
f->vertex[0] = e->endpts[0];
f->vertex[1] = e->endpts[1];
SWAP( s, f->edge[1], f->edge[2] );
}
f->vertex[2] = p;
}
#define SWAP(t,x,y) { t = x; x = y; y = t; }
CleanUp: 在 AddOne 之后调用 CleanUp 之前，新的凸包已经构建完成了，所有的面、边和一个新顶点都已经相互链接，
并且正确地链接到了旧的结构上。但是，这个锥体是通过可见区域边界上 newface 字段"粘贴"在旧结构上的。
而且，锥体内部的旧结构需要被删除。CleanUp 的目的就是"清理"这些数据结构，使其精确且仅表示新的凸包，为下轮迭代做准备。
这项任务并没有想象中那么简单。我们将它分为三个部分，如 Code4.19 所示，分别清理顶点、边和面列表，但处理这三者的顺序很重要。
最容易的是删除那些被标记为 f->visible 的面。要删除的边已经记录在 e->delete 中，
两个相邻的面都是可见的。要删除的顶点比较难搞，因为这些顶点在新凸包上没有关联的边。
void CleanUp(void) {
CleanEdges();
CleanFaces();
CleanVertices();
}
我们先描述 CleanFaces，参见 Code4.20，它直接删除所有被标记为 visible 的面。
从新添加的点可以看到这些面，因此它们都位于新凸包内部。有一点需要注意，通常都是从头开始遍历列表所有元素，在 do-while 循环中再次遇到头时停止。
但是如果，面列表的前两个元素 \(A\) 和 \(B\) 都是 visible 的，就都应该被删除。从 f=faces 开始，
删除元素 f=A 后\(f\) 被置为 \(B\)，宏DELETE 也会修改 faces 指向 B。
如果仍使用标准的循环终止条件 while(f != faces) 可能会提前终止循环。
void CleanFaces( void ) {
tFace f; /* Primary pointer into face list. */
tFace t; /* Temporary pointer, for deleting. */
while ( faces && faces->visible ) {
f = faces;
DELETE( faces, f );
}
f = faces->next;
do {
if ( f->visible ) {
t = f;
f = f->next;
DELETE( faces, t );
}
else f = f->next;
} while ( f != faces );
}
先检查链表的头部是否需要删除，若是则删除之，直到链表头是不可见的或者链表为空，才开始常规循环，可以避开这个问题。
同样的策略也用于 CleanEdges 和 CleanVertices 中的删除操作。
新添加的锥体总是连接到可见区域的边界上。对于这些边界上的每一条边，CleanEdges 将 newface
复制到对应的 adjface 字段中，如 Code4.21 所示。之所以在 CleanFaces 之前调用
CleanEdges，是因为我们需要访问相邻面的 visible 字段来决定覆盖哪一个。
因此，旧的面必须存在才能正确地整合新的面。同时，CleanEdges 删除了所有先前被标记为删除的边，就是那些由 AddOne 标记的边。
void CleanEdges( void ) {
tEdge e; /* Primary index into edge list. */
tEdge t; /* Temporary edge pointer. */
/* Integrate the newfaces into the data structure. */
/* Check every edge. */
e = edges;
do {
if ( e->newface ) {
if ( e->adjface[0]->visible )
e->adjface[0] = e->newface;
else e->adjface[1] = e->newface;
e->newface = NULL;
}
e = e->next;
} while ( e != edges );
/* Delete any edges marked for deletion. */
while ( edges && edges->delete ) {
e = edges;
DELETE( edges, e );
}
e = edges->next;
do {
if ( e->delete ) {
t = e;
e = e->next;
DELETE( edges, t );
}
else e = e->next;
} while ( e != edges );
}
要删除的顶点并没有被任何先前调用的程序标记。由于我们先调用了 CleanEdges，就可以做出如下推断: 如果一个顶点没有关联的边，那么它将严格位于可见区域的内部。
因为那些内部边现在都已经被删除了。如 Code4.22 中的函数 CleanVertices，遍历一遍边列表，将每个边的端点都通过字段
v->onhull 标记为在凸包上。然后，在一个循环中删除所有那些已经处理过但不在凸包上的点。最后，重置顶点记录中的各种标志位。
void CleanVertices(void) {
tEdge e;
tVertex v, t;
/* Mark all vertices incident to some undeleted edge as on the hull. */
e = edges;
do {
e->endpts[0]->onhull = e->endpts[1]->onhull = ONHULL;
e = e->next;
} while ( e != edges );
/* Delete all vertices that have been processed but are not on the hull. */
while ( vertices && vertices->mark && !vertices->onhull ) {
v = vertices;
DELETE( vertices, v );
}
v = vertices->next;
do {
if ( v->mark && !v->onhull ) {
t = v;
v = v->next;
DELETE( vertices, t );
}
else v = v->next;
} while ( v != vertices );
/* Reset flags. */
v = vertices;
do {
v->duplicate = NULL;
v->onhull = !ONHULL;
v = v->next;
} while ( v != vertices );
}
至此，代码的介绍就结束了。显然，看似简单直接的算法与实现代码之间还是有很大距离的。我们将在接下来的三个小节中继续讨论一些更现实的问题。
4.3.3. 检查(Checks)
指望这样复杂的程序在首次实现时就能正确运行是不现实的。我已省略了调试打印日志，以免干扰读者，这些日志可以通过命令行标志开启。 代码中未展示的另一部分或许更值得讨论，一致性检查。同样通过命令行标志，我们可以调用一些函数来彻底梳理数据结构，检查在各种已知属性在一切正常时是否成立。 目前使用的一组检查如下:
- 面的方向(Face orientations): 检查每条边的两个端点在相邻于该边的两个面中是否以相反的顺序出现。
- 凸(Convexity): 检查凸包的每个面与凸包的每个顶点是否构成非负体积。
- 欧拉关系(Euler's relations): 检查是否满足 \(F = 2V - 4\) 以及 \(2E = 3V\)。
这些测试在每次迭代后运行。它们非常慢，但如果能通过所有这些检查，就能让人对程序的正确性产生一定的信心。
4.3.4 性能(Performance)
该程序本质上是二次的，但其性能随数据的不同会有很大差异。在此我们给出两种极端情况的数据，立方体内均匀分布的随机点，以及球面附近均匀分布的随机点。 图 4.14 和 4.15 显示了 \(n = 10,000\) 时的示例。立方体中的大多数点最终并不在凸包（hull）上，而球面附近的点中有很大一部分是凸包的一部分。 在图 4.14 中，凸包有 124 个顶点，因此 10,000 个点中有 9,876 个是内部点。图 4.15 中的凸包有 2,356 个顶点，其余 7,644 个点位于距离球面 2% 的球半径的范围内。 球面点是由长度 \(r = 100\) 的随机向量生成的，其端点随后被舍入为整数坐标，这些截断向量的长度中约有四分之三超过 99。
图 4.16 中描述了 \(n\) 取不同值时，两种情况的计算时间，横坐标表示点的数量，纵坐标是计算时间，单位为秒。 图中的曲线是在 Silicon Graphics 133 MHz Indy 工作站上测试得到的。超线性增长在球体曲线中显而易见，而在立方体曲线中几乎难以察觉。
4.3.5 体积溢出(Volume Overflow)
刚才展示的代码中，所有的几何计算都集中在一处: 体积计算(volume computation)。我们坚持使用整数坐标来表示点，以确保计算结算是正确的。 但现在我们不得不面对一个令人不快的现实，即使使用整数运算来计算体积，也不能保证给出正确的结果，因为存在溢出的可能。 在当前大多数的机器上，有符号整数使用 32 位，可以表示从 \(-2^{31} = -2147483648\) 到 \(2^{31} - 1 = 2147483647\) 的数字，大约二十亿，\(\pm 2.1 \times 10^9\)。 当计算(例如加法或乘法)超出这些界限时，C 程序会继续运行而不报错。这与除零不同，整数溢出是不会被检测的。相反，这 32 位仅仅被解释为一个普通的有符号整数， 这通常意味着稍微超过 \(2^{31} - 1\) 的数字会"回绕"变成负整数。
在大多数情况下，这并没有影响，因为我们几乎不会用到十分大的数字。但是体积计算会将三个坐标相乘。我们把方程(\(\ref{f4.6}\)) 完全展开，有
$$ \begin{equation}\label{f4.8} \begin{array}{l} &-b_x c_y d_z + a_x c_y d_z + b_x c_x d_z - a_y c_x d_z - a_x b_y d_z \\ &+a_y b_x d_z + b_y c_z d_x - a_z c_y d_x - b_z c_x d_y + a_z c_x d_y \\ &+a_x b_z d_y - a_z b_x d_y - b_y c_z d_x + a_y c_z d_x + b_z c_y d_x \\ &-a_z c_y d_x - a_y b_z d_x + a_z b_y d_x + a_x b_y c_z - a_y b_x c_z \\ &-a_x b_z c_y + a_z b_x c_y + a_y b_z c_x - a_z b_y c_x. \end{array} \end{equation} $$该计算的通用项是 \(abc\)，其中 \(a, b, c\) 分别是各点的三个坐标之一。
让我们探索一下这个计算的“安全范围”。Because of the many terms, the freedom of compilers to reorganize the computation, and the possible cancellations,of even incorrect calculations, this is not an easy question to answer. 我能够使计算出错的最小示例使用的坐标仅为 \(\pm 512\)。 从由 \((1, 1, 1)\)、\((1, -1, -1)\)、\((-1, 1, -1)\) 和 \((-1, -1, 1)\) 定义的正四面体 \(T\) 开始，它由以原点为中心的立方体的四个顶点构成。按比例缩放常数 \(c\)， 该四面体的体积为 \(16c^3\)。取 \(c = 2^9 = 512\)，体积为 \(2^{3(9)+4} = 2^{31}\)。因此，
$$ \begin{vmatrix} 512 & 512 & 512 & 1 \\ 512 & -512 & -512 & 1 \\ -512 & 512 & -512 & 1 \\ -512 & -512 & 512 & 1 \end{vmatrix} = 2^{31} = 2147483648 $$但是按照方程(\(\ref{f4.8}\)) 计算得到的结果是 \(-2147483648 = -2^{31}\)
如此小的坐标值就会计算出错，这严重限制了代码的实用性。幸运的是，在现代机器上有一种方法，可以在不付出额外努力的情况下扩展计算的安全范围。
由于大多数机器为 double 分配 64 位，其中超过 50 位用于尾数(mantissa)。因此，计算可以用浮点数更准确地执行。
特别是，上述在整数运算中计算错误的例子，当计算使用浮点算术时，会被正确计算。
但是，使用 double 只是将精度问题转移到了别处。例如，四个点 \((3, 0, 0), (0, 3, 0), (0, 0, 3), (1, 1, 1)\) 是共面的。
第四个点是前三个点确定的三角形的质心。将这些点按比例缩放 \(c\)，得到体积的行列式:
当 \(c = 200001 \approx 2 \times 10^5\) 时，方程 (4.10) 中所有变量都为 double 的计算结果约为 \(16^{21}\)。
原因是计算中的一些中间项很大:
Code4.16 并没有计算方程(\(\ref{f4.8}\)) 的体积，而是计算了一个更高效的公式，参见第 1 章方程(1.12)。
这里的效率是以乘法的数量来衡量的，乘法在大多数机器上比加法或减法更耗时。函数 VolumeSign 几乎总是能正确返回，因为 cancellations 使得没有哪一项需要超过 54 位
但是，这仅仅是将"崩溃点"推得更远了一点。仍然存在三个坐标差值的相乘的情况。当 \(c = 80000001 \approx 8 \times 10^6\) 时，按照 double 计算，
得到的体积为 \(-1.16453 \times 10^{27}\) 而不是 0。一些中间结果能达到 \((10^9)^3 = 10^{27}\)。同样，这超过了可以用 54 位尾数精确表示的 \(10^{16}\)。
我不知道 Code4.16 开始产生错误结果的确切坐标值，但我猜测在大多数机器上，坐标值约为 \(10^6\) 时都会给出精确的结果(练习 4.3.6[12])。
关于 VolumeSign 代码的最后一个要点是，它返回的仅仅是体积的符号，而不是体积本身。
对于可见性测试，这已经足够了。更重要的是，将正确的 double 体积转换为 int 会因为类型转换而返回错误的结果。
这里所面临的基础问题，并没有简单的解决方案。这个问题被称为稳健计算robust computation。以下是几种应对策略：
- 报告算术溢出(Report arithmetic overflows)。C++ 允许定义一个数字的类型，当溢出发生时会报告。其他语言也会报告溢出。这并没有扩展代码的范围，但至少用户会知道它失败了。
- 使用更高精度的算术(Use higher precision arithmetic)。现在提供 64 位整数计算的机器，可以将体积计算的范围扩展到更舒适的水平。
- 使用大数(Use bignums)。一些语言，如 LISP 和 Mathematica，使用任意精度算术，通常称为"大数(bignums)"。这个问题在这些语言中消失了， 但是它们通常不是那么容易与其他应用程序结合使用的。最近，也出现了一些针对几何计算的，任意精度表达式包(Yap 1997)。其中LEDA库可能是功能最强并且广泛使用的一个(Mehlhorn & Näher 1995)。
- 使用特殊的行列式代码(Incorporate special determinant code)。精确计算行列式的迫切需求，驱动了大量关于这个主题的研究。 这类算法最近的一个例子是 Clarkson 的方法（1992），它仅用比坐标所需的位数多几位的数来来计算行列式的符号。 该算法专注于获得正确的符号，而不是试图找到行列式的确切值。这避免了需要更多位数的坐标乘法。
第 1 章中的计算面积所面临的所有问题，在体积计算中都出现了，并且更恶劣。
尽管如此，在 AreaSign 函数中使用 double 是有意义的。
代码 4.23 中的 AreaSign 函数在本书随附的代码中广泛使用，只要需要区域符号的地方都会见到它。
不过，这对于练习 1.6.8[5] 中的质心计算来说是不够的。注意这里将整数强制转换为 double，以便乘法具有更多可用的位数。
我们将在第 7.2 节回到这个问题。
int AreaSign( tPointi a, tPointi b, tPointi c ) {
double area2;
area2 = ( b[0] - a[0] ) * (double)( c[1] - a[1] ) -
( c[0] - a[0] ) * (double)( b[1] - a[1] );
/* The area should be an integer. */
if ( area2 > 0.5 ) return 1;
else if ( area2 < -0.5 ) return -1;
else return 0;
}
4.3.6 练习(Exercises)
- [Programming] 探索 chull.c Explore chull.c。学习如何使用 chull.c 及其相关例程。主要程序有三个，chull、sphere 和 cube。
sphere n 输出球体表面附近的 n 个随机点。cube n 输出立方体内部的 n 个随机点。chull 从标准输入读取点并输出它们的凸包。
sphere 或 cube 的输出可以直接通过管道传递给 chull:
sphere 100 | chull。 有关输入输出格式约定和其他相关信息的详细信息，请参阅首条注释。虽然 chull 生成 Postscript 输出，但可以轻松修改以用于其他图形显示。 - [Programming] 测量时间复杂度 Measure time complexity。通过测量 chull 在 sphere 和 cube 生成的随机数据上的执行时间来测量其时间复杂度。
您可以使用 Unix 函数 time。请参阅
man time。确保您只测量 chull 的执行时间，而不是生成点的时间。 将您机器上的执行时间与图 4.16 中所示的执行时间进行比较。 - [Programming] Profile。使用 Unix 性能分析工具 profiling 工具分析 chull 程序运行时间最长的部分。 使用 -p 标志编译然后运行实用程序 prof。请参阅手册页。
- [Programming] 加速 chull Speed up chull。David Dobkin 通过各种改进，在某些情况下将我的代码速度提高了五倍。提出一些改进建议并实现它们。
- [Programming] 拆分体积计算 Distributed volume computation。如果将体积计算视为面积乘以高度，则可以通过计算每个面 \(f\) 的面积法向量 \(a\)， 然后将从面到 \(p\) 的向量与 \(a\) 点积来计算四面体的高度，其中 \(p\) 是添加到凸包的点，从而节省一些时间。实现此更改并查看它能使代码速度提高多少
- 可见区域Visibility region。证明可见区域(从点 \(p\) 可见的 \(Q\) 区域)是连通的(参见练习 4.2.3[3])。 证明可见区域的边界边构成一个简单的环(与图 4.7 中的情况相反)。基于此性质提出代码改进建议。
- 数据结构Criticize data structures。指出你能想到的所有数据结构的弱点，并提出替代方案。
- 一致性检查Consistency checks。是否存在一种情况，数据结构可能不正确，但无法被 4.3.3 节中讨论的一致性检查检测到。设计一个测试方法来捕获这种情况。
- 具有多个顶点的面Faces with many vertices。设计一个允许面具有任意数量顶点的数据结构。
- 不同的点Distinct points。当并非所有输入点都不同时，代码是否有效？
- [Programming] 删除共面点Deleting coplanar points。对凸包数据结构进行后处理，删除不必要的共面点。
- [Programming] 体积范围Volume range。对于一台有 \(L\) 位浮点尾数(mantissas)的机器，找到一个整数 \(m\)，
使得所有在 [-m, +m] 范围内的顶点坐标，函数
VolumeSign的结果都是正确的。 - [Programming] 体积和双精度浮点数Volume and doubles。找到一个示例，使得在你的机器上 VolumeSign 的双精度浮点数计算不正确，并且坐标的绝对值要尽量小。
- [Programming] 代码分析Break the code。找到一组非共面的点，使得 chull.c 的输出不正确，但所有体积计算都正确。通知作者。