Plots a filled voronoi diagram, using the given points and their colors.
  The first parameter must be an array-like or a scipy.stats.Voronoi object.

def voronoi_filled(points_or_voronoi, colors, show_points=False,
                   padding_fraction=0.05, cmap=None, ax=None, alpha=None,
                   edgecolor=None):
  '''Plots a filled voronoi diagram, using the given points and their colors.
  The first parameter must be an array-like or a scipy.stats.Voronoi object.
  '''
  from scipy.spatial import Voronoi  # Late import; scipy is optional

  # Disambiguate the first parameter
  if isinstance(points_or_voronoi, Voronoi):
    vor = points_or_voronoi
  else:
    points = np.asanyarray(points_or_voronoi)
    assert points.shape[1] == 2, 'Input points must be 2D'
    vor = Voronoi(points)

  # Borrowed from http://nbviewer.ipython.org/gist/pv/8037100
  regions = []
  vertices = vor.vertices.tolist()

  center = vor.points.mean(axis=0)
  radius = vor.points.ptp().max()*2

  # Construct a map containing all ridges for a given point
  all_ridges = {}
  for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
    all_ridges.setdefault(p1, []).append((p2, v1, v2))
    all_ridges.setdefault(p2, []).append((p1, v1, v2))

  # Reconstruct infinite regions
  for p1, region in enumerate(vor.point_region):
    verts = vor.regions[region]
    if all(v >= 0 for v in verts):
      # finite region
      regions.append(verts)
      continue

    # reconstruct a non-finite region
    ridges = all_ridges[p1]
    new_region = [v for v in verts if v >= 0]

    for p2, v1, v2 in ridges:
      if v2 < 0:
        v1, v2 = v2, v1
      if v1 >= 0:
        # finite ridge: already in the region
        continue

      # Compute the missing endpoint of an infinite ridge
      t = vor.points[p2] - vor.points[p1]  # tangent
      t /= np.linalg.norm(t)
      n = np.array([-t[1], t[0]])  # normal

      midpoint = vor.points[[p1, p2]].mean(axis=0)
      direction = np.sign((midpoint - center).dot(n)) * n
      far_point = vor.vertices[v2] + direction * radius

      new_region.append(len(vertices))
      vertices.append(far_point.tolist())

    # sort region counterclockwise
    vs = np.asarray([vertices[v] for v in new_region])
    vs -= vs.mean(axis=0)
    angle_order = np.argsort(np.arctan2(vs[:,1], vs[:,0]))
    new_region = np.array(new_region)[angle_order]

    # finish
    regions.append(new_region)
  vertices = np.asarray(vertices)

  # Plot colored polygons
  if ax is None:
    ax = plt.gca()
  polys = PatchCollection([Polygon(vertices[region]) for region in regions],
                          cmap=cmap, alpha=alpha, edgecolor=edgecolor)
  polys.set_array(np.asanyarray(colors))
  ax.add_collection(polys)

  if show_points:
    ax.plot(vor.points[:,0], vor.points[:,1], 'ko')

  # Zoom to a reasonable scale.
  pad = padding_fraction * (vor.max_bound - vor.min_bound)
  mins = vor.min_bound - pad
  maxes = vor.max_bound + pad
  ax.set_xlim(mins[0], maxes[0])
  ax.set_ylim(mins[1], maxes[1])
  return polys