Back to the main page.
Bug 3055 - ft_prepare_headmodel 'dipoli' output does not contain the vol.mat field
| Status | ASSIGNED |
| Reported | 2016-01-29 17:10:00 +0100 |
| Modified | 2016-02-01 13:32:03 +0100 |
| Product: | FieldTrip |
| Component: | core |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Linux |
| Importance: | P5 normal |
| Assigned to: | Robert Oostenveld |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
Diego Lozano Soldevilla - 2016-01-29 17:10:09 +0100
I'm preparing a pipeline for Coimbra workshop to show how to use a standard electrodes/MRIs to do EEG-source reconstruction. I noticed that the 'dipoli' method produce a volume conduction model without the 'vol.mat' field that later produces an error on ft_prepare_leadfield: grid = ft_prepare_leadfield(cfg) using gradiometers specified in the configuration converting units from 'mm' to 'cm' converting units from 'mm' to 'cm' creating dipole grid based on automatic 3D grid with specified resolution using gradiometers specified in the configuration creating dipole grid with 1 cm resolution 1907 dipoles inside, 5199 dipoles outside brain making tight grid 1907 dipoles inside, 1873 dipoles outside brain the call to "ft_prepare_sourcemodel" took 1 seconds and required the additional allocation of an estimated 0 MB computing leadfield Error using eeg_leadfieldb (line 53) there is no BEM system matrix present Error in ft_compute_leadfield (line 431) lf = eeg_leadfieldb(dippos, sens.elecpos, headmodel); Error in ft_prepare_leadfield (line 219) grid.leadfield{thisindx} = ft_compute_leadfield(grid.pos(thisindx,:), sens, headmodel, 'reducerank', cfg.reducerank, 'normalize', cfg.normalize, 'normalizeparam', cfg.normalizeparam, 'backproject The issue is that the template dipoli volume model has the 'vol.mat'(fieldtrip/template/headmodel/standard_bem.mat') and no error is produced. There's a strange error computing the dipoli method that does not produce a matlab error: Fatal error in dipoli: during computation of B-matrix; vertex 916 of /tmp/tp0436fb35_ca6c_4b05_b160_a94786e27dc4_1.tri touches triangle 2962 of /tmp/tp0436fb35_ca6c_4b05_b160_a94786e27dc4_2.tri Warning: an error ocurred while running dipoli > In ft_headmodel_dipoli at 211 In ft_prepare_headmodel at 278 Error using fread Invalid file identifier. Use fopen to generate a valid file identifier. Warning: File '/tmp/tp5316e4d2_d1a3_4402_b1f9_bbd9a98a8249.ama' not found. > In ft_headmodel_dipoli at 219 In ft_prepare_headmodel at 278 the call to "ft_prepare_headmodel" took 3 seconds and required the additional allocation of an estimated 0 MB Please find below the whole pipeline to reproduce the error: ftdir = fileparts(which('ft_defaults')); mri = ft_read_mri(fullfile(ftdir,'template','headmodel','standard_mri.mat')); % used later to plot dipoles cfg = []; cfg.resolution = 1; cfg.xrange = [-100 100]; cfg.yrange = [-140 140]; cfg.zrange = [-80 120]; mris = ft_volumereslice(cfg, mri); mris = ft_convert_units(mris, 'cm'); % segmentation cfg = []; cfg.output = {'brain','skull','scalp'}; seg = ft_volumesegment(cfg,mri); % triangulation cfg = []; cfg.tissue = {'brain','skull','scalp'}; cfg.numvertices = [3000 2000 1000]; cfg.sourceunits = seg.unit; bnd = ft_prepare_mesh(cfg,seg); % headmodel cfg = []; cfg.method = 'dipoli'; vol = ft_prepare_headmodel(cfg,bnd); vol.cond = [1 1/20 1].*0.33; % brain, skull, scalp % visualization figure; ft_plot_mesh(vol.bnd(1),'facecolor','none'); ft_plot_mesh(vol.bnd(2),'facecolor','none'); ft_plot_mesh(vol.bnd(3),'facecolor','none'); % get 64 electrodes (biosemi64 customized layout) label={'Fp1','AF7','AF3','F1','TP9','F5h','F7','FT7','FC5h','PO9','FC1','C1','C3','C5','T7','TP7','CP5','CP3','CP1','P1','I1','P5h','P7','P9','PO7','PO3','O1','Iz','Oz','POz','Pz','CPz','Fpz','Fp2','AF8','AF4','AFz','Fz','F2','TP10','F6h','F8','FT8','FC6h','PO10','FC2','FCz','Cz','C2','C4','C6','T8','TP8','CP6','CP4','CP2','P2','I2','P6h','P8','P10','PO8','PO4','O2'} elec = ft_read_sens(fullfile(ftdir,'template','electrode','standard_1005.elc')); [sel1,sel2] = match_str(elec.label,label); elec.chanpos = elec.chanpos(sel1,:); elec.elecpos = elec.elecpos(sel1,:); elec.label = label; elec.label % check alignment figure; ft_plot_sens(elec,'style','sk'); hold on; ft_plot_mesh(vol.bnd(1),'facealpha', 0.85, 'edgecolor', 'none', 'facecolor', [0.65 0.65 0.65]); %scalp figure; ft_plot_sens(elec,'style','sk'); hold on; ft_plot_mesh(vol.bnd(2),'facealpha', 0.85, 'edgecolor', 'none', 'facecolor', [0.65 0.65 0.65]); %skull figure; ft_plot_sens(elec,'style','sk'); hold on; ft_plot_mesh(vol.bnd(3),'facealpha', 0.85, 'edgecolor', 'none', 'facecolor', [0.65 0.65 0.65]); %brain cfg = []; cfg.grad = elec; cfg.headmodel = vol; cfg.reducerank = 2; cfg.channel = 'all'; cfg.grid.resolution = 1; % use a 3-D grid with a 1 cm resolution cfg.grid.unit = 'cm'; grid = ft_prepare_leadfield(cfg); vol2 = ft_read_vol(fullfile(ftdir,'template','headmodel','standard_bem.mat')); cfg.headmodel = vol2; grid = ft_prepare_leadfield(cfg);
Robert Oostenveld - 2016-01-30 09:54:05 +0100
I am working through your code, replicating the problem. Let me comment along the way: you should not use cfg.sourceunits in ft_prepare_mesh. I rather suggest that you do >> seg = ft_convert_units(seg, 'cm'); It is actually weird that your MRI is (explicit) in cm, but that the output of ft_volumesegment is in mm. That is due to the SPM template (which is internally used) being in mm, but better would be to retain the output in the units of the input.
Robert Oostenveld - 2016-01-30 10:07:20 +0100
(In reply to Robert Oostenveld from comment #1) I looked it up, and there is indeed no explicit unit handling at the end of the function. I have updated the code, so that the output has the same units as the input. mac011> svn commit ft_volumesegment.m Sending ft_volumesegment.m Transmitting file data . Committed revision 11154.
Robert Oostenveld - 2016-01-30 10:20:05 +0100
(In reply to Robert Oostenveld from comment #2) I was wrong, the output was already consistent with the input. I removed the unnecessary line of code. mac011> svn commit ft_volumesegment.m Sending ft_volumesegment.m Transmitting file data . Committed revision 11155.
Robert Oostenveld - 2016-01-30 10:35:03 +0100
I can confirm the problem. I changed the code such that it is not a warning but an error. That makes much more sense. ... Bmat countdown: 5110 (pass 2) Bmat countdown: 5100 (pass 2) Bmat countdown: 5090 (pass 2) Fatal error in dipoli: during computation of B-matrix; vertex 916 of /private/tmp/tp1073ac04_28a4_447d_a830_a28fbae00a99_1.tri touches triangle 811 of /private/tmp/tp1073ac04_28a4_447d_a830_a28fbae00a99_2.tri Error using ft_headmodel_dipoli (line 211) an error ocurred while running the dipoli executable - please look at the screen output Error in ft_prepare_headmodel (line 278) headmodel = ft_headmodel_dipoli(geometry, 'conductivity', cfg.conductivity, 'isolatedsource', cfg.isolatedsource); ------- mac011> svn commit forward/ft_headmodel_dipoli.m Sending forward/ft_headmodel_dipoli.m Transmitting file data . Committed revision 11156.
Robert Oostenveld - 2016-01-30 10:55:08 +0100
The actual error (mesh not being BEM compatible) can be visualized with % ft_plot_mesh(vol.bnd(1),'facecolor','none'); % skin ft_plot_mesh(vol.bnd(2),'facecolor','none'); % skull ft_plot_mesh(vol.bnd(3),'facecolor','white','edgecolor','none'); % brain hold on ft_plot_mesh(bnd(3).pnt(916,:), 'vertexmarker', '*', 'vertexcolor', 'r'); tmp = vol.bnd(2); tmp.tri = tmp.tri(811,:); ft_plot_mesh(tmp,'facecolor','g'); % skull The problem is at the base of the skull, where it touches the brain.
Robert Oostenveld - 2016-01-30 11:02:58 +0100
Created attachment 769 screen shot showing hole (in axial slice) If you do seg.anatomy = mri.anatomy cfg = [] cfg.funparameter = 'skull' ft_sourceplot(cfg, seg) you can see at the bottom that the skull is not closed.
Diego Lozano Soldevilla - 2016-01-30 18:10:51 +0100
(In reply to Robert Oostenveld from comment #6) Hi Robert, Thanks a lot for your prompt reply! Good catch: I see the skull problem. I've been playing around with the segmentation cfgs (brainsmooth, scalpsmooth, brainthreshold, scalpthreshold) but I was unable to solve the issue. However, I checked the stardard segementation file and it's fine although there's no cfg info to know how to reproduce it. In addition, the format it's a bit peculiar: seg = ft_read_mri(fullfile(ftdir,'template','headmodel','standard_seg.mat')); seg = dim: [181 217 181] anatomy: [181x217x181 uint8] transform: [4x4 double] hdr: [1x1 struct] unit: 'mm' coordsys: 'spm' seg: [181x217x181 double] This seg representation produces an error in ft_prepare_mesh. To solve it, I did the following and it seems I got proper meshes: brain = zeros(size(seg.seg)); bri = find(seg.seg==3); brain(bri)=1; skull = zeros(size(seg.seg)); ski = find(seg.seg==2); skull(ski)=1; scalp = zeros(size(seg.seg)); sci = find(seg.seg==1); scalp(sci)=1; seg.brain = brain; seg.scalp = scalp; seg.skull = skull; seg = rmfield(seg,'seg'); cfg = []; cfg.tissue = {'brain','skull','scalp'}; cfg.numvertices = [3000 2000 1000]; cfg.sourceunits = seg.unit; bnd = ft_prepare_mesh(cfg,seg); figure; ft_plot_mesh(bnd(1),'facecolor','none'); ft_plot_mesh(bnd(2),'facecolor','none'); ft_plot_mesh(bnd(3),'facecolor','none'); Unfortunately, I don't know if dipoli method works after that because I'm testing in windows and there's no dipoli mex file compiled (my remote connection to my windows machine doesn't work). On monday I'll test this further. If you know the parameters to get a proper segmentation with the standard 'colin.mri', let me know. I've planned to do an example script (wiki) to show a basic dipole fitting and beamforming on EEG data using template MRI/electrodes. Lots of people are using these approach but I cannot get all the details from some papers that used fieldtrip.