[ Delaunay Triangulation ] Delaunay 三角剖分

我高中的時候在morris的網站第一次看到這個東西，那個時候就很想試著寫出來看看，但那時候英文很不好所以不會找資料就沒有成功。直到今年的NPSC，聽說台大出了一題找Voronoi Diagram的題目沒有人寫出來，於是我就想挑戰看看，寒假開始之後就認真研讀相關資料，終於在前幾天完成了隨機增量法以及分治法的模板。

先來說說甚麼是三角剖分 Triangulation 吧！假設平面上散布很多不重複的點，現在用一些線段將這些點連成許多的三角形，滿足所有三角形的內部不會有其他點，而且三角形的數量越多越好，這就是一個三角剖分。

Delaunay Triangulation是一種特殊的三角剖分，滿足所有的三角形其外接圓中不會出現其他的點，這樣的性質會使得所有三角形中最小的角盡量大，而且最重要的是它是 Voronoi Diagram 的對偶圖，Voronoi Diagram 直接求的話很容易有浮點誤差而且程式也不好寫，因此從Delaunay Triangulation直接轉換過去會更加方便。

還有一個很重要的點，透過Delaunay Triangulation 的邊集合可以找出平面歐幾里得距離最小生成樹。

這裡我先給出點的資料結構，方便接下來的講解：

	template<typename T>
	struct point{
	T x,y;
	point(){}
	point(const T&x,const T&y):x(x),y(y){}
	point operator+(const point &b)const{
	return point(x+b.x,y+b.y);
	}
	point operator-(const point &b)const{
	return point(x-b.x,y-b.y);
	}
	point operator*(const T &b)const{
	return point(x*b,y*b);
	}
	point operator/(const T &b)const{
	return point(x/b,y/b);
	}
	bool operator==(const point &b)const{
	return x==b.x&&y==b.y;
	}
	T dot(const point &b)const{
	return x*b.x+y*b.y;
	}
	T cross(const point &b)const{
	return x*b.y-y*b.x;
	}
	T abs2()const{//向量長度的平方
	return dot(*this);
	}
	};

view raw point2D.cpp hosted with ❤ by GitHub

一般來說找出 Delaunay Triangulation 的其中一個操作是要判斷一個點是不是在一個三角形的外接圓內，但直接去求外接圓圓心會有浮點誤差，因此有一個作法是將所有點投影到三維的拋物線曲面中，也就是下圖的碗狀曲面，會發現把一個圓投影上去之後，該圓會落在一個平面上，圓內的點會在平面下面，圓外的點會在平面上面，這樣就可以利用一些處理三維幾何的技術就可以避免浮點誤差的判斷點是否在外接圓內了。

範例程式碼中點在圓內輸出1，圓外輸出-1，圓上輸出0：

	template<class T>
	struct point3D{
	T x, y, z;
	point3D(){}
	point3D(const T &x, const T &y, const T &z):
	x(x), y(y), z(z) {}
	point3D(const point<T> &p){
	x=p.x, y=p.y, z=x*x+y*y;
	}
	point3D operator-(const point3D &b)const{
	return point3D(x-b.x, y-b.y, z-b.z);
	}
	T dot(const point3D &b){
	return x*b.x+y*b.y+z*b.z;
	}
	point3D cross(const point3D &b)const{
	return point3D(y*b.z-z*b.y, z*b.x-x*b.z, x*b.y-y*b.x);
	}
	};
	template<class T>
	int inCircle(const point<T> &a, point<T> b, point<T> c, const point<T> &p){
	if((b-a).cross(c-a)<0) swap(b, c);
	point3D<T> a3(a), b3(b), c3(c), p3(p);
	b3=b3-a3, c3=c3-a3, p3=p3-a3;
	T res = p3.dot(b3.cross(c3));
	return res<0 ? 1 : (res>0 ? -1 : 0);
	}

view raw inCircle.cpp hosted with ❤ by GitHub

當然實作時並不會真正的寫一個三維點的資料結構，這樣程式碼會比較長。

