def processAlgorithm(self, progress):
        layer = dataobjects.getObjectFromUri(self.getParameterValue(self.INPUT))

        buf = self.getParameterValue(self.BUFFER)

        writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(
            layer.pendingFields().toList(), QGis.WKBPolygon, layer.crs())

        inFeat = QgsFeature()
        outFeat = QgsFeature()
        extent = layer.extent()
        extraX = extent.height() * (buf / 100.0)
        extraY = extent.width() * (buf / 100.0)
        height = extent.height()
        width = extent.width()
        c = voronoi.Context()
        pts = []
        ptDict = {}
        ptNdx = -1

        features = vector.features(layer)
        for inFeat in features:
            geom = QgsGeometry(inFeat.geometry())
            point = geom.asPoint()
            x = point.x() - extent.xMinimum()
            y = point.y() - extent.yMinimum()
            pts.append((x, y))
            ptNdx += 1
            ptDict[ptNdx] = inFeat.id()

        if len(pts) < 3:
            raise GeoAlgorithmExecutionException(
                self.tr('Input file should contain at least 3 points. Choose '
                        'another file and try again.'))

        uniqueSet = Set(item for item in pts)
        ids = [pts.index(item) for item in uniqueSet]
        sl = voronoi.SiteList([voronoi.Site(i[0], i[1], sitenum=j) for (j,
                              i) in enumerate(uniqueSet)])
        voronoi.voronoi(sl, c)
        inFeat = QgsFeature()

        current = 0
        total = 100.0 / float(len(c.polygons))

        for (site, edges) in c.polygons.iteritems():
            request = QgsFeatureRequest().setFilterFid(ptDict[ids[site]])
            inFeat = layer.getFeatures(request).next()
            lines = self.clip_voronoi(edges, c, width, height, extent, extraX, extraY)

            geom = QgsGeometry.fromMultiPoint(lines)
            geom = QgsGeometry(geom.convexHull())
            outFeat.setGeometry(geom)
            outFeat.setAttributes(inFeat.attributes())
            writer.addFeature(outFeat)

            current += 1
            progress.setPercentage(int(current * total))

        del writer

  def delaunay_triangulation( self ):
    import voronoi
    from sets import Set
    vprovider = self.vlayer.dataProvider()

    fields = QgsFields()
    fields.append( QgsField( "POINTA", QVariant.Double ) )
    fields.append( QgsField( "POINTB", QVariant.Double ) )
    fields.append( QgsField( "POINTC", QVariant.Double ) )

    writer = QgsVectorFileWriter( self.myName, self.myEncoding, fields,
                                  QGis.WKBPolygon, vprovider.crs() )
    inFeat = QgsFeature()
    c = voronoi.Context()
    pts = []
    ptDict = {}
    ptNdx = -1
    fit = vprovider.getFeatures()
    while fit.nextFeature( inFeat ):
      geom = QgsGeometry( inFeat.geometry() )
      point = geom.asPoint()
      x = point.x()
      y = point.y()
      pts.append( ( x, y ) )
      ptNdx +=1
      ptDict[ptNdx] = inFeat.id()
    if len(pts) < 3:
      return False
    uniqueSet = Set( item for item in pts )
    ids = [ pts.index( item ) for item in uniqueSet ]
    sl = voronoi.SiteList( [ voronoi.Site( *i ) for i in uniqueSet ] )
    c.triangulate = True
    voronoi.voronoi( sl, c )
    triangles = c.triangles
    feat = QgsFeature()
    nFeat = len( triangles )
    nElement = 0
    self.emit( SIGNAL( "runStatus( PyQt_PyObject )" ), 0 )
    self.emit( SIGNAL( "runRange( PyQt_PyObject )" ), ( 0, nFeat ) )
    for triangle in triangles:
      indicies = list( triangle )
      indicies.append( indicies[ 0 ] )
      polygon = []
      step = 0
      for index in indicies:
        vprovider.getFeatures( QgsFeatureRequest().setFilterFid( ptDict[ ids[ index ] ] ) ).nextFeature( inFeat )
        geom = QgsGeometry( inFeat.geometry() )
        point = QgsPoint( geom.asPoint() )
        polygon.append( point )
        if step <= 3: feat.setAttribute( step, QVariant( ids[ index ] ) )
        step += 1
      geometry = QgsGeometry().fromPolygon( [ polygon ] )
      feat.setGeometry( geometry )
      writer.addFeature( feat )
      nElement += 1
      self.emit( SIGNAL( "runStatus( PyQt_PyObject )" ), nElement )
    del writer
    return True

 def processAlgorithm(self, progress):
     vlayer = QGisLayers.getObjectFromUri(self.getParameterValue(Delaunay.INPUT))
     vprovider = vlayer.dataProvider()
     allAttrs = vprovider.attributeIndexes()
     vprovider.select( allAttrs )
     fields = {
               0 : QgsField( "POINTA", QVariant.Double ),
               1 : QgsField( "POINTB", QVariant.Double ),
               2 : QgsField( "POINTC", QVariant.Double ) }
     writer = self.getOutputFromName(Delaunay.OUTPUT).getVectorWriter(fields, QGis.WKBPolygon, vprovider.crs() )
     inFeat = QgsFeature()
     c = voronoi.Context()
     pts = []
     ptDict = {}
     ptNdx = -1
     while vprovider.nextFeature(inFeat):
       geom = QgsGeometry(inFeat.geometry())
       point = geom.asPoint()
       x = point.x()
       y = point.y()
       pts.append((x, y))
       ptNdx +=1
       ptDict[ptNdx] = inFeat.id()
     if len(pts) < 3:
       return False
     uniqueSet = Set(item for item in pts)
     ids = [pts.index(item) for item in uniqueSet]
     sl = voronoi.SiteList([voronoi.Site(*i) for i in uniqueSet])
     c.triangulate = True
     voronoi.voronoi(sl, c)
     triangles = c.triangles
     feat = QgsFeature()
     nFeat = len( triangles )
     nElement = 0
     for triangle in triangles:
       indicies = list(triangle)
       indicies.append(indicies[0])
       polygon = []
       step = 0
       for index in indicies:
         vprovider.featureAtId(ptDict[ids[index]], inFeat, True,  allAttrs)
         geom = QgsGeometry(inFeat.geometry())
         point = QgsPoint(geom.asPoint())
         polygon.append(point)
         if step <= 3: feat.addAttribute(step, QVariant(ids[index]))
         step += 1
       geometry = QgsGeometry().fromPolygon([polygon])
       feat.setGeometry(geometry)
       writer.addFeature(feat)
       nElement += 1
       progress.setPercentage(nElement/nFeat * 100)
     del writer

 def voronoi_polygons( self ):
   vprovider = self.vlayer.dataProvider()
   allAttrs = vprovider.attributeIndexes()
   vprovider.select( allAttrs )
   writer = QgsVectorFileWriter( self.myName, self.myEncoding, vprovider.fields(),
                                 QGis.WKBPolygon, vprovider.crs() )
   inFeat = QgsFeature()
   outFeat = QgsFeature()
   extent = self.vlayer.extent()
   extraX = extent.height() * ( self.myParam / 100.00 )
   extraY = extent.width() * ( self.myParam / 100.00 )
   height = extent.height()
   width = extent.width()
   c = voronoi.Context()
   pts = []
   ptDict = {}
   ptNdx = -1
   while vprovider.nextFeature( inFeat ):
     geom = QgsGeometry( inFeat.geometry() )
     point = geom.asPoint()
     x = point.x() - extent.xMinimum()
     y = point.y() - extent.yMinimum()
     pts.append( ( x, y ) )
     ptNdx +=1
     ptDict[ ptNdx ] = inFeat.id()
   self.vlayer = None
   if len( pts ) < 3:
     return False
   uniqueSet = Set( item for item in pts )
   ids = [ pts.index( item ) for item in uniqueSet ]
   sl = voronoi.SiteList( [ voronoi.Site( i[ 0 ], i[ 1 ], sitenum = j ) for j, i in enumerate( uniqueSet ) ] )
   voronoi.voronoi( sl, c )
   inFeat = QgsFeature()
   nFeat = len( c.polygons )
   nElement = 0
   self.emit( SIGNAL( "runStatus( PyQt_PyObject )" ), 0 )
   self.emit( SIGNAL( "runRange( PyQt_PyObject )" ), ( 0, nFeat ) )
   for site, edges in c.polygons.iteritems():
     vprovider.featureAtId( ptDict[ ids[ site ] ], inFeat, True,  allAttrs )
     lines = self.clip_voronoi( edges, c, width, height, extent, extraX, extraY )
     geom = QgsGeometry.fromMultiPoint( lines )
     geom = QgsGeometry( geom.convexHull() )
     outFeat.setGeometry( geom )
     outFeat.setAttributeMap( inFeat.attributeMap() )
     writer.addFeature( outFeat )
     nElement += 1
     self.emit( SIGNAL( "runStatus( PyQt_PyObject )" ), nElement )
   del writer
   return True

 def display(network):
     """Display with matplotlib"""
     import matplotlib
     import matplotlib.pyplot as plt
     try:
         from voronoi import voronoi
     except:
         voronoi = None        
     fig = plt.figure(figsize=(10,10))
     axes = fig.add_subplot(1,1,1)
     # Draw samples
     x, y = samples[:,0], samples[:,1]
     plt.scatter(x, y, s=1.0, color='b', alpha=0.1, zorder=1)
     # Draw network
     x, y = network.nodes[...,0], network.nodes[...,1]
     if len(network.nodes.shape) > 2:
         for i in range(network.nodes.shape[0]):
             plt.plot (x[i,:], y[i,:], 'k', alpha=0.85, lw=1.5, zorder=2)
         for i in range(network.nodes.shape[1]):
             plt.plot (x[:,i], y[:,i], 'k', alpha=0.85, lw=1.5, zorder=2)
     else:
         plt.plot (x, y, 'k', alpha=0.85, lw=1.5, zorder=2)
     plt.scatter (x, y, s=50, c='w', edgecolors='k', zorder=3)
     if voronoi is not None:
         segments = voronoi(x.ravel(), y.ravel())
         lines = matplotlib.collections.LineCollection(segments, color='0.65')
         axes.add_collection(lines)
     plt.axis([0,1,0,1])
     plt.xticks([]), plt.yticks([])
     plt.show()

 def processAlgorithm(self, progress):
     vlayer = QGisLayers.getObjectFromUri(self.getParameterValue(VoronoiPolygons.INPUT))
     vprovider = vlayer.dataProvider()
     allAttrs = vprovider.attributeIndexes()
     vprovider.select( allAttrs )
     writer = self.getOutputFromName(VoronoiPolygons.OUTPUT).getVectorWriter(vprovider.fields(), QGis.WKBPolygon, vprovider.crs() )
     inFeat = QgsFeature()
     outFeat = QgsFeature()
     extent = vlayer.extent()
     height = extent.height()
     width = extent.width()
     c = voronoi.Context()
     pts = []
     ptDict = {}
     ptNdx = -1
     while vprovider.nextFeature(inFeat):
       geom = QgsGeometry(inFeat.geometry())
       point = geom.asPoint()
       x = point.x()-extent.xMinimum()
       y = point.y()-extent.yMinimum()
       pts.append((x, y))
       ptNdx +=1
       ptDict[ptNdx] = inFeat.id()
     #self.vlayer = None
     if len(pts) < 3:
       return False
     uniqueSet = Set(item for item in pts)
     ids = [pts.index(item) for item in uniqueSet]
     sl = voronoi.SiteList([voronoi.Site(i[0], i[1], sitenum=j) for j, i in enumerate(uniqueSet)])
     voronoi.voronoi(sl, c)
     inFeat = QgsFeature()
     nFeat = len(c.polygons)
     nElement = 0
     for site, edges in c.polygons.iteritems():
       vprovider.featureAtId(ptDict[ids[site]], inFeat, True,  allAttrs)
       lines = self.clip_voronoi(edges, c, width, height, extent, 0, 0)
       geom = QgsGeometry.fromMultiPoint(lines)
       geom = QgsGeometry(geom.convexHull())
       outFeat.setGeometry(geom)
       outFeat.setAttributeMap(inFeat.attributeMap())
       writer.addFeature(outFeat)
       nElement += 1
       progress.setPercentage(nElement/nFeat * 100)
     del writer
     return True

