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

	//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->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