Delaunay Triangulation 我所知道的有兩種，分別是隨機增量法以及分治法。隨機增量法再點隨機加入的時候複雜度是$\ord{n\log n}$，分治法則是穩定$\ord{n\log n}$，但我的實作分治法的效能是隨機增量法的三倍快。由於兩種做法都已經有詳細教學了這裡就不用文字敘述。

1. 隨機增量法教學在《Computational Geometry - Algorithms and Applications》書中有提到，該書有中文版
2. 分治法教學請點該網站連結

這裡附上程式碼，順帶一提morris的分治法好像在某些測資會因為三點共線爛掉

隨機增量法：

	template<class T>
	class Delaunay{
	typedef point<T> PT;
	struct face{
	int a, b, c, tag;
	vector<int> DAG;
	face(int a,int b,int c):a(a),b(b),c(c),tag(0){}
	};
	struct edge{
	int v, fid, pre, next;
	edge(int v,int fid,int pre,int next):
	v(v),fid(fid),pre(pre),next(next){}
	};
	vector<PT> S;
	vector<face> F;
	vector<edge> E;
	int onLeft(int a, int b, int p){
	if(a<3&&b<3) return (a+1)%3==b;
	PT ba = S[b]-S[a], pa = S[p]-S[a];
	if(a<3) ba = S[a]*-1, pa = S[p]-S[b];
	if(b<3) ba = S[b];
	if(p<3) pa = S[p];
	T res = ba.cross(pa);
	return res>0 ? 1 : (res<0 ? -1 : 0);
	}
	int inTriangle(int p, int fid){
	int a = E[F[fid].a].v, b = E[F[fid].b].v;
	int c = E[F[fid].c].v, ab = onLeft(a,b,p);
	int bc = onLeft(b,c,p), ca = onLeft(c,a,p);
	if(!ab&&bc*ca>=0) return F[fid].b;
	if(!bc&&ab*ca>=0) return F[fid].c;
	if(!ca&&ab*bc>=0) return F[fid].a;
	return ab==bc&&bc==ca ? F[fid].a : -1;
	}
	int timeTag;
	int pointIn(int p, int fid = 0){
	if(F[fid].tag==timeTag) return -1;
	F[fid].tag = timeTag;
	int eid = inTriangle(p, fid);
	if(eid<0||F[fid].DAG.empty()) return eid;
	for(int v:F[fid].DAG){
	int res = pointIn(p, v);
	if(res>=0) return res;
	}
	return -2; // error!
	}
	T det3(T a11,T a12,T a13,T a21,T a22,T a23,T a31,T a32,T a33){
	return a11*(a22*a33-a32*a23)-a12*(a21*a33-a31*a23)+a13*(a21*a32-a31*a22);
	}
	int inCircle(const PT &a, const PT &b, const PT &c, const PT &p){
	T as = a.abs2(), bs = b.abs2(), cs = c.abs2(), ps = p.abs2();
	T res = a.x * det3(b.y,bs,1,c.y,cs,1,p.y,ps,1)
	-a.y * det3(b.x,bs,1,c.x,cs,1,p.x,ps,1)
	+as * det3(b.x,b.y,1,c.x,c.y,1,p.x,p.y,1)
	-det3(b.x,b.y,bs,c.x,c.y,cs,p.x,p.y,ps);
	return res<0 ? 1 : (res>0 ? -1 : 0);
	}
	void addFlipFace(int pid, int eid){
	F[E[eid].fid].DAG.emplace_back(F.size());
	F[E[eid^1].fid].DAG.emplace_back(F.size());
	F.emplace_back(E[eid].next, E.size(), E[eid^1].pre);
	int a = F.back().a, c = F.back().c;
	E.emplace_back(pid, F.size()-1, a, c);
	E[a].pre = c;
	E[a].next = E[c].pre = F.back().b;
	E[a].fid = E[c].fid = F.size()-1;
	E[c].next = a;
	}
	void legalizeEdge(int pid,int eid){
	if(E[eid^1].fid<0) return;
	int i=E[eid].v, j=E[eid^1].v, k=E[E[eid^1].next].v;
	if(i>2&&j>2&&k<3) return;
	bool notLegal;
	if(i<3) notLegal = onLeft(pid,j,k)==1;
	else if(j<3) notLegal = onLeft(i,pid,k)==1;
	else notLegal = inCircle(S[pid], S[i], S[j], S[k])==1;
	if(notLegal){
	addFlipFace(k, eid);
	addFlipFace(pid, eid^1);
	E[eid].fid = E[eid^1].fid = -3;
	legalizeEdge(pid, E[eid^1].pre);
	legalizeEdge(pid, E[eid^1].next);
	}
	}
	public:
	void init(){
	S = { PT(-1,-1), PT(1,0), PT(0,1) };
	F = { face(0,2,4) };
	E = { edge(1,0,4,2), edge(0,-1,3,5), edge(2,0,0,4),
	edge(1,-1,5,1), edge(0,0,2,0), edge(2,-1,1,3) };
	timeTag = 0;
	}
	void add(const PT& p){
	S.emplace_back(p);
	int eid = (++timeTag, pointIn(S.size()-1));
	static vector<int> fe; fe.clear();
	fe.emplace_back(E[eid].next);
	fe.emplace_back(E[eid].pre);
	if(!onLeft(E[eid^1].v, E[eid].v, S.size()-1)){
	fe.emplace_back(E[eid^1].next);
	fe.emplace_back(E[eid^1].pre);
	E[eid].fid = E[eid^1].fid = -4;
	}else fe.emplace_back(eid);
	int fn = fe.size(), pid = S.size()-1;
	for(int e:fe){
	F[E[e].fid].DAG.emplace_back(F.size());
	E[e].fid = F.size();
	E[e].next = E.size();
	F.emplace_back(e,E.size(),0);
	E.emplace_back(pid,F.size()-1,e,0);
	E.emplace_back(E[e].v,0,0,0);
	}
	for(int i=0,j=fn-1; i<fn; j=i++){
	int last = F[E[fe[j]].fid].b^1;
	F[E[fe[i]].fid].c = last;
	E[fe[i]].pre = last;
	E[E[fe[i]].next].next = last;
	E[last].pre = E[fe[i]].next;
	E[last].next = fe[i];
	E[last].fid = E[fe[i]].fid;
	}
	for(int e:fe) legalizeEdge(pid, e);
	}
	vector<pair<int,int>> getEdge(){
	vector<pair<int,int>> res;
	for(const auto &f:F){
	if(f.DAG.empty()){
	int a = E[f.a].v,b = E[f.b].v,c = E[f.c].v;
	if(a>2&&b>2) res.emplace_back(a-3, b-3);
	if(b>2&&c>2) res.emplace_back(b-3, c-3);
	if(c>2&&a>2) res.emplace_back(c-3, a-3);
	}
	}
	return res;
	}
	};

