Dipole estimation by minimum error search
[Vdipole,Jdipole,err,iter] = vb_job_dipole(proj_root, dipole_parm)
--- Input
proj_root : data directory
dipole_parm : parameter structure
- Input MEG file
dipole_parm.megfile = MEG file (relative path from proj_root)
- Initial dipole position
dipole_parm.Vinit = Initial dipole position [Ndipole x 3]
number of dipoles in minimum error search is given by the size of Vinit
- Error function
dipole_parm.errfunc = Error function mode for dipole search
= 0 : 'vb_dipole_error': for external current (Biot-Savart Eq.)
= 1 : 'vb_dipole_error_sarvas': for cortical current (Sarvas Eq.)
= 2 : 'vb_dipole_error_magdipole': for coil current (magnetic dipole)
- Optional parameter
dipole_parm.dipolefile = Output file (relative path from proj_root)
if this field is not exist or empty, no output file is produced.
dipole_parm.twin_meg = Time window for analysis
if this field is not exist or empty, whole time period is used.
dipole_parm.trial_id = trial number to estimate
if this field is not exist or empty, all trials are averaged.
- Optimization setting (Optional)
dipole_parm.MaxIter = 5000; % max iteration number in optimization
dipole_parm.TolX = 1.e-8; % Xの終了許容値
dipole_parm.TolFun = 1.e-10; % 関数値の終了許容値
dipole_parm.mode = 1;
= 0; % Use 'fminsearch' in Matlab
= 1; % Use 'fminunc' in Optimization Toolbox
= 2; % Use 'fminunc1'
--- Output
Vdipole : Dipole position [Ndipole x 3]
Jdipole : dipole current [Ndipole x 3 x Time]
Jdipole(n ,k ,t) = current at n-th dipole ,time 't' and k-th direction
err : normalized error
iter : iteration number
2007-7-6 made by M.Sato
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [Vdipole,Jdipole,err,iter] = vb_job_dipole(proj_root, dipole_parm) 0002 % Dipole estimation by minimum error search 0003 % [Vdipole,Jdipole,err,iter] = vb_job_dipole(proj_root, dipole_parm) 0004 % --- Input 0005 % proj_root : data directory 0006 % dipole_parm : parameter structure 0007 % - Input MEG file 0008 % dipole_parm.megfile = MEG file (relative path from proj_root) 0009 % - Initial dipole position 0010 % dipole_parm.Vinit = Initial dipole position [Ndipole x 3] 0011 % number of dipoles in minimum error search is given by the size of Vinit 0012 % - Error function 0013 % dipole_parm.errfunc = Error function mode for dipole search 0014 % = 0 : 'vb_dipole_error': for external current (Biot-Savart Eq.) 0015 % = 1 : 'vb_dipole_error_sarvas': for cortical current (Sarvas Eq.) 0016 % = 2 : 'vb_dipole_error_magdipole': for coil current (magnetic dipole) 0017 % - Optional parameter 0018 % dipole_parm.dipolefile = Output file (relative path from proj_root) 0019 % if this field is not exist or empty, no output file is produced. 0020 % dipole_parm.twin_meg = Time window for analysis 0021 % if this field is not exist or empty, whole time period is used. 0022 % dipole_parm.trial_id = trial number to estimate 0023 % if this field is not exist or empty, all trials are averaged. 0024 % - Optimization setting (Optional) 0025 % dipole_parm.MaxIter = 5000; % max iteration number in optimization 0026 % dipole_parm.TolX = 1.e-8; % Xの終了許容値 0027 % dipole_parm.TolFun = 1.e-10; % 関数値の終了許容値 0028 % dipole_parm.mode = 1; 0029 % = 0; % Use 'fminsearch' in Matlab 0030 % = 1; % Use 'fminunc' in Optimization Toolbox 0031 % = 2; % Use 'fminunc1' 0032 % --- Output 0033 % Vdipole : Dipole position [Ndipole x 3] 0034 % Jdipole : dipole current [Ndipole x 3 x Time] 0035 % Jdipole(n ,k ,t) = current at n-th dipole ,time 't' and k-th direction 0036 % err : normalized error 0037 % iter : iteration number 0038 % 0039 % 2007-7-6 made by M.Sato 0040 % 0041 % Copyright (C) 2011, ATR All Rights Reserved. 0042 % License : New BSD License(see VBMEG_LICENSE.txt) 0043 0044 proj_root = vb_rm_trailing_slash(proj_root); 0045 0046 % --- Input file 0047 megfile = [proj_root '/' dipole_parm.megfile]; 0048 0049 % Error function name for dipole search 0050 % err_mode 0051 % = 0 : 'vb_dipole_error': for external current (Biot-Savart Eq.) 0052 % = 1 : 'vb_dipole_error_sarvas': for cortical current (Sarvas Eq.) 0053 % = 2 : 'vb_dipole_error_magdipole': for external current (magnetic dipole) 0054 err_mode = dipole_parm.errfunc; 0055 0056 switch err_mode, 0057 case 0 0058 errfunc = 'vb_dipole_error'; 0059 case 1 0060 errfunc = 'vb_dipole_error_sarvas'; 0061 case 2 0062 errfunc = 'vb_dipole_error_magdipole'; 0063 otherwise 0064 error('error function mode is incorrect') 0065 end 0066 0067 % Vinit : initial dipole position (Ndipole x 3) 0068 Vinit = dipole_parm.Vinit; 0069 0070 % Load MEG data 0071 bexp = vb_load_meg_data(megfile); 0072 [N, T, M] = size(bexp); 0073 0074 % set time window & trial number to estimate 0075 if isfield(dipole_parm,'trial_id') & ~isempty(dipole_parm.trial_id) 0076 trial_id = dipole_parm.trial_id; 0077 else 0078 trial_id = 1:M; 0079 end 0080 0081 if isfield(dipole_parm,'twin_meg') & ~isempty(dipole_parm.twin_meg) 0082 t_begin = dipole_parm.twin_meg(1); 0083 t_end = dipole_parm.twin_meg(2); 0084 else 0085 t_begin = 1; 0086 t_end = T; 0087 end 0088 0089 % MEG data 0090 bexp = bexp(:, t_begin:t_end, trial_id); 0091 bexp = sum(bexp,3)/size(bexp,3); 0092 0093 % Load sensor info 0094 [pick,Qpick,Wsensor,V0] = vb_load_sensor(megfile); 0095 0096 % Change center of coordinate 0097 pick = [pick(:,1)-V0(1), pick(:,2)-V0(2), pick(:,3)-V0(3)]; 0098 0099 % Optimization method 0100 if isfield(dipole_parm,'mode') & ~isempty(dipole_parm.mode) 0101 opt_mode = dipole_parm.mode; 0102 else 0103 opt_mode = vb_tool_check('Optimization'); 0104 % = 0; % Use 'fminsearch' 0105 % = 1; % Use 'fminunc' in Optimization Toolbox 0106 end 0107 0108 % Dipole position search 0109 [Vdipole, err, iter] = ... 0110 vb_dipole_search(Vinit,bexp,pick,Qpick,Wsensor,errfunc,dipole_parm,opt_mode); 0111 0112 % Dipole current time cource for fixed position dipole 0113 [Jdipole, err] = ... 0114 vb_dipole_current(Vdipole, bexp, pick, Qpick, Wsensor, [0 0 0], err_mode); 0115 0116 % Right-hand SPM (m) coordinate 0117 Vdipole = [ Vdipole(:,1)+V0(1), Vdipole(:,2)+V0(2), Vdipole(:,3)+V0(3)]; 0118 0119 % Save result 0120 % --- Output file 0121 if isfield(dipole_parm,'dipolefile') & ~isempty(dipole_parm.dipolefile) 0122 dipolefile = [proj_root '/' dipole_parm.dipolefile]; 0123 vb_save(dipolefile, 'Vdipole','Jdipole','dipole_parm','err'); 0124 end