Delaunay Triangulation算法学习

基于Voronoi分割的Delaunay Triangulation算法的openCV实现

 

Delaunay三角剖分是前苏联数学家 Delaunay在 1934年提出的:对于任意给定的平面点集 ,只存在着唯一的一种三角剖分方法 ,满足所谓的“ 最大 — 最小角 ” 优化准则 ,即所有最小内角之和最大 ,这就是 Delaunay三角剖分。这种剖分方法遵循“ 最小角最大 ” 和“ 空外接圆 ” 准则

“ 最小角最大 ” 准则是在不出现奇异性的情况下,Delaunay三角剖分最小角之和均大于任何非 Delaunay剖分所形成三角形最小角之和 ,三角形的最小内角之和最大 ,从而使得划分的三角形不会出现某个内角过小的情况 ,比较有利于有限元的后续计算。“ 空外接圆 ” 准则是 Delaunay三角剖分中任意三角形的外接圆内不包括其他结点。因此 ,在各种二维三角剖分中 ,只有 Delaunay三角剖分才同时满足全局和局部最优。

“learning openCV”中提到的最简单的一种计算Delaunay算法（原文为chapter9 中的delaunay triangulation）P301，我翻译半天都感觉翻译模糊不清，不知哪位有中文版的把这部分算法给贴出来？

 

我主要参考就是《平面域中的 Delaunay三角算法》中的三角网生长算法，我觉得好像learning openCV中提到的计算delaunay，是用的这个算法，不过我翻译了半天也感觉模糊不清。索性我就用这个算法。也希望大家把中文版的那个算法贴出来好吗？我们一起学习

 

下面是我调试的代码，贴出来大家一起学习，还有那篇论文我也贴出来，是中文的，好理解的一点。

PS:自己调试了代码好久，后来在OpenCV2.0\samples\c\delaunay.c,找到源代码，白白浪费了自己好长时间，我把代码重新整理下放到下面，也希望大家调试代码有问题时，参看人家源代码，帮助很大的。大家也可以进行设置断点进行delaunay学习

 

//adapted by the example of leanring opencv by crazy_007
//adapted by the OpenCV2.0\samples\c\delaunay.c
//2010-4-22

#include <cv.h>
#include <highgui.h>


//initialization conventence function for delaunay subdivsion
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 ); //rect sets the bounds

return subdiv;
}

//draw subdiv point
void draw_subdiv_point( IplImage* img, CvPoint2D32f fp, CvScalar color )
{
cvCircle( img, cvPoint(cvRound(fp.x), cvRound(fp.y)), 3, color, CV_FILLED, 8, 0 );
}

//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);
dst_pt = cvSubdiv2DEdgeDst(edge);

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
void draw_subdiv( IplImage* img, CvSubdiv2D* subdiv,
CvScalar delaunay_color, CvScalar voronoi_color )
{
CvSeqReader  reader;
int i, total = subdiv->edges->total;
int elem_size = subdiv->edges->elem_size;

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

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

if( CV_IS_SET_ELEM( edge ))
{
draw_subdiv_edge( img, (CvSubdiv2DEdge)edge + 1, voronoi_color );
draw_subdiv_edge( img, (CvSubdiv2DEdge)edge, delaunay_color );
}

CV_NEXT_SEQ_ELEM( elem_size, reader );
}
}

//use an external point to locate an edge or vertex or step around the edges of a delaunay tirangle
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 );

if( e0 )
{
e = e0;
do// always s edges--this is a triangulation,after all.
{
draw_subdiv_edge( img, e, active_color );
e = cvSubdiv2DGetEdge(e,CV_NEXT_AROUND_LEFT);
}
while( e != e0 );
}

draw_subdiv_point( img, fp, active_color );
}

//draw the voronoi facet
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));
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 );

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


int main(int argc,char** argv[])
{

//strorage and structure for delaunay subdivsion
CvRect        rect = { 0, 0, 600, 600 };   //Our outer bounding box
CvMemStorage* storage;                     //Storage for the Delaunay subdivsion
CvSubdiv2D*   subdiv;                      //The subdivision itself
int i;
IplImage* img;                   
CvScalar active_facet_color, delaunay_color, voronoi_color, bkgnd_color;


cvNamedWindow("Delaunay",1);

//initialize
storage = cvCreateMemStorage(0);
subdiv = init_delaunay( storage, rect);             
img=cvCreateImage(cvSize(rect.width,rect.height),8,3);
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);


for(i = 0; i < 100; i++ )
{
CvPoint2D32f fp;               //This is our point holder
// However you want to set points
fp = cvPoint2D32f( (float)(rand()%(rect.width-10)+5),(float)(rand()%(rect.height-10)+5));

locate_point(subdiv,fp,img,active_facet_color);
cvShowImage("Delaunay",img);

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

cvSubdivDelaunay2DInsert( subdiv, fp );
cvCalcSubdivVoronoi2D( subdiv );  // Fill out Voronoi data in subdiv
cvSet(img,bkgnd_color,0);
draw_subdiv(img,subdiv,delaunay_color,voronoi_color);
cvShowImage("Delaunay",img);

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

cvSet(img,bkgnd_color,0);
paint_voronoi(subdiv,img);
cvShowImage("Delaunay",img);

cvWaitKey(0);

cvReleaseMemStorage( &storage );
cvReleaseImage(&img);
cvDestroyWindow("Delaunay");

return 1;
}