ip_dipolefitting.m

function [dipoles_str, session] = ip_dipolefitting(EEG, eloc, subject_name, session_name, comp_index, constr, warpback, varargin)

% Usage:
%   >> dipoles_str = ip_dipolefitting(EEG, sensor_file, subject_name,
%   session_name, comp_index, warpback);
%
% Inputs:
%   EEG - EEGLAB data structure
%   sensor_file - sensor file name
%   subject_name - subject name
%   session_name - session name
%   comp_index - component indices
%   warpback - warping structures
%
% Optional keywords:
%
%   sensor_name:  sensor name (default: [subject_name '_' session_name '.sensors'])
%   ss_name : sourcespace name (default: [subject_name '_sourcespace.dip'])
%
% Outputs:
%   dipoles_str - dipole structure as in dipfit

% Author: Zeynep Akalin Acar, SCCN, 2008

% Copyright (C) 2007 Zeynep Akalin Acar, SCCN, zeynep@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


% session_name should be entered with its folder
% model_name should be entered with its folder
% if flag1 == 1 mri_based, no warping back

% for old type of sensor structure, this may be removed in future
%eloc = readlocs(sensor_file);

% load sensor index

sensor_name = [subject_name '_' session_name '.sensors'];
ss_name = [subject_name '_sourcespace.dip'];


for i = 1:2:length(varargin) % for each Keyword
      Keyword = varargin{i};
      Value = varargin{i+1};

      if ~isstr(Keyword)
         fprintf('keywords must be strings')
         return
      end

      if strcmp(Keyword,'sensor_name')
         if ~isstr(Value)
            fprintf('sensor_name must be a string');
            return
         else
             sensor_name = Value;
         end
      elseif strcmp(Keyword,'ss_name')
         if ~isstr(Value)
            fprintf('ss_name must be a string');
            return
         else
             ss_name = Value;
         end
      
      end
end

a = dir(sensor_name);
if size(a,1) == 0
    b = [subject_name '_' session_name '_sensorindex'];
    if size(b,1) == 1
        ind_n = [subject_name '_' session_name '_sensorindex'];
        sens_index = load(ind_n);
    else
        error('sensor file is not found!')
    end
else
    se = load(sensor_name,'-mat');
    sens_index = se.ind;
end

% load LFM
load([session_name '_LFM']);
avg = mean(LFM, 1);    LFM = LFM - repmat(avg, size(LFM,1), 1);

disp('Dipole fitting...')

[A2, session, LFM2, ind_fp, ind_eeg, elocn] = eloc2eeglab_r(EEG, session_name, LFM, sens_index, eloc);
%elocn2 = elocn(ind_fp); % yeni warp edilmis electrod noktalari
% for the inverse problem solution
ss = load(ss_name);
spacing = max(abs(ss(1,1:3)-ss(2,1:3)))/2;

% check session type
if ~isfield(session, 'type') ||  ~strcmp(session.type, 'fem')
    % BEM
    if session.model.mod>0
        if ~isfield(session.model, 'iinv')
            session.model = bem_load_model_matrix(session.model,'iinv');
        end
        if ~isfield(session.model, 'dmt')
            session.model = bem_load_model_matrix(session.model,'dmt');
        end
    end
    %load vol 
    [vol, sens] = session2vol(session);
    sens.label = sens.label(ind_fp);
    sens.pnt = sens.pnt(ind_fp,:);
else
    %FEM
    vol = session.vol;
    sens = session.sens;
    sens.label='A1';
   
    metufem('setup',session.vol.mesh_name,'','')
    metufem('setrf',session.sens.rf)
end

%if isfield(session.model.mesh,'transform')
%    if length(session.model.mesh.transform) == 3
%        ss(:,1:3) = ss(:,1:3) - ones(size(ss,1),1) * session.model.mesh.transform;
%    end
%end
%tr=session.model.mesh.transform;

