Difference between revisions of "FreeSurferToFieldTrip"

From WikiMEG
Jump to: navigation, search
(to CTF coord sys)
(BEM surfaces: add help functions)
 
(4 intermediate revisions by one user not shown)
Line 19: Line 19:
 
= In MATLAB =
 
= In MATLAB =
  
 +
== Load data ==
 
Read the Nifti
 
Read the Nifti
  
 
<pre>t1 = ft_read_mri('fs/anattms_50/mri/T1.mgz.nii')</pre>
 
<pre>t1 = ft_read_mri('fs/anattms_50/mri/T1.mgz.nii')</pre>
 +
 +
Set the voxel to head space transform, using the file generated previously
 +
 +
<pre>t1.transform = load('fs/anattms_50/transforms/t1-vox2ras-tkr.xfm')</pre>
  
 
Read the BEM surfaces (this function is given below)
 
Read the BEM surfaces (this function is given below)
Line 34: Line 39:
 
plot_bem_surfs(bems);
 
plot_bem_surfs(bems);
 
</pre>
 
</pre>
 +
 +
== 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.  
 
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.  
Line 47: Line 54:
  
 
<pre>bems = update_bems_coord_sys(bems, t1.transform/t1.transformorig);</pre>
 
<pre>bems = update_bems_coord_sys(bems, t1.transform/t1.transformorig);</pre>
 +
 +
Now verify again everything's aligned
 +
 +
<pre>
 +
ft_determine_coordsys(t1, 'interactive', 'no')
 +
hold on
 +
plot_bem_surfs(bems);
 +
</pre>
 +
 +
== building the FieldTrip head model ==
 +
 +
= Helper functions =
 +
 +
== BEM surfaces ==
 +
 +
<pre>
 +
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;
 +
</pre>
 +
 +
<pre>
 +
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
 +
</pre>
 +
 +
<pre>
 +
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
 +
</pre>

Latest revision as of 16: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