added fast robust predicates

python/framework/framework.pyx

@@ -5,6 +5,7 @@ from math import pi as PI
 
 from gl cimport *
 include "polydecomp.pxi"
+include "seidel.pxi"
 
 cdef extern from 'math.h':
     double cos(double)

python/framework/polydecomp.pxi

@@ -4,15 +4,18 @@
 ##
 from sys import float_info
 
-def makeCCW(list poly):
+cdef extern from 'predicates.h':
+    double orient2d(double *pa, double *pb, double *pc)
+    
+def make_ccw(list poly):
     cdef int br = 0
     # find bottom right point
     for i from 1 <= i < len(poly):
         if poly[i][1] < poly[br][1] or (poly[i][1] == poly[br][1] and poly[i][0] > poly[br][0]):
-            br = i
+            br = i       
     # reverse poly if clockwise
     if not left(at(poly, br - 1), at(poly, br), at(poly, br + 1)):
-        poly.reverse()
+        poly.reverse()    
 
 cpdef list decompose_poly(list poly, list polys):
 
@@ -48,8 +51,8 @@ cpdef list decompose_poly(list poly, list polys):
 
             # if there are no vertices to connect to, choose a point in the middle
             if lower_index == (upper_index + 1) % len(poly):
-                p[0] = (lowerInt[0] + upperInt[0]) / 2
-                p[1] = (lowerInt[1] + upperInt[1]) / 2
+                p[0] = (lowerInt[0] + upperInt[0]) * 0.5
+                p[1] = (lowerInt[1] + upperInt[1]) * 0.5
 
                 if i < upper_index:
                     lower_poly.extend(poly[i:upper_index+1])
@@ -99,26 +102,19 @@ cpdef list decompose_poly(list poly, list polys):
             else:
                 decompose_poly(upper_poly, polys)
                 decompose_poly(lower_poly, polys)
-            
             return
 
     polys.append(poly)
     
 cdef list intersection(list p1, list p2, list q1, list q2):
-    cdef list i = []
-    cdef float a1, b1, c1, a2, b2, c2, det
-    a1 = p2[1] - p1[1]
-    b1 = p1[0] - p2[0]
-    c1 = a1 * p1[0] + b1 * p1[1]
-    a2 = q2[1] - q1[1]
-    b2 = q1[0] - q2[0]
-    c2 = a2 * q1[0] + b2 * q1[1]
-    det = a1 * b2 - a2 * b1
-    if not eq(det, 0): 
-        # lines are not parallel
-        i.append((b2 * c1 - b1 * c2) / det)
-        i.append((a1 * c2 - a2 * c1) / det)
-    return i
+    cdef double pqx, pqy, bax, bay, t
+    pqx = p1[0] - p2[0]
+    pqy = p1[1] - p2[1]
+    t = pqy*(q1[0]-p2[0]) - pqx*(q1[1]-p2[1])
+    t /= pqx*(q2[1]-q1[1]) - pqy*(q2[0]-q1[0])
+    bax = t*(q2[0]-q1[0]) + q1[0]
+    bay = t*(q2[1]-q1[1]) + q1[1]
+    return [bax, bay]
 
 cdef bool eq(float a, float b):
     return abs(a - b) <= 1e-8
@@ -130,16 +126,28 @@ cdef float area(list a, list b, list c):
     return (((b[0] - a[0])*(c[1] - a[1]))-((c[0] - a[0])*(b[1] - a[1])))
 
 cdef bool left(list a, list b, list c):
-    return area(a, b, c) > 0
+    cdef double *x = [a[0], a[1]]
+    cdef double *y = [b[0], b[1]]
+    cdef double *z = [c[0], c[1]]
+    return orient2d(x, y, z) > 0.0
 
 cdef bool leftOn(list a, list b, list c):
-    return area(a, b, c) >= 0
+    cdef double *x = [a[0], a[1]]
+    cdef double *y = [b[0], b[1]]
+    cdef double *z = [c[0], c[1]]
+    return orient2d(x, y, z) >= 0.0
 
 cdef bool right(list a, list b, list c):
-    return area(a, b, c) < 0
+    cdef double *x = [a[0], a[1]]
+    cdef double *y = [b[0], b[1]]
+    cdef double *z = [c[0], c[1]]
+    return orient2d(x, y, z) < 0.0
 
 cdef bool rightOn(list a, list b, list c):
