Delaunay三角剖分算法

点集的三角剖分（Triangulation），对数值分析（比如有限元分析）以及图形学来说，都是极为重要的一项预处理技术。
尤其是Delaunay三角剖分，由于其独特性，关于点集的很多种几何图都和Delaunay三角剖分相关，如Voronoi图，EMST树，Gabriel图等。

Delaunay剖分所具备的优异特性：

1.最接近：以最近的三点形成三角形，且各线段(三角形的边)皆不相交。
2.唯一性：不论从区域何处开始构建，最终都将得到一致的结果。
3.最优性：任意两个相邻三角形形成的凸四边形的对角线如果可以互换的话，那么两个三角形六个内角中最小的角度不会变大。
4.最规则：如果将三角网中的每个三角形的最小角进行升序排列，则Delaunay三角网的排列得到的数值最大。
5.区域性：新增、删除、移动某一个顶点时只会影响临近的三角形。
6.具有凸多边形的外壳：三角网最外层的边界形成一个凸多边形的外壳。

概念及定义
二维实数域(二维平面)上的三角剖分

定义1：假设V是二维实数域上的有限点集，边e是由点集中的点作为端点构成的封闭线段, E为e的集合。
那么该点集V的一个三角剖分T=(V,E)是一个平面图G，该平面图满足条件：
1.除了端点，平面图中的边不包含点集中的任何点。
2.没有相交边。
3.平面图中所有的面都是三角面，且所有三角面的合集就是点集V的凸包。
那什么是Delaunay三角剖分呢?不过是一种特殊的三角剖分罢了。从Delaunay边说起。

Delaunay边
定义2：假设E中的一条边e（两个端点为a,b），e若满足下列条件，则称之为Delaunay边：
存在一个圆经过a,b两点，圆内不含点集V中任何的点，这一特性又称空圆特性。

Delaunay三角剖分
定义3：如果点集V的一个三角剖分T只包含Delaunay边，那么该三角剖分称为Delaunay三角剖分。

定义4：假设T为V的任一三角剖分，则T是V的一个Delaunay三角剖分，当前仅当T中的每个三角形的外接圆的内部不包含V中任何的点。

定义5：V的Voronoi图是由多边形区域的集合（有些区域可能不是闭合的），该区域仅含点集中的一个点v，区域中的任何位置到v的距离都比该位置到点集中其它所有点的距离短。

由Voronoi图和Delaunay三角剖分的关系，可以引出另一个Delaunay三角剖分的定义：
定义6：将Voronoi图相邻区域（共边的区域）中的点连接起来构成的图，称为Delaunay三角剖分。
如下图：

概念部分到此，下面看看怎么求Delaunay三角剖分。

计算Delaunay三角剖分

问题1：计算二维Delaunay三角剖分
问题输入：二维实数域上的点集V
问题输出：Delaunay三角剖分DT = (V, E).

算法

由不同的定义对应有不同的算法。用得较多的是基于定义3或4的算法。
目前常用的算法又分为好几种，被不同的家伙发现。什么扫描线法（Sweepline），随机增量法（Incremental），分治法（Divide and Conquer）

c++实现：（转载）

//adapted by the example of leanring opencv by crazy_007
//adapted by the OpenCV2.0samplescdelaunay.c
//2010-4-22
#include "StdAfx.h"
#include <cv.h>
#include <highgui.h>
void draw_subdiv_edge( IplImage* img, CvSubdiv2DEdge edge, CvScalar color ); //为了查看代码方便，不然需调整调用函数顺序

/* the script demostrates iterative construction of
delaunay triangulation and voronoi tesselation */

CvSubdiv2D* init_delaunay( CvMemStorage* storage,
CvRect rect )
{
CvSubdiv2D* subdiv;

subdiv = cvCreateSubdiv2D( CV_SEQ_KIND_SUBDIV2D, sizeof(*subdiv),
sizeof(CvSubdiv2DPoint),
sizeof(CvQuadEdge2D),
storage );
cvInitSubdivDelaunay2D( subdiv, rect );//矩形确定的边界

return subdiv;
}

