Create ROIs at random locations with fixed number of voxels

Great that you managed !

If you add this constraint (each cluster must be a connected component) then I think you need to check nilearn/regions/rena_clustering.py especially the recursive_neighbor_agglomeration() function.

The matrix “connectivity” is what you need to enforce this local connectedness. Then the aggregation scheme has to be tailored with your constraints.

Good luck