seidel

python/framework/framework.pyx

@@ -5,13 +5,10 @@ from math import pi as PI
 
 from gl cimport *
 
-#from triangulator import Point
-
-include "triangulator.pyx"
-
 cdef extern from 'math.h':
     double cos(double)
     double sin(double)
+    double sqrt(double)
 
 SEGMENTS = 25
 INCREMENT = 2.0 * PI / SEGMENTS
@@ -61,11 +58,6 @@ from glfw cimport *
 
 import sys
 
-cdef extern from 'math.h':
-    double cos(double)
-    double sin(double)
-    double sqrt(double)
-
 # Keyboard callback wrapper
 kbd_callback_method = None
 
@@ -78,9 +70,6 @@ cdef class Game:
     title = "Poly2Tri"
     
     def __init__(self, window_width, window_height):
-        
-        points = [Point(100,100), Point(-100,100), Point(-100,-100), Point(100,-100)]
-        seidel = Triangulator(points)
      
         glfwInit()
         

python/framework/triangulator.pyx

@@ -43,145 +43,9 @@ cdef extern from 'math.h':
     double floor(double)
     double sqrt(double)
 
-class Triangulator:
+class Point:
-  
-    def __init__(self, points):
-        self.polygons = []
-        self.edge_list = self.init_edges(points)
-        self.trapezoids = []
-        self.trapezoidal_map = TrapezoidalMap()
-        bounding_box = self.trapezoidal_map.bounding_box(self.edge_list)
-        self.query_graph = QueryGraph(Sink(bounding_box))
-        self.xmono_poly = []
-        self.process()
-  
-    def trapezoidMap(self): 
-        return self.trapezoidal_map.map
     
-    # Build the trapezoidal map and query graph
+    def __cinit__(self, float x, float y):
-    def process(self):
-        for e in self.edge_list:
-            traps = self.query_graph.follow_edge(e)
-            for t in traps:
-                try:
-                    self.trapezoidal_map.map.remove(t)
-                except:
-                    pass
-            for t in traps:
-                tlist = []
-                cp = t.contains(e.p)
-                cq = t.contains(e.q)
-                if cp and cq:
-                    tlist = self.trapezoidal_map.case1(t, e)
-                    self.query_graph.case1(t.sink, e, tlist)
-                elif cp and not cq:
-                    tlist = self.trapezoidal_map.case2(t, e) 
-                    self.query_graph.case2(t.sink, e, tlist)
-                elif not cp and not cq:
-                    tlist = self.trapezoidal_map.case3(t, e)
-                    self.query_graph.case3(t.sink, e, tlist)
-                else:
-                    tlist = self.trapezoidal_map.case4(t, e)
-                    self.query_graph.case4(t.sink, e, tlist)
-                
-                # Add new trapezoids to map
-                for t in tlist:
-                  self.trapezoidal_map.map.append(t)
-          
-            self.trapezoidal_map.clear()
-
-        # Mark outside trapezoids
-        for t in self.trapezoidal_map.map:
-            self.mark_outside(t)
-        
-        # Collect interior trapezoids
-        for t in self.trapezoidal_map.map:
-            if t.inside():
-                self.trapezoids.append(t)
-                t.add_points()
-
-        self.create_mountains()
-
-    def mono_polies(self):
-        polies = []
-        for x in self.xmono_poly:
-            polies.append(x.monoPoly)
-        return polies
-  
-    def create_mountains(self):
-        for s in self.edge_list:          
-            if len(s.mpoints) > 0:
-                mountain = MonotoneMountain()
-                k = merge_sort(s.mpoints)
-                points = [s.p] + k + [s.q]
-                for p in points:
-                    mountain.append(p)
-                mountain.process()
-                for t in mountain.triangles:
-                    self.polygons.append(t)
-                self.xmono_poly.append(mountain)
-  
-    def mark_outside(self, t):
-        if t.top is self.bounding_box.top or t.bottom is self.bounding_box.bottom:
-            t.trimNeighbors()
-  
-    def init_edges(self, points):
-        edges = []
-        for i in range(len(points)-1):
-            edges.append(Edge(points[i], points[i+1]))
-        edges.append(Edge(points[0], points[-1]))
-        return self.order_edges(edges)
-  
-    def order_edges(self, edges):
-        segs = []
-        for s in edges:
-            p = self.shearTransform(s.p)
-            q = self.shearTransform(s.q)
-            if p.x > q.x: segs.append(Edge(q, p))
-            elif p.x < q.x: segs.append(Edge(p, q))
-        shuffle(segs)
-        return segs
-
-    def shearTransform(self, point):
-        return Point(point.x + 1e-4 * point.y, point.y)
- 
-cdef list merge_sort(list l):
-    cdef list lleft, lright
-    cdef int p1, p2, p
-    if len(l)>1 :
-        lleft = merge_sort(l[:len(l)/2])
-        lright = merge_sort(l[len(l)/2:])
-        p1, p2, p = 0, 0, 0
-        while p1<len(lleft) and p2<len(lright):
-            if lleft[p1].x < lright[p2].x:
-                l[p]=lleft[p1]
-                p+=1
-                p1+=1
-            else:
-                l[p]=lright[p2]
-                p+=1
-                p2+=1
-        if p1<len(lleft):l[p:]=lleft[p1:]
-        elif p2<len(lright):l[p:]=lright[p2:]
-        else : print "internal error"
-    return l
-  
-cdef class Point:
-  
-    cdef float x, y
-    
-    next = None
-    prev = None
-    edge = None
-    edges = []
-    
-    property x:
-        def __get__(self): return self.x
-    
-    property y:
-        def __get__(self): return self.y
-        
-    def __init__(self, float x, float y):
         self.x = x
         self.y = y
 
