Create 3D blobs with fixed sizes in a predefined subspace

My question is similar to this post (and possibly also related to this one).

I would like to create random 3D blobs in a predefined subspace. Specifically, my subspace is defined by a binary 3D mask image that 'lives' inside a bigger 3D cube. This mask image again defines the space where the blobs are supposed to 'live' in.

There are a few constraints that should be met:

1.) It should be possible to define the number of blobs.

2.) It should be possible to define the size of each blob as a number of voxels. In other words, each blob should take up a defined fraction of my subspace.

3.) Blobs should not be allowed to touch or overlap (that would violate the first rule because it decreases the defined number of blobs)

4.) Each blob should be a set of connected voxels, i.e. each voxel must be directly connected to another voxel.

I came up with code to obtain a mask_coords array that holds all coordinates defining the subspace. Now I need an algorithm that assigns some of these coordinates to distinct blobs meeting my constraints:

import numpy as np
from nilearn import masking
from nilearn.image.resampling import coord_transform
from nilearn.plotting import plot_roi
from nilearn.image import new_img_like
from nilearn.datasets import fetch_icbm152_brain_gm_mask

# get a grey matter mask image from nilearn. This is a binary 3D image that
# defines our subspace in our bigger 3D cube. 
mask_img = fetch_icbm152_brain_gm_mask()

## Get the x,y,z coordinates of our mask image. ###############################

# get mask data array and affine
mask_data,affine = masking._load_mask_img(mask_img)

# get data indices for all '1' voxels inside mask data array
mask_data_indices = np.asarray(np.nonzero(mask_data)).T.tolist()

# return coordinates for those '1' voxels. This gives us all coordinates that
# correspond to our subspace
mask_coords = np.asarray(list(zip(*mask_data_indices)))

mask_coords = coord_transform(mask_coords[0],
                              mask_coords[1],
                              mask_coords[2],
                              affine)

mask_coords = np.asarray(mask_coords).T

Here's also some code that we can use to come back to a 3D numpy array and that we can use to plot our 3D mask image:

## We can convert the mask coordinates back to data array #####################

mask_indices = coord_transform(mask_coords[:,0],
                               mask_coords[:,1],
                               mask_coords[:,2],
                               np.linalg.inv(affine))

mask_indices = np.asarray(mask_indices,dtype=int).T

empty_array = np.zeros(mask_data.shape)

for indices in mask_indices.tolist():
    empty_array[indices[0],indices[1],indices[2]] = True
    
mask_img_from_indices = new_img_like(mask_img,empty_array)

# show results
plot_roi(mask_img_from_indices,bg_img=None,draw_cross=False)