def delaunay(X):
    """Return a Delaunay triangulation of the specified Coords.

    While the Coords are 3d, only the first 2 components are used.

    Returns a TriSurface with the Delaunay trinagulation in the x-y plane.
    """
    from voronoi import voronoi
    return TriSurface(X,voronoi(X[:,:2]).triangles)

    def computeWalls(self):
        nodes = self.graph.nodes()

        # Poached from voronoi.computeVoronoiDiagram
        # http://stackoverflow.com/questions/9441007/how-can-i-get-a-dictionary-of-cells-from-this-voronoi-diagram-data
        site_list = SiteList(nodes)
        context = Context()
        voronoi(site_list, context)
        verts = context.vertices

        walls = [None for _ in context.edges]

        for i, v1, v2 in context.edges:
            path = context.bisectors[i]
            if all(v != -1 and not self.pointOutsideBounds(*verts[v])
                   for v in (v1, v2)):
                # edge has a vertex at either end, easy
                walls[i] = WallCrossing((Coords(*verts[v1]), Coords(*verts[v2])), path)
                continue
            if self.pointOutsideBounds(*verts[v1]):
                v1 = -1
            elif self.pointOutsideBounds(*verts[v2]):
                v2 = -1
            # apparently v1, v2 go left to right
            # calculate an edge point using slope = -a/b
            a, b, _ = context.lines[i]
            p0 = Coords(*verts[v1 if v1 != -1 else v2])
            if self.pointOutsideBounds(*p0):
                continue
            # need to handle case where b is 0
            p1 = Coords(*constrain(p0, (-a / b) if b else (-a * float('inf')),
                        self.dims, rightward=(v1 != -1)))
            walls[i] = WallCrossing((p0, p1), path)

        wall_dict = dict()
        for site, edge_list in context.polygons.items():
            wall_dict[nodes[site]] = [walls[i].wall for i, _, _ in edge_list
                                      if walls[i] is not None]

        return [wall for wall in walls if wall is not None], wall_dict

    def compute(self,points):
        import voronoi
        self.mcontext  = voronoi.voronoi(points)

        self.medges    = self.buildtriforedges()

        self.mvertices = self.buildtriforvertices(self.medges)

        self.mexts     = self.buildextseqsforvertices(self.mvertices)

       # centers = buildcentersfortri(context)

        self.mpolygons = self.buildpolygonsfromexts(self.mpoints,self.mexts)

