I am trying to implement an image stippling algorithm in python, and want to vectorize calculating the density (average luminance) of labelled image regions (Voronoi cells). Currently I'm able to do so using a loop, but this is too computationally intensive for large numbers of regions. How can I vectorize this operation?
import numpy as np
from skimage import io
from scipy.interpolate import griddata
number_of_points = 1000
img = io.imread('https://www.kindpng.com/picc/m/111-1114964_house-icon-png-old-house-easy-drawing-transparent.png', as_gray=True)
height, width = img.shape
# generate random points
rng = np.random.default_rng()
points = rng.random((number_of_points,2)) * [width, height]
# calculate labelled regions
grid_x, grid_y = np.mgrid[0:width, 0:height]
labels = griddata(points, np.arange(number_of_points), (grid_x, grid_y), method='nearest')
# calculate density per region (mean of grayscale values of pixels in each region)
point_idxs = np.arange(len(points))
density = [np.mean(img[labels.T==i]) for i in point_idxs] # <-- this is the bottleneck