How to correctly triangulate a polygon in C++

Question

I'm working on triangulating an object (ultimately, I want to implement a Delaunay triangulation but the triangulation doesn't work even before legalizing edges, so I would like to focus on a simple triangulation first). I'm including the relevant code below. I'm implementing a triangulation method similar to the one described in "Computation Geometry: Algorithms and Application Third Edition" by Mark de Berg, among others. The pseudocode given is below (I'll remove it if need be): Note: I modified the pseudo code by creating a bordering triangle instead of using the "lexicographically highest point of P" because I wasn't too sure how to define p_-1 and p_-2 as the textbook says to define them "symbolically" rather than defining exact units (It's certainly possible that I simply misunderstood what it was trying to say, to be fair). Also, the legalizing isn't part of my code (yet) as that is necessary for the Delaunay Triangulation, but I want to make sure a simple triangulation works as intended.

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.

Giovanna

Please, review your blocks, I think you have unpaired {. And then reformat again your code with the right indentation. – Alonsoalonzo

@Alonsoalonzo this should be fixed now. Thank you. – Giovanna

@MaxVollmer I included the entire code at the bottom of the post. I was able to run this and reproduce the problem so I hope it's sufficient. Thank you. – Giovanna

Thank you, but this doesn't look like a minimal example. For example, I don't think the keyboard and mouse handling code is needed to reproduce your issue. Please try to reduce your code to the absolute minimum required to reproduce the issue. Otherwise your bug just drowns in noise and anyone who wants to help has to work through verbose code that isn't relevant. Thank you! – Nonconformist

@MaxVollmer My bad. I removed the controls and now the main method just draws a square and triangulates it. I also removed some more irrelevant methods. Thanks! – Giovanna

What kind of algorithm are you implementing? Ear Clipping or something more complicated? – Dogberry

@Giovanna What is exactly the algorithm you are trying to implement? Why are you shuffling the points after sorting them? What is that hardcoded value (20) in the triangulate method? – Inspired

@SergioMonteleone I'm trying to implement the algorithm as defined by the textbook that I mentioned. I was thinking of including the pseudo-code from the textbook but I wasn't sure if I should or not due to copyright issues or whatever else. As for shuffling, it suggests creating a random permutation of the points before adding them. I tried both without and with shuffling and my results were similar from what I could tell, so I'll likely remove the shuffle unless I see a significant difference. The 20 was from some other post regarding how to create a triangle around an entire polygon. – Giovanna

About "symbolically p0,p1,p2", Shewchuk uses a "ghost vertex". See Fig 3.6 at Lecture Notes on Delaunay Mesh Generation. This ghost vertex requires special deal, as opposed to use a big bounding triangle, but saves a bit of memory and avoids knowing in advance the bounding box. – Alonsoalonzo

Answer 1

Your code is quite sub-optimal, but that doesn't matter now (you're learning, right?. I'll focus on triangulation issues.

EDITED: You initialize yMin and yMax members of Triangulation at sort() and use them later for the "big enclosing" triangle. If you decide not to use "sort()" you'll use initialized values. Put some defaults on it.

Sorting the points is not needed for building the triangulation. You use it just to find the BB, too much effort, and at the end shuffle them, again too much effort.

The main issue (in your first post, before you edited it) I see is the way of finding if a point is inside a triangle, on its boundary, or outside of it.
Triangle::isOnTriangle() is horrible. You calculate several crossproducts and return '0' (inside triangle) is none of them equals '0'.
You may argue that you know in advance that the point is not outside because you tested it before by Triangle::contains(), but this function is also horrible, not that much though.

The best (or at least easiest and most used) way to find the relative position of a point to a line is

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

res is a positive value if {px,py} is at right of {x1,y1} to {x2,y2} line. Negative if at left and zero if it's on the line.
Two important things here:

• a) Swapping the direction (i.e. the line is {x2,y2} to {x1,y1}) changes the sign of res.
• b) Telling if res is really zero is not easy due to numerical issues, as with any other float-precision representation.

For a) you must be sure that all of triangles have the same orientation (or the former expression will be used wrongly). You can take extra-care on the order of points you pass to the triangle ctor. Or you can add a "orientTri" function that sets them. Currently your bounding triangle is clockwise-order. The most common order (also used in OpenGL) is counter-clockwise; but you can choose the one you like, just be aware of it.