@@ -203,10 +67,10 @@ cdef class Point:
     def __div__(self, float a):
         return Point(self.x / a, self.y / a)
     
-    def cross(self, Point p):
+    def cross(self, p):
         return self.x * p.y - self.y * p.x
     
-    def dot(self, Point p):
+    def dot(self, p):
         return self.x * p.x + self.y * p.y
     
     def length(self):
@@ -215,38 +79,22 @@ cdef class Point:
     def normalize(self):
         return self / self.length() 
         
-    def less(self, Point p):
+    def less(self, p):
         return self.x < p.x
-  
-    '''
-    # Sort along y axis
-    def greater(self, p):
-        if y < p.y: 
-          return True
-        elif y > p.y:
-          return False
-        else:
-          if x < p.x:
-            return True
-          else:
-            return False
-    '''
     
     def not_equal(self, p):
         return not (p.x == self.x and p.y == self.y)
         
     def clone(self):
         return Point(self.x, self.y)
-  
-cdef class Edge:
 
-    cdef Point p, q
+cdef class Edge:
-    cdef bool above, below
-    cdef float slope, b
  
+    cdef object above, below
+    cdef float slope, b
     mpoints = []
     
-    def __init__(self, Point p, Point q):
+    def __cinit__(self, p, q):
         self.p = p
         self.q = q
         self.slope = (q.y - p.y)/(q.x - p.x)
@@ -264,12 +112,12 @@ cdef class Edge:
     property below:
         def __get__(self): return self.below
         
-    cdef bool is_above(self, Point point):
+    cdef bool is_above(self, point):
         return (floor(point.y) < floor(self.slope * point.x + self.b))
-    cdef bool is_below(self, Point point):
+    cdef bool is_below(self, point):
         return (floor(point.y) > floor(self.slope * point.x + self.b))
 
-    cdef float intersect(self, Point c, Point d):
+    cdef float intersect(self, c, d):
         cdef float a1, a2, a3, a4, t
         cdef Point a, b
         a = self.p
@@ -284,24 +132,18 @@ cdef class Edge:
                 return a + ((b - a) * t)
         return 0.0
         
-    cdef float signed_area(self, Point a, Point b, Point c):
+    cdef float signed_area(self, a, b, c):
         return (a.x - c.x) * (b.y - c.y) - (a.y - c.y) * (b.x - c.x)
-
-cdef Point line_intersect(Edge e, float x):
-    cdef float y =  e.slope * x + e.b
-    return Point(x, y)
         
 cdef class Trapezoid:
-
+  
-    cdef:
+    cdef Edge top, bottom
-        Point left_point, right_point
+    cdef Trapezoid upper_left, lower_left
-        Edge top, bottom
+    cdef Trapezoid upper_right, lower_right
-        Trapezoid upper_left, lower_left
+    cdef object sink
-        Trapezoid upper_right, lower_right
+    cdef bool inside
-        bool inside
-        object sink
     
-    def __init__(self, Point left_point, Point right_point, Edge top, Edge bottom):
+    def __cinit__(self, left_point, right_point, Edge top, Edge bottom):
         self.left_point = left_point
         self.right_point = right_point
         self.top = top
@@ -313,6 +155,9 @@ cdef class Trapezoid:
         self.inside = True
         self.sink = None
     
+    property inside:
+        def __get__(self): return self.inside
+        
     property top:
         def __get__(self): return self.top
         
@@ -321,9 +166,11 @@ cdef class Trapezoid:
         
     property left_point:
         def __get__(self): return self.left_point
-    
+        def __set__(self, lp): self.left_point = lp
+        
     property right_point:
         def __get__(self): return self.right_point
+        def __set__(self, rp): self.right_point = rp
         
     property sink:
         def __get__(self): return self.sink
@@ -375,7 +222,7 @@ cdef class Trapezoid:
             if self.upper_right != None: self.upper_right.trim_neighbors()
             if self.lower_right != None: self.lower_right.trim_neighbors()
   
-    def contains(self, Point point):
+    def contains(self, point):
         return (point.x > self.left_point.x and point.x < self.right_point.x and 
                 self.top.is_above(point) and self.bottom.is_below(point))
   
@@ -397,6 +244,140 @@ cdef class Trapezoid:
         if self.right_point != self.top.q: 
             self.top.mpoints.append(self.right_point.clone)
 
