bem_generate_eeg_transfer_matrix.m

% bem_generate_eeg_transfer_matrix() - Generates the EEG transfer matrix.
%
% Usage:
%   >> session = bem_generate_eeg_transfer_matrix(session);
%
% Inputs:
%   session - session structure generated by BEM_CREATE_SESSION
%
% Outputs:
%   session - session object; it may have been modified while
%             loading the model matrices.
%
% Notes: The matrix is created using the name given in the session structure
%       with "tmte" extension for the EEG transfer matrix. It can be loaded
%       using bem_load_transfer_matrix().
%
% Author: Zeynep Akalin Acar, SCCN, 2007

% Copyright (C) 2007 Zeynep Akalin Acar, SCCN, zeynep@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


function session = bem_generate_eeg_transfer_matrix(session)

% check session
if ~isempty(find(isfield(session, {'name', 'model', 'Smatrix'}) == 0,1))
    error('BEM:bem_generate_eeg_transfer_matrix:session','%s','Invalid session');
end
model = session.model;
if ~isempty(find(isfield(model, {'name', 'mesh'}) == 0,1))
    error('BEM:bem_generate_eeg_transfer_matrix:model','%s','Invalid model');
end
mesh = model.mesh;
if ~isempty(find(isfield(mesh, {'name','num_nodes'}) == 0,1))
    error('BEM:bem_generate_eeg_transfer_matrix:mesh','%s','Invalid mesh');
end
if ~isfield(model, 'cmt')
    session.model = bem_load_model_matrix(model, 'cmt');
    model = session.model;
end

if session.num_electrodes < 50
    hh=waitbar(0,'calculating rows of the transfer matrix');
    parfor k = 1:session.num_electrodes
        waitbar(k/session.num_electrodes);
        disp(sprintf('BEM:status: solving for electrode %d', k));
        b = zeros(mesh.num_nodes,1);
        ind = find(session.Smatrix(:,1) == k);
        b(session.Smatrix(ind,2))=session.Smatrix(ind,3);
        x = gmres(model.cmt', b, 100);
        tmte(:,k) = x;
    end
    close(hh)
else
    for k = 1:session.num_electrodes
        b = zeros(mesh.num_nodes,1);
        ind = find(session.Smatrix(:,1) == k);
        b(session.Smatrix(ind,2))=session.Smatrix(ind,3);
        bi(k,:) = b;
    end
    iD=inv(model.cmt); tmte = bi*iD; tmte=tmte';
end       

fn = sprintf('%s.tmte', session.name);
save(fn, 'tmte', '-ascii', '-double');