//draw subdiv point
void draw_subdiv_point( IplImage* img, CvPoint2D32f fp, CvScalar color )
{
cvCircle( img, cvPoint(cvRound(fp.x), cvRound(fp.y)), 5, color, CV_FILLED, 8, 0 );//画fp为圆心，5为半径的实心圆表示delaunay顶点
}
//use an external point to locate an edge or vertex or step around the edges of a delaunay tirangle
//画出delaunay 顶点
void locate_point( CvSubdiv2D* subdiv, CvPoint2D32f fp, IplImage* img,
CvScalar active_color )
{
CvSubdiv2DEdge e;
CvSubdiv2DEdge e0 = 0;
CvSubdiv2DPoint* p = 0;

cvSubdiv2DLocate( subdiv, fp, &e0, &p );//使用一个外部的点定位边缘或顶点
//该函数填充三角形的边缘和顶点或者填充该点所处在的Voronoi面
if( e0 )
{
e = e0;
do
{
draw_subdiv_edge( img, e, active_color );//调用下面函数：画出红色直线
e = cvSubdiv2DGetEdge(e,CV_NEXT_AROUND_LEFT);//遍历Delaunay图：返回左区域的下一条的边缘
}
while( e != e0 );
}

draw_subdiv_point( img, fp, active_color );//调用上面的函数：实现在该点处画半径为5的圆
}




/************* 分割线 ************************************/
//draw subdiv edge
void draw_subdiv_edge( IplImage* img, CvSubdiv2DEdge edge, CvScalar color )
{
CvSubdiv2DPoint* org_pt;
CvSubdiv2DPoint* dst_pt;
CvPoint2D32f org;
CvPoint2D32f dst;
CvPoint iorg, idst;

org_pt = cvSubdiv2DEdgeOrg(edge);//Delaunay或者Voronoi边缘的原始点
dst_pt = cvSubdiv2DEdgeDst(edge);//Delaunay或者Voronoi边缘的终点

if( org_pt && dst_pt )
{
org = org_pt->pt;
dst = dst_pt->pt;

iorg = cvPoint( cvRound( org.x ), cvRound( org.y ));
idst = cvPoint( cvRound( dst.x ), cvRound( dst.y ));

cvLine( img, iorg, idst, color, 1, CV_AA, 0 );//画红色直线
}
}

//draw subdiv:遍历所有的Delaunay边
void draw_subdiv( IplImage* img, CvSubdiv2D* subdiv,
CvScalar delaunay_color, CvScalar voronoi_color )
{
CvSeqReader reader;//使用cvSeqReader逐步遍历边：获得细分结构
int i, total = subdiv->edges->total;
int elem_size = subdiv->edges->elem_size;

cvStartReadSeq( (CvSeq*)(subdiv->edges), &reader, 0 );//initialize reader of the sequence

for( i = 0; i < total; i++ )//total是边数目
{
CvQuadEdge2D* edge = (CvQuadEdge2D*)(reader.ptr);
//CvQuadEdge2D平面划分中的Quad-edge(四方边缘结构):四个边缘 (e, eRot(红色) 以及它们的逆（绿色）
if( CV_IS_SET_ELEM( edge ))
{
draw_subdiv_edge( img, (CvSubdiv2DEdge)edge + 1, voronoi_color ); //不知如何理解(CvSubdiv2DEdge)edge + 1
//书中P346：voronoi_edge=(CvSubdiv2DEdge)edge + 1
//直接采用数组位移法进行各种边的对应的(即edge+1),
//cvSubdiv2DRotateEdge((CvSubdiv2DEdge)edge,1)=(CvSubdiv2DEdge)edge+1
//参考网址为：http://tech.ddvip.com/2007-12/119897724239781.html
draw_subdiv_edge( img, (CvSubdiv2DEdge)edge, delaunay_color ); //调用上面的子函数
}

CV_NEXT_SEQ_ELEM( elem_size, reader );
}
}



