Index: voronoi.py =================================================================== --- voronoi.py (revision 13853) +++ voronoi.py (working copy) @@ -110,41 +110,49 @@ 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 -# self.extra_edges = [] # list of additional vertex 2-tubles (x,y) based on bounded voronoi tesselation -# self.set_bounds(None) -# self.use_bound = False + self.extra_edges = [] # list of additional vertex 2-tubles (x,y) based on bounded voronoi tesselation + self.set_bounds(None) + self.use_bound = False self.xmin = self.ymin = self.xmax = self.ymax = None + self.poly_vertices = [] # AJS def circle(self,x,y,rad): pass def clip_line(self,edge,lid,rid): - pass # here is where I will create false verticies if # the voronoi line extends to infinity... # the extra verticies will be added to the # extra edges list as 2-tuples -# a,b,c = edge.a,edge.b,edge.c -# if lid == -1: -# x = self.xMin -# y = (c-a*x) / b -# if y < self.yMin or y > self.yMax: -# if y < self.yMin: y = self.yMin -# elif y > self.yMax: y = self.yMax -# x = (c-b*y) / a -# self.extra_edges.append((x,y)) -# lid = -(len(self.extra_edges)-1) -# if rid == -1: -# x = self.xMax -# y = (c-a*x) / b -# if y < self.yMin or y > self.yMax: -# if y < self.yMin: y = self.yMin -# elif y > self.yMax: y = self.yMax -# x = (c-b*y) / a -# self.extra_edges.append((x,y)) -# rid = -(len(self.extra_edges)-1) -# print lid,rid -# return (lid,rid) + a,b,c = edge.a,edge.b,edge.c + if lid == -1: + x = self.xmin + y = (c-a*x) / b + if y < self.ymin or y > self.ymax: + if y < self.ymin: y = self.ymin + elif y > self.ymax: y = self.ymax + x = (c-b*y) / a + self.extra_edges.append((x,y)) + ang0 = math.atan2((y - edge.reg[0].y),(x - edge.reg[0].x)) # AJS Add false vertices to list of polygon vertices + ang1 = math.atan2((y - edge.reg[1].y),(x - edge.reg[1].x)) # AJS + self.poly_vertices.append((edge.reg[0].sitenum,ang0,x,y)) # AJS + self.poly_vertices.append((edge.reg[1].sitenum,ang1,x,y)) # AJS + lid = -(len(self.extra_edges)-1) + if rid == -1: + x = self.xmax + y = (c-a*x) / b + if y < self.ymin or y > self.ymax: + if y < self.ymin: y = self.ymin + elif y > self.ymax: y = self.ymax + x = (c-b*y) / a + self.extra_edges.append((x,y)) + ang0 = math.atan2((y - edge.reg[0].y),(x - edge.reg[0].x)) # AJS As above + ang1 = math.atan2((y - edge.reg[1].y),(x - edge.reg[1].x)) # AJS + self.poly_vertices.append((edge.reg[0].sitenum,ang0,x,y)) # AJS + self.poly_vertices.append((edge.reg[1].sitenum,ang1,x,y)) # AJS + rid = -(len(self.extra_edges)-1) + print lid,rid + return (lid,rid) def line(self,x0,y0,x1,y1): pass @@ -159,7 +167,10 @@ elif(self.doPrint): print "s %f %f" % (s.x, s.y) - def outVertex(self,s): + def outVertex(self,s,bot,top,mid,botang,topang,midang): # AJS + self.poly_vertices.append((bot,botang,s.x,s.y)) # AJS + self.poly_vertices.append((top,topang,s.x,s.y)) # AJS + self.poly_vertices.append((mid,midang,s.x,s.y)) # AJS self.vertices.append((s.x,s.y)) if s.x < self.xmin: self.xmin = s.x elif s.x > self.xmax: self.xmax = s.x @@ -198,8 +209,8 @@ if edge.ep[Edge.RE] is not None: sitenumR = edge.ep[Edge.RE].sitenum -# if sitenumL == -1 or sitenumR == -1 and self.use_bound: -# sitenumL,sitenumR = self.clip_line(edge,sitenumL,sitenumR) + if (sitenumL == -1 or sitenumR == -1) and self.use_bound: # AJS (Logic error - OR condition need braces) + sitenumL,sitenumR = self.clip_line(edge,sitenumL,sitenumR) # AJS self.edges.append((edge.edgenum,sitenumL,sitenumR)) if(not self.triangulate): @@ -299,7 +310,10 @@ # couldn't do this earlier since we didn't know when it would be processed v = lbnd.vertex siteList.setSiteNumber(v) - context.outVertex(v) + botang = math.atan2((v.y - bot.y),(v.x - bot.x)) # AJS Calculate angles from vertex to source layer points + topang = math.atan2((v.y - top.y),(v.x - top.x)) # AJS + midang = math.atan2((v.y - mid.y),(v.x - mid.x)) # AJS + context.outVertex(v,bot.sitenum,top.sitenum,mid.sitenum,botang,topang,midang) # AJS # set the endpoint of the left and right Halfedge to be this vector if lbnd.edge.setEndpoint(lbnd.pm,v): @@ -784,9 +798,10 @@ """ siteList = SiteList( points ) context = Context() + context.use_bound = True # AJS context.set_bounds( siteList ) voronoi( siteList, context ) - return ( context.vertices, context.lines, context.edges, (context.xmin,context.ymin,context.xmax,context.ymax)) + return ( context.vertices, context.lines, context.edges, (context.xmin,context.ymin,context.xmax,context.ymax),context.poly_vertices) # AJS #------------------------------------------------------------------ def computeDelaunayTriangulation( points ): Index: doGeometry.py =================================================================== --- doGeometry.py (revision 13853) +++ doGeometry.py (working copy) @@ -457,8 +457,61 @@ del writer return True + def voronoi_polygons( self ): + import voronoi + + ptDict = {} + ptNdx = -1 + vprovider = self.vlayer.dataProvider() + allAttrs = vprovider.attributeIndexes() + vprovider.select( allAttrs ) + fields = { + 0 : QgsField( "PointId", QVariant.Int ), + 1 : QgsField( "PointName", QVariant.String )} + writer = QgsVectorFileWriter( self.myName, self.myEncoding, + fields, QGis.WKBPolygon, vprovider.crs() ) + inFeat = QgsFeature() + points = [] + while vprovider.nextFeature( inFeat ): + inGeom = QgsGeometry( inFeat.geometry() ) + point = inGeom.asPoint() + points.append( point ) + ptNdx +=1 + ptDict[ptNdx] = inFeat.id() + vprovider.rewind() + vprovider.select( allAttrs ) + vertices,equations,edges,bounds,polygon_vertices = voronoi.computeVoronoiDiagram(points) + feat = QgsFeature() + nFeat = len( vertices ) + nElement = 0 + self.emit( SIGNAL( "runStatus(PyQt_PyObject)" ), 0 ) + self.emit( SIGNAL( "runRange(PyQt_PyObject)" ), ( 0, nFeat ) ) + + polygon_vertices.sort() + cur = -1 + for i,v in enumerate(polygon_vertices): + if (not(v[0] == cur) and i > 0) or (i == len(polygon_vertices) - 1): + polygon.append(curStart) + feat.addAttribute( 0, QVariant( ptDict[ cur ] ) ) + geometry = QgsGeometry().fromPolygon( [ polygon ] ) + feat.setGeometry( geometry ) + writer.addFeature( feat ) + nElement += 1 + self.emit( SIGNAL( "runStatus(PyQt_PyObject)" ), nElement ) + if (not(v[0] == cur)) : + cur = v[0] + curStart = QgsPoint(v[2],v[3]) + polygon = [] + polygon.append( QgsPoint(v[2],v[3]) ) + + del writer + return True + def delaunay_triangulation( self ): import voronoi + + ptDict = {} + ptNdx = -1 vprovider = self.vlayer.dataProvider() allAttrs = vprovider.attributeIndexes() vprovider.select( allAttrs ) @@ -474,6 +527,8 @@ inGeom = QgsGeometry( inFeat.geometry() ) point = inGeom.asPoint() points.append( point ) + ptNdx +=1 + ptDict[ptNdx] = inFeat.id() vprovider.rewind() vprovider.select( allAttrs ) triangles = voronoi.computeDelaunayTriangulation( points ) @@ -488,7 +543,7 @@ polygon = [] step = 0 for index in indicies: - vprovider.featureAtId( index, inFeat, True, allAttrs ) + vprovider.featureAtId( ptDict[index], inFeat, True, allAttrs ) geom = QgsGeometry( inFeat.geometry() ) point = QgsPoint( geom.asPoint() ) polygon.append( point )