Efficient algorithm for generating a triangle mesh from a concave polygon

Question

I am working on a project that involves generating meshes of triangles from potentially concave polygons. I am far more interested in a performant solution than an optimal solution. Most of what I've come across involves using OpenGL (not an option) or focuses on optimally or approximately optimally decomposing into convex polygons, with O(n³) or O(n⁴). Ear clipping seems like the most straightforward solution and is supposed to be O(n²):

for /* ever */ {
    for each(vertex) {
        candidate := Triangle{previousVertex, currentVertex, nextVertex}
        if !candidate.IsInterior() || candidate.ContainsAnyOtherVertex() {
            continue
        }

        result.Push(candidate)
        vertex = vertex.Except(currentVertex)
    }

    if len(vertex) == 3 {
        result.Push(Triangle(vertex))
        break
    }
}

I don't see how this can be O(n²). As I see it, this will have to involve three nested loops, each proportional to N: the two explicit loops above, and candidate.ContainsAnyOtherVertex(). Am I missing something? Is there some way to detect whether the candidate triangle contains any of the other remaining vertices that does not involve iterating over all the remaining vertices?

Alternatively, is there a different algorithm I should use? I'd be fine with one that decomposes into convex polygons, but A) I don't care if the solution is optimal, and B) I am not interested in reading through CGAL's source code or academic literature.

Answer 1

Ear clipping is O(n²) if you do it this way:

• First, for each vertex with an angle < 180 degrees, count the number of vertexes inside its angle. (O(n²))
• You can clip + remove any such vertex with 0 count. You will do this O(n) times. When you remove a vetex:
  • First remove it from the counts of any triangles that it's in (O(n)); and
  • Count the other vertexes in the at-most-2 new triangles you form by clipping (also O(n)).

Answer 2

For anyone else having trouble with this, here is a condensed version of what I'm using:

package mesh

import "sort"

type Vec []float32

type Face []Vec

type FaceTri struct {
    face  Face
    index [3]int
}

func triangulate(face Face, offset int) [][3]int {
    normal := face.Normal()

    index := make([]int, len(face))
    convex := map[int]bool{}
    reflex := map[int]bool{}
    ear := map[int]bool{}

    // Mark convex and reflex
    for i := range index {
        index[i] = i
        tri := face.tri(i, i+1, i+2)

        // Skip colinear vertices
        if tri.Area().Len() == 0 {
            continue
        }

        if tri.Area().Dot(normal) > 0 {
            convex[tri.index[1]] = true
        } else {
            reflex[tri.index[1]] = true
        }
    }

    // Mark ears
    for i := range convex {
        if isEar(face.tri(i-1, i, i+1), reflex) {
            ear[i] = true
        }
    }

    var tris [][3]int
    for len(reflex) > 0 || len(index) > 3 {
        ni := len(index)

        // Next ear
        var j int
        for j = range ear {
            break
        }

        // Find J in the list of unclipped vertices and get 2 neighbors to each side
        x := sort.SearchInts(index, j)
        h, i, k, l := (ni+x-2)%ni, (ni+x-1)%ni, (x+1)%ni, (x+2)%ni
        h, i, k, l = index[h], index[i], index[k], index[l]

        // Clip (I,J,K)
        index = append(index[:x], index[x+1:]...)
        tris = append(tris, [3]int{offset + i, offset + j, offset + k})
        delete(ear, j)
        delete(convex, j)

        // Update prior vertex
        update(face.tri(h, i, k), normal, reflex, convex, ear)

        // Update later vertex
        if h != l {
            update(face.tri(i, k, l), normal, reflex, convex, ear)
        }
    }

    tris = append(tris, [3]int{
        offset + index[0],
        offset + index[1],
        offset + index[2],
    })
    return tris
}

func update(tri *FaceTri, faceNormal Vec, reflex, convex, ear map[int]bool) {
    idx := tri.index[1]
    wasConvex, wasEar := convex[idx], ear[idx]

    isConvex := wasConvex || tri.Area().Dot(faceNormal) > 0
    if !wasConvex && isConvex {
        convex[idx] = true
        delete(reflex, idx)
    }

    if !wasEar && isConvex && isEar(tri, reflex) {
        ear[idx] = true
    }
}

func isEar(tri *FaceTri, reflex map[int]bool) bool {
    // It is sufficient to only check reflex vertices - a convex vertex
    // cannot be contained without a reflex vertex also being contained
    for j := range reflex {
        if tri.HasIndex(j) {
            continue
        }

        if tri.ContainsPoint(tri.face[j]) {
            return false
        }
    }

    return true
}

// Euclidean length of V
func (v Vec) Len() float32

// Dot product of A and B
func (a Vec) Dot(b Vec) float32

// Calculates the normal vector to the face - for concave polygons, this
// is the summation of the normals of each vertex (normalized)
func (Face) Normal() Vec

// Tri returns a FaceTri for I, J, K, modulo len(f)
func (Face) tri(i, j, k int) *FaceTri

// Area returns the cross product of the vector from v1 to v0 and the vector
// from v1 to v2
func (*FaceTri) Area() Vec

// Returns true if v is in f.index
func (f *FaceTri) HasIndex(v int) bool

// Returns true if P is contained within the described triangle
func (*FaceTri) ContainsPoint(p Vec) bool

This is based on the description from https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf. I think it is essentially the same as what Matt described.

For the sake of clarity and brevity, I omitted the contents of some methods, and cases where I reused results in a way that reduces readability. IMO the only interesting part I left out is ContainsPoint. In 2D, figuring out if triangle ABC contains point P is straight forward. In 3D, less so, since P is not necessarily co-planar with ABC:

V = vector from B to A
U = vector from B to C
Q = vector from B to P

M = MatrixFromColumns(V, U, V x U)
X = (inverse of M) * Q

Q === (V * X[0]) + (U * X[1]) + ((V x U) * X[2])

If X[2] is zero, then Q has a component out of the plane of the triangle

P is in ABC if and only if X[0] and X[1] are positive and X[0] + X[1] <= 1