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

Delaunay 三角网格学习相关推荐

VTK修炼之道46：图形基本操作进阶_三角网格体积、表面积、测地距离、包围盒
1.基本图形操作意义图形处理,比如图形平滑.多分辨率分析.特征提取等都离不开一些基本的图形操作.掌握这些基本的图形操作有助于理解和深入学习图形处理和分析方法. VTK中提供了多种图形的基本操作,其中 ...
bullet物理引擎与OpenGL结合导入3D模型进行碰撞检测以及画三角网格的坑
原文作者:aircraft 原文链接:https://www.cnblogs.com/DOMLX/p/11681069.html 一.初始化世界以及模型 /// 冲突配置包含内存的默认设置,冲突设置. ...
三角网格去噪-DNF_Net论文解读
三角网格去噪-DNF_Net论文解读用三维扫描仪或者深度相机获取的原始三维模型,都不可避免的含有噪声.这些含有噪声的三维模型需要做去噪处理,以便后续应用.三维模型常见的表示形式是点云或者三角网格,对 ...
Matlab论文插图绘制模板第67期—三角网格图(Trimesh)
在之前的文章中,分享了Matlab网格曲面图的绘制模板: 以及一些特殊形式的网格曲面图: 这一次,再来分享一种特殊的网格曲面图:三角网格图. 先来看一下成品效果: 特别提示:Matlab论文插图绘制模 ...
ITK:转换三角网格为二进制图像
ITK:转换三角网格为二进制图像内容提要输出结果输入输出 C++实现代码内容提要将三角形的itk :: Mesh转换为二进制的itk :: Image 输出结果输入输出 C++实现代码
数字几何处理作业1：编程实现三角网格上高斯曲率和平均曲率的计算编程部分
三.编程 1.代码用的是中国科大傅孝明老师的框架:框架下载及配置运行 (1)在哪儿添加代码梳理框架的结构后,在MeshViewerWidget.中添加求解曲率的函数,并在MainViewerWid ...
Siggraph三角网格变形之拉普拉斯变换
三角网格变形一直是CAGD相关领域的重点,刚上研究生的时候,感觉有点神奇.而且一上来导师就给我发了一篇基于格林坐标的自由变形的相关paper,让我看,外文文献,看了n多天,第一次看外文文献,啥也没看懂 ...
双边网格学习一：2021cvpr《Ultra-High-Definition Image Dehazing via Multi-Guided Bilateral Learning》
目录背景方法总结背景众所周知,目前深度学习去雾方法速度方面一直是个痛点,难以达到效果和速度兼容.基于此,本论文设计了快速的去雾网络,网络基于双边网格,能够通过提取一种双边网格的数据结构对原 ...
三角网格模型及基于RBF隐曲面方程求解的曲面重建
资料来源:径向基函数和神经网络技术在逆向工程中的应用研究(博士论文:王宏涛) RBF神经网络模型 RBF神经网络起源于数值分析中多变量插值的RBF方法,1988年Broomhead等人首先将该算法应用 ...

Delaunay 三角网格学习

Delaunay 三角网格学习相关推荐

最新文章

热门文章