Delaunay 2D算法

整个算法的流程如下所示：

Incremental Delaunay Triangulation.    1 Construct an initial triangle, which is large.    2 Randomly generate a point p in the unit square.    3 InsertVertex(p).    4 Repeat step 2 and 3

整个算法的关键在于insertVertex函数

InsertVertex(p)    1 pFace = LocatePoint(p)    2 FaceSplit( pFace )    3 Legalize three edges of the original pFace.

inCircle函数

用于判断P是否位于三角形[vi,vj,vk]的外接圆中
若位于三角形[vi,vj,vk]中,则：

inCircle(pi,pj,pk,p)=det⎛⎝⎜⎜⎜⎜⎜pixpkxpjxpxpiypkypjypyp2ix+p2iyp2kx+p2kyp2jx+p2jyp2x+p2y1111⎞⎠⎟⎟⎟⎟⎟>0

locatePoint函数

首先，三角形的面积计算如下：
Given a triangle [v0,v1,v2] with vi = (xi,yi), the area is given by

S(v0,v1,v2)=12det⎛⎝⎜⎜⎜x0x1x2y0y1y2111⎞⎠⎟⎟⎟

对于一个三角形[v0,v1,v2]，p是同一平面上一点，重心坐标

αk=S(p,vk+1,vk+2)S(v0,v1,v2)

(α0,α1,α2)称为P在三角形[v0,v1,v2]上的重心坐标
可见，每个αk对应三角形的一条有向边[vi,vj]，若αk为负数，则p在[vi,vj]外面

Face * LocatePoint( Point p)    1 Arbitrarily choose the initial face [vi,vj,vk]    2 Compute the barycentric coordinates of p with respect to current    face.    3 If αi,αj,αk are non-negative, then return the current face.    4 Suppose αi is negative, get the face adjacent to the current face    sharing edge [vj,vk], denote as F ˜    5 If F ˜ is empty, return NULL. The point is outside the whole range.    6 Set current face to be F ˜, repeat through step 2

faceSplit函数

将p所在的face分割

edgeSwap函数

legalizeEdge函数

Suppose e = [v0,v1], the vertex against v is v2, compute the    circum circle c through v,v0,v1.    2 If v2 is outside c, then return false.    3 EdgeSwap(e)    4 Recursive call LegalizeEdge(v,[v1,v2]);    5 Recursive call LegalizeEdge(v,[v0,v2]);    6 return true.

运算结果：

调用函数：

1000个随机点：