% EEG_LAPLAC - gives the laplacian for the data contained in EEG
% for the values 1:stepl:end with electrodes position given by EEG.chanlocs. 
% Step must be a number between 1 and length(EEG.data) 
% The method used for computig the laplacian is the described by Perrin et al.
% 1989 using a spline interpolation between the electrodes. This function 
% requires the functions G2 and GM_1 to be in the path or in the same dir. 
% Usage:
%        >>  [laplac time] = eeg_laplac(EEG, step);
%        >>  laplac = eeg_laplac(EEG, step);
% Inputs:
%   EEG     - A EEG structure containing the EEG.data an the EEG.chanlocs
%   step    - The step for subsampling the data
%    
% Outputs:
%   laplac     - A matrix containing the value of laplacian for the data
%   time       - Time spent in calculation  
%
% See also: DEL2MAP
%
% Author: Juan Sebastian Gonzalez, DFI / Universidad de Chile
% 
% Message from Jrgen Kayser, one of the main author of the CSD (Current
% Source Density) toolbox: "I do not find any differences between his Laplacian 
% solution and our CSD transform (i.e., apart from the lack of a smoothing
% constant, the choice of the m constant, and the specific spherical locations 
% assigned to any given electrode). Even more important would be a warning that
% the flexibility constant m is fixed at 3 (2 = max flexible, 3 = medium, 
% 4=medium rigid, 5+=increasingly less flexible splines), which is a debatable 
% choice and not recommended by us. The choice of the m constant (and the 
% regulization or smoothing constant) should depend on the objective of the
% study."

% Copyright (C)2006  Juan Sebastian Gonzalez, DFI / Universidad de Chile
%                                          
% 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 [laplac, time] = eeg_laplac(EEG, paso)

if nargin < 2
   help eeg_laplac;
   return;
end

tic
n = length(EEG.data(:,1));
a = zeros(n,3);
g = zeros(n);
for i=1:n
    a(i,1)=EEG.chanlocs(i).sph_phi;
    a(i,2)=EEG.chanlocs(i).sph_theta;
end
V = zeros(1,n+1);
aux = a*pi/180;
matriz_g = zeros(n+1);
for i=1:n
    matriz_g(i,1:n) = g2(aux(i,1)*ones(n,1),aux(i,2)*ones(n,1),aux(:,1),aux(:,2));
    g(i,:) = gm_1(aux(i,1)*ones(n,1),aux(i,2)*ones(n,1),aux(:,1),aux(:,2));
end
matriz_g(:,n+1) = 1;
matriz_g(n+1,:) = 1;
inv_m = inv(matriz_g);
k = 1;
largo = length(EEG.data(1,:));
laplac = zeros(n,round((largo-1)/paso));
for j=1:paso:largo
    % ciclo para sacar los anglos phi y theta en grados    
    for i=1:n
        V(i) = EEG.data(i,j);
    end
    % calculo la constante a partir de la matriz invertida 
    constant = inv_m*V';
    % una vez que tengo las constantes, calculo el laplaciano segun la formula
    % de perrin 89 corregida
    for i=1:n
        laplac(i,k) = g(i,:)*constant(1:n);
    end
    k = k + 1;
end
time = toc;
%---------------------------------------------------------------------
% internal functions g2 gm_1
function valor=g2(r1,r2,r3,r4)
%g entrega la sumatoria hasta 7 con m=3 de la ec. (2) de la technical note del
%Electrical Geodesics, salvo que esta esta bien

x=cos(r1).*cos(r2).*cos(r3).*cos(r4) + ...
    cos(r1).*sin(r2).*cos(r3).*sin(r4) + ...
    sin(r1).*sin(r3);
L = length(x);
for i=1:L
    if x(i) > 1
        x(i) = 1;
    
    elseif x(i) < -1
        x(i) = -1;
        
    end
end

P = zeros(25,length(x));
for n=1:25
    aux = legendre(n,x); % calculo los polinomios de legendre de orden n %
    P(n,:) = aux(1,:); % tomo el que me sirve %
end

n = 1:25;

for i=1:25
    aux2(i,:) = ((2*n(i)+1)./(n(i).*(n(i)+1)).^3)*ones(1,L);
end
    valor = 1/4/pi*sum(aux2.*P);
 %-----------------------------------------------------   
function valor=gm_1(r1,r2,r3,r4)
%g entrega la sumatoria hasta 7 con m=3 de la ec. (2) de la technical note del
%Electrical Geodesics

x=cos(r1).*cos(r2).*cos(r3).*cos(r4) + ...
    cos(r1).*sin(r2).*cos(r3).*sin(r4) + ...
    sin(r1).*sin(r3);
L = length(x);
for i=1:L
    if x(i) > 1
        x(i) = 1;
    
    elseif x(i) < -1
        x(i) = -1;
        
    end
end
P = zeros(25,length(x));
for n=1:25
    aux=legendre(n,x);
    P(n,:)=aux(1,:);
end

n=1:25;
%el problema es calcular esto para cada columna...sin hacer un for
%respuesta, hacer for de 7 vueltas jajaj
for i=1:25
    aux2(i,:)=(2*n(i)+1)/(n(i)*(n(i)+1))^2*ones(1,L);
end

    valor=1/4/pi*sum(aux2.*P);
    


