function [ad, direction, linecoords] = comp_ad(px,py,a,b)
% This function computes the absolute (Euclidean) distance of a point P
% (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.
%
% Additional output arguments are direction (1 = above the line, 0 = on
% the line, -1 = below the line) and linecoords (structure with the fields
% x and y for x0 and y0).
%
%                                                           Roland Pfister
%                                                           16.11.2011
%
% Background:
%
% The Euclidean distance d between P and the line defined by a and b can
% be computed as a function of the x coordinate of a point G(x,y) on the
% line (general case). To avoid confusion, I will denote px and py as p1
% and p2, respectively.
%
%    d(x) = sqrt( (p1-x)^2 + (p2-y)^2 )       | ()^2
%
%    d(x)^2 = (p1-x)^2 + (p2-y)^2             | solve brackets
%
%    d(x)^2 = p1^2 - 2*p1*x + x^2 + p2^2 - 2*p2*a*x - ...
%             2*p2*b + 2*a*x*b + a^2*x^2 + b^2
%
%    Rearranging to bring all terms including x's to the front...
%
%    d(x)^2 = (a^2+1)*x^2 + (2*a*b - 2*p1 - 2*p2*a)*x + ...
%             p1^2 + p2^2 - 2*p2*b + b^2
%
%    Deriving this latter function helps to find the minimum
%    (absolute) distance between P and the line (note that none
%    of the terms in the 2nd line of the above equation contains
%    any x's):
%
%    d(x)^2' = 2*(a^2+1)*x + 2*(a*b - p1 - p2*a)
%
%    d(x)^2' =(!)= 0
%
% The x coordinate corresponding to the global minimum is denoted as x0:
%
%    x0 = (p1 + p2*a - a*b) / (a^2+1)
%
x0 = (px + py*a - a*b) ./ (a^2 + 1);
y0 = a*x0 + b;
linecoords.x = x0;
linecoords.y = y0;

% Finally, ad is computed as the simple Euclidean distance between the two
% points P(px,py) and G(x0,y0):
%
ad = sqrt((px-x0).^2 + (py-y0).^2);

% Compute direction.
direction = sign(py - (a*px+b));



