The procedure is as follows:
- compute the results we want to morph
- select reference spaces (in our case the subject’s individual T1 and “fsaverage”, which comes with mne or FreeSurfer). Note that we only use the “brain”, excluding skull and other surrounding parts. A brain only image can be obtained from the FreeSurfer segmentation.
- morph subject brain to average brain
- apply the resulting transform to the volumetric grid obtained from the lcmv beamformer
pip install mne dipy nibabel nilearn
But first let’s import some packages:
import numpy as np import mne from mne.datasets import sample from mne.beamformer import make_lcmv, apply_lcmv from dipy.align.imwarp import SymmetricDiffeomorphicRegistration from dipy.align.metrics import CCMetric from dipy.align.reslice import reslice from dipy.align.transforms import (TranslationTransform3D, RigidTransform3D, AffineTransform3D) from dipy.align.imaffine import (transform_centers_of_mass, MutualInformationMetric, AffineRegistration) import nibabel as nib from nilearn.plotting import plot_anat from nilearn.image import index_img, resample_img
# Setup paths data_path = sample.data_path() raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif' event_fname = data_path + '/MEG/sample/sample_audvis_raw-eve.fif' fname_fwd = data_path + '/MEG/sample/sample_audvis-meg-vol-7-fwd.fif' # Get epochs event_id, tmin, tmax = [1, 2], -0.2, 0.5 # Setup for reading the raw data raw = mne.io.read_raw_fif(raw_fname, preload=True) raw.info['bads'] = ['MEG 2443', 'EEG 053'] # 2 bad channels events = mne.read_events(event_fname) # Set up pick list: gradiometers and magnetometers, excluding bad channels picks = mne.pick_types(raw.info, meg=True, eeg=False, stim=True, eog=True, exclude='bads') # Pick the channels of interest raw.pick_channels([raw.ch_names[pick] for pick in picks]) # Re-normalize our empty-room projectors, so they are fine after subselection raw.info.normalize_proj() # Read epochs proj = False # already applied epochs = mne.Epochs(raw, events, event_id, tmin, tmax, baseline=(None, 0), preload=True, proj=proj, reject=dict(grad=4000e-13, mag=4e-12, eog=150e-6)) evoked = epochs.average() # Read regularized noise covariance and compute regularized data covariance noise_cov = mne.compute_covariance(epochs, tmin=tmin, tmax=0, method='shrunk') data_cov = mne.compute_covariance(epochs, tmin=0.04, tmax=0.15, method='shrunk') ############################################################################### # Read forward model forward = mne.read_forward_solution(fname_fwd) # Compute weights of free orientation (vector) beamformer with weight # normalization (neural activity index, NAI). Providing a noise covariance # matrix enables whitening of the data and forward solution. Source orientation # is optimized by setting pick_ori to 'max-power'. # weight_norm can also be set to 'unit-noise-gain'. Source orientation can also # be 'normal' (but only when using a surface-based source space) or None, # which computes a vector beamfomer. Note, however, that not all combinations # of orientation selection and weight normalization are implemented yet. filters = make_lcmv(evoked.info, forward, data_cov, reg=0.05, noise_cov=noise_cov, pick_ori='max-power', weight_norm='nai') # Apply this spatial filter to the evoked data. stc = apply_lcmv(evoked, filters, max_ori_out='signed') # take absolute values for plotting stc.data[:, :] = np.abs(stc.data) # Save result in stc files stc.save('lcmv-vol') # select time window (tmin, tmax) stc.crop(0.0, 0.0) # Save result in a 4D nifti file mne.save_stc_as_volume('lcmv_inverse.nii.gz', stc, forward['src'], mri_resolution=True)
Note, that we set “mri_resolution” to “true”. This helps us later on in morphing the grid, since the (non-linear) transform will be estimated based on the anatomical MRI(s).
Next, we will prepare our data for the morph. This includes loading the respective volumes, defining a new voxel size (note, that because this is an example, a voxel size was chosen towards optimizing performance, rather than resolution), and normalizing our two reference T1s. Furthermore we will define the number of iterations used for each optimizations step (explained below). In order to keep things comprehensible, we will from now on call the fsaverage T1 image Static, because it is the volume we want to move to. The subject’s T1 in turn will be called Moving, since it will be moved towards Static:
# isometric voxel size of 3 mm - reconsider for real data voxel_size = (3., 3., 3.) # number of iterations for each optimisation level niter_affine = [100, 100, 10] niter_sdr = [5, 5, 3] # load lcmv inverse img_vol = nib.load('lcmv_inverse.nii.gz') # reslice lcmv inverse img_vol_res, img_vol_res_affine = reslice(img_vol.get_data(), img_vol.affine, img_vol.header.get_zooms()[:3], voxel_size) img_vol = nib.Nifti1Image(img_vol_res, img_vol_res_affine) img_vol_first = index_img(img_vol, 0) # for plotting # select first time volume and output ndarray img = img_vol.dataobj[:, :, :, 0] # load subject brain (Moving) t1_fname = data_path + '/subjects/sample/mri/brain.mgz' t1_m_img = nib.load(t1_fname) # reslice Moving t1_m_img_res, t1_m_img_res_affine = reslice(t1_m_img.get_data(), t1_m_img.affine, t1_m_img.header.get_zooms()[:3], voxel_size) t1_m_img = nib.Nifti1Image(t1_m_img_res, t1_m_img_res_affine) # get Moving to world transform t1_m_grid2world = t1_m_img.affine # output Moving as ndarray t1_m = t1_m_img.dataobj[:, :, :] # normalize values t1_m = t1_m.astype('float') / t1_m.max() # load fsaverage brain (Static) t1_fname = data_path + '/subjects/fsaverage/mri/brain.mgz' t1_s_img = nib.load(t1_fname) # reslice Static t1_s_img_res, t1_s_img_res_affine = reslice(t1_s_img.get_data(), t1_s_img.affine, t1_s_img.header.get_zooms()[:3], voxel_size) t1_s_img = nib.Nifti1Image(t1_s_img_res, t1_s_img_res_affine) # get Static to world transform t1_s_grid2world = t1_s_img.affine # output Static as ndarray t1_s = t1_s_img.dataobj[:, :, :] # normalize values t1_s = t1_s.astype('float') / t1_s.max()
Now we can compute the transform. First we will use an affine transformation, to get both images pre-aligned as good as possible. Based on that result, the non-linear transform will follow.
Step 1 – affine transform:
# compute center of mass c_of_mass = transform_centers_of_mass(t1_s, t1_s_grid2world, t1_m, t1_m_grid2world) nbins = 32 # prepare affine registration affreg = AffineRegistration(metric=MutualInformationMetric(nbins, None), level_iters=niter_affine, sigmas=[3.0, 1.0, 0.0], factors=[4, 2, 1]) # translation translation = affreg.optimize(t1_s, t1_m, TranslationTransform3D(), None, t1_s_grid2world, t1_m_grid2world, starting_affine=c_of_mass.affine) # rigid body transform (translation + rotation) rigid = affreg.optimize(t1_s, t1_m, RigidTransform3D(), None, t1_s_grid2world, t1_m_grid2world, starting_affine=translation.affine) # affine transform (translation + rotation + scaling) affine = affreg.optimize(t1_s, t1_m, AffineTransform3D(), None, t1_s_grid2world, t1_m_grid2world, starting_affine=rigid.affine) # apply affine transformation t1_m_affine = affine.transform(t1_m)
We shift the center of mass of Moving to Static. Afterwards, we compute a rigid body transform, that is preserving the distance between each set of points. Hence the object will be rotated, translated and reflected to match Static even better. In a third step, we will finally compute a full affine transform, that additionally includes scaling ans shearing.
We use the mutual information as an optimization metric as common in the field (Mattes et al., 2003).
Our two brains should overlap quite nicely by now. However, we still have the problem of “individuality”. Even though our brains overlap very well by now, structural differences across different individuals are still not accounted.
Step 2 – non-linear transform:
# set up Symmetric Diffeomorphic Registration (metric, iterations per level) sdr = SymmetricDiffeomorphicRegistration(CCMetric(3), niter_sdr) # compute mapping mapping = sdr.optimize(t1_s, t1_m_affine) # morph img data img_sdr_affine = mapping.transform(affine.transform(img))
The result is a displacement map, telling exactly which voxel in Moving corresponds to which voxel in Static and how Moving has to be moved to match Static. We will use the cross-correlation optimization metric as described in Avants et al. (2009).
Step 3 – plotting, saving, happiness:
# make transformed ndarray a nifti img_vol_sdr_affine = nib.Nifti1Image(img_sdr_affine, affine=t1_s_grid2world) # save morphed grid nib.save(img_vol_sdr_affine, 'lcmv_inverse_fsavg.nii.gz') # plot result imgs = [img_vol_first, img_vol_first, nib.load('lcmv_inverse_fsavg.nii.gz')] t1_imgs = [t1_m_img, t1_s_img, t1_s_img] titles = ['subject brain', 'fsaverage brian', 'fsaverage brian\nmorphed source'] saveas = ['lcmv_subject.png', 'lcmv_fsaverage.png', 'lcmv_fsaverage_morphed.png'] for img, t1_img, title, fname in zip(imgs, t1_imgs, titles, saveas): display = plot_anat(t1_img, display_mode='x', title=title) display.add_overlay(img) display.savefig(fname) display.close()
The result of the above code will vary slightly. In the example 3 PNG files will be saved in order to visualize the transform. For reasons of eye candy I transformed the results into GIF images:
Here you see how the subject’s anatomical T1 and the result of the source analysis overlap nicely:
When overlaying the FreeSurfer average brain with the very same overlay, then we see that it does not fit. The spaces are off. The same would happen, if we tried to apply the result-space of subject A to the brain of subject B. Hence in order to make results comparable, they must fit a common space (in this case fsaverage):
After estimating the volumetric morph based on the anatomical scans, we can now transform our beamformer results to the same space. By applying first the affine transform and afterwards the non-linear transform, we shift our individual results into a common space (fsaverage):
If we would do this for a set of subjects, we could compute any statistics as usual. Because anatomical (/ functional) areas overlap now between subjects.
The next step would be to provide this example as an official example as part of the MNE-Python documentation.
Afterwards it’s time to go into the API that the will make the workflow accessible from higher order functions. Thereby it could be linked directly to the existing surface morph, such that the user experiences no difference in the actual handling of the function. The magic will be done “under the hood”.
Alright! Have fun trying it out and stay tuned 😉
Avants, B. B., Epstein, C. L., Grossman, M., & Gee, J. C. (2009). Symmetric Diffeomorphic Image Registration with Cross- Correlation: Evaluating Automated Labeling of Elderly and Neurodegenerative Brain, 12(1), 26-41.
Mattes, D., Haynor, D. R., Vesselle, H., Lewellen, T. K., & Eubank, W. (2003). PET-CT image registration in the chest using free-form deformations. IEEE transactions on medical imaging, 22(1), 120-128.