def VoronoiLineEdges(PointsMap):
  Sitepts = []
  pts = {}

  #print CurrentDate, PointsMap[PointsMap.keys()[0]]
  for grid, stn in PointsMap.items():

    x=float(stn[0])
    y=float(stn[1])
    station=grid
    #station.extend( stn[3:])
    #print x,y,station

    pts[ (x,y) ]=station
  
  stncounts=len(pts.keys())
  #print stncounts, "points"
  
  site_points=[]
  for pt in pts.keys():
    Sitepts.append(voronoi.Site(pt[0],pt[1]))
    site_points.append( (pt[0],pt[1]) )
      

  #print "Calculating Voronoi Lattice",
  
  siteList = voronoi.SiteList(Sitepts)
  context  = voronoi.Context()
  voronoi.Edge.EDGE_NUM=0  
  voronoi.voronoi(siteList,context)

  vertices=context.vertices
  lines=context.lines
  edges=context.edges
  triangles=context.triangles
  has_edge=context.has_edge

  return vertices, lines, edges, has_edge

    def learn(network, samples, epochs=25000, sigma=(10, 0.01), lrate=(0.5,0.005)):
        network.learn(samples, epochs)

        fig = plt.figure(figsize=(10,10))
        axes = fig.add_subplot(1,1,1)
        # Draw samples
        x,y = samples[:,0], samples[:,1]
        plt.scatter(x, y, s=1.0, color='b', alpha=0.1, zorder=1)
        # Draw network
        x,y = network.codebook[...,0], network.codebook[...,1]
        plt.scatter (x, y, s=50, c='w', edgecolors='k', zorder=3)
        if voronoi is not None:
            segments = voronoi(x.ravel(),y.ravel())
            lines = matplotlib.collections.LineCollection(segments, color='0.65')
            axes.add_collection(lines)
        plt.axis([0,1,0,1])
        plt.xticks([]), plt.yticks([])
        plt.show()

