function planehandle = plothalfplane(A,b,xbounds,ybounds,fighandle)
% plothalfplane - 2-d plotting tool that plots a half plane in 2-d, of the form A*x >= b
% usage: plothalfplane(A,b,xbounds,ybounds)
% usage: plothalfplane(A,b,xbounds,ybounds,fighandle)
% usage: planehandle = plothalfplane(A,b,xbounds,ybounds)
% usage: planehandle = plothalfplane(A,b,xbounds,ybounds,fighandle)
%
% the indicated half plane will always be a single color, 'blue', at a
% transparency of 0.25. These parameters and any others can be adjusted by
% accessing the graphics handle in planehandle.
%
% arguments: (input)
%   A - 1x2 vector, contains the coefficients of the bounding
%           relationship, of the form A(1)*x + A(2)*y >= b
%   b - constant
%   xbounds - vector of length 2 defining the lower and upper bounds for x
%   ybounds - vector of length 2 defining the lower and upper bounds for y
%
%   fighandle (optional) - scalar figure handle to plot on. if no figure
%     handle is supplied, the plothalfplane adds the indicated 2-d half
%     plane to the current figure, or if no figure currently exists, it
%     creates a new figure.
%
% arguments: (output)
%   planehandle - graphics handle to the plane as generated (not returned
%           if no output arguments are not requested from plothalfplane)

% was a figure handle supplied?
if (nargin < 5) || isempty(fighandle)
  fighandle = gcf;
end

% check the other arguments
if (numel(A) ~= 2) || ~isnumeric(A)
  error('A must be a numeric vector of length 2')
end

if ~isscalar(b) || ~isnumeric(b)
  error('b must be a numeric scalar')
end

if (nargin < 2) || isempty(xbounds)
  xbounds = xlim;
elseif  (numel(xbounds) ~= 2) || ~isnumeric(xbounds)
  error('If supplied, xbounds must be a numeric vector of length 2')
end

if (nargin < 3) || isempty(ybounds)
  ybounds = ylim;
elseif (numel(ybounds) ~= 2) || ~isnumeric(ybounds)
  error('If supplied, ybounds must be a numeric vector of length 2')
end

% just insure they are sorted
xbounds = sort(xbounds(:));
ybounds = sort(ybounds(:));

% there are 4 vertices for the box
pxy = [xbounds(1),ybounds(1);xbounds(2),ybounds(1);...
  xbounds(2),ybounds(2);xbounds(1),ybounds(2)];

% which vertices of the bounding box still remain. explicitly use > for this test
indremain = pxy*A(:) > b;

% treat as binary, converting to decimal
remaindec = [1 2 4 8]*indremain;

% the remaining vertices tell us which edges of the bounding box are crossed
% there are multiple possibilities we must consider
switch remaindec
  case 0
    % the plane does not see any part of the bounding box
    pxyplot = [];
  case 15
    % the plane sees the entire bounding box, so it is entirely covered
    pxyplot = pxy;
  case {1,14}
    % vertex 1, or vertices [2,3,4] are seen, so edges [1 2], [1 4] are crossed
    pxyplot = [pxy(indremain,:);xintercept(A,b,ybounds(1));yintercept(A,b,xbounds(1))];

  case {2,13}
    % vertex 2, or vertices [1,3,4] are seen, so edges [1 2], [2 3] are crossed
    pxyplot = [pxy(indremain,:);xintercept(A,b,ybounds(1));yintercept(A,b,xbounds(2))];

  case {4,11}
    % vertex 3, or vertices [1,2,4] are seen, so edges [2 3], [3 4] are crossed
    pxyplot = [pxy(indremain,:);xintercept(A,b,ybounds(2));yintercept(A,b,xbounds(2))];
    
  case {8,7}
    % vertex 4, or vertices [1,2,3] are seen, so edges [1 4], [3 4] are crossed
    pxyplot = [pxy(indremain,:);xintercept(A,b,ybounds(1));yintercept(A,b,xbounds(2))];
    
  case {3,12}
    % vertices [1 2], or vertices [3 4] are seen, so edges [2 3], [1 4] are crossed
    pxyplot = [pxy(indremain,:);yintercept(A,b,xbounds(1));yintercept(A,b,xbounds(2))];
    
  case {6,9}
    % vertices [2 3], or vertices [1 4] are seen, so edges [1 2], [3 4] are crossed
    pxyplot = [pxy(indremain,:);xintercept(A,b,ybounds(1));xintercept(A,b,ybounds(2))];

  otherwise
    % the cases for 5 or 10 are impossible to see

end

% At this point, the rows of pxyplot represent a patch. a convex hull will
% essentially sort them for us in a good order.
E = convhull(pxyplot(:,1),pxyplot(:,2));

planehandle = patch(pxyplot(E,1),pxyplot(E,2),'b');
planehandle.FaceAlpha = 0.1;

xlim(xbounds)
ylim(ybounds)

if nargout < 1
  clear planehandle;
end



function xyp = xintercept(A,b,yp)
% find the x-intercept of the line A(1)*x + A(2)*y = b, where it crosses
% a horizontal line at y = yp

% no need to worry about the case where A(1) == 0, since then an intersection
% with this edge would not have been found.
xp = (b - A(2)*yp)/A(1);
xyp = [xp,yp];

function xyp = yintercept(A,b,xp)
% find the y-intercept of the line A(1)*x + A(2)*y = b, where it crosses
% the vertical line at x = xp

% no need to worry about the case where A(2) == 0, since then an intersection
% with this edge would not have been found.
yp = (b - A(1)*xp)/A(2);
xyp = [xp,yp];



