Delaunay Triangulation

DelaunayDivideAndConquer.cpp

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;
	}
};

DelaunayRandom.cpp

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;
	}
};

inCircle.cpp

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);
}

point2D.cpp

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);
	}
};