From: Alvin Penner Date: Thu, 21 Jan 2010 01:07:47 +0000 (-0500) Subject: Generate Voronoi extension. Code by Steven Fortune (C) and Bill Simons (Python). X-Git-Url: https://git.tokkee.org/?a=commitdiff_plain;h=d9b1d836390f0891fb0d18f7894a385cd1b38020;p=inkscape.git Generate Voronoi extension. Code by Steven Fortune (C) and Bill Simons (Python). --- diff --git a/share/extensions/Makefile.am b/share/extensions/Makefile.am index 1650923e0..33dfcc25c 100644 --- a/share/extensions/Makefile.am +++ b/share/extensions/Makefile.am @@ -61,6 +61,7 @@ extensions = \ fractalize.py \ funcplot.py \ gears.py\ + generate_voronoi.py \ gimp_xcf.py \ grid_cartesian.py \ grid_polar.py \ @@ -129,6 +130,7 @@ extensions = \ txt2svg.pl \ uniconv-ext.py \ uniconv_output.py \ + voronoi.py \ web-set-att.py \ web-transmit-att.py \ whirl.py \ @@ -190,6 +192,7 @@ modules = \ fractalize.inx \ funcplot.inx \ gears.inx\ + generate_voronoi.inx \ gimp_xcf.inx \ grid_cartesian.inx \ grid_polar.inx \ diff --git a/share/extensions/generate_voronoi.inx b/share/extensions/generate_voronoi.inx new file mode 100644 index 000000000..4a04ffe1d --- /dev/null +++ b/share/extensions/generate_voronoi.inx @@ -0,0 +1,21 @@ + + + <_name>Voronoi Pattern + com.vaxxine.generate.voronoi + generate_voronoi.py + voronoi.py + inkex.py + <_param name="title1" type="description">Generate a random pattern of Voronoi cells. The pattern will be accessible in the Fill and Stroke dialog. You must select an object or a group. + <_param name="title2" type="description">If border is zero, the pattern will be discontinuous at the edges. Use a positive border, preferably greater than the cell size, to produce a smooth join of the pattern at the edges. Use a negative border to reduce the size of the pattern and get an empty border. + 10 + 0 + + all + + + + + + diff --git a/share/extensions/generate_voronoi.py b/share/extensions/generate_voronoi.py new file mode 100644 index 000000000..3359685fc --- /dev/null +++ b/share/extensions/generate_voronoi.py @@ -0,0 +1,187 @@ +#!/usr/bin/env python +""" +Copyright (C) 2010 Alvin Penner, penner@vaxxine.com + +- Voronoi Diagram algorithm and C code by Steven Fortune, 1987, http://ect.bell-labs.com/who/sjf/ +- Python translation to file voronoi.py by Bill Simons, 2005, http://www.oxfish.com/ + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +""" +import random, inkex, simplestyle, gettext, voronoi +_ = gettext.gettext + +try: + from subprocess import Popen, PIPE +except: + inkex.errormsg(_("Failed to import the subprocess module. Please report this as a bug at : https://bugs.launchpad.net/inkscape.")) + inkex.errormsg("Python version is : " + str(inkex.sys.version_info)) + exit() + +def clip_line(x1, y1, x2, y2, w, h): + if x1 < 0 and x2 < 0: + return [0, 0, 0, 0] + if x1 > w and x2 > w: + return [0, 0, 0, 0] + if x1 < 0: + y1 = (y1*x2 - y2*x1)/(x2 - x1) + x1 = 0 + if x2 < 0: + y2 = (y1*x2 - y2*x1)/(x2 - x1) + x2 = 0 + if x1 > w: + y1 = y1 + (w - x1)*(y2 - y1)/(x2 - x1) + x1 = w + if x2 > w: + y2 = y1 + (w - x1)*(y2 - y1)/(x2 - x1) + x2 = w + if y1 < 0 and y2 < 0: + return [0, 0, 0, 0] + if y1 > h and y2 > h: + return [0, 0, 0, 0] + if x1 == x2 and y1 == y2: + return [0, 0, 0, 0] + if y1 < 0: + x1 = (x1*y2 - x2*y1)/(y2 - y1) + y1 = 0 + if y2 < 0: + x2 = (x1*y2 - x2*y1)/(y2 - y1) + y2 = 0 + if y1 > h: + x1 = x1 + (h - y1)*(x2 - x1)/(y2 - y1) + y1 = h + if y2 > h: + x2 = x1 + (h - y1)*(x2 - x1)/(y2 - y1) + y2 = h + return [x1, y1, x2, y2] + +class Pattern(inkex.Effect): + def __init__(self): + inkex.Effect.__init__(self) + self.OptionParser.add_option("--size", + action="store", type="int", + dest="size", default=10, + help="Average size of cell (px)") + self.OptionParser.add_option("--border", + action="store", type="int", + dest="border", default=0, + help="Size of Border (px)") + + def effect(self): + if not self.options.ids: + inkex.errormsg(_("Please select an object")) + exit() + q = {'x':0,'y':0,'width':0,'height':0} # query the bounding box of ids[0] + for query in q.keys(): + p = Popen('inkscape --query-%s --query-id=%s "%s"' % (query, self.options.ids[0], self.args[-1]), shell=True, stdout=PIPE, stderr=PIPE) + rc = p.wait() + q[query] = float(p.stdout.read()) + defs = self.xpathSingle('/svg:svg//svg:defs') + pattern = inkex.etree.SubElement(defs ,inkex.addNS('pattern','svg')) + pattern.set('id', 'Voronoi' + str(random.randint(1, 9999))) + pattern.set('width', str(q['width'])) + pattern.set('height', str(q['height'])) + pattern.set('patternTransform', 'translate(%s,%s)' % (q['x'], q['y'])) + pattern.set('patternUnits', 'userSpaceOnUse') + + # generate random pattern of points + c = voronoi.Context() + pts = [] + b = float(self.options.border) # width of border + for i in range(int(q['width']*q['height']/self.options.size/self.options.size)): + x = random.random()*q['width'] + y = random.random()*q['height'] + if b > 0: # duplicate border area + pts.append(voronoi.Site(x, y)) + if x < b: + pts.append(voronoi.Site(x + q['width'], y)) + if y < b: + pts.append(voronoi.Site(x + q['width'], y + q['height'])) + if y > q['height'] - b: + pts.append(voronoi.Site(x + q['width'], y - q['height'])) + if x > q['width'] - b: + pts.append(voronoi.Site(x - q['width'], y)) + if y < b: + pts.append(voronoi.Site(x - q['width'], y + q['height'])) + if y > q['height'] - b: + pts.append(voronoi.Site(x - q['width'], y - q['height'])) + if y < b: + pts.append(voronoi.Site(x, y + q['height'])) + if y > q['height'] - b: + pts.append(voronoi.Site(x, y - q['height'])) + elif x > -b and y > -b and x < q['width'] + b and y < q['height'] + b: + pts.append(voronoi.Site(x, y)) # leave border area blank + # dot = inkex.etree.SubElement(pattern, inkex.addNS('rect','svg')) + # dot.set('x', str(x-1)) + # dot.set('y', str(y-1)) + # dot.set('width', '2') + # dot.set('height', '2') + if len(pts) < 3: + inkex.errormsg("Please choose a larger object, or smaller cell size") + exit() + + # plot Voronoi diagram + sl = voronoi.SiteList(pts) + voronoi.voronoi(sl, c) + for edge in c.edges: + if edge[1] >= 0 and edge[2] >= 0: # two vertices + [x1, y1, x2, y2] = clip_line(c.vertices[edge[1]][0], c.vertices[edge[1]][1], c.vertices[edge[2]][0], c.vertices[edge[2]][1], q['width'], q['height']) + elif edge[1] >= 0: # only one vertex + if c.lines[edge[0]][1] == 0: # vertical line + xtemp = c.lines[edge[0]][2]/c.lines[edge[0]][0] + if c.vertices[edge[1]][1] > q['height']/2: + ytemp = q['height'] + else: + ytemp = 0 + else: + xtemp = q['width'] + ytemp = (c.lines[edge[0]][2] - q['width']*c.lines[edge[0]][0])/c.lines[edge[0]][1] + [x1, y1, x2, y2] = clip_line(c.vertices[edge[1]][0], c.vertices[edge[1]][1], xtemp, ytemp, q['width'], q['height']) + elif edge[2] >= 0: # only one vertex + if c.lines[edge[0]][1] == 0: # vertical line + xtemp = c.lines[edge[0]][2]/c.lines[edge[0]][0] + if c.vertices[edge[2]][1] > q['height']/2: + ytemp = q['height'] + else: + ytemp = 0 + else: + xtemp = 0 + ytemp = c.lines[edge[0]][2]/c.lines[edge[0]][1] + [x1, y1, x2, y2] = clip_line(xtemp, ytemp, c.vertices[edge[2]][0], c.vertices[edge[2]][1], q['width'], q['height']) + if x1 or x2 or y1 or y2: + path = 'M %f,%f %f,%f' % (x1, y1, x2, y2) + attribs = {'d': path, 'style': 'stroke:#000000'} + inkex.etree.SubElement(pattern, inkex.addNS('path', 'svg'), attribs) + + # link selected object to pattern + obj = self.selected[self.options.ids[0]] + style = {} + if obj.attrib.has_key('style'): + style = simplestyle.parseStyle(obj.attrib['style']) + style['fill'] = 'url(#%s)' % pattern.get('id') + obj.attrib['style'] = simplestyle.formatStyle(style) + if obj.tag == inkex.addNS('g', 'svg'): + for node in obj: + style = {} + if node.attrib.has_key('style'): + style = simplestyle.parseStyle(node.attrib['style']) + style['fill'] = 'url(#%s)' % pattern.get('id') + node.attrib['style'] = simplestyle.formatStyle(style) + +if __name__ == '__main__': + e = Pattern() + e.affect() + +# vim: expandtab shiftwidth=4 tabstop=8 softtabstop=4 encoding=utf-8 textwidth=99 diff --git a/share/extensions/voronoi.py b/share/extensions/voronoi.py new file mode 100644 index 000000000..be15fbe13 --- /dev/null +++ b/share/extensions/voronoi.py @@ -0,0 +1,789 @@ +############################################################################# +# +# Voronoi diagram calculator/ Delaunay triangulator +# Translated to Python by Bill Simons +# September, 2005 +# +# Calculate Delaunay triangulation or the Voronoi polygons for a set of +# 2D input points. +# +# Derived from code bearing the following notice: +# +# The author of this software is Steven Fortune. Copyright (c) 1994 by AT&T +# Bell Laboratories. +# Permission to use, copy, modify, and distribute this software for any +# purpose without fee is hereby granted, provided that this entire notice +# is included in all copies of any software which is or includes a copy +# or modification of this software and in all copies of the supporting +# documentation for such software. +# THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED +# WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY +# REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY +# OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. +# +# Comments were incorporated from Shane O'Sullivan's translation of the +# original code into C++ (http://mapviewer.skynet.ie/voronoi.html) +# +# Steve Fortune's homepage: http://netlib.bell-labs.com/cm/cs/who/sjf/index.html +# +############################################################################# + +def usage(): + print """ +voronoi - compute Voronoi diagram or Delaunay triangulation + +voronoi [-t -p -d] [filename] + +Voronoi reads from filename (or standard input if no filename given) for a set +of points in the plane and writes either the Voronoi diagram or the Delaunay +triangulation to the standard output. Each input line should consist of two +real numbers, separated by white space. + +If option -t is present, the Delaunay triangulation is produced. +Each output line is a triple i j k, which are the indices of the three points +in a Delaunay triangle. Points are numbered starting at 0. + +If option -t is not present, the Voronoi diagram is produced. +There are four output record types. + +s a b indicates that an input point at coordinates a b was seen. +l a b c indicates a line with equation ax + by = c. +v a b indicates a vertex at a b. +e l v1 v2 indicates a Voronoi segment which is a subsegment of line number l + with endpoints numbered v1 and v2. If v1 or v2 is -1, the line + extends to infinity. + +Other options include: + +d Print debugging info + +p Produce output suitable for input to plot (1), rather than the forms + described above. + +On unsorted data uniformly distributed in the unit square, voronoi uses about +20n+140 bytes of storage. + +AUTHOR +Steve J. Fortune (1987) A Sweepline Algorithm for Voronoi Diagrams, +Algorithmica 2, 153-174. +""" + +############################################################################# +# +# For programmatic use two functions are available: +# +# computeVoronoiDiagram(points) +# +# Takes a list of point objects (which must have x and y fields). +# Returns a 3-tuple of: +# +# (1) a list of 2-tuples, which are the x,y coordinates of the +# Voronoi diagram vertices +# (2) a list of 3-tuples (a,b,c) which are the equations of the +# lines in the Voronoi diagram: a*x + b*y = c +# (3) a list of 3-tuples, (l, v1, v2) representing edges of the +# Voronoi diagram. l is the index of the line, v1 and v2 are +# the indices of the vetices at the end of the edge. If +# v1 or v2 is -1, the line extends to infinity. +# +# computeDelaunayTriangulation(points): +# +# Takes a list of point objects (which must have x and y fields). +# Returns a list of 3-tuples: the indices of the points that form a +# Delaunay triangle. +# +############################################################################# +import math +import sys +import getopt +TOLERANCE = 1e-9 +BIG_FLOAT = 1e38 + +#------------------------------------------------------------------ +class Context(object): + def __init__(self): + self.doPrint = 0 + self.debug = 0 + self.plot = 0 + self.triangulate = False + self.vertices = [] # list of vertex 2-tuples: (x,y) + self.lines = [] # equation of line 3-tuple (a b c), for the equation of the line a*x+b*y = c + self.edges = [] # edge 3-tuple: (line index, vertex 1 index, vertex 2 index) if either vertex index is -1, the edge extends to infiinity + self.triangles = [] # 3-tuple of vertex indices + + def circle(self,x,y,rad): + pass + + def clip_line(self,edge): + pass + + def line(self,x0,y0,x1,y1): + pass + + def outSite(self,s): + if(self.debug): + print "site (%d) at %f %f" % (s.sitenum, s.x, s.y) + elif(self.triangulate): + pass + elif(self.plot): + self.circle (s.x, s.y, cradius) + elif(self.doPrint): + print "s %f %f" % (s.x, s.y) + + def outVertex(self,s): + self.vertices.append((s.x,s.y)) + if(self.debug): + print "vertex(%d) at %f %f" % (s.sitenum, s.x, s.y) + elif(self.triangulate): + pass + elif(self.doPrint and not self.plot): + print "v %f %f" % (s.x,s.y) + + def outTriple(self,s1,s2,s3): + self.triangles.append((s1.sitenum, s2.sitenum, s3.sitenum)) + if(self.debug): + print "circle through left=%d right=%d bottom=%d" % (s1.sitenum, s2.sitenum, s3.sitenum) + elif(self.triangulate and self.doPrint and not self.plot): + print "%d %d %d" % (s1.sitenum, s2.sitenum, s3.sitenum) + + def outBisector(self,edge): + self.lines.append((edge.a, edge.b, edge.c)) + if(self.debug): + 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) + elif(self.triangulate): + if(self.plot): + self.line(edge.reg[0].x, edge.reg[0].y, edge.reg[1].x, edge.reg[1].y) + elif(self.doPrint and not self.plot): + print "l %f %f %f" % (edge.a, edge.b, edge.c) + + def outEdge(self,edge): + sitenumL = -1 + if edge.ep[Edge.LE] is not None: + sitenumL = edge.ep[Edge.LE].sitenum + sitenumR = -1 + if edge.ep[Edge.RE] is not None: + sitenumR = edge.ep[Edge.RE].sitenum + self.edges.append((edge.edgenum,sitenumL,sitenumR)) + if(not self.triangulate): + if self.plot: + self.clip_line(edge) + elif(self.doPrint): + print "e %d" % edge.edgenum, + print " %d " % sitenumL, + print "%d" % sitenumR + +#------------------------------------------------------------------ +def voronoi(siteList,context): + edgeList = EdgeList(siteList.xmin,siteList.xmax,len(siteList)) + priorityQ = PriorityQueue(siteList.ymin,siteList.ymax,len(siteList)) + siteIter = siteList.iterator() + + bottomsite = siteIter.next() + context.outSite(bottomsite) + newsite = siteIter.next() + minpt = Site(-BIG_FLOAT,-BIG_FLOAT) + while True: + if not priorityQ.isEmpty(): + minpt = priorityQ.getMinPt() + + if (newsite and (priorityQ.isEmpty() or cmp(newsite,minpt) < 0)): + # newsite is smallest - this is a site event + context.outSite(newsite) + + # get first Halfedge to the LEFT and RIGHT of the new site + lbnd = edgeList.leftbnd(newsite) + rbnd = lbnd.right + + # if this halfedge has no edge, bot = bottom site (whatever that is) + # create a new edge that bisects + bot = lbnd.rightreg(bottomsite) + edge = Edge.bisect(bot,newsite) + context.outBisector(edge) + + # create a new Halfedge, setting its pm field to 0 and insert + # this new bisector edge between the left and right vectors in + # a linked list + bisector = Halfedge(edge,Edge.LE) + edgeList.insert(lbnd,bisector) + + # if the new bisector intersects with the left edge, remove + # the left edge's vertex, and put in the new one + p = lbnd.intersect(bisector) + if p is not None: + priorityQ.delete(lbnd) + priorityQ.insert(lbnd,p,newsite.distance(p)) + + # create a new Halfedge, setting its pm field to 1 + # insert the new Halfedge to the right of the original bisector + lbnd = bisector + bisector = Halfedge(edge,Edge.RE) + edgeList.insert(lbnd,bisector) + + # if this new bisector intersects with the right Halfedge + p = bisector.intersect(rbnd) + if p is not None: + # push the Halfedge into the ordered linked list of vertices + priorityQ.insert(bisector,p,newsite.distance(p)) + + newsite = siteIter.next() + + elif not priorityQ.isEmpty(): + # intersection is smallest - this is a vector (circle) event + + # pop the Halfedge with the lowest vector off the ordered list of + # vectors. Get the Halfedge to the left and right of the above HE + # and also the Halfedge to the right of the right HE + lbnd = priorityQ.popMinHalfedge() + llbnd = lbnd.left + rbnd = lbnd.right + rrbnd = rbnd.right + + # get the Site to the left of the left HE and to the right of + # the right HE which it bisects + bot = lbnd.leftreg(bottomsite) + top = rbnd.rightreg(bottomsite) + + # output the triple of sites, stating that a circle goes through them + mid = lbnd.rightreg(bottomsite) + context.outTriple(bot,top,mid) + + # get the vertex that caused this event and set the vertex number + # couldn't do this earlier since we didn't know when it would be processed + v = lbnd.vertex + siteList.setSiteNumber(v) + context.outVertex(v) + + # set the endpoint of the left and right Halfedge to be this vector + if lbnd.edge.setEndpoint(lbnd.pm,v): + context.outEdge(lbnd.edge) + + if rbnd.edge.setEndpoint(rbnd.pm,v): + context.outEdge(rbnd.edge) + + + # delete the lowest HE, remove all vertex events to do with the + # right HE and delete the right HE + edgeList.delete(lbnd) + priorityQ.delete(rbnd) + edgeList.delete(rbnd) + + + # if the site to the left of the event is higher than the Site + # to the right of it, then swap them and set 'pm' to RIGHT + pm = Edge.LE + if bot.y > top.y: + bot,top = top,bot + pm = Edge.RE + + # Create an Edge (or line) that is between the two Sites. This + # creates the formula of the line, and assigns a line number to it + edge = Edge.bisect(bot, top) + context.outBisector(edge) + + # create a HE from the edge + bisector = Halfedge(edge, pm) + + # insert the new bisector to the right of the left HE + # set one endpoint to the new edge to be the vector point 'v' + # If the site to the left of this bisector is higher than the right + # Site, then this endpoint is put in position 0; otherwise in pos 1 + edgeList.insert(llbnd, bisector) + if edge.setEndpoint(Edge.RE - pm, v): + context.outEdge(edge) + + # if left HE and the new bisector don't intersect, then delete + # the left HE, and reinsert it + p = llbnd.intersect(bisector) + if p is not None: + priorityQ.delete(llbnd); + priorityQ.insert(llbnd, p, bot.distance(p)) + + # if right HE and the new bisector don't intersect, then reinsert it + p = bisector.intersect(rrbnd) + if p is not None: + priorityQ.insert(bisector, p, bot.distance(p)) + else: + break + + he = edgeList.leftend.right + while he is not edgeList.rightend: + context.outEdge(he.edge) + he = he.right + +#------------------------------------------------------------------ +def isEqual(a,b,relativeError=TOLERANCE): + # is nearly equal to within the allowed relative error + norm = max(abs(a),abs(b)) + return (norm < relativeError) or (abs(a - b) < (relativeError * norm)) + +#------------------------------------------------------------------ +class Site(object): + def __init__(self,x=0.0,y=0.0,sitenum=0): + self.x = x + self.y = y + self.sitenum = sitenum + + def dump(self): + print "Site #%d (%g, %g)" % (self.sitenum,self.x,self.y) + + def __cmp__(self,other): + if self.y < other.y: + return -1 + elif self.y > other.y: + return 1 + elif self.x < other.x: + return -1 + elif self.x > other.x: + return 1 + else: + return 0 + + def distance(self,other): + dx = self.x - other.x + dy = self.y - other.y + return math.sqrt(dx*dx + dy*dy) + +#------------------------------------------------------------------ +class Edge(object): + LE = 0 + RE = 1 + EDGE_NUM = 0 + DELETED = {} # marker value + + def __init__(self): + self.a = 0.0 + self.b = 0.0 + self.c = 0.0 + self.ep = [None,None] + self.reg = [None,None] + self.edgenum = 0 + + def dump(self): + print "(#%d a=%g, b=%g, c=%g)" % (self.edgenum,self.a,self.b,self.c) + print "ep",self.ep + print "reg",self.reg + + def setEndpoint(self, lrFlag, site): + self.ep[lrFlag] = site + if self.ep[Edge.RE - lrFlag] is None: + return False + return True + + @staticmethod + def bisect(s1,s2): + newedge = Edge() + newedge.reg[0] = s1 # store the sites that this edge is bisecting + newedge.reg[1] = s2 + + # to begin with, there are no endpoints on the bisector - it goes to infinity + # ep[0] and ep[1] are None + + # get the difference in x dist between the sites + dx = float(s2.x - s1.x) + dy = float(s2.y - s1.y) + adx = abs(dx) # make sure that the difference in positive + ady = abs(dy) + + # get the slope of the line + newedge.c = float(s1.x * dx + s1.y * dy + (dx*dx + dy*dy)*0.5) + if adx > ady : + # set formula of line, with x fixed to 1 + newedge.a = 1.0 + newedge.b = dy/dx + newedge.c /= dx + else: + # set formula of line, with y fixed to 1 + newedge.b = 1.0 + newedge.a = dx/dy + newedge.c /= dy + + newedge.edgenum = Edge.EDGE_NUM + Edge.EDGE_NUM += 1 + return newedge + + +#------------------------------------------------------------------ +class Halfedge(object): + def __init__(self,edge=None,pm=Edge.LE): + self.left = None # left Halfedge in the edge list + self.right = None # right Halfedge in the edge list + self.qnext = None # priority queue linked list pointer + self.edge = edge # edge list Edge + self.pm = pm + self.vertex = None # Site() + self.ystar = BIG_FLOAT + + def dump(self): + print "Halfedge--------------------------" + print "left: ", self.left + print "right: ", self.right + print "edge: ", self.edge + print "pm: ", self.pm + print "vertex: ", + if self.vertex: self.vertex.dump() + else: print "None" + print "ystar: ", self.ystar + + + def __cmp__(self,other): + if self.ystar > other.ystar: + return 1 + elif self.ystar < other.ystar: + return -1 + elif self.vertex.x > other.vertex.x: + return 1 + elif self.vertex.x < other.vertex.x: + return -1 + else: + return 0 + + def leftreg(self,default): + if not self.edge: + return default + elif self.pm == Edge.LE: + return self.edge.reg[Edge.LE] + else: + return self.edge.reg[Edge.RE] + + def rightreg(self,default): + if not self.edge: + return default + elif self.pm == Edge.LE: + return self.edge.reg[Edge.RE] + else: + return self.edge.reg[Edge.LE] + + + # returns True if p is to right of halfedge self + def isPointRightOf(self,pt): + e = self.edge + topsite = e.reg[1] + right_of_site = pt.x > topsite.x + + if(right_of_site and self.pm == Edge.LE): + return True + + if(not right_of_site and self.pm == Edge.RE): + return False + + if(e.a == 1.0): + dyp = pt.y - topsite.y + dxp = pt.x - topsite.x + fast = 0; + if ((not right_of_site and e.b < 0.0) or (right_of_site and e.b >= 0.0)): + above = dyp >= e.b * dxp + fast = above + else: + above = pt.x + pt.y * e.b > e.c + if(e.b < 0.0): + above = not above + if (not above): + fast = 1 + if (not fast): + dxs = topsite.x - (e.reg[0]).x + above = e.b * (dxp*dxp - dyp*dyp) < dxs*dyp*(1.0+2.0*dxp/dxs + e.b*e.b) + if(e.b < 0.0): + above = not above + else: # e.b == 1.0 + yl = e.c - e.a * pt.x + t1 = pt.y - yl + t2 = pt.x - topsite.x + t3 = yl - topsite.y + above = t1*t1 > t2*t2 + t3*t3 + + if(self.pm==Edge.LE): + return above + else: + return not above + + #-------------------------- + # create a new site where the Halfedges el1 and el2 intersect + def intersect(self,other): + e1 = self.edge + e2 = other.edge + if (e1 is None) or (e2 is None): + return None + + # if the two edges bisect the same parent return None + if e1.reg[1] is e2.reg[1]: + return None + + d = e1.a * e2.b - e1.b * e2.a + if isEqual(d,0.0): + return None + + xint = (e1.c*e2.b - e2.c*e1.b) / d + yint = (e2.c*e1.a - e1.c*e2.a) / d + if(cmp(e1.reg[1],e2.reg[1]) < 0): + he = self + e = e1 + else: + he = other + e = e2 + + rightOfSite = xint >= e.reg[1].x + if((rightOfSite and he.pm == Edge.LE) or + (not rightOfSite and he.pm == Edge.RE)): + return None + + # create a new site at the point of intersection - this is a new + # vector event waiting to happen + return Site(xint,yint) + + + +#------------------------------------------------------------------ +class EdgeList(object): + def __init__(self,xmin,xmax,nsites): + if xmin > xmax: xmin,xmax = xmax,xmin + self.hashsize = int(2*math.sqrt(nsites+4)) + + self.xmin = xmin + self.deltax = float(xmax - xmin) + self.hash = [None]*self.hashsize + + self.leftend = Halfedge() + self.rightend = Halfedge() + self.leftend.right = self.rightend + self.rightend.left = self.leftend + self.hash[0] = self.leftend + self.hash[-1] = self.rightend + + def insert(self,left,he): + he.left = left + he.right = left.right + left.right.left = he + left.right = he + + def delete(self,he): + he.left.right = he.right + he.right.left = he.left + he.edge = Edge.DELETED + + # Get entry from hash table, pruning any deleted nodes + def gethash(self,b): + if(b < 0 or b >= self.hashsize): + return None + he = self.hash[b] + if he is None or he.edge is not Edge.DELETED: + return he + + # Hash table points to deleted half edge. Patch as necessary. + self.hash[b] = None + return None + + def leftbnd(self,pt): + # Use hash table to get close to desired halfedge + bucket = int(((pt.x - self.xmin)/self.deltax * self.hashsize)) + + if(bucket < 0): + bucket =0; + + if(bucket >=self.hashsize): + bucket = self.hashsize-1 + + he = self.gethash(bucket) + if(he is None): + i = 1 + while True: + he = self.gethash(bucket-i) + if (he is not None): break; + he = self.gethash(bucket+i) + if (he is not None): break; + i += 1 + + # Now search linear list of halfedges for the corect one + if (he is self.leftend) or (he is not self.rightend and he.isPointRightOf(pt)): + he = he.right + while he is not self.rightend and he.isPointRightOf(pt): + he = he.right + he = he.left; + else: + he = he.left + while (he is not self.leftend and not he.isPointRightOf(pt)): + he = he.left + + # Update hash table and reference counts + if(bucket > 0 and bucket < self.hashsize-1): + self.hash[bucket] = he + return he + + +#------------------------------------------------------------------ +class PriorityQueue(object): + def __init__(self,ymin,ymax,nsites): + self.ymin = ymin + self.deltay = ymax - ymin + self.hashsize = int(4 * math.sqrt(nsites)) + self.count = 0 + self.minidx = 0 + self.hash = [] + for i in range(self.hashsize): + self.hash.append(Halfedge()) + + def __len__(self): + return self.count + + def isEmpty(self): + return self.count == 0 + + def insert(self,he,site,offset): + he.vertex = site + he.ystar = site.y + offset + last = self.hash[self.getBucket(he)] + next = last.qnext + while((next is not None) and cmp(he,next) > 0): + last = next + next = last.qnext + he.qnext = last.qnext + last.qnext = he + self.count += 1 + + def delete(self,he): + if (he.vertex is not None): + last = self.hash[self.getBucket(he)] + while last.qnext is not he: + last = last.qnext + last.qnext = he.qnext + self.count -= 1 + he.vertex = None + + def getBucket(self,he): + bucket = int(((he.ystar - self.ymin) / self.deltay) * self.hashsize) + if bucket < 0: bucket = 0 + if bucket >= self.hashsize: bucket = self.hashsize-1 + if bucket < self.minidx: self.minidx = bucket + return bucket + + def getMinPt(self): + while(self.hash[self.minidx].qnext is None): + self.minidx += 1 + he = self.hash[self.minidx].qnext + x = he.vertex.x + y = he.ystar + return Site(x,y) + + def popMinHalfedge(self): + curr = self.hash[self.minidx].qnext + self.hash[self.minidx].qnext = curr.qnext + self.count -= 1 + return curr + + +#------------------------------------------------------------------ +class SiteList(object): + def __init__(self,pointList): + self.__sites = [] + self.__sitenum = 0 + + self.__xmin = pointList[0].x + self.__ymin = pointList[0].y + self.__xmax = pointList[0].x + self.__ymax = pointList[0].y + for i,pt in enumerate(pointList): + self.__sites.append(Site(pt.x,pt.y,i)) + if pt.x < self.__xmin: self.__xmin = pt.x + if pt.y < self.__ymin: self.__ymin = pt.y + if pt.x > self.__xmax: self.__xmax = pt.x + if pt.y > self.__ymax: self.__ymax = pt.y + self.__sites.sort() + + def setSiteNumber(self,site): + site.sitenum = self.__sitenum + self.__sitenum += 1 + + class Iterator(object): + def __init__(this,lst): this.generator = (s for s in lst) + def __iter__(this): return this + def next(this): + try: + return this.generator.next() + except StopIteration: + return None + + def iterator(self): + return SiteList.Iterator(self.__sites) + + def __iter__(self): + return SiteList.Iterator(self.__sites) + + def __len__(self): + return len(self.__sites) + + def _getxmin(self): return self.__xmin + def _getymin(self): return self.__ymin + def _getxmax(self): return self.__xmax + def _getymax(self): return self.__ymax + xmin = property(_getxmin) + ymin = property(_getymin) + xmax = property(_getxmax) + ymax = property(_getymax) + + +#------------------------------------------------------------------ +def computeVoronoiDiagram(points): + """ Takes a list of point objects (which must have x and y fields). + Returns a 3-tuple of: + + (1) a list of 2-tuples, which are the x,y coordinates of the + Voronoi diagram vertices + (2) a list of 3-tuples (a,b,c) which are the equations of the + lines in the Voronoi diagram: a*x + b*y = c + (3) a list of 3-tuples, (l, v1, v2) representing edges of the + Voronoi diagram. l is the index of the line, v1 and v2 are + the indices of the vetices at the end of the edge. If + v1 or v2 is -1, the line extends to infinity. + """ + siteList = SiteList(points) + context = Context() + voronoi(siteList,context) + return (context.vertices,context.lines,context.edges) + +#------------------------------------------------------------------ +def computeDelaunayTriangulation(points): + """ Takes a list of point objects (which must have x and y fields). + Returns a list of 3-tuples: the indices of the points that form a + Delaunay triangle. + """ + siteList = SiteList(points) + context = Context() + context.triangulate = true + voronoi(siteList,context) + return context.triangles + +#----------------------------------------------------------------------------- +if __name__=="__main__": + try: + optlist,args = getopt.getopt(sys.argv[1:],"thdp") + except getopt.GetoptError: + usage() + sys.exit(2) + + doHelp = 0 + c = Context() + c.doPrint = 1 + for opt in optlist: + if opt[0] == "-d": c.debug = 1 + if opt[0] == "-p": c.plot = 1 + if opt[0] == "-t": c.triangulate = 1 + if opt[0] == "-h": doHelp = 1 + + if not doHelp: + pts = [] + fp = sys.stdin + if len(args) > 0: + fp = open(args[0],'r') + for line in fp: + fld = line.split() + x = float(fld[0]) + y = float(fld[1]) + pts.append(Site(x,y)) + if len(args) > 0: fp.close() + + if doHelp or len(pts) == 0: + usage() + sys.exit(2) + + sl = SiteList(pts) + voronoi(sl,c) +