%% Copyright (c) 2011, INRA
%% 2007-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{edge} =} clipLine (@var{line}, @var{box})
%% Clip a line with a box.
%% 
%%   @var{line} is a straight line given as a 4 element row vector: [x0 y0 dx dy],
%%   with (x0 y0) being a point of the line and (dx dy) a direction vector,
%%   @var{box} is the clipping box, given by its extreme coordinates: 
%%   [xmin xmax ymin ymax].
%%   The result is given as an edge, defined by the coordinates of its 2
%%   extreme points: [x1 y1 x2 y2].
%%   If line does not intersect the box, [NaN NaN NaN NaN] is returned.
%%   
%%   Function works also if @var{line} is a Nx4 array, if @var{box} is a Nx4 array, or
%%   if both @var{line} and @var{box} are Nx4 arrays. In these cases, @var{edge} is a Nx4
%%   array.
%%
%%   Example:
%%
%% @example
%%   line = [30 40 10 0];
%%   box = [0 100 0 100];
%%   res = clipLine(line, box)
%%   res = 
%%       0 40 100 40
%% @end example
%%
%% @seealso{lines2d, boxes2d, edges2d, clipEdge, clipRay}
%% @end deftypefn

function edge = clipLine(lin, bb, varargin)

  % adjust size of two input arguments
  if size(lin, 1)==1
      lin = repmat(lin, size(bb, 1), 1);
  elseif size(bb, 1)==1
      bb = repmat(bb, size(lin, 1), 1);
  elseif size(lin, 1) ~= size(bb, 1)
      error('bad sizes for input');
  end

  % allocate memory
  nbLines = size(lin, 1);
  edge    = zeros(nbLines, 4);

  % main loop on lines
  for i=1:nbLines
      % extract limits of the box
      xmin = bb(i, 1);
      xmax = bb(i, 2);
      ymin = bb(i, 3);
      ymax = bb(i, 4);
      
      % use direction vector for box edges similar to direction vector of the
      % line in order to reduce computation errors
      delta = hypot(lin(i,3), lin(i,4));
      
      
	  % compute intersection with each edge of the box
      
      % lower edge
	  px1 = intersectLines(lin(i,:), [xmin ymin delta 0]);
      % right edge
	  px2 = intersectLines(lin(i,:), [xmax ymin 0 delta]);
      % upper edge
	  py1 = intersectLines(lin(i,:), [xmax ymax -delta 0]);
      % left edge
	  py2 = intersectLines(lin(i,:), [xmin ymax 0 -delta]);
      
      % remove undefined intersections (case of lines parallel to box edges)
      points = [px1 ; px2 ; py1 ; py2];
      points = points(isfinite(points(:,1)), :);
	
      % sort points according to their position on the line
      pos = linePosition(points, lin(i,:));
      [pos inds] = sort(pos); %#ok<ASGLU>
      points = points(inds, :);
      
      % create clipped edge by using the two points in the middle
      ind = size(points, 1)/2;
      inter1 = points(ind,:);
      inter2 = points(ind+1,:);
      edge(i, 1:4) = [inter1 inter2];
      
      % check that middle point of the edge is contained in the box
      midX = mean(edge(i, [1 3]));
      xOk = xmin <= midX && midX <= xmax;
      midY = mean(edge(i, [2 4]));
      yOk = ymin <= midY && midY <= ymax;
      
      % if one of the bounding condition is not met, set edge to NaN
      if ~(xOk && yOk)
          edge (i,:) = NaN;
      end
  end
endfunction

%!demo
%!   lin = [30 40 10 0];
%!   bb = [0 100 0 100];
%!   res = clipLine(lin, bb)
%!
%! drawBox(bb,'color','k');
%! line(lin([1 3]),lin([2 4]),'color','b');
%! line(res([1 3]),res([2 4]),'color','r','linewidth',2);
%! axis tight
%! v = axis ();
%! axis(v+[0 10 -10 0])

%!test % inside, to the right % inside, to the left% outside
%! bb = [0 100 0 100];
%! lin = [30 40 10 0];
%! edge = [0 40 100 40];
%! assert (edge, clipLine(lin, bb), 1e-6);
%! lin = [30 40 -10 0];
%! edge = [100 40 0 40];
%! assert (edge, clipLine(lin, bb), 1e-6);
%! lin = [30 140 10 0];
%! assert (sum(isnan(clipLine(lin, bb)))==4);

%!test % inside, upward % inside, downward % outside
%!  bb = [0 100 0 100];
%!  lin = [30 40 0 10];
%!  edge = [30 0 30 100];
%!  assert (edge, clipLine(lin, bb), 1e-6);
%!  lin = [30 40 0 -10];
%!  edge = [30 100 30 0];
%!  assert (edge, clipLine(lin, bb), 1e-6);
%!  lin = [140 30 0 10];
%!  assert (sum(isnan(clipLine(lin, bb)))==4);

%!test % inside, top right corner% inside, down right corner % outside
%!  bb = [0 100 0 100];
%!  lin = [80 30 10 10];
%!  edge = [50 0 100 50];
%!  assert (edge, clipLine(lin, bb), 1e-6);
%!  lin = [20 70 10 10];
%!  edge = [0 50 50 100];
%!  assert (edge, clipLine(lin, bb), 1e-6);
%!  lin = [140 -30 10 10];
%!  assert (sum(isnan(clipLine(lin, bb)))==4);
%!  lin = [-40 130 10 10];
%!  assert (sum(isnan(clipLine(lin, bb)))==4);

%!test %multilines % inside, top right corner
%!  bb = [0 100 0 100];
%!  lin = [...
%!      80 30 10 10; ...
%!      20 70 10 10; ...
%!      140 -30 10 10; ...
%!      -40 130 10 10];
%!  edge = [...
%!      50 0 100 50; ...
%!      0 50 50 100; ...
%!      NaN NaN NaN NaN; ...
%!      NaN NaN NaN NaN; ...
%!      ];
%!  clipped = clipLine(lin, bb);
%!  assert (4, size(clipped, 1));
%!  assert (edge(1:2, :), clipped(1:2, :), 1e-6);
%!  assert (sum(isnan(clipped(3,:)))==4);
%!  assert (sum(isnan(clipped(4,:)))==4);

%!test % test clipping of horizontal lines % inside, to the right
%!  bb = [-1 1 -1 1]*1e10;
%!  lin = [3 0 1 2];
%!  D = 1e10;
%!  edge = [3-D/2 -D 3+D/2 D];
%!  clipped = clipLine(lin, bb);
%!  assert (edge, clipped);

%!test % inside, to the right
%!  bb = [-1 1 -1 1]*100;
%!  lin = [3 0 1*1e10 2*1e10];
%!  D = 100;
%!  edge = [3-D/2 -D 3+D/2 D];
%!  clipped = clipLine(lin, bb);
%!  assert (edge, clipped, 1e-6);
