From 3d3a82e865445f15862c5f4804e6351159a79a37 Mon Sep 17 00:00:00 2001 From: Jan-Mathijs Schoffelen <1517611+schoffelen@users.noreply.github.com> Date: Mon, 11 Nov 2024 16:45:10 +0100 Subject: [PATCH] ENH - allow volumetric source space to be read (#2464) * ENH - allow volumetric source space to be read * ENH - test function --- fileio/ft_read_headshape.m | 132 ++++++++++++++++++++----------------- fileio/private/pos2dim.m | 2 +- fileio/private/pos2dim3d.m | 10 +-- test/test_pull2464.m | 13 ++++ 4 files changed, 89 insertions(+), 68 deletions(-) create mode 100644 test/test_pull2464.m diff --git a/fileio/ft_read_headshape.m b/fileio/ft_read_headshape.m index b7bb46f496..74183ec894 100644 --- a/fileio/ft_read_headshape.m +++ b/fileio/ft_read_headshape.m @@ -440,74 +440,82 @@ case 'mne_source' % read the source space from an MNE file ft_hastoolbox('mne', 1); + + try + % a fif-file can also contain a source space that is volumetric, in which case the below function call will fail (due to the add_geom + % being specified as true, but the file does not contain triangulation information. strictly speaking, the fif-file then + % does not represent a headshape, but as a service to the casual user, let's support the reading of this type of file + src = mne_read_source_spaces(filename, 1); + catch + src = mne_read_source_spaces(filename); + end - src = mne_read_source_spaces(filename, 1); + if numel(src)==2 + if ~isempty(annotationfile) + ft_hastoolbox('freesurfer', 1); + if numel(annotationfile)~=2 + ft_error('two annotationfiles expected, one for each hemisphere'); + end + for k = 1:numel(annotationfile) + [v{k}, label{k}, c(k)] = read_annotation(annotationfile{k}, 1); + end - if ~isempty(annotationfile) - ft_hastoolbox('freesurfer', 1); - if numel(annotationfile)~=2 - ft_error('two annotationfiles expected, one for each hemisphere'); - end - for k = 1:numel(annotationfile) - [v{k}, label{k}, c(k)] = read_annotation(annotationfile{k}, 1); + % match the annotations with the src structures + if src(1).np == numel(label{1}) && src(2).np == numel(label{2}) + src(1).labelindx = label{1}; + src(2).labelindx = label{2}; + elseif src(1).np == numel(label{2}) && src(1).np == numel(label{1}) + src(1).labelindx = label{2}; + src(2).labelindx = label{1}; + else + ft_warning('incompatible annotation with triangulations, not using annotation information'); + end + if ~isequal(c(1),c(2)) + ft_error('the annotation tables differ, expecting equal tables for the hemispheres'); + end + c = c(1); end - % match the annotations with the src structures - if src(1).np == numel(label{1}) && src(2).np == numel(label{2}) - src(1).labelindx = label{1}; - src(2).labelindx = label{2}; - elseif src(1).np == numel(label{2}) && src(1).np == numel(label{1}) - src(1).labelindx = label{2}; - src(2).labelindx = label{1}; - else - ft_warning('incompatible annotation with triangulations, not using annotation information'); + shape = []; + % only keep the points that are in use + inuse1 = src(1).inuse==1; + inuse2 = src(2).inuse==1; + shape.pos=[src(1).rr(inuse1,:); src(2).rr(inuse2,:)]; + + % only keep the triangles that are in use; these have to be renumbered + newtri1 = src(1).use_tris; + newtri2 = src(2).use_tris; + for i=1:numel(src(1).vertno) + newtri1(newtri1==src(1).vertno(i)) = i; + end + for i=1:numel(src(2).vertno) + newtri2(newtri2==src(2).vertno(i)) = i; end - if ~isequal(c(1),c(2)) - ft_error('the annotation tables differ, expecting equal tables for the hemispheres'); + shape.tri = [newtri1; newtri2 + numel(src(1).vertno)]; + if isfield(src(1), 'use_tri_area') + shape.area = [src(1).use_tri_area(:); src(2).use_tri_area(:)]; end - c = c(1); - end - - shape = []; - % only keep the points that are in use - inuse1 = src(1).inuse==1; - inuse2 = src(2).inuse==1; - shape.pos=[src(1).rr(inuse1,:); src(2).rr(inuse2,:)]; - - % only keep the triangles that are in use; these have to be renumbered - newtri1 = src(1).use_tris; - newtri2 = src(2).use_tris; - for i=1:numel(src(1).vertno) - newtri1(newtri1==src(1).vertno(i)) = i; - end - for i=1:numel(src(2).vertno) - newtri2(newtri2==src(2).vertno(i)) = i; - end - shape.tri = [newtri1; newtri2 + numel(src(1).vertno)]; - if isfield(src(1), 'use_tri_area') - shape.area = [src(1).use_tri_area(:); src(2).use_tri_area(:)]; - end - if isfield(src(1), 'use_tri_nn') - shape.nn = [src(1).use_tri_nn; src(2).use_tri_nn]; - end - shape.orig.pos = [src(1).rr; src(2).rr]; - shape.orig.tri = [src(1).tris; src(2).tris + src(1).np]; - shape.orig.inuse = [src(1).inuse src(2).inuse]'; - shape.orig.nn = [src(1).nn; src(2).nn]; - if isfield(src(1), 'labelindx') - shape.orig.labelindx = [src(1).labelindx;src(2).labelindx]; - shape.labelindx = [src(1).labelindx(inuse1); src(2).labelindx(inuse2)]; - % ulabelindx = unique(c.table(:,5)); - % for k = 1:c.numEntries - % % the values are really high (apart from the 0), so I guess it's safe to start - % % numbering from 1 - % shape.orig.labelindx(shape.orig.labelindx==ulabelindx(k)) = k; - % shape.labelindx(shape.labelindx==ulabelindx(k)) = k; - % end - % FIXME the above screws up the interpretation of the labels, because the color table is not sorted - shape.label = c.struct_names; - shape.annotation = c.orig_tab; % to be able to recover which one - shape.ctable = c.table; + if isfield(src(1), 'use_tri_nn') + shape.nn = [src(1).use_tri_nn; src(2).use_tri_nn]; + end + shape.orig.pos = [src(1).rr; src(2).rr]; + shape.orig.tri = [src(1).tris; src(2).tris + src(1).np]; + shape.orig.inuse = [src(1).inuse src(2).inuse]'; + shape.orig.nn = [src(1).nn; src(2).nn]; + if isfield(src(1), 'labelindx') + shape.orig.labelindx = [src(1).labelindx;src(2).labelindx]; + shape.labelindx = [src(1).labelindx(inuse1); src(2).labelindx(inuse2)]; + shape.label = c.struct_names; + shape.annotation = c.orig_tab; % to be able to recover which one + shape.ctable = c.table; + end + else + ft_warning('the fif-file did not seem to contain triangulation information, probably you try to read a volumetric source space'); + + shape = []; + shape.pos = src.rr; + shape.inside = src.inuse(:)>0; + shape.dim = pos2dim3d(src.rr); end case {'neuromag_fif' 'neuromag_mne'} diff --git a/fileio/private/pos2dim.m b/fileio/private/pos2dim.m index 328601d8e1..918df4dd3c 100644 --- a/fileio/private/pos2dim.m +++ b/fileio/private/pos2dim.m @@ -27,7 +27,7 @@ dim = nan(1, 3); [tmp, ind] = max(dpos, [], 2); dim(1) = find(tmp>1.5, 1, 'first'); -dpos = dpos(dim:dim:npos-1, :); +dpos = dpos(dim(1):dim(1):npos-1, :); [tmp, ind] = max(dpos(:, setdiff(1:3, ind(dim(1)))), [], 2); dim(2) = find(tmp>1.1*min(tmp), 1, 'first'); % this threshold seems to work on what I tried out dim(3) = npos./prod(dim(1:2)); diff --git a/fileio/private/pos2dim3d.m b/fileio/private/pos2dim3d.m index 757b14f84a..453bf3016c 100644 --- a/fileio/private/pos2dim3d.m +++ b/fileio/private/pos2dim3d.m @@ -17,22 +17,22 @@ % Copyright (C) 2009, Jan-Mathijs Schoffelen -if nargin==1 && ~isstruct(pos), +if nargin==1 && ~isstruct(pos) dimold = zeros(0,2); -elseif isstruct(pos), +elseif isstruct(pos) % the input is a FieldTrip data structure dimord = pos.dimord; dimtok = tokenize(dimord, '_'); for i = 1:length(dimtok) - if strcmp(dimtok{i},'pos'), + if strcmp(dimtok{i},'pos') dimold(i,1) = size(pos.pos,1); else - dimold(i,1) = numel(getfield(pos, dimtok{i})); + dimold(i,1) = numel(pos.(dimtok{i})); end end pos = pos.pos; else - if size(pos,1)~=dimold(1), + if size(pos,1)~=dimold(1) ft_error('the first element in the second input should be equal to the number of positions'); end end diff --git a/test/test_pull2464.m b/test/test_pull2464.m new file mode 100644 index 0000000000..e1429f7cf0 --- /dev/null +++ b/test/test_pull2464.m @@ -0,0 +1,13 @@ +function test_pull2464 + +% MEM 1gb +% WALLTIME 00:10:00 +% DEPENDENCY ft_read_headshape +% DATA private + +%% + +datadir = dccnpath('/home/common/matlab/fieldtrip/data/test/'); +fname = fullfile(datadir, 'pull2464.fif'); +hs = ft_read_headshape(fname, 'format', 'mne_source'); +