dipole search Dipole estimation by minimum error search Masa-aki Sato 2010-1-15 Copyright (C) 2011, ATR All Rights Reserved. License : New BSD License(see VBMEG_LICENSE.txt)
0001 % dipole search 0002 % Dipole estimation by minimum error search 0003 % 0004 % Masa-aki Sato 2010-1-15 0005 % 0006 % Copyright (C) 2011, ATR All Rights Reserved. 0007 % License : New BSD License(see VBMEG_LICENSE.txt) 0008 0009 clear all 0010 0011 % Number of diple 0012 Ndipole = 1; 0013 % # of trial initial position for diple search 0014 Ntrial = 10; 0015 0016 % Directory of data file 0017 root = 'D:\' ; 0018 root = getenv('MATHOME') ; 0019 proj_root = [root '/NH_SEF/proj_200']; 0020 0021 % MEG data file name (relative path from proj_root) 0022 parm.megfile = 'NH070604a1000_ave.meg.mat'; 0023 % Output file (relative path from proj_root) 0024 parm.dipolefile = 'dipole.mat'; 0025 0026 % Time window for dipole fit (sample number) 0027 parm.twin_meg = [100 200]; 0028 0029 % Error function mode for dipole search 0030 parm.errfunc = 1; 0031 % = 0 : for External source current (Biot-Savart Eq.) 0032 % = 1 : for Cortical source current (Sarvas Eq.) 0033 % = 2 : for External coil current (magnetic dipole) 0034 0035 % Initial dipole position parameter 0036 Xmax = 0.07; % Max coordinate for initial position [m] 0037 Xmin = -0.07; % Min coordinate for initial position [m] 0038 0039 % Initilize min error 0040 errmin = realmax; 0041 0042 if isfield(parm,'dipolefile') && ~isempty(parm.dipolefile) 0043 dipolefile = [proj_root '/' parm.dipolefile]; 0044 else 0045 dipolefile = []; 0046 end 0047 0048 parm.dipolefile = []; 0049 parm.mode = 2; 0050 0051 % Dipole search by changing initial position 0052 for n=1:Ntrial 0053 % Vinit is defined where origin is center of sphere model 0054 parm.coord_mode = 0; 0055 parm.Vinit = (Xmax - Xmin)*rand(Ndipole,3) + Xmin; 0056 0057 [X,J,err] = vb_job_dipole(proj_root, parm); 0058 0059 fprintf('Error = %f\n',err) 0060 % check min error 0061 if err < errmin, 0062 errmin = err; 0063 Vdipole = X; 0064 end; 0065 end 0066 0067 % Vinit is defined in SPM-right-m coordinate 0068 parm.coord_mode = 1; 0069 parm.Vinit = Vdipole; 0070 0071 % Final dipole search from best result 0072 [Vdipole,Jdipole,err] = vb_job_dipole(proj_root, parm); 0073 0074 for n=1:Ndipole 0075 fprintf('X = [%g , %g, %g]\n',Vdipole(n,:)); 0076 end 0077 0078 % Save result 0079 % --- Output file 0080 if ~isempty(dipolefile) 0081 vb_fsave(dipolefile, 'Vdipole','Jdipole','parm','err'); 0082 end 0083 0084 % [Vdipole,Jdipole,err,iter] = vb_job_dipole(proj_root, parm) 0085 % --- Input 0086 % proj_root : data directory 0087 % parm : parameter structure 0088 % - Input MEG file 0089 % parm.megfile = MEG file (relative path from proj_root) 0090 % - Initial dipole position 0091 % parm.Vinit = Initial dipole position [Ndipole x 3] 0092 % number of dipoles in minimum error search is given by the size of Vinit 0093 % - Error function 0094 % parm.errfunc = Error function mode for dipole search 0095 % = 0 : 'vb_dipole_error': for external current (Biot-Savart Eq.) 0096 % = 1 : 'vb_dipole_error_sarvas': for cortical current (Sarvas Eq.) 0097 % = 2 : 'vb_dipole_error_magdipole': for coil current (magnetic dipole) 0098 % - Optional parameter 0099 % parm.dipolefile = Output file (relative path from proj_root) 0100 % if this field is not exist or empty, no output file is produced. 0101 % parm.twin_meg = Time window for analysis 0102 % if this field is not exist or empty, whole time period is used. 0103 % parm.trial_id = trial number to estimate 0104 % if this field is not exist or empty, all trials are averaged. 0105 % - Optimization setting (Optional) 0106 % parm.MaxIter = 5000; % max iteration number in optimization 0107 % parm.TolX = 1.e-8; % Xの終了許容値 0108 % parm.TolFun = 1.e-10; % 関数値の終了許容値 0109 % parm.mode = 1; 0110 % = 0; % Use 'fminsearch' in Matlab 0111 % = 1; % Use 'fminunc' in Optimization Toolbox 0112 % = 2; % Use 'fminunc1' 0113 % --- Output 0114 % Vdipole : Dipole position [Ndipole x 3] 0115 % Jdipole : dipole current [Ndipole x 3 x Time] 0116 % Jdipole(n ,k ,t) = current at n-th dipole ,time 't' and k-th direction 0117 % err : normalized error 0118 % iter : iteration number