Revision as of 16:17, 5 May 2015

The FreeSurfer software provides a good automatic segmentation of T1 images. It includes utilities to generate BEM surfaces that can be used to construct a head model for FieldTrip.

The following assumes you've run FreeSurfer on your T1.

Post-recon steps in the shell

Generate BEM surfaces (brain, inner/outer skull & outer skin)

SUBJECTS_DIR=~/mris-tms/fs SUBJECT=anattms_50 mne_watershed_bem

Convert T1 to nifti

mri_convert fs/anattms_50/mri/T1.mgz fs/anattms_50/mri/T1.mgz.nii

Obtain transform from T1 voxels to BEM surface space

mri_info --vox2raw-tkr fs/anattms_50/mri/T1.mgz > fs/anattms_50/transforms/t1-vox2ras-tkr.xfm

In MATLAB

Load data

Read the Nifti

t1 = ft_read_mri('fs/anattms_50/mri/T1.mgz.nii')

Set the voxel to head space transform, using the file generated previously

t1.transform = load('fs/anattms_50/transforms/t1-vox2ras-tkr.xfm')

Read the BEM surfaces (this function is given below)

bems = read_bem_surfs('fs/', 'anattms_50');

Plot the two together

ft_determine_coordsys(t1, 'interactive', 'no')
hold on
plot_bem_surfs(bems);

To CTF coordinates

At this point, the T1 and BEM surfaces are in the same coord sys but FieldTrip doesn't know which sys this is, and more importantly, a registration is required with the sensors.

Identify fiducial points (however you like; I use a simple imagesc based function https://gist.github.com/maedoc/97f7d4d4e30c9123e748)

fids = gofish(t1);
t1 = ft_volumerealign(fids, t1);

Now, the T1 has a CTF based coord system, but the BEM surfaces need to have their vertices updated (again, this is defined below):

bems = update_bems_coord_sys(bems, t1.transform/t1.transformorig);

Now verify again everything's aligned

ft_determine_coordsys(t1, 'interactive', 'no')
hold on
plot_bem_surfs(bems);

building the FieldTrip head model

Helper functions

BEM surfaces

update_bems_coord_sys

function bems = update_bems_coord_sys(bems, xfm)

surfs = {'brain', 'inner_skull', 'outer_skull', 'outer_skin'};

for i=1:length(surfs)
	si = surfs{i};
	v = bems.(si).vertices;
	v = [v ones(size(v, 1), 1)];
	v = v*xfm';
	bems.(si).vertices = v(:, 1:3);
end