Back to the main page.

Bug 2556 - ft_sourceparcellate

Status ASSIGNED
Reported 2014-05-01 23:27:00 +0200
Modified 2014-06-27 18:39:29 +0200
Product: FieldTrip
Component: core
Version: unspecified
Hardware: Macintosh
Operating System: Mac OS
Importance: P5 normal
Assigned to: Robert Oostenveld
URL:
Tags:
Depends on:
Blocks:
See also: http://bugzilla.fcdonders.nl/show_bug.cgi?id=2579

Raghavan Gopalakrishnan - 2014-05-01 23:27:52 +0200

I used the following function to read the AAL parcellation aal = ft_read_atlas('/Users/gopalar/Documents/MATLAB/fieldtrip-20140430/template/atlas/aal/ROI_MNI_V4.nii'); aal = dim: [91 109 91] hdr: [1x1 struct] transform: [4x4 double] unit: 'mm' tissue: [91x109x91 double] tissuelabel: {1x116 cell} coordsys: 'mni' I suppose 'aal' should be consisent with the format as given in ft_datatype_parcellation. But, looks different. According to definition, I should be able to use output of ft_read_atlas in ft_sourcepacellate as a parcellation. Then, I tried using the aal to generate parcels for my source structure cfg.method = 'mean'; cfg.parcellation = 'tissue'; cfg.parameter ={'pow'}; parcel = ft_sourceparcellate(cfg, source , aal); First problem, ft_sourceparcellate does not recognize aal as a datatype parcellation (at line 77) Second, line 177 required source.inside to be a non-empty matrix. Otherwise, gives an error Not sure if there is a bug or mistake in my understanding.


Robert Oostenveld - 2014-05-02 09:48:54 +0200

Hi Raghavan, > I suppose 'aal' should be consisent with the format as given in ft_datatype_parcellation. The output of ft_read_atlas should be consistent with either the description in ft_datatype_parcellate (surface), or with ft_datatype_segmentation (volume). In this case it should be the latter, as AAL is a volumetric atlas. > Then, I tried using the aal to generate parcels for my source structure that is indeed how it is supposed to work. Note that the function is rather new and has not been used/tested very much yet. I can confirm the first error for line 77 Error using ft_checkdata (line 368) This function requires parcellation data as input. I will work on this. For the second error I would expect that line 78 would have added a non-empty boolean inside vector. Could you check this?


Robert Oostenveld - 2014-05-02 09:52:27 +0200

A small correction to my previous comment. A segmentation must be a volume, i.e. ft_datatype_segmentation is a special class of ft_datatype_volume (which is always a regular 3d grid). A parcellation however is not necessarily a surface. In general ft_datatype_parcellation is a structure according to ft_datatype_source, which can describe a regular or irregular geometry, such as a surface, but also a volume. The source structure is more flexible in the anatomical description, and allows more dimensions (such as time) to be represented as well.


Robert Oostenveld - 2014-05-02 10:00:15 +0200

I fixed the problem on line 77, it was a simple omission in ft_checkdata. mac011> svn commit test utilities/ Adding test/test_bug2556.m Sending utilities/ft_checkdata.m Transmitting file data .. Committed revision 9464. >> aal = ft_read_atlas(filename) Rescaling NIFTI: slope = 1, intercept = 0 aal = dim: [91 109 91] hdr: [1x1 struct] transform: [4x4 double] unit: 'mm' tissue: [91x109x91 double] tissuelabel: {1x116 cell} coordsys: 'mni' >> aal = ft_checkdata(aal, 'datatype', 'parcellation', 'parcellationstyle', 'indexed') aal = dim: [91 109 91] hdr: [1x1 struct] unit: 'mm' tissue: [91x109x91 double] tissuelabel: {1x116 cell} coordsys: 'mni' pos: [902629x3 double]


Raghavan Gopalakrishnan - 2014-05-05 18:30:59 +0200

Robert, Line 77 problem is now eliminated Line 79 adds a boolean vector. My source looks like this time: [1x1500 double] pos: [8196x3 double] inside: [8196x1 double] outside: [1x0 double] method: 'average' avg: [1x1 struct] cfg: [1x1 struct] When, I try to parcellate, there seems to be multiple problems, First one at line 177, then at 289. The parcellation is much bigger than number of sources available (I am using a cortical sheet, I think I sent you some example source structure before). Please help Thanks, Raghavan


Robert Oostenveld - 2014-05-07 15:07:47 +0200

