1 #############################################################################
2 #
3 # Voronoi diagram calculator/ Delaunay triangulator
4 # Translated to Python by Bill Simons
5 # September, 2005
6 #
7 # Calculate Delaunay triangulation or the Voronoi polygons for a set of
8 # 2D input points.
9 #
10 # Derived from code bearing the following notice:
11 #
12 # The author of this software is Steven Fortune. Copyright (c) 1994 by AT&T
13 # Bell Laboratories.
14 # Permission to use, copy, modify, and distribute this software for any
15 # purpose without fee is hereby granted, provided that this entire notice
16 # is included in all copies of any software which is or includes a copy
17 # or modification of this software and in all copies of the supporting
18 # documentation for such software.
19 # THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
20 # WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
21 # REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
22 # OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
23 #
24 # Comments were incorporated from Shane O'Sullivan's translation of the
25 # original code into C++ (http://mapviewer.skynet.ie/voronoi.html)
26 #
27 # Steve Fortune's homepage: http://netlib.bell-labs.com/cm/cs/who/sjf/index.html
28 #
29 #############################################################################
31 def usage():
32 print """
33 voronoi - compute Voronoi diagram or Delaunay triangulation
35 voronoi [-t -p -d] [filename]
37 Voronoi reads from filename (or standard input if no filename given) for a set
38 of points in the plane and writes either the Voronoi diagram or the Delaunay
39 triangulation to the standard output. Each input line should consist of two
40 real numbers, separated by white space.
42 If option -t is present, the Delaunay triangulation is produced.
43 Each output line is a triple i j k, which are the indices of the three points
44 in a Delaunay triangle. Points are numbered starting at 0.
46 If option -t is not present, the Voronoi diagram is produced.
47 There are four output record types.
49 s a b indicates that an input point at coordinates a b was seen.
50 l a b c indicates a line with equation ax + by = c.
51 v a b indicates a vertex at a b.
52 e l v1 v2 indicates a Voronoi segment which is a subsegment of line number l
53 with endpoints numbered v1 and v2. If v1 or v2 is -1, the line
54 extends to infinity.
56 Other options include:
58 d Print debugging info
60 p Produce output suitable for input to plot (1), rather than the forms
61 described above.
63 On unsorted data uniformly distributed in the unit square, voronoi uses about
64 20n+140 bytes of storage.
66 AUTHOR
67 Steve J. Fortune (1987) A Sweepline Algorithm for Voronoi Diagrams,
68 Algorithmica 2, 153-174.
69 """
71 #############################################################################
72 #
73 # For programmatic use two functions are available:
74 #
75 # computeVoronoiDiagram(points)
76 #
77 # Takes a list of point objects (which must have x and y fields).
78 # Returns a 3-tuple of:
79 #
80 # (1) a list of 2-tuples, which are the x,y coordinates of the
81 # Voronoi diagram vertices
82 # (2) a list of 3-tuples (a,b,c) which are the equations of the
83 # lines in the Voronoi diagram: a*x + b*y = c
84 # (3) a list of 3-tuples, (l, v1, v2) representing edges of the
85 # Voronoi diagram. l is the index of the line, v1 and v2 are
86 # the indices of the vetices at the end of the edge. If
87 # v1 or v2 is -1, the line extends to infinity.
88 #
89 # computeDelaunayTriangulation(points):
90 #
91 # Takes a list of point objects (which must have x and y fields).
92 # Returns a list of 3-tuples: the indices of the points that form a
93 # Delaunay triangle.
94 #
95 #############################################################################
96 import math
97 import sys
98 import getopt
99 TOLERANCE = 1e-9
100 BIG_FLOAT = 1e38
102 #------------------------------------------------------------------
103 class Context(object):
104 def __init__(self):
105 self.doPrint = 0
106 self.debug = 0
107 self.plot = 0
108 self.triangulate = False
109 self.vertices = [] # list of vertex 2-tuples: (x,y)
110 self.lines = [] # equation of line 3-tuple (a b c), for the equation of the line a*x+b*y = c
111 self.edges = [] # edge 3-tuple: (line index, vertex 1 index, vertex 2 index) if either vertex index is -1, the edge extends to infiinity
112 self.triangles = [] # 3-tuple of vertex indices
114 def circle(self,x,y,rad):
115 pass
117 def clip_line(self,edge):
118 pass
120 def line(self,x0,y0,x1,y1):
121 pass
123 def outSite(self,s):
124 if(self.debug):
125 print "site (%d) at %f %f" % (s.sitenum, s.x, s.y)
126 elif(self.triangulate):
127 pass
128 elif(self.plot):
129 self.circle (s.x, s.y, cradius)
130 elif(self.doPrint):
131 print "s %f %f" % (s.x, s.y)
133 def outVertex(self,s):
134 self.vertices.append((s.x,s.y))
135 if(self.debug):
136 print "vertex(%d) at %f %f" % (s.sitenum, s.x, s.y)
137 elif(self.triangulate):
138 pass
139 elif(self.doPrint and not self.plot):
140 print "v %f %f" % (s.x,s.y)
142 def outTriple(self,s1,s2,s3):
143 self.triangles.append((s1.sitenum, s2.sitenum, s3.sitenum))
144 if(self.debug):
145 print "circle through left=%d right=%d bottom=%d" % (s1.sitenum, s2.sitenum, s3.sitenum)
146 elif(self.triangulate and self.doPrint and not self.plot):
147 print "%d %d %d" % (s1.sitenum, s2.sitenum, s3.sitenum)
149 def outBisector(self,edge):
150 self.lines.append((edge.a, edge.b, edge.c))
151 if(self.debug):
152 print "line(%d) %gx+%gy=%g, bisecting %d %d" % (edge.edgenum, edge.a, edge.b, edge.c, edge.reg[0].sitenum, edge.reg[1].sitenum)
153 elif(self.triangulate):
154 if(self.plot):
155 self.line(edge.reg[0].x, edge.reg[0].y, edge.reg[1].x, edge.reg[1].y)
156 elif(self.doPrint and not self.plot):
157 print "l %f %f %f" % (edge.a, edge.b, edge.c)
159 def outEdge(self,edge):
160 sitenumL = -1
161 if edge.ep[Edge.LE] is not None:
162 sitenumL = edge.ep[Edge.LE].sitenum
163 sitenumR = -1
164 if edge.ep[Edge.RE] is not None:
165 sitenumR = edge.ep[Edge.RE].sitenum
166 self.edges.append((edge.edgenum,sitenumL,sitenumR))
167 if(not self.triangulate):
168 if self.plot:
169 self.clip_line(edge)
170 elif(self.doPrint):
171 print "e %d" % edge.edgenum,
172 print " %d " % sitenumL,
173 print "%d" % sitenumR
175 #------------------------------------------------------------------
176 def voronoi(siteList,context):
177 edgeList = EdgeList(siteList.xmin,siteList.xmax,len(siteList))
178 priorityQ = PriorityQueue(siteList.ymin,siteList.ymax,len(siteList))
179 siteIter = siteList.iterator()
181 bottomsite = siteIter.next()
182 context.outSite(bottomsite)
183 newsite = siteIter.next()
184 minpt = Site(-BIG_FLOAT,-BIG_FLOAT)
185 while True:
186 if not priorityQ.isEmpty():
187 minpt = priorityQ.getMinPt()
189 if (newsite and (priorityQ.isEmpty() or cmp(newsite,minpt) < 0)):
190 # newsite is smallest - this is a site event
191 context.outSite(newsite)
193 # get first Halfedge to the LEFT and RIGHT of the new site
194 lbnd = edgeList.leftbnd(newsite)
195 rbnd = lbnd.right
197 # if this halfedge has no edge, bot = bottom site (whatever that is)
198 # create a new edge that bisects
199 bot = lbnd.rightreg(bottomsite)
200 edge = Edge.bisect(bot,newsite)
201 context.outBisector(edge)
203 # create a new Halfedge, setting its pm field to 0 and insert
204 # this new bisector edge between the left and right vectors in
205 # a linked list
206 bisector = Halfedge(edge,Edge.LE)
207 edgeList.insert(lbnd,bisector)
209 # if the new bisector intersects with the left edge, remove
210 # the left edge's vertex, and put in the new one
211 p = lbnd.intersect(bisector)
212 if p is not None:
213 priorityQ.delete(lbnd)
214 priorityQ.insert(lbnd,p,newsite.distance(p))
216 # create a new Halfedge, setting its pm field to 1
217 # insert the new Halfedge to the right of the original bisector
218 lbnd = bisector
219 bisector = Halfedge(edge,Edge.RE)
220 edgeList.insert(lbnd,bisector)
222 # if this new bisector intersects with the right Halfedge
223 p = bisector.intersect(rbnd)
224 if p is not None:
225 # push the Halfedge into the ordered linked list of vertices
226 priorityQ.insert(bisector,p,newsite.distance(p))
228 newsite = siteIter.next()
230 elif not priorityQ.isEmpty():
231 # intersection is smallest - this is a vector (circle) event
233 # pop the Halfedge with the lowest vector off the ordered list of
234 # vectors. Get the Halfedge to the left and right of the above HE
235 # and also the Halfedge to the right of the right HE
236 lbnd = priorityQ.popMinHalfedge()
237 llbnd = lbnd.left
238 rbnd = lbnd.right
239 rrbnd = rbnd.right
241 # get the Site to the left of the left HE and to the right of
242 # the right HE which it bisects
243 bot = lbnd.leftreg(bottomsite)
244 top = rbnd.rightreg(bottomsite)
246 # output the triple of sites, stating that a circle goes through them
247 mid = lbnd.rightreg(bottomsite)
248 context.outTriple(bot,top,mid)
250 # get the vertex that caused this event and set the vertex number
251 # couldn't do this earlier since we didn't know when it would be processed
252 v = lbnd.vertex
253 siteList.setSiteNumber(v)
254 context.outVertex(v)
256 # set the endpoint of the left and right Halfedge to be this vector
257 if lbnd.edge.setEndpoint(lbnd.pm,v):
258 context.outEdge(lbnd.edge)
260 if rbnd.edge.setEndpoint(rbnd.pm,v):
261 context.outEdge(rbnd.edge)
264 # delete the lowest HE, remove all vertex events to do with the
265 # right HE and delete the right HE
266 edgeList.delete(lbnd)
267 priorityQ.delete(rbnd)
268 edgeList.delete(rbnd)
271 # if the site to the left of the event is higher than the Site
272 # to the right of it, then swap them and set 'pm' to RIGHT
273 pm = Edge.LE
274 if bot.y > top.y:
275 bot,top = top,bot
276 pm = Edge.RE
278 # Create an Edge (or line) that is between the two Sites. This
279 # creates the formula of the line, and assigns a line number to it
280 edge = Edge.bisect(bot, top)
281 context.outBisector(edge)
283 # create a HE from the edge
284 bisector = Halfedge(edge, pm)
286 # insert the new bisector to the right of the left HE
287 # set one endpoint to the new edge to be the vector point 'v'
288 # If the site to the left of this bisector is higher than the right
289 # Site, then this endpoint is put in position 0; otherwise in pos 1
290 edgeList.insert(llbnd, bisector)
291 if edge.setEndpoint(Edge.RE - pm, v):
292 context.outEdge(edge)
294 # if left HE and the new bisector don't intersect, then delete
295 # the left HE, and reinsert it
296 p = llbnd.intersect(bisector)
297 if p is not None:
298 priorityQ.delete(llbnd);
299 priorityQ.insert(llbnd, p, bot.distance(p))
301 # if right HE and the new bisector don't intersect, then reinsert it
302 p = bisector.intersect(rrbnd)
303 if p is not None:
304 priorityQ.insert(bisector, p, bot.distance(p))
305 else:
306 break
308 he = edgeList.leftend.right
309 while he is not edgeList.rightend:
310 context.outEdge(he.edge)
311 he = he.right
313 #------------------------------------------------------------------
314 def isEqual(a,b,relativeError=TOLERANCE):
315 # is nearly equal to within the allowed relative error
316 norm = max(abs(a),abs(b))
317 return (norm < relativeError) or (abs(a - b) < (relativeError * norm))
319 #------------------------------------------------------------------
320 class Site(object):
321 def __init__(self,x=0.0,y=0.0,sitenum=0):
322 self.x = x
323 self.y = y
324 self.sitenum = sitenum
326 def dump(self):
327 print "Site #%d (%g, %g)" % (self.sitenum,self.x,self.y)
329 def __cmp__(self,other):
330 if self.y < other.y:
331 return -1
332 elif self.y > other.y:
333 return 1
334 elif self.x < other.x:
335 return -1
336 elif self.x > other.x:
337 return 1
338 else:
339 return 0
341 def distance(self,other):
342 dx = self.x - other.x
343 dy = self.y - other.y
344 return math.sqrt(dx*dx + dy*dy)
346 #------------------------------------------------------------------
347 class Edge(object):
348 LE = 0
349 RE = 1
350 EDGE_NUM = 0
351 DELETED = {} # marker value
353 def __init__(self):
354 self.a = 0.0
355 self.b = 0.0
356 self.c = 0.0
357 self.ep = [None,None]
358 self.reg = [None,None]
359 self.edgenum = 0
361 def dump(self):
362 print "(#%d a=%g, b=%g, c=%g)" % (self.edgenum,self.a,self.b,self.c)
363 print "ep",self.ep
364 print "reg",self.reg
366 def setEndpoint(self, lrFlag, site):
367 self.ep[lrFlag] = site
368 if self.ep[Edge.RE - lrFlag] is None:
369 return False
370 return True
372 @staticmethod
373 def bisect(s1,s2):
374 newedge = Edge()
375 newedge.reg[0] = s1 # store the sites that this edge is bisecting
376 newedge.reg[1] = s2
378 # to begin with, there are no endpoints on the bisector - it goes to infinity
379 # ep[0] and ep[1] are None
381 # get the difference in x dist between the sites
382 dx = float(s2.x - s1.x)
383 dy = float(s2.y - s1.y)
384 adx = abs(dx) # make sure that the difference in positive
385 ady = abs(dy)
387 # get the slope of the line
388 newedge.c = float(s1.x * dx + s1.y * dy + (dx*dx + dy*dy)*0.5)
389 if adx > ady :
390 # set formula of line, with x fixed to 1
391 newedge.a = 1.0
392 newedge.b = dy/dx
393 newedge.c /= dx
394 else:
395 # set formula of line, with y fixed to 1
396 newedge.b = 1.0
397 newedge.a = dx/dy
398 newedge.c /= dy
400 newedge.edgenum = Edge.EDGE_NUM
401 Edge.EDGE_NUM += 1
402 return newedge
405 #------------------------------------------------------------------
406 class Halfedge(object):
407 def __init__(self,edge=None,pm=Edge.LE):
408 self.left = None # left Halfedge in the edge list
409 self.right = None # right Halfedge in the edge list
410 self.qnext = None # priority queue linked list pointer
411 self.edge = edge # edge list Edge
412 self.pm = pm
413 self.vertex = None # Site()
414 self.ystar = BIG_FLOAT
416 def dump(self):
417 print "Halfedge--------------------------"
418 print "left: ", self.left
419 print "right: ", self.right
420 print "edge: ", self.edge
421 print "pm: ", self.pm
422 print "vertex: ",
423 if self.vertex: self.vertex.dump()
424 else: print "None"
425 print "ystar: ", self.ystar
428 def __cmp__(self,other):
429 if self.ystar > other.ystar:
430 return 1
431 elif self.ystar < other.ystar:
432 return -1
433 elif self.vertex.x > other.vertex.x:
434 return 1
435 elif self.vertex.x < other.vertex.x:
436 return -1
437 else:
438 return 0
440 def leftreg(self,default):
441 if not self.edge:
442 return default
443 elif self.pm == Edge.LE:
444 return self.edge.reg[Edge.LE]
445 else:
446 return self.edge.reg[Edge.RE]
448 def rightreg(self,default):
449 if not self.edge:
450 return default
451 elif self.pm == Edge.LE:
452 return self.edge.reg[Edge.RE]
453 else:
454 return self.edge.reg[Edge.LE]
457 # returns True if p is to right of halfedge self
458 def isPointRightOf(self,pt):
459 e = self.edge
460 topsite = e.reg[1]
461 right_of_site = pt.x > topsite.x
463 if(right_of_site and self.pm == Edge.LE):
464 return True
466 if(not right_of_site and self.pm == Edge.RE):
467 return False
469 if(e.a == 1.0):
470 dyp = pt.y - topsite.y
471 dxp = pt.x - topsite.x
472 fast = 0;
473 if ((not right_of_site and e.b < 0.0) or (right_of_site and e.b >= 0.0)):
474 above = dyp >= e.b * dxp
475 fast = above
476 else:
477 above = pt.x + pt.y * e.b > e.c
478 if(e.b < 0.0):
479 above = not above
480 if (not above):
481 fast = 1
482 if (not fast):
483 dxs = topsite.x - (e.reg[0]).x
484 above = e.b * (dxp*dxp - dyp*dyp) < dxs*dyp*(1.0+2.0*dxp/dxs + e.b*e.b)
485 if(e.b < 0.0):
486 above = not above
487 else: # e.b == 1.0
488 yl = e.c - e.a * pt.x
489 t1 = pt.y - yl
490 t2 = pt.x - topsite.x
491 t3 = yl - topsite.y
492 above = t1*t1 > t2*t2 + t3*t3
494 if(self.pm==Edge.LE):
495 return above
496 else:
497 return not above
499 #--------------------------
500 # create a new site where the Halfedges el1 and el2 intersect
501 def intersect(self,other):
502 e1 = self.edge
503 e2 = other.edge
504 if (e1 is None) or (e2 is None):
505 return None
507 # if the two edges bisect the same parent return None
508 if e1.reg[1] is e2.reg[1]:
509 return None
511 d = e1.a * e2.b - e1.b * e2.a
512 if isEqual(d,0.0):
513 return None
515 xint = (e1.c*e2.b - e2.c*e1.b) / d
516 yint = (e2.c*e1.a - e1.c*e2.a) / d
517 if(cmp(e1.reg[1],e2.reg[1]) < 0):
518 he = self
519 e = e1
520 else:
521 he = other
522 e = e2
524 rightOfSite = xint >= e.reg[1].x
525 if((rightOfSite and he.pm == Edge.LE) or
526 (not rightOfSite and he.pm == Edge.RE)):
527 return None
529 # create a new site at the point of intersection - this is a new
530 # vector event waiting to happen
531 return Site(xint,yint)
535 #------------------------------------------------------------------
536 class EdgeList(object):
537 def __init__(self,xmin,xmax,nsites):
538 if xmin > xmax: xmin,xmax = xmax,xmin
539 self.hashsize = int(2*math.sqrt(nsites+4))
541 self.xmin = xmin
542 self.deltax = float(xmax - xmin)
543 self.hash = [None]*self.hashsize
545 self.leftend = Halfedge()
546 self.rightend = Halfedge()
547 self.leftend.right = self.rightend
548 self.rightend.left = self.leftend
549 self.hash[0] = self.leftend
550 self.hash[-1] = self.rightend
552 def insert(self,left,he):
553 he.left = left
554 he.right = left.right
555 left.right.left = he
556 left.right = he
558 def delete(self,he):
559 he.left.right = he.right
560 he.right.left = he.left
561 he.edge = Edge.DELETED
563 # Get entry from hash table, pruning any deleted nodes
564 def gethash(self,b):
565 if(b < 0 or b >= self.hashsize):
566 return None
567 he = self.hash[b]
568 if he is None or he.edge is not Edge.DELETED:
569 return he
571 # Hash table points to deleted half edge. Patch as necessary.
572 self.hash[b] = None
573 return None
575 def leftbnd(self,pt):
576 # Use hash table to get close to desired halfedge
577 bucket = int(((pt.x - self.xmin)/self.deltax * self.hashsize))
579 if(bucket < 0):
580 bucket =0;
582 if(bucket >=self.hashsize):
583 bucket = self.hashsize-1
585 he = self.gethash(bucket)
586 if(he is None):
587 i = 1
588 while True:
589 he = self.gethash(bucket-i)
590 if (he is not None): break;
591 he = self.gethash(bucket+i)
592 if (he is not None): break;
593 i += 1
595 # Now search linear list of halfedges for the corect one
596 if (he is self.leftend) or (he is not self.rightend and he.isPointRightOf(pt)):
597 he = he.right
598 while he is not self.rightend and he.isPointRightOf(pt):
599 he = he.right
600 he = he.left;
601 else:
602 he = he.left
603 while (he is not self.leftend and not he.isPointRightOf(pt)):
604 he = he.left
606 # Update hash table and reference counts
607 if(bucket > 0 and bucket < self.hashsize-1):
608 self.hash[bucket] = he
609 return he
612 #------------------------------------------------------------------
613 class PriorityQueue(object):
614 def __init__(self,ymin,ymax,nsites):
615 self.ymin = ymin
616 self.deltay = ymax - ymin
617 self.hashsize = int(4 * math.sqrt(nsites))
618 self.count = 0
619 self.minidx = 0
620 self.hash = []
621 for i in range(self.hashsize):
622 self.hash.append(Halfedge())
624 def __len__(self):
625 return self.count
627 def isEmpty(self):
628 return self.count == 0
630 def insert(self,he,site,offset):
631 he.vertex = site
632 he.ystar = site.y + offset
633 last = self.hash[self.getBucket(he)]
634 next = last.qnext
635 while((next is not None) and cmp(he,next) > 0):
636 last = next
637 next = last.qnext
638 he.qnext = last.qnext
639 last.qnext = he
640 self.count += 1
642 def delete(self,he):
643 if (he.vertex is not None):
644 last = self.hash[self.getBucket(he)]
645 while last.qnext is not he:
646 last = last.qnext
647 last.qnext = he.qnext
648 self.count -= 1
649 he.vertex = None
651 def getBucket(self,he):
652 bucket = int(((he.ystar - self.ymin) / self.deltay) * self.hashsize)
653 if bucket < 0: bucket = 0
654 if bucket >= self.hashsize: bucket = self.hashsize-1
655 if bucket < self.minidx: self.minidx = bucket
656 return bucket
658 def getMinPt(self):
659 while(self.hash[self.minidx].qnext is None):
660 self.minidx += 1
661 he = self.hash[self.minidx].qnext
662 x = he.vertex.x
663 y = he.ystar
664 return Site(x,y)
666 def popMinHalfedge(self):
667 curr = self.hash[self.minidx].qnext
668 self.hash[self.minidx].qnext = curr.qnext
669 self.count -= 1
670 return curr
673 #------------------------------------------------------------------
674 class SiteList(object):
675 def __init__(self,pointList):
676 self.__sites = []
677 self.__sitenum = 0
679 self.__xmin = pointList[0].x
680 self.__ymin = pointList[0].y
681 self.__xmax = pointList[0].x
682 self.__ymax = pointList[0].y
683 for i,pt in enumerate(pointList):
684 self.__sites.append(Site(pt.x,pt.y,i))
685 if pt.x < self.__xmin: self.__xmin = pt.x
686 if pt.y < self.__ymin: self.__ymin = pt.y
687 if pt.x > self.__xmax: self.__xmax = pt.x
688 if pt.y > self.__ymax: self.__ymax = pt.y
689 self.__sites.sort()
691 def setSiteNumber(self,site):
692 site.sitenum = self.__sitenum
693 self.__sitenum += 1
695 class Iterator(object):
696 def __init__(this,lst): this.generator = (s for s in lst)
697 def __iter__(this): return this
698 def next(this):
699 try:
700 return this.generator.next()
701 except StopIteration:
702 return None
704 def iterator(self):
705 return SiteList.Iterator(self.__sites)
707 def __iter__(self):
708 return SiteList.Iterator(self.__sites)
710 def __len__(self):
711 return len(self.__sites)
713 def _getxmin(self): return self.__xmin
714 def _getymin(self): return self.__ymin
715 def _getxmax(self): return self.__xmax
716 def _getymax(self): return self.__ymax
717 xmin = property(_getxmin)
718 ymin = property(_getymin)
719 xmax = property(_getxmax)
720 ymax = property(_getymax)
723 #------------------------------------------------------------------
724 def computeVoronoiDiagram(points):
725 """ Takes a list of point objects (which must have x and y fields).
726 Returns a 3-tuple of:
728 (1) a list of 2-tuples, which are the x,y coordinates of the
729 Voronoi diagram vertices
730 (2) a list of 3-tuples (a,b,c) which are the equations of the
731 lines in the Voronoi diagram: a*x + b*y = c
732 (3) a list of 3-tuples, (l, v1, v2) representing edges of the
733 Voronoi diagram. l is the index of the line, v1 and v2 are
734 the indices of the vetices at the end of the edge. If
735 v1 or v2 is -1, the line extends to infinity.
736 """
737 siteList = SiteList(points)
738 context = Context()
739 voronoi(siteList,context)
740 return (context.vertices,context.lines,context.edges)
742 #------------------------------------------------------------------
743 def computeDelaunayTriangulation(points):
744 """ Takes a list of point objects (which must have x and y fields).
745 Returns a list of 3-tuples: the indices of the points that form a
746 Delaunay triangle.
747 """
748 siteList = SiteList(points)
749 context = Context()
750 context.triangulate = true
751 voronoi(siteList,context)
752 return context.triangles
754 #-----------------------------------------------------------------------------
755 if __name__=="__main__":
756 try:
757 optlist,args = getopt.getopt(sys.argv[1:],"thdp")
758 except getopt.GetoptError:
759 usage()
760 sys.exit(2)
762 doHelp = 0
763 c = Context()
764 c.doPrint = 1
765 for opt in optlist:
766 if opt[0] == "-d": c.debug = 1
767 if opt[0] == "-p": c.plot = 1
768 if opt[0] == "-t": c.triangulate = 1
769 if opt[0] == "-h": doHelp = 1
771 if not doHelp:
772 pts = []
773 fp = sys.stdin
774 if len(args) > 0:
775 fp = open(args[0],'r')
776 for line in fp:
777 fld = line.split()
778 x = float(fld[0])
779 y = float(fld[1])
780 pts.append(Site(x,y))
781 if len(args) > 0: fp.close()
783 if doHelp or len(pts) == 0:
784 usage()
785 sys.exit(2)
787 sl = SiteList(pts)
788 voronoi(sl,c)