function auc = comp_auc(px,py,a,b)
% This function computes areas under the curve (AUCs) for trajectory data
% with the 2D-coordinates px and py and a line, defined by slope a and
% intercept b. The arguments px and py may be scalars or vectors of the
% same length. In the latter case, ad will be a vector of distances.
%
% AUCs are computed by dividing the area in triangular and square-shaped
% pieces and adding the up across the trajectory. The output is a vector
% representing the cumulative AUC across course of the trajectory.
%
%                                                           Roland Pfister
%                                                           16.11.2011
%                                                           03.08.2016(v2)
% Get distances from line and matching line coordinates.
[ad, direction, linecoords] = comp_ad(px,py,a,b);

cum_auc = ones(1,numel(ad));
cum_auc(1,1) = 0;
for i = 2:numel(ad)
    % Only increase cum_auc if coordinates change.
    if linecoords.x(i)==linecoords.x(i-1) && ...
       linecoords.y(i)==linecoords.y(i-1)
       cum_auc(1,i) = cum_auc(1,i-1);
    else
        % Compute distance on line.
        ldist = ((linecoords.x(i)-linecoords.x(i-1))^2 + ...
                (linecoords.y(i)-linecoords.y(i-1))^2)^0.5;
        % Check for line flips
        if direction(i) == direction(i-1)
            % If no flip is detected: Increment cum_auc (1) by 
            % rectangle (2) and triangle (3).
            cum_auc(1,i) = cum_auc(1,i-1) + direction(i)  * ( ... % (1)
                ldist * min(ad(i),ad(i-1))                  + ... % (2)
                0.5 * sign(ad(i)-ad(i-1)) * (ad(i)-ad(i-1)) * ... % (3)
                ldist );           
        else
            % Add triangle or trapezoid (verschraenktes Trapez, Wiki)
            % if trajectory hits or crosses ideal line.
            if direction(i)==0
                cum_auc(1,i) = cum_auc(1,i-1) + direction(i-1) * ...
                    (ldist * (ad(i-1)-ad(i)) / 2 );
            elseif direction(i-1)==0
                cum_auc(1,i) = cum_auc(1,i-1) + direction(i) * ... 
                    (ldist * (ad(i)-ad(i-1)) / 2 );            
            else
                cum_auc(1,i) = cum_auc(1,i-1) + direction(i) * ... 
                    (ldist * (ad(i)-ad(i-1)) / 2 );
            end;
            
        end;
    
    end;
    
   
end;

auc = cum_auc;

end

