Dear Yekaterina,
I think you are right. I made a fix as in the attached version. I hope
it doesn't break anything.
Thanks,
Vladimir
On Wed, Feb 7, 2018 at 3:36 PM, Lyzhko Ekaterina <[log in to unmask]> wrote:
> Dear Vladimir
>
> I have a question about function spm_eeg_lgainmat.m
> The line '159 save(D);' doesn't save 'scale' of forward model
> Maybe if line '179 D.inv{val}.forward = forward;' will be before 159 it will better?
>
> Best regards,
> Ekaterina Lyzhko
>
>
function [L,D] = spm_eeg_lgainmat(D,Is,channels)
% Load or compute if necessary a gain matrix
% FORMAT [L,D] = spm_eeg_lgainmat(D,Is,channels)
% D - Data structure
% Is - indices of vertices
%
% L - Lead-field or gain matrix L(:,Is)
%__________________________________________________________________________
% Copyright (C) 2008-2017 Wellcome Trust Centre for Neuroimaging
% Karl Friston
% $Id: spm_eeg_lgainmat.m 7255 2018-02-07 22:06:07Z vladimir $
SVNrev = '$Rev: 7255 $';
%-Get gain or lead-field matrix
%--------------------------------------------------------------------------
val = D.val;
forward = D.inv{val}.forward;
for ind = 1:numel(forward)
modality = forward(ind).modality;
%-Channels
%----------------------------------------------------------------------
if isequal(modality, 'MEG')
chanind = D.indchantype({'MEG', 'MEGPLANAR'}, 'GOOD');
else
chanind = D.indchantype(modality, 'GOOD');
end
if ~isempty(chanind)
forward(ind).channels = D.chanlabels(chanind);
else
error(['No good ' modality ' channels were found.']);
end
end
if nargin < 3
channels = [forward(:).channels];
end
try
fname = D.inv{val}.gainmat;
G = load(fullfile(D.path, fname)); % Relative path
label = G.label;
G = G.G;
if numel(label) ~= size(G, 1) || ~all(ismember(channels, label))
error('Gain matrix has an incorrect number of channels.');
end
catch
spm('sFnBanner', mfilename, SVNrev);
spm('Pointer', 'Watch');
G = {};
label = {};
for ind = 1:numel(forward)
%-Create a new lead-field matrix
%==================================================================
%-Head Geometry (create tesselation file)
%------------------------------------------------------------------
vert = forward(ind).mesh.vert;
face = forward(ind).mesh.face;
%-Normals
%------------------------------------------------------------------
norm = spm_mesh_normals(struct('faces',face,'vertices',vert),true);
vol = forward(ind).vol;
if ischar(vol)
vol = ft_read_vol(vol);
end
modality = forward(ind).modality;
if isfield(forward, 'siunits') && forward(ind).siunits
units = D.units(D.indchannel(forward(ind).channels));
sens = forward(ind).sensors;
siunits = isempty(strmatch('unknown', units));
else
siunits = false;
sens = D.inv{val}.datareg(ind).sensors;
end
%-Forward computation
%------------------------------------------------------------------
[vol, sens] = ft_prepare_vol_sens(vol, sens, 'channel', forward(ind).channels);
nvert = size(vert, 1);
spm_progress_bar('Init', nvert, ['Computing ' modality ' leadfields']);
if nvert > 100, Ibar = floor(linspace(1, nvert,100));
else Ibar = [1:nvert]; end
if ~isequal(ft_voltype(vol), 'interpolate')
Gxyz = zeros(length(forward(ind).channels), 3*nvert);
for i = 1:nvert
if siunits
Gxyz(:, (3*i - 2):(3*i)) = ft_compute_leadfield(vert(i, :), sens, vol,...
'dipoleunit', 'nA*m', 'chanunit', units);
else
Gxyz(:, (3*i - 2):(3*i)) = ft_compute_leadfield(vert(i, :), sens, vol);
end
if any(Ibar == i)
spm_progress_bar('Set', i);
end
end
else
if siunits
Gxyz = ft_compute_leadfield(vert, sens, vol, 'dipoleunit', 'nA*m', 'chanunit', units);
else
Gxyz = ft_compute_leadfield(vert, sens, vol);
end
end
spm_progress_bar('Clear');
spm_progress_bar('Init', nvert, ['Orienting ' modality ' leadfields']);
G{ind} = zeros(size(Gxyz, 1), size(Gxyz, 2)/3);
for i = 1:nvert
G{ind}(:, i) = Gxyz(:, (3*i- 2):(3*i))*norm(i, :)';
if ismember(i,Ibar)
spm_progress_bar('Set', i);
end
end
spm_progress_bar('Clear');
%-Condition the scaling of the lead-field
%------------------------------------------------------------------
[Gs, scale] = spm_cond_units(G{ind});
if siunits && abs(log10(scale))>2
warning(['Scaling expected to be 1 for SI units, actual scaling ' num2str(scale)]);
G{ind} = Gs;
else
scale = 1;
end
label = [label; forward(ind).channels(:)];
forward(ind).scale = scale;
end
if numel(G) > 1
G = cat(1, G{:});
else
G = G{1};
end
%-Save
%----------------------------------------------------------------------
D.inv{val}.gainmat = ['SPMgainmatrix_' spm_file(D.fname, 'basename') '_' num2str(val) '.mat'];
save(fullfile(D.path, D.inv{val}.gainmat), 'G', 'label', spm_get_defaults('mat.format'));
D.inv{val}.forward = forward;
save(D);
spm('Pointer', 'Arrow');
end
[sel1, sel2] = spm_match_str(channels, label);
if length(sel2) ~= numel(channels)
error('Did not find a match for all the requested channels.');
end
L = sparse(G(sel2, :));
%-Retain selected sources if necessary
%--------------------------------------------------------------------------
if nargin > 1 && ~isempty(Is)
L = L(:,Is);
end
D.inv{val}.forward = forward;
|