0001 function [basis, basis_parm] = ...
0002 vb_job_leadfield_meg_harmonics(proj_root, basis_parm)
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024 [proj_root, basis_parm] = inner_check_arguments(proj_root, basis_parm);
0025
0026
0027
0028
0029
0030 [Basis_mode, pick, Qpick, Wsensor, V0, V, xx] = ...
0031 vb_lf_preprocess(proj_root, basis_parm);
0032
0033
0034 head_file = [proj_root filesep basis_parm.head_file];
0035 load([head_file], 'Vhead', 'Fhead', 'XXhead');
0036
0037
0038 if basis_parm.DEBUG == ON,
0039
0040 Hcenter = [0 0 0];
0041 else
0042 Hcenter = mean(Vhead);
0043 end
0044
0045
0046 Vhead = [Vhead(:,1)-Hcenter(1), Vhead(:,2)-Hcenter(2), Vhead(:,3)-Hcenter(3)];
0047 V = [ V(:,1)-Hcenter(1), V(:,2)-Hcenter(2), V(:,3)-Hcenter(3)];
0048 pick = [ pick(:,1)-Hcenter(1), pick(:,2)-Hcenter(2), pick(:,3)-Hcenter(3)];
0049
0050
0051 [V,xx] = vb_current_vector(V, xx, Basis_mode);
0052
0053 fprintf('# of vertices = %d\n',size(Vhead,1));
0054
0055
0056 meg_file = [proj_root filesep basis_parm.meg_file];
0057 MEGinfo = vb_load_measurement_info(meg_file);
0058
0059 if ~isfield(MEGinfo,'SPHinfo') | basis_parm.Recalc == ON,
0060
0061
0062 Rmax = vb_average_radius(Hcenter, Vhead);
0063
0064
0065 func_order = basis_parm.func_order;
0066
0067
0068 SPHinfo.harmo_coef = ...
0069 vb_spherical_harmo_coef(Vhead, XXhead, func_order, Rmax, pick, Qpick);
0070
0071
0072 SPHinfo.func_order = func_order;
0073 SPHinfo.radius = Rmax;
0074 SPHinfo.head_file = basis_parm.head_file;
0075 MEGinfo.SPHinfo = SPHinfo;
0076
0077 vb_save(meg_file, 'MEGinfo');
0078 else
0079
0080 SPHinfo = MEGinfo.SPHinfo;
0081 end
0082
0083
0084 fprintf('--- MEG Spherical Harmonics expansion ');
0085 N = SPHinfo.func_order;
0086 Rmax = SPHinfo.radius;
0087 A = SPHinfo.harmo_coef;
0088
0089 BB = vb_spherical_harmo_magnetic_field(V, xx, N, Rmax, pick, Qpick, A);
0090
0091 basis = BB * Wsensor';
0092
0093
0094
0095
0096
0097
0098
0099 function [proj_root, basis_parm] = inner_check_arguments(proj_root, basis_parm)
0100 if ~isfield(basis_parm, 'func_order'),
0101
0102 basis_parm.func_order = 35;
0103 end
0104
0105
0106
0107 return;
0108
0109
0110
0111