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())
m = cm.ScalarMappable(cmap=cmap,norm=norm)

# plot histogram 

enter image description here enter image description here

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.masking import apply_mask
from nilearn.masking import unmask
from nilearn.plotting import plot_stat_map

# download freely available mask image and statistical image from nilearn
mask_img = fetch_icbm152_brain_gm_mask()
stat_img = fetch_localizer_calculation_task(n_subjects=1).cmaps[0]

# mask statistical image to grey matter only
stat_img_resampled = resample_to_img(stat_img,mask_img)
stat_img_data_masked = apply_mask(stat_img_resampled,mask_img)
stat_img_masked = unmask(stat_img_data_masked,mask_img)

# plot the statistical image

# this is the data array
stat_img_masked_data = stat_img_masked.get_fdata()

enter image description here

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

region_img,_ = connected_regions(stat_img_masked,


enter image description here

  • $\begingroup$ "Same spatial structure" is too vague to support any useful answer. Could you clarify what that means, preferably in a quantitative manner? $\endgroup$ – whuber Apr 8 at 15:47
  • $\begingroup$ I updated the question and tried to make it more clear what I mean $\endgroup$ – Johannes Wiesner Apr 12 at 12:10
  • $\begingroup$ Is my understanding correct that the question can be reduced to: How do I simulate spatially correlated data? $\endgroup$ – Frans Rodenburg Apr 12 at 12:43
  • $\begingroup$ Yes, you can! I would like to simulate statistical MRI images from a given one. $\endgroup$ – Johannes Wiesner Apr 12 at 15:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.