%% Copyright (c) 2011, INRA
%% 2003-2011, David Legland <david.legland@grignon.inra.fr>
%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
%%
%% All rights reserved.
%% (simplified BSD License)
%%
%% Redistribution and use in source and binary forms, with or without
%% modification, are permitted provided that the following conditions are met:
%%
%% 1. Redistributions of source code must retain the above copyright notice, this
%%    list of conditions and the following disclaimer.
%%     
%% 2. Redistributions in binary form must reproduce the above copyright notice, 
%%    this list of conditions and the following disclaimer in the documentation
%%    and/or other materials provided with the distribution.
%%
%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%% POSSIBILITY OF SUCH DAMAGE.
%%
%% The views and conclusions contained in the software and documentation are
%% those of the authors and should not be interpreted as representing official
%% policies, either expressed or implied, of copyright holder.

%% -*- texinfo -*-
%% @deftypefn {Function File} {@var{b} = } isPointOnEdge (@var{point}, @var{edge})
%% @deftypefnx {Function File} {@var{b} = } isPointOnEdge (@var{point}, @var{edge}, @var{tol})
%% @deftypefnx {Function File} {@var{b} = } isPointOnEdge (@var{point}, @var{edgearray})
%% @deftypefnx {Function File} {@var{b} = } isPointOnEdge (@var{pointarray}, @var{edgearray})
%% Test if a point belongs to an edge.
%
%   with @var{point} being [xp yp], and @var{edge} being [x1 y1 x2 y2], returns TRUE if
%   the point is located on the edge, and FALSE otherwise.
%
%   Specify an optilonal tolerance value @var{tol}. The tolerance is given as a
%   fraction of the norm of the edge direction vector. Default is 1e-14. 
%
%   When one of the inputs has several rows, return the result of the test
%   for each element of the array tested against the single parameter.
%
%   When both @var{pointarray} and @var{edgearray} have the same number of rows,
%   returns a column vector with the same number of rows.
%   When the number of rows are different and both greater than 1, returns
%   a Np-by-Ne matrix of booleans, containing the result for each couple of
%   point and edge.
%
%   @seealso{edges2d, points2d, isPointOnLine}
%% @end deftypefn

