////// Preprocess ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Concatenate all the subjects
fmri_data = np.concatenate(dataset.func, axis=3)
// Apply a small amount of Gaussian smoothing: in the case of ICA it is
// important as it introduces a spatial model that ICA lacks and greatly
// reduces the high-frequency signal
from scipy import ndimage
for image in fmri_data.T:
// This works efficiently because image is a view on fmri_data
image[...] = ndimage.gaussian_filter(image, 1.5)
// Take the mean along axis 3: the direction of time
mean_img = np.mean(fmri_data, axis=3)
// Mask non brain areas
from nisl import masking
mask = masking.compute_epi_mask(mean_img)
data_masked = fmri_data[mask]
////// Apply ICA //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
from sklearn.decomposition import FastICA
n_components = 20
ica = FastICA(n_components=n_components, random_state=42)
components_masked = ica.fit(data_masked).transform(data_masked)
// We normalize the estimated components, for thresholding to make sens
components_masked -= components_masked.mean(axis=0)
components_masked /= components_masked.std(axis=0)
// Threshold
components_masked[np.abs(components_masked) < .5] = 0
// Now we inverting the masking operation, to go back to a full 3D
// representation
(x, y, z) = mean_img.shape
components = np.zeros((x, y, z, n_components))
components[mask] = components_masked
// Using a masked array is important to have transparency in the figures
components = np.ma.masked_equal(components, 0, copy=False)
After Change
fmri_data = data_masked
// Take the mean along axis 3: the direction of time
mean_img = masker.inverse_transform(fmri_data.mean(axis=-1))
////// Apply ICA //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////