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三⾓剖分。
te定义4：假设T为V的任⼀三⾓剖分，则T是V的⼀个Delaunay三⾓剖分，当前仅当T中的每个三⾓形的外接圆的内部不包含V中任何的点。智商怎么提高
定义5：V的Voronoi图是由多边形区域的集合（有些区域可能不是闭合的），该区域仅含点集中的⼀个点v，区域中的任何位置到v的距离都⽐该位置到点集中其它所有点的距离短。
由Voronoi图和Delaunay三⾓剖分的关系，可以引出另⼀个Delaunay三⾓剖分的定义：
定义6：将Voronoi图相邻区域（共边的区域）中的点连接起来构成的图，称为Delaunay三⾓剖分。
如下图：
概念部分到此，下⾯看看怎么求Delaunay三⾓剖分。
piab计算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.0\samples\c\delaunay.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 teslation */
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顶点}
//u 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 )
{
adwinCvSubdiv2DEdge 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* org_pt;
CvSubdiv2DPoint* dst_pt;
CvPoint2D32f org;
stuffy
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 quence
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+1sick什么意思
//参考⽹址为：/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);
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.\n"
"To stop the process, press any key\n\n");
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;
cvSubdivDelaunay2DInrt( 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);
cvReleaMemStorage( &storage );
cvReleaImage(&img);
cvDestroyWindow( win );
}
int main( int argc, char** argv )
{
run();
return 0;
}