def voronoi_polygons(points):
    """Return series of polygons for given coordinates. """
    c = voronoi(points)
    # For each point find triangles (vertices) of a cell
    point_in_triangles = defaultdict(set)
    for t_ind, ps in enumerate(c.triangles):
        for p in ps:
            point_in_triangles[p].add(t_ind)
            
    # Vertex connectivity graph
    vertex_graph = defaultdict(set)
    for e_ind, (_, r, l) in enumerate(c.edges):
        vertex_graph[r].add(l)
        vertex_graph[l].add(r)
    
    # Calculate cells
    def cell(point):
        if point not in point_in_triangles:
            return None
        vertices = set(point_in_triangles[point]) # copy
        v_cell = [vertices.pop()]
#        vertices.add(-1)  # Simulate infinity :-)
        while vertices:
            neighbours = vertex_graph[v_cell[-1]] & vertices
            if not neighbours:
                break
            v_cell.append(neighbours.pop())
            vertices.discard(v_cell[-1])
        return v_cell
        
    # Finally, produces polygons
    polygons = []
    for i, p in enumerate(points):
        vertices = []
        point_cell = cell(i)
        for j in point_cell:
            vertices.append(c.vertices[j])
        polygons.append(tuple(vertices))
    return tuple(polygons)

def quantize(data, cutoff, numRegions=8, tolerance=0.0001):

    newData = data[::4, ::4]  # Downsample data to speed up by a factor of 16
    # Next remove clouds poorly represented data
    cleanData = newData[np.logical_and(newData > 0, newData < cutoff)]

    tess = voronoi(cleanData, numRegions, tolerance, data)

    thresh = []
    for val in tess.values():
        thresh.append(np.min(val))
    thresh = np.sort(thresh)

    (rows, cols) = data.shape
    qData = np.ndarray((rows, cols), dtype='uint8')
    qData.fill(255)

    for val, t in enumerate(thresh):
        qData[data > t] = val + 1

    qData[data == 0] = 0  # Put back the land values
    qData[data == -1] = 255  # Put back the cloud values

    return qData

    def effect(self):
        if not self.options.ids:
            inkex.errormsg(_("Please select an object"))
            exit()
        scale = self.unittouu('1px')            # convert to document units
        self.options.size *= scale
        self.options.border *= scale
        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] = scale*float(p.stdout.read())
        mat = simpletransform.composeParents(self.selected[self.options.ids[0]], [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
        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'] - mat[0][2], q['y'] - mat[1][2]))
        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)
        path = ""
        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 %.3f,%.3f %.3f,%.3f ' % (x1, y1, x2, y2)

        patternstyle = {'stroke': '#000000', 'stroke-width': str(scale)}
        attribs = {'d': path, 'style': simplestyle.formatStyle(patternstyle)}
        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)
