patch2825.diff

changes as patch (untested) - Jürgen Fischer, 2010-06-30 11:55 AM

Download (8.58 KB)

View differences:

voronoi.py (working copy)
110 110
        self.lines     = []    # equation of line 3-tuple (a b c), for the equation of the line a*x+b*y = c  
111 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 112
        self.triangles = []    # 3-tuple of vertex indices
113
#        self.extra_edges = []  # list of additional vertex 2-tubles (x,y) based on bounded voronoi tesselation
114
#        self.set_bounds(None)
115
#        self.use_bound = False
113
        self.extra_edges = []  # list of additional vertex 2-tubles (x,y) based on bounded voronoi tesselation
114
        self.set_bounds(None)
115
        self.use_bound = False
116 116
        self.xmin = self.ymin = self.xmax = self.ymax = None
117
        self.poly_vertices = [] # AJS
117 118

  
118 119
    def circle(self,x,y,rad):
119 120
        pass
120 121

  
121 122
    def clip_line(self,edge,lid,rid):
122
        pass
123 123
        # here is where I will create false verticies if
124 124
        # the voronoi line extends to infinity...
125 125
        # the extra verticies will be added to the
126 126
        # extra edges list as 2-tuples
127
#        a,b,c = edge.a,edge.b,edge.c
128
#        if lid == -1:
129
#            x = self.xMin
130
#            y = (c-a*x) / b
131
#            if y < self.yMin or y > self.yMax:
132
#                if y < self.yMin: y = self.yMin
133
#                elif y > self.yMax: y = self.yMax
134
#                x = (c-b*y) / a
135
#            self.extra_edges.append((x,y))
136
#            lid = -(len(self.extra_edges)-1)
137
#        if rid == -1:
138
#            x = self.xMax
139
#            y = (c-a*x) / b
140
#            if y < self.yMin or y > self.yMax:
141
#                if y < self.yMin: y = self.yMin
142
#                elif y > self.yMax: y = self.yMax
143
#                x = (c-b*y) / a
144
#            self.extra_edges.append((x,y))
145
#            rid = -(len(self.extra_edges)-1)
146
#        print lid,rid
147
#        return (lid,rid)
127
        a,b,c = edge.a,edge.b,edge.c
128
        if lid == -1:
129
            x = self.xmin
130
            y = (c-a*x) / b
131
            if y < self.ymin or y > self.ymax:
132
                if y < self.ymin: y = self.ymin
133
                elif y > self.ymax: y = self.ymax
134
                x = (c-b*y) / a
135
            self.extra_edges.append((x,y))
136
            ang0 = math.atan2((y - edge.reg[0].y),(x - edge.reg[0].x))    # AJS Add false vertices to list of polygon vertices
137
            ang1 = math.atan2((y - edge.reg[1].y),(x - edge.reg[1].x))    # AJS  
138
            self.poly_vertices.append((edge.reg[0].sitenum,ang0,x,y))	# AJS 
139
            self.poly_vertices.append((edge.reg[1].sitenum,ang1,x,y))	# AJS
140
            lid = -(len(self.extra_edges)-1)
141
        if rid == -1:
142
            x = self.xmax
143
            y = (c-a*x) / b
144
            if y < self.ymin or y > self.ymax:
145
                if y < self.ymin: y = self.ymin
146
                elif y > self.ymax: y = self.ymax
147
                x = (c-b*y) / a
148
            self.extra_edges.append((x,y))
149
            ang0 = math.atan2((y - edge.reg[0].y),(x - edge.reg[0].x))    # AJS As above
150
            ang1 = math.atan2((y - edge.reg[1].y),(x - edge.reg[1].x))    # AJS
151
            self.poly_vertices.append((edge.reg[0].sitenum,ang0,x,y))	# AJS
152
            self.poly_vertices.append((edge.reg[1].sitenum,ang1,x,y))	# AJS
153
            rid = -(len(self.extra_edges)-1)
154
        print lid,rid
155
        return (lid,rid)
148 156

  
149 157
    def line(self,x0,y0,x1,y1):
150 158
        pass
......
159 167
        elif(self.doPrint):
160 168
            print "s %f %f" % (s.x, s.y)
161 169

  
162
    def outVertex(self,s):
170
    def outVertex(self,s,bot,top,mid,botang,topang,midang):	# AJS
171
        self.poly_vertices.append((bot,botang,s.x,s.y))	# AJS
172
        self.poly_vertices.append((top,topang,s.x,s.y))	# AJS
173
        self.poly_vertices.append((mid,midang,s.x,s.y))	# AJS
163 174
        self.vertices.append((s.x,s.y))
164 175
        if s.x < self.xmin: self.xmin = s.x
165 176
        elif s.x > self.xmax: self.xmax = s.x
......
198 209
        if edge.ep[Edge.RE] is not None:
199 210
            sitenumR = edge.ep[Edge.RE].sitenum
200 211
            
201
#        if sitenumL == -1 or sitenumR == -1 and self.use_bound:
202
#            sitenumL,sitenumR = self.clip_line(edge,sitenumL,sitenumR)
212
        if (sitenumL == -1 or sitenumR == -1) and self.use_bound:	# AJS (Logic error - OR condition need braces)
