Moved everything into a package that can be easly installed.
Moved everything into a package that can be easly installed.

--- a/test-association.py
+++ b/test-association.py
@@ -3,7 +3,7 @@
 import random
 import time
 
-import unittests
+from voronoi import unittests
 
 # Config 
 width=1000

--- a/test-random.py
+++ b/test-random.py
@@ -3,13 +3,13 @@
 import random
 import time
 
-import unittests
+from voronoi import unittests
 
 # Config 
 width=1000
 height=1000
 npoints=1000
-#seed=1341469547
+#seed=1
 seed=int(time.time())
 
 # Voronoi

--- a/test-verticaledge.py
+++ b/test-verticaledge.py
@@ -3,7 +3,7 @@
 import numpy
 import time
 
-import unittests
+from voronoi import unittests
 
 # Config 
 width=500

--- a/test-verticalmirror.py
+++ b/test-verticalmirror.py
@@ -3,7 +3,7 @@
 import random
 import time
 
-import unittests
+from voronoi import unittests
 
 # Config 
 width=500

file:a/tools.py (deleted)
--- a/tools.py
+++ /dev/null
@@ -1,57 +1,1 @@
-def img_coords(edge,width,height):
-    intersections=[]
 
-    # Intersection with north and south edges (horizontal)
-    if edge.direction[1]!=0:
-        # North (y=height)
-        t=(height-edge.r[1])/edge.direction[1]
-        x=edge.value_x(t)
-        if 0<x and x<width:
-            intersections.append(t)
-        # South (y=0)
-        t=-edge.r[1]/edge.direction[1]
-        x=edge.value_x(t)
-        if 0<x and x<width:
-            intersections.append(t)
-
-    # Intersections with east and west (vertical)
-    if edge.direction[0]!=0:
-        # East (x=height)
-        t=(width-edge.r[0])/edge.direction[0]
-        y=edge.value_y(t)
-        if 0<y and y<height:
-            intersections.append(t)
-        # West (x=0)
-        t=-edge.r[0]/edge.direction[0]
-        y=edge.value_y(t)
-        if 0<y and y<height:
-            intersections.append(t)
-
-    if len(intersections)==0:
-        return None # Line misses image completely
-
-    # Figure out the range of values the line takes inside the box
-    mint=min(intersections)
-    maxt=max(intersections)
-
-    # Determine what range of values the line defined as now takes
-    if edge.start==None: # Make sure to clip the edge if it is infinite
-        edge.start=mint # This is when the box is left
-
-    if edge.end==None: # Make sure to clip the edge if it is infinite
-        edge.end=maxt # This is when the box is left
-
-    p1t=edge.start
-    p2t=edge.end
-
-    if p2t<mint or p1t>maxt: # This edge lies outside of the box
-        return None
-
-    # Determine the points that the edge takes after being clipped by the box
-    p1t=max(p1t,mint)
-    p2t=min(p2t,maxt)
-
-    start=edge.value(p1t)
-    end=edge.value(p2t)
-    return ((int(start[0]),int(start[1])),(int(end[0]),int(end[1])))
-

file:a/unittests.py (deleted)
--- a/unittests.py
+++ /dev/null
@@ -1,74 +1,1 @@
-import numpy
-import voronoi
-import cv,cv2
-import tools
-import random
 
-# Useful functions
-def toint(t2):
-    return (int(t2[0]),int(t2[1]))
-
-# Run a test on a set of points. The test does the following:
-#  1) Creates a CV window named "testwin"
-#  2) Draw the points into it as tiny circles
-#  3) Draws the edges calculated by voronoi
-#  4) Waits for a keystroke
-# The testwin is created based on passed width and height
-def run_test(points,width,height):
-    # Init open CV
-    cv.NamedWindow('testwin')
-
-    img = numpy.zeros((height,width, 1), numpy.uint8)
-    
-    for p in points:
-        cv2.circle(img,toint(p),2,255,-1)
-
-    v=voronoi.Voronoi(points,width,height)
-    v.img=img
-    for edge in v.edges:
-        edge=tools.img_coords(edge,width,height)
-        if edge==None:
-            continue
-
-        cv2.line(img,edge[0],edge[1],255)
-
-    # Show results
-    cv2.imshow("testwin",img)
-    cv2.waitKey(0)
-
-# Test the ability to get the edges surrounding a point. It does the following:
-#  1) Creates a CV window named "testwin"
-#  2) Draw the points into it as tiny circles
-#  3) Draws the edges calculated by voronoi for the specified point
-#  4) Waits for a keystroke
-# The testwin is created based on passed width and height
-def run_test_association(points,plot_points,width,height):
-    # Init open CV
-    cv.NamedWindow('testwin')
-
-    img = numpy.zeros((height,width, 1), numpy.uint8)
-    
-    v=voronoi.Voronoi(points,width,height)
-
-    for point in plot_points:
-        polygon=[]
-    
-        for edge in v.sites[point]:
-            edge=tools.img_coords(edge,width,height)
-            if edge==None:
-                continue
-
-            polygon.append(edge[0])
-            polygon.append(edge[1])
-
-        color=random.randint(150,255)
-        cv2.fillPoly(img,numpy.array([polygon],'int32'),color)
-
-    # Stuff
-    for p in points:
-        cv2.circle(img,toint(p),2,255,-1)
-
-    # Show results
-    cv2.imshow("testwin",img)
-    cv2.waitKey(0)
-

