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.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 # 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 plot_stat_map(stat_img_masked) # this is the data array stat_img_masked_data = stat_img_masked.get_fdata()
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, mask_img=mask_img, smoothing_fwhm=None, min_region_size=100) plot_prob_atlas(region_img, view_type='filled_contours', cut_coords=(-3,-25,11))