#.........这里部分代码省略.........

    def processAlgorithm(self, progress):
        layer = dataobjects.getObjectFromUri(
            self.getParameterValue(self.INPUT))

        fields = [QgsField('POINTA', QVariant.Double, '', 24, 15),
                  QgsField('POINTB', QVariant.Double, '', 24, 15),
                  QgsField('POINTC', QVariant.Double, '', 24, 15)]

        writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(fields,
                                                                     QGis.WKBPolygon, layer.crs())

        pts = []
        ptDict = {}
        ptNdx = -1
        c = voronoi.Context()
        features = vector.features(layer)
        total = 100.0 / len(features)
        for current, inFeat in enumerate(features):
            geom = QgsGeometry(inFeat.geometry())
            point = geom.asPoint()
            x = point.x()
            y = point.y()
            pts.append((x, y))
            ptNdx += 1
            ptDict[ptNdx] = inFeat.id()
            progress.setPercentage(int(current * total))

        if len(pts) < 3:
            raise GeoAlgorithmExecutionException(
                self.tr('Input file should contain at least 3 points. Choose '
                        'another file and try again.'))

        uniqueSet = set(item for item in pts)
        ids = [pts.index(item) for item in uniqueSet]
        sl = voronoi.SiteList([voronoi.Site(*i) for i in uniqueSet])
        c.triangulate = True
        voronoi.voronoi(sl, c)
        triangles = c.triangles
        feat = QgsFeature()

        total = 100.0 / len(triangles)
        for current, triangle in enumerate(triangles):
            indicies = list(triangle)
            indicies.append(indicies[0])
            polygon = []
            attrs = []
            step = 0
            for index in indicies:
                request = QgsFeatureRequest().setFilterFid(ptDict[ids[index]])
                inFeat = layer.getFeatures(request).next()
                geom = QgsGeometry(inFeat.geometry())
                point = QgsPoint(geom.asPoint())
                polygon.append(point)
                if step <= 3:
                    attrs.append(ids[index])
                step += 1
            feat.setAttributes(attrs)
            geometry = QgsGeometry().fromPolygon([polygon])
            feat.setGeometry(geometry)
            writer.addFeature(feat)
            progress.setPercentage(int(current * total))

        del writer

    # Draw contour map of potential
    im = grid[gk].contourf(X,Y,V,clines,cmap=my_cmap)

    # Draw circles defining state boundaries
    cirA = plt.Circle((-3,0), radius=1.0,  ec='k', fill=False, lw=0.5)
    cirB = plt.Circle((3,0), radius=1.0,  ec='k', fill=False, lw=0.5)
    grid[gk].add_patch(cirA)
    grid[gk].add_patch(cirB)

    if k == 0:
        cen = centers[:ncenters,:]
    else:
        cen = centers[ncenters:,:]

    grid[gk].plot(cen[:,0],cen[:,1],'k.-',lw=1,ms=2,marker='o')

    segments = voronoi.voronoi(cen[:,0],cen[:,1])
    lines = mpl.collections.LineCollection(segments, color='k',lw=0.5)
    grid[gk].add_collection(lines)

    grid[gk].axis([-5,5,-5,5])

