vb_dipole_current

 Dipole current time cource for dipoles at fixed position
  [J, err, B] = vb_dipole_current(V, bexp, pick, Qpick, Wsensor, V0, mode)
 --- Output
  J   : dipole current [NV, 3 ,Time]
  J(n ,k ,t) =  current at n-th dipole ,time 't' and k-th direction
 err  : Reconstruction error
 B(m,t) : reconstructed magnetic field at m-th channel and time 't' 
 --- INPUT
 V      : position of dipoles (NV x 3)
 V(n,:) : n-th dipole position (3-D coordinate)

 bexp(m,t)   : 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)
 --- Optional input
 V0   : Center of spherical model (1 x 3)
 mode = 0: Biot-Savart      , L = 3: independent current direction
      = 1: Sarvas           , L = 2: independent current direction
      = 2: Magnetic dipole  , L = 3: independent current direction
      = 3: EEG              , L = 3: independent current direction

 2006-12-16 made by M.Sato
 2007-4-3 made by M.Sato
 2007-7-9 modified by M.Sato

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function  [J, err, B] = vb_dipole_current(V, bexp, pick, Qpick, Wsensor, V0, mode)
0002 % Dipole current time cource for dipoles at fixed position
0003 %  [J, err, B] = vb_dipole_current(V, bexp, pick, Qpick, Wsensor, V0, mode)
0004 % --- Output
0005 %  J   : dipole current [NV, 3 ,Time]
0006 %  J(n ,k ,t) =  current at n-th dipole ,time 't' and k-th direction
0007 % err  : Reconstruction error
0008 % B(m,t) : reconstructed magnetic field at m-th channel and time 't'
0009 % --- INPUT
0010 % V      : position of dipoles (NV x 3)
0011 % V(n,:) : n-th dipole position (3-D coordinate)
0012 %
0013 % bexp(m,t)   : 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 % --- Optional input
0021 % V0   : Center of spherical model (1 x 3)
0022 % mode = 0: Biot-Savart      , L = 3: independent current direction
0023 %      = 1: Sarvas           , L = 2: independent current direction
0024 %      = 2: Magnetic dipole  , L = 3: independent current direction
0025 %      = 3: EEG              , L = 3: independent current direction
0026 %
0027 % 2006-12-16 made by M.Sato
0028 % 2007-4-3 made by M.Sato
0029 % 2007-7-9 modified by M.Sato
0030 %
0031 % Copyright (C) 2011, ATR All Rights Reserved.
0032 % License : New BSD License(see VBMEG_LICENSE.txt)
0033 
0034 if ~exist('V0','var')|isempty(V0), V0 = [0 0 0]; end;
0035 if ~exist('mode','var'), mode = 0; end;
0036 
0037 T  = size(bexp,2);
0038 NV = length(V(:))/3;        % Number of dipole position
0039 V  = reshape(V,[NV 3]);
0040 
0041 switch    mode
0042 case    0
0043     % Magnetic field by Biot-Savart
0044     L  = 3;
0045 case    1
0046     % Magnetic field by Sarvas
0047     L  = 2;
0048 case    2
0049     % Magnetic field by magnetic dipole
0050     L  = 3;
0051 case    3
0052     % EEG 3-shell Sphere model
0053     L  = 3;
0054 end
0055 
0056 % L = 2 :
0057 %   xx(n,:)      = 1st orthogonal vectors to V(n,:)
0058 %   xx(n + NV,:) = 2nd orthogonal vectors to V(n,:)
0059 % L = 3 :
0060 %    xx(n       ,:) =  [1 0 0]
0061 %    xx(n + NV  ,:) =  [0 1 0]
0062 %    xx(n + 2*NV,:) =  [0 0 1]
0063 
0064 % xx : current direction
0065 [V,xx] = vb_current_vector(V,[],L);
0066 
0067 % G : Lead Field (NV*L, Nsensor)
0068 G = vb_dipole_basis(V, xx, pick, Qpick, Wsensor, V0, mode);
0069 
0070 GGinv = pinv( G*G' ); % (NV*L x NV*L) matrix inverse
0071 
0072 % Diploe current
0073 JJ = (GGinv * G) * bexp;    % (NV*L x T)
0074 
0075 % reconstructed field
0076 B  = G' * JJ;
0077 
0078 JJ = reshape(JJ, [NV, L, T]);
0079 
0080 if L == 2,
0081     J  = zeros(NV,3,T);
0082     ix = 1:NV;
0083     
0084     for k = 1:3
0085         J(:,k,:) = JJ(ix,1,:) .* repmat(xx(ix,k), [1, 1, T]) ...
0086                  + JJ(ix,2,:) .* repmat(xx(ix + NV,k), [1, 1, T]) ;
0087     end
0088 else
0089     J = JJ;
0090 end
0091 
0092 if nargout == 1, return; end;
0093 
0094 % Reconstruction error
0095 
0096 % MEG covariance matrix
0097 BB    = bexp*bexp';    % BB(Nsensor,Nsensor)
0098 Bnorm = sum(diag(BB));
0099 
0100 % Error = ( B - G'*J )^2 = B^2  - ( B'* G'* GGinv * G * B )
0101 err  = 1 - sum(sum(G.*(GGinv*G*BB)))/Bnorm;
0102 
0103 return
0104 
0105 % --- DEBUG
0106 
0107 err2  = sum(sum((bexp-B).^2))/Bnorm;
0108 delta = err - err2