vb_job_leadfield_extra_eeg

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

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