grid[2].cax.colorbar(im)
grid[2].cax.toggle_label(True)

    
#plt.show()
plt.savefig('potential_string.eps',dpi=600,format='eps',bbox_inches='tight')

    num_points = min(ii.size,500)
    jj = ii[:num_points]

    if k % 2 == 0:
        cc = '0.8'
    else:
        cc = '0.4'

    grid[0].plot(crd[jj,0],crd[jj,1],c=cc,mec=cc,mfc=cc,marker='.',ms=3,ls='None',zorder=1)


# Plot string
grid[0].plot(centers[:,0],centers[:,1],ms=3,marker='o',color='k',ls='None',zorder=2)

# Plot Voronoi Cells
segments = voronoi.voronoi(cen_rep[:,0],cen_rep[:,1])
lines = mpl.collections.LineCollection(segments, color='k',lw=1.0,zorder=3)
grid[0].add_collection(lines)

grid[0].set_xlabel('X')
grid[0].set_ylabel('Y')


grid[0].axis([-1,1,0,1])
grid[0].cax.colorbar(im)
grid[0].cax.toggle_label(True)

fig.set_size_inches(fsize)
plt.tight_layout()
plt.savefig('potential_string.eps',dpi=600,format='eps',bbox_inches='tight')
plt.show()

def generate_map(): 
 

	screen.fill((0, 0, 0))

	points = generate_random_points(num_points, width, height, buf)


	#for x, y in points:
	#	pygame.draw.circle(screen, WHITE, (x,y), 2, 1)

	voronoi_context = voronoi(points)

	voronoi_point_dict = {}
	point_to_segment_dict = {}
	segments = []
	vertices = []

	top_l =  Point(0,0)
	top_r = Point(width,0)
	bottom_l = Point(0, height)
	bottom_r = Point(width, height) 

	top = Line(top_l, top_r) 
	left = Line(top_l, bottom_l) 
	right = Line(top_r, bottom_r) 
	bottom = Line(bottom_l, bottom_r) 

	boundaries = [top, right, bottom, left]

	for edge in voronoi_context.edges:
		il, i1, i2 = edge # index of line, index of vertex 1, index of vertex 2

		line_color = RED 

		vert1 = None
		vert2 = None
		print_line = True

		if i1 is not -1 and i2 is not -1:
			vert1 = voronoi_context.vertices[i1]
			vert2 = voronoi_context.vertices[i2]

		else:
			line_point = None

			if i1 is -1:
				line_p = voronoi_context.vertices[i2]
			if i2 is -1: 
				line_p = voronoi_context.vertices[i1]

			line_point = Point(line_p[0], line_p[1])
			line = voronoi_context.lines[il] 

			p1 = None
			p2 = None
			if line[1] == 0:
				p1 = line_point
				p2 = Point(line[0]/line[2], 1)
			else: 
				p1 = Point(0, line[2]/line[1])
				p2 = line_point

			l = Line(p1, p2)

			top_intersect = l.intersection(top)
			bottom_intersect = l.intersection(bottom)
			right_intersect = l.intersection(right)
			left_intersect = l.intersection(left)

			distances = []

			top_dist = None
			bottom_dist = None
			right_dist = None
			left_dist = None

			if len(top_intersect) != 0: 
				top_dist = abs(line_point.distance(top_intersect[0]))
				distances.append(top_dist)
			if len(bottom_intersect) != 0 : 
				bottom_dist = abs(line_point.distance(bottom_intersect[0]))
				distances.append(bottom_dist)
			if len(right_intersect) != 0:
				right_dist = abs(line_point.distance(right_intersect[0]))
				distances.append(right_dist)
			if len(left_intersect) != 0: 
				left_dist = abs(line_point.distance(left_intersect[0]))
				distances.append(left_dist)

			vert1 = line_p 
			v2 = None

			if top_dist == min(distances):
				v2 = top_intersect[0]
			elif bottom_dist == min(distances):
				v2 = bottom_intersect[0]
			elif right_dist == min(distances):
				v2 = right_intersect[0]
			elif left_dist == min(distances):
