diff --git a/run_all.py b/run_all.py
index 4d047f7..feace41 100755
--- a/run_all.py
+++ b/run_all.py
@@ -159,23 +159,26 @@ else:
     hull_edges = list(utils.tour(hull))
     LOGN( "\t\tHull of",len(hull_edges),"edges" )
 
-    LOGN( "\tRemove triangles that have at least one edge in common with the convex hull" )
-    def adjoin_hull(triangle):
-        """Return True if the given triangle has at least one edge that is in the set hull_edges."""
-        for (p,q) in utils.tour(list(triangle)):
-            if (p,q) in hull_edges or (q,p) in hull_edges:
-                return True
-        return False
+    LOGN( "\tRemove triangles that are not sub-parts of the Penrose tiling" )
+    # def adjoin_hull(triangle):
+    #     """Return True if the given triangle has at least one edge that is in the set hull_edges."""
+    #     for (p,q) in utils.tour(list(triangle)):
+    #         if (p,q) in hull_edges or (q,p) in hull_edges:
+    #             return True
+    #     return False
 
     def acute_triangle(triangle):
         """Return True if the center of the circumcircle of the given triangle lies inside the triangle.
            That is if the triangle is acute."""
         return triangulation.in_triangle( triangulation.circumcircle(triangle)[0], triangle )
 
+    # FIXME at depth 3, some triangles have an edge in the convex hull...
     # Filter out edges that are in hull_edges
-    tri_nohull = list(filter_if_not( adjoin_hull, triangles ))
+    # tri_nohull = list(filter_if_not( adjoin_hull, triangles ))
     # Filter out triangles that are obtuse
-    triangulated = list(filter_if_not( acute_triangle, tri_nohull ))
+    # triangulated = list(filter_if_not( acute_triangle, tri_nohull ))
+
+    triangulated = list(filter_if_not( acute_triangle, triangles ))
     LOGN( "\t\tRemoved", len(triangles)-len(triangulated), "triangles" )
 
     triangulation_edges = triangulation.edges_of( triangulated )
@@ -189,15 +192,18 @@ else:
 ########################################################################
 
 if args.voronoi:
-    voronoi_centers = utils.load_points(args.voronoi)
+    # voronoi_centers = utils.load_points(args.voronoi)
+    pass
 
 else:
-    LOGN( "Compute the nodes of the Voronoï diagram" )
-    voronoi_centers = voronoi.centers(triangulated)
+    # LOGN( "Compute the nodes of the Voronoï diagram" )
+    voronoi_graph = voronoi.dual( triangulated )
+    voronoi_edges = voronoi.edges_of( voronoi_graph )
+    voronoi_centers = voronoi_graph.keys()
 
-    with open("d%i_voronoi_centers.points" % depth, "w") as fd:
-        for p in voronoi_centers:
-            fd.write( "%f %f\n" % (p[0],p[1]) )
+    # with open("d%i_voronoi_centers.points" % depth, "w") as fd:
+    #     for p in voronoi_centers:
+    #         fd.write( "%f %f\n" % (p[0],p[1]) )
 
 
 ########################################################################