-    return area(a, b, c) <= 0
+    cdef double *x = [a[0], a[1]]
+    cdef double *y = [b[0], b[1]]
+    cdef double *z = [c[0], c[1]]
+    return orient2d(x, y, z) <= 0.0
 
 cdef float sqdist(list a, list b):
     cdef float dx = b[0] - a[0]

python/framework/predicates.h

@@ -0,0 +1,4 @@
+double orient2d(pa, pb, pc);
+double *pa;
+double *pb;
+double *pc;

python/framework/seidel.pxi

@@ -0,0 +1,636 @@
+#
+# 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
+
+##
+## Based on Raimund Seidel'e paper "A simple and fast incremental randomized
+## algorithm for computing trapezoidal decompositions and for triangulating polygons"
+## (Ported from poly2tri)
+##
+
+# Shear transform. May effect numerical robustness
+SHEAR = 1e-6
+
+cdef extern from 'math.h':
+    double atan2(double, double)
+    
+cdef extern from 'predicates.h':
+    double orient2d(double *pa, double *pb, double *pc)
+    
+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 neq(self, other):
+        return other.x != self.x or other.y != self.y
+        
+    def clone(self):
+        return Point(self.x, self.y)
+
+class Edge(object):
+    
+    def __init__(self, p, q):
+        self.p = p
+        self.q = q
+        self.slope = (q.y - p.y) / (q.x - p.x)
+        self.b = p.y - (p.x * self.slope)
+        self.above, self.below = None, None
+        self.mpoints = []
+        self.mpoints.append(p)
+        self.mpoints.append(q)
+    
+    ##
+    ## NOTE Rounding accuracy significantly effects numerical robustness!!!
+    ##
+    
+    def is_above(self, point):
+        cdef double *a = [self.p.x, self.p.y]
+        cdef double *b = [self.q.x, self.q.y]
+        cdef double *c = [point.x, point.y]
+        return orient2d(a, b, c) < 0
+        
+    def is_below(self, point):
+        cdef double *a = [self.p.x, self.p.y]
+        cdef double *b = [self.q.x, self.q.y]
+        cdef double *c = [point.x, point.y]
+        return orient2d(a, b, c) > 0
+
+    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.0 and a2 != 0.0 and (a1 * a2) < 0.0:
+            a3 = self.signed_area(c, d, a)
+            a4 = a3 + a2 - a1
+            if a3 * a4 < 0.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
+        self.key = hash(self)
+        
+    def update_left(self, ul, ll):
+        self.upper_left = ul
+        if ul != None: ul.upper_right = self
+        self.lower_left = ll
+        if ll != None: ll.lower_right = self
+          
+    def update_right(self, ur, lr):
+        self.upper_right = ur
+        if ur != None: ur.upper_left = self
+        self.lower_right = lr
+        if lr != None: lr.lower_left = self 
+          
+    def update_left_right(self, ul, ll, ur, lr):
+        self.upper_left = ul
+        if ul != None: ul.upper_right = self
+        self.lower_left = ll
+        if ll != None: ll.lower_right = self
+        self.upper_right = ur
+        if ur != None: ur.upper_left = self
+        self.lower_right = lr
+        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):
+        v1 = line_intersect(self.top, self.left_point.x)
+        v2 = line_intersect(self.bottom, self.left_point.x)
+        v3 = line_intersect(self.bottom, self.right_point.x)
+        v4 = line_intersect(self.top, self.right_point.x)
+        return v1, v2, v3, v4
+  
+    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 x, y
+    
+class Triangulator(object):
+     
+    def __init__(self, poly_line):
+        assert len(poly_line) > 3, "Number of points must be > 3"
+        self.polygons = []
+        self.trapezoids = []
+        self.xmono_poly = []
+        self.edge_list = self.init_edges(poly_line)
+        self.trapezoidal_map = TrapezoidalMap()
+        self.bounding_box = self.trapezoidal_map.bounding_box(self.edge_list)
+        self.query_graph = QueryGraph(isink(self.bounding_box))
+            
+        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:
+                # Remove old trapezods
+                del self.trapezoidal_map.map[t.key]
+                # Bisect old trapezoids and create new
+                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[t.key] = t
+            self.trapezoidal_map.clear()
+                    
+        # Mark outside trapezoids w/ depth-first search
+        for k, t in self.trapezoidal_map.map.items():
+            self.mark_outside(t)
+            
+        # Collect interior trapezoids
+        for k, t in self.trapezoidal_map.map.items():
+            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) > 2:
+                mountain = MonotoneMountain()
+                points = merge_sort(edge.mpoints)
+                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):
+        edge_list = []
+        size = len(points)
+        for i in range(size):
+            j = i + 1 if i < size-1 else 0
+            p = points[i][0], points[i][1]
+            q = points[j][0], points[j][1]
+            edge_list.append((p, q))
+        return self.order_edges(edge_list)
+  
+    def order_edges(self, edge_list):
+        edges = []
+        for e in edge_list:
+            p = shear_transform(e[0])
+            q = shear_transform(e[1])
+            if p.x > q.x: 
+                edges.append(Edge(q, p))
+            else: 
+                edges.append(Edge(p, q))
+        # Randomized incremental algorithm
+        shuffle(edges)
+        return edges
+
+def shear_transform(point):
+    return Point(point[0] + SHEAR * point[1], point[1])
+ 
+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):
+
+    def __init__(self):
+        self.map = {}
+        self.margin = 50.0
+        self.bcross = None
+        self.tcross = None
+
+    def clear(self):
+        self.bcross = None
+        self.tcross = None
+        
+    def case1(self, t, e):
+        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, t, 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))
+        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, t, 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 = []
+        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, t, e):
+        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 e in edges:
+            if e.p.x > max.x: max = Point(e.p.x + margin, max.y)
+            if e.p.y > max.y: max = Point(max.x, e.p.y + margin)
+            if e.q.x > max.x: max = Point(e.q.x + margin, max.y)
+            if e.q.y > max.y: max = Point(max.x, e.q.y + margin)
+            if e.p.x < min.x: min = Point(e.p.x - margin, min.y)
+            if e.p.y < min.y: min = Point(min.x, e.p.y - margin)
+            if e.q.x < min.x: min = Point(e.q.x - margin, min.y)
+            if e.q.y < min.y: min = Point(min.x, e.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 = top.p
+        right = top.q
+        trap = Trapezoid(left, right, top, bottom)
+        self.map[trap.key] = trap
+        return trap
+        
+class Node(object):
+
+    def __init__(self, lchild, rchild):
+        self.parent_list = []
+        self.lchild = lchild
+        self.rchild = rchild
+        if lchild != None: 
+            lchild.parent_list.append(self)
+        if rchild != None: 
+            rchild.parent_list.append(self)
+        
+    def replace(self, node):
+        for parent in node.parent_list:
+            if parent.lchild is node: 
+                parent.lchild = self
+            else: 
+                parent.rchild = 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 is None: 
+        return Sink(trapezoid)
+    return trapezoid.sink
+    
+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.rchild.locate(edge)
+        return self.lchild.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.rchild.locate(edge)
+        if self.edge.is_below(edge.p): 
+            return self.lchild.locate(edge)
+        if edge.slope < self.edge.slope: 
+            return self.rchild.locate(edge)
+        return self.lchild.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)]
+        while(edge.q.x > trapezoids[-1].right_point.x):
+            if edge.is_above(trapezoids[-1].right_point):
+                trapezoids.append(trapezoids[-1].upper_right)
+            else:
+                trapezoids.append(trapezoids[-1].lower_right)
+        return trapezoids
+  
+    def replace(self, sink, node):
+        if sink.parent_list:
+            node.replace(sink)
+        else:
+            self.head = node
+
+    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 is 0: 
+            self.head = point
+            self.size = 1
+        elif self.size is 1:
+            if point.neq(self.head):
+                self.tail = point
+                self.tail.prev = self.head
+                self.head.next = self.tail
+                self.size = 2
+        else:
+            if point.neq(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 != self.tail:
+            a = self.angle(p)
+            if a >= PI_SLOP or a <= -PI_SLOP or a == 0: 
+                self.remove(p)
+            elif self.is_convex(p): 
+                self.convex_points.append(p)
+            p = p.next
+        self.triangulate()
+
+    def triangulate(self):
+        while self.convex_points:
+            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