function [channels,overall_amplitude] = eeg_amplitudearea(EEG, channels, resrate, wstart, wend)
% EEG_AMPLITUDEAREA - Resamples an ERP average using spline interpolation 
%                       at a new sample rate (resrate) in Hz to get the exact limits 
%                       of the window of integration. Finely samples the window 
%                       and adds together very narrow rectangles capped by 
%                       right-angled triangles under the window. Output is in uV. 
%                       Trade-off between speed and number of resamples and number of 
%                        channels selected occurs.
% Usage:
%      >> [channels, amplitude] = eeg_amplitudearea(EEG,channels, resrate, wstart, wend);
% Inputs:
%   EEG         - EEGLAB data struct containing a (3-D) epoched data matrix 
%   channels    - vector of channel indices
%   resrate     - resampling rate for window of integration in Hz
%   wstart      - start of window of integration in ms post stimulus-onset
%   wend        - end of window of integration in ms post stimulus-onset
%
% Outputs:
%   channels    - a vector of channel indices.
%   amplitude   - 1-dimensional array in uV for the channels
%
% Example
%    >> [channels, amplitude] = eeg_amplitudearea(EEG,[12 18 25 29], 2000, 90.52, 120.52);
%
% Author: Tom Campbell, Helsinki Collegium for Advanced Studies, Biomag Laboratory, 
%         Engineering Centre, Helsinki University Central Hospital Helsinki Brain 
%         Research Centre (tom.campbell@helsinki.fi) Spartam nanctus es: Hanc exorna. 
%         Combined with AMPLITUDEAREA_MSUV by Darren Weber, UCSF 28/1/05
%         Retested and debugged Tom Campbell 2/2/05
%         Reconceived, factored somewhat, tested and debugged Tom Campbell 13:24 23.3.2005

if wstart > wend
    error ('ERROR: wstart must be greater than wend')
else
    [channels, overall_amplitude] = eeg_amplitudearea_msuV (EEG,channels, resrate, wstart, wend);
    overall_amplitude = overall_amplitude/(wend - wstart);
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [channels, overall_area] = eeg_amplitudearea_msuV (EEG, channels, resrate, wstart, wend)

%if ndim(EEG.data) ~= 3
%  error('EEG.data must be 3-D data epochs');
%end
erp = mean(EEG.data,3);

[tmp ind1] =min( abs( EEG.times - wstart ) ); % closest time index to wstart
[tmp ind2] =min( abs( EEG.times - wend ) ); % closest time index to wend
restep = 1/resrate;

if EEG.times(ind1) > wstart 
    ind1= ind1 -1;
end

if EEG.times(ind2) < wend 
    ind2= ind2 +1;
end    

for x= ind1:ind2
    t = (x -ind1)+1;
    tim(t) = EEG.times(x);
end

tr = 1;
timr(tr) = wstart;
while timr(tr) < wend
    tr = tr + 1;
    timr(tr) = timr(tr-1)+ restep;
end

for x = 1:size(channels,2)
    channel = channels(x);
    %resamples
    rerp(x, 1:tr) = spline(tim(:),erp(channel, ind1:ind2), timr(1:tr));
    pent = timr(tr - 1);
    overall_area(x) = 0;
    for y = 1:(tr -1)
        v1 =  rerp(x,(y));
        v2 =  rerp(x,(y+1));
        if ((v1 > 0) && (v2 < 0)) || ((v1 < 0) && (v2 > 0))
            if (y == (tr-1)) && (timr(y+1)> wend)
                area1 = zero_crossing_truncated(v1, v2, restep, wend, pent);
            else    
                area1 = zero_crossing(v1, v2, restep);
            end
        else
            if( y == (tr-1)) && (timr(y+1)> wend)
                area1 = rect_tri_truncated(v1, v2, restep,wend,pent);
            else
                area1 = rect_tri(v1, v2, restep);
            end
        end
        overall_area(x) = overall_area(x) + area1;
    end
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [area] = zero_crossing(v1,v2,step)
    if (v1 > v2)
        T1 = v1;
        T2 = v2;
    else
        T1 = v2;
        T2 = v1;
    end
    tantheta = (abs(T1)+ abs(T2))/step;
    if (v1 > v2)
        %decline
        z = abs(T1)/tantheta;
        tr1= abs(T1)*(z/2);
        tr2= abs(T2)*((step-z)/2);
    else
        %incline
        z = abs(T2)/tantheta;
        tr2= abs(T2)*(z/2);
        tr1= abs(T1)*((step-z)/2);
    end
    [area] = (tr1 - tr2);
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [area] = zero_crossing_truncated(v1,v2,step,wend,pent)
    if (v1 > v2)
        T1 = v1;
        T2 = v2;
    else
        T1 = v2;
        T2 = v1;
    end
    tantheta = (abs(T1)+ abs(T2))/step;
    s = wend - pent;
    if (v1 > v2)
        z = abs(T1)/tantheta
        if s < z   
            %decline,truncated before zerocrossing
            t1 = tantheta * s;
            r1 = abs(T1)-abs(t1);
            tr1= abs(t1)*(s/2);
            tr2= 0;
            rect1 = r1*s;
            rect2 = 0;
        else
            %decline,truncated after zerocrossing
            t2= tantheta*(s-z);
            tr1= abs(T1)*(z/2);
            tr2 = abs(t2)*((s-z)/2);
            rect1 = 0;
            rect2 = 0;
        end    
    else
        z = abs(T2)/tantheta;
        if s < z
            %incline,truncated before zerocrossing
            t2 = tantheta * s;
            r2 = abs(T2)-abs(t2);
            tr1= 0;
            tr2= abs(t2)*(s/2);
            rect1 = 0;
            rect2 = r2*s;
        else
            %incline,truncated after zerocrossing
            t1= tantheta*(s-z);
            tr1 = abs(t1)*((s-z)/2);
            tr2 = abs(T2) * (z/2);
            rect1 = 0;
            rect2 = 0;
        end
    end

[area] = ((rect1 + tr1) - (rect2 + tr2));
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [area] = rect_tri(v1,v2,step)
    if (abs(v1) > abs(v2))
        T = abs(v1)-abs(v2);
        R = abs(v2);
    else
        T = abs(v2)-abs(v1);
        R = abs(v1);
    end
    rect = R*step;
    tri = T*(step/2);
    if v1 > 0
        area = 1* (rect+tri);
    else
        area = -1 * (rect+tri);
    end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [area] = rect_tri_truncated(v1,v2,step,wend,pent)
    if (abs(v1) > abs(v2))
        T = abs(v1)-abs(v2);
        R = abs(v2);
    else
        T = abs(v2)-abs(v1);
        R = abs(v1);
    end
    
    tantheta = abs(T)/step;
    s = wend -pent;
    
    if (v1>0)
        if v1 >v2
        %positive decline
            t = tantheta*s;
            e = abs(T)-abs(t);
            rect = s*R;
            exrect = s*e;
            tri = (s/2)*R;
        else
        %positive incline
            t = tantheta*s;
            rect = s*R;
            exrect = 0;
            tri = (s/2)*R;
        end
    else
       if v1 >v2
        %negative decline
            t = tantheta*s;
            rect = s*R;
            exrect = 0;
            tri = (s/2)*R;
        else
        %negative incline 
            t = tantheta*s;
            e = abs(T)-abs(t);
            rect = s*R;
            exrect = s*e;
            tri = (s/2)*R;
        end
    end
    tri = T*(step/2);
    if v1 > 0
        area = 1* (rect+exrect+tri);
    else
        area = -1 * (rect+exrect+tri);
    end
return