+cdef Point line_intersect(Edge e, float x):
+    cdef float y =  e.slope * x + e.b
+    return Point(x, y)
+    
+class Triangulator:
+  
+    def __init__(self, poly_line):
+        self.polygons = []
+        self.edge_list = self.init_edges(poly_line)
+        self.trapezoids = []
+        self.trapezoidal_map = TrapezoidalMap()
+        self.bounding_box = self.trapezoidal_map.bounding_box(self.edge_list)
+        self.query_graph = QueryGraph(isink(self.bounding_box))
+        self.xmono_poly = []
+        
+        self.process()
+  
+    def trapezoidMap(self): 
+        return self.trapezoidal_map.map
+    
+    # Build the trapezoidal map and query graph
+    def process(self):
+    
+        for e in self.edge_list:
+            traps = self.query_graph.follow_edge(e)
+            for t in traps:
+                try:
+                    self.trapezoidal_map.map.remove(t)
+                except:
+                    pass
+            for t in traps:
+                tlist = []
+                cp = t.contains(e.p)
+                cq = t.contains(e.q)
+                if cp and cq:
+                    tlist = self.trapezoidal_map.case1(t, e)
+                    self.query_graph.case1(t.sink, e, tlist)
+                elif cp and not cq:
+                    tlist = self.trapezoidal_map.case2(t, e) 
+                    self.query_graph.case2(t.sink, e, tlist)
+                elif not cp and not cq:
+                    tlist = self.trapezoidal_map.case3(t, e)
+                    self.query_graph.case3(t.sink, e, tlist)
+                else:
+                    tlist = self.trapezoidal_map.case4(t, e)
+                    self.query_graph.case4(t.sink, e, tlist)
+                # Add new trapezoids to map
+                for t in tlist:
+                  self.trapezoidal_map.map.append(t)
+            self.trapezoidal_map.clear()
+        
+        # Mark outside trapezoids
+        for t in self.trapezoidal_map.map:
+            self.mark_outside(t)
+        
+        # Collect interior trapezoids
+        for t in self.trapezoidal_map.map:
+            if t.inside:
+                self.trapezoids.append(t)
+                t.add_points()
+
+        self.create_mountains()
+
+    def mono_polies(self):
+        polies = []
+        for x in self.xmono_poly:
+            polies.append(x.monoPoly)
+        return polies
+  
+    def create_mountains(self):
+        for s in self.edge_list:          
+            if len(s.mpoints) > 0:
+                mountain = MonotoneMountain()
+                print s.mpoints
+                k = merge_sort(s.mpoints)
+                points = [s.p] + k + [s.q]
+                for p in points:
+                    mountain.append(p)
+                mountain.process()
+                for t in mountain.triangles:
+                    self.polygons.append(t)
+                self.xmono_poly.append(mountain)
+  
+    def mark_outside(self, t):
+        if t.top is self.bounding_box.top or t.bottom is self.bounding_box.bottom:
+            t.trim_neighbors()
+  
+    def init_edges(self, points):
+        edges = []
+        for i in range(len(points)-1):
+            p = Point(points[i][0], points[i][1])
+            q = Point(points[i+1][0], points[i+1][1])
+            edges.append(Edge(p, q))
+        p = Point(points[0][0], points[0][1])
+        q = Point(points[-1][0], points[-1][1])
+        edges.append(Edge(p, q))
+        return self.order_edges(edges)
+  
+    def order_edges(self, edges):
+        segs = []
+        for s in edges:
+            p = self.shearTransform(s.p)
+            q = self.shearTransform(s.q)
+            if p.x > q.x: 
+                segs.append(Edge(q, p))
+            elif p.x < q.x: 
+                segs.append(Edge(p, q))
+        shuffle(segs)
+        return segs
+
+    def shearTransform(self, point):
+        return Point(point.x + 1e-4 * point.y, point.y)
+ 
+cdef list merge_sort(l):
+    cdef list lleft, lright
+    cdef int p1, p2, p
+    if len(l)>1 :
+        lleft = merge_sort(l[:len(l)/2])
+        lright = merge_sort(l[len(l)/2:])
+        p1, p2, p = 0, 0, 0
+        while p1<len(lleft) and p2<len(lright):
+            if lleft[p1].x < lright[p2].x:
+                l[p]=lleft[p1]
+                p+=1
+                p1+=1
+            else:
+                l[p]=lright[p2]
+                p+=1
+                p2+=1
+        if p1<len(lleft):l[p:]=lleft[p1:]
+        elif p2<len(lright):l[p:]=lright[p2:]
+        else : print "internal error"
+    return l
+
 class TrapezoidalMap:
 
     map = []
@@ -408,7 +389,7 @@ class TrapezoidalMap:
         self.bcross = None
         self.tcross = None
 
-    def case1(self, t, e):
+    def case1(self, Trapezoid t, Edge e):
         trapezoids = []
         trapezoids.append(Trapezoid(t.left_point, e.p, t.top, t.bottom))
         trapezoids.append(Trapezoid(e.p, e.q, t.top, e))
@@ -420,7 +401,7 @@ class TrapezoidalMap:
         trapezoids[3].update_right(t.upper_right, t.lower_right)
         return trapezoids
 
-    def case2(self, t, e):
+    def case2(self, Trapezoid t, Edge e):
         rp = e.q if e.q.x == t.right_point.x else t.right_point
         trapezoids = []
         trapezoids.append(Trapezoid(t.left_point, e.p, t.top, t.bottom))
@@ -435,7 +416,7 @@ class TrapezoidalMap:
         e.below = trapezoids[2]
         return trapezoids
   
-    def case3(self, t, e):
+    def case3(self, Trapezoid t, Edge e):
         lp = e.p if e.p.x == t.left_point.x  else t.left_point
         rp = e.q if e.q.x == t.right_point.x else t.right_point
         trapezoids = []
@@ -459,7 +440,7 @@ class TrapezoidalMap:
         e.below = trapezoids[1]
         return trapezoids
 
-    def case4(self, t, e):
+    def case4(self, Trapezoid t, Edge e):
         lp = e.p if e.p.x == t.left_point.x else t.left_point
         trapezoids = []
         if self.tcross is t.top:
@@ -498,61 +479,81 @@ class TrapezoidalMap:
         right = top.q
         return Trapezoid(left, right, top, bottom)
 
-class Node:
+cdef class Node:
 
-    parent_list = []
+    cdef Node left, right
+    cdef object parent_list
     
-    def __init__(self, left, right):
+    def __init__(self, Node left, Node right):
+        self.parent_list = []
         self.left = left
         self.right = right
-        if left != None: left.parent_list.append(self)
+        if left != None: 
-        if right != None: right.parent_list.append(self)
+            left.parent_list.append(self)
-  
+        if right != None: 
-    def replace(self, node):
+            right.parent_list.append(self)
+    
+    property left:
+        def __get__(self): return self.left
+        def __set__(self, Node left): self.left = left
+    
+    property right:
+        def __get__(self): return self.right
+        def __set__(self, Node right): self.right = right
+    
+    property parent_list:
+        def __get__(self): return self.parent_list
+        
+    def replace(self, Node node):
         for parent in node.parent_list:
             if parent.left is node: 
                 parent.left = self
             else: 
                 parent.right = self 