@@ -238,8 +244,9 @@ uberplot.scatter_segments( ax, penrose_segments, edgecolor=tcol, alpha=0.9, line
 # triangulation
 uberplot.plot_segments( ax, triangulation_edges, edgecolor="green", alpha=0.2, linewidth=1 )
 
-# Voronoï centers
-uberplot.scatter_points( ax, voronoi_centers, edgecolor="none", facecolor="green", linewidth=0 )
+# Voronoï
+uberplot.scatter_points( ax, voronoi_centers, edgecolor="magenta", facecolor="white", s=200, alpha=0.5 )
+uberplot.plot_segments( ax, voronoi_edges, edgecolor="magenta", alpha=0.2, linewidth=1 )
 
 ax.set_aspect('equal')
 
diff --git a/uberplot.py b/uberplot.py
index 59805b8..9264231 100644
--- a/uberplot.py
+++ b/uberplot.py
@@ -39,7 +39,8 @@ def scatter_segments( ax, segments, **kwargs  ):
 def scatter_points( ax, points, **kwargs  ):
     x = [i[0] for i in points]
     y = [i[1] for i in points]
-    ax.scatter( x,y, s=20, marker='o', **kwargs) 
+    # ax.scatter( x,y, s=20, marker='o', **kwargs) 
+    ax.scatter( x,y, marker='o', **kwargs) 
 
 
 if __name__=="__main__":
diff --git a/voronoi.py b/voronoi.py
new file mode 100644
index 0000000..fdf979d
--- /dev/null
+++ b/voronoi.py
@@ -0,0 +1,118 @@
+#/usr/bin/env python
+
+from utils import tour,LOG,LOGN,x,y
+import triangulation
+
+def nodes( triangles ):
+    """Compute the locations of the centers of all the circumscribed circles of the given triangles"""
+    for triangle in triangles:
+        (cx,cy),r = triangulation.circumcircle(triangle)
+        yield (cx,cy)
+
+
+def edge_in( edge, edges ):
+    """Return True if the given edge (or its symetric) is in the given polygon."""
+    n,m = edge
+    assert( len(n) == 2 )
+    assert( len(m) == 2 )
+    return (n,m) in edges or (m,n) in edges
+
+
+def neighbours( triangle, polygons ):
+    """Returns the set of triangles in candidates that have an edge in common with the given triangle."""
+    for polygon in polygons:
+        if polygon == triangle:
+            continue
+
+        # Convert list of points to list of edges, because we want to find EDGE neighbours.
+        edges_poly = list(tour(list(polygon )))
+        assert( len(list(edges_poly)) > 0 )
+        edges_tri  = list(tour(list(triangle)))
+        assert( len(list(edges_tri)) > 0 )
+
+        # If at least one of the edge that are in availables polygons are also in the given triangle.
+        # Beware the symetric edges.
+        # if any( edge_in( edge, edges_tri ) for edge in edges_poly ):
+        #     yield polygon
+        for edge in edges_poly:
+            if edge_in( edge, edges_tri ):
+                yield polygon
+                break
+
+
+def dual( triangles ):
+    graph = {}
+
+    def add_edge( current, neighbor ):
+        if current in graph:
+            if neighbor not in graph[current]:
+                graph[current].append( neighbor )
+        else:
+            graph[current] = [ neighbor ]
+
+    for triangle in triangles:
+        assert( len(triangle) == 3 )
+        current_node = triangulation.circumcircle(triangle)[0]
+        assert( len(current_node) == 2 )
+
+        for neighbor_triangle in neighbours( triangle, triangles ):
+            neighbor_node = triangulation.circumcircle(neighbor_triangle)[0]
+            assert( len(neighbor_node) == 2 )
+
+            add_edge(  current_node, neighbor_node )
+            add_edge( neighbor_node,  current_node )
+
+    return graph
+
+
+def edges_of( graph ):
+    # edges = set()
+    edges = []
+    for k in graph:
+        for n in graph[k]:
+            if k != n and (k,n) not in edges and (n,k) not in edges:
+                # edges.add( (k,n) )
+                edges.append( (k,n) )
+
+    return edges
+
+
+if __name__ == "__main__":
+    import sys
+    import random
+    import utils
+    import uberplot
+    import triangulation
+    import matplotlib.pyplot as plot
+
+    if len(sys.argv) > 1:
+        scale = 100
+        nb = int(sys.argv[1])
+        points = [ (scale*random.random(),scale*random.random()) for i in range(nb)]
+    else:
+        points = [
+                (0,40),
+                (100,60),
+                (40,0),
+                (50,100),
+                (90,10),
+                # (50,50),
+                ]
+
+    fig = plot.figure()
+
+    triangles = triangulation.delaunay_bowyer_watson( points )
+    delaunay_edges = triangulation.edges_of( triangles )
+
+    voronoi_graph = dual( triangles )
+    voronoi_edges = edges_of( voronoi_graph )
+    print voronoi_edges
+
+    ax = fig.add_subplot(111)
+    ax.set_aspect('equal')
+    uberplot.scatter_segments( ax, delaunay_edges, facecolor = "blue" )
+    uberplot.plot_segments( ax, delaunay_edges, edgecolor = "blue" )
+    uberplot.scatter_segments( ax, voronoi_edges, facecolor = "red" )
+    uberplot.plot_segments( ax, voronoi_edges, edgecolor = "red" )
+    plot.show()
+