(In reply to Raghavan Gopalakrishnan from comment #4) You will indeed have problems, because the source and parcellation are not anatomically consistent. Better would be to start with the same anatomical model for the source reconstruction as you want to use for the parcellation. See http://fieldtrip.fcdonders.nl/example/create_single-subject_grids_in_individual_head_space_that_are_all_aligned_in_mni_space, except that here you would want to start with dipoles defined according to the voxel locations in the AAL atlas. I added the following around line 81 % ensure that the source and the parcellation are anatomically consistent if ~isequal(source.pos, parcellation.pos) error('the source positions are not consistent with the parcellation, please use FT_SOURCEINTERPOLATE'); end


Robert Oostenveld - 2014-05-07 15:09:50 +0200

Sorry, I should have added this: it should be possible to do sourceAAL = ft_sourceinterpolate(cfg, source, aal) to interpolate the source reconstruction results on the AAL atlas. Note that this does require them to be coregistered, i.e. the source model (sheet) should be expressed in MNI space, like the AAL atlas.


Raghavan Gopalakrishnan - 2014-05-19 20:55:32 +0200

ft_sourceparcellate still looks for avg.pow structure. The output of sourcegrandaverage no longer has avg.pow field, instead it has only pow field When I used the following structure atlassource_timelock_stim = time: [1x1500 double] pow: [902629x1500 double] pos: [902629x3 double] dim: [91 109 91] inside: [1x902629 double] outside: [1x0 double] coordsys: 'mni' unit: 'mm' cfg: [1x1 struct] gave the following error Warning: the input does not contain an avg or trial field > In ft_checkdata>fixsource at 1327 In ft_checkdata at 715 In ft_sourceparcellate at 79 Warning: parameter "pow" cannot be parcellated > In ft_sourceparcellate at 149 Error using ft_sourceparcellate (line 162) there are no source parameters that can be parcellated When the structure had pow given as 'avg.pow' it works.


Raghavan Gopalakrishnan - 2014-05-19 21:02:33 +0200

The output of ft_sourceanalysis has the pow field inside avg. So, currently, ft_sourceparcellate is compatible with output of ft_sourceanalysis, but not the output of ft_sourcegrandaverage. Can pow field have a standard location either inside of outside avg?


Robert Oostenveld - 2014-05-20 19:11:29 +0200

(In reply to Raghavan Gopalakrishnan from comment #8) the future standard is source.pow, the old standard is source.avg.pow. Migrating all code and documentation and in that process trying to maintain compatibility with users' scripts cannot be done on the flick of a switch. After trying out various things, I managed to reproduce the problem you describe. It appears that specifying cfg.parameter='pow' in ft_sourceparcellate does not work as expected. I have made some changes, please also see the second part of the test function test/test_bug2556.m and try to come with similar small sections of code if it still fails for you. It is a lot of work to pinpoint the specific problem based on a general statement of something not working. A specific example allows me to reproduce it more quickly and understand the specific issue. mac011> svn commit private/getdimord.m ft_sourceparcellate.m test/test_bug2556.m Sending ft_sourceparcellate.m Sending private/getdimord.m Sending test/test_bug2556.m Transmitting file data ... Committed revision 9557.


Raghavan Gopalakrishnan - 2014-06-25 23:16:50 +0200

(In reply to Robert Oostenveld from comment #6) Robert, You mentioned that source should be expressed in mni space inorder to use ft_sourceinterpolate since aal is in mni space. My sources are in neuromag space, and I am not sure how to convert coordsys at this stage. Using ft_convert_coordsys on sources does not seem to work since there is no transform matrix associated that defines neuromag to mni space. Is there a easy way to tackle this problem? On another note, is 'spm' the same as 'mni' space? Thanks, Raghavan


Robert Oostenveld - 2014-06-26 10:46:07 +0200

(In reply to Raghavan Gopalakrishnan from comment #10) spm and min are indeed the same. ft_volumenormalize can be used to determine the spatial deformation between some anatomical MRI and apply that transformation to the anatomical, and to any functional data that is in coregister with the anatomy. That means that the transformation is done volumetrically after source estimation. We consider it more optimal (and more flexible) to make a source grid that is expressed in single subject MEG space, but where each source location maps directly into MNI space. http://fieldtrip.fcdonders.nl/example/create_single-subject_grids_in_individual_head_space_that_are_all_aligned_in_mni_space This allows you to update the source positions from MEG to MNI, keeping the source estimates as they are (i.e. no interpolation). There is also this example script that was provided by a Neuromag user: http://fieldtrip.fcdonders.nl/example/read_neuromag_mri_and_create_single-subject_grids_in_individual_head_space_that_are_all_aligned_in_mni_space


Raghavan Gopalakrishnan - 2014-06-27 18:39:29 +0200

Robert, Thank you vey much. I require a new functionality in ft_sourceparcellate (or may be ft_sourcedescriptives or another function). I understand ft_sourcedescriptives with cfg.projectmom='yes' constraints the sources in the direction of largest eigenvector. And when the results are used with ft_sourceparcellate, I can compute mean (or some function) of the constrained sources with in each parcel. In addition, I would like to sign flip each of the source inside a parcel to align with its dominant direction, before taking the mean (or some another function). Is is possible to implement this in FT? Or may be its already possible to do it and am not sure how to? I am including a bit of code here of what I need. % find direction of largest eigen [Up,Sp,Vp]=svd(GridOrient,'econ'); % compute MNE [U,S,V] = svd(Leadfields,'econ'); s=diag(S); ss = s ./ (s.^2 + lambda); for par=1: N % number of parcels vtx= % vertices (index) of dipoles in that parcel mne_parcel = V(vtx,:)*(diag(ss) * U' * dat; % mne for that parcel mne_parcel_sign_flip=mean(diag(sign(Up(vtx,1))*mne_parcel,1); end Please advise. Thanks, Raghavan