eeg_interp.m

% EEG_INTERP - interpolate data channels
%
% Usage: EEGOUT = eeg_interp(EEG, badchans, method, t_range);
%
% Inputs: 
%     EEG      - EEGLAB dataset
%     badchans - [integer array] indices of channels to interpolate.
%                For instance, these channels might be bad.
%                [chanlocs structure] channel location structure containing
%                either locations of channels to interpolate or a full
%                channel structure (missing channels in the current 
%                dataset are interpolated).
%     method   - [string] method used for interpolation (default is 'spherical').
%                'invdist'/'v4' uses inverse distance on the scalp
%                'spherical' uses superfast spherical interpolation. 
%                'sphericalKang' uses Kang et al. 2015 parameters 
%                'spacetime' uses griddata3 to interpolate both in space 
%                and time (very slow and cannot be interrupted).
%     t_range  - [integer array with just two elements] time interval of the
%                badchans which should be interpolated. First element is
%                the start time and the second element is the end time.
%     params   - [1x3 array] Parameters [lambda m n]. Defaults are
%                'spherical'     -> lambda=0, m = 4, n =7 (Perin et al., 1989)
%                'sphericalKang' -> lambda=1e-8, m = 3 , n = 50 (Kang et al., 2015)
%                if the 'params' input is used, then a spherical
%                interpolation is used with these parameters (method is
%                ignored). See also https://github.com/sccn/eeglab/issues/500
% Output: 
%     EEGOUT   - data set with bad electrode data replaced by
%                interpolated data
%
% Note: See also the sphericalSplineInterpolate function of the clean_rawdata
%       EEGLAB plugin (in the private folder) which approximate Legendre 
%       polynomials with up to 500 iteration.                  .
%
% References:
% Perrin, F., Pernier, J., Bertrand, O., & Echallier, J. F. (1989). Spherical 
% splines for scalp potential and current density mapping. Electroencephalography 
% and clinical neurophysiology, 72(2), 184–187. https://doi.org/10.1016/0013-4694(89)90180-6
%
% Kang, S.S., Lano, T.J. & Sponheim, S.R. Distortions in EEG interregional phase 
% synchrony by spherical spline interpolation: causes and remedies. Neuropsychiatr 
% Electrophysiol 1, 9 (2015). https://doi.org/10.1186/s40810-015-0009-5
%
% Author: Arnaud Delorme, CERCO, CNRS, Mai 2006-

