介绍
想必研究网格细分技术的同学们不会不知道Jonathan Richard Shewchuk，凭借作品Triangle获得了2003年数值计算软件威尔金森奖。先拿下他的图来镇下楼（图片出自：点击打开链接，点击打开链接）：
他的作品Triangle用于生成二维网格和Delaunay三角形。
分析
基础概念
先介绍一下他所定义的几个结构体：
/* The vertex data structure. Each vertex is actually an array of REALs. */
/* The number of REALs is unknown until runtime. An integer boundary */
/* marker, and sometimes a pointer to a triangle, is appended after the */
/* REALs. */
typedef REAL *vertex;
/* The triangle data structure. Each triangle contains three pointers to */
/* adjoining triangles, plus three pointers to vertices, plus three */
/* pointers to subsegments (declared below; these pointers are usually */
/* `dummysub'). It may or may not also contain user-defined attributes */
/* and/or a floating-point "area constraint." It may also contain extra */
/* pointers for nodes, when the user asks for high-order elements. */
/* Because the size and structure of a `triangle' is not decided until */
/* runtime, I haven't simply declared the type `triangle' as a struct. */
typedef REAL **triangle; /* Really: typedef triangle *triangle */
/* An oriented triangle: includes a pointer to a triangle and orientation. */
/* The orientation denotes an edge of the triangle. Hence, there are */
/* three possible orientations. By convention, each edge always points */
/* counterclockwise about the corresponding triangle. */
struct otri {
triangle *tri;
int orient; /* Ranges from 0 to 2. */
};
下面再来介绍几个宏：
/* Fast lookup arrays to speed some of the mesh manipulation primitives. */
int plus1mod3[3] = {1, 2, 0};
int minus1mod3[3] = {2, 0, 1};
/* These primitives determine or set the origin, destination, or apex of a */
/* triangle. */
#define org(otri, vertexptr) \
vertexptr = (vertex) (otri).tri[plus1mod3[(otri).orient] + 3]
#define dest(otri, vertexptr) \
vertexptr = (vertex) (otri).tri[minus1mod3[(otri).orient] + 3]
#define apex(otri, vertexptr) \
vertexptr = (vertex) (otri).tri[(otri).orient + 3]
#define setorg(otri, vertexptr) \
(otri).tri[plus1mod3[(otri).orient] + 3] = (triangle) vertexptr
#define setdest(otri, vertexptr) \
(otri).tri[minus1mod3[(otri).orient] + 3] = (triangle) vertexptr
#define setapex(otri, vertexptr) \
(otri).tri[(otri).orient + 3] = (triangle) vertexptr
/* Bond two triangles together. */
#define bond(otri1, otri2) \
(otri1).tri[(otri1).orient] = encode(otri2); \
(otri2).tri[(otri2).orient] = encode(otri1)
三角形数据结构
在Jonathan Richard Shewchuk的论文[1]中更一步讲到了有向三角形的设计思想，它这个设计是参考及改进的Leonidas Guibas 和Jorge Stolfi的研究[2]中的四边结构。JRS引入了假三角形（Ghost Triangle）的概念。这里引述一下他论文中的原文和图片：
To preserve the elegance of Guibas and Stolfi’s presentation of the divide-and-conquer algorithm, each triangulation is surrounded with a layer of “ghost” triangles, one triangle per convex hull edge. The ghost triangles are connected to each other in a ring about a “vertex at infinity” (really just a null pointer). A single edge is represented by two ghost triangles, as illustrated in Figure 2.
每个三角形都包围有一层假三角形，一条凸包边有一个假三角形，这个假三角形中的一个顶点是指向无穷远的点，在Triangle工程中设置这个指向这个无穷远点的指针为NULL。
分治法中三角形的几何信息和拓扑信息的操作
Triangle工程中分治法三角形剖分迭代函数divconqrecurse的源码如下：
/*****************************************************************************/
/* */
/* divconqrecurse() Recursively form a Delaunay triangulation by the */
/* divide-and-conquer method. */
/* */
/* Recursively breaks down the problem into smaller pieces, which are */
/* knitted together by mergehulls(). The base cases (problems of two or */
/* three vertices) are handled specially here. */
/* */
/* On completion, `farleft' and `farright' are bounding triangles such that */
/* the origin of `farleft' is the leftmost vertex (breaking ties by */
/* choosing the highest leftmost vertex), and the destination of */
/* `farright' is the rightmost vertex (breaking ties by choosing the */
/* lowest rightmost vertex). */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void divconqrecurse(struct mesh *m, struct behavior *b, vertex *sortarray,
int vertices, int axis,
struct otri *farleft, struct otri *farright)
#else /* not ANSI_DECLARATORS */
void divconqrecurse(m, b, sortarray, vertices, axis, farleft, farright)
struct mesh *m;
struct behavior *b;
vertex *sortarray;
int vertices;
int axis;
struct otri *farleft;
struct otri *farright;
#endif /* not ANSI_DECLARATORS */
{
struct otri midtri, tri1, tri2, tri3;
struct otri innerleft, innerright;
REAL area;
int divider;
if (b->verbose > 2) {
printf(" Triangulating %d vertices.\n", vertices);
}
if (vertices == 2) {
/* The triangulation of two vertices is an edge. An edge is */
/* represented by two bounding triangles. */
maketriangle(m, b, farleft);
setorg(*farleft, sortarray[0]);
setdest(*farleft, sortarray[1]);
/* The apex is intentionally left NULL. */
maketriangle(m, b, farright);
setorg(*farright, sortarray[1]);
setdest(*farright, sortarray[0]);
/* The apex is intentionally left NULL. */
bond(*farleft, *farright);
lprevself(*farleft);
lnextself(*farright);
bond(*farleft, *farright);
lprevself(*farleft);
lnextself(*farright);
bond(*farleft, *farright);
if (b->verbose > 2) {
printf(" Creating ");
printtriangle(m, b, farleft);
printf(" Creating ");
printtriangle(m, b, farright);
}
/* Ensure that the origin of `farleft' is sortarray[0]. */
lprev(*farright, *farleft);
return;
} else if (vertices == 3) {
/* The triangulation of three vertices is either a triangle (with */
/* three bounding triangles) or two edges (with four bounding */
/* triangles). In either case, four triangles are created. */
maketriangle(m, b, &midtri);
maketriangle(m, b, &tri1);
maketriangle(m, b, &tri2);
maketriangle(m, b, &tri3);
area = counterclockwise(m, b, sortarray[0], sortarray[1], sortarray[2]);
if (area == 0.0) {
/* Three collinear vertices; the triangulation is two edges. */
setorg(midtri, sortarray[0]);
setdest(midtri, sortarray[1]);
setorg(tri1, sortarray[1]);
setdest(tri1, sortarray[0]);
setorg(tri2, sortarray[2]);
setdest(tri2, sortarray[1]);
setorg(tri3, sortarray[1]);
setdest(tri3, sortarray[2]);
/* All apices are intentionally left NULL. */
bond(midtri, tri1);
bond(tri2, tri3);
lnextself(midtri);
lprevself(tri1);
lnextself(tri2);
lprevself(tri3);
bond(midtri, tri3);
bond(tri1, tri2);
lnextself(midtri);
lprevself(tri1);
lnextself(tri2);
lprevself(tri3);
bond(midtri, tri1);
bond(tri2, tri3);
/* Ensure that the origin of `farleft' is sortarray[0]. */
otricopy(tri1, *farleft);
/* Ensure that the destination of `farright' is sortarray[2]. */
otricopy(tri2, *farright);
} else {
/* The three vertices are not collinear; the triangulation is one */
/* triangle, namely `midtri'. */
setorg(midtri, sortarray[0]);
setdest(tri1, sortarray[0]);
setorg(tri3, sortarray[0]);
/* Apices of tri1, tri2, and tri3 are left NULL. */
if (area > 0.0) {
/* The vertices are in counterclockwise order. */
setdest(midtri, sortarray[1]);
setorg(tri1, sortarray[1]);
setdest(tri2, sortarray[1]);
setapex(midtri, sortarray[2]);
setorg(tri2, sortarray[2]);
setdest(tri3, sortarray[2]);
} else {
/* The vertices are in clockwise order. */
setdest(midtri, sortarray[2]);
setorg(tri1, sortarray[2]);
setdest(tri2, sortarray[2]);
setapex(midtri, sortarray[1]);
setorg(tri2, sortarray[1]);
setdest(tri3, sortarray[1]);
}
/* The topology does not depend on how the vertices are ordered. */
bond(midtri, tri1);
lnextself(midtri);
bond(midtri, tri2);
lnextself(midtri);
bond(midtri, tri3);
lprevself(tri1);
lnextself(tri2);
bond(tri1, tri2);
lprevself(tri1);
lprevself(tri3);
bond(tri1, tri3);
lnextself(tri2);
lprevself(tri3);
bond(tri2, tri3);
/* Ensure that the origin of `farleft' is sortarray[0]. */
otricopy(tri1, *farleft);
/* Ensure that the destination of `farright' is sortarray[2]. */
if (area > 0.0) {
otricopy(tri2, *farright);
} else {
lnext(*farleft, *farright);
}
}
if (b->verbose > 2) {
printf(" Creating ");
printtriangle(m, b, &midtri);
printf(" Creating ");
printtriangle(m, b, &tri1);
printf(" Creating ");
printtriangle(m, b, &tri2);
printf(" Creating ");
printtriangle(m, b, &tri3);
}
return;
} else {
/* Split the vertices in half. */
divider = vertices >> 1;
/* Recursively triangulate each half. */
divconqrecurse(m, b, sortarray, divider, 1 - axis, farleft, &innerleft);
divconqrecurse(m, b, &sortarray[divider], vertices - divider, 1 - axis,
&innerright, farright);
if (b->verbose > 1) {
printf(" Joining triangulations with %d and %d vertices.\n", divider,
vertices - divider);
}
/* Merge the two triangulations into one. */
mergehulls(m, b, farleft, &innerleft, &innerright, farright, axis);
}
}
onext = lprev * sym;
opre = sym * lnext;
dnext = sym * lpre;
dpre = lnext * sym;
rnext = sym * lnext * sym;
rpre = sym * lpre * sym;
a: org, b: dest, c: apex;
abc表示三角形，ab表示三角形的邻接边，*表示其他未标识的顶点。
sym(abc) = ba*;
lnext(abc) = bca;
lpre(abc) = cab;
onext(abc) = ac*;
opre(abc) = a*b;
dnext(abc) = *ba;
dpre(abc) = cb*;
rnext(abc) = *a*;
rpre(abc) = b**;
所以分析方法divconqrecurse：（假设sortarray[0]: a, sortarray[1]: b, sortarray[2]:c）
i. 当顶点数为2时
Step1:
farleft: ab* lprev: *ab lprev: b*a
farright: ba* lnext: a*b lnext: *ba
Step2:
farright: *ba
lprev(farright, farleft) ->farleft: a*b
ii. 当顶点数为3且三点共线时：
Step1:
midtri: ab* lnext: b*a lnext: *ab
tri1: ba* lprev: *ba lprev: a*b
tri2: cb* lnext: b*c lnext: *cb
tri3: bc* lprev: *bc lprev: c*b
Step2:
otricopy(tri1, farleft) ->farleft: a*b
otricopy(tri2, farright) ->farright: *cb
iii. 当顶点为3且为逆序时：
Step1:
midtri: abc lnext: bca lnext: cab
tri1: ba* lprev: *ba lprev: a*b
tri2: cb* lnext: b*c lnext: *cb
tri3: ac* lprev: *ac lprev: c*a
Step2:
otricopy(tri1, farleft) ->farleft: a*b
otricopy(tri2, farright) –>farright: *cb
iv. 当顶点为3且顺序时：(保证组装三角形时仍为逆时针)
Step1:
Midtri: acb lnext: cba lnext: bac
Tri1: ca* lprev: *ca lprev: a*c
Tri2: bc* lnext: c*b lnext: *bc
Tri3: ab* lprev: *ab lprev: b*a
Step2:
otricopy(tri1, farleft) ->farleft: a*c
lnext(farleft, farright) ->farright: *ca
symmetric与flip
symmetric是对称的意思，它与flip是一个对立的概念。如何对立呢？我先来讲一下direction和orientation这两个概念。direction可以这样子看，一条有向线段，从它的起始点位置到它重点位置的方向，这个方向是direction；orientation则更直观，表示逆时针顺序或者顺时针顺序。定义一条有direction和orientation的边，symmetric是指与direction相反，orientation相同的边；而flipped是指与orientation相反，direction相同的边。如下所示示意图：
left与right, interior与exterior
l其实代表左面的意思，操作宏rnext中的r代表的是右面的意思。左右面的概念我给出一张图大家就能明白了：
左边的边表示是使用参数化方程来表示的：P = Q(S, T)；右边的边表示是使用隐式方程来表示的：Q(S, N) = 0。知道边的切线方向T，要求得边的法向N，将T逆时针旋转90°即得N。如右图所示，假如存在点P使得Q(S, N) < 0，则点P在边的内侧（interior）；如果存在点P’使得Q(S, N) > 0，则点P’在边的外侧（exterior）；如果存在点P’’使得Q(S, N) = 0，则点P’’在边上。对应的，如左图所示，边的内侧与边的右侧相同。
一个逆时针环绕的三角形，三角形本身是三角形所有边的左面，三角形外部是三角形所有边的右面。
文献：
[1] Jonathan Richard Shewchuk 《Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator》
[2] LEONIDAS GUIBAS and JORGE STOLFI 《Primitives for the Manipulation of General Subdivisions and the Computation of Voronoi Diagrams》
860
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
