Delaunay 三角网格学习

  本文是为《Mastering Opencv...》第七章准备的，他要使用主动外观模型AMM和POSIT算法做一个头部3D姿态估计。图像上的特征点由人工标定，但由于头部运动，比如张嘴，会导致外观形状的扭曲，即特征带点坐标变化，但相对位置几乎不变。因此我们要将这些特征点映射到一个稳定的模型上。我们采用Delaunay三角网格。【As the shape we are looking for might be distorted, such as an open mouth for instance, we are required to map our texture back to a mean shape】

 本文内容将首先介绍Delaunay相关概念，然后分别给出OPENCV实现、c代码实现【改写自网上，结果似乎有问题】

一、Delaunay三角刨分算法简介

1.三角刨分定义

  三角刨分：假设V是二维实数域上的有限点集，边e是由点集合V中的点作为端点构成的封闭线段，E为e的集合。那么，点集V的一个三角刨分T=(V,E)是一个平面图G，该平面图满足条件：

a.除了端点，平面图中的边不包含点集中的任何点

b.没有相交边

c.平面图中所有的面都是三角面 ，且所有三角面的合集是散点集V的凸包【用不严谨的话来讲，给定二维平面上的点集，凸包就是将最外层的点连接起来构成的凸多边型，它能包含点集中所有点的】

2.Delaunay三角刨分定义

 Delaunay三角刨分是一种特殊的三角刨分：

a.Delaunay边：E中的一条边e(两个端点为a,b)，满足条件：存在一个圆经过a,b两点，圆内(在圆内、圆上最多三点共圆)不包含点集V中其他任何点

b.Delaunay刨分：如果点集V的一个三角刨分T只包含Delaunay边，那么该三角刨分称为Delaunay刨分

3.Delaunay刨分算法

主要有3种：

a.Lawson算法，首先建立一个大的三角形或多边形，把所有数据点包围起来 ，向其中插入一点，该点与包含他的三角形三个点相连，形成3个新 三角形，然后对其逐一进行空外接圆检测，同时用Lawson设计的局部优化过程LOP进行优化，即通过交换对角线的方法来保证所形成的三角网为Delaunay三角网。

b,Bowyer-Watson算法，The Bowyer–Watson algorithm is an incremental algorithm. It works by adding points, one at a time, to a valid Delaunay triangulation of a subset of the desired points. After every insertion, any triangles whose circumcircles contain the new point are deleted, leaving a star-shaped polygonal hole which is then re-triangulated using the new point.

具体可查文献：

c.生长法，待会下文会给出，就是执行效率比较慢：

(1)随机选择一个起始点A，然后选择一个离这个点距离最近的点B，构成初始边，加入边表；

(2)在剩余点中选择一个点作为第三个点C，使得角ACB最大，新生成两个边AC和BC加入边表，并把三角形ABC作为第一个三角形加入三角形表中；

(3)从边表中取出一条边DE，边的两个端点是E和D，设其已在三角形DEF中；边DE把平面分成两个半平面，在剩余的离散点中寻找一个离散点G，它与点F不在边DE的同一个半平面中，且角DGE最大，把新边DG和EG加入边表，把新三角形DGE加入三角形表中；

(4)如果剩余的离散点表中还有点，则转至(3)，否则算法结束。

二、具体实现：

测试的三点集数据，存在txt中，第一行分别表示行数和列列数：

1.opencv版本：

// My_Triangulation.cpp : 定义控制台应用程序的入口点。

#include "stdafx.h"

#include iostream 

#include "cv.h"

#include "highgui.h"

#include opencv2/legacy/legacy.hpp 

using namespace std;

