Triangulation of convex shapes (2D)

To draw freehand shapes in PyQt I wanted to sample the user’s mouse drag and calculate a triangulated polygon for drawing per-triangle.

Basic triangulation proved cumbersome with concave shapes, so here’s a more advanced method (complete class at the bottom):

I create a polygon that contains a list of points, each point is connected to the previous (and 0 is connected to the last), so the points are sorted as a sequence that form a closed polygon.

Then for every point I attempt to create a triangle between that point and the next two points. There are two problems when creating that triangle for a concave polygon. First the new triangle may be outside the polygon because the angle between the edges is more than 180 degrees, second the edge may intersect other edge of the polygon:

To avoid the first issue it is easiest to take the center of the new edge and see if that point lies in the polygon, for that the polygon requires and intersectPoint method; I made mine using the winding number algorithm described in A Winding Number and Point-in-Polygon Algorithm by D. G. Alciatore, R. Miranda

UPDATE: there is a problem with this method obviously, a better way to check whether the new edge lies outside the polygon is by checking whether the new edge creates a triangle with a different winding direction.

First we need the normal of the input polygon by taking the cross product of each connected edge-pair. In 2D the cross product only can be computed for the Z component. X and Y are always zero. So I am assuming you have a 2D cross product which returns a single float.
sign( (P0-P1) X (P2-P1) )
Doing this for each edge pair and adding the results together and taking the sign of that again gives us either 1 or -1. This determines whether the polygon winds clockwise or counterclockwise (which is not important).

Then for the new edge we can take the cross product with one of the adjacent points of the starting point. So consider the new edge consisting of P0 and PE, we have the adjacent point P1.
(P1-P0) X (PE-P0) should give the same sign as the polygon, if not, the edge is outside the polygon (if you’re getting invalid results, consider taking the other adjacent point).

Now you can ignore the next code block and continue reading.

    def intersectPoint(self, in_pt):
        '''
        winding number algorithm by D. G. Alciatore, R. Miranda
        http://www.engr.colostate.edu/~dga/dga/papers/point_in_polygon.pdf

        @param in_pt: Vec[2]: point to check
        '''
        w = 0
        ptlen = len(self.points)
        for i in range(ptlen):
            pt0, pt1 = self.points[i]-in_pt, self.points[(i+1)%ptlen]-in_pt

            if pt0[1] == pt1[1]: #parallel
                continue

            #up or down
            mod = ( pt0[1] < pt1[1] )*2-1

            #start/end on the axis
            if 0 in (pt0[1],pt1[1]):
                mod *= 0.5
            #line crossing X axis
            if ( pt0[1] <= 0 and pt1[1] >= 0 ) or ( pt0[1] >= 0 and pt1[1] <= 0 ):
                #line on the right of the pivot
                if pt0[0] >= 0 and pt1[0] >= 0:
                    w += mod
                    continue
                #get intersection point X
                xpery = (pt1[0]-pt0[0])/(pt1[1]-pt0[1])
                ix = xpery*-pt0[1]+pt0[0]
                #intersection on the right of the pivot
                if ix >= 0:
                    w += mod
        return w

To solve the self-intersection problem I will check the new line against every other line in the polygon, skipping lines that can not be the problem because they share one of the points with the new line.

This exclusion is given as parameter because not all line intersection checks are with a line of the same polygon. Excluding these lines however is necessary to avoid false positives.

def intersectLine(self, in_start, in_end, in_exclude = []):
        '''
        @param excluderange: list of int: excludes points with
        indices matching the given list, defaults to empty list
        '''
        ptlen = len(self.points)
        for i in range(ptlen):
            if i in in_exclude:
                continue
            tmp = Vmath.utils.lineLineIntersect2D(self.points[i], \
                                            self.points[(i+1)%ptlen], \
                                            in_start, in_end)
            if tmp[0]:
                return tmp
        return [False]

This method references lineLineIntersect2D from Vmath.utils. The line to line intersection formula comes from euclideanspace.com by Martin Baker. This is only a quick python implementation and it can be optimized, there are also some added checks to avoid division by 0 and related errors.

def lineLineIntersect2D(pt0, pt1, pt2, pt3):
    '''
    formula source from
    http://www.euclideanspace.com/maths/geometry/elements/intersection/twod/index.htm
    '''
    a,c = 0,0
    e = (pt1[0]-pt0[0])
    b = (pt1[0]*pt0[1] - pt0[0]*pt1[1])
    if e != 0:
        a = (pt1[1]-pt0[1]) / e #step size of line 1
        b /= e
    f = (pt3[0]-pt2[0])
    d = (pt3[0]*pt2[1] - pt2[0]*pt3[1])
    if f != 0:
        c = (pt3[1]-pt2[1]) / f #step size of line 2
        d /=  f

    g = a-c

    if g == 0: #lines are parallel
        if e == 0:
            out = Vec( pt0[0], (a*d-b*c) )
        else:
            out = Vec( pt2[0], (a*d-b*c) )
    else:
        out = Vec( (d-b) / g, (a*d-b*c) / g )

    if out[0] >= min(pt0[0],pt1[0]) and out[0] <= max(pt0[0],pt1[0]) and \
        out[0] >= min(pt2[0],pt3[0]) and out[0] <= max(pt2[0],pt3[0]) and \
        out[1] >= min(pt0[1],pt1[1]) and out[1] <= max(pt0[1],pt1[1]) and \
        out[1] >= min(pt2[1],pt3[1]) and out[1] <= max(pt2[1],pt3[1]):
            return [True, out]
    return [False, out]