#.........这里部分代码省略.........

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import matplotlib
import numpy as np
from voronoi import voronoi
import matplotlib.pyplot as plt

X, Y = np.meshgrid(np.linspace(-0.1, 1.1, 25), np.linspace(-0.1, 1.1, 25))
X = X.ravel() + np.random.uniform(-0.025, 0.025, X.size)
Y = Y.ravel() + np.random.uniform(-0.025, 0.025, Y.size)
cells, triangles, circles = voronoi(X, Y)

fig = plt.figure(figsize=(8, 6))
axes = plt.subplot(111, aspect=1)
for cell in cells:
    codes = (
        [matplotlib.path.Path.MOVETO]
        + [matplotlib.path.Path.LINETO] * (len(cell) - 2)
        + [matplotlib.path.Path.CLOSEPOLY]
    )
    path = matplotlib.path.Path(cell, codes)
    color = np.random.uniform(0.4, 0.9, 3)
    patch = matplotlib.patches.PathPatch(path, facecolor=color, edgecolor="w", zorder=-1)
    axes.add_patch(patch)

plt.axis([0, 1, 0, 1])
plt.xticks([]), plt.yticks([])
plt.show()

def VoronoiPolygons(PointsMap, BoundingBox="W", PlotMap=False):
  global PlotIt
  global WorldRange
  global WorldRanges

  if type(BoundingBox)==type([]):
    if len(BoundingBox)==4:
      WorldRange=BoundingBox
    else:
      return "Error in Bounding Box"
  else:
    WorldRange=WorldRanges[BoundingBox]
    
  PlotIt=PlotMap

  currenttime=time.time()
  Sitepts = []
  pts = {}

  #print CurrentDate, PointsMap[PointsMap.keys()[0]]
  for grid, stn in PointsMap.items():

    x=float(stn[0])
    y=float(stn[1])
    station=grid
    #station.extend( stn[3:])
    #print x,y,station

    pts[ (x,y) ]=station
  
  stncounts=len(pts.keys())
  #print stncounts, "points"
  
  site_points=[]
  for pt in pts.keys():
    Sitepts.append(voronoi.Site(pt[0],pt[1]))
    site_points.append( (pt[0],pt[1]) )
      

  #print "Calculating Voronoi Lattice",
  
  siteList = voronoi.SiteList(Sitepts)
  context  = voronoi.Context()
  voronoi.Edge.EDGE_NUM=0   
  voronoi.voronoi(siteList,context)

  vertices=context.vertices
  lines=context.lines
  edges=context.edges
  
  #print edges

  #For Faster Access
  edge_dic={}
  for edge in edges:
    edge_dic[edge[0]]=edge[1:]
  
  triangles=context.triangles
  has_edge=context.has_edge

  voronoi_lattice={}

  m_range={}
  m_range["max_x"]=-9999999999
  m_range["max_y"]=-9999999999
  m_range["min_x"]=9999999999
  m_range["min_y"]=9999999999

  #Get the range!!
  for pnt in site_points:
    m_range=update_maxmin(m_range, pnt[0], pnt[1])

  #print "Getting the Polygons"

  prev_percent=0
  for station, ls in has_edge.items():


    voronoi_lattice[station]={}
    voronoi_lattice[station]["coordinate"]=site_points[station]
    voronoi_lattice[station]["info"]=pts[ site_points[station] ]

    polygon=[]
 
    prev_extreme=[]
    Verbose=True
    if Verbose: 
      current_percent=int(station/float(stncounts)*100)
      if current_percent!=prev_percent:
        #print station,"/", stncounts, current_percent, "% Done" 
        timeelapse=time.time()-currenttime
        #print station, timeelapse
        currenttime=time.time()

      prev_percent=current_percent

    #For every lines that the station owns
    for l in ls:
      e=edge_dic[l]

#.........这里部分代码省略.........