int _tmain(int argc, _TCHAR* argv[])

 FILE* in;

 IplImage *dest;

 int rows,cols,i,j,a,b,total,count;

 CvMemStorage* storage;

 CvSubdiv2D* subdiv;

 vector CvPoint points;

 CvSeqReader reader;

 CvPoint buf[3];

 //读入文本数据

 in = fopen("data.txt","r");

 fscanf(in,"%d%d", rows, cols);

 CvRect rect = { 0, 0, 640, 480 };

 storage = cvCreateMemStorage(0);

 subdiv = cvCreateSubdivDelaunay2D(rect,storage);

 count = 0;

 //初始化坐标

 for(i=0;i rows;i++){

 for(j=0;j cols/2;j++){

 fscanf(in,"%d%d", a, 

 points.push_back(cvPoint(a,b));

 //iterate through points inserting them in the subdivision

 for(int i=0;i points.size();i++){

 float x = points.at(i).x;

 float y = points.at(i).y;

 //printf("%f,%f\n",x,y);

 CvPoint2D32f floatingPoint = cvPoint2D32f(x, y);

 cvSubdivDelaunay2DInsert( subdiv, floatingPoint );

 //draw image and iterating through all the triangles

 dest = cvCreateImage(cvSize(640,480),IPL_DEPTH_8U,3);

 total = subdiv- edges- total;

 int elem_size = subdiv- edges- elem_size;

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

 CvNextEdgeType triangleDirections[2] = {CV_NEXT_AROUND_LEFT,CV_NEXT_AROUND_RIGHT};

 for(int tdi = 0;tdi tdi++){

 CvNextEdgeType triangleDirection = triangleDirections[tdi];

 for(i = 0; i total; i++){

 CvQuadEdge2D* edge = (CvQuadEdge2D*)(reader.ptr);

 if(CV_IS_SET_ELEM(edge)){

 CvSubdiv2DEdge t = (CvSubdiv2DEdge)edge;

 int shouldPaint=1;

 for(j=0;j j++){

 CvSubdiv2DPoint* pt = cvSubdiv2DEdgeOrg(t);

 if(!pt) break;

 buf[j] = cvPoint(cvRound(pt- pt.x), cvRound(pt- pt.y));

 t = cvSubdiv2DGetEdge(t, triangleDirection);

 if((pt- pt.x 0)||(pt- pt.x dest- width))

 shouldPaint=0;

 if((pt- pt.y 0)||(pt- pt.y dest- height))

 shouldPaint=0;

 if(shouldPaint){

 cvLine(dest,buf[0],buf[1],CV_RGB(0,255,0),1,8,0);

 cvLine(dest,buf[1],buf[2],CV_RGB(0,255,0),1,8,0);

 cvLine(dest,buf[2],buf[0],CV_RGB(0,255,0),1,8,0);

 printf("%d,%d,%d,%d,%d,%d\n",buf[0].x,buf[0].y,buf[1].x,buf[1].y,buf[2].x,buf[2].y);

 count++;

 CV_NEXT_SEQ_ELEM(elem_size, reader);

 printf("%d",count);

 cvSaveImage("test.jpg",dest);

 cvNamedWindow("三角刨分",0);

 cvShowImage("三角刨分",dest);

 cvWaitKey(0);

 return 0;

实验结果：

2.生长法：

Delaunay三角网有一个特性，每个三角网形成的外接圆都不包含其他参考点。利用这一个性质，我们可以直接构成Delaunay三角网：

(1).建立第一个三角形

a.判断用来建立TIN的总脚点数，小于3则报错退出。

b.第一点的选择：链表的第一节点，命名为Pt1;

c.第二点的选择，命名为Pt2，条件1：非Pt1点；条件2：B.距Pt1最近；

d.第三点的选择，命名为Pt3，条件1：非Pt1，Pt2点；条件2：与Pt1，Pt2点组成的三角形的外接圆内无其他节点；条件3：与Pt1，Pt2组成的三角形中的角∠Pt1Pt3Pt2最大。

e.生成三边，加入边表

 f.生成第一个三角形，组建三角形表
(2).扩展TIN
a.从边表头取一边，要求：该边flag标志为假（只在一个三角形中）

b.从点链表中搜索一点，要求：条件1：与边中的Pixel3在边的异侧；条件2：该点与边组成的三角形的外接圆内无其他点；条件3：满足上面两条件的点中角Pt1Pt3Pt2最大的点为Pt3。

c.判断新生成的边，如果边表中没有，则加入边表尾，设定标志flag为假，如果有，则设定该边flag为真。

d.将生成的三角形假如三角形表

e.设定选中的边的标志flag为真

f.转至a，直至边表边的标志flag全部为真。

代码，主体部分改编自csdn论坛某篇帖子，都是链表操作。。。

// Triangulation_noopencv.cpp : 定义控制台应用程序的入口点。

#include "stdafx.h"

#include "cv.h"

#include "highgui.h"

struct Pixel //散点数据

 double x,y;

 bool flag;

struct List //数据链表

 Pixel *pixel;

 List *next;

struct Line //三角形边

 Pixel *pixel1; //三角形边一端点

 Pixel *pixel2; //三角形边另一端点

 Pixel *pixel3; //三角形边所对顶点

 bool flag;

struct Linelist //三角形边表

 Line *line;

 Linelist *next;

struct Triangle //三角形表

 Line *line1;

 Line *line2;

 Line *line3;

 Triangle *next;


Triangle *CreateDelaunayTIN(List *list);

double calc_dist(double x1,double y1,double x2,double y2);

Pixel circumcenter (Pixel p0 , Pixel p1 , Pixel p2);//求外接圆圆心

bool vaild_triangulaion(List *list,double radius,Pixel center,Pixel *a,Pixel *b ,Pixel *c);

int _tmain(int argc, _TCHAR* argv[])

 FILE* in;

 int rows,cols,a,b,i,j,count;

 List *my_list,*p;

 IplImage *result;

 my_list = (List *)malloc(sizeof(List));

 my_list- pixel = NULL;

 my_list- next = NULL;

 p = my_list;

 count = 0;

 result = cvCreateImage(cvSize(600,480),IPL_DEPTH_8U,3);

 //读入文本数据

 in = fopen("data.txt","r");

 fscanf(in,"%d%d", rows, cols);

 //初始化坐标,并存入List链表中

 for(i=0;i rows;i++){

 for(j=0;j cols/2;j++){

 fscanf(in,"%d%d", a, 

 Pixel* data = (Pixel *)malloc(sizeof(Pixel));

 data- x = a;

 data- y = b;

 if (my_list- pixel == NULL)

 p- pixel = data;

 p- next = NULL;

 }else{

 List *q = (List *)malloc(sizeof(List));

 q- pixel = data;

 q- next = NULL;

 p- next =q;

 p = p- next;

 //显示图像

 CvPoint pt1,pt2,pt3;

 Triangle *head = CreateDelaunayTIN(my_list);

 Triangle *hp = head;

 while (hp) //fuck,每个三角形存三条边，六个顶点，尼玛有病是不？

 pt1.x = hp- line1- pixel1- pt1.y = hp- line1- pixel1- 

 pt2.x = hp- line1- pixel2- pt2.y = hp- line1- pixel2- 

 //寻找第三个点

 if (hp- line2- pixel1- x ==pt1.x hp- line2- pixel1- y ==pt1.y)

 pt3.x = hp- line2- pixel2- 

 pt3.y = hp- line2- pixel2- 

 }else if (hp- line2- pixel2- x ==pt1.x hp- line2- pixel2- y ==pt1.y)

 pt3.x = hp- line2- pixel1- 

 pt3.y = hp- line2- pixel1- 

 }else if (hp- line2- pixel1- x ==pt2.x hp- line2- pixel1- y ==pt2.y)

 pt3.x = hp- line2- pixel2- 

 pt3.y = hp- line2- pixel2- 

 }else

 pt3.x = hp- line2- pixel1- 

 pt3.y = hp- line2- pixel1- 

 printf("%d,%d,%d,%d,%d,%d\n",pt1.x,pt1.y,pt2.x,pt2.y,pt3.x,pt3.y);

 cvLine(result,pt1,pt2,CV_RGB(0,255,0),1,8,0);

 cvLine(result,pt2,pt3,CV_RGB(0,255,0),1,8,0);

 cvLine(result,pt3,pt1,CV_RGB(0,255,0),1,8,0);

 hp = hp- next;

 count++;

 printf("%d",count);

 cvNamedWindow("三角刨分",0);

 cvShowImage("三角刨分",result);

 cvSaveImage("result.jpg",result);

 cvWaitKey(0);

 return 0;

double calc_dist(double x1,double y1,double x2,double y2)

 return (x1-x2)*(x1-x2)+(y1-y2)*(y1-y2);

//外接圆圆心

Pixel circumcenter (Pixel *p0 , Pixel *p1 , Pixel *p2) 

 Pixel ret;

 double a1=p1- x-p0- x,b1 = p1- y - p0- y,c1 = (sqrt(a1) + sqrt(b1)) / 2; 

 double a2=p2- x-p0- x,b2 = p2- y - p0- y,c2 = (sqrt(a2) + sqrt(b2)) / 2; 

 double d = a1 * b2 - a2 * b1; 

 ret.x = p0- x + (c1 * b2 - c2 * b1) / d; 

 ret.y = p0- y + (a1 * c2 - a2 * c1) / d; 

 return ret;

bool vaild_triangulaion(List *list,double radius,Pixel center,Pixel *a,Pixel *b ,Pixel *c)

 bool result = true;

 List *p;

 p=list;

 while (p)

 if (p- pixel!=a p- pixel!=b p- pixel!=c)

 if ((calc_dist(p- pixel- x,p- pixel- y,center.x,center.y) - radius ) = 0)

 //外接圆内有其他点

 result = false;

 p = p- next;

 return result;


 if(calc_dist(pt1- x,pt1- y,node- pixel- x,node- pixel- y) 

 calc_dist(pt1- x,pt1- y,pt2- x,pt2- y)){

 pt2=node- pixel;

 node=node- next;

 node=list- next;

 pt3=NULL;

 while(node!=NULL)

 if(node- pixel==pt1 || node- pixel==pt2){

 node=node- next;

 continue;

 if(pt3==NULL)

 pt3=node- pixel;

 }else //pt3根据，pt1,pt2组成的边，计算a^2+b^2-c^2/(2*a*b)，取最小的pt3

 //余弦定理，让选pt3，使得∠pt1pt3pt2最大

 if((pow(calc_dist(pt1- x,pt1- y,node- pixel- x,node- pixel- y),2)+pow(calc_dist(pt2- x,pt2- y,node- pixel- x,node- pixel- y),2)

 - pow(calc_dist(pt1- x,pt1- y,pt2- x,pt2- y),2))/(2*calc_dist(pt1- x,pt1- y,node- pixel- x,node- pixel- y)*calc_dist(pt2- x,pt2- y,node- pixel- x,node- pixel- y))

 (pow(calc_dist(pt1- x,pt1- y,pt3- x,pt3- y),2)+pow(calc_dist(pt2- x,pt2- y,pt3- x,pt3- y),2)

 -pow(calc_dist(pt1- x,pt1- y,pt2- x,pt2- y),2))/(2*calc_dist(pt1- x,pt1- y,pt3- x,pt3- y)*calc_dist(pt2- x,pt2- y,pt3- x,pt3- y))){


 Pixel temp = circumcenter(pt1,pt2,node- pixel); //求外接圆圆心

 double radius = calc_dist(temp.x,temp.y,pt1- x,pt1- 

 //遍历所有结点，如果pt3与Pt1，Pt2点组成的三角形的外接圆内无其他节点；

 if (vaild_triangulaion(list,radius,temp,pt1,pt2,node- pixel))

 pt3=node- pixel;

 node=node- next; 

 //LineList

 Linelist *linehead,*linenode,*linelast;

 Line *ln1,*ln2,*ln3;

 linenode=new Linelist;

 linenode- line=new Line;

 linenode- line- pixel1=pt1;

 linenode- line- pixel2=pt2;

 linenode- line- pixel3=pt3;

 linenode- line- flag=false;

 linenode- next=NULL; //第一个三角形边

 linehead=linelast=linenode;

 ln1=linenode- line;

 linenode=new Linelist;

 linenode- line=new Line;

 linenode- line- pixel1=pt2;

 linenode- line- pixel2=pt3;

 linenode- line- pixel3=pt1;

 linenode- line- flag=false;

 linenode- next=NULL; //构造第二个三角形边

 linelast- next=linenode;

 linelast=linenode;

 ln2=linenode- line;

 linenode=new Linelist;

 linenode- line=new Line;

 linenode- line- pixel1=pt3;

 linenode- line- pixel2=pt1;

 linenode- line- pixel3=pt2;

 linenode- line- flag=false;

 linenode- next=NULL; //构造第三个三角形边

 linelast- next=linenode;

 linelast=linenode;

 ln3=linenode- line;

 //first Triangle

 Triangle *tglhead,*tglnode,*tgllast;

 tglnode=new Triangle;

 tglnode- line1=ln1;

 tglnode- line2=ln2;

 tglnode- line3=ln3;

 tglnode- next=NULL;

 tglhead=tgllast=tglnode;

 //expend tin;

 Linelist *linetmp,*linetemp;

 List *pixeltmp;

 double x1,y1,x2,y2,x3,y3;

 linetmp = linehead;

 while(linetmp!=NULL)

 if(linetmp- line- flag==true) //从边表中取出一条边，该边只在个三角形中

 linetmp=linetmp- next;

 continue;

 ln1=linetmp- line;

 pt1=linetmp- line- pixel1;

 pt2=linetmp- line- pixel2;

 x1=linetmp- line- pixel1- 

 y1=linetmp- line- pixel1- 

 x2=linetmp- line- pixel2- 

 y2=linetmp- line- pixel2- 

 x3=linetmp- line- pixel3- //该边对面的顶点

 y3=linetmp- line- pixel3- 

 pixeltmp=list;

 pt3=NULL;

 pt3=NULL;

 while(pixeltmp!=NULL)

 if(pixeltmp- pixel==pt1 || pixeltmp- pixel==pt2) //从点表中取出一点

 pixeltmp=pixeltmp- next;

 continue;

 if(((y2-y1)*pixeltmp- pixel- x+(x1-x2)*pixeltmp- pixel- y+(x2*y1-x1*y2))*((y2-y1)*x3+(x1-x2)*y3+(x2*y1-x1*y2)) =0)//边对应顶点的异侧判断

 pixeltmp=pixeltmp- next;

 continue;

 if(pt3==NULL)pt3=pixeltmp- pixel;

 else

 //余弦定理，让选pt3，使得∠pt1pt3pt2最大

 if((pow(calc_dist(pt1- x,pt1- y,pixeltmp- pixel- x,pixeltmp- pixel- y),2)+pow(calc_dist(pt2- x,pt2- y,pixeltmp- pixel- x,pixeltmp- pixel- y),2)-pow(calc_dist(pt1- x,pt1- y,pt2- x,pt2- y),2))/(2*calc_dist(pt1- x,pt1- y,pixeltmp- pixel- x,pixeltmp- pixel- y)*calc_dist(pt2- x,pt2- y,pixeltmp- pixel- x,pixeltmp- pixel- y))

 (pow(calc_dist(pt1- x,pt1- y,pt3- x,pt3- y),2)+pow(calc_dist(pt2- x,pt2- y,pt3- x,pt3- y),2)-pow(calc_dist(pt1- x,pt1- y,pt2- x,pt2- y),2))/(2*calc_dist(pt1- x,pt1- y,pt3- x,pt3- y)*calc_dist(pt2- x,pt2- y,pt3- x,pt3- y)))

 Pixel temp = circumcenter(pt1,pt2,pixeltmp- pixel); //求外接圆圆心

 double radius = calc_dist(temp.x,temp.y,pt1- x,pt1- 

 //遍历所有结点，如果pt3与Pt1，Pt2点组成的三角形的外接圆内无其他节点；

 if (vaild_triangulaion(list,radius,temp,pt1,pt2,pixeltmp- pixel))

 pt3=pixeltmp- pixel;

 pixeltmp=pixeltmp- next;

 if(pt3!=NULL)

 linetemp=linehead;

 flag=false;

 while(linetemp!=NULL) //判断新生成的边

 if((pt1==linetemp- line- pixel1 pt3==linetemp- line- pixel2)

 || (pt3==linetemp- line- pixel1 pt1==linetemp- line- pixel2))

 linetemp- line- flag=true;

 flag=true;

 break;

 linetemp=linetemp- next;

 if(!flag) //在边表末尾插入新的边pt1pt3

 linenode=new Linelist;

 linenode- line=new Line;

 linenode- line- pixel1=pt3;

 linenode- line- pixel2=pt1;

 linenode- line- pixel3=pt2;

 linenode- line- flag=false;

 linenode- next=NULL;

 linelast- next=linenode;

 linelast=linenode;

 ln2=linenode- line;

 linetemp=linehead;

 flag=false;

 while(linetemp!=NULL)

 if((pt2==linetemp- line- pixel1 pt3==linetemp- line- pixel2)

 || (pt3==linetemp- line- pixel1 pt2==linetemp- line- pixel2))

 linetemp- line- flag=true;

 flag=true;

 break;

 linetemp=linetemp- next;

 if(!flag) //在边表末尾插入新的边pt2pt3

 linenode=new Linelist;

 linenode- line=new Line;

 linenode- line- pixel1=pt2;

 linenode- line- pixel2=pt3;

 linenode- line- pixel3=pt1;

 linenode- line- flag=false;

 linenode- next=NULL;

 linelast- next=linenode;

 linelast=linenode;

 ln3=linenode- line;

 tglnode=new Triangle;

 tglnode- line1=ln1;

 tglnode- line2=ln2;

 tglnode- line3=ln3;

 tglnode- next=NULL;

 tgllast- next=tglnode; //在三角形表插入新的三角形

 tgllast=tglnode;

 linetmp- line- flag=true;

 linetmp=linetmp- next;

 TIN=tglhead;

 return TIN;

实验结果：

居然画成这样，我也很无语。。。但是生长法还是很好实现的，具体论文什么的，自己去中国智网搜吧

另外，用opencv做三角刨分的，这个哥们做的很详细：http://blog.csdn.net/raby_gyl/article/details/17409717

百度文库里，还有一种用java实现代码：http://wenku.baidu.com/view/1f0d15147375a417866f8f70.html

经典算法详解（10）图中有多少个三角形 题目：请说出下面图形中包含多少个三角形？请用一个程序完成计算。 C++版本 1 #include 3 using namespace std; 5 const char NO_POINT = 0 7 //任意的一条线 8 const char *map[]...