Converting Voronoi Diagram Cell Areas to Lists of Pixel Coordinates

Question

I am using Voronoi diagrams for image processing (procedurally generated stippling). In order to do this I need to create a list (cells) of a list (coords_within_cell) of tuples (x,y pixel locations).

I have developed a couple brute-force algorithms to accomplish this (see below), but they are too slow to process more than ~10 points. The scipy spatial utilities seem to be more than 1000x more efficient. Because of this, I would like to use scipy to generate the Voronoi diagram: https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Voronoi.html

Using scipy to generate the Voronoi diagram is fairly simple but unfortunately I cannot figure out how to convert the cell areas into pixel coordinates. What is the best way to do this?

I found a related question, but it has no answers and it was deleted: https://web.archive.org/web/20200120151304/https://stackoverflow.com/questions/57703129/converting-a-voronoi-diagram-into-bitmap

Brute Force Algorithm 1 (too slow)

import math
import random
from PIL import Image 

def distance(x1, y1, x2, y2):
    return math.hypot(x2 - x1, y2 - y1)

# define the size of the x and y bounds
screen_width = 1260
screen_height = 1260

# define the number of points that should be used
number_of_points = 16

# randomly generate a list of n points within the given x and y bounds
point_x_coordinates = random.sample(range(0, screen_width), number_of_points)
point_y_coordinates = random.sample(range(0, screen_height), number_of_points)
points = list(zip(point_x_coordinates, point_y_coordinates))

# each point needs to have a corresponding list of pixels
point_pixels = []
for i in range(len(points)):
    point_pixels.append([]) 

# for each pixel within bounds, determine which point it is closest to and add it to the corresponding list in point_pixels
for pixel_y_coordinate in range(screen_height):
    for pixel_x_coordinate in  range(screen_width):
        distance_to_closest_point = float('inf')
        closest_point_index = 1

        for point_index, point in enumerate(points):
            distance_to_point = distance(pixel_x_coordinate, pixel_y_coordinate, point[0], point[1])
            if(distance_to_point < distance_to_closest_point):
                closest_point_index = point_index
                distance_to_closest_point = distance_to_point
        
        point_pixels[closest_point_index].append((pixel_x_coordinate, pixel_y_coordinate))

# each point needs to have a corresponding centroid
point_pixels_centroid = []

for pixel_group in point_pixels:
    x_sum = 0
    y_sum = 0
    for pixel in pixel_group:
        x_sum += pixel[0]
        y_sum += pixel[1]
    
    x_average = x_sum / len(pixel_group)
    y_average = y_sum / len(pixel_group)

    point_pixels_centroid.append((round(x_average), round(y_average)))



# display the resulting voronoi diagram
display_voronoi = Image.new("RGB", (screen_width, screen_height), "white")
for pixel_group in point_pixels:
    rgb = random.sample(range(0, 255), 3)
    for pixel in pixel_group:
        display_voronoi.putpixel( pixel, (rgb[0], rgb[1], rgb[2], 255) )

for centroid in point_pixels_centroid:
    print(centroid)
    display_voronoi.putpixel( centroid, (1, 1, 1, 255) )

display_voronoi.show()

Brute Force Algorithm 2 (also too slow): Based on this concept.

import math
import random
from PIL import Image 

def distance(x1, y1, x2, y2):
    return math.hypot(x2 - x1, y2 - y1)

# define the size of the x and y bounds
screen_width = 500
screen_height = 500

# define the number of points that should be used
number_of_points = 4

# randomly generate a list of n points within the given x and y bounds
point_x_coordinates = random.sample(range(0, screen_width), number_of_points)
point_y_coordinates = random.sample(range(0, screen_height), number_of_points)
points = list(zip(point_x_coordinates, point_y_coordinates))

# each point needs to have a corresponding list of pixels
point_pixels = []
for i in range(len(points)):
    point_pixels.append([]) 

# for each pixel within bounds, determine which point it is closest to and add it to the corresponding list in point_pixels
# do this by continuously growing circles outwards from the points
# if circles overlap then whoever was their first claims the location

# keep track of whether pixels have been used or not
# this is done via a 2D list of booleans
is_drawn_on = []
for i in range(screen_width):
    is_drawn_on.append([]) 
    for j in range(screen_height):
        is_drawn_on[i].append(False)

circles_are_growing = True
radius = 1
while(circles_are_growing):
    circles_are_growing = False
    for point_index, point in enumerate(points):
        for i in range(point[0] - radius, point[0] + radius):
            for j in range(point[1] - radius, point[1] + radius):
                # print(str(i)+" vs "+str(len(is_drawn_on)))
                if(i >= 0 and i < len(is_drawn_on)):
                    if(j >= 0 and j < len(is_drawn_on[i])):
                        if(not is_drawn_on[i][j] and distance(i, j, point[0], point[1]) <= radius):
                            point_pixels[point_index].append((i, j))
                            circles_are_growing = True
                            is_drawn_on[i][j] = True
    radius += 1

# each point needs to have a corresponding centroid
point_pixels_centroid = []

