计算机视觉:三角剖分算法
算法解析
参考论文Triangulate得到对应的伪代码
subroutine triangulate
input : vertex list
output : triangle list
initialize the triangle list
determine the supertriangle
add supertriangle vertices to the end of the vertex list
add the supertriangle to the triangle list
for each sample point in the vertex list
initialize the edge buffer
for each triangle currently in the triangle list
calculate the triangle circumcircle center and radius
if the point lies in the triangle circumcircle then
add the three triangle edges to the edge buffer
remove the triangle from the triangle list
endif
endfor
delete all doubly specified edges from the edge buffer
this leaves the edges of the enclosing polygon only
add to the triangle list all triangles formed between the point
and the edges of the enclosing polygon
endfor
remove any triangles from the triangle list that use the supertriangle vertices
remove the supertriangle vertices from the vertex list
end
同时,参考了这篇博客对算法进行了优化,提高了运行速度
input: 顶点列表(vertices) //vertices为外部生成的随机或乱序顶点列表
output:已确定的三角形列表(triangles)
初始化顶点列表
创建索引列表(indices = new Array(vertices.length)) //indices数组中的值为0,1,2,3,......,vertices.length-1
基于vertices中的顶点x坐标对indices进行sort //sort后的indices值顺序为顶点坐标x从小到大排序(也可对y坐标,本例中针对x坐标)
确定超级三角形
将超级三角形保存至未确定三角形列表(temp triangles)
将超级三角形push到triangles列表
遍历基于indices顺序的vertices中每一个点 //基于indices后,则顶点则是由x从小到大出现
初始化边缓存数组(edge buffer)
遍历temp triangles中的每一个三角形
计算该三角形的圆心和半径
如果该点在外接圆的右侧
则该三角形为Delaunay三角形,保存到triangles
并在temp里去除掉
跳过
如果该点在外接圆外(即也不是外接圆右侧)
则该三角形为不确定 //后面会在问题中讨论
跳过
如果该点在外接圆内
则该三角形不为Delaunay三角形
将三边保存至edge buffer
在temp中去除掉该三角形
对edge buffer进行去重
将edge buffer中的边与当前的点进行组合成若干三角形并保存至temp triangles中
将triangles与temp triangles进行合并
除去与超级三角形有关的三角形
end
算法实现
// 使用到的自定义的类
// 二维图像上的一个点
class Point{
public:
double x, y;
Point(double a, double b) {
x = a;
y = b;
}
Point(){
x = 0;
y = 0;
}
bool operator==(const Point& other) {
return x == other.x && y == other.y;
}
};
// 二维图像上的一个点对
struct PointPair {
Point s, t;
PointPair(Point a, Point b){
s.x = a.x;
s.y = a.y;
t.x = b.x;
t.y = b.y;
}
};
// 二维图像上的边,用点对表示
typedef PointPair Edge;
// 二维图像上的一个三角形
struct Triangle {
Point a,b,c,center;
double r;
Triangle(Point p1, Point p2, Point p3, int mode = 0) {
// 计算外接圆半径(利用向量法求外接圆圆心和半径)
if(mode == 1) {
double A = pow(p1.x, 2) + pow(p1.y, 2);
double B = pow(p2.x, 2) + pow(p2.y, 2);
double C = pow(p3.x, 2) + pow(p3.y, 2);
double G = (p3.y-p2.y) * p1.x + (p1.y-p3.y) * p2.x + (p2.y-p1.y) * p3.x;
center.x = ((B-C) * p1.y + (C-A) * p2.y + (A-B) * p3.y) / (2*G);
center.y = ((C-B) * p1.x + (A-C) * p2.x + (B-A) * p3.x) / (2*G);
r = sqrt(pow(p1.x-center.x, 2) + pow(p1.y-center.y, 2));
}
else {
center.x = 0;
center.y = 0;
r = 0;
}
a.x = p1.x;
a.y = p1.y;
b.x = p2.x;
b.y = p2.y;
c.x = p3.x;
c.y = p3.y;
}
};
// 三角剖分算法
void Delaunay() {
vector<Point> vertex;
vector<Triangle> temp_triangles, triangles;
double minX, maxX, minY, maxY;
double maxTriS, maxTriH;
unsigned char green[3] = {0,255,0};
unsigned char red[3] = {255,0,0};
for(int i = 0; i < pointPairList.size(); i++) {
vertex.push_back(pointPairList[i].s);
}
// 初始化一个大三角形,将所有标记点包含起来
sort(vertex.begin(), vertex.end(), [](const Point& a, const Point& b) { return a.y < b.y; });
minY = vertex[0].y;
maxY = vertex[vertex.size()-1].y;
sort(vertex.begin(), vertex.end(), [](const Point& a, const Point& b) { return a.x < b.x; });
minX = vertex[0].x;
maxX = vertex[vertex.size()-1].x;
maxTriS = (maxX - minX)/2;
maxTriH = (maxY - minY)*2;
Point right(maxX+maxTriS+15, maxY+5);
Point left(minX-maxTriS-15, maxY+5);
Point top(minX+maxTriS, maxY-maxTriH);
Triangle maxTriangle(right, left, top, 1);
temp_triangles.push_back(maxTriangle);
triangles.push_back(maxTriangle);
// 三角剖分算法核心代码
for(int i = 0; i < vertex.size(); i++) {
vector<Edge> buff, buff_res;
for(vector<Triangle>::iterator p = temp_triangles.begin(); p != temp_triangles.end(); ){
// method 1: faster
if(getDistance(vertex[i], p->center) > p->r && vertex[i].x - p->center.x > p->r) {
triangles.push_back(*p);
p = temp_triangles.erase(p);
continue;
}
else if(getDistance(vertex[i], p->center) > p->r && vertex[i].x - p->center.x <= p->r) {
p++;
continue;
}
else {
buff.push_back(Edge(p->a, p->b));
buff.push_back(Edge(p->b, p->c));
buff.push_back(Edge(p->a, p->c));
p = temp_triangles.erase(p);
}
// method 2: slower
/*if(getDistance(vertex[i], p->center) < p->r) {
buff.push_back(Edge(p->a, p->b));
buff.push_back(Edge(p->b, p->c));
buff.push_back(Edge(p->a, p->c));
p = temp_triangles.erase(p);
} else {
p++;
}*/
}
int* pos = new int[buff.size()];
memset(pos, 0, sizeof(int) * buff.size());
for(int j = 0; j < buff.size()-1; j++) {
if(pos[j] != 0) continue;
for(int k = j+1; k < buff.size(); k++) {
if((buff[j].s == buff[k].s && buff[j].t == buff[k].t) ||
(buff[j].s == buff[k].t && buff[j].t == buff[k].s)){
pos[k] = 1;
pos[j] = 1;
}
}
}
for(int j = 0; j < buff.size(); j++){
if(pos[j] == 0) {
buff_res.push_back(buff[j]);
}
}
for(int j = 0; j < buff_res.size(); j++) {
Triangle t(buff_res[j].s, buff_res[j].t, vertex[i], 1);
temp_triangles.push_back(t);
}
delete[] pos;
}
// 将上述得到的临时三角形存下
for(int i = 0; i < temp_triangles.size(); i++) {
triangles.push_back(temp_triangles[i]);
}
// 删去所有与外部三角形的端点相关的三角形
for(vector<Triangle>::iterator p = triangles.begin(); p != triangles.end(); p++) {
if((*p).a == right || (*p).a == left || (*p).a == top ||
(*p).b == right || (*p).b == left || (*p).b == top ||
(*p).c == right || (*p).c == left || (*p).c == top
){
continue;
}
else {
res_triangle.push_back((*p));
}
}
}