-            self.parent_list.append(parent)
+        self.parent_list += node.parent_list
-
+    
-class Sink(Node):
+cdef class Sink(Node):
-
+    
-    def __new__(cls, trapezoid):
+    cdef Trapezoid trapezoid
-        if trapezoid.sink != None: 
+    
-            return trapezoid.sink
-        return Sink(trapezoid)
-            
     def __init__(self, trapezoid):
         self.trapezoid = trapezoid
-        Node.__init__(self, None, None)
+        super(Sink, self).__init__(None, None)
         trapezoid.sink = self
   
+    property trapezoid:
+        def __get__(self): return self.trapezoid
+        
     def locate(self, e): 
         return self
-        
-class XNode(Node):
 
-    def __init__(self, point, lchild, rchild):
+cdef Sink isink(Trapezoid trapezoid):
-        Node.__init__(self, lchild, rchild)
+    if trapezoid.sink != None: 
-        self.point = point
+        return trapezoid.sink
-        self.lchild = lchild
+    return Sink(trapezoid)
-        self.rchild = rchild
     
-    def locate(self, e): 
+cdef class XNode(Node):
+
+    cdef Point point
+    
+    def __init__(self, Point point, Node lchild, Node rchild):
+        super(XNode, self).__init__(lchild, rchild)
+        self.point = point
+    
+    def locate(self, Edge e): 
         if e.p.x >= self.point.x: 
             return self.right.locate(e)
         return self.left.locate(e)
 
-class YNode(Node):
+cdef class YNode(Node):
 
-    def __init__(self, edge, lchild, rchild):
+    cdef Edge edge
-        Node.__init__(self, lchild, rchild)
+    
+    def __init__(self, Edge edge, Node lchild, Node rchild):
+        super(YNode, self).__init__(lchild, rchild)
         self.edge = edge
-        self.lchild = lchild
-        self.rchild = rchild
         
-    def locate(self, e):
+    def locate(self, Edge e):
         if self.edge.is_above(e.p): 
             return self.right.locate(e)
         elif self.edge.is_below(e.p): 
@@ -560,53 +561,54 @@ class YNode(Node):
         else:
             if e.slope < self.edge.slope: 
                 return self.right.locate(e)
-        return self.left.locate(e)
+            else:
+                return self.left.locate(e)
             
-class QueryGraph:
+cdef class QueryGraph:
 
-    head = None
+    cdef Node head
     
-    def __init__(self, head):
+    def __init__(self, Node head):
         self.head = head
         
-    def locate(self, e):
+    def locate(self, Edge e):
         return self.head.locate(e).trapezoid
   
-    def follow_edge(self, e):
+    def follow_edge(self, Edge e):
         trapezoids = [self.locate(e)]
-        j = 0
+        cdef int j = 0
         while(e.q.x > trapezoids[j].right_point.x):
-            if e > trapezoids[j].right_point:
+            if e.is_above(trapezoids[j].right_point):
                 trapezoids.append(trapezoids[j].upper_right)
             else:
-                trapezoids .append(trapezoids[j].lower_right)
+                trapezoids.append(trapezoids[j].lower_right)
             j += 1
         return trapezoids
   
-    def replace(self, sink, node):
+    def replace(self, Sink sink, Node node):
         if not sink.parent_list:
           self.head = node
         else:
           node.replace(sink)
 
-    def case1(self, sink, e, tlist):
+    def case1(self, Sink sink, Edge e, tlist):
-        yNode = YNode(e, Sink(tlist[1]), Sink(tlist[2]))
+        cdef Node yNode = YNode(e, isink(tlist[1]), isink(tlist[2]))
-        qNode = XNode(e.q, yNode, Sink(tlist[3]))
+        cdef Node qNode = XNode(e.q, yNode, isink(tlist[3]))
-        pNode = XNode(e.p, Sink(tlist[0]), qNode)
+        cdef Node pNode = XNode(e.p, isink(tlist[0]), qNode)
         self.replace(sink, pNode)
   
-    def case2(self, sink, e, tlist):
+    def case2(self, Sink sink, Edge e, tlist):
-        yNode = YNode(e, Sink(tlist[1]), Sink(tlist[2]))
+        yNode = YNode(e, isink(tlist[1]), isink(tlist[2]))
-        pNode = XNode(e.p, Sink(tlist[0]), yNode)
+        pNode = XNode(e.p, isink(tlist[0]), yNode)
         self.replace(sink, pNode)
   
-    def case3(self, sink, e, tlist):
+    def case3(self, Sink sink, Edge e, tlist):
-        yNode = YNode(e, Sink(tlist[0]), Sink(tlist[1]))
+        yNode = YNode(e, isink(tlist[0]), isink(tlist[1]))
         self.replace(sink, yNode)
 
-    def case4(self, sink, e, tlist):
+    def case4(self, Sink sink, Edge e, tlist):
-        yNode = YNode(e, Sink(tlist[0]), Sink(tlist[1]))
+        yNode = YNode(e, isink(tlist[0]), isink(tlist[1]))
-        qNode = XNode(e.q, yNode, Sink(tlist[2]))
+        qNode = XNode(e.q, yNode, isink(tlist[2]))
         self.replace(sink, qNode)
 
 cdef float PI_SLOP = 3.1
@@ -624,7 +626,8 @@ cdef class MonotoneMountain:
 
     def __init__(self):
         self.size = 0
-        self.tail, self.head = None
+        self.tail = None
+        self.head = None
         self.positive = False
         self.convex_points = []
         self.mono_poly = []

python/poly2tri.py

@@ -1,20 +1,48 @@
 #!/usr/bin/env python2.6
