如何在C++中正确地对多边形进行三角剖分

Question

我工作的三角物体(最终,我想实现一个Delaunay三角,但三角甚至不合法化的边缘前工作,所以我想专注于一个简单的三角测量第一).我在下面列出了相关代码.我正在实施类似于Mark de Berg的"Computation Geometry:Algorithms and Application Third Edition"中描述的三角测量方法等.给出的伪代码如下(如果需要,我将删除它): 注意:我通过创建一个边界三角形来修改伪代码,而不是使用"按字典顺序排列的最高点P",因为我不太确定如何定义p _-1和p _-2,因为教科书说它"象征性地"定义它们而不是定义确切的单位(我当然可能误解了它试图说的是公平的).此外,合法化不是我的代码的一部分(但是)因为这对于Delaunay三角测量是必要的,但我想确保简单的三角测量按预期工作.

问题是我得到了一些诸如此类的三角测量 其中蓝线来自原始多边形.

其中一些行没有被添加,因为它们是点p0,p1和p2的三角形的一部分,我没有在findSmallest方法中添加它们.然而,如果我也添加这些三角形,我得到这样的东西:(注意p0,p1和p2超出了图片的范围).来自原始多边形的一些线(此时为绿色)仍未添加到三角测量中.我不确定我哪里出错了.

我希望代码是清楚的,但我将解释一些方法/结构以防万一:

TriPoint

是Point结构的子代.

p0, p1, p2

是三个点形成多边形周围的边界三角形.我从这篇文章中得到了这个想法.

contains(Point p) 

如果某个点位于三角形内或其中一条边上,则返回true.

findCommonTriangle(TriPoint *a, TriPoint *b, Triangle *t) 

沿着边缘ab返回t的三角形事件.(我没有使用Edges来计算三角测量因此我决定以这种方式得到事件三角形).

isOnTriangle(Point s)

在三角形abc上调用,如果该点位于边缘ab,则返回1;如果该点位于边缘bc,则返回2;如果该点位于边缘cd,则返回3.如果它在三角形内,则返回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是使用此处的代码创建的.

Answer 1

您的代码相当次优，但这现在并不重要（您正在学习，对吗？我将重点关注三角测量问题。

编辑：您初始化at的yMin和成员，并稍后将它们用于“大封闭”三角形。如果您决定不使用“sort()”，您将使用初始化值。对其进行一些默认设置。yMaxTriangulationsort()

构建三角测量不需要对点进行排序。你用它只是为了找到 BB，太费力了，最后洗牌，又太费力了。

我看到的主要问题（在您的第一篇文章中，在编辑之前）是查找点是否位于三角形内部、其边界上或外部的方法。
Triangle::isOnTriangle()太可怕了。您计算了几个crossproducts 并返回“0”（三角形内），因为它们都不等于“0”。
你可能会说你事先知道该点不在外面，因为你之前通过 测试过它Triangle::contains()，但这个函数也很糟糕，但没有那么糟糕。

查找点与线的相对位置的最佳（或至少最简单和最常用）方法是

res = (y2 - y1)*(px - x1) - (x2 - x1)*(py - y1)

res{px,py}如果位于线的右侧，则为正值{x1,y1} to {x2,y2}。如果在左边则为负，如果在线上则为零。
这里有两件重要的事情：

• a) 交换方向（即线为{x2,y2} to {x1,y1}）会改变 的符号res。
• b)res与任何其他浮点精度表示一样，由于数值问题，判断是否真的为零并不容易。

对于a），您必须确保所有三角形具有相同的方向（否则前一个表达式将被错误地使用）。您可以特别注意传递给三角形向量的点的顺序。或者您可以添加一个“orientTri”函数来设置它们。目前，您的边界三角形是顺时针顺序的。最常见的顺序（也在 OpenGL 中使用）是逆时针；但你可以选择你喜欢的，只要注意这一点即可。

对于 b) 将浮点数与“0”进行比较并不是一个好主意。在某些场景下你可以使用if std::abs(value) < someEpsilon. 但特别是对于三角测量来说，这还不够。几何计算中稳健性问题的课堂示例很好地解释了为什么您的计算必须“稳健”。Shewchuk Robust Predicates是一个非常好的解决方案。

一旦解决了这两个问题，“三角形中的点”问题就可以这样处理：

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”条件（外接圆中没有其他三角形包含此新点）。可以按照书中或我发布的链接中所示的方式完成。同样，您需要强大的谓词。

打乱点的插入顺序可以提高定位新点所在三角形的性能。根据用于定位部件的方法，这可能是正确的。您使用三角形的层次结构，因此数据是否排序不会产生影响。
顺便说一句，维护每个添加/删除/更改的三角形的层次结构对于 CPU 和 RAM 来说是昂贵的。当您获得网格划分经验后，也许您可​​能会找到其他方法。

由于没有层次结构，“步行到点”方法（谷歌搜索）似乎通过随机输入更快。但使用缓存（最后构建的三角形）效率要高得多。

祝您网格化好运。开始和调试都很困难，细节决定成败。