Error for dipole estimation (Biot-Savart)
err = vb_dipole_error(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(m,n) = n-th coil weight for m-th sensor channel
basis(channel,dipole) = Wsensor * basis(coil,dipole)
双極子磁場計算
x0 の双極子 J0 が 観測点 X に作る磁場の Q 方向射影
B0 = Q*(JJ x Xd)/|Xd|^3 , Xd = X - x0
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(V, BB, pick, Qpick, Wsensor) 0002 % Error for dipole estimation (Biot-Savart) 0003 % err = vb_dipole_error(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(m,n) = n-th coil weight for m-th sensor channel 0019 % basis(channel,dipole) = Wsensor * basis(coil,dipole) 0020 % 0021 % 双極子磁場計算 0022 % x0 の双極子 J0 が 観測点 X に作る磁場の Q 方向射影 0023 % B0 = Q*(JJ x Xd)/|Xd|^3 , Xd = X - x0 0024 % 0025 % 2006-12-16 made by M.Sato 0026 % 2007-4-3 made by M.Sato 0027 % 0028 % Copyright (C) 2011, ATR All Rights Reserved. 0029 % License : New BSD License(see VBMEG_LICENSE.txt) 0030 0031 0032 NV = length(V(:))/3; % Number of dipoles 0033 V = reshape(V,[NV 3]); 0034 0035 % 3 current directions 0036 vec1 = zeros(NV,3); vec1(:,1) = 1; 0037 vec2 = zeros(NV,3); vec2(:,2) = 1; 0038 vec3 = zeros(NV,3); vec3(:,3) = 1; 0039 0040 % magnetic field by Biot-Savart 0041 G1 = vb_dipole_magnetic(V, vec1, pick, Qpick ); 0042 G2 = vb_dipole_magnetic(V, vec2, pick, Qpick ); 0043 G3 = vb_dipole_magnetic(V, vec3, pick, Qpick ); 0044 0045 % G : Lead Field (NV*L, Nsensor) 0046 G = [G1; G2; G3]* Wsensor'; 0047 0048 Bnorm = sum(diag(BB)); 0049 0050 GGinv = pinv( G*G' ); % (NV*L x NV*L) matrix inverse 0051 0052 % Reconstruction error 0053 % Error = ( B - G'*J )^2 = B^2 - ( B'* G'* GGinv * G * B ) 0054 err = 1 - sum(sum(G.*(GGinv*G*BB)))/Bnorm; 0055