-from framework import Game
+from framework import Game, draw_polygon, reset_zoom
- 
+
+from seidel import Triangulator
+
 class Poly2Tri(Game):
 
     #Screen size
     screen_size = 800.0, 600.0
 
     def __init__(self):
-        super(Poly2Tri, self).__init__(*self.screen_size)
+        super(Poly2Tri, self).__init__(*self.screen_size)       
+        
+        # Load point set
+        file_name = "../data/star.dat"
+        points = self.load_points(file_name)
+        
+        # Triangulate
+        t1 = self.time
+        seidel = Triangulator(points)
+        self.triangles = seidel.triangles()
+        dt = self.time - t1
+        print "time = %f  , num triangles = %d" % (dt, len(self.triangles))
+            
         self.main_loop()
         
     def update(self):
         pass
         
     def render(self):
-        pass
+        reset_zoom(1.0, (0,0), self.screen_size)
+        red = 255, 0, 0
+        for t in self.triangles:
+            draw_polygon(t, red)
+            
+    def load_points(self, file_name):
+        infile = open(file_name, "r")
+        points = []
+        while infile:
+            line = infile.readline()
+            s = line.split()
+            if len(s) == 0:
+                break
+            points.append((float(s[0]), float(s[1])))
+        return points
         
 if __name__ == '__main__':
     demo = Poly2Tri()

python/seidel.py

