def voronoi_edges(shape: List[int], radius: int, ncells: int,
                  flat_faces: bool = True):
    r"""
    Create an image of the edges in a Voronoi tessellation

    Parameters
    ----------
    shape : array_like
        The size of the image to generate in [Nx, Ny, Nz] where Ni is the
        number of voxels in each direction.

    radius : scalar
        The radius to which Voronoi edges should be dilated in the final image.

    ncells : scalar
        The number of Voronoi cells to include in the tesselation.

    flat_faces : Boolean
        Whether the Voronoi edges should lie on the boundary of the
        image (True), or if edges outside the image should be removed (False).

    Returns
    -------
    image : ND-array
        A boolean array with ``True`` values denoting the pore space

    """
    print(78*'―')
    print('voronoi_edges: Generating', ncells, 'cells')
    shape = sp.array(shape)
    if sp.size(shape) == 1:
        shape = sp.full((3, ), int(shape))
    im = sp.zeros(shape, dtype=bool)
    base_pts = sp.rand(ncells, 3)*shape
    if flat_faces:
        # Reflect base points
        Nx, Ny, Nz = shape
        orig_pts = base_pts
        base_pts = sp.vstack((base_pts,
                              [-1, 1, 1] * orig_pts + [2.0*Nx, 0, 0]))
        base_pts = sp.vstack((base_pts,
                              [1, -1, 1] * orig_pts + [0, 2.0*Ny, 0]))
        base_pts = sp.vstack((base_pts,
                              [1, 1, -1] * orig_pts + [0, 0, 2.0*Nz]))
        base_pts = sp.vstack((base_pts, [-1, 1, 1] * orig_pts))
        base_pts = sp.vstack((base_pts, [1, -1, 1] * orig_pts))
        base_pts = sp.vstack((base_pts, [1, 1, -1] * orig_pts))
    vor = sptl.Voronoi(points=base_pts)
    vor.vertices = sp.around(vor.vertices)
    vor.vertices *= (sp.array(im.shape)-1) / sp.array(im.shape)
    vor.edges = _get_Voronoi_edges(vor)
    for row in vor.edges:
        pts = vor.vertices[row].astype(int)
        if sp.all(pts >= 0) and sp.all(pts < im.shape):
            line_pts = line_segment(pts[0], pts[1])
            im[tuple(line_pts)] = True
    im = spim.distance_transform_edt(~im) > radius
    return im