file:a/voronoi.py (deleted)
--- a/voronoi.py
+++ /dev/null
@@ -1,500 +1,1 @@
-#!/usr/bin/env python
 
-import sys
-import random
-from math import *
-from heapq import *
-
-def sign(x):
-    if x>0:
-        return 1
-    elif x<0:
-        return -1
-    else:
-        return 0
-
-class Event:
-    TYPE_SITE=0
-    TYPE_CIRCLE=1
-    def __init__(self,type):
-        self.type=type
-
-class Parabola:
-    def __init__(self,point):
-        self.point=point
-
-    def value(self,x,ypp):
-        if self.point[1]==ypp: # Vertical edge!
-            return None
-
-        return ((self.point[0]-x)**2-ypp**2+self.point[1]**2)/(2*(self.point[1]-ypp))
-
-    def debug(self,ypp):
-        print '((%f-x)**2-%f**2+%f**2)/(2(%f-%f))'%(self.point[0],ypp,self.point[1],self.point[1], ypp)
-
-class Edge:
-    def __init__(self,r,lp,rp):
-        self.r=r
-        self.start=0
-        self.end=None
-        self.lp=lp
-        self.rp=rp
-
-        self.normal=(rp[0]-lp[0],rp[1]-lp[1])
-        A=sqrt(self.normal[0]**2+self.normal[1]**2)
-        self.normal=(self.normal[0]/A,self.normal[1]/A)
-        self.direction=(-self.normal[1],self.normal[0])
-
-    def flip(self):
-        edge=Edge(self.r,self.rp,self.lp)
-        if self.start==None: edge.end=None
-        else: edge.end=-self.start
-        if self.end==None: edge.start=None
-        else: edge.start=-self.end
-
-        return edge
-
-    def value(self,t):
-        return (self.value_x(t),self.value_y(t))
-
-    def value_x(self,t):
-        if t==None:
-            return None
-        return self.direction[0]*t+self.r[0]
-
-    def value_y(self,t):
-        if t==None:
-            return None
-        return self.direction[1]*t+self.r[1]
-
-    def intersects_with(self,edge):
-        x=self.r[0]-edge.r[0]
-        y=self.r[1]-edge.r[1]
-        
-        a=-self.direction[0]
-        b=edge.direction[0]
-        c=-self.direction[1]
-        d=edge.direction[1]
-
-        det=a*d-b*c
-
-        if det==0:
-            return None
-
-        t=(d*x-b*y)/det
-        s=(-c*x+a*y)/det
-
-        return (t,s)
-
-
-class Node:
-    def __init__(self):
-        self.left=None
-        self.right=None
-        self.parent=None
-
-    def clone(self):
-        node=Node()
-        node.parabola=self.parabola
-        return node
-        
-    def set_left(self,left):
-        self.left=left
-        left.parent=self
-
-    def set_right(self,right):
-        self.right=right
-        right.parent=self
-
-    def is_leaf(self):
-        return (self.left==None and self.right==None)
-
-    def find_leaf(self,point):
-        if self.is_leaf():
-            return self
-
-        x=None
-
-        if self.edge.lp[1]==self.edge.rp[1]:
-            b=2*(self.edge.rp[0]-self.edge.lp[0])
-            c=(self.edge.lp[0]**2+self.edge.lp[1]**2-self.edge.rp[0]**2-self.edge.rp[1]**2)
-            x=(self.edge.lp[0]+self.edge.rp[0])/2
-        elif self.edge.rp[1]==point[1]:
-            x=self.edge.rp[0]
-        elif self.edge.lp[1]==point[1]:
-            x=self.edge.lp[0]
-        elif self.edge.lp[1]!=self.edge.rp[1]:
-            r1=1/(self.edge.lp[1]-point[1])
-            r2=1/(self.edge.rp[1]-point[1])
-            a=r1-r2
-            b=2*(self.edge.rp[0]*r2-self.edge.lp[0]*r1)
-            c=r1*(self.edge.lp[0]**2+self.edge.lp[1]**2-point[1]**2)-r2*(self.edge.rp[0]**2+self.edge.rp[1]**2-point[1]**2)
-            
-            x1=(-b-sqrt(b**2-4*a*c))/(2*a)
-            x2=(-b+sqrt(b**2-4*a*c))/(2*a)
-
-            x=x1
-
-        if x==None:
-            print 'WARNING: Unhandled case ',self.edge.rp[1],self.edge.lp[1],point[1]
-        if point[0]<x:
-            return self.left.find_leaf(point)
-        else:
-            return self.right.find_leaf(point)
-
-        return None
-
-    def is_left_node(self):
-        if self.parent==None:
-            return False
-
-        return self.parent.left==self
-
-    def is_right_node(self):
-        if self.parent==None:
-            return False
-
-        return self.parent.right==self
-
-    def first_left_parent(self):
-        if self.parent==None:
-            return None
-
-        if self.parent.left==self:
-            return self.parent
-
-        return self.parent.first_left_parent()
-
-    def first_right_parent(self):
-        if self.parent==None:
-            return None
-
-        if self.parent.right==self:
-            return self.parent
-
-        return self.parent.first_right_parent()
-
-    def left_leaf(self):
-        if self.left.is_leaf():
-            return self.left
-        else:
-            return self.left.left_leaf()
-
-    def right_leaf(self):
-        if self.right.is_leaf():
-            return self.right
-        else:
-            return self.right.right_leaf()
-
-    def prev_parabola(self):
-        lp=self.first_right_parent()
-        if lp==None:
-            return None
-        if lp.left.is_leaf():
-            return lp.left
-        return lp.left.right_leaf()
-
-    def next_parabola(self):
-        rp=self.first_left_parent()
-        if rp==None:
-            return None
-        if rp.right.is_leaf():
-            return rp.right
-        return rp.right.left_leaf()
-
-    def remove(self,reason):
-        if reason==self.left:
-            if self.is_left_node():
-                self.parent.set_left(self.right)
-            else:
-                self.parent.set_right(self.right)
-        else:
-            if self.is_left_node():
-                self.parent.set_left(self.left)
-            else:
-                self.parent.set_right(self.left)
-
-    def dump_tree(self,path):
-        fh=open(path,'w')
-
-        fh.write('digraph G {\n')
-        name='root'
-        self.dump_node(fh,name)
-        fh.write('}\n')
-        fh.close()
-
-    def dump_node(self,fh,name):
-        if self.is_leaf():
-            label='%s\\n(%.1f,%.1f)'%(self,self.parabola.point[0],self.parabola.point[1])
-            fh.write('\t%s [label="%s"];\n'%(name,label))
-        else:
-            label='%s\\np (%.1f,%.1f)\\nd (%.1f,%.1f)'%(self,self.edge.r[0],self.edge.r[1],self.edge.direction[0],self.edge.direction[1])
-            fh.write('\t%s [label="%s"];\n'%(name,label))
-
-            self.left.dump_node(fh,name+'l')
-            self.right.dump_node(fh,name+'r')
-
-            fh.write('\t%s -> %s;\n'%(name,name+'l'));
-            fh.write('\t%s -> %s;\n'%(name,name+'r'));
-            
-    def dump_beachline(self,ypp):
-        if self.is_leaf():
-            return str(self.parabola.point)
-        else:
-            str1=self.left.dump_beachline(ypp)
-            str2=self.right.dump_beachline(ypp)
-
-            r1=1/(self.edge.lp[1]-ypp)
-            r2=1/(self.edge.rp[1]-ypp)
-            a=r1-r2
-            b=2*(self.edge.rp[0]*r2-self.edge.lp[0]*r1)
-            c=r1*(self.edge.lp[0]**2+self.edge.lp[1]**2-ypp**2)-r2*(self.edge.rp[0]**2+self.edge.rp[1]**2-ypp**2)
-            
-            x1=(-b-sqrt(b**2-4*a*c))/(2*a)
-            x2=(-b+sqrt(b**2-4*a*c))/(2*a)
-            
-            x=x1
-
-            return '%s |%f| %s'%(str1,x,str2)
-
-class Voronoi:
-    def __init__(self,points,width,height):
-        self.width=width
-        self.height=height
-
-        # Results are stored here
-        self.edges=[]
-        self.sites={}
-
-        # Sort the points in x first
-        self.points=[]
-        for point in points:
-            insert=len(self.points)
-            for idx in range(0,len(self.points)):
-                if self.points[idx][0]>point[0]:
-                    insert=idx
-                    break
-            self.points.insert(insert,point)
-            self.sites[point]=[]
-
-        # Run
-        self.run()
-
-    def run(self):
-        self.root=None
-        self.internal_edges=[]
-
-        self.events=[]
-        for point in self.points:
-            e=Event(Event.TYPE_SITE)
-            e.point=point
-            heappush(self.events,((point[1],e)))
-
-        done=False
-        while not done:
-            if len(self.events)==0:
-                done=True
-                continue
-            
-            e_info=heappop(self.events)
-            self.sweepline=e_info[0]
-            e=e_info[1]
-#            print 'SWEEPLINE',self.sweepline
-            if e.type==Event.TYPE_SITE:
-                self.add_parabola(e.point)
-            elif e.type==Event.TYPE_CIRCLE:
-                self.remove_parabola(e.parabola_node)
-    
-        self.finish_edges()
-
-    def add_parabola(self,point):
-#        print 'add_parabola',point
-        par=Parabola(point)
- 
-        if self.root==None:
-            node=Node()
-            node.parabola=par
-            self.root=node
-            return
-
-        if self.root.is_leaf():
-            leaf=self.root
-        else:
-            leaf=self.root.find_leaf(point)
-#        print 'splice into',leaf.parabola.point
-        if hasattr(leaf,'circle_event'):
-            self.remove_event(leaf.circle_event)
-            del leaf.circle_event
-
-        self.insert_parabola(leaf,par)
-#        self.root.dump_tree('add_parabola_%s.dot'%str(point))
-
-    def insert_parabola(self,node,parabola):
-        if parabola.point[1]!=node.parabola.point[1]:
-            intersection=(parabola.point[0],node.parabola.value(parabola.point[0],parabola.point[1]))
-
-            newnode=Node()
-            newnode.parabola=parabola
-
-            curleaf1=node.clone()
-            curleaf2=node.clone()
-
-            edge1=Edge(intersection,node.parabola.point,parabola.point)
-            edge2=Edge(intersection,parabola.point,node.parabola.point)
-            self.internal_edges.append((edge1,edge2))
-
-            edge_node2=Node()
-            edge_node2.set_left(newnode)
-            edge_node2.set_right(curleaf2)
-            edge_node2.edge=edge2
-        
-            node.set_left(curleaf1)
-            node.set_right(edge_node2)
-            node.edge=edge1
-
-            # Check circle events
-            if curleaf1.prev_parabola()!=None and curleaf1.next_parabola()!=None:
-                self.circle_check(curleaf1)
-            if curleaf2.prev_parabola()!=None and curleaf2.next_parabola()!=None:
-                self.circle_check(curleaf2)
-        else:
-            intersection=((parabola.point[0]+node.parabola.point[0])/2,
-                          (parabola.point[1]+node.parabola.point[1])/2)
-
-            newnode=Node()
-            newnode.parabola=parabola
-            
-            curleaf=node.clone()
-
-            edge=Edge(intersection,node.parabola.point,parabola.point)
-            edge.start=None
-            self.internal_edges.append(edge)
-
-            node.set_left(curleaf)
-            node.set_right(newnode)
-            node.edge=edge
-
-
-    def circle_check(self,circle_node):
-#        print 'circle_check'
-        if circle_node.prev_parabola()==None or circle_node.next_parabola()==None:
-            return
-        if circle_node.prev_parabola().parabola==circle_node.next_parabola().parabola:
-            return
-#        print 'circles',circle_node.prev_parabola().parabola.point,circle_node.parabola.point,circle_node.next_parabola().parabola.point
-
-        edge_to_end1=circle_node.first_left_parent().edge
-        edge_to_end2=circle_node.first_right_parent().edge
-        t,s=edge_to_end1.intersects_with(edge_to_end2)
-        
-#        print 't,s',t,s
-        if (t<0 and edge_to_end1.start!=None) or (s<0 and edge_to_end2.start!=None):
-            return
-
-        if t==s and t==0:
-            return
-
-        center=edge_to_end1.value(t)
-        R=sqrt((circle_node.parabola.point[0]-center[0])**2+(circle_node.parabola.point[1]-center[1])**2)
-
-        ytrig=R+center[1]
-
-        e=Event(Event.TYPE_CIRCLE)
-        e.center=center
-        e.t=t
-        e.s=s
-        e.parabola_node=circle_node
-        circle_node.circle_event=e
-        
-        heappush(self.events,(ytrig,e))
-
-    def remove_parabola(self,parabola_node):
-#        print 'remove_parabola',parabola_node,parabola_node.parabola.point
-        prev_parabola_node=parabola_node.prev_parabola()
-        next_parabola_node=parabola_node.next_parabola()
-
-        if hasattr(prev_parabola_node,'circle_event'):
-            self.remove_event(prev_parabola_node.circle_event)
-            del prev_parabola_node.circle_event
-
-        if hasattr(next_parabola_node,'circle_event'):
-            self.remove_event(next_parabola_node.circle_event)
-            del next_parabola_node.circle_event
- 
-        edge_to_end1=parabola_node.first_left_parent().edge
-        edge_to_end2=parabola_node.first_right_parent().edge
-        
-        edge_to_end1.end=parabola_node.circle_event.t
-        edge_to_end2.end=parabola_node.circle_event.s
-
-        # Remove this arch
-        parabola_node.parent.remove(parabola_node)
-
-        # Connect the two parabolas that are left
-        rp=prev_parabola_node.first_left_parent()
-        rp.edge=Edge(parabola_node.circle_event.center,prev_parabola_node.parabola.point,next_parabola_node.parabola.point)
-#        print 'fromto',prev_parabola_node.parabola.point,next_parabola_node.parabola.point,rp.edge.direction
-        self.internal_edges.append(rp.edge)
-
-        self.circle_check(prev_parabola_node)
-        self.circle_check(next_parabola_node)
-
-#        self.root.dump_tree('remove_parabola_%s.dot'%str(parabola_node))
-
-    def finish_edges(self):
-        for edge in self.internal_edges:
-            if type(edge)==tuple:
-                tmpedges=edge
-                edge=tmpedges[0]
-                if tmpedges[1].end==None:
-                    edge.start=None
-                else:
-                    edge.start=-tmpedges[1].end
-
-            self.edges.append(edge)
-            self.sites[edge.lp].append(edge)
-            self.sites[edge.rp].append(edge)
-            
-        new_sites={}
-        for site in self.sites:
-            edges=self.sites[site]
-            new_edges=[]
-            # Step 1: Correctly orient edges such that lp is always the site
-            for edge in edges:
-                if edge.lp!=site:
-                    edge=edge.flip()
-                new_edges.append(edge)
-
-            # Step 2: Find the first edge (aka one that starts at inifinity, or 
-            # random if no such edge exists).
-            first_edge=new_edges[0]
-            for edge in new_edges:
-                if edge.start==None:
-                    first_edge=edge
-                    break
-
-			# Step 3: Find subsequent edges
-            self.sites[site]=[first_edge]
-            new_edges.remove(first_edge)
-            last_end=first_edge.value(first_edge.end)
-            while len(new_edges)!=0:
-                for edge in new_edges:
-                    start=edge.value(edge.start)
-                    if (start[0]-last_end[0])**2+(start[1]-last_end[1])**2<sys.float_info.epsilon:
-                        last_end=edge.value(edge.end)
-                        new_edges.remove(edge)
-                        self.sites[site].append(edge)
-                        break
-
-
-
-            new_sites[site]=new_edges
-            
- 
-    def remove_event(self,event):
-        for entry in self.events:
-            if event==entry[1]:
-                self.events.remove(entry)
-                break
-        heapify(self.events)
-