@@ -0,0 +1,622 @@
+#
+# Poly2Tri
+# Copyright (c) 2009, Mason Green
+# http://code.google.com/p/poly2tri/
+# 
+# All rights reserved.
+# 
+# Redistribution and use in source and binary forms, with or without modification,
+# are permitted provided that the following conditions are met:
+# 
+# Redistributions of source code must retain the above copyright notice,
+# self list of conditions and the following disclaimer.
+# Redistributions in binary form must reproduce the above copyright notice,
+# self list of conditions and the following disclaimer in the documentation
+# and/or other materials provided with the distribution.
+# Neither the name of Poly2Tri nor the names of its contributors may be
+# used to endorse or promote products derived from self software without specific
+# prior written permission.
+# 
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+# 
+from random import shuffle
+from math import atan2, floor
+
+###
+### Based on Raimund Seidel'e paper "A simple and fast incremental randomized
+### algorithm for computing trapezoidal decompositions and for triangulating polygons"
+### (Ported from poly2tri)
+
+class Point(object):
+    
+    def __init__(self, x, y):
+        self.x = x
+        self.y = y
+        self.next, self.prev = None, None
+
+    def __sub__(self, other):
+        if isinstance(other, Point):
+            return Point(self.x - other.x, self.y - other.y) 
+        else:
+            return Point(self.x - other, self.y - other)
+    
+    def __add__(self, other):
+        if isinstance(other, Point):
+            return Point(self.x + other.x, self.y + other.y) 
+        else:
+            return Point(self.x + other, self.y + other)
+    
+    def __mul__(self, f):
+        return Point(self.x * f, self.y * f)
+    
+    def __div__(self, a):
+        return Point(self.x / a, self.y / a)
+    
+    def cross(self, p):
+        return self.x * p.y - self.y * p.x
+    
+    def dot(self, p):
+        return self.x * p.x + self.y * p.y
+    
+    def length(self):
+        return sqrt(self.x * self.x + self.y * self.y)
+        
+    def normalize(self):
+        return self / self.length() 
+        
+    def less(self, p):
+        return self.x < p.x
+    
+    def not_equal(self, p):
+        return not (p.x == self.x and p.y == self.y)
+        
+    def clone(self):
+        return Point(self.x, self.y)
+
+class Edge(object):
+
+    mpoints = []
+    above, below = None, None
+    
+    def __init__(self, p, q):
+        self.p = p
+        self.q = q
+        self.slope = 0.0 if (q.x - p.x) == 0 else (q.y - p.y)/(q.x - p.x)
+        self.b = p.y - (p.x * self.slope)
+          
+    def is_above(self, point):
+        return (floor(point.y) < floor(self.slope * point.x + self.b))
+        
+    def is_below(self, point):
+        return (floor(point.y) > floor(self.slope * point.x + self.b))
+
+    def intersect(self, c, d):
+        a = self.p
+        b = self.q
+        a1 = self.signed_area(a, b, d)
+        a2 = self.signed_area(a, b, c)
+        if  a1 != 0 and a2 != 0 and (a1 * a2) < 0:
+            a3 = self.signed_area(c, d, a)
+            a4 = a3 + a2 - a1
+            if a3 * a4 < 0:
+                t = a3 / (a3 - a4)
+                return a + ((b - a) * t)
+        return 0.0
+        
+    def signed_area(self, a, b, c):
+        return (a.x - c.x) * (b.y - c.y) - (a.y - c.y) * (b.x - c.x)
+        
+class Trapezoid(object):
+    
+    def __init__(self, left_point, right_point, top, bottom):
+        self.left_point = left_point
+        self.right_point = right_point
+        self.top = top
+        self.bottom = bottom
+        self.upper_left = None
+        self.upper_right = None
+        self.lower_left = None
+        self.lower_right = None
+        self.inside = True
+        self.sink = None
+        
+    def update_left(self, ul, ll):
+        self.upper_left = ul
+        self.lower_left = ll
+        if ul != None: ul.upper_right = self
+        if ll != None: ll.lower_right = self  
+  
+    def update_right(self, ur, lr):
+        self.upper_right = ur
+        self.lower_right = lr
+        if ur != None: ur.upper_left = self
+        if lr != None: lr.lower_left = self   
+
+    def update_left_right(self, ul, ll, ur, lr):
+        self.upper_left = ul
+        self.lower_left = ll 
+        self.upper_right = ur 
+        self.lower_right = lr   
+        if ul != None: ul.upper_right = self
+        if ll != None: ll.lower_right = self
+        if ur != None: ur.upper_left = self
+        if lr != None: lr.lower_left = self  
+        
+    def trim_neighbors(self):
+        if self.inside:
+            self.inside = False
+            if self.upper_left != None: self.upper_left.trim_neighbors()
+            if self.lower_left != None: self.lower_left.trim_neighbors()
+            if self.upper_right != None: self.upper_right.trim_neighbors()
+            if self.lower_right != None: self.lower_right.trim_neighbors()
+  
+    def contains(self, point):
+        return (point.x > self.left_point.x and point.x < self.right_point.x and 
+                self.top.is_above(point) and self.bottom.is_below(point))
+  
+    def vertices(self):
+        verts = []
+        verts.append(line_intersect(self.top, self.left_point.x))
+        verts.append(line_intersect(self.bottom, self.left_point.x))
+        verts.append(line_intersect(self.bottom, self.right_point.x))
+        verts.append(line_intersect(self.top, self.right_point.x))
+        return verts
+  
+    def add_points(self):
+        if self.left_point != self.bottom.p: 
+            self.bottom.mpoints.append(self.left_point.clone())
+        if self.right_point != self.bottom.q: 
+            self.bottom.mpoints.append(self.right_point.clone())
+        if self.left_point != self.top.p: 
+            self.top.mpoints.append(self.left_point.clone())
+        if self.right_point != self.top.q: 
+            self.top.mpoints.append(self.right_point.clone())
+
+def line_intersect(edge, x):
+    y =  edge.slope * x + edge.b
+    return Point(x, y)
+    
+class Triangulator(object):
+  
+    def __init__(self, poly_line):
+    
+        self.polygons = []
+        self.edge_list = self.init_edges(poly_line)
+        self.trapezoids = []
+        self.trapezoidal_map = TrapezoidalMap()
+        self.bounding_box = self.trapezoidal_map.bounding_box(self.edge_list)
+        self.query_graph = QueryGraph(isink(self.bounding_box))
+        self.xmono_poly = []
+        self.process()        
+  
+    def triangles(self):
+        triangles = []
+        for p in self.polygons:
+            verts = []
+            for v in p:
+                verts.append((v.x, v.y))
+            triangles.append(verts)
+        return triangles
+            
+    def trapezoid_map(self): 
+        return self.trapezoidal_map.map
+    
+    # Build the trapezoidal map and query graph
+    def process(self):
+        for edge in self.edge_list:
+            traps = self.query_graph.follow_edge(edge)
+            for t in traps:
+                try:
+                    self.trapezoidal_map.map.remove(t)
+                except:
+                    pass
+            for t in traps:
+                cp = t.contains(edge.p)
+                cq = t.contains(edge.q)
+                if cp and cq:
+                    tlist = self.trapezoidal_map.case1(t, edge)
+                    self.query_graph.case1(t.sink, edge, tlist)
+                elif cp and not cq:
+                    tlist = self.trapezoidal_map.case2(t, edge) 
+                    self.query_graph.case2(t.sink, edge, tlist)
+                elif not cp and not cq:
+                    tlist = self.trapezoidal_map.case3(t, edge)
+                    self.query_graph.case3(t.sink, edge, tlist)
+                else:
+                    tlist = self.trapezoidal_map.case4(t, edge)
+                    self.query_graph.case4(t.sink, edge, tlist)
+                # Add new trapezoids to map
+                for t in tlist:
+                  self.trapezoidal_map.map.append(t)
+            self.trapezoidal_map.clear()
+        
+        # Mark outside trapezoids
+        for t in self.trapezoidal_map.map:
+            self.mark_outside(t)
+        
+        # Collect interior trapezoids
+        for t in self.trapezoidal_map.map:
+            if t.inside:
+                self.trapezoids.append(t)
+                t.add_points()
+
+        # Generate the triangles
+        self.create_mountains()
+
+    def mono_polies(self):
+        polies = []
+        for x in self.xmono_poly:
+            polies.append(x.monoPoly)
+        return polies
+  
+    def create_mountains(self):
+        for edge in self.edge_list:          
+            if len(edge.mpoints) > 0:
+                mountain = MonotoneMountain()
+                k = merge_sort(edge.mpoints)
+                points = [edge.p] + k + [edge.q]
+                for p in points:
+                    mountain.add(p)
+                mountain.process()
+                for t in mountain.triangles:
+                    self.polygons.append(t)
+                self.xmono_poly.append(mountain)
+  
+    def mark_outside(self, t):
+        if t.top is self.bounding_box.top or t.bottom is self.bounding_box.bottom:
+            t.trim_neighbors()
+  
+    def init_edges(self, points):
+        edges = []
+        for i in range(len(points)-1):
+            p = Point(points[i][0], points[i][1])
+            q = Point(points[i+1][0], points[i+1][1])
+            edges.append(Edge(p, q))
+        p = Point(points[0][0], points[0][1])
+        q = Point(points[-1][0], points[-1][1])
+        edges.append(Edge(p, q))
+        return self.order_edges(edges)
+  
+    def order_edges(self, edge_list):
+        edges = []
+        for e in edge_list:
+            p = self.shear_transform(e.p)
+            q = self.shear_transform(e.q)
+            if p.x > q.x: 
+                edges.append(Edge(q, p))
+            elif p.x < q.x: 
+                edges.append(Edge(p, q))
+        shuffle(edges)
+        return edges
+
+    def shear_transform(self, point):
+        return Point(point.x + 1e-4 * point.y, point.y)
+ 
+def merge_sort(l):
+    if len(l)>1 :
+        lleft = merge_sort(l[:len(l)/2])
+        lright = merge_sort(l[len(l)/2:])
+        p1, p2, p = 0, 0, 0
+        while p1<len(lleft) and p2<len(lright):
+            if lleft[p1].x < lright[p2].x:
+                l[p]=lleft[p1]
+                p+=1
+                p1+=1
+            else:
+                l[p]=lright[p2]
+                p+=1
+                p2+=1
+        if p1<len(lleft):l[p:]=lleft[p1:]
+        elif p2<len(lright):l[p:]=lright[p2:]
+        else : print "internal error"
+    return l
+
+class TrapezoidalMap(object):
+
+    map = []
+    margin = 50
+    bcross = None
+    tcross = None
+
+    def clear(self):
+        self.bcross = None
+        self.tcross = None
+        map = []
+
+    def case1(self, trapezoid, edge):
+        t = trapezoid; e = edge
+        trapezoids = []
+        trapezoids.append(Trapezoid(t.left_point, e.p, t.top, t.bottom))
+        trapezoids.append(Trapezoid(e.p, e.q, t.top, e))
+        trapezoids.append(Trapezoid(e.p, e.q, e, t.bottom))
+        trapezoids.append(Trapezoid(e.q, t.right_point, t.top, t.bottom))
+        trapezoids[0].update_left(t.upper_left, t.lower_left)
+        trapezoids[1].update_left_right(trapezoids[0], None, trapezoids[3], None)
+        trapezoids[2].update_left_right(None, trapezoids[0], None, trapezoids[3])
+        trapezoids[3].update_right(t.upper_right, t.lower_right)
+        return trapezoids
+
+    def case2(self, trapezoid, edge):
+        t = trapezoid; e = edge
+        rp = e.q if e.q.x == t.right_point.x else t.right_point
+        trapezoids = []
+        trapezoids.append(Trapezoid(t.left_point, e.p, t.top, t.bottom))
+        trapezoids.append(Trapezoid(e.p, rp, t.top, e))
+        trapezoids.append(Trapezoid(e.p, rp, e, t.bottom))
+        trapezoids[0].update_left(t.upper_left, t.lower_left)
+        trapezoids[1].update_left_right(trapezoids[0], None, t.upper_right, None)
+        trapezoids[2].update_left_right(None, trapezoids[0], None, t.lower_right)
+        self.bcross = t.bottom
+        self.tcross = t.top
+        e.above = trapezoids[1]
+        e.below = trapezoids[2]
+        return trapezoids
+  
+    def case3(self, trapezoid, edge):
+        t = trapezoid; e = edge
+        lp = e.p if e.p.x == t.left_point.x  else t.left_point
+        rp = e.q if e.q.x == t.right_point.x else t.right_point
+        trapezoids = []
+        if self.tcross is t.top:
+            trapezoids.append(t.upper_left)
+            trapezoids[0].update_right(t.upper_right, None)
+            trapezoids[0].right_point = rp
+        else:
+            trapezoids.append(Trapezoid(lp, rp, t.top, e))
+            trapezoids[0].update_left_right(t.upper_left, e.above, t.upper_right, None)
+        if self.bcross is t.bottom:
+            trapezoids.append(t.lower_left)
+            trapezoids[1].update_right(None, t.lower_right)
+            trapezoids[1].right_point = rp
+        else:
+            trapezoids.append(Trapezoid(lp, rp, e, t.bottom))
+            trapezoids[1].update_left_right(e.below, t.lower_left, None, t.lower_right)
+        self.bcross = t.bottom
+        self.tcross = t.top
+        e.above = trapezoids[0]
+        e.below = trapezoids[1]
+        return trapezoids
+
+    def case4(self, trapezoid, edge):
+        t = trapezoid; e = edge
+        lp = e.p if e.p.x == t.left_point.x else t.left_point
+        trapezoids = []
+        if self.tcross is t.top:
+            trapezoids.append(t.upper_left)
+            trapezoids[0].right_point = e.q
+        else:
+            trapezoids.append(Trapezoid(lp, e.q, t.top, e))
+            trapezoids[0].update_left(t.upper_left, e.above)
+        if self.bcross is t.bottom:
+            trapezoids.append(t.lower_left)
+            trapezoids[1].right_point = e.q
+        else:
+            trapezoids.append(Trapezoid(lp, e.q, e, t.bottom))
+            trapezoids[1].update_left(e.below, t.lower_left)
+        trapezoids.append(Trapezoid(e.q, t.right_point, t.top, t.bottom))
+        trapezoids[2].update_left_right(trapezoids[0], trapezoids[1], t.upper_right, t.lower_right)
+        return trapezoids
+  
+    def bounding_box(self, edges): 
+        margin = self.margin
+        max = edges[0].p + margin
+        min = edges[0].q - margin
+        for edge in edges:
+            if edge.p.x > max.x: max = Point(edge.p.x + margin, max.y)
+            if edge.p.y > max.y: max = Point(max.x, edge.p.y + margin)
+            if edge.q.x > max.x: max = Point(edge.q.x + margin, max.y)
+            if edge.q.y > max.y: max = Point(max.x, edge.q.y + margin)
+            if edge.p.x < min.x: min = Point(edge.p.x - margin, min.y)
+            if edge.p.y < min.y: min = Point(min.x, edge.p.y - margin)
+            if edge.q.x < min.x: min = Point(edge.q.x - margin, min.y)
+            if edge.q.y < min.y: min = Point(min.x, edge.q.y - margin)
+        top = Edge(Point(min.x, max.y), Point(max.x, max.y))
+        bottom = Edge(Point(min.x, min.y), Point(max.x, min.y))
+        left = bottom.p
+        right = top.q
+        return Trapezoid(left, right, top, bottom)
+
+class Node(object):
+
+    def __init__(self, left, right):
+        self.parent_list = []
+        self.left = left
+        self.right = right
+        if left != None: 
+            left.parent_list.append(self)
+        if right != None: 
+            right.parent_list.append(self)
+        
+    def replace(self, node):
+        for parent in node.parent_list:
+            if parent.left is node: 
+                parent.left = self
+            else: 
+                parent.right = self 
+        self.parent_list += node.parent_list
+    
+class Sink(Node):
+        
+    def __init__(self, trapezoid):
+        super(Sink, self).__init__(None, None)
+        self.trapezoid = trapezoid
+        trapezoid.sink = self
+        
+    def locate(self, edge): 
+        return self
+
+def isink(trapezoid):
+    if trapezoid.sink != None: 
+        return trapezoid.sink
+    return Sink(trapezoid)
+    
+class XNode(Node):
+
+    def __init__(self, point, lchild, rchild):
+        super(XNode, self).__init__(lchild, rchild)
+        self.point = point
+    
+    def locate(self, edge): 
+        if edge.p.x >= self.point.x: 
+            return self.right.locate(edge)
+        return self.left.locate(edge)
+
+class YNode(Node):
+    
+    def __init__(self, edge, lchild, rchild):
+        super(YNode, self).__init__(lchild, rchild)
+        self.edge = edge
+        
+    def locate(self, edge):
+        if self.edge.is_above(edge.p): 
+            return self.right.locate(edge)
+        elif self.edge.is_below(edge.p): 
+            return self.left.locate(edge)
+        else:
+            if edge.slope < self.edge.slope: 
+                return self.right.locate(edge)
+            else:
+                return self.left.locate(edge)
+            
+class QueryGraph:
+    
+    def __init__(self, head):
+        self.head = head
+        
+    def locate(self, edge):
+        return self.head.locate(edge).trapezoid
+  
+    def follow_edge(self, edge):
+        trapezoids = [self.locate(edge)]
+        j = 0
+        while(edge.q.x > trapezoids[j].right_point.x):
+            if edge.is_above(trapezoids[j].right_point):
+                trapezoids.append(trapezoids[j].upper_right)
+            else:
+                trapezoids.append(trapezoids[j].lower_right)
+            j += 1
+        return trapezoids
+  
+    def replace(self, sink, node):
+        if not sink.parent_list:
+          self.head = node
+        else:
+          node.replace(sink)
+
+    def case1(self, sink, edge, tlist):
+        yNode = YNode(edge, isink(tlist[1]), isink(tlist[2]))
+        qNode = XNode(edge.q, yNode, isink(tlist[3]))
+        pNode = XNode(edge.p, isink(tlist[0]), qNode)
+        self.replace(sink, pNode)
+  
+    def case2(self, sink, edge, tlist):
+        yNode = YNode(edge, isink(tlist[1]), isink(tlist[2]))
+        pNode = XNode(edge.p, isink(tlist[0]), yNode)
+        self.replace(sink, pNode)
+  
+    def case3(self, sink, edge, tlist):
+        yNode = YNode(edge, isink(tlist[0]), isink(tlist[1]))
+        self.replace(sink, yNode)
+
+    def case4(self, sink, edge, tlist):
+        yNode = YNode(edge, isink(tlist[0]), isink(tlist[1]))
+        qNode = XNode(edge.q, yNode, isink(tlist[2]))
+        self.replace(sink, qNode)
+
+
+PI_SLOP = 3.1
+
+class MonotoneMountain:
+
+    def __init__(self):
+        self.size = 0
+        self.tail = None
+        self.head = None
+        self.positive = False
+        self.convex_points = []
+        self.mono_poly = []
+        self.triangles = []
+        self.convex_polies = []
+        
+    def add(self, point):
+        if self.size == 0: 
+            self.head = point
+            self.size += 1
+        elif self.size == 1:
+            if point.not_equal(self.head):
+                self.tail = point
+                self.tail.prev = self.head
+                self.head.next = self.tail
+                self.size += 1
+        else:
+            if point.not_equal(self.tail):
+                self.tail.next = point
+                point.prev = self.tail
+                self.tail = point
+                self.size += 1
+
+    def remove(self, point):
+        next = point.next
+        prev = point.prev
+        point.prev.next = next
+        point.next.prev = prev
+        self.size -= 1
+
+    def process(self):
+        self.positive = self.angle_sign()
+        self.gen_mono_poly()
+        p = self.head.next
+        while p is not self.tail:
+            a = self.angle(p)
+            if a >= PI_SLOP or a <= -PI_SLOP: 
+                self.remove(p)
+            elif self.is_convex(p): 
+                self.convex_points.append(p)
+            p = p.next
+        self.triangulate()
+
+    def triangulate(self):
+        while len(self.convex_points) > 0:
+            ear = self.convex_points.pop(0)
+            a = ear.prev
+            b = ear
+            c = ear.next
+            triangle = [a, b, c]
+            self.triangles.append(triangle)
+            self.remove(ear)
+            if self.valid(a): self.convex_points.append(a)
+            if self.valid(c): self.convex_points.append(c)
+        assert(self.size <= 3, "Triangulation bug, please report")
+
+    def valid(self, p):
+        return p != self.head and p != self.tail and self.is_convex(p)
+
+    def gen_mono_poly(self): 
+        p = self.head
+        while(p != None):
+            self.mono_poly.append(p)
+            p = p.next
+
+    def angle(self, p):
+        a = p.next - p
+        b = p.prev - p
+        return atan2(a.cross(b), a.dot(b))
+
+    def angle_sign(self):
+        a = self.head.next - self.head
+        b = self.tail - self.head
+        return atan2(a.cross(b), a.dot(b)) >= 0
+
+    def is_convex(self, p):
+        if self.positive != (self.angle(p) >= 0): return False
+        return True