Latest revision as of 17:20, 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

function bem = read_bem_surfs(subjects_dir, subject)

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

for i=1:length(surfs)
	if ~strcmp(surfs{i}, 'head')
		bname = [subject '_' surfs{i} '_surface'];
		fname = fullfile(subjects_dir, subject, 'bem', 'watershed', bname);

	% not used for the moment, this is a high resolution surface (~300k vert)
	% whereas outer_skin is good enough
	else
		bname = 'lh.seghead';
		fname = fullfile(subjects_dir, subject, 'surf', bname);
	end
	[v, f] = freesurfer_read_surf(fname);
	bem.(surfs{i}).vertices = v;
	bem.(surfs{i}).faces = f;
end

bem.surfs = surfs;

function ps = plot_bem_surfs(bems)
%
% bems = read_bem_surfs(..)
% ps = plot_bem_surfs(bems);

surfs = {'brain', 'inner_skull', 'outer_skull', 'outer_skin'};
colors = {[1 1 1] [1 0 0] [0 1 0] [0 0 1]};

for i=1:length(surfs)
	si = surfs{i};
	p = patch(bems.(si));
	set(p, 'edgealpha', 0.0);
	set(p, 'facealpha', 0.2);
	set(p, 'facecolor', colors{i});
	ps.(si) = p;
end

axis equal

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