for compi = comp_index
    compi
    Vdata = A2(:, compi);
    Vdata = Vdata - mean(Vdata);
    [pos_bin, griderror] = Grid_dipole2(LFM2, Vdata, ss);
    [~,I1] = sort(griderror);
    pos_bin_grid = pos_bin(I1(1),:);
    
    % Calculate pos_bin_grid on a finer grid for FEM
    if isfield(session, 'type')
        if strcmp(session.type, 'fem')
            % compute on fine grid 
            ss1 = new_ss(pos_bin_grid, spacing);
            LFMp = metufem('pot', ss1','interp');
            avg = mean(LFMp, 1);    LFMp = LFMp - repmat(avg, size(LFMp,1), 1);
            [pos_bin, griderror] = Grid_dipole2(LFMp, Vdata, ss1);
            [Y1,I1] = sort(griderror);
            pos_bin_grid = pos_bin(I1(1),:);
        end
    end
    [dip_bin] = dipole_fit(struct('pos',pos_bin_grid), sens, vol, Vdata, 'constr', constr,'maxiter',500);

    % for multiple dipole source localization - necessary for bilateral dipoles !!!!
    %for i=1:size(dip_bin.pos,1)
    %    dip_bin.mom(i,:) = dip_bin.mom((i-1)*3+1:i*3);
    %end
    %dip_bin.mom = mom;

    lf = ft_compute_leadfield(dip_bin.pos, sens, vol);
    mom = pinv(lf)*Vdata;
    dif = Vdata - lf*mom;
    
    % relative residual variance
    fval = sum(dif(:).^2) / sum(Vdata(:).^2);

    
    if length(warpback) > 0 % warp back dipole locations
        % warp back
        for i = 1:size(dip_bin.pos,1)
            dip_bin.pos(i,:) = dip_bin.pos(i,:) + warp_lm(dip_bin.pos(i,:), warpback.A, warpback.W, warpback.LMd);
        end
    end
    
    
    % check this !!! (before or after warping???)
    %if isfield(session.model.mesh,'transform')
    %    if length(session.model.mesh.transform) == 3
    %        dip_bin.pos = dip_bin.pos + ones(size(dip_bin.pos,1),1)*session.model.mesh.transform;
    %    end
    %end
    
    dipoles_str(compi).posxyz = dip_bin.pos;
    dipoles_str(compi).momxyz = dip_bin.mom';
    dipoles_str(compi).rv = fval;
end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [pos,griderror] = Grid_dipole2(lf, Vdata, BrSS)
% compute grid error of individual dipole-sets in the LFM
% The LFM consists of sets of three orthagonal dipoles occupying
% the same location. 

% Zeynep Akalin Acar, 2008

nchans = length(Vdata);
ndip = size(lf,2) / 3;
m = size(BrSS,2)-6;
griderror = zeros(ndip,1);
pos = zeros(ndip,3);

s = 1;
for d = 1:ndip;
    i = [d d+ndip d+ndip+ndip];
    griderror(d) = sum(((eye(nchans)-lf(:,i)*pinv(lf(:,i)))*Vdata).^2);
    pos(d,:) = BrSS(d,1+m:3+m);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [pos,griderror] = Grid_dipole2sym(lf, Vdata, BrSS)
% compute grid error of individual dipole-sets in the symmetric LFM
% The symmetric LFM consists of sets of six orthagonal dipoles
% for the same symmetric pair. 

% Zeynep Akalin Acar, 2011

nchans = length(Vdata);
ndip = size(lf,2) / 6;
m = size(BrSS,2)-6;
griderror = zeros(ndip,1);
pos = zeros(ndip,3);

s = 1;
for d = 1:ndip;
    i = [d d+ndip d+2*ndip d+3*ndip d+4*ndip d+5*ndip];
    griderror(d) = sum(((eye(nchans)-lf(:,i)*pinv(lf(:,i)))*Vdata).^2);
    pos(d,:) = BrSS(d,1+m:3+m);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rw] = warp_lm(r,A,W,p)
% performs warp transformation with linear 3D RFB see Ermer's Thesis
rw = r * A(1:3,1:3) + repmat(A(4,:), size(r,1), 1);
for i = 1 : size(p,1)
    U = sqrt(sum((r - repmat(p(i,:), size(r,1),1)).^2, 2));  
    rw = rw + U * W(i,:);
