diff --git a/triangulation.py b/triangulation.py
new file mode 100644
index 0000000..4fcba28
--- /dev/null
+++ b/triangulation.py
@@ -0,0 +1,338 @@
+
+import sys
+import math
+from utils import tour,LOG,LOGN
+from itertools import ifilterfalse as filter_if_not
+
+# Based on http://paulbourke.net/papers/triangulate/
+# Efficient Triangulation Algorithm Suitable for Terrain Modelling
+#     An Algorithm for Interpolating Irregularly-Spaced Data
+#     with Applications in Terrain Modelling
+# Written by Paul Bourke
+# Presented at Pan Pacific Computer Conference, Beijing, China.
+# January 1989
+
+def x( point ):
+    return point[0]
+
+def y( point ):
+    return point[1]
+
+def mid( xy, pa, pb ):
+    return ( xy(pa) + xy(pb) ) / 2.0
+
+def middle( pa, pb ):
+    return mid(x,pa,pb),mid(y,pa,pb)
+
+def mtan( pa, pb ):
+    return -1 * ( x(pa) - x(pb) ) / ( y(pa) - y(pb) )
+
+
+class CoincidentPointsError(Exception):
+    """Coincident points"""
+    pass
+
+def circumcircle( triangle, epsilon = sys.float_info.epsilon ):
+    """Compute the circumscribed circle of a triangle and 
+    Return a 2-tuple: ( (center_x, center_y), radius )"""
+
+    assert( len(triangle) == 3 )
+    p0,p1,p2 = triangle
+    assert( len(p0) == 2 )
+    assert( len(p1) == 2 )
+    assert( len(p2) == 2 )
+
+    dy01 = abs( y(p0) - y(p1) )
+    dy12 = abs( y(p1) - y(p2) )
+
+    if dy01 < epsilon and dy12 < epsilon:
+        # coincident points
+        raise CoincidentPointsError
+
+    elif dy01 < epsilon:
+        m12 = mtan( p2,p1 )
+        mx12,my12 = middle( p1, p2 )
+        cx = mid( x, p1, p0 )
+        cy = m12 * (cx - mx12) + my12
+
+    elif dy12 < epsilon:
+        m01 = mtan( p1, p0 )
+        mx01,my01 = middle( p0, p1 )
+        cx = mid( x, p2, p1 )
+        cy = m01 * ( cx - mx01 ) + my01
+
+    else:
+        m01 =  mtan( p1, p0 )
+        m12 =  mtan( p2, p1 )
+        mx01,my01 = middle( p0, p1 )
+        mx12,my12 = middle( p1, p2 )
+        cx = ( m01 * mx01 - m12 * mx12 + my12 - my01 ) / ( m01 - m12 )
+        if dy01 > dy12:
+            cy = m01 * ( cx - mx01 ) + my01
+        else:
+            cy = m12 * ( cx - mx12 ) + my12
+
+    dx1 = x(p1) - cx
+    dy1 = y(p1) - cy
+    r = math.sqrt(dx1**2 + dy1**2)
+
+    return (cx,cy),r
+
+
+def in_circle( p, center, radius, epsilon  = sys.float_info.epsilon ):
+    """Return True if the given point p is in the circumscribe circle of the given triangle"""
+
+    assert( len(p) == 2 )
+    cx,cy = center
+
+    dxp = x(p) - cx
+    dyp = y(p) - cy
+    dr = math.sqrt(dxp**2 + dyp**2)
+
+    if (dr - radius) <= epsilon:
+        return True
+    else:
+        return False
+
+
+def in_circumcircle( p, triangle, epsilon  = sys.float_info.epsilon ):
+    """Return True if the given point p is in the circumscribe circle of the given triangle"""
+
+    assert( len(p) == 2 )
+    (cx,cy),r = circumcircle( triangle, epsilon )
+    return in_circle( p, (cx,cy), r, epsilon )
+
+
+def bounds( vertices ):
+    # find vertices set bounds
+    xmin = x(vertices[0])
+    ymin = y(vertices[0])
+    xmax = xmin
+    ymax = ymin
+
+    # we do not use min(vertices,key=x) because it would iterate 4 times over the list, instead of just one
+    for v in vertices:
+        xmin = min(x(v),xmin)
+        xmax = max(x(v),xmax)
+        ymin = min(y(v),ymin)
+        ymax = max(y(v),ymax)
+    return (xmin,ymin),(xmax,ymax)
+
+
+def edges_of( triangulation ):
+    edges = []
+    for t in triangulation:
+        for e in utils.tour(list(t)):
+            edges.append( e )
+    return edges
+
+
+def delaunay_bowyer_watson( points, epsilon = sys.float_info.epsilon, supert=20, do_plot = True ):
+
+    if do_plot and len(points) > 10:
+        print "WARNING it is a bad idea to plot each steps of a triangulation of many points"
+        return []
+
+    # sort points first on the x-axis, then on the y-axis
+    vertices = sorted( points )
+
+    (xmin,ymin),(xmax,ymax) = bounds( vertices )
+
+    dx = xmax - xmin
+    dy = ymax - ymin
+    dmax = max( dx, dy )
+    xmid = (xmax + xmin) / 2.0
+    ymid = (ymax + ymin ) / 2.0
+
+
+    # compute the super triangle, that encompasses all the vertices
+    supertri = ( (xmid-supert*dmax, ymid-dmax ),
+                 (xmid, ymid+supert*dmax),
+                 (xmid+supert*dmax, ymid-dmax) )
+
+    LOGN( "super-triangle",supertri )
+
+    # it is the first triangle of the list
+    triangles = [ supertri ]
+
+    completed = { supertri: False }
+
+    # The predicate returns true if at least one of the vertices
+    # is also found in the supertriangle
+    def match_supertriangle( tri ):
+        if tri[0] in supertri or \
+           tri[1] in supertri or \
+           tri[2] in supertri:
+            return True
+
+    # insert vertices one by one
+    it=0
+    for vi,vertex in enumerate(vertices):
+        LOGN( "\tvertex",vertex )
+        assert( len(vertex) == 2 )
+
+        if do_plot:
+            fig = plot.figure()
+            ax = fig.add_subplot(111)
+            scatter_x = [ p[0] for p in vertices[:vi]+list(supertri)]
+            scatter_y = [ p[1] for p in vertices[:vi]+list(supertri)]
+            ax.scatter( scatter_x,scatter_y, s=30, marker='o', facecolor="black")
+            ax.scatter( vertex[0],vertex[1], s=30, marker='o', facecolor="red")
+            uberplot.plot_segments( ax, edges_of(triangles), edgecolor = "blue", alpha=0.3, linestyle='dashed' )
+
+        # All the triangles whose circumcircle encloses the point to be added are identified,
+        # the outside edges of those triangles form an enclosing polygon.
+
+        # forget previous candidate polygon's edges
+        enclosing = []
+
+        removed = []
+        for triangle in triangles:
+            LOGN( "\t\ttriangle",triangle )
+            assert( len(triangle) == 3 )
+
+            # if completed has a key, test it, else return False
+            if completed.get( triangle, False ):
+                LOGN( "\t\t\tAlready completed" )
+                # if do_plot:
+                    # uberplot.plot_segments( ax, tour(list(triangle)), edgecolor = "magenta", alpha=1, lw=1, linestyle='dotted' )
+                continue
+
+            LOGN( "\t\t\tCircumcircle" ) 
+            assert( triangle[0] != triangle[1] and triangle[1] != triangle [2] and triangle[2] != triangle[0] )
+            center,radius = circumcircle( triangle, epsilon )
+
+            # if it match Delaunay's conditions
+            if x(center) < x(vertex) and math.sqrt((x(vertex)-x(center))**2) > radius:
+                LOGN( "\t\t\tMatch Delaunay, mark as completed" ) 
+                completed[triangle] = True
+
+            # if the current vertex is inside the circumscribe circle of the current triangle
+            # add the current triangle's edges to the candidate polygon
+            if in_circle( vertex, center, radius, epsilon ):
+                LOGN( "\t\t\tIn circumcircle, add to enclosing polygon",triangle )
+                if do_plot:
+                    # if not match_supertriangle( triangle ):
+                    circ = plot.Circle(center, radius, facecolor='yellow', edgecolor="orange", alpha=0.1)
+                    ax.add_patch(circ)
+
+                for p0,p1 in tour(list(triangle)):
+                    # then add this edge to the polygon enclosing the vertex
+                    enclosing.append( (p0,p1) )
+                # and remove the corresponding triangle from the current triangulation
+                removed.append( triangle )
+                completed.pop(triangle,None)
+
+            elif do_plot:
+                # if not match_supertriangle( triangle ):
+                circ = plot.Circle(center, radius, facecolor='lightgrey', edgecolor="grey", alpha=0.1)
+                ax.add_patch(circ)
+
+
+        # end for triangle in triangles
+
+        # The triangles in the enclosing polygon are deleted and
+        # new triangles are formed between the point to be added and
+        # each outside edge of the enclosing polygon. 
+
+        # actually remove triangles
+        for triangle in removed:
+            # if do_plot:
+                # if not match_supertriangle( triangle ):
+                # uberplot.plot_segments( ax, tour(list(triangle)), edgecolor = "orange", alpha=0.3, lw=2 )
+
+            triangles.remove(triangle)
+
+
+        # remove duplicated edges
+        # this leaves the edges of the enclosing polygon only,
+        # because enclosing edges are only in a single triangle,
+        # but edges inside the polygon are at least in two triangles.
+        # duplicated = []
+        # for i,ei in enumerate(enclosing):
+        #     for j,ej in enumerate(enclosing,i+1):
+        #         if (ei[0] == ej[1] and ei[1] == ej[0]) or (ei[0] == ej[0] and ei[1] == ej[1]):
+        #             duplicated.append( ei )
+        # for e in duplicated:
+        #     enclosing.remove(e)
+        hull = []
+        for i,(p0,p1) in enumerate(enclosing):
+            if (p0,p1) not in enclosing[i+1:] and (p1,p0) not in enclosing:
+                hull.append((p0,p1))
+
+        if do_plot:
+            uberplot.plot_segments( ax, hull, edgecolor = "red", alpha=1, lw=1, linestyle='solid' )
+
+
+        # create new triangles using the current vertex and the enclosing hull
+        # All candidates should be arranged in clockwise order!
+        LOGN( "\t\tCreate new triangles" )
+        for p0,p1 in hull:
+            assert( p0 != p1 )
+            # if p0 != vertex and p1 != vertex:
+            # triangle = tuple(sorted([p0,p1,vertex]))
+            triangle = tuple([p0,p1,vertex])
+            LOGN("\t\t\tNew triangle",triangle)
+            triangles.append( triangle )
+            completed[triangle] = False
+
+            if do_plot: # linestyle = ['solid' | 'dashed' | 'dashdot' | 'dotted']
+                uberplot.plot_segments( ax, tour(list(triangle)), edgecolor = "green", alpha=0.3, linestyle='solid' )
+
+
+        with open("triangulation_%i.dat" % it, 'w') as fd:
+            for triangle in triangles:
+                for edge in tour(list(triangle)):
+                    coords = tuple([coord for point in edge for coord in point])
+                    fd.write( "%f %f %f %f\n" % coords )
+
+        if do_plot:
+            # ax.set_ylim([-100,200])
+            # ax.set_xlim([-100,200])
+            plot.savefig("triangulation_%i.png" % it, dpi=300)
+            plot.close()
+
+        it+=1
+
+    # end for vertex in vertices
+
+
+    # Remove triangles that have at least one of the supertriangle vertices
+    LOGN( "\tRemove super-triangles" ) 
+
+    # filter out elements for which the predicate is False
+    # here: *keep* elements that *do not* have a common vertex
+    triangulation = filter_if_not( match_supertriangle, triangles )
+
+    return triangulation
+
+
+
+if __name__ == "__main__":
+    import random
+    import utils
+    import uberplot
+    import matplotlib.pyplot as plot
+    from matplotlib.path import Path
+    import matplotlib.patches as patches
+
+    scale = 100
+    nb = 10
+    points = [ (scale*random.random(),scale*random.random()) for i in range(nb)]
+    # points = [
+    #         (0,40),
+    #         (100,60),
+    #         (40,0),
+    #         (50,100),
+    #         (90,10),
+    #         ]
+
+    triangles = delaunay_bowyer_watson( points, epsilon=10e-4, supert=3 )
+
+    edges = edges_of( triangles )
+
+    fig = plot.figure()
+    ax = fig.add_subplot(111)
+    uberplot.scatter_segments( ax, edges, facecolor = "red" )
+    uberplot.plot_segments( ax, edges, edgecolor = "blue", alpha=0.2 )
+    plot.show()