Below is the final class, the triangulate method is the most interesting of course, it will duplicate all internal points and for every three points attempt to create a triangle. When succeeding the intermediate point is removed from the target points as it is now cutoff from the rest of the polygon by a new edge.

The ptmap variable contains the original point indices for self.points querying (i in range(ptlen) will contain correct data for the points duplicate but the intersection functions will use self.points and require the pmap). The triangle indices that are generated relate to self.points using ptmap because the internal copy (points) is mutated during the process and deleted at the end.

'''
Created on Nov 20, 2012

@author: TC
'''
import Vmath.utils
from Vmath.vec import Vec
reload(Vmath.utils)

class Polygon2D(object):
    '''
    Polygon contains a set of points
    and their triangulation data
    as well as intersection support
    for points, lines and other polygons
    '''

    def __init__(self, in_points):
        self.points = in_points
        self.triangulate()

    def triangulate(self):
        '''
        self.triangles contains ids of the
        points array for drawing purposes; it
        triangulates concave polygons properly
        '''
        self.triangles = []
        points = self.points[:]
        ptlen = len(points)
        ptmap = range(ptlen)
        prevptlen = ptlen
        i = 0
        while ptlen > 2:
            if ptlen == 3:
                self.triangles.extend(ptmap[0:3])
                break
            pt0, pt1 = points[i], points[(i+2)%ptlen]
            #exclude adjacent edges
            exclude = set([ptmap[i],
            (ptmap[i]-1)%len(self.points),
            (ptmap[i]+1)%len(self.points),
            (ptmap[(i+2)%ptlen])%len(self.points),
            (ptmap[(i+2)%ptlen]-1)%len(self.points),
            (ptmap[(i+2)%ptlen]+1)%len(self.points)])
            #if new line does not cross any outer edges
            #and midpoint of line is inside this poly
            if not self.intersectLine(pt0, pt1, exclude)[0] \
                and self.intersectPoint( (pt1-pt0)*0.5+pt0 ) != 0:
                self.triangles.extend([ptmap[i], ptmap[(i+1)%ptlen], ptmap[(i+2)%ptlen]])
                points.pop((i+1)%ptlen)
                ptmap.pop((i+1)%ptlen)
                ptlen-=1
            i+=1
            #wrap around
            if i > ptlen-1:
                #if no change since previous wrap, bail out of infinite loop
                if ptlen == prevptlen:
                    raise ValueError("Polygon could not be triangulated, it contains self intersection")
                    break
                i = 0
                prevptlen = ptlen

    def addPoints(self, in_points):
        self.points.extend(in_points)
        self.triangulate()

    def setPoints(self, in_points):
        self.points = in_points[:]
        self.triangulate()

    def intersectLine(self, in_start, in_end, in_exclude = []):
        '''
        @param excluderange: list of int: excludes points with
        indices matching the given list, defaults to empty list
        '''
        ptlen = len(self.points)
        for i in range(ptlen):
            if i in in_exclude:
                continue
            tmp = Vmath.utils.lineLineIntersect2D(self.points[i], \
                                            self.points[(i+1)%ptlen], \
                                            in_start, in_end)
            if tmp[0]:
                return tmp
        return [False]

    def intersectPoint(self, in_pt):
        '''
        winding number algorithm by D. G. Alciatore, R. Miranda
        http://www.engr.colostate.edu/~dga/dga/papers/point_in_polygon.pdf

        @param in_pt: Vec[2]: point to check
        '''
        w = 0
        ptlen = len(self.points)
        for i in range(ptlen):
            pt0, pt1 = self.points[i]-in_pt, self.points[(i+1)%ptlen]-in_pt

            if pt0[1] == pt1[1]: #parallel
                continue

            #up or down
            mod = ( pt0[1] < pt1[1] )*2-1

            #start/end on the axis
            if 0 in (pt0[1],pt1[1]):
                mod *= 0.5
            #line crossing X axis
            if ( pt0[1] <= 0 and pt1[1] >= 0 ) or ( pt0[1] >= 0 and pt1[1] <= 0 ):
                #line on the right of the pivot
                if pt0[0] >= 0 and pt1[0] >= 0:
                    w += mod
                    continue
                #get intersection point X
                xpery = (pt1[0]-pt0[0])/(pt1[1]-pt0[1])
                ix = xpery*-pt0[1]+pt0[0]
                #intersection on the right of the pivot
                if ix >= 0:
                    w += mod
        return w

And for any Maya users out there, here’s a quick example that draws curve for each triangle which I used to debug it:

from Vmath.vec import Vec
import Vmath.polygon
reload(Vmath.polygon)
from Vmath.polygon import Polygon2D

poly = Polygon2D([
Vec(1,0),
Vec(2.8,0),
Vec(3,2),
Vec(2,4),
Vec(0,3),
Vec(0,2),
Vec(1,2.2),
Vec(2.45,1),
Vec(0,1)])
from maya import cmds
for i in range(0, len(poly.triangles), 3):
    p0 = (poly.points[ poly.triangles[i] ][0], 0, poly.points[ poly.triangles[i] ][1])
    cmds.curve(d=1, p=[p0,
    (poly.points[ poly.triangles[i+1] ][0], 0,
    poly.points[ poly.triangles[i+1] ][1]),
    (poly.points[ poly.triangles[i+2] ][0], 0,
    poly.points[ poly.triangles[i+2] ][1]), p0], k=[0,1,2,3])

Leave a Reply

Your email address will not be published. Required fields are marked *