vb_job_leadfield_extra_eeg

 EEG Lead field matrix for extra dipole using BEM
 
 --- Syntax
 vb_job_leadfield_extra_eeg(extra_basis_parm)
 vb_job_leadfield_extra_eeg(proj_root, extra_basis_parm)    [old style]

 --- Input
 Fields of extra_basis_parm
 .basis_file: Leadfield matrix is saved as this file
 .eeg_file  : to load sensor position 'pick'
 .mps_file  : Extra dipole model file (.mps.mat)
 .head_file : Head shell model file for BEM
 .sigma     : Conductivity  from innermost to outermost
 .Recalc    : (Optional) re-calculation flag for solid angle matrix 'Omega'
   if 'Omega' file exist, time consuming 'Omega' calculation is not done
   if Recalc==ON, 'Omega' is re-calculated even of 'Omega' file exist

 --- Output
 Basis file
 basis( Norient * Nvertex , Nsensor)
       basis( n, k )   : k-th sensor field for dipole current at n

 --- History
 2011-09-06 Masa-aki Sato

0001 function vb_job_leadfield_extra_eeg(varargin)
0002 % EEG Lead field matrix for extra dipole using BEM
0003 %
0004 % --- Syntax
0005 % vb_job_leadfield_extra_eeg(extra_basis_parm)
0006 % vb_job_leadfield_extra_eeg(proj_root, extra_basis_parm)    [old style]
0007 %
0008 % --- Input
0009 % Fields of extra_basis_parm
0010 % .basis_file: Leadfield matrix is saved as this file
0011 % .eeg_file  : to load sensor position 'pick'
0012 % .mps_file  : Extra dipole model file (.mps.mat)
0013 % .head_file : Head shell model file for BEM
0014 % .sigma     : Conductivity  from innermost to outermost
0015 % .Recalc    : (Optional) re-calculation flag for solid angle matrix 'Omega'
0016 %   if 'Omega' file exist, time consuming 'Omega' calculation is not done
0017 %   if Recalc==ON, 'Omega' is re-calculated even of 'Omega' file exist
0018 %
0019 % --- Output
0020 % Basis file
0021 % basis( Norient * Nvertex , Nsensor)
0022 %       basis( n, k )   : k-th sensor field for dipole current at n
0023 %
0024 % --- History
0025 % 2011-09-06 Masa-aki Sato
0026 
0027 if length(varargin) == 1
0028   proj_root = [];
0029   extra_basis_parm = varargin{1};
0030 elseif length(varargin) == 2
0031   proj_root = varargin{1};
0032   extra_basis_parm = varargin{2};
0033 end
0034 
0035 % Filename
0036 mpsfile = fullfile(proj_root, extra_basis_parm.mps_file);
0037 eegfile = fullfile(proj_root, extra_basis_parm.eeg_file);
0038 
0039 % Prepare extra dipole information
0040 load(mpsfile,'Pointlist');
0041 extra_pos = [];
0042 extra_direction = [];
0043 
0044 %%%% ---- This part may be changed ---- %%%%
0045 for i=1:length(Pointlist)
0046   % Three dipoles at each position
0047   extra_pos = [extra_pos; Pointlist{i}.point];
0048   extra_pos = [extra_pos; Pointlist{i}.point];
0049   extra_pos = [extra_pos; Pointlist{i}.point];
0050   extra_direction = [extra_direction; 1 0 0];
0051   extra_direction = [extra_direction; 0 1 0];
0052   extra_direction = [extra_direction; 0 0 1];
0053 end
0054 
0055 % Load EEG sensor
0056 [pick] = vb_load_sensor(eegfile);
0057 
0058 % sigma  : Conductivity  from innermost to outermost
0059 sigma = extra_basis_parm.sigma;
0060 
0061 if iscell(extra_basis_parm.head_file)
0062     headfile = extra_basis_parm.head_file;
0063     joinfile = [vb_join_cell_file_name(headfile) '.head.mat'];
0064     joinpath = fullfile(proj_root, joinfile);
0065     
0066     if ~exist(joinpath,'file')
0067         sigma=[];layer=[];mode = 3; % MAX mode for radius
0068         vb_head_join_files(headfile,proj_root,joinpath, sigma,layer,mode);
0069     end
0070     extra_basis_parm.head_file = joinfile;
0071 end
0072 
0073 new_headfile = fullfile(proj_root, extra_basis_parm.head_file);
0074 
0075 omega_file   = [vb_get_basename(new_headfile,'.head.mat') '.omega.mat'];
0076 
0077 fprintf('Load %s\n',new_headfile)
0078 load(new_headfile);
0079 % Vhead : (# of vertex) x 3
0080 % Fhead : (# of patch) x 3
0081 % Nvertex : (# of surface) x 2
0082 % Npatch  : (# of surface) x 2
0083 % Sigma   : 1 x (# of surface)
0084 
0085 fprintf('--- BEM for EEG\n')
0086 fprintf('# of vertices = %d\n',size(Vhead,1));
0087 
0088 if ~exist('Nvertex', 'var'), Nvertex = [1 length(Vhead)]; end;
0089 if ~exist('Npatch', 'var'),  Npatch  = [1 length(Fhead)]; end;
0090 
0091 Nsurf = size(Nvertex,1);
0092 fprintf('# of surfaces = %d\n',Nsurf);
0093 
0094 BEM.sigma   = [sigma(1:Nsurf) 0];
0095 BEM.Nvertex = Nvertex;
0096 BEM.Npatch  = Npatch;
0097 
0098 if exist(omega_file,'file')
0099     fprintf('Load %s\n',omega_file)
0100     load(omega_file)
0101 end
0102 
0103 if  ~exist('Omega','var') || extra_basis_parm.Recalc == ON,
0104     tic;
0105     fprintf('--- Solid angle calculation \n');
0106         
0107     % (Normal vectors (XXhead) are outwarded here)
0108     [Omega, Sout] = vb_solid_angle_grk(Vhead,Fhead,XXhead);
0109 
0110     vb_ptime(toc);
0111     vb_fsave(omega_file,'Omega', 'Sout');
0112     fprintf('Saved omega file: %s\n',omega_file);
0113 end
0114 
0115 tic;
0116 
0117 % Multilayer model correction
0118 fprintf('--- Multisurface correction\n');
0119 [Omega] = vb_bem_matrix(Omega,BEM);
0120 
0121 % Self-weighting factor
0122 %   : Add integrating result of the basis function on the triangular surface.
0123 fprintf('--- Auto solid angle calculation\n');
0124 Omega = vb_solid_auto_grk(Omega,Fhead,Sout);
0125 
0126 vb_ptime(toc);
0127 
0128 fprintf('--- EEG potential (BEM) \n');
0129     
0130 % Electric potential
0131 basis = vb_bem_eeg( Omega, Sout, Vhead, Fhead, ...
0132         extra_pos, extra_direction, pick, BEM );
0133 basis = basis';
0134 
0135 % Save data
0136 vb_fsave(fullfile(proj_root, extra_basis_parm.basis_file),'extra_basis_parm','basis', ...
0137          'Pointlist');