vb_prepare_bem

 Preparation for leadfield BEM calculation

  [BEM, Vhead, Fhead, Omega, Sout, SPHinfo] = vb_prepare_bem(basis_parm)

 --- Output
  BEM      : Parameters for EEG-BEM 
  Vhead    : Vertex position of the surface model (NV x 3)
  Fhead    : Patch data of the surface model (NP x 3)
  Sout     : Normal of boudary surface (NV x 3
  Omega    : Solid angle matrix        
  SPHinfo  : Parameters for Spherical Harmonics

 --- 'BEM' should be saved in headfile for EEG-BEM case
      BEM = []  for MEG cases
 Nsurf         : number of boundary surface = 3 for EEG
 BEM.sigma     : ( 1 x Nsurf)
               = conductivity of shell region (各領域の伝導率 )
               = [sigma(1), ..., sigma(Nsurf)]
                 from innermost to outermost

 BEM.R         : ( 1 x Nsurf)
               = Radius of each surface (for Sphere model EEG case)
               = [R(1), ..., R(Nsurf)]
 
 BEM.Nvertex   : ( Nsurf x 2 )
               = start & end vertex index of each surface
               = [ start_id(1)     end_id(1)     ; 
                             ...                 ;
                   start_id(Nsurf) end_id(Nsurf) ]
               = [開始インデックス, 終了インデックス]

 BEM.Npatch    : ( Nsurf x 2 )
               = start & end patch index of each surface
               = [開始インデックス, 終了インデックス]

 written by M. Sato  2005-8-8
 modified by T. Sako 2006-8-23
 modified by T. Sako 2006-12-12 - for EEG 3-shell model

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function    [BEM, Vhead, Fhead, Omega, Sout, SPHinfo] = vb_prepare_bem(basis_parm)
0002 % Preparation for leadfield BEM calculation
0003 %
0004 %  [BEM, Vhead, Fhead, Omega, Sout, SPHinfo] = vb_prepare_bem(basis_parm)
0005 %
0006 % --- Output
0007 %  BEM      : Parameters for EEG-BEM
0008 %  Vhead    : Vertex position of the surface model (NV x 3)
0009 %  Fhead    : Patch data of the surface model (NP x 3)
0010 %  Sout     : Normal of boudary surface (NV x 3
0011 %  Omega    : Solid angle matrix
0012 %  SPHinfo  : Parameters for Spherical Harmonics
0013 %
0014 % --- 'BEM' should be saved in headfile for EEG-BEM case
0015 %      BEM = []  for MEG cases
0016 % Nsurf         : number of boundary surface = 3 for EEG
0017 % BEM.sigma     : ( 1 x Nsurf)
0018 %               = conductivity of shell region (各領域の伝導率 )
0019 %               = [sigma(1), ..., sigma(Nsurf)]
0020 %                 from innermost to outermost
0021 %
0022 % BEM.R         : ( 1 x Nsurf)
0023 %               = Radius of each surface (for Sphere model EEG case)
0024 %               = [R(1), ..., R(Nsurf)]
0025 %
0026 % BEM.Nvertex   : ( Nsurf x 2 )
0027 %               = start & end vertex index of each surface
0028 %               = [ start_id(1)     end_id(1)     ;
0029 %                             ...                 ;
0030 %                   start_id(Nsurf) end_id(Nsurf) ]
0031 %               = [開始インデックス, 終了インデックス]
0032 %
0033 % BEM.Npatch    : ( Nsurf x 2 )
0034 %               = start & end patch index of each surface
0035 %               = [開始インデックス, 終了インデックス]
0036 %
0037 % written by M. Sato  2005-8-8
0038 % modified by T. Sako 2006-8-23
0039 % modified by T. Sako 2006-12-12 - for EEG 3-shell model
0040 %
0041 % Copyright (C) 2011, ATR All Rights Reserved.
0042 % License : New BSD License(see VBMEG_LICENSE.txt)
0043 
0044 % --- Boundary surface information
0045 %  Vhead : Vertex position of the surface model (NV x 3)
0046 % XXhead : Normal vector at vertices (NV x 3)
0047 %  Fhead : Patch index of the surface model (NP x 3)
0048 
0049 bem_mode = basis_parm.bem_mode;
0050 
0051 flg_joined_head = false;
0052 if isfield(basis_parm, 'head_file') && ~iscell(basis_parm.head_file)
0053   headfile = [basis_parm.proj_root filesep basis_parm.head_file];
0054 else
0055   headfile = [];
0056   flg_joined_head = true;
0057 end
0058 
0059 global vbmeg_inst;
0060 const = vbmeg_inst.const;
0061 % const.BASIS_MEG_SPHERE;   %  MEG Sphere model (Sarvas)
0062 % const.BASIS_MEG_BEM;      %  MEG BEM
0063 % const.BASIS_MEG_HARMONICS;%  MEG Spherical harmonics expansion
0064 % const.BASIS_EEG_SPHERE;   %  EEG 3-shell Sphere model
0065 % const.BASIS_EEG_BEM;      %  EEG 3-shell BEM
0066 
0067 switch bem_mode
0068 case const.BASIS_MEG_SPHERE,
0069     %  MEG Sphere model (Sarvas)
0070     BEM   = [];
0071     Vhead = [];
0072     Fhead = [];
0073     Omega = [];
0074     Sout  = [];
0075 
0076     SPHinfo = [];
0077 
0078 case const.BASIS_MEG_BEM, 
0079     %  MEG BEM 1 shell model
0080     
0081     % load head file
0082     load([headfile]);
0083     fprintf('# of vertices = %d\n',size(Vhead,1));
0084     
0085     if  ~exist('Omega','var') | basis_parm.Recalc == ON,
0086           tic;
0087           fprintf('--- Solid angle calculation ');
0088     
0089           % 境界面係数行列計算
0090           % (Normal vectors (XXhead) are outwarded here)
0091           [Omega, Sout]= vb_solid_angle_grk(Vhead,Fhead,XXhead);
0092     
0093           fprintf('%f[sec]\n',toc);
0094           vb_save([headfile],'Omega','Sout');
0095     end
0096     
0097     BEM = [];
0098     SPHinfo = [];
0099 
0100 case const.BASIS_MEG_HARMONICS, 
0101     %  MEG Spherical harmonics expansion
0102     
0103     % load head file
0104     load([headfile], 'Vhead', 'Fhead', 'XXhead');
0105     
0106     fprintf('# of vertices = %d\n',size(Vhead,1));
0107     
0108     % load meg file
0109     meg_file = [basis_parm.proj_root filesep basis_parm.meg_file];
0110     [pick, Qpick, CoilWeight, Vcenter] = vb_load_sensor(meg_file);
0111     MEGinfo = vb_load_meg_info(meg_file);
0112     
0113     if ~isfield(MEGinfo,'SPHinfo') | basis_parm.Recalc == ON,
0114        
0115         % Center of Head
0116         if basis_parm.DEBUG == ON,
0117             % Setting for DEBUG to check Leadfield calculation
0118             Hcenter = [0 0 0];
0119         else
0120             Hcenter = mean(Vhead);
0121         end
0122 
0123         % average radius from center of head
0124 %         Rmax = vb_average_radius(Hcenter, Vhead);
0125 
0126         % Change center of coordinate
0127         Vhead = [Vhead(:,1)-Hcenter(1),  ...
0128                  Vhead(:,2)-Hcenter(2), Vhead(:,3)-Hcenter(3)];
0129         pick  = [ pick(:,1)-Hcenter(1),   ...
0130                   pick(:,2)-Hcenter(2), pick(:,3)-Hcenter(3)];
0131 
0132     % test
0133         Rmax = vb_average_radius(Hcenter, Vhead);
0134 
0135     % function order
0136         func_order = basis_parm.func_order;
0137         
0138         % Calculate spherical harmonics coefficient
0139         SPHinfo.harmo_coef = ...
0140             vb_spherical_harmo_coef(Vhead, XXhead, func_order, Rmax, pick, Qpick);
0141 
0142         % Save SPHinfo to MEGinfo
0143         SPHinfo.func_order = func_order;
0144         SPHinfo.radius     = Rmax;
0145         MEGinfo.SPHinfo    = SPHinfo;
0146         vb_save(meg_file, 'MEGinfo');
0147     else
0148         % Saved result
0149         SPHinfo = MEGinfo.SPHinfo;
0150     end
0151     
0152     BEM   = [];
0153     Omega = [];
0154     Sout  = [];
0155 
0156 case const.BASIS_EEG_SPHERE, 
0157   % EEG 3-shell Sphere model
0158 
0159   % R    : Relative radii of sphere from innermost to outermost
0160   % sigma: Conductivity  from innermost to outermost
0161     
0162     BEM.R     = basis_parm.radius;
0163     BEM.sigma = basis_parm.sigma;
0164     
0165     Vhead = [];
0166     Fhead = [];
0167     Omega = [];
0168     Sout  = [];
0169     SPHinfo = [];
0170   
0171 case const.BASIS_EEG_BEM, 
0172   %  EEG BEM 3-shell model
0173     
0174   % check fields
0175   required_fields = {'head_file', 'sigma'};
0176   optional_fields = {'LayerTag'};
0177 
0178   for i = 1:length(required_fields)
0179     if ~isfield(basis_parm, required_fields{i})
0180       error('incomplete basis_parm : does not exist : %s', required_fields{i});
0181     end
0182   end
0183   
0184   for i = 1:length(optional_fields)
0185     if ~isfield(basis_parm, optional_fields{i})
0186       warning('incomplete basis_parm : does not exist : %s', optional_fields{i});
0187       basis_parm.LayerTag = char('CSF', 'Skull', 'Scalp');
0188     end
0189   end
0190 
0191   % load head file
0192   new_headfile = vb_util_make_joined_head_file( ...
0193     basis_parm.proj_root, basis_parm.head_file, basis_parm.sigma, basis_parm.LayerTag);
0194   
0195     load([new_headfile]);
0196     fprintf('# of vertices = %d\n',size(Vhead,1));
0197   
0198     if ~exist('BEM','var'), 
0199     % make BEM from information of headfile
0200     necessary_info = {'R', 'Sigma', 'Nvertex', 'Npatch'};
0201     if exist('R', 'var') && exist('Sigma', 'var') && exist('Nvertex', 'var') && exist('Npatch', 'var')
0202       BEM.R       = R;
0203       BEM.sigma   = [Sigma 0];  % for vb_bem_matrix
0204       BEM.Nvertex = Nvertex;
0205       BEM.Npatch  = Npatch;
0206     elseif  ~flg_joined_head
0207       if ~exist('Vhead', 'var') || ~exist('Fhead', 'var')
0208         error('Head file must have Vhead and Fhead field.(%s)\n', new_headfile);
0209       end
0210       
0211       % single shell head
0212       if ~exist('R', 'var')             BEM.R = calc_average_radius(Vhead); end;
0213       if ~exist('Sigma', 'var')     BEM.sigma = [1.0 0];                    end;
0214       if ~exist('Nvertex', 'var') BEM.Nvertex = [1 length(Vhead)];          end;
0215       if ~exist('Npatch', 'var')   BEM.Npatch = [1 length(Fhead)];          end;
0216     else
0217       % output error information
0218       for i = 1:length(necessary_info)
0219         if ~exist(necessary_info{i}, 'var')
0220           warning('The headfile must have <%s>.\n', necessary_info{i});
0221         end
0222       end
0223       error('Making BEM is incomplete. Check --> %s', headfile);
0224     end
0225   end
0226     
0227   if  ~exist('Omega','var') | basis_parm.Recalc == ON,
0228     tic;
0229     fprintf('--- Solid angle calculation ');
0230     
0231     % 境界面係数行列計算
0232     % (Normal vectors (XXhead) are outwarded here)
0233     [Omega, Sout]= vb_solid_angle_grk(Vhead,Fhead,XXhead);
0234     
0235     fprintf('%f[sec]\n',toc);
0236       
0237     omega_file = vb_util_make_omega_file(new_headfile, Omega, Sout);
0238   end
0239   SPHinfo = [];
0240 end
0241 
0242 % --- END OF FILE --- %