Quad-Edge Data Structure and Library
Quad-Edge 的动机
我们希望能够创造各种类型的几何体,有些简单并且高度对称的多边形体比如柏拉图立体(The Platonic Solids)。
或者是更普遍的阿基米德多面体(Archimedean Solids)
或者甚至是更普遍的均匀多面体( Uniform Polyhedra)
另外有一些非常好的收集在 Polyhedra and Polytopes
所有这些都可以通过truncations(截断每个点),stellations(在每个面上建一个角柱),或者其它在柏拉图固体(Platonic solids)上的操作来实现。编写程序来实现这些多面体将会是一件非常有意思的事情,甚至可以动画化它们之间的转换(3-D morphs)。但是我们要用什么样的数据结构来完成这些事情呢?
First Try:对多边形进行列举。 首先,我们能想到的最简单的数据结构可能是对多边形进行简单的列举,每一个多边形多存储它的点坐标(有冗余!!)。在C++中可以这样写:
struct Vert(double x,y,z); //点的位置
Vert tri[NFACES][3]; //三角形的列表,每一个有三个点
使用这个数据结构,面f的顶点将会是tri[f][i] for i = 0,1,2;上述方案适用于三角模型,其中每个面(多边形)具有三个面,但显然可以将其推广到具有n面的面的模型。使用这个数据结构,进行像顶点截断(vertex trunction)的操作将会是非常笨拙的,因为对于一个给定的点你需要扎到所有的邻接点。为此,你需要在面列表里搜索具有相同坐标的其他顶点。–inefficient and inelegant.
Second Try:点和面列表。 一个更好的替换将会是对点进行分别存储,并且让每个面成为点的指针:
struct Vert{double x,y,z;}; //顶点坐标
Vert vert[NVERTS]; //顶点列表
struct Tri{Vert *p,*q,*r}; //三角形拥有三个顶点指针
Tri tri[NFACES]; //三角形面的列表
同样,这只是三角形面的代码。第二种方法减少了冗余,但是给定点的邻接点将会是非常费时的(O(NFACES)),因为你需要搜索整个面列表。
上面的两种数据结构在存储几何信息上表示点的几何信息(geometric information)比如顶点坐标都很好,但是它们缺少记录点,边和面之间连接性(connectedness)和邻接性(adjacenies)的拓扑信息(topological information)。第一种方法没有存储拓扑信息,第二种方法只存储了从面到点的指针。
我们可以做的更好。为了这样做我们需要存储更多的拓扑信息,使得我们可以在常数时间内找到与给定的顶点/边/面紧邻的顶点/边/面。
什么是Quad-Edge
多边体的一个特别优雅的数据结构是quad-edge数据结构,它是由 Guibas and Stolfi发明的。它并不能表示所有的多面体,只能表示流形多面体
(
m
a
n
i
f
o
l
d
s
)
(manifolds)
(manifolds) (每个点邻域的面在拓扑上和圆圈
(
d
i
s
k
)
(disk)
(disk)等同。边被两个平面共享)。
在quad-edge数据结构中,有点,边和面的类,但是边发挥了主要的作用。边完整的存储了拓扑信息。面和顶点存储的拓扑信息都和边存储的信息冗余。形象地说,边形成了骨骼,而点和面是可选的装饰品,悬挂在边外。顶点拥有大多数的几何信息。
我们现在描述quad-edge的实现。我们将在较高的API层次上进行描述。完整的代码可以在这里找到。
怎样使用Quad-Edge Library
Edge
Edge类表示了一条有向边。给定Edge* e
,你可以直接找到邻接点,面,和边,以及指向相反方向的对称边(“symmetric edge”)。这些操作都很快。因为边是有向的,并且我们经常想象我们是从外面来观察一个物体,所以我们可以说出边的起点终点以及边的左右面。我们总结下面的界面(详细的看cell/edge.hh)。
-
Edge *e - directed edge
-
Vertex *e->Org() - vertex at the origin of e
-
Vertex *e->Dest() - vertex at the destination of e
-
Face *e->Left() - face on the left of e
-
Face *e->Right() - face on the right of e
在下面的函数中,"next"是指相邻面或顶点以逆时针方式(CCW)表示的下一个。 -
Edge *e->Rnext() - next edge around right face,with same right face
-
Edge *e->Lnext() - next edge around left face,with same left face
-
Edge *e->Onext() - next edge around origin,with same origin
-
Edge *e->Dnext() - next edge around dest,with same dest
在下面的函数中,"prev"是指在相邻面或顶点以顺时针方式(CW)表示的下一个。
- Edge *e->Rprev() - prev edge around right face,with same right face
- Edge *e->Lprev() - prev edge around left face,with same left face
- Edge *e->Oprev() - prev edge around origin,with same origin
- Edge *e->Dprev() - prev edge around dest,with same dest
建议花一点时间在图中验证一下,Lnext和Rnext如何使用绕面旋转,而Onext和Dnext如何使用绕顶点旋转。
下面的成员函数为每条边返回一个唯一的整型ID
- unsigned int e->getID() - id# of this edge
在debug的时候,你可能需要打印ID
使用Lnext
,我们可以以逆时针的方式对边estart
左边的面的周围的边进行遍历。
void leftFromEdge(Edge *estart){
Edge *e = estart;
do{
<do something with edge e>
e = e->Lnext();
}while(e!=estart);
}
一个面的边的数量(面的degree或者valence)通常是3或者更大,但是有时在数据结构的构建的过程中,有1个或2个边的面会很有用,这在几何上对应于环或者退化的条形多边形(sliver polygons)。
同样的,以逆时针方式访问边estart
原点周围的边可以被写为:
void orgFromEdge(Edge *estart){
Edge *e = estart;
do{
<do something with edge e>
e = e->Onext();
}while(e!=estart);
}
点的度是三或者更大,但是像面一样,在构建的过程中我们通常需要度是1或者2的点。
因为访问一个面周围的边非常普遍,我们实现了一些迭代器类来简化代码。使用迭代器,leftFromEdge
的替代方式是:
void edgesOfFace(Face *face) {
// visit edges of face in ccw order; edges have face on the left
FaceEdgeIterator faceEdges(face);
Edge *edge;
while ((edge = faceEdges.next()) != 0)
<do something with edge e>
}
//orgFromEdge的一种替换方式是
void edgesOfVertex(Vertex *vert) {
// visit edges of vertex in ccw order; edges have vert as origin
VertexEdgeIterator vertEdges(vert);
Edge *edge;
while ((edge = vertEdges.next()) != 0)
<do something with edge e>
}
Duality
你也许主要到在上面点和面的处理方式非常相似。这不是意外,因为Guibas and Stolfi在设计quad-edge的时候使用了Duality(对偶)。
多面体的对偶是通过将边旋转90度,将顶点替换为面,将面替换为顶点而产生的多面体。新的顶点位置可以视为旧面的质心。 例如,立方体的对偶是八面体,反之亦然。十二面体和二十面体也是彼此的对偶。 四面体和它自身的一个旋转互为对偶。quad-edge数据结构之所以得名,是因为通过将四倍有向边存储在一起,对偶性在较低层次得以建立。
- Edge *e - directed edge
- Edge *e->Sym() - edge pointing opposite to e
- Edge *e->Rot() - dual edge pointing to the left of e
- Edge *e->InvRot - dual edge pointing to the right of e
你也许会用到Sym,也不不会用到Rot或者InvRot。在内部,我们的实现只存储每个边的4个基本信息(origin vertex,left face,Onext and quadedge index),剩下的邻接操作(Dest,Right,Lnext,Rprev,Dnext…)是从Sym和Rot派生的。对Voronoi图和Delaunay三角剖分,对偶也非常有用,但这是另一门课程(计算几何)。
Vertex
顶点处存储的信息存储了一点拓扑信息(指向顶点输出边之一的指针),以及几何信息( ( x , y , z ) 坐 标 (x,y,z)坐标 (x,y,z)坐标)和可选的属性(color,normals等等)。你将会用到的公共接口总结如下(完整代码见cell/vertex.hh)
class Vertex {
Vec3 pos; // (x,y,z) position of this vertex
const void *data; // data pointer for this vertex
Edge *getEdge(); // an outgoing edge of this vertex
unsigned int getID(); // id# of this vertex (unique for this Cell)
};
这里,Vec3是由三个double类型组成的数组,Vec3来自Simple Vector Library(SVL),他的文档在这里。
在结构体对于每一个顶点有一个data指针,它可以被方便扩展。你可以在里面存储任意4bytes信息,或者是指向额外内存的指针。(比如说:颜色,法线和纹理坐标)。
Face
每个面存储一点拓扑信息,一个指向这个面逆时针方向的边之一的指针,另外有一些额外的属性(颜色等)。你将会用到的公共接口总结如下(代码见cell/face.hh):
class Face{
Edge *getEdge();//这个面的一个逆时针方向的边
const void *data; //当前面的data指针
unsigned int getID();//这个面的ID
}
这里的data指针和vertices的相似。
Cell , and Euler Operators
一个cell是单个多边形,包含了一组点,边和面。你将会用到的最多的代码是下面的4个:
class Cell{
Edge *makeVertexEdge(Vertex *v,Face *left,Face *right);
Edge *makeFaceEdge(Face *f,Vertex *org,Vertex *dest);
void killVertexEdge(Edge *e);
void killFaceEdge(Edge *e);
}
它们叫做欧拉算子,因为这些公式保持了genus为0的多面体的欧拉公式
V
−
E
+
F
=
2
V-E+F=2
V−E+F=2,如果拓扑结构在函数调用前是合法的多面体,那么在函数调用后依然是合法的多面体。注意到这些函数更新了拓扑结构,但是它们使用了Vertex和Face的默认的构造器,所以生成的新的点的位置是(0,0,0)–你将需要自己设定它们。
同样,这些并不强制边的线性和面的平面性。拥有度是1或2的点或者面是可以接受的,比如:
给定Cell *c
,函数的调用会做以下的事情:
c->makeVertexEdge(v,left,right)
对vertex v进行切分,产生一个都位于left和right面之间的新的边和点。这个新的边的左边是left面,右边是right面,v在它的origin点,并且向新的点在它的destination。这个新的边会被返回;新的点通过Dest()函数很容易被找到。这个新的点和边将被存储在和Cell C
相关联的集合中。如果left和right不和vertex v相邻,那么就会产生报错信息,并且发生core dump。c->makeFaceEdge(f,org,dest)
是上面函数的对偶。它分割平面f,在org和dest两个点之间创建一个新的边和一个新的面。新的边使用org作为它的origin点并且使用dest作为它的destination,f作为它左边的平面,新产生的平面在它的右边。新生成的边会作为返回值被返回。新生成的面可以很容易的通过返回值的Right()函数进行查找。新产生的边和面被存储在和Cell C
相关联的集合中。如果org和dest不是和f相连的,那么就会产生报错信息,同时发生core dump。c->KillVertexEdge(e)
与makeVertexEdge
相反。它移除边e
和点e->Dest()
。因此,c->killVertexEdge(c->makeVertexEdge(v,left,right)
将不会做任何事情。c->KillFaceEdge(e)
与makeFaceEdge
相反。它移除边e
和面e->Right()
。因此,c->KillFaceEdge(c->makeFaceEdge(f,org,dest))
将不会做任何事情。
不使用这些函数来构建quad-edge数据结构也是可以的,通过使用lower-level的函数可以达成。但是它会是数倍困难,并且容易出错,所以我们并不推荐这么做。
为了debug或者display的目的,你会经常需要遍历一个cell的所有的点,边和面。
遍历Cell *c
的所有的边
CellVertexIterator cellVerts(c);
Vertex *v;
while ((v = cellVerts.next()) != 0)
<do something with vertex v>
遍历Cell *c
的所有的面:
CellFaceIterator cellFaces(c);
Face *f;
while ((f = cellFaces.next()) != 0)
<do something with face f>
因此,使用OpenGL来绘制所有的面,使用随机颜色:
#include <stdlib.h>
double frand() {return (double)rand()/RAND_MAX;}
CellFaceIterator cellFaces(c);
Face *f;
while ((f = cellFaces.next()) != 0) {
glColor3f(frand(), frand(), frand()};
glBegin(GL_POLYGON);
FaceEdgeIterator faceEdges(f);
Edge *edge;
while ((edge = faceEdges.next()) != 0)
glVertex3dv(&edge->Org()->pos[0]);
glEnd();
}
在上面的代码中,edge->Org()
是面的当前边的origin,并且&...->pos[0]
取得它的x坐标的地址。我们可以使用Dest()
来等价代替。步过一个cell的(无向)边将会是更复杂,因为我们已经完成了设置。注意有向边的边数是无向边的两倍。上面的代码对所有的有向边遍历了一次,所以它访问了每个无向边两次。但是对于wireframe drawing和其他的目的你希望这些操作只在无向边上运行一次。不对点进行标记或者使用额外的链表队列的一个比较聪明的方法是访问每个无向边两次,但是利用我们访问顶点时是使用addresses来访问的事实,其中的一个地址会大于另一个:
CellFaceIterator cellFaces(c);
Face *f;
while ((f = cellFaces.next()) != 0) {
// visit each face of cell c
FaceEdgeIterator faceEdges(f);
Edge *edge;
while ((edge = faceEdges.next()) != 0) {
// visit each edge of face f
// if edge passes the following, its Sym will not,
// and vice-versa
if (edge->Org() < edge->Dest())
<do something with edge>
}
}
OBJ File I/O
到目前为止我们看到了怎样修改多面体,但是首先如何创建一个?最简单的方法是从文件中直接读一个。我们有读取和写入OBJ文件格式的代码。这种格式的完整的文档在这里获取,这里有一些更复杂的OBJ格式。我们读取和写入的格式如下:
# comment
v x y z
f v1 v2 ... vn
其中第i条以v开头的行定义了顶点i,x,y,z是浮点数字。每一个f开头的行定义一个n边的面,其中vj以点索引的形式出现。举例来说,下面的代码定义了一个四面体:
v -1 -1 -1
v 1 1 -1
v -1 1 1
v 1 -1 1
f 2 3 4
f 1 4 3
f 1 3 2
f 1 2 4
其中顶点1(v1)在坐标(-1,-1,-1),面1由顶点v2,v3,v4构成。为方便起见,从外面观察时face是逆时针方向的。平面提的OBJ文件在这个目录中。
为了读取一个OBJ文件,使用函数Cell *objReadCell(char *filename)
,写一个文件,使用objWriteCell(Cell *cell,char *filename)
。前者在出错时返回NULL。
A Complete Quad-Edge Program
为了将所有的部分整合到一起,我们现在展示一个程序,它将读入一个cube文件,然后将每条边分成两部分。
Edge *split(Edge *e) {
// split edge e, putting new vertex at midpoint
// get Cell pointer from vertex (Edges don't have one)
Cell *c = e->Org()->getCell();
// split, creating new edge and vertex (sets topology only)
Edge *enew = c->makeVertexEdge(e->Org(), e->Left(), e->Right());
// At this point enew->Dest()==e->Org(),
// and enew->Dest(), the new vertex, is between enew and e.
// You might want to check the defn of makeVertexEdge to
// convince yourself of this.
// position new vertex at midpoint (note use of Vec3::operator+)
enew->Dest()->pos = .5*(enew->Org()->pos + e->Dest()->pos);
return enew; // return new edge
}
void splitAll(Cell *c) {
{
// first, set the splitme bits of all original edges
CellFaceIterator cellFaces(c);
Face *f;
while ((f = cellFaces.next()) != 0) {
// visit each face of cell c
FaceEdgeIterator faceEdges(f);
Edge *edge;
while ((edge = faceEdges.next()) != 0) {
// visit each edge of face f
int splitme = edge->Org() < edge->Dest();
// splitme = 0 or 1
// my Sym's bit will be the complement of mine
edge->data = splitme; // set bit
}
}
}
{
// go through again, splitting marked edges
// need to construct a new iterator, hence the {}'s
CellFaceIterator cellFaces(c);
Face *f;
while ((f = cellFaces.next()) != 0) {
// visit each face of cell c
FaceEdgeIterator faceEdges(f);
Edge *edge;
while ((edge = faceEdges.next()) != 0) {
// visit each edge of face f
// if its "splitme" bit set then split it
if ((int)edge->data) {
Edge *enew = split(edge);
// clear splitme bits on two sub-edges and
// their Syms to avoid recursive splitting
edge->data = 0;
edge->Sym()->data = 0;
enew->data = 0;
enew->Sym()->data = 0;
}
}
}
}
}
void main() {
Cell *c;
c = objReadCell("cube.obj"); // read cube.obj
if (!c) exit(1);
splitAll(c); // split all edges
}