Question: I'm trying to create some artistic "plots" like the ones below:  The color of the regions do not really matter, what I'm trying to achieve is the variable "thickness" of the edges along the Voronoi regions (espescially, how they look like a bigger rounded blob where they meet in corners, and thinner at their middle point). I've tried by "painting manually" each pixel based on the minimum distance to each centroid (each associated with a color): n_centroids = 10 centroids = [(random.randint(0, h), random.randint(0, w)) for _ in range(n_centroids)] colors = np.array([np.random.choice(range(256), size=3) for _ in range(n_centroids)]) / 255 for x, y in it.product(range(h), range(w)): distances = np.sqrt([(x - c[0])**2 + (y - c[1])**2 for c in centroids]) centroid_i = np.argmin(distances) img[x, y] = colors[centroid_i] plt.imshow(img, cmap='gray')   Or by scipy.spatial.Voronoi, that also gives me the vertices points, although I still can't see how I can draw a line through them with the desired variable thickness. from scipy.spatial import Voronoi, voronoi_plot_2d # make up data points points = [(random.randint(0, 10), random.randint(0, 10)) for _ in range(10)] # add 4 distant dummy points points = np.append(points, [[999,999], [-999,999], [999,-999], [-999,-999]], axis = 0) # compute Voronoi tesselation vor = Voronoi(points) # plot voronoi_plot_2d(vor) # colorize for region in vor.regions: if not -1 in region: polygon = [vor.vertices[i] for i in region] plt.fill(*zip(*polygon)) # fix the range of axes plt.xlim([-2,12]), plt.ylim([-2,12]) plt.show()   Edit: I've managed to get a somewhat satisfying result via erosion + corner smoothing (via median filter as suggested in the comments) on each individual region, then drawing it into a black background. res = np.zeros((h,w,3)) for color in colors: region = (img == color)[:,:,0] region = region.astype(np.uint8) * 255 region = sg.medfilt2d(region, 15) # smooth corners # make edges from eroding regions region = cv2.erode(region, np.ones((3, 3), np.uint8)) region = region.astype(bool) res[region] = color plt.imshow(res)   But as you can see the "stretched" line along the boundaries/edges of the regions is not quite there. Any other suggestions?

Answer: This is what @JohanC suggestion looks like. IMO, it looks much better than my attempt with Bezier curves. However, there appears to be a small problem with the RoundedPolygon class, as there are sometimes small defects at the corners (e.g. between blue and purple in the image below). Edit: I fixed the RoundedPolygon class.  #!/usr/bin/env python # coding: utf-8 """ https://stackoverflow.com/questions/72061965/create-voronoi-art-with-rounded-region-edges """ import numpy as np import matplotlib.pyplot as plt from matplotlib import patches, path from scipy.spatial import Voronoi, voronoi_plot_2d def shrink(polygon, pad): center = np.mean(polygon, axis=0) resized = np.zeros_like(polygon) for ii, point in enumerate(polygon): vector = point - center unit_vector = vector / np.linalg.norm(vector) resized[ii] = point - pad * unit_vector return resized class RoundedPolygon(patches.PathPatch): # https://stackoverflow.com/a/66279687/2912349 def __init__(self, xy, pad, **kwargs): p = path.Path(*self.__round(xy=xy, pad=pad)) super().__init__(path=p, **kwargs) def __round(self, xy, pad): n = len(xy) for i in range(0, n): x0, x1, x2 = np.atleast_1d(xy[i - 1], xy[i], xy[(i + 1) % n]) d01, d12 = x1 - x0, x2 - x1 l01, l12 = np.linalg.norm(d01), np.linalg.norm(d12) u01, u12 = d01 / l01, d12 / l12 x00 = x0 + min(pad, 0.5 * l01) * u01 x01 = x1 - min(pad, 0.5 * l01) * u01 x10 = x1 + min(pad, 0.5 * l12) * u12 x11 = x2 - min(pad, 0.5 * l12) * u12 if i == 0: verts = [x00, x01, x1, x10] else: verts += [x01, x1, x10] codes = [path.Path.MOVETO] + n*[path.Path.LINETO, path.Path.CURVE3, path.Path.CURVE3] verts[0] = verts[-1] return np.atleast_1d(verts, codes) if __name__ == '__main__': # make up data points n = 100 max_x = 20 max_y = 10 points = np.c_[np.random.uniform(0, max_x, size=n), np.random.uniform(0, max_y, size=n)] # add 4 distant dummy points points = np.append(points, [[2 * max_x, 2 * max_y], [ -max_x, 2 * max_y], [2 * max_x, -max_y], [ -max_x, -max_y]], axis = 0) # compute Voronoi tesselation vor = Voronoi(points) fig, ax = plt.subplots(figsize=(max_x, max_y)) for region in vor.regions: if region and (not -1 in region): polygon = np.array([vor.vertices[i] for i in region]) resized = shrink(polygon, 0.15) ax.add_patch(RoundedPolygon(resized, 0.2, color=plt.cm.Reds(0.5 + 0.5*np.random.rand()))) ax.axis([0, max_x, 0, max_y]) ax.axis('off') ax.set_facecolor('black') ax.add_artist(ax.patch) ax.patch.set_zorder(-1) plt.show()