end    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [A2, session, LFM2, ind_fp, ind_eeg, elocn] = eloc2eeglab_r(EEG, session_name, LFM, elp_index, eloc);
% realistic data icin

A = EEG.icawinv;

% neglect the FID electrodes
y = elp_index;
N = length(y);
for i = 1:N
    elocn(i).labels = eloc(y(i)).labels;
    elocn(i).X = eloc(y(i)).X;
    elocn(i).Y = eloc(y(i)).Y;
    elocn(i).Z = eloc(y(i)).Z;
 %   elocn(i).type = eloc(y(i)).type;
end

Nel2 = length(elocn);
Neeg = length(EEG.chanlocs);

Mel2 = zeros(1, Nel2);
Meeg = zeros(1, Neeg);

% Mel2(i) = index of EEG.chanlocs that correspond to electrode i of eloc2
% Meeg(i) = index of elocn that correspond to electrode i of EEG.chanlocs
clear Mel2 Meeg
for i = 1:Nel2
    for j = 1:Neeg
        if strcmp(elocn(i).labels, EEG.chanlocs(j).labels)
            Mel2(i) = j;
            Meeg(j) = i;
            continue;
        end
    end
end

ind_fp = find(Mel2>0); % index for the FP outputs (LFM, TM, session)
ind_eeg = find(Meeg>0); % index for the EEG structure (ICs)

if ~isfield(EEG.etc,'nft')
    EEG.etc.nft=[];
end
if ~isfield(EEG.etc.nft,'session')
    % load the session and the transfer matrix
    ssave = load([session_name '.session'], '-MAT');
    if isfield(ssave, 'model_name') 
        % BEM
        model = load_model([ssave.model_name, '.model']);
        session = bem_create_session(ssave.name, model, ssave.Smatrix);
        session = bem_load_transfer_matrix(session, 'tmte');
        % update the transfer matrix of the session wrt ind_fp
        session.tmte = session.tmte(:,ind_fp);
        session.num_electrodes = length(ind_fp);
        if session.model.mod>0
            session.model = bem_load_model_matrix(session.model,'iinv');
            session.model = bem_load_model_matrix(session.model,'dmt');
        end
   else
        session = ssave.session;
        session.sens.rf = session.sens.rf(:,ind_fp);
        session.sens.pnt = session.sens.pnt(ind_fp,:);    
    end
else
    session = EEG.etc.nft.session;
end




% update LFM wrt. ind_fp
LFM2 = LFM(ind_fp,:);

% update EEG.icawinv wrt. ind_eeg

A2 = A(ind_eeg,:);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function model = load_model(file)
msave = load(file, '-MAT');
mesh = bem_load_mesh(msave.mesh_name);
model = bem_create_model(msave.name, mesh, msave.cond, msave.mod);
session = [];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ss1 = new_ss(pos_bin_grid, gm)
% find a fine grid 1with 1 mm spacing +-5mm around pos_bin_grid

spacing = 1; % 1 mm spacing
%gm = 4; 
ma = pos_bin_grid + [gm gm gm];
mi = pos_bin_grid - [gm gm gm];
md = ((ma - mi)/spacing);

n = 0;
for i = 1:md(1)+1
    a(i) = mi(1)+spacing*(i-1);
    for j = 1:md(2)+1
        b(j) = mi(2)+spacing*(j-1);
        for k = 1:md(3)+1
            c(k) = mi(3)+spacing*(k-1);
            n = n+1;
            so(n,1:3) = [a(i) b(j) c(k)];
        end
    end
end
ne = size(so,1);    
ss1 = zeros(ne*3,6);
ss1(1:ne,1:3) = so;
ss1(ne+1:2*ne,1:3) = so;
ss1(2*ne+1:3*ne,1:3) = so;
ss1(1:ne,4) = 1;
ss1(ne+1:2*ne,5) = 1;
ss1(2*ne+1:3*ne,6) = 1;