error for dipole estimation (Sarvas) err = vb_dipole_error_sarvas(V, BB, pick, Qpick, Wsensor) --- Output err : Reconstruction error --- INPUT V : dipole position (Ndipole x 3) V(n,:) : dipole position (3-D coordinate) at n-th vertex MEG covariance matrix BB = bexp*bexp'; % BB(Nsensor,Nsensor) bexp : MEG-data (Nsensor ,Time) pick(k, 1:3) : Sensor coil position : コイル位置, Qpick(k, 1:3) : Sensor coil direction : コイル方向 Wsensor : (sensor channel) x (sensor coil) : (Nsensor x Npick) basis(channel,dipole) = Wsensor * basis(coil,dipole) 2006-12-16 made by M.Sato 2007-4-3 made by M.Sato Copyright (C) 2011, ATR All Rights Reserved. License : New BSD License(see VBMEG_LICENSE.txt)
0001 function err = vb_dipole_error_sarvas(V, BB, pick, Qpick, Wsensor) 0002 % error for dipole estimation (Sarvas) 0003 % err = vb_dipole_error_sarvas(V, BB, pick, Qpick, Wsensor) 0004 % --- Output 0005 % err : Reconstruction error 0006 % 0007 % --- INPUT 0008 % V : dipole position (Ndipole x 3) 0009 % V(n,:) : dipole position (3-D coordinate) at n-th vertex 0010 % 0011 % MEG covariance matrix 0012 % BB = bexp*bexp'; % BB(Nsensor,Nsensor) 0013 % bexp : MEG-data (Nsensor ,Time) 0014 % 0015 % pick(k, 1:3) : Sensor coil position : コイル位置, 0016 % Qpick(k, 1:3) : Sensor coil direction : コイル方向 0017 % 0018 % Wsensor : (sensor channel) x (sensor coil) : (Nsensor x Npick) 0019 % 0020 % basis(channel,dipole) = Wsensor * basis(coil,dipole) 0021 % 0022 % 2006-12-16 made by M.Sato 0023 % 2007-4-3 made by M.Sato 0024 % 0025 % Copyright (C) 2011, ATR All Rights Reserved. 0026 % License : New BSD License(see VBMEG_LICENSE.txt) 0027 0028 NV = length(V(:))/3; % Number of dipoles 0029 V = reshape(V,[NV 3]); 0030 Npick = size(pick,1); % Number of sensor coil 0031 0032 %% Transform into spherical coordinate 0033 [phi,theta,r] = cart2sph(V(:,1),V(:,2),V(:,3)); 0034 0035 % Calculate unit dipole direction vector 0036 [x1,y1,z1] = sph2cart(phi+0.5*pi, zeros(NV,1), 1); 0037 [x2,y2,z2] = sph2cart( phi, theta+0.5*pi, 1); 0038 0039 vec1 = [x1,y1,z1]; 0040 vec2 = [x2,y2,z2]; 0041 %vec3 = V./repmat(sqrt(sum(V.^2, 2)), [NV 3]); 0042 0043 % calculate gradient of magnetic field by Sarvas equation 0044 % G : Lead Field (Ndipole, Nsensor) 0045 % L = 2 (電流方向) 0046 L = 2; 0047 G = zeros(NV*L , Npick); 0048 0049 for i=1:NV 0050 G(i,:) = vb_sarvas_sensor(V(i,:),vec1(i,:),pick, Qpick)'; 0051 G(i+NV,:) = vb_sarvas_sensor(V(i,:),vec2(i,:),pick, Qpick)'; 0052 end 0053 0054 G = G * Wsensor'; 0055 0056 GGinv = pinv( G*G' ); % (NV*L x NV*L) matrix inverse 0057 0058 % Reconstruction error 0059 % Error = ( B - G'*J )^2 = B^2 - ( B'* G'* GGinv * G * B ) 0060 Bnorm = sum(diag(BB)); 0061 err = 1 - sum(sum(G.*(GGinv*G*BB)))/Bnorm; 0062