% Copyright (C) Arnaud Delorme, CERCO, 2006, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function EEG = eeg_interp(ORIEEG, bad_elec, method, t_range, params)

    if nargin < 2
        help eeg_interp;
        return;
    end
    EEG = ORIEEG;
    
    if nargin < 3
        disp('Using spherical interpolation');
        method = 'spherical';
    end
    
    if nargin < 4
        t_range = [ORIEEG.xmin ORIEEG.xmax];
    end
    if nargin < 5
        if strcmpi(method, 'spherical')
            params = [0 4 7];
        elseif strcmpi(method, 'sphericalKang')
            params = [1e-8 3 50];
        elseif strcmpi(method, 'sphericalCRD')
            params = [1e-5 4 500];
        end
    else
        if length(params) ~=3 
            error('parameters array must have 3 elements')
        end
        method = 'spherical';
    end

    % check channel structure
    tmplocs = ORIEEG.chanlocs;
    if isempty(tmplocs) || isempty([tmplocs.X])
        error('Interpolation require channel location');
    end
    
    if isstruct(bad_elec)
       
        % remove duplicate channel labels
        % -------------------------------
        allLabels = { bad_elec.labels };
        if length(unique(allLabels)) ~= length(allLabels)
            [~,uniqueInd] = unique(allLabels);
            bad_elec = bad_elec(uniqueInd);
        end

        % add missing channels in interpolation structure
        % -----------------------------------------------
        lab1 = { bad_elec.labels };
        tmpchanlocs = EEG.chanlocs;
        lab2 = { tmpchanlocs.labels };
        [~, tmpchan] = setdiff_bc( lab2, lab1);
        tmpchan = sort(tmpchan);
        
        % From 'bad_elec' using only fields present on EEG.chanlocs
        fields = fieldnames(bad_elec);
        [~, indx1] = setxor(fields,fieldnames(EEG.chanlocs)); clear tmp;
        if ~isempty(indx1)
            bad_elec = rmfield(bad_elec,fields(indx1));
            fields = fieldnames(bad_elec);
        end
        
        if ~isempty(tmpchan)
            newchanlocs = [];
            for index = 1:length(fields)
                if isfield(bad_elec, fields{index})
                    for cind = 1:length(tmpchan)
                        fieldval = getfield( EEG.chanlocs, { tmpchan(cind) },  fields{index});
                        newchanlocs = setfield(newchanlocs, { cind }, fields{index}, fieldval);
                    end
                end
            end
            newchanlocs(end+1:end+length(bad_elec)) = bad_elec;
            bad_elec = newchanlocs;
        end
        if length(EEG.chanlocs) == length(bad_elec), return; end
        
        lab1 = { bad_elec.labels };
        tmpchanlocs = EEG.chanlocs;
        lab2 = { tmpchanlocs.labels };
        [~, badchans] = setdiff_bc( lab1, lab2);
        fprintf('Interpolating %d channels...\n', length(badchans));
        if isempty(badchans), return; end
        goodchans      = sort(setdiff(1:length(bad_elec), badchans));
       
        % re-order good channels
        % ----------------------
        [~, tmp2, neworder] = intersect_bc( lab1, lab2 );
        [~, ordertmp2] = sort(tmp2);
        neworder = neworder(ordertmp2);
        EEG.data = EEG.data(neworder, :, :);

        % looking at channels for ICA
        % ---------------------------
        %[tmp sorti] = sort(neworder);
        %{ EEG.chanlocs(EEG.icachansind).labels; bad_elec(goodchans(sorti(EEG.icachansind))).labels }
        
        % update EEG dataset (add blank channels)
        % ---------------------------------------
        if ~isempty(EEG.icasphere)
            
            [~, sorti] = sort(neworder);
            EEG.icachansind = sorti(EEG.icachansind);
            EEG.icachansind = goodchans(EEG.icachansind);
            EEG.chaninfo.icachansind = EEG.icachansind;
            
            % TESTING SORTING
            %icachansind = [ 3 4 5 7 8]
            %data = round(rand(8,10)*10)
            %neworder = shuffle(1:8)
            %data2 = data(neworder,:)
            %icachansind2 = sorti(icachansind)
            %data(icachansind,:)
            %data2(icachansind2,:)
        end
        % { EEG.chanlocs(neworder).labels; bad_elec(sort(goodchans)).labels }
        %tmpdata                  = zeros(length(bad_elec), size(EEG.data,2), size(EEG.data,3));
        %tmpdata(goodchans, :, :) = EEG.data;
        
        % looking at the data
        % -------------------
        %tmp1 = mattocell(EEG.data(sorti,1));
        %tmp2 = mattocell(tmpdata(goodchans,1));
        %{ EEG.chanlocs.labels; bad_elec(goodchans).labels; tmp1{:}; tmp2{:} }
        %EEG.data      = tmpdata;
        
        EEG.chanlocs  = bad_elec;

    else
        badchans  = bad_elec;
        goodchans = setdiff_bc(1:EEG.nbchan, badchans);
        if strcmpi(method, 'sphericalfast')
            EEG.data(badchans,:) = [];
            EEG.nbchan = length(goodchans);
        else
            oldelocs  = EEG.chanlocs;
            EEG       = pop_select(EEG, 'nochannel', badchans);
            EEG.chanlocs = oldelocs;
            disp('Interpolating missing channels...');
        end
    end

    % find non-empty good channels
    % ----------------------------
    origoodchans = goodchans;
    chanlocs     = EEG.chanlocs;
    nonemptychans = find(~cellfun('isempty', { chanlocs.theta }));
    [tmp indgood ] = intersect_bc(goodchans, nonemptychans);
    goodchans = goodchans( sort(indgood) );
    datachans = getdatachans(goodchans,badchans);
    badchans  = intersect_bc(badchans, nonemptychans);
    if isempty(badchans), return; end
    
    % scan data points
    % ----------------
    if strcmpi(method, 'spherical') || strcmpi(method, 'sphericalfast') || strcmpi(method, 'sphericalKang')
        % get theta, rad of electrodes
        % ----------------------------
        tmpgoodlocs = EEG.chanlocs(goodchans);
        xelec = [ tmpgoodlocs.X ];
        yelec = [ tmpgoodlocs.Y ];
        zelec = [ tmpgoodlocs.Z ];
        rad = sqrt(xelec.^2+yelec.^2+zelec.^2);
        xelec = xelec./rad;
        yelec = yelec./rad;
        zelec = zelec./rad;
        tmpbadlocs = EEG.chanlocs(badchans);
        xbad = [ tmpbadlocs.X ];
        ybad = [ tmpbadlocs.Y ];
        zbad = [ tmpbadlocs.Z ];
        rad = sqrt(xbad.^2+ybad.^2+zbad.^2);
        xbad = xbad./rad;
        ybad = ybad./rad;
        zbad = zbad./rad;
        if isempty(xbad)
            error('Trying to interpolate a channel without coordinates assigned to it');
        end

        EEG.data = reshape(EEG.data, EEG.nbchan, EEG.pnts*EEG.trials);
        %[tmp1 tmp2 tmp3 tmpchans] = spheric_spline_old( xelec, yelec, zelec, EEG.data(goodchans,1));
        %max(tmpchans(:,1)), std(tmpchans(:,1)), 
        %[tmp1 tmp2 tmp3 EEG.data(badchans,:)] = spheric_spline( xelec, yelec, zelec, xbad, ybad, zbad, EEG.data(goodchans,:));
        [tmp1 tmp2 tmp3 badchansdata] = spheric_spline( xelec, yelec, zelec, xbad, ybad, zbad, EEG.data(datachans,:), params);
        %max(EEG.data(goodchans,1)), std(EEG.data(goodchans,1))
        %max(EEG.data(badchans,1)), std(EEG.data(badchans,1))
        EEG.data = reshape(EEG.data, EEG.nbchan, EEG.pnts, EEG.trials);
    elseif strcmpi(method, 'spacetime') % 3D interpolation, works but x10 times slower
        disp('Warning: if processing epoch data, epoch boundary are ignored...');
        disp('3-D interpolation, this can take a long (long) time...');
        tmpgoodlocs = EEG.chanlocs(goodchans);
        tmpbadlocs = EEG.chanlocs(badchans);
        [xbad ,ybad]  = pol2cart([tmpbadlocs.theta],[tmpbadlocs.radius]);
        [xgood,ygood] = pol2cart([tmpgoodlocs.theta],[tmpgoodlocs.radius]);
        pnts = size(EEG.data,2)*size(EEG.data,3);
        zgood = 1:pnts;
        zgood = repmat(zgood, [length(xgood) 1]);    
        zgood = reshape(zgood,numel(zgood),1);
        xgood = repmat(xgood, [1 pnts]); xgood = reshape(xgood,numel(xgood),1);
        ygood = repmat(ygood, [1 pnts]); ygood = reshape(ygood,numel(ygood),1);
        tmpdata = reshape(EEG.data(datachans,:,:), numel(EEG.data(datachans,:,:)),1);
        zbad = 1:pnts;
        zbad = repmat(zbad, [length(xbad) 1]);     
        zbad = reshape(zbad,numel(zbad),1);
        xbad = repmat(xbad, [1 pnts]); xbad = reshape(xbad,numel(xbad),1);
        ybad = repmat(ybad, [1 pnts]); ybad = reshape(ybad,numel(ybad),1);
        badchansdata = griddatan([ygood, xgood, zgood], tmpdata,[ybad, xbad, zbad], 'nearest'); % interpolate data                                            
    else 
        % get theta, rad of electrodes
        % ----------------------------
        tmpchanlocs = EEG.chanlocs;
        [xbad ,ybad]  = pol2cart([tmpchanlocs( badchans).theta],[tmpchanlocs( badchans).radius]);
        [xgood,ygood] = pol2cart([tmpchanlocs(goodchans).theta],[tmpchanlocs(goodchans).radius]);

        fprintf('Points (/%d):', size(EEG.data,2)*size(EEG.data,3));
        badchansdata = zeros(length(badchans), size(EEG.data,2)*size(EEG.data,3));
        
        for t=1:(size(EEG.data,2)*size(EEG.data,3)) % scan data points
            if mod(t,100) == 0, fprintf('%d ', t); end
            if mod(t,1000) == 0, fprintf('\n'); end;          
        
            %for c = 1:length(badchans)
            %   [h EEG.data(badchans(c),t)]= topoplot(EEG.data(goodchans,t),EEG.chanlocs(goodchans),'noplot', ...
            %        [EEG.chanlocs( badchans(c)).radius EEG.chanlocs( badchans(c)).theta]);
            %end
            tmpdata = reshape(EEG.data, size(EEG.data,1), size(EEG.data,2)*size(EEG.data,3) );
            if strcmpi(method, 'invdist'), method = 'v4'; end
            [Xi,Yi,badchansdata(:,t)] = griddata(ygood, xgood , double(tmpdata(datachans,t)'),...
                                                    ybad, xbad, method); % interpolate data                                            
        end
        fprintf('\n');
    end
    
    tmpdata               = zeros(length(bad_elec), EEG.pnts, EEG.trials);
    tmpdata(origoodchans, :,:) = EEG.data;
    %if input data are epoched reshape badchansdata for Octave compatibility...
    if length(size(tmpdata))==3
        badchansdata = reshape(badchansdata,length(badchans),size(tmpdata,2),size(tmpdata,3));
    end

    if ~isempty(t_range)
        t_start = t_range(1);
        t_end   = t_range(2);
        if t_start~=ORIEEG.xmin || t_end~=ORIEEG.xmax
            if length(size(tmpdata))==2
                times_2b_ignored = [1:floor(t_start*ORIEEG.srate), floor(t_end*ORIEEG.srate):floor(ORIEEG.xmax*ORIEEG.srate)];
                badchansdata(:, times_2b_ignored) = ORIEEG.data(badchans, times_2b_ignored);
            end
        end
    end

    tmpdata(badchans,:,:) = badchansdata;
    EEG.data = tmpdata;
    EEG.nbchan = size(EEG.data,1);
    if ~strcmpi(method, 'sphericalfast')
        EEG = eeg_checkset(EEG);
    end
    
% get data channels
% -----------------
function datachans = getdatachans(goodchans, badchans);
      datachans = goodchans;
      badchans  = sort(badchans);
      for index = length(badchans):-1:1
          datachans(find(datachans > badchans(index))) = datachans(find(datachans > badchans(index)))-1;
      end
        
% -----------------
% spherical splines
% -----------------
% function [x, y, z, Res] = spheric_spline_old( xelec, yelec, zelec, values)
% 
% SPHERERES = 20;
% [x,y,z] = sphere(SPHERERES);
% x(1:(length(x)-1)/2,:) = []; x = [ x(:)' ];
% y(1:(length(y)-1)/2,:) = []; y = [ y(:)' ];
% z(1:(length(z)-1)/2,:) = []; z = [ z(:)' ];
% 
% Gelec = computeg(xelec,yelec,zelec,xelec,yelec,zelec);
% Gsph  = computeg(x,y,z,xelec,yelec,zelec);
% 
% % equations are 
% % Gelec*C + C0  = Potential (C unknown)
% % Sum(c_i) = 0
% % so 
% %             [c_1]
% %      *      [c_2]
% %             [c_ ]
% %    xelec    [c_n]
% % [x x x x x]         [potential_1]
% % [x x x x x]         [potential_ ]
% % [x x x x x]       = [potential_ ]
% % [x x x x x]         [potential_4]
% % [1 1 1 1 1]         [0]
% 
% % compute solution for parameters C
% % ---------------------------------
% meanvalues = mean(values); 
% values = values - meanvalues; % make mean zero
% C = pinv([Gelec;ones(1,length(Gelec))]) * [values(:);0];
% 
% % apply results
% % -------------
% Res = zeros(1,size(Gsph,1));
% for j = 1:size(Gsph,1)
%     Res(j) = sum(C .* Gsph(j,:)');
% end
% Res = Res + meanvalues;
% Res = reshape(Res, length(x(:)),1);

function [xbad, ybad, zbad, allres] = spheric_spline( xelec, yelec, zelec, xbad, ybad, zbad, values, params)

newchans = length(xbad);
numpoints = size(values,2);

%SPHERERES = 20;
%[x,y,z] = sphere(SPHERERES);
%x(1:(length(x)-1)/2,:) = []; xbad = [ x(:)'];
%y(1:(length(x)-1)/2,:) = []; ybad = [ y(:)'];
%z(1:(length(x)-1)/2,:) = []; zbad = [ z(:)'];

Gelec = computeg(xelec,yelec,zelec,xelec,yelec,zelec,params);
Gsph  = computeg(xbad,ybad,zbad,xelec,yelec,zelec,params);

% compute solution for parameters C
% ---------------------------------
meanvalues = mean(values); 
values = values - repmat(meanvalues, [size(values,1) 1]); % make mean zero

values = [values;zeros(1,numpoints)];
lambda = params(1);
C = pinv([Gelec+eye(size(Gelec))*lambda;ones(1,length(Gelec))]) * values;
%C = pinv([Gelec;ones(1,length(Gelec))]) * values;

clear values;
allres = zeros(newchans, numpoints);

% apply results
% -------------
for j = 1:size(Gsph,1)
    allres(j,:) = sum(C .* repmat(Gsph(j,:)', [1 size(C,2)]));        
end
allres = allres + repmat(meanvalues, [size(allres,1) 1]);

% compute G function
% ------------------
function g = computeg(x,y,z,xelec,yelec,zelec, params)

unitmat = ones(length(x(:)),length(xelec));
EI = unitmat - sqrt((repmat(x(:),1,length(xelec)) - repmat(xelec,length(x(:)),1)).^2 +... 
                (repmat(y(:),1,length(xelec)) - repmat(yelec,length(x(:)),1)).^2 +...
                (repmat(z(:),1,length(xelec)) - repmat(zelec,length(x(:)),1)).^2);

g = zeros(length(x(:)),length(xelec));
m = params(2);
maxn = params(3);

for n = 1:maxn % 200
    if ismatlab
        L = legendre(n,EI);
    else % Octave legendre function cannot process 2-D matrices
        for icol = 1:size(EI,2)
            tmpL = legendre(n,EI(:,icol));
            if icol == 1, L = zeros([ size(tmpL) size(EI,2)]); end
            L(:,:,icol) = tmpL;
        end
    end
    g = g + ((2*n+1)/(n^m*(n+1)^m))*squeeze(L(1,:,:));
end
g = g/(4*pi);