介绍
想必研究网格细分技术的同学们不会不知道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. */};第一个定义的是顶点，一个二维或者三维数组；第二个定义的是三角形，别看它定义的类型是double型的二维指针，其实它这里采用了4字节对齐技术，这个triangle里面不仅仅存的是组成三角形的三个顶点，它的存储空间布局是这样的：前三个空间存储的是指向它的邻接三角形的指针，接下来的三个空间存储的是指向顶点的指针，如果存在用户输入的边界边数据，那么接下来的三个空间存储的是指向边界边的指针，还可能会包含任意数量的用户定义的浮点型属性选项，和一个可选的面积限制；第三个定义的带有orientation的triangle，这个orientation的值区间为[0, 2]，用来确定这个三角形的顶点org,dest,apex，所表示的边(顶点org到dest的有向边)以及与这条边邻接的三角形。 下面再来介绍几个宏：/* 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)这几个宏和数组的定义可帮助我们确定在有向三角形otri中三个顶点对应的位置以及两个三角形时如何连接在一起的。当orientation为0时，三角形顶点org,dest,apex在数组中的索引位置为：4,5,3；当orientation为1时，索引位置为：5,3,4；当orientation为2时，索引位置为：3,4,5。再如第三段代码宏定义bond，分别在数组索引位置为orientation的地址处存储指向以顶点org到dest的边为邻接的三角形的指针。三角形数据结构
在Jonathan Richard Shewchuk的论文[1]中更一步讲到了有向三角形的设计思想，它这个设计是参考及改进的Leonidas Guibas 和Jorge Stolfi的研究[2]中的四边结构。JRS引入了假三角形（Ghost Triangle）的概念。这里引述一下他论文中的原文和图片：
To preserve the elegance of Guibas and Stolfi’s presentationof the divide-and-conquer algorithm, each triangulationis surrounded with a layer of “ghost” triangles, one triangleper convex hull edge. The ghost triangles are connected toeach other in a ring about a “vertex at infinity” (really justa null pointer). A single edge is represented by two ghosttriangles, 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_DECLARATORSvoid 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); }}这其中三角形的几何信息的操作宏有：setorg, setdest和setapex，拓扑信息的操作宏有：sym, lnext/lprev, onext/oprev, dnext/oprev, rnext/rprev。在论文[2]中对这些操作有比较详细的介绍。这些几何信息的操作函数比较简单，就是设置顶点org,dest和apex，在这里就不做解释了。下面还是重点讲下拓扑信息的操作宏。这几个宏之间的关系式如下：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**;也许这就是一个人无法抗拒的命运，有你、有我、也有他。
