Error for dipole estimation (magnetic dipole model)
(Note: it is different from current dipole)
err = vb_dipole_error_magdipole(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)
Magnetic field produced by magnetic dipole
Coil current can be approximated by magnetic dipole.
磁気双極子磁場 (コイルに流れる電流の近似式)
x0 の磁気双極子 J0 が 観測点 X に作る磁場の Q 方向射影
B = 3(Q*Xd)(J0*Xd)/|Xd|^5 - (Q*J0)/|Xd|^3 , Xd = X - x0
2007-7-9 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_magdipole(V, BB, pick, Qpick, Wsensor) 0002 % Error for dipole estimation (magnetic dipole model) 0003 % (Note: it is different from current dipole) 0004 % err = vb_dipole_error_magdipole(V, BB, pick, Qpick, Wsensor) 0005 % --- Output 0006 % err : Reconstruction error 0007 % 0008 % --- INPUT 0009 % V : dipole position (Ndipole x 3) 0010 % V(n,:) : dipole position (3-D coordinate) at n-th vertex 0011 % 0012 % MEG covariance matrix 0013 % BB = bexp*bexp'; % BB(Nsensor,Nsensor) 0014 % bexp : MEG-data (Nsensor ,Time) 0015 % 0016 % pick(k, 1:3) : Sensor coil position : コイル位置, 0017 % Qpick(k, 1:3) : Sensor coil direction : コイル方向 0018 % 0019 % Wsensor(m,n) = n-th coil weight for m-th sensor channel 0020 % basis(channel,dipole) = Wsensor * basis(coil,dipole) 0021 % 0022 % Magnetic field produced by magnetic dipole 0023 % Coil current can be approximated by magnetic dipole. 0024 % 0025 % 磁気双極子磁場 (コイルに流れる電流の近似式) 0026 % x0 の磁気双極子 J0 が 観測点 X に作る磁場の Q 方向射影 0027 % 0028 % B = 3(Q*Xd)(J0*Xd)/|Xd|^5 - (Q*J0)/|Xd|^3 , Xd = X - x0 0029 % 0030 % 2007-7-9 made by M.Sato 0031 % 0032 % Copyright (C) 2011, ATR All Rights Reserved. 0033 % License : New BSD License(see VBMEG_LICENSE.txt) 0034 0035 0036 NV = length(V(:))/3; % Number of dipoles 0037 V = reshape(V,[NV 3]); 0038 0039 % 3 current directions 0040 vec1 = zeros(NV,3); vec1(:,1) = 1; 0041 vec2 = zeros(NV,3); vec2(:,2) = 1; 0042 vec3 = zeros(NV,3); vec3(:,3) = 1; 0043 0044 % magnetic field by Biot-Savart 0045 G1 = vb_dipole_field_magdipole(V, vec1, pick, Qpick ); 0046 G2 = vb_dipole_field_magdipole(V, vec2, pick, Qpick ); 0047 G3 = vb_dipole_field_magdipole(V, vec3, pick, Qpick ); 0048 0049 % G : Lead Field (NV*L, Nsensor) 0050 G = [G1; G2; G3]* Wsensor'; 0051 0052 Bnorm = sum(diag(BB)); 0053 0054 GGinv = pinv( G*G' ); % (NV*L x NV*L) matrix inverse 0055 0056 % Reconstruction error 0057 % Error = ( B - G'*J )^2 = B^2 - ( B'* G'* GGinv * G * B ) 0058 err = 1 - sum(sum(G.*(GGinv*G*BB)))/Bnorm; 0059