c++ - 如何在 C++ 中正确三角化多边形

Question

我正在对一个对象进行三角剖分(最终，我想实现 Delaunay 三角剖分，但三角剖分甚至在边缘合法化之前也不起作用，所以我想首先关注一个简​​单的三角剖分)。我在下面包含了相关代码。我正在实现一种类似于 Mark de Berg 等人在“计算几何:算法和应用程序第三版”中描述的三角测量方法。给出的伪代码如下(如果需要我会删除它): 注意:我通过创建一个边界三角形而不是使用“P 的字典序最高点”来修改伪代码，因为我不太确定如何定义 p_-1 和 p_-2 正如教科书所说“象征性地”定义它们而不是定义确切的单位(公平地说，我很可能只是误解了它想说的话)。此外，合法化不是我的代码的一部分(目前)，因为这对于 Delaunay 三角剖分是必要的，但我想确保简单的三角剖分按预期工作。

问题是我得到了一些三角测量，例如 其中蓝线来自原始多边形。

其中一些线没有被添加，因为它们是我没有在 findSmallest 方法中添加的点 p0、p1 和 p2 的三角形的一部分。然而，如果我也添加这些三角形，我会得到这样的结果: (注意p0、p1、p2不在图片范围内)。原始多边形中的一些线(这次为绿色)仍未添加到三角测量中。我不确定哪里出错了。

我希望代码清晰，但为了以防万一，我将解释一些方法/结构:

TriPoint

是 Point 结构的子结构。

p0, p1, p2

三个点围绕多边形形成边界三角形。我的想法来自 this post .

contains(Point p) 

如果点在三角形内或其中一条边上，则返回真。

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 是使用来自 here 的代码创建的.

Answer 1

您的代码非常次优，但现在这并不重要(您正在学习，对吧？。我将专注于三角剖分问题。

已编辑:您初始化 yMin和 yMax Triangulation的成员在 sort()稍后将它们用于“大封闭”三角形。如果您决定不使用“sort()”，您将使用初始化值。设置一些默认值。

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

主要问题(在你的第一篇文章中，在你编辑它之前)我看到的是找到一个点是否在三角形内部的方法，在它的边界上，或在它之外。
Triangle::isOnTriangle()太可怕了。你算了几个crossproduct s 并返回“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)你必须确保所有的三角形都具有相同的方向(否则前一个表达式会被错误使用)。您可以特别注意传递给三角形 ctor 的点的顺序。或者您可以添加一个设置它们的“orientTri”函数。目前你的边界三角形是顺时针顺序的。最常见的顺序(也用于 OpenGL)是逆时针方向；但是你可以选择你喜欢的那个，只是要注意它。

对于 b) 将 float 与 '0' 进行比较不是一个好主意。在某些场景中，您可以使用 if std::abs(value) < someEpsilon .但特别是三角测量还不够。 Classroom Examples ofRobustness Problems in Geometric Computations很好地解释了为什么您的计算必须“稳健”。 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 中是昂贵的。也许您稍后会找到其他方法，当您获得网格划分的经验时。

在没有层次结构的情况下，“步行到点”方法(谷歌搜索)似乎使用随机输入更快。但是使用缓存(最后构建的三角形)要有效得多。

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