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?

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()