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.

具体可查文献：

• Bowyer, Adrian (1981). "Computing Dirichlet tessellations". Comput. J. 24 (2): 162–166. doi:10.1093/comjnl/24.2.162. edit
• Watson, David F. (1981). "Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes". Comput. J. 24 (2): 167–172.doi:10.1093/comjnl/24.2.167. edit

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,&b);points.push_back(cvPoint(a,b));}}//iterate through points inserting them in the subdivisionfor(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 trianglesdest = 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<2;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<3;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,&b);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->x;pt1.y = hp->line1->pixel1->y;pt2.x = hp->line1->pixel2->x;pt2.y = hp->line1->pixel2->y;//寻找第三个点if (hp->line2->pixel1->x ==pt1.x&&hp->line2->pixel1->y ==pt1.y){pt3.x = hp->line2->pixel2->x;pt3.y = hp->line2->pixel2->y;}else if (hp->line2->pixel2->x ==pt1.x&&hp->line2->pixel2->y ==pt1.y){pt3.x = hp->line2->pixel1->x;pt3.y = hp->line2->pixel1->y;}else if (hp->line2->pixel1->x ==pt2.x&&hp->line2->pixel1->y ==pt2.y){pt3.x = hp->line2->pixel2->x;pt3.y = hp->line2->pixel2->y;}else{pt3.x = hp->line2->pixel1->x;pt3.y = hp->line2->pixel1->y;}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;
}Triangle *CreateDelaunayTIN(List *list)
{//组建第一个三角形List *node;Pixel *pt1,*pt2,*pt3;bool flag;Triangle *TIN;pt1=list->pixel;pt2=list->next->pixel;node=list->next->next;while(node!=NULL)			//找p2点，使得p2与p1组成的边最短{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->y);//遍历所有结点，如果pt3与Pt1，Pt2点组成的三角形的外接圆内无其他节点；if (vaild_triangulaion(list,radius,temp,pt1,pt2,node->pixel)){pt3=node->pixel;}}}node=node->next;	}//LineListLinelist *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 TriangleTriangle *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->x;y1=linetmp->line->pixel1->y;x2=linetmp->line->pixel2->x;y2=linetmp->line->pixel2->y;x3=linetmp->line->pixel3->x;//该边对面的顶点y3=linetmp->line->pixel3->y;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->y);//遍历所有结点，如果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