123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182 |
- 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;
|