function b = isPointOnEdge(point, edge, varargin)

  % extract computation tolerance
  tol = 1e-14;
  if ~isempty(varargin)
      tol = varargin{1};
  end

  % number of edges and of points
  Np = size(point, 1);
  Ne = size(edge, 1);

  % adapt size of inputs if needed, and extract elements for computation
  if Np == Ne
      % When the number of points and edges is the same, the one-to-one test
      % will be computed, so there is no need to repeat matrices
      dx = edge(:,3) - edge(:,1);
      dy = edge(:,4) - edge(:,2);
      lx = point(:,1) - edge(:,1);
      ly = point(:,2) - edge(:,2);
      
  elseif Np == 1
      % one point, several edges
      dx = edge(:, 3) - edge(:, 1);
      dy = edge(:, 4) - edge(:, 2);
      lx = point(ones(Ne, 1), 1) - edge(:, 1);
      ly = point(ones(Ne, 1), 2) - edge(:, 2);
      
  elseif Ne == 1
      % several points, one edge
      dx = (edge(3) - edge(1)) * ones(Np, 1);
      dy = (edge(4) - edge(2)) * ones(Np, 1);
      lx = point(:, 1) - edge(1);
      ly = point(:, 2) - edge(2);

  else
      % Np points and Ne edges:
      % Create an array for each parameter, so that the result will be a
      % Np-by-Ne matrix of booleans (requires more memory, and uses repmat)

      x0 = repmat(edge(:, 1)', Np, 1);
      y0 = repmat(edge(:, 2)', Np, 1);
      dx = repmat(edge(:,3)', Np,  1) - x0;
      dy = repmat(edge(:,4)', Np,  1) - y0;
      
      lx = repmat(point(:, 1), 1, Ne) - x0;
      ly = repmat(point(:, 2), 1, Ne) - y0;
  end

  % test if point is located on supporting line
  b1 = (abs(lx.*dy - ly.*dx) ./ hypot(dx, dy)) < tol;

  % compute position of point with respect to edge bounds
  % use different tests depending on line angle
  ind     = abs(dx) > abs(dy);
  t       = zeros(size(dx));
  t(ind)  = lx( ind) ./ dx( ind);
  t(~ind) = ly(~ind) ./ dy(~ind);

  % check if point is located between edge bounds
  b = t >- tol & t-1 < tol & b1;

endfunction

%!shared points, vertices, edges

%!demo
%!   % create a point array
%!   points = [10 10;15 10; 30 10];
%!   % create an edge array
%!   vertices = [10 10;20 10;20 20;10 20];
%!   edges = [vertices vertices([2:end 1], :)];
%!
%!   % Test one point and one edge
%!   isPointOnEdge(points(1,:), edges(1,:))
%!   isPointOnEdge(points(3,:), edges(1,:))

%!demo
%!   % create a point array
%!   points = [10 10;15 10; 30 10];
%!   % create an edge array
%!   vertices = [10 10;20 10;20 20;10 20];
%!   edges = [vertices vertices([2:end 1], :)];
%!
%!   % Test one point and several edges
%!   isPointOnEdge(points(1,:), edges)'

%!demo
%!   % create a point array
%!   points = [10 10;15 10; 30 10];
%!   % create an edge array
%!   vertices = [10 10;20 10;20 20;10 20];
%!   edges = [vertices vertices([2:end 1], :)];
%!
%!   % Test several points and one edge
%!   isPointOnEdge(points, edges(1,:))'

%!demo
%!   % create a point array
%!   points = [10 10;15 10; 30 10];
%!   % create an edge array
%!   vertices = [10 10;20 10;20 20;10 20];
%!   edges = [vertices vertices([2:end 1], :)];
%!
%!   % Test N points and N edges
%!   isPointOnEdge(points, edges(1:3,:))'

%!demo
%!   % create a point array
%!   points = [10 10;15 10; 30 10];
%!   % create an edge array
%!   vertices = [10 10;20 10;20 20;10 20];
%!   edges = [vertices vertices([2:end 1], :)];
%!
%!   % Test NP points and NE edges
%!   isPointOnEdge(points, edges)

%!test
%!  p1 = [10 20];
%!  p2 = [80 20];
%!  edge = [p1 p2];
%!  p0 = [10 20];
%!  assert (isPointOnEdge(p0, edge));
%!  p0 = [80 20];
%!  assert (isPointOnEdge(p0, edge));
%!  p0 = [50 20];
%!  assert (isPointOnEdge(p0, edge));
%!  p0 = [9.99 20];
%!  assert (!isPointOnEdge(p0, edge));
%!  p0 = [80.01 20];
%!  assert (!isPointOnEdge(p0, edge));
%!  p0 = [50 21];
%!  assert (!isPointOnEdge(p0, edge));
%!  p0 = [79 19];
%!  assert (!isPointOnEdge(p0, edge));

%!test
%!  p1 = [20 10];
%!  p2 = [20 80];
%!  edge = [p1 p2];
%!  p0 = [20 10];
%!  assert (isPointOnEdge(p0, edge));
%!  p0 = [20 80];
%!  assert (isPointOnEdge(p0, edge));
%!  p0 = [20 50];
%!  assert (isPointOnEdge(p0, edge));
%!  p0 = [20 9.99];
%!  assert (!isPointOnEdge(p0, edge));
%!  p0 = [20 80.01];
%!  assert (!isPointOnEdge(p0, edge));
%!  p0 = [21 50];
%!  assert (!isPointOnEdge(p0, edge));
%!  p0 = [19 79];
%!  assert (!isPointOnEdge(p0, edge));

%!test
%!  p1 = [10 20];
%!  p2 = [60 70];
%!  edge = [p1 p2];
%!  assert (isPointOnEdge(p1, edge));
%!  assert (isPointOnEdge(p2, edge));
%!  p0 = [11 21];
%!  assert (isPointOnEdge(p0, edge));
%!  p0 = [59 69];
%!  assert (isPointOnEdge(p0, edge));
%!  p0 = [9.99 19.99];
%!  assert (!isPointOnEdge(p0, edge));
%!  p0 = [60.01 70.01];
%!  assert (!isPointOnEdge(p0, edge));
%!  p0 = [30 50.01];
%!  assert (!isPointOnEdge(p0, edge));

%!test
%!  edge = [10 20 80 20; 20 10 20 80; 20 10 60 70];
%!  p0 = [20 20];
%!  assert ([true ; true ; false], isPointOnEdge(p0, edge));

%!test
%!  k = 1e15;
%!  p1 = [10 20]*k;
%!  p2 = [60 70]*k;
%!  edge = [p1 p2];
%!  assert (isPointOnEdge(p1, edge));
%!  assert (isPointOnEdge(p2, edge));
%!  p0 = [11 21]*k;
%!  assert (isPointOnEdge(p0, edge));
%!  p0 = [59 69]*k;
%!  assert (isPointOnEdge(p0, edge));
%!  p0 = [9.99 19.99]*k;
%!  assert (!isPointOnEdge(p0, edge));
%!  p0 = [60.01 70.01]*k;
%!  assert (!isPointOnEdge(p0, edge));
%!  p0 = [30 50.01]*k;
%!  assert (!isPointOnEdge(p0, edge));

%!test
%!  k = 1e-10;
%!  p1 = [10 20]*k;
%!  p2 = [60 70]*k;
%!  edge = [p1 p2];
%!  assert (isPointOnEdge(p1, edge));
%!  assert (isPointOnEdge(p2, edge));
%!  p0 = [11 21]*k;
%!  assert (isPointOnEdge(p0, edge));
%!  p0 = [59 69]*k;
%!  assert (isPointOnEdge(p0, edge));
%!  p0 = [9.99 19.99]*k;
%!  assert (!isPointOnEdge(p0, edge));
%!  p0 = [60.01 70.01]*k;
%!  assert (!isPointOnEdge(p0, edge));
%!  p0 = [30 50.01]*k;
%!  assert (!isPointOnEdge(p0, edge));

%!test
%!  p1 = [10 20];
%!  p2 = [80 20];
%!  edge = [p1 p2];
%!  p0 = [10 20; 80 20; 50 20;50 21];
%!  exp = [true;true;true;false];
%!  assert (exp, isPointOnEdge(p0, edge));

%!test
%!  p1 = [10 20];
%!  p2 = [80 20];
%!  edge = [p1 p2];
%!  p0 = [40 20];
%!  exp = [true;true;true;true];
%!  assert (exp, isPointOnEdge(p0, [edge;edge;edge;edge]));

%!test
%!  edge1 = [10 20 80 20];
%!  edge2 = [30 10 30 80];
%!  edges = [edge1; edge2];
%!  p0 = [40 20;30 90];
%!  exp = [true;false];
%!  assert (exp, isPointOnEdge(p0, edges));
