我正在进行物体三角测量（最终，我想实现Delaunay三角测量，但在合法化边缘之前，三角测量无法正常工作，因此我想先专注于简单的三角测量）。我在下面包含了相关代码。我正在实现类似于Mark de Berg等人在“计算几何：算法与应用第三版”中描述的三角测量方法。下面是伪代码（如果需要，我将删除它）：
注意：我通过创建一个边框三角形来修改了伪代码，而不是使用“P中字典序最高的点”，因为我不太确定如何定义p-1和p-2，因为教科书上说要“符号化”定义它们，而不是定义确切的单位（公平地说，我可能误解了它的意思）。此外，合法化还不是我的代码的一部分（尚未），因为这对于Delaunay三角测量是必要的，但我想确保简单的三角测量按预期工作。
注意：我通过创建一个边框三角形来修改了伪代码，而不是使用“P中字典序最高的点”，因为我不太确定如何定义p-1和p-2，因为教科书上说要“符号化”定义它们，而不是定义确切的单位（公平地说，我可能误解了它的意思）。此外，合法化还不是我的代码的一部分（尚未），因为这对于Delaunay三角测量是必要的，但我想确保简单的三角测量按预期工作。
The issue is that I get some triangulations such as 
where the blue lines are from the original polygon.
Some of these lines don't get added because they are part of the triangles of points p0, p1, and p2 which I don't add in the findSmallest method. Yet, if I add those triangles as well, I get something like this:
(Note p0, p1, and p2 are beyond the scope of the picture). Some of the lines from the original polygon (this time in green) still aren't added to the triangulation. I'm not sure where I'm going wrong.
I hope the code is clear but I'm going to explain some of the methods/structs just in case:
TriPoint
is a child of the Point struct.
p0, p1, p2
are three points form a bounding triangle around the polygon. I got the idea from this post.
contains(Point p)
returns true if a point is within the triangle or on one of the edges.
findCommonTriangle(TriPoint *a, TriPoint *b, Triangle *t)
returns the triangle incident to t along edge ab. (I'm not using Edges to calculate the triangulation so I decided to get the incident triangle this way).
isOnTriangle(Point s)
is called on a Triangle abc, and returns 1 if the point is on edge ab, 2 if the pointis on edge bc, 3 if the point is on edge cd. If it is within the triangle, it returns 0.
The code for the triangulation itself is located below:
#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;
}
Note: predicates.cpp and predicates.h were created using code from here.
