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));