--- /dev/null
+++ b/voronoi/__init__.py

--- /dev/null
+++ b/voronoi/fortune.py
@@ -1,1 +1,500 @@
-
+#!/usr/bin/env python
+
+import sys
+import random
+from math import *
+from heapq import *
+
+def sign(x):
+    if x>0:
+        return 1
+    elif x<0:
+        return -1
+    else:
+        return 0
+
+class Event:
+    TYPE_SITE=0
+    TYPE_CIRCLE=1
+    def __init__(self,type):
+        self.type=type
+
+class Parabola:
+    def __init__(self,point):
+        self.point=point
+
+    def value(self,x,ypp):
+        if self.point[1]==ypp: # Vertical edge!
+            return None
+
+        return ((self.point[0]-x)**2-ypp**2+self.point[1]**2)/(2*(self.point[1]-ypp))
+
+    def debug(self,ypp):
+        print '((%f-x)**2-%f**2+%f**2)/(2(%f-%f))'%(self.point[0],ypp,self.point[1],self.point[1], ypp)
+
+class Edge:
+    def __init__(self,r,lp,rp):
+        self.r=r
+        self.start=0
+        self.end=None
+        self.lp=lp
+        self.rp=rp
+
+        self.normal=(rp[0]-lp[0],rp[1]-lp[1])
+        A=sqrt(self.normal[0]**2+self.normal[1]**2)
+        self.normal=(self.normal[0]/A,self.normal[1]/A)
+        self.direction=(-self.normal[1],self.normal[0])
+
+    def flip(self):
+        edge=Edge(self.r,self.rp,self.lp)
+        if self.start==None: edge.end=None
+        else: edge.end=-self.start
+        if self.end==None: edge.start=None
+        else: edge.start=-self.end
+
+        return edge
+
+    def value(self,t):
+        return (self.value_x(t),self.value_y(t))
+
+    def value_x(self,t):
+        if t==None:
+            return None
+        return self.direction[0]*t+self.r[0]
+
+    def value_y(self,t):
+        if t==None:
+            return None
+        return self.direction[1]*t+self.r[1]
+
+    def intersects_with(self,edge):
+        x=self.r[0]-edge.r[0]
+        y=self.r[1]-edge.r[1]
+        
+        a=-self.direction[0]
+        b=edge.direction[0]
+        c=-self.direction[1]
+        d=edge.direction[1]
+
+        det=a*d-b*c
+
+        if det==0:
+            return None
+
+        t=(d*x-b*y)/det
+        s=(-c*x+a*y)/det
+
+        return (t,s)
+
+
+class Node:
+    def __init__(self):
+        self.left=None
+        self.right=None
+        self.parent=None
+
+    def clone(self):
+        node=Node()
+        node.parabola=self.parabola
+        return node
+        
+    def set_left(self,left):
+        self.left=left
+        left.parent=self
+
+    def set_right(self,right):
+        self.right=right
+        right.parent=self
+
+    def is_leaf(self):
+        return (self.left==None and self.right==None)
+
+    def find_leaf(self,point):
+        if self.is_leaf():
+            return self
+
+        x=None
+
+        if self.edge.lp[1]==self.edge.rp[1]:
+            b=2*(self.edge.rp[0]-self.edge.lp[0])
+            c=(self.edge.lp[0]**2+self.edge.lp[1]**2-self.edge.rp[0]**2-self.edge.rp[1]**2)
+            x=(self.edge.lp[0]+self.edge.rp[0])/2
+        elif self.edge.rp[1]==point[1]:
+            x=self.edge.rp[0]
+        elif self.edge.lp[1]==point[1]:
+            x=self.edge.lp[0]
+        elif self.edge.lp[1]!=self.edge.rp[1]:
+            r1=1/(self.edge.lp[1]-point[1])
+            r2=1/(self.edge.rp[1]-point[1])
+            a=r1-r2
+            b=2*(self.edge.rp[0]*r2-self.edge.lp[0]*r1)
+            c=r1*(self.edge.lp[0]**2+self.edge.lp[1]**2-point[1]**2)-r2*(self.edge.rp[0]**2+self.edge.rp[1]**2-point[1]**2)
+            
+            x1=(-b-sqrt(b**2-4*a*c))/(2*a)
+            x2=(-b+sqrt(b**2-4*a*c))/(2*a)
+
+            x=x1
+
+        if x==None:
+            print 'WARNING: Unhandled case ',self.edge.rp[1],self.edge.lp[1],point[1]
+        if point[0]<x:
+            return self.left.find_leaf(point)
+        else:
+            return self.right.find_leaf(point)
+
+        return None
+
+    def is_left_node(self):
+        if self.parent==None:
+            return False
+
+        return self.parent.left==self
+
+    def is_right_node(self):
+        if self.parent==None:
+            return False
+
+        return self.parent.right==self
+
+    def first_left_parent(self):
+        if self.parent==None:
+            return None
+
+        if self.parent.left==self:
+            return self.parent
+
+        return self.parent.first_left_parent()
+
+    def first_right_parent(self):
+        if self.parent==None:
+            return None
+
+        if self.parent.right==self:
+            return self.parent
+
+        return self.parent.first_right_parent()
+
+    def left_leaf(self):
+        if self.left.is_leaf():
+            return self.left
+        else:
+            return self.left.left_leaf()
+
+    def right_leaf(self):
+        if self.right.is_leaf():
+            return self.right
+        else:
+            return self.right.right_leaf()
+
+    def prev_parabola(self):
+        lp=self.first_right_parent()
+        if lp==None:
+            return None
+        if lp.left.is_leaf():
+            return lp.left
+        return lp.left.right_leaf()
+
+    def next_parabola(self):
+        rp=self.first_left_parent()
+        if rp==None:
+            return None
+        if rp.right.is_leaf():
+            return rp.right
+        return rp.right.left_leaf()
+
+    def remove(self,reason):
+        if reason==self.left:
+            if self.is_left_node():
+                self.parent.set_left(self.right)
+            else:
+                self.parent.set_right(self.right)
+        else:
+            if self.is_left_node():
+                self.parent.set_left(self.left)
+            else:
+                self.parent.set_right(self.left)
+
+    def dump_tree(self,path):
+        fh=open(path,'w')
+
+        fh.write('digraph G {\n')
+        name='root'
+        self.dump_node(fh,name)
+        fh.write('}\n')
+        fh.close()
+
+    def dump_node(self,fh,name):
+        if self.is_leaf():
+            label='%s\\n(%.1f,%.1f)'%(self,self.parabola.point[0],self.parabola.point[1])
+            fh.write('\t%s [label="%s"];\n'%(name,label))
+        else:
+            label='%s\\np (%.1f,%.1f)\\nd (%.1f,%.1f)'%(self,self.edge.r[0],self.edge.r[1],self.edge.direction[0],self.edge.direction[1])
+            fh.write('\t%s [label="%s"];\n'%(name,label))
+
+            self.left.dump_node(fh,name+'l')
+            self.right.dump_node(fh,name+'r')
+
+            fh.write('\t%s -> %s;\n'%(name,name+'l'));
+            fh.write('\t%s -> %s;\n'%(name,name+'r'));
+            
+    def dump_beachline(self,ypp):
+        if self.is_leaf():
+            return str(self.parabola.point)
+        else:
+            str1=self.left.dump_beachline(ypp)
+            str2=self.right.dump_beachline(ypp)
+
+            r1=1/(self.edge.lp[1]-ypp)
+            r2=1/(self.edge.rp[1]-ypp)
+            a=r1-r2
+            b=2*(self.edge.rp[0]*r2-self.edge.lp[0]*r1)
+            c=r1*(self.edge.lp[0]**2+self.edge.lp[1]**2-ypp**2)-r2*(self.edge.rp[0]**2+self.edge.rp[1]**2-ypp**2)
+            
+            x1=(-b-sqrt(b**2-4*a*c))/(2*a)
+            x2=(-b+sqrt(b**2-4*a*c))/(2*a)
+            
+            x=x1
+
+            return '%s |%f| %s'%(str1,x,str2)
+
+class Voronoi:
+    def __init__(self,points,width,height):
+        self.width=width
+        self.height=height
+
+        # Results are stored here
+        self.edges=[]
+        self.sites={}
+
+        # Sort the points in x first
+        self.points=[]
+        for point in points:
+            insert=len(self.points)
+            for idx in range(0,len(self.points)):
+                if self.points[idx][0]>point[0]:
+                    insert=idx
+                    break
+            self.points.insert(insert,point)
+            self.sites[point]=[]
+
+        # Run
+        self.run()
+
+    def run(self):
+        self.root=None
+        self.internal_edges=[]
+
+        self.events=[]
+        for point in self.points:
+            e=Event(Event.TYPE_SITE)
+            e.point=point
+            heappush(self.events,((point[1],e)))
+
+        done=False
+        while not done:
+            if len(self.events)==0:
+                done=True
+                continue
+            
+            e_info=heappop(self.events)
+            self.sweepline=e_info[0]
+            e=e_info[1]
+#            print 'SWEEPLINE',self.sweepline
+            if e.type==Event.TYPE_SITE:
+                self.add_parabola(e.point)
+            elif e.type==Event.TYPE_CIRCLE:
+                self.remove_parabola(e.parabola_node)
+    
+        self.finish_edges()
+
+    def add_parabola(self,point):
+#        print 'add_parabola',point
+        par=Parabola(point)
+ 
+        if self.root==None:
+            node=Node()
+            node.parabola=par
+            self.root=node
+            return
+
+        if self.root.is_leaf():
+            leaf=self.root
+        else:
+            leaf=self.root.find_leaf(point)
+#        print 'splice into',leaf.parabola.point
+        if hasattr(leaf,'circle_event'):
+            self.remove_event(leaf.circle_event)
+            del leaf.circle_event
+
+        self.insert_parabola(leaf,par)
+#        self.root.dump_tree('add_parabola_%s.dot'%str(point))
+
+    def insert_parabola(self,node,parabola):
+        if parabola.point[1]!=node.parabola.point[1]:
+            intersection=(parabola.point[0],node.parabola.value(parabola.point[0],parabola.point[1]))
+
+            newnode=Node()
+            newnode.parabola=parabola
+
+            curleaf1=node.clone()
+            curleaf2=node.clone()
+
+            edge1=Edge(intersection,node.parabola.point,parabola.point)
+            edge2=Edge(intersection,parabola.point,node.parabola.point)
+            self.internal_edges.append((edge1,edge2))
+
+            edge_node2=Node()
+            edge_node2.set_left(newnode)
+            edge_node2.set_right(curleaf2)
+            edge_node2.edge=edge2
+        
+            node.set_left(curleaf1)
+            node.set_right(edge_node2)
+            node.edge=edge1
+
+            # Check circle events
+            if curleaf1.prev_parabola()!=None and curleaf1.next_parabola()!=None:
+                self.circle_check(curleaf1)
+            if curleaf2.prev_parabola()!=None and curleaf2.next_parabola()!=None:
+                self.circle_check(curleaf2)
+        else:
+            intersection=((parabola.point[0]+node.parabola.point[0])/2,
+                          (parabola.point[1]+node.parabola.point[1])/2)
+
+            newnode=Node()
+            newnode.parabola=parabola
+            
+            curleaf=node.clone()
+
+            edge=Edge(intersection,node.parabola.point,parabola.point)
+            edge.start=None
+            self.internal_edges.append(edge)
+
+            node.set_left(curleaf)
+            node.set_right(newnode)
+            node.edge=edge
+
+
+    def circle_check(self,circle_node):
+#        print 'circle_check'
+        if circle_node.prev_parabola()==None or circle_node.next_parabola()==None:
+            return
+        if circle_node.prev_parabola().parabola==circle_node.next_parabola().parabola:
+            return
+#        print 'circles',circle_node.prev_parabola().parabola.point,circle_node.parabola.point,circle_node.next_parabola().parabola.point
+
+        edge_to_end1=circle_node.first_left_parent().edge
+        edge_to_end2=circle_node.first_right_parent().edge
+        t,s=edge_to_end1.intersects_with(edge_to_end2)
+        
+#        print 't,s',t,s
+        if (t<0 and edge_to_end1.start!=None) or (s<0 and edge_to_end2.start!=None):
+            return
+
+        if t==s and t==0:
+            return
+
+        center=edge_to_end1.value(t)
+        R=sqrt((circle_node.parabola.point[0]-center[0])**2+(circle_node.parabola.point[1]-center[1])**2)
+
+        ytrig=R+center[1]
+
+        e=Event(Event.TYPE_CIRCLE)
+        e.center=center
+        e.t=t
+        e.s=s
+        e.parabola_node=circle_node
+        circle_node.circle_event=e
+        
+        heappush(self.events,(ytrig,e))
+
+    def remove_parabola(self,parabola_node):
+#        print 'remove_parabola',parabola_node,parabola_node.parabola.point
+        prev_parabola_node=parabola_node.prev_parabola()
+        next_parabola_node=parabola_node.next_parabola()
+
+        if hasattr(prev_parabola_node,'circle_event'):
+            self.remove_event(prev_parabola_node.circle_event)
+            del prev_parabola_node.circle_event
+
+        if hasattr(next_parabola_node,'circle_event'):
+            self.remove_event(next_parabola_node.circle_event)
+            del next_parabola_node.circle_event
+ 
+        edge_to_end1=parabola_node.first_left_parent().edge
+        edge_to_end2=parabola_node.first_right_parent().edge
+        
+        edge_to_end1.end=parabola_node.circle_event.t
+        edge_to_end2.end=parabola_node.circle_event.s
+
+        # Remove this arch
+        parabola_node.parent.remove(parabola_node)
+
+        # Connect the two parabolas that are left
+        rp=prev_parabola_node.first_left_parent()
+        rp.edge=Edge(parabola_node.circle_event.center,prev_parabola_node.parabola.point,next_parabola_node.parabola.point)
+#        print 'fromto',prev_parabola_node.parabola.point,next_parabola_node.parabola.point,rp.edge.direction
+        self.internal_edges.append(rp.edge)
+
+        self.circle_check(prev_parabola_node)
+        self.circle_check(next_parabola_node)
+
+#        self.root.dump_tree('remove_parabola_%s.dot'%str(parabola_node))
+
+    def finish_edges(self):
+        for edge in self.internal_edges:
+            if type(edge)==tuple:
+                tmpedges=edge
+                edge=tmpedges[0]
+                if tmpedges[1].end==None:
+                    edge.start=None
+                else:
+                    edge.start=-tmpedges[1].end
+
+            self.edges.append(edge)
+            self.sites[edge.lp].append(edge)
+            self.sites[edge.rp].append(edge)
+            
+        new_sites={}
+        for site in self.sites:
+            edges=self.sites[site]
+            new_edges=[]
+            # Step 1: Correctly orient edges such that lp is always the site
+            for edge in edges:
+                if edge.lp!=site:
+                    edge=edge.flip()
+                new_edges.append(edge)
+
+            # Step 2: Find the first edge (aka one that starts at inifinity, or 
+            # random if no such edge exists).
+            first_edge=new_edges[0]
+            for edge in new_edges:
+                if edge.start==None:
+                    first_edge=edge
+                    break
+
+			# Step 3: Find subsequent edges
+            self.sites[site]=[first_edge]
+            new_edges.remove(first_edge)
+            last_end=first_edge.value(first_edge.end)
+            while len(new_edges)!=0:
+                for edge in new_edges:
+                    start=edge.value(edge.start)
+                    if (start[0]-last_end[0])**2+(start[1]-last_end[1])**2<sys.float_info.epsilon:
+                        last_end=edge.value(edge.end)
+                        new_edges.remove(edge)
+                        self.sites[site].append(edge)
+                        break
+
+
+
+            new_sites[site]=new_edges
+            
+ 
+    def remove_event(self,event):
+        for entry in self.events:
+            if event==entry[1]:
+                self.events.remove(entry)
+                break
+        heapify(self.events)
+

