vb_dipole_error_magdipole

 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