for pixel_group in point_pixels:
    x_sum = 0
    y_sum = 0
    for pixel in pixel_group:
        x_sum += pixel[0]
        y_sum += pixel[1]
    
    x_average = x_sum / len(pixel_group)
    y_average = y_sum / len(pixel_group)

    point_pixels_centroid.append((round(x_average), round(y_average)))



# display the resulting voronoi diagram
display_voronoi = Image.new("RGB", (screen_width, screen_height), "white")
for pixel_group in point_pixels:
    rgb = random.sample(range(0, 255), 3)
    for pixel in pixel_group:
        display_voronoi.putpixel( pixel, (rgb[0], rgb[1], rgb[2], 255) )

for centroid in point_pixels_centroid:
    print(centroid)
    display_voronoi.putpixel( centroid, (1, 1, 1, 255) )

display_voronoi.show()

Answer 1

Rather than build and interrogate the Voronoi diagram directly, it is easier to build and query a standard search tree. Below is my modification of your code using scipy.spatial.KDTree to determine the closest point for each pixel location followed by an image of the result (a 500x500 image with 500 Voronoi points).

The code is still a little slow but now scales well in the number of Voronoi points. This could be faster if you avoid building the list of pixel locations for each Voronoi cell and instead just directly set data in the image.

The fastest solution may involve building the Voronoi diagam and walking through it one pixel at a time and associating the closest Voronoi cell, looking at neighboring Voronoi cells when needed (since the previous pixel gives a very good guess for finding the Voronoi cell for the next pixel). But that will involve writing a lot more code that using the KDTree naively like this and probably not yield huge gains: the slow part of the code at this point is building all the per-pixel arrays/data which can be cleaned up independently.

import math
import random
from PIL import Image 
from scipy import spatial
import numpy as np

# define the size of the x and y bounds
screen_width = 500
screen_height = 500

# define the number of points that should be used
number_of_points = 500

# randomly generate a list of n points within the given x and y bounds
point_x_coordinates = random.sample(range(0, screen_width), number_of_points)
point_y_coordinates = random.sample(range(0, screen_height), number_of_points)
points = list(zip(point_x_coordinates, point_y_coordinates))

# each point needs to have a corresponding list of pixels
point_pixels = []
for i in range(len(points)):
    point_pixels.append([]) 

# build a search tree
tree = spatial.KDTree(points)

# build a list of pixed coordinates to query
pixel_coordinates = np.zeros((screen_height*screen_width, 2));
i = 0
for pixel_y_coordinate in range(screen_height):
    for pixel_x_coordinate in  range(screen_width):
        pixel_coordinates[i] = np.array([pixel_x_coordinate, pixel_y_coordinate])
        i = i+1

# for each pixel within bounds, determine which point it is closest to and add it to the corresponding list in point_pixels
[distances, indices] = tree.query(pixel_coordinates)

i = 0
for pixel_y_coordinate in range(screen_height):
    for pixel_x_coordinate in  range(screen_width):
        point_pixels[indices[i]].append((pixel_x_coordinate, pixel_y_coordinate))
        i = i+1

# each point needs to have a corresponding centroid
point_pixels_centroid = []

for pixel_group in point_pixels:
    x_sum = 0
    y_sum = 0
    for pixel in pixel_group:
        x_sum += pixel[0]
        y_sum += pixel[1]

    x_average = x_sum / max(len(pixel_group),1)
    y_average = y_sum / max(len(pixel_group),1)

    point_pixels_centroid.append((round(x_average), round(y_average)))


# display the resulting voronoi diagram
display_voronoi = Image.new("RGB", (screen_width, screen_height), "white")
for pixel_group in point_pixels:
    rgb = random.sample(range(0, 255), 3)
    for pixel in pixel_group:
        display_voronoi.putpixel( pixel, (rgb[0], rgb[1], rgb[2], 255) )

for centroid in point_pixels_centroid:
    #print(centroid)
    display_voronoi.putpixel( centroid, (1, 1, 1, 255) )

#display_voronoi.show()
display_voronoi.save("test.png")

Answer 2

scipy.interpolate.griddata does exactly that, and more or less by the same method as in Alex's answer

import numpy as np
from numpy.random import default_rng
from scipy.interpolate import griddata

# define the size of the x and y bounds
screen_width = 1260
screen_height = 1260

# define the number of points that should be used
number_of_points = 16

# randomly generate a list of n points within the given x and y bounds
rng = default_rng()
points = rng.random((number_of_points,2)) * [screen_width, screen_height]

grid_x, grid_y = np.mgrid[0:screen_width, 0:screen_height]

labels = griddata(points, np.arange(number_of_points), (grid_x, grid_y), method='nearest')

Then, you can use np.where(labels==10) to get the coordinates of all the pixels that belong to cell #10.

Or your can use all the machinery in scipy.ndimage to measure various properties of labeld regions. For instance the center of gravity.

If you want to display colorful cells:

from matplotlib.pyplot import imsave
rgb = rng.integers(0, 255, size=(number_of_points,3))
rgb_labels = rgb[labels]
imsave('test.png', rgb_labels)