vb_dipole_error

 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