Difference between revisions of "FreeSurferToFieldTrip"
(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.
Contents
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