function [ do_intersect, are_parallel, x, alpha, beta ] = intersect_half_lines( a, u, b, v, epsilon )
  
%  function [ do_intersect, are_parallel, x, alpha, beta ] = intersect_half_lines( a, u, b, v, epsilon )
%
% Find the intersection (if any) of two half lines:
%    x = a + alpha * u,   alpha > 0
%    x = b + beta  * v,   beta  > 0
% 
%   Inputs: 
%
% The 2-by-1 vectors a and b are the starting point of each half-line.
% The 2-by-1 vectors u and v are the directions of each half-line.
%
% Optionally, epsilon is a very small value for equality tests
% ( default 1e-10).
%
%
%   Outputs:
%
% do_intersect and are_parallel are binary flags (0 or 1).
%
% The 2-by-1 vector x is defined only when do_intersect == 1 and
% are_parallel == 0 otherwise a [] value is returned.
%
% --- G. Lathoud 2005
  
  
  if nargin < 4
    error( [ mfilename ' needs at least 4 input arguments.' ] );
  end
  
  if ~exist( 'epsilon', 'var' )
    epsilon = 1e-10;
  end
  
  % Make sure we have 2 x 1 column vectors
  % ( throwing errors makes it easier to debug higher-level code. )
  if ~all( size( a ) == [ 2 1 ] )
    error( [ mfilename ': input argument "u" must be a 2 x 1 vector.' ] );
  end
  
  if ~all( size( u ) == [ 2 1 ] )
    error( [ mfilename ': input argument "u" must be a 2 x 1 vector.' ] );
  end
  
  if ~all( size( b ) == [ 2 1 ] )
    error( [ mfilename ': input argument "b" must be a 2 x 1 vector.' ] );
  end
  
  if ~all( size( v ) == [ 2 1 ] )
    error( [ mfilename ': input argument "v" must be a 2 x 1 vector.' ] );
  end
  
    
  delta = b - a;

  % Potential intersection point
  %      x = a + alpha * u
  % and  x = b + beta * v
  %
  % In case that there is no intersection, we return [].
  
  x = [];
  alpha = [];
  beta = [];
  
  % Are they parallel?
  
  DELTA = det( [ u, -v ] );
  
  are_parallel   = ( abs( DELTA ) < epsilon );
  
  
  do_intersect   = 0;
  
  if are_parallel
    
    % Check whether the supporting lines intersect
    do_intersect = abs( det( [ u, delta ] ) ) < epsilon;
    
    if do_intersect
      
      % Check whether the parallel half-lines intersect, which amounts
      % to exclude only one case: 
      %
      % "opposite directions u,v" and "the middle (a+b)/2 does not
      % belong to at least one of the half-lines."
      
      middle = ( a + b ) / 2;
      
      do_intersect = ~( ( dot( u, v ) < epsilon ) & ( ( dot( middle - a, u ) < epsilon ) | ( dot( middle - b, v ) < epsilon ) ) );
      
    end
    
  else
    
    % Not parallel: find the intersection point of the supporting lines
    alpha = det( [ delta, -v ] ) / DELTA;
    beta  = det( [ u, delta  ] ) / DELTA;
    
    % Check whether the intersection point belongs to both half-lines
    do_intersect = ( alpha > -epsilon ) & ( beta > -epsilon );
  
    if do_intersect
      % "mean" to reduce numerical imprecision errors
      x = mean( [ a + alpha * u,    b + beta * v ], 2 );
    end
    
  end