# Generate data that matches a frequency distribution while preserving the original spatial structure

I am dealing with a 3D array containing values representing the "importance" of each voxel. For my analysis, I would like to synthesize n new arrays from my original array to have a comparison condition. The values in the 3D array can be clustered into spatially distinct regions, meaning that there is a positive relationship between the spatial proximity and value similarity between two voxels. Just randomly shuffling the values to create a synthetic array would create an array without the spatial structure of the original array thus leading to an "unfair" comparison. Instead, I am looking for a way to create a synthetic 3D array from my original array that has the same value distribution while somehow conserving the non-random spatial distribution of the voxel values.

Here's some code to visualize my problem:

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate
from matplotlib import cm

rng = np.random.RandomState(41)

# Make this bigger to generate a more dense grid
N = 10

# Generate data
x,y,z,d = rng.random((4,20))
xi = np.linspace(x.min(),x.max(),N),
yi = np.linspace(y.min(),y.max(),N),
zi = np.linspace(z.min(),z.max(),N)
xi,yi,zi = np.meshgrid(xi,yi,zi)
rbf = scipy.interpolate.Rbf(x,y,z,d,function='linear')
di = rbf(xi,yi,zi)

# plot voxels
fig = plt.figure()
ax = fig.gca(projection='3d')
cmap = plt.get_cmap('hsv')
norm= plt.Normalize(di.min(),di.max())
ax.voxels(np.ones_like(di),facecolors=cmap(norm(di)))
m = cm.ScalarMappable(cmap=cmap,norm=norm)
m.set_array([])
plt.colorbar(m)
plt.show()

# plot histogram
plt.hist(di.ravel())



As you can see, the values in this cube data are not randomly distributed but tend to spatially cluster together (I guess 'local entropy' is the right keyword here?). Now is there a way to create n other cubes from this one given my constraints? Note: The cube above is just a simplified example for my question. In reality, I would like to apply this to a statistical brain image. Here's some code to get such a data array:

from nilearn.datasets import fetch_icbm152_brain_gm_mask,fetch_localizer_calculation_task
from nilearn.image import resample_to_img
from nilearn.plotting import plot_stat_map

# mask statistical image to grey matter only

# plot the statistical image

# this is the data array


We can cluster this image using a random walker segmentation provided by nilearn

# we can cluster the masked statistical image using a random walker
# segmentation algorithm
# NOTE: This can take a while!
from nilearn.regions import connected_regions
from nilearn.plotting import plot_prob_atlas