For b) the comparison of a float with '0' is not a good idea. In some scenaries you can use if std::abs(value) < someEpsilon. But specially with triangulations that's not enough. Classroom Examples of Robustness Problems in Geometric Computations explains very well why your calculations must be "robust". Shewchuk Robust Predicates is a very good solution.

Once you have those two subjects solved the problem of "point in triangle" can be handled as this:

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
}

For the rest of the process of triangulation some advices:

Because you want a Delaunay triangulation, every time you add a new point you must check the "inCircle" condition (no other triangle whose circumcircle contains this new point). It can be done as shown in the book or in the links I posted. Again, you need robust predicates.

Shuffling the order of insertion of the points may improve performance for locating the triangle where the new point lies. This is may be true depending on the method used for the locating part. You use a hierarchy of triangles, so if the data is sorted or not does not affect.
By the way, maintaining the hierarchy structure for each triangle added/removed/changed is expensive in CPU and RAM. Perhaps you may find another ways later, when you get experience with meshing.

With no hierarchy, the "walk to point" method (google for it) seems faster with a randomized input. But using a cache (the last triangle built) is far more efficent.

Good luck with meshing. It's hard to begin with and to debug, the devil lives in the details.

Alonsoalonzo

I'd like to ask why yMin or yMax may cause a failure when I initalize them in the sort method? I don't actually use them until the triangulate method, which I call after sorting. And I wasn't sure if the sorting was necessary, so I'll likely remove it unless I see a significant improvement (which I have yet to see). I'm going to be adjusting my method then for determining if a point is within the triangle, so thank you! – Giovanna

You're right about yMin/yMax. I've edited my answer. – Alonsoalonzo

Answer 2

Apart from the issues already pointed out in Ripi2's answer, I'd like to suggest the following:

1. random_shuffle:

I see you do not initialize the random generator (by calling srand in your main function). That means you will always use the same pseudorandom sequence, hence executing the program several times will lead to the exact same result. I tested your code, and shuffling has indeed an impact on the code. After initializing the random generator correctly, you can see that the triangulation changes every time, producing different triangles.

2. comparisons:

In your code you compare points and triangles by comparing pointers. That means, as an example, that two points will be equal if and only if they are exactly the same struct in memory. Two point structs with the same coordinates will be considered different points. I am unsure if this is what you want to obtain or not so I would suggest to think about that.

3. triangle around polygon:

Aside from the hardcoded value (20), I don't understand why this code should produce a valid triangulation. You create several triangles over the polygon, but they all share one of those 3 fixed points outside the triangle.

Here is a picture taken reducing the hardcoded parameter to 2, to fit all the triangles in the viewport:

Same code, different point order (after initializing srand with time(0)):

I don't have access to the pseudocode of the algorithm, but I would suggest to edit your answer to briefly describe it, just for the sake of clarity.

Good luck with your implementation :)

Inspired

I added the pseudo code (I can remove it if necessary). I was more worried about making sure that the polygon was triangulated correctly and that all edges were included. I was planning on dealing with the outside points later. For example, if a given leaf node has 2 points of which are p0, p1, or p2, I can be certain that I don't need to include it. If only 1 point of the triangle is of those 3 points, then I would have to take more care into making sure relevant edges were included. If I could get that to work, I'd be satisfied but it relies on a proper triangulation to start. – Giovanna

About duplicating points, bad idea, because you will get 0-length edges and a mess in triangles connections. If you find a duplicate, just skip it. – Alonsoalonzo

@Giovanna Triangles that have p0,p1,p2 as one or two points are handled exactly the same as the rest. There are no shorteners. Just skip them at displaying stage. – Alonsoalonzo

@Alonsoalonzo Not displaying those triangles relies on the triangulation to be correct. I find that even after using orient2d() from Shewchuk, my result of the triangulation is essentially the same. I'm updating the minimal and complete example from my original post with the new code. I'm thinking part of my error may have been when I oriented triangle in orientTri(). Also, I call orient2d() on (b, a, p) to determine p's location relative to line segment AB (similarly for the other 2 segments), and I'm not sure if that's correct either but res1, res2, and res3 would all be positive otherwise. – Giovanna

@Toj19. I told you meshing is hard to begin with. Now you have corrected some parts of your code my best advice is that you use a debugger step by step to see what's happening. For example, a wrong assigment, wrong sign, messed hierachy, some vertex missed, etc. – Alonsoalonzo

@Alonsoalonzo Ok. Thanks for the help! – Giovanna