file:b/voronoi/tools.py (new)
--- /dev/null
+++ b/voronoi/tools.py
@@ -1,1 +1,57 @@
+def img_coords(edge,width,height):
+    intersections=[]
 
+    # Intersection with north and south edges (horizontal)
+    if edge.direction[1]!=0:
+        # North (y=height)
+        t=(height-edge.r[1])/edge.direction[1]
+        x=edge.value_x(t)
+        if 0<x and x<width:
+            intersections.append(t)
+        # South (y=0)
+        t=-edge.r[1]/edge.direction[1]
+        x=edge.value_x(t)
+        if 0<x and x<width:
+            intersections.append(t)
+
+    # Intersections with east and west (vertical)
+    if edge.direction[0]!=0:
+        # East (x=height)
+        t=(width-edge.r[0])/edge.direction[0]
+        y=edge.value_y(t)
+        if 0<y and y<height:
+            intersections.append(t)
+        # West (x=0)
+        t=-edge.r[0]/edge.direction[0]
+        y=edge.value_y(t)
+        if 0<y and y<height:
+            intersections.append(t)
+
+    if len(intersections)==0:
+        return None # Line misses image completely
+
+    # Figure out the range of values the line takes inside the box
+    mint=min(intersections)
+    maxt=max(intersections)
+
+    # Determine what range of values the line defined as now takes
+    if edge.start==None: # Make sure to clip the edge if it is infinite
+        edge.start=mint # This is when the box is left
+
+    if edge.end==None: # Make sure to clip the edge if it is infinite
+        edge.end=maxt # This is when the box is left
+
+    p1t=edge.start
+    p2t=edge.end
+
+    if p2t<mint or p1t>maxt: # This edge lies outside of the box
+        return None
+
+    # Determine the points that the edge takes after being clipped by the box
+    p1t=max(p1t,mint)
+    p2t=min(p2t,maxt)
+
+    start=edge.value(p1t)
+    end=edge.value(p2t)
+    return ((int(start[0]),int(start[1])),(int(end[0]),int(end[1])))
+

