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

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])