view raw DelaunayRandom.cpp hosted with ❤ by GitHub

分治法：

	template<class T>
	class Delaunay{
	struct PT:public point<T>{
	int g[2];
	PT(const point<T> &p):
	point<T>(p){ g[0]=g[1]=-1; }
	};
	static bool cmp(const PT &a,const PT &b){
	return a.x<b.x||(a.x==b.x&&a.y<b.y);
	}
	struct edge{
	int v,g[2];
	edge(int v,int g0,int g1):
	v(v){g[0]=g0,g[1]=g1;}
	};
	vector<PT> S;
	vector<edge> E;
	bool convex(int &from,int to,T LR){
	for(int i=0;i<2;++i){
	int c = E[S[from].g[i]].v;
	auto A=S[from]-S[to], B=S[c]-S[to];
	T v = A.cross(B)*LR;
	if(v>0||(v==0&&B.abs2()<A.abs2()))
	return from = c, true;
	}
	return false;
	}
	void addEdge(int v,int g0,int g1){
	E.emplace_back(v,g0,g1);
	E[E.back().g[0]].g[1] = E.size()-1;
	E[E.back().g[1]].g[0] = E.size()-1;
	}
	void climb(int &p, int e, int n, int nl, int nr, int LR){
	for(int i=E[e].g[LR]; (S[nr]-S[nl]).cross(S[E[i].v]-S[n])>0;){
	if(inCircle(S[E[i].v],S[nl],S[nr],S[E[E[i].g[LR]].v])>=0)
	{ p = i; break; }
	for(int j=0;j<4;++j)
	E[E[i^j/2].g[j%2^1]].g[j%2] = E[i^j/2].g[j%2];
	int j=i; i=E[i].g[LR];
	E[j].g[0]=E[j].g[1]=E[j^1].g[0]=E[j^1].g[1]=-1;
	}
	}
	T det3(T a11,T a12,T a13,T a21,T a22,T a23,T a31,T a32,T a33){
	return a11*(a22*a33-a32*a23)-a12*(a21*a33-a31*a23)+a13*(a21*a32-a31*a22);
	}
	int inCircle(const PT &a, const PT &b, const PT &c, const PT &p){
	T as = a.abs2(), bs = b.abs2(), cs = c.abs2(), ps = p.abs2();
	T res = a.x * det3(b.y,bs,1,c.y,cs,1,p.y,ps,1)
	-a.y * det3(b.x,bs,1,c.x,cs,1,p.x,ps,1)
	+as * det3(b.x,b.y,1,c.x,c.y,1,p.x,p.y,1)
	-det3(b.x,b.y,bs,c.x,c.y,cs,p.x,p.y,ps);
	return res<0 ? 1 : (res>0 ? -1 : 0);
	}
	void divide(int l, int r){
	if(l>=r)return;
	if(l+1==r){
	int A=S[l].g[0]=S[l].g[1]=E.size();
	E.emplace_back(r,A,A);
	int B=S[r].g[0]=S[r].g[1]=E.size();
	E.emplace_back(l,B,B);
	return;
	}
	int mid = (l+r)/2;
	divide(l,mid), divide(mid+1, r);
	int nl = mid, nr = mid+1;
	for(;;){
	if(convex(nl,nr,1)) continue;
	if(S[nr].g[0]!=-1&&convex(nr,nl,-1)) continue;
	break;
	}
	addEdge(nr,S[nl].g[0],S[nl].g[1]);
	S[nl].g[1] = E.size()-1;
	if(S[nr].g[0]==-1){
	addEdge(nl,E.size(),E.size());
	S[nr].g[1] = E.size()-1;
	}else addEdge(nl,S[nr].g[0],S[nr].g[1]);
	S[nr].g[0] = E.size()-1;
	int cl = nl, cr = nr;
	for(;;){
	int pl=-1, pr=-1, side;
	climb(pl,E.size()-2,nl,nl,nr,1);
	climb(pr,E.size()-1,nr,nl,nr,0);
	if(pl==-1&&pr==-1) break;
	if(pl==-1||pr==-1) side = pl==-1;
	else side=inCircle(S[E[pl].v],S[nl],S[nr],S[E[pr].v])<=0;
	if(side){
	nr = E[pr].v;
	addEdge(nr,E.size()-2,E[E.size()-2].g[1]);
	addEdge(nl,E[pr^1].g[0],pr^1);
	}else{
	nl = E[pl].v;
	addEdge(nr,pl^1,E[pl^1].g[1]);
	addEdge(nl,E[E.size()-2].g[0],E.size()-2);
	}
	}
	if(cl==nl&&cr==nr) return;//Collinearity
	S[nl].g[0] = E.size()-2;
	S[nr].g[1] = E.size()-1;
	}
	public:
	void solve(const vector<point<T>> &P){
	S.clear(), E.clear();
	for(const auto &p:P) S.emplace_back(p);
	sort(S.begin(),S.end(),cmp);
	divide(0,int(S.size())-1);
	}
	vector<pair<int,int>> getEdge(){
	vector<pair<int,int>> res;
	for(size_t i=0;i<E.size();i+=2)
	if(E[i].g[0]!=-1)
	res.emplace_back(E[i].v,E[i^1].v);
	return res;
	}
	};

view raw DelaunayDivideAndConquer.cpp hosted with ❤ by GitHub