/************* 分割线 ************************************/
//draw the voronoi facet:遍历Voronoi面
void draw_subdiv_facet( IplImage* img, CvSubdiv2DEdge edge )
{
CvSubdiv2DEdge t = edge;
int i, count = 0;
CvPoint* buf = 0;

// count number of edges in facet
do
{
count++;
t = cvSubdiv2DGetEdge( t, CV_NEXT_AROUND_LEFT );//返回左区域的下一条的边缘
} while (t != edge );

buf = (CvPoint*)malloc( count * sizeof(buf[0]));

// gather points
t = edge;
for( i = 0; i < count; i++ )
{
CvSubdiv2DPoint* pt = cvSubdiv2DEdgeOrg( t );//获得边缘的起点
if( !pt ) break;
buf[i] = cvPoint( cvRound(pt->pt.x), cvRound(pt->pt.y));//点记录在buf中
t = cvSubdiv2DGetEdge( t, CV_NEXT_AROUND_LEFT );
}

if( i == count )
{
CvSubdiv2DPoint* pt = cvSubdiv2DEdgeDst( cvSubdiv2DRotateEdge( edge, 1 ));//获得边缘的终点
cvFillConvexPoly( img, buf, count, CV_RGB(rand()&255,rand()&255,rand()&255), CV_AA, 0 );//一次只能画一个多边形，而且只能画凸多边形
cvPolyLine( img, &buf, &count, 1, 1, CV_RGB(0,0,0), 1, CV_AA, 0);//一次调用中绘制多个多边形
draw_subdiv_point( img, pt->pt, CV_RGB(0,0,0));//画圆
}
free( buf );
}
//draw & paint voronoi graph
void paint_voronoi( CvSubdiv2D* subdiv, IplImage* img )
{
CvSeqReader reader;
int i, total = subdiv->edges->total;
int elem_size = subdiv->edges->elem_size;

cvCalcSubdivVoronoi2D( subdiv );//计算Voronoi图表的细胞结构

cvStartReadSeq( (CvSeq*)(subdiv->edges), &reader, 0 );

for( i = 0; i < total; i++ )
{
CvQuadEdge2D* edge = (CvQuadEdge2D*)(reader.ptr);

if( CV_IS_SET_ELEM( edge ))
{
CvSubdiv2DEdge e = (CvSubdiv2DEdge)edge;
// left
draw_subdiv_facet( img, cvSubdiv2DRotateEdge( e, 1 ));//调用上面的子函数

// right
draw_subdiv_facet( img, cvSubdiv2DRotateEdge( e, 3 ));
}

CV_NEXT_SEQ_ELEM( elem_size, reader );
}
}
/************* 分割线 ************************************/

void run(void)
{
char win[] = "source";
int i;
CvRect rect = { 0, 0, 600, 600 };//外界边界矩形大小
CvMemStorage* storage;//为delaunay申请内存空间
CvSubdiv2D* subdiv;//细分
IplImage* img;
CvScalar active_facet_color, delaunay_color, voronoi_color, bkgnd_color;

active_facet_color = CV_RGB( 255, 0, 0 );
delaunay_color = CV_RGB( 0,0,0);
voronoi_color = CV_RGB(0, 180, 0);
bkgnd_color = CV_RGB(255,255,255);

img = cvCreateImage( cvSize(rect.width,rect.height), 8, 3 );//创建白色背景的图像
cvSet( img, bkgnd_color, 0 );

cvNamedWindow( win, 1 );//窗口名字为"source"

storage = cvCreateMemStorage(0);//初始化内存空间
subdiv = init_delaunay( storage, rect );//initialization convenience function for delaunay subdivision：调用子函数确定矩形边界

printf("Delaunay triangulation will be build now interactively. "
"To stop the process, press any key ");

for( i = 0; i < 20; i++ )
{
CvPoint2D32f fp = cvPoint2D32f( (float)(rand()%(rect.width-10)+5),
(float)(rand()%(rect.height-10)+5));//This is our point holder

locate_point( subdiv, fp, img, active_facet_color );//调用函数
cvShowImage( win, img );

if( cvWaitKey( 100 ) >= 0 ) //等待600ms
break;

cvSubdivDelaunay2DInsert( subdiv, fp );//向Delaunay三角测量中插入一个点
cvCalcSubdivVoronoi2D( subdiv );//计算Voronoi图表的细胞结构
cvSet( img, bkgnd_color, 0 ); // 给一个对象全部元素赋值:void cvSet( CvArr* arr, CvScalar value, const CvArr* mask=NULL );
draw_subdiv( img, subdiv, delaunay_color, voronoi_color );//调用子函数：cvSeqReader逐步遍历边来获得细分结构
cvShowImage( win, img );

if( cvWaitKey(100) >= 0 )
break;
}

cvSet( img, bkgnd_color, 0 );
paint_voronoi( subdiv, img );//调用子函数：
cvShowImage( win, img );

cvWaitKey(0);

cvReleaseMemStorage( &storage );
cvReleaseImage(&img);
cvDestroyWindow( win );
}

int main( int argc, char** argv )
{
run();
return 0;
}