title: Scipy: Voronoi diagram cells appear to spawn from nowhere
body: <p>I am researching <a href="https://en.wikipedia.org/wiki/Lloyd%27s_algorithm" rel="nofollow noreferrer">Lloyd iteration</a>, an iterative algorithm that distributes points in a space. In each iteration, Lloyd's algorithm builds a Voronoi map with each input point in its own Voronoi cell, then centers each point in its Voronoi cell.</p>
  
  <p>I'm seeing some strange behavior with Scipy's Voronoi implementation though: it seems that certain points are materializing out of nowhere in some iterations. The image below captures this behavior. If you look closely, you'll see two new points and cells appear in the center of the map a few iterations in:</p>
  
  <p><a href="https://i.stack.imgur.com/UNM9X.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/UNM9X.gif" alt="enter image description here"></a></p>
  
  <p>Here's the code used to generate the Voronoi distributions:</p>
  
  <pre><code>from scipy.spatial import Voronoi, voronoi_plot_2d
  import matplotlib.pyplot as plt
  import numpy as np
  import umap, os
  
  def find_centroid(verts):
    '''Return the centroid of a polygon described by `verts`'''
    area = 0
    x = 0
    y = 0
    for i in range(len(verts)-1):
      step = (verts[i, 0] * verts[i+1, 1]) - (verts[i+1, 0] * verts[i, 1])
      area += step
      x += (verts[i, 0] + verts[i+1, 0]) * step
      y += (verts[i, 1] + verts[i+1, 1]) * step
    if area == 0: area += 0.01
    return np.array([  (1/(3*area))*x,  (1/(3*area))*y  ])
  
  def lloyd_iterate(X):
    voronoi = Voronoi(X, qhull_options='Qbb Qc Qx')
    centroids = []
    for i in voronoi.regions:
      region = voronoi.vertices[i + [i[0]]]
      centroids.append( find_centroid( region ) )
    return np.array(centroids)
  
  def plot(X, name):
    '''Plot the Voronoi map of 2D numpy array X'''
    v = Voronoi(X, qhull_options='Qbb Qc Qx')
    plot = voronoi_plot_2d(v, show_vertices=False, line_colors='y', line_alpha=0.5, point_size=5)
    plot.set_figheight(14)
    plot.set_figwidth(20)
    plt.axis([-10, 10, -10, 10])
    if not os.path.exists('plots'): os.makedirs('plots')
    plot.savefig( 'plots/' + str(name) + '.png' )
  
  # get 1000 observations in two dimensions and plot their Voronoi map
  X = np.random.rand(1000, 4)
  X = umap.UMAP().fit_transform(X)
  plot(X, 0)
  
  # run several iterations, plotting each result
  for i in range(20):
    X = lloyd_iterate(X)
    plot(X, i)
  </code></pre>
  
  <p>Am I overlooking something, or is there something funny going on here? Any insight others can offer would be greatly appreciated. </p>