--- /dev/null
+++ b/voronoi/unittests.py
@@ -1,1 +1,74 @@
+import numpy
+import fortune
+import cv,cv2
+import tools
+import random
 
+# Useful functions
+def toint(t2):
+    return (int(t2[0]),int(t2[1]))
+
+# Run a test on a set of points. The test does the following:
+#  1) Creates a CV window named "testwin"
+#  2) Draw the points into it as tiny circles
+#  3) Draws the edges calculated by voronoi
+#  4) Waits for a keystroke
+# The testwin is created based on passed width and height
+def run_test(points,width,height):
+    # Init open CV
+    cv.NamedWindow('testwin')
+
+    img = numpy.zeros((height,width, 1), numpy.uint8)
+    
+    for p in points:
+        cv2.circle(img,toint(p),2,255,-1)
+
+    v=fortune.Voronoi(points,width,height)
+    v.img=img
+    for edge in v.edges:
+        edge=tools.img_coords(edge,width,height)
+        if edge==None:
+            continue
+
+        cv2.line(img,edge[0],edge[1],255)
+
+    # Show results
+    cv2.imshow("testwin",img)
+    cv2.waitKey(0)
+
+# Test the ability to get the edges surrounding a point. It does the following:
+#  1) Creates a CV window named "testwin"
+#  2) Draw the points into it as tiny circles
+#  3) Draws the edges calculated by voronoi for the specified point
+#  4) Waits for a keystroke
+# The testwin is created based on passed width and height
+def run_test_association(points,plot_points,width,height):
+    # Init open CV
+    cv.NamedWindow('testwin')
+
+    img = numpy.zeros((height,width, 1), numpy.uint8)
+    
+    v=fortune.Voronoi(points,width,height)
+
+    for point in plot_points:
+        polygon=[]
+    
+        for edge in v.sites[point]:
+            edge=tools.img_coords(edge,width,height)
+            if edge==None:
+                continue
+
+            polygon.append(edge[0])
+            polygon.append(edge[1])
+
+        color=random.randint(150,255)
+        cv2.fillPoly(img,numpy.array([polygon],'int32'),color)
+
+    # Stuff
+    for p in points:
+        cv2.circle(img,toint(p),2,255,-1)
+
+    # Show results
+    cv2.imshow("testwin",img)
+    cv2.waitKey(0)
+