如何在 C++ 中正确三角化多边形
How to correctly triangulate a polygon in C++
我正在对一个对象进行三角剖分（最终，我想实现 Delaunay 三角剖分，但即使在边合法化之前三角剖分也不起作用，所以我想首先关注简单的三角剖分）。我在下面包含了相关代码。我正在实施一种类似于 Mark de Berg 等人在 "Computation Geometry: Algorithms and Application Third Edition" 中描述的三角测量方法。给出的伪代码如下（如果需要我会删除它）：
注意：我通过创建一个边界三角形而不是使用 "lexicographically highest point of P" 来修改伪代码，因为我不太确定如何定义 p-1 和 p-2 正如教科书所说定义它们 "symbolically" 而不是定义确切的单位（公平地说，我当然可能只是误解了它想说的内容）。此外，合法化不是我的代码的一部分（目前），因为这对于 Delaunay 三角剖分是必要的，但我想确保简单的三角剖分按预期工作。
问题是我得到了一些三角剖分，例如 ，其中蓝线来自原始多边形。
其中一些线没有被添加，因为它们是我没有在 findSmallest 方法中添加的点 p0、p1 和 p2 的三角形的一部分。然而，如果我也添加这些三角形，我会得到类似这样的结果：（注意 p0、p1 和 p2 超出了图片的范围）。原始多边形中的一些线（这次为绿色）仍未添加到三角测量中。我不确定我哪里出错了。
我希望代码清晰，但我将解释一些 methods/structs 以防万一：
TriPoint
是 Point 结构的子结构。
p0, p1, p2
三个点围绕多边形形成边界三角形。我的想法来自 this post。
contains(Point p)
returns 如果点在三角形内或其中一条边上则为真。
findCommonTriangle(TriPoint *a, TriPoint *b, Triangle *t)
returns 沿边 ab 入射到 t 的三角形。 （我没有使用 Edges 来计算三角剖分，所以我决定以这种方式获得入射三角形）。
isOnTriangle(Point s)
在三角形 abc 上调用，returns 如果点在边 ab 上则为 1，如果点在边 bc 上为 2，如果点在边 cd 上为 3。如果在三角形内，则 returns 0.
三角剖分本身的代码位于下方：
#include <GL\glew.h>
#include <GL\freeglut.h>
#include <iostream>
#include <array>
#include <vector>
#include "predicates.h"
struct Point {
float x, y;
Point() { }
Point(float a, float b) {
x = a;
y = b;
}
};
struct Triangle;
struct Triangulation;
std::vector<Triangulation *> triangulations;
struct TriPoint : Point {
std::vector<Triangle *> triangles;
TriPoint() { };
int index;
TriPoint(Point a) {
x = a.x;
y = a.y;
}
TriPoint(float x, float y) : Point(x, y) {};
void removeTriangle(Triangle *t) {
for (size_t i = 0; i < triangles.size(); i++) {
if (triangles[i] == t) {
triangles.erase(triangles.begin() + i);
}
}
}
void addTriangle(Triangle *t) {
triangles.push_back(t);
}
};
double pointInLine(Point *a, Point *b, Point *p) {
REAL *A, *B, *P;
A = new REAL[2];
B = new REAL[2];
P = new REAL[2];
A[0] = a->x;
A[1] = a->y;
B[0] = b->x;
B[1] = b->y;
P[0] = p->x;
P[1] = p->y;
double orient = orient2d(A, B, P);
delete(A);
delete(B);
delete(P);
return orient;
}
struct Triangle {
TriPoint *a, *b, *c;
std::vector<Triangle *> children;
Triangle() { };
Triangle(TriPoint *x, TriPoint *y, TriPoint *z) {
a = x;
b = y;
c = z;
orientTri();
x->addTriangle(this);
y->addTriangle(this);
z->addTriangle(this);
}
bool hasChildren() {
return children.size() != 0;
}
void draw() {
glBegin(GL_LINE_STRIP);
glVertex2f(a->x, a->y);
glVertex2f(b->x, b->y);
glVertex2f(c->x, c->y);
glVertex2f(a->x, a->y);
glEnd();
}
bool contains(Point s) {
float as_x = s.x - a->x;
float as_y = s.y - a->y;
bool s_ab = (b->x - a->x)*as_y - (b->y - a->y)*as_x > 0;
if ((c->x - a->x)*as_y - (c->y - a->y)*as_x > 0 == s_ab) return false;
if ((c->x - b->x)*(s.y - b->y) - (c->y - b->y)*(s.x - b->x) > 0 != s_ab) return false;
return true;
}
int isOnTriangle(Point p) {
//Return -1 if outside
//Returns 1 if on AB
//Returns 2 if on BC
//Returns 3 if on CA
//Returns 4 if on B
//Returns 5 if on C
//Returns 6 if on A
double res1 = pointInLine(b, a, &p);
double res2 = pointInLine(c, b, &p);
double res3 = pointInLine(a, c, &p);
/*If triangles are counter-clockwise oriented then a point is inside
the triangle if the three 'res' are < 0, at left of each oriented edge
*/
if (res1 > 0 || res2 > 0 || res3 > 0)
return -1; //outside
if (res1 < 0) {
if (res2 < 0) {
if (res3 < 0) {
return 0; //inside
} else { //res3 == 0
return 3; //on edge3
}
} else { //res2 == 0
if (res3 == 0) {
return 5; //is point shared by edge2 and edge3
}
return 2; //on edge2
}
} else { //res1 == 0
if (res2 == 0) {
return 4; //is point shared by edge1 and edge2
} else if (res3 == 0) {
return 6; //is point shared by edge1 and 3
}
return 1; //on edge 1
}
}
TriPoint *getThirdPoint(TriPoint *x, TriPoint *y) {
if (a != x && a != y)
return a;
if (b != x && b != y)
return b;
return c;
}
bool hasPoint(TriPoint *p) {
return a == p || b == p || c == p;
}
void orientTri() {
REAL *A, *B, *C;
A = new REAL[2];
B = new REAL[2];
C = new REAL[2];
A[0] = a->x;
A[1] = a->y;
B[0] = b->x;
B[1] = b->y;
C[0] = c->x;
C[1] = c->y;
double orientation = orient2d(A, B, C);
if (orientation < 0) {
TriPoint *temp = a;
a = b;
b = temp;
}
delete(A);
delete(B);
delete(C);
}
};
struct Poly {
std::vector<Point> points;
bool selected = false;
};
Triangle *findCommonTriangle(TriPoint *a, TriPoint *b, Triangle *t) {
//Returns a triangle shared by a and b incident to t
for (Triangle *aTri : a->triangles) {
for (Triangle *bTri : b->triangles) {
if (aTri == bTri && aTri != t) {
return aTri;
}
}
}
return NULL;
}
struct Triangulation {
std::vector<Point> points;
std::vector<Triangle *> triangles;
float xMin = 9999;
float xMax = 0;
float yMin;
float yMax;
Triangulation() { };
Triangulation(Poly p) {
points = p.points;
sort();
triangulate();
}
void draw() {
for (Triangle *t : triangles) {
t->draw();
}
}
void sort() {
//Sort by y-value in ascending order.
//If y-values are equal, sort by x in ascending order.
for (size_t i = 0; i < points.size() - 1; i++) {
if (points[i].x < xMin) {
xMin = points[i].x;
}
if (points[i].x > xMax) {
xMax = points[i].x;
}
int index = i;
for (size_t j = i; j < points.size(); j++) {
if (points[index].y > points[j].y) {
index = j;
} else if (points[index].y == points[j].y) {
if (points[index].x > points[j].x) {
index = j;
}
}
}
std::swap(points[i], points[index]);
}
yMin = points[0].y;
yMax = points[points.size() - 1].y;
std::random_shuffle(points.begin(), points.end());
}
void triangulate() {
Triangle *root;
float dx = xMax - xMin;
float dy = yMax - yMin;
float deltaMax = std::max(dx, dy);
float midx = (xMin + xMax) / 2.f;
float midy = (yMin + yMax) / 2.f;
TriPoint *p0;
TriPoint *p1;
TriPoint *p2;
p0 = new TriPoint(midx - 2 * deltaMax, midy - deltaMax);
p1 = new TriPoint(midx, midy + 2 * deltaMax);
p2 = new TriPoint(midx + 2 * deltaMax, midy - deltaMax);
p0->index = 0;
p1->index = -1;
p2->index = -2;
root = new Triangle(p0, p1, p2);
for (size_t i = 0; i < points.size(); i++) {
TriPoint *p = new TriPoint(points[i]);
p->index = i + 1;
Triangle *temp = root;
double in;
while (temp->hasChildren()) {
for (size_t j = 0; j < temp->children.size(); j++) {
in = temp->children[j]->isOnTriangle(points[i]);
if (in >= 0) {
temp = temp->children[j];
break;
}
}
}
in = temp->isOnTriangle(points[i]);
if (in > 0 ) { //Boundary
if (in == 1) { //AB
Triangle *other = findCommonTriangle(temp->a, temp->b, temp);
TriPoint *l = other->getThirdPoint(temp->a, temp->b);
l->removeTriangle(other);
temp->a->removeTriangle(other);
temp->b->removeTriangle(other);
temp->a->removeTriangle(temp);
temp->b->removeTriangle(temp);
temp->c->removeTriangle(temp);
Triangle *n1 = new Triangle(temp->a, p, temp->c);
Triangle *n2 = new Triangle(temp->b, temp->c, p);
Triangle *n3 = new Triangle(temp->a, l, p);
Triangle *n4 = new Triangle(temp->b, p, l);
temp->children.push_back(n1);
temp->children.push_back(n2);
other->children.push_back(n3);
other->children.push_back(n4);
} else if (in == 2) { //BC
Triangle *other = findCommonTriangle(temp->b, temp->c, temp);
TriPoint *l = other->getThirdPoint(temp->b, temp->c);
l->removeTriangle(other);
temp->b->removeTriangle(other);
temp->c->removeTriangle(other);
temp->a->removeTriangle(temp);
temp->b->removeTriangle(temp);
temp->c->removeTriangle(temp);
Triangle *n1 = new Triangle(temp->a, p, temp->c);
Triangle *n2 = new Triangle(temp->b, temp->a, p);
Triangle *n3 = new Triangle(temp->c, p, l);
Triangle *n4 = new Triangle(temp->b, l, p);
temp->children.push_back(n1);
temp->children.push_back(n2);
other->children.push_back(n3);
other->children.push_back(n4);
} else if (in == 3) { //CA
Triangle *other = findCommonTriangle(temp->a, temp->c, temp);
TriPoint *l = other->getThirdPoint(temp->a, temp->c);
l->removeTriangle(other);
temp->a->removeTriangle(other);
temp->c->removeTriangle(other);
temp->a->removeTriangle(temp);
temp->b->removeTriangle(temp);
temp->c->removeTriangle(temp);
Triangle *n1 = new Triangle(temp->b, temp->c, p);
Triangle *n2 = new Triangle(temp->a, temp->b, p);
Triangle *n3 = new Triangle(temp->c, l, p);
Triangle *n4 = new Triangle(temp->a, p, l);
temp->children.push_back(n1);
temp->children.push_back(n2);
other->children.push_back(n3);
other->children.push_back(n4);
}
} else { //Point is inside of triangle
Triangle *t1 = new Triangle(temp->a, temp->b, p);
Triangle *t2 = new Triangle(temp->b, temp->c, p);
Triangle *t3 = new Triangle(temp->c, temp->a, p);
temp->a->removeTriangle(temp);
temp->b->removeTriangle(temp);
temp->c->removeTriangle(temp);
temp->children.push_back(t1);
temp->children.push_back(t2);
temp->children.push_back(t3);
}
} //Triangulation done
findSmallest(root, p0, p1, p2);
triangulations.push_back(this);
}
void findSmallest(Triangle *root, TriPoint *p0, TriPoint *p1, TriPoint *p2) {
bool include = true; //Controls drawing triangles with p0, p1, and p2
if (root->hasChildren()) {
for (Triangle *t : root->children) {
findSmallest(t, p0, p1, p2);
}
} else {
int i0 = root->hasPoint(p0);
int i1 = root->hasPoint(p1);
int i2 = root->hasPoint(p2);
if ((!i0 && !i1 && !i2) || include) {
triangles.push_back(root);
}
}
}
};
Poly polygon;
void changeViewPort(int w, int h)
{
glViewport(0, 0, w, h);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glOrtho(0, glutGet(GLUT_WINDOW_WIDTH), 0, glutGet(GLUT_WINDOW_HEIGHT), -1, 1);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glTranslatef(0.375, 0.375, 0.0);
}
void render() {
glClear(GL_COLOR_BUFFER_BIT);
glLineWidth(2.5);
changeViewPort(glutGet(GLUT_WINDOW_WIDTH), glutGet(GLUT_WINDOW_HEIGHT));
glColor3f(0, 0, 1); //Blue
glBegin(GL_LINE_STRIP);
for (size_t j = 0; j < polygon.points.size(); j++) {
std::vector<Point> ps = polygon.points;
Point p1 = ps[j];
glVertex2i(p1.x, p1.y);
}
glVertex2i(polygon.points[0].x, polygon.points[0].y);
glEnd();
glColor3f(1, 0, 1);
for (Triangulation *t : triangulations) {
t->draw();
}
glutSwapBuffers();
}
int main(int argc, char* argv[]) {
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
glutInitWindowSize(800, 600);
glutCreateWindow("Stack Overflow Question");
glutReshapeFunc(changeViewPort);
glutDisplayFunc(render);
exactinit();
polygon.points.push_back(*new Point(300.0f, 300.0f));
polygon.points.push_back(*new Point(300.0f, 400.0f));
polygon.points.push_back(*new Point(400.0f, 400.0f));
polygon.points.push_back(*new Point(400.0f, 300.0f));
Triangulation t = *(new Triangulation(polygon));
glutMainLoop();
return 0;
}
注意：predicates.cpp 和 predicates.h 是使用来自 here 的代码创建的。
您的代码不是最理想的，但现在这并不重要（您正在学习，对吧？。我将专注于三角剖分问题。
已编辑：您在 sort() 初始化 Triangulation 的 yMin 和 yMax 成员，稍后将它们用于 "big enclosing" 三角形.如果您决定不使用 "sort()"，您将使用初始化值。设置一些默认值。
构建三角剖分不需要对点进行排序。你用它只是为了找到BB，太费力了，最后洗牌，又太费力了。
主要问题（在你的第一个post，在你编辑它之前）我看到的是找到一个点是否在里面的方法一个三角形，在它的边界上，或者在它的外面。
Triangle::isOnTriangle() 太可怕了。您计算了几个 crossproduct 和 return '0'（三角形内）是 none 其中的 none 等于 '0'。
你可能会争辩说你提前知道点不在外面，因为你之前用Triangle::contains()测试过，但这个功能也很糟糕，虽然不是那么多。
找到点与线的相对位置的最佳（或者至少是最简单和最常用的）方法是
res = (y2 - y1)*(px - x1) - (x2 - x1)*(py - y1)
如果 {px,py} 位于 {x1,y1} to {x2,y2} 行的右侧，则 res 为正值。如果在左边则为负，如果在直线上则为零。
这里有两件重要的事情：
- a) 交换方向（即线是
{x2,y2} to {x1,y1}）改变 res. 的符号
- b) 由于数值问题，判断
res 是否真的为零并不容易，就像任何其他浮点精度表示一样。
对于a）你必须确保所有的三角形都具有相同的方向（否则前一个表达式会被错误使用）。您可以特别注意传递给三角形 ctor 的点的顺序。或者您可以添加一个 "orientTri" 函数来设置它们。目前你的边界三角形是顺时针顺序的。最常见的顺序（也用于 OpenGL）是逆时针方向；但是你可以选择你喜欢的那个，注意就好了。
对于 b) 浮点数与 '0' 的比较不是一个好主意。在某些场景中，您可以使用 if std::abs(value) < someEpsilon。但特别是三角测量还不够。 课堂范例
几何计算中的鲁棒性问题很好地解释了为什么你的计算必须是"robust"。 Shewchuk Robust Predicates 是一个很好的解决方案。
一旦你解决了这两个问题，"point in triangle"的问题就可以这样处理了：
double pointInLine(line *L, point *p)
{
//returns positive, negative or exactly zero value
return robustCalculus(L, p);
}
int pointInTriangle(triangle *t, point *p)
{
double res1 = pointInLine(t->edge1, p);
double res2 = pointInLine(t->edge2, p);
double res3 = pointInLine(t->edge3, p);
/*If triangles are counter-clockwise oriented then a point is inside
the triangle if the three 'res' are < 0, at left of each oriented edge
*/
if ( res1 > 0 || res2 > 0 || res3 > 0 )
return -1; //outside
if ( res1 < 0 )
if ( res2 < 0 )
if ( res3 < 0 )
return 0; //inside
else //res3 == 0
return 3; //on edge3
else //res2 == 0
{ if ( res3 == 0 )
return 5; //is point shared by edge2 and edge3
return 2; //on edge2
}
else
... test for other edges or points
}
对于剩下的三角剖分过程一些建议：
因为您需要 Delaunay 三角剖分，所以每次添加新点时都必须检查 "inCircle" 条件（没有其他三角形的外接圆包含这个新点）。可以按照书中或我 post 编辑的链接中的说明来完成。同样，您需要 robust 谓词。
打乱点的插入顺序可能会提高定位新点所在三角形的性能。根据用于定位部件的方法，这可能是正确的。您使用三角形的层次结构，因此数据是否排序不影响。
顺便说一句，维护每个三角形 added/removed/changed 的层次结构在 CPU 和 RAM 中是昂贵的。也许您稍后会找到其他方法，当您获得网格划分的经验时。
在没有层次结构的情况下，"walk to point" 方法（google 方法）在随机输入时似乎更快。但是使用缓存（最后构建的三角形）要有效得多。
祝网格化顺利。开始和调试都很难，细节决定成败。
除了Ripi2的回答中已经指出的问题，我想提出以下建议：
1. random_shuffle:
我看到你没有初始化随机生成器（通过在主函数中调用 srand）。这意味着您将始终使用相同的伪随机序列，因此多次执行该程序将导致完全相同的结果。我测试了你的代码，shuffling确实对代码有影响。正确初始化随机生成器后，可以看到三角剖分每次都在变化，产生不同的三角形。
2。比较：
在您的代码中，您通过比较 指针 来比较点和三角形。这意味着，例如，当且仅当它们在内存中是完全相同的结构时，两个点才相等。具有相同坐标的两个点结构将被视为不同的点。我不确定这是否是您想要获得的，所以我建议您考虑一下。
3。围绕多边形的三角形：
除了硬编码值 (20)，我不明白为什么这段代码应该产生有效的三角剖分。您在多边形上创建了几个三角形，但它们 所有 共享三角形外的那 3 个固定点之一。
这是将硬编码参数减少到 2 以适应视口中所有三角形的图片：
相同的代码，不同的点顺序（在用 time(0) 初始化 srand 之后）：
我无法访问该算法的伪代码，但为了清楚起见，我建议编辑您的答案以对其进行简要描述。
祝你实施顺利:)
我正在对一个对象进行三角剖分（最终，我想实现 Delaunay 三角剖分，但即使在边合法化之前三角剖分也不起作用，所以我想首先关注简单的三角剖分）。我在下面包含了相关代码。我正在实施一种类似于 Mark de Berg 等人在 "Computation Geometry: Algorithms and Application Third Edition" 中描述的三角测量方法。给出的伪代码如下（如果需要我会删除它）：
问题是我得到了一些三角剖分，例如 其中一些线没有被添加，因为它们是我没有在 findSmallest 方法中添加的点 p0、p1 和 p2 的三角形的一部分。然而，如果我也添加这些三角形，我会得到类似这样的结果： 我希望代码清晰，但我将解释一些 methods/structs 以防万一： 是 Point 结构的子结构。 三个点围绕多边形形成边界三角形。我的想法来自 this post。 returns 如果点在三角形内或其中一条边上则为真。 returns 沿边 ab 入射到 t 的三角形。 （我没有使用 Edges 来计算三角剖分，所以我决定以这种方式获得入射三角形）。 在三角形 abc 上调用，returns 如果点在边 ab 上则为 1，如果点在边 bc 上为 2，如果点在边 cd 上为 3。如果在三角形内，则 returns 0. 三角剖分本身的代码位于下方： 注意：predicates.cpp 和 predicates.h 是使用来自 here 的代码创建的。TriPoint
p0, p1, p2
contains(Point p)
findCommonTriangle(TriPoint *a, TriPoint *b, Triangle *t)
isOnTriangle(Point s)
#include <GL\glew.h>
#include <GL\freeglut.h>
#include <iostream>
#include <array>
#include <vector>
#include "predicates.h"
struct Point {
float x, y;
Point() { }
Point(float a, float b) {
x = a;
y = b;
}
};
struct Triangle;
struct Triangulation;
std::vector<Triangulation *> triangulations;
struct TriPoint : Point {
std::vector<Triangle *> triangles;
TriPoint() { };
int index;
TriPoint(Point a) {
x = a.x;
y = a.y;
}
TriPoint(float x, float y) : Point(x, y) {};
void removeTriangle(Triangle *t) {
for (size_t i = 0; i < triangles.size(); i++) {
if (triangles[i] == t) {
triangles.erase(triangles.begin() + i);
}
}
}
void addTriangle(Triangle *t) {
triangles.push_back(t);
}
};
double pointInLine(Point *a, Point *b, Point *p) {
REAL *A, *B, *P;
A = new REAL[2];
B = new REAL[2];
P = new REAL[2];
A[0] = a->x;
A[1] = a->y;
B[0] = b->x;
B[1] = b->y;
P[0] = p->x;
P[1] = p->y;
double orient = orient2d(A, B, P);
delete(A);
delete(B);
delete(P);
return orient;
}
struct Triangle {
TriPoint *a, *b, *c;
std::vector<Triangle *> children;
Triangle() { };
Triangle(TriPoint *x, TriPoint *y, TriPoint *z) {
a = x;
b = y;
c = z;
orientTri();
x->addTriangle(this);
y->addTriangle(this);
z->addTriangle(this);
}
bool hasChildren() {
return children.size() != 0;
}
void draw() {
glBegin(GL_LINE_STRIP);
glVertex2f(a->x, a->y);
glVertex2f(b->x, b->y);
glVertex2f(c->x, c->y);
glVertex2f(a->x, a->y);
glEnd();
}
bool contains(Point s) {
float as_x = s.x - a->x;
float as_y = s.y - a->y;
bool s_ab = (b->x - a->x)*as_y - (b->y - a->y)*as_x > 0;
if ((c->x - a->x)*as_y - (c->y - a->y)*as_x > 0 == s_ab) return false;
if ((c->x - b->x)*(s.y - b->y) - (c->y - b->y)*(s.x - b->x) > 0 != s_ab) return false;
return true;
}
int isOnTriangle(Point p) {
//Return -1 if outside
//Returns 1 if on AB
//Returns 2 if on BC
//Returns 3 if on CA
//Returns 4 if on B
//Returns 5 if on C
//Returns 6 if on A
double res1 = pointInLine(b, a, &p);
double res2 = pointInLine(c, b, &p);
double res3 = pointInLine(a, c, &p);
/*If triangles are counter-clockwise oriented then a point is inside
the triangle if the three 'res' are < 0, at left of each oriented edge
*/
if (res1 > 0 || res2 > 0 || res3 > 0)
return -1; //outside
if (res1 < 0) {
if (res2 < 0) {
if (res3 < 0) {
return 0; //inside
} else { //res3 == 0
return 3; //on edge3
}
} else { //res2 == 0
if (res3 == 0) {
return 5; //is point shared by edge2 and edge3
}
return 2; //on edge2
}
} else { //res1 == 0
if (res2 == 0) {
return 4; //is point shared by edge1 and edge2
} else if (res3 == 0) {
return 6; //is point shared by edge1 and 3
}
return 1; //on edge 1
}
}
TriPoint *getThirdPoint(TriPoint *x, TriPoint *y) {
if (a != x && a != y)
return a;
if (b != x && b != y)
return b;
return c;
}
bool hasPoint(TriPoint *p) {
return a == p || b == p || c == p;
}
void orientTri() {
REAL *A, *B, *C;
A = new REAL[2];
B = new REAL[2];
C = new REAL[2];
A[0] = a->x;
A[1] = a->y;
B[0] = b->x;
B[1] = b->y;
C[0] = c->x;
C[1] = c->y;
double orientation = orient2d(A, B, C);
if (orientation < 0) {
TriPoint *temp = a;
a = b;
b = temp;
}
delete(A);
delete(B);
delete(C);
}
};
struct Poly {
std::vector<Point> points;
bool selected = false;
};
Triangle *findCommonTriangle(TriPoint *a, TriPoint *b, Triangle *t) {
//Returns a triangle shared by a and b incident to t
for (Triangle *aTri : a->triangles) {
for (Triangle *bTri : b->triangles) {
if (aTri == bTri && aTri != t) {
return aTri;
}
}
}
return NULL;
}
struct Triangulation {
std::vector<Point> points;
std::vector<Triangle *> triangles;
float xMin = 9999;
float xMax = 0;
float yMin;
float yMax;
Triangulation() { };
Triangulation(Poly p) {
points = p.points;
sort();
triangulate();
}
void draw() {
for (Triangle *t : triangles) {
t->draw();
}
}
void sort() {
//Sort by y-value in ascending order.
//If y-values are equal, sort by x in ascending order.
for (size_t i = 0; i < points.size() - 1; i++) {
if (points[i].x < xMin) {
xMin = points[i].x;
}
if (points[i].x > xMax) {
xMax = points[i].x;
}
int index = i;
for (size_t j = i; j < points.size(); j++) {
if (points[index].y > points[j].y) {
index = j;
} else if (points[index].y == points[j].y) {
if (points[index].x > points[j].x) {
index = j;
}
}
}
std::swap(points[i], points[index]);
}
yMin = points[0].y;
yMax = points[points.size() - 1].y;
std::random_shuffle(points.begin(), points.end());
}
void triangulate() {
Triangle *root;
float dx = xMax - xMin;
float dy = yMax - yMin;
float deltaMax = std::max(dx, dy);
float midx = (xMin + xMax) / 2.f;
float midy = (yMin + yMax) / 2.f;
TriPoint *p0;
TriPoint *p1;
TriPoint *p2;
p0 = new TriPoint(midx - 2 * deltaMax, midy - deltaMax);
p1 = new TriPoint(midx, midy + 2 * deltaMax);
p2 = new TriPoint(midx + 2 * deltaMax, midy - deltaMax);
p0->index = 0;
p1->index = -1;
p2->index = -2;
root = new Triangle(p0, p1, p2);
for (size_t i = 0; i < points.size(); i++) {
TriPoint *p = new TriPoint(points[i]);
p->index = i + 1;
Triangle *temp = root;
double in;
while (temp->hasChildren()) {
for (size_t j = 0; j < temp->children.size(); j++) {
in = temp->children[j]->isOnTriangle(points[i]);
if (in >= 0) {
temp = temp->children[j];
break;
}
}
}
in = temp->isOnTriangle(points[i]);
if (in > 0 ) { //Boundary
if (in == 1) { //AB
Triangle *other = findCommonTriangle(temp->a, temp->b, temp);
TriPoint *l = other->getThirdPoint(temp->a, temp->b);
l->removeTriangle(other);
temp->a->removeTriangle(other);
temp->b->removeTriangle(other);
temp->a->removeTriangle(temp);
temp->b->removeTriangle(temp);
temp->c->removeTriangle(temp);
Triangle *n1 = new Triangle(temp->a, p, temp->c);
Triangle *n2 = new Triangle(temp->b, temp->c, p);
Triangle *n3 = new Triangle(temp->a, l, p);
Triangle *n4 = new Triangle(temp->b, p, l);
temp->children.push_back(n1);
temp->children.push_back(n2);
other->children.push_back(n3);
other->children.push_back(n4);
} else if (in == 2) { //BC
Triangle *other = findCommonTriangle(temp->b, temp->c, temp);
TriPoint *l = other->getThirdPoint(temp->b, temp->c);
l->removeTriangle(other);
temp->b->removeTriangle(other);
temp->c->removeTriangle(other);
temp->a->removeTriangle(temp);
temp->b->removeTriangle(temp);
temp->c->removeTriangle(temp);
Triangle *n1 = new Triangle(temp->a, p, temp->c);
Triangle *n2 = new Triangle(temp->b, temp->a, p);
Triangle *n3 = new Triangle(temp->c, p, l);
Triangle *n4 = new Triangle(temp->b, l, p);
temp->children.push_back(n1);
temp->children.push_back(n2);
other->children.push_back(n3);
other->children.push_back(n4);
} else if (in == 3) { //CA
Triangle *other = findCommonTriangle(temp->a, temp->c, temp);
TriPoint *l = other->getThirdPoint(temp->a, temp->c);
l->removeTriangle(other);
temp->a->removeTriangle(other);
temp->c->removeTriangle(other);
temp->a->removeTriangle(temp);
temp->b->removeTriangle(temp);
temp->c->removeTriangle(temp);
Triangle *n1 = new Triangle(temp->b, temp->c, p);
Triangle *n2 = new Triangle(temp->a, temp->b, p);
Triangle *n3 = new Triangle(temp->c, l, p);
Triangle *n4 = new Triangle(temp->a, p, l);
temp->children.push_back(n1);
temp->children.push_back(n2);
other->children.push_back(n3);
other->children.push_back(n4);
}
} else { //Point is inside of triangle
Triangle *t1 = new Triangle(temp->a, temp->b, p);
Triangle *t2 = new Triangle(temp->b, temp->c, p);
Triangle *t3 = new Triangle(temp->c, temp->a, p);
temp->a->removeTriangle(temp);
temp->b->removeTriangle(temp);
temp->c->removeTriangle(temp);
temp->children.push_back(t1);
temp->children.push_back(t2);
temp->children.push_back(t3);
}
} //Triangulation done
findSmallest(root, p0, p1, p2);
triangulations.push_back(this);
}
void findSmallest(Triangle *root, TriPoint *p0, TriPoint *p1, TriPoint *p2) {
bool include = true; //Controls drawing triangles with p0, p1, and p2
if (root->hasChildren()) {
for (Triangle *t : root->children) {
findSmallest(t, p0, p1, p2);
}
} else {
int i0 = root->hasPoint(p0);
int i1 = root->hasPoint(p1);
int i2 = root->hasPoint(p2);
if ((!i0 && !i1 && !i2) || include) {
triangles.push_back(root);
}
}
}
};
Poly polygon;
void changeViewPort(int w, int h)
{
glViewport(0, 0, w, h);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glOrtho(0, glutGet(GLUT_WINDOW_WIDTH), 0, glutGet(GLUT_WINDOW_HEIGHT), -1, 1);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glTranslatef(0.375, 0.375, 0.0);
}
void render() {
glClear(GL_COLOR_BUFFER_BIT);
glLineWidth(2.5);
changeViewPort(glutGet(GLUT_WINDOW_WIDTH), glutGet(GLUT_WINDOW_HEIGHT));
glColor3f(0, 0, 1); //Blue
glBegin(GL_LINE_STRIP);
for (size_t j = 0; j < polygon.points.size(); j++) {
std::vector<Point> ps = polygon.points;
Point p1 = ps[j];
glVertex2i(p1.x, p1.y);
}
glVertex2i(polygon.points[0].x, polygon.points[0].y);
glEnd();
glColor3f(1, 0, 1);
for (Triangulation *t : triangulations) {
t->draw();
}
glutSwapBuffers();
}
int main(int argc, char* argv[]) {
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
glutInitWindowSize(800, 600);
glutCreateWindow("Stack Overflow Question");
glutReshapeFunc(changeViewPort);
glutDisplayFunc(render);
exactinit();
polygon.points.push_back(*new Point(300.0f, 300.0f));
polygon.points.push_back(*new Point(300.0f, 400.0f));
polygon.points.push_back(*new Point(400.0f, 400.0f));
polygon.points.push_back(*new Point(400.0f, 300.0f));
Triangulation t = *(new Triangulation(polygon));
glutMainLoop();
return 0;
}
您的代码不是最理想的，但现在这并不重要（您正在学习，对吧？。我将专注于三角剖分问题。
已编辑：您在 sort() 初始化 Triangulation 的 yMin 和 yMax 成员，稍后将它们用于 "big enclosing" 三角形.如果您决定不使用 "sort()"，您将使用初始化值。设置一些默认值。
构建三角剖分不需要对点进行排序。你用它只是为了找到BB，太费力了，最后洗牌，又太费力了。
主要问题（在你的第一个post，在你编辑它之前）我看到的是找到一个点是否在里面的方法一个三角形，在它的边界上，或者在它的外面。
Triangle::isOnTriangle() 太可怕了。您计算了几个 crossproduct 和 return '0'（三角形内）是 none 其中的 none 等于 '0'。
你可能会争辩说你提前知道点不在外面，因为你之前用Triangle::contains()测试过，但这个功能也很糟糕，虽然不是那么多。
找到点与线的相对位置的最佳（或者至少是最简单和最常用的）方法是
res = (y2 - y1)*(px - x1) - (x2 - x1)*(py - y1)
{px,py} 位于 {x1,y1} to {x2,y2} 行的右侧，则 res 为正值。如果在左边则为负，如果在直线上则为零。
这里有两件重要的事情：
- a) 交换方向（即线是
{x2,y2} to {x1,y1}）改变res. 的符号
- b) 由于数值问题，判断
res是否真的为零并不容易，就像任何其他浮点精度表示一样。
对于a）你必须确保所有的三角形都具有相同的方向（否则前一个表达式会被错误使用）。您可以特别注意传递给三角形 ctor 的点的顺序。或者您可以添加一个 "orientTri" 函数来设置它们。目前你的边界三角形是顺时针顺序的。最常见的顺序（也用于 OpenGL）是逆时针方向；但是你可以选择你喜欢的那个，注意就好了。
对于 b) 浮点数与 '0' 的比较不是一个好主意。在某些场景中，您可以使用 if std::abs(value) < someEpsilon。但特别是三角测量还不够。 课堂范例
几何计算中的鲁棒性问题很好地解释了为什么你的计算必须是"robust"。 Shewchuk Robust Predicates 是一个很好的解决方案。
一旦你解决了这两个问题，"point in triangle"的问题就可以这样处理了：
double pointInLine(line *L, point *p)
{
//returns positive, negative or exactly zero value
return robustCalculus(L, p);
}
int pointInTriangle(triangle *t, point *p)
{
double res1 = pointInLine(t->edge1, p);
double res2 = pointInLine(t->edge2, p);
double res3 = pointInLine(t->edge3, p);
/*If triangles are counter-clockwise oriented then a point is inside
the triangle if the three 'res' are < 0, at left of each oriented edge
*/
if ( res1 > 0 || res2 > 0 || res3 > 0 )
return -1; //outside
if ( res1 < 0 )
if ( res2 < 0 )
if ( res3 < 0 )
return 0; //inside
else //res3 == 0
return 3; //on edge3
else //res2 == 0
{ if ( res3 == 0 )
return 5; //is point shared by edge2 and edge3
return 2; //on edge2
}
else
... test for other edges or points
}
对于剩下的三角剖分过程一些建议：
因为您需要 Delaunay 三角剖分，所以每次添加新点时都必须检查 "inCircle" 条件（没有其他三角形的外接圆包含这个新点）。可以按照书中或我 post 编辑的链接中的说明来完成。同样，您需要 robust 谓词。
打乱点的插入顺序可能会提高定位新点所在三角形的性能。根据用于定位部件的方法，这可能是正确的。您使用三角形的层次结构，因此数据是否排序不影响。
顺便说一句，维护每个三角形 added/removed/changed 的层次结构在 CPU 和 RAM 中是昂贵的。也许您稍后会找到其他方法，当您获得网格划分的经验时。
在没有层次结构的情况下，"walk to point" 方法（google 方法）在随机输入时似乎更快。但是使用缓存（最后构建的三角形）要有效得多。
祝网格化顺利。开始和调试都很难，细节决定成败。
除了Ripi2的回答中已经指出的问题，我想提出以下建议：
1. random_shuffle:
我看到你没有初始化随机生成器（通过在主函数中调用 srand）。这意味着您将始终使用相同的伪随机序列，因此多次执行该程序将导致完全相同的结果。我测试了你的代码，shuffling确实对代码有影响。正确初始化随机生成器后，可以看到三角剖分每次都在变化，产生不同的三角形。
2。比较：
在您的代码中，您通过比较 指针 来比较点和三角形。这意味着，例如，当且仅当它们在内存中是完全相同的结构时，两个点才相等。具有相同坐标的两个点结构将被视为不同的点。我不确定这是否是您想要获得的，所以我建议您考虑一下。
3。围绕多边形的三角形：
除了硬编码值 (20)，我不明白为什么这段代码应该产生有效的三角剖分。您在多边形上创建了几个三角形，但它们 所有 共享三角形外的那 3 个固定点之一。
这是将硬编码参数减少到 2 以适应视口中所有三角形的图片：
相同的代码，不同的点顺序（在用 time(0) 初始化 srand 之后）：
我无法访问该算法的伪代码，但为了清楚起见，我建议编辑您的答案以对其进行简要描述。
祝你实施顺利:)