213
            sitenumL,sitenumR = self.clip_line(edge,sitenumL,sitenumR)	# AJS
203 214

  
204 215
        self.edges.append((edge.edgenum,sitenumL,sitenumR))
205 216
        if(not self.triangulate):
......
299 310
            # couldn't do this earlier since we didn't know when it would be processed
300 311
            v = lbnd.vertex                 
301 312
            siteList.setSiteNumber(v)
302
            context.outVertex(v)
313
            botang = math.atan2((v.y - bot.y),(v.x - bot.x))    # AJS Calculate angles from vertex to source layer points
314
            topang = math.atan2((v.y - top.y),(v.x - top.x))    # AJS
315
            midang = math.atan2((v.y - mid.y),(v.x - mid.x))    # AJS
316
            context.outVertex(v,bot.sitenum,top.sitenum,mid.sitenum,botang,topang,midang)	# AJS
303 317
            
304 318
            # set the endpoint of the left and right Halfedge to be this vector
305 319
            if lbnd.edge.setEndpoint(lbnd.pm,v):
......
784 798
    """
785 799
    siteList = SiteList( points )
786 800
    context  = Context()
801
    context.use_bound = True	# AJS
787 802
    context.set_bounds( siteList )
788 803
    voronoi( siteList, context )
789
    return ( context.vertices, context.lines, context.edges, (context.xmin,context.ymin,context.xmax,context.ymax))
804
    return ( context.vertices, context.lines, context.edges, (context.xmin,context.ymin,context.xmax,context.ymax),context.poly_vertices) # AJS
790 805

  
791 806
#------------------------------------------------------------------
792 807
def computeDelaunayTriangulation( points ):
doGeometry.py (working copy)
457 457
    del writer
458 458
    return True
459 459
    
460
  def voronoi_polygons( self ):
461
    import voronoi
462
    
463
    ptDict = {} 
464
    ptNdx = -1 
465
    vprovider = self.vlayer.dataProvider()
466
    allAttrs = vprovider.attributeIndexes()
467
    vprovider.select( allAttrs )
468
    fields = {
469
    0 : QgsField( "PointId", QVariant.Int ),
470
    1 : QgsField( "PointName", QVariant.String )}
471
    writer = QgsVectorFileWriter( self.myName, self.myEncoding,
472
    fields, QGis.WKBPolygon, vprovider.crs() )
473
    inFeat = QgsFeature()
474
    points = []
475
    while vprovider.nextFeature( inFeat ):
476
      inGeom = QgsGeometry( inFeat.geometry() )
477
      point = inGeom.asPoint()
478
      points.append( point )
479
      ptNdx +=1  
480
      ptDict[ptNdx] = inFeat.id() 
481
    vprovider.rewind()
482
    vprovider.select( allAttrs )
483
    vertices,equations,edges,bounds,polygon_vertices = voronoi.computeVoronoiDiagram(points)
484
    feat = QgsFeature()
485
    nFeat = len( vertices )
486
    nElement = 0
487
    self.emit( SIGNAL( "runStatus(PyQt_PyObject)" ), 0 )
488
    self.emit( SIGNAL( "runRange(PyQt_PyObject)" ), ( 0, nFeat ) )
489

  
490
    polygon_vertices.sort()
491
    cur = -1
492
    for i,v in enumerate(polygon_vertices):
493
      if (not(v[0] == cur) and i > 0) or (i == len(polygon_vertices) - 1):
494
        polygon.append(curStart)
495
        feat.addAttribute( 0, QVariant( ptDict[ cur ] ) )
496
        geometry = QgsGeometry().fromPolygon( [ polygon ] )
497
        feat.setGeometry( geometry )      
498
        writer.addFeature( feat )
499
        nElement += 1
500
        self.emit( SIGNAL( "runStatus(PyQt_PyObject)" ),  nElement )
501
      if (not(v[0] == cur)) :
502
        cur = v[0]
503
        curStart = QgsPoint(v[2],v[3])
504
        polygon = []
505
      polygon.append( QgsPoint(v[2],v[3]) )
506

  
507
    del writer
508
    return True
509

  
460 510
  def delaunay_triangulation( self ):
461 511
    import voronoi
512
    
513
    ptDict = {}
514
    ptNdx = -1
462 515
    vprovider = self.vlayer.dataProvider()
463 516
    allAttrs = vprovider.attributeIndexes()
464 517
    vprovider.select( allAttrs )
......
474 527
      inGeom = QgsGeometry( inFeat.geometry() )
475 528
      point = inGeom.asPoint()
476 529
      points.append( point )
530
      ptNdx +=1
531
      ptDict[ptNdx] = inFeat.id()
477 532
    vprovider.rewind()
478 533
    vprovider.select( allAttrs )
479 534
    triangles = voronoi.computeDelaunayTriangulation( points )
......
488 543
      polygon = []
489 544
      step = 0
490 545
      for index in indicies:
491
        vprovider.featureAtId( index, inFeat, True,  allAttrs )
546
        vprovider.featureAtId( ptDict[index], inFeat, True,  allAttrs )
492 547
        geom = QgsGeometry( inFeat.geometry() )
493 548
        point = QgsPoint( geom.asPoint() )
494 549
        polygon.append( point )