function vispoly2mask(x,y,M,N)
%VISPOLY2MASK Visualize the algorithm used by poly2mask.
%   VISPOLY2MASK(X,Y,M,N) displays the output image of poly2mask, and then
%   superimposes:
%      - The 5-by-5 subpixel resolution grid along the polygon edges
%      - The original polygon
%      - The horizontal edges of the modified polygon
%
%   Example
%   -------
%       x = [10 35 40 10];
%       y = [10 40 15 10];
%       vispoly2mask(x,y,50,50)

%   $Revision$
%   Copyright 2006 The MathWorks, Inc.

bw = poly2mask(x,y,M,N);

I = ones(size(bw));
I(bw) = 0.8;
imshow(I,'InitialMag','fit')

hold on
plot(x,y,'b','LineStyle','--','Marker','*')

xq = (round(5*(x - 0.5) + 1) - 1)/5 + 0.5;
yq = (round(5*(y - 0.5) + 1) - 1)/5 + 0.5;
plot(xq,yq,'r','LineStyle','--','Marker','*')

[xe,ye] = poly2edgelist(x,y,5);

for k = 1:length(xe)
    xk = floor((xe(k) - 1)/5) + 1;
    yk = floor((ye(k) - 1)/5) + 1;

    % Plot subpixel grid.
    [xxsub,yysub] = meshgrid(linspace(xk-.5, xk+.5, 6), ...
                             linspace(yk-.5, yk+.5, 6));
    plot(xxsub,yysub,'Color',[.6 .6 .6]);
    plot(xxsub',yysub','Color',[.6 .6 .6]);

    % Plot pixel border
    xx = [xk-.5 xk+.5 xk+.5 xk-.5 xk-.5];
    yy = [yk-.5 yk-.5 yk+.5 yk+.5 yk-.5];
    plot(xx,yy,'k','LineWidth',2);
    
    % Plot edge
    xedge = (xe(k) - 1)/5 + 0.5;
    yedge = (ye(k) - 1)/5 + 0.5;
    plot([xedge xedge+0.2], [yedge yedge], 'r', 'LineWidth', 4)
end

axis on
box on
hold off

function [xe, ye] = poly2edgelist(x,y,scale)
%POLY2EDGELIST Computes list of horizontal edges from polygon.
%   [XE,YE] = POLY2EDGELIST(X,Y) rescales the polygon represented
%   by vectors X and Y by a factor of 5, and then quantizes them to
%   integer locations.  It then creates a new polygon by filling in
%   new vertices on the integer grid in between the original scaled
%   vertices.  The new polygon has only horizontal and vertical
%   edges.  Finally, it computes the locations of all horizontal
%   edges, returning the result in the vectors XE and YE.
%
%   [XE,YE] = POLY2EDGELIST(X,Y,SCALE) rescales by a factor of
%   SCALE instead of 5.  SCALE must be a positive odd integer.
%
%   See also EDGELIST2MASK, POLY2MASK, ROIPOLY.

%   Copyright 1993-2003 The MathWorks, Inc.
%   $Revision: 1.1.4.2 $  $Date: 2003/01/11 00:06:17 $

if nargin < 3
    scale = 5;
end

% Scale and quantize (x,y) locations to the higher resolution grid.
x = round(scale*(x - 0.5) + 1);
y = round(scale*(y - 0.5) + 1);

num_segments = length(x) - 1;
x_segments = cell(num_segments,1);
y_segments = cell(num_segments,1);
for k = 1:num_segments
    [x_segments{k},y_segments{k}] = intline(x(k),x(k+1),y(k),y(k+1));
end

% Concatenate segment vertices.
x = cat(1,x_segments{:});
y = cat(1,y_segments{:});

% Horizontal edges are located where the x-value changes.
d = diff(x);
edge_indices = find(d);
xe = x(edge_indices);

% Wherever the diff is negative, the x-coordinate should be x-1 instead of
% x.
shift = find(d(edge_indices) < 0);
xe(shift) = xe(shift) - 1;

% In order for the result to be the same no matter which direction we are
% tracing the polynomial, the y-value for a diagonal transition has to be
% biased the same way no matter what.  We'll always chooser the smaller
% y-value associated with diagonal transitions.
ye = min(y(edge_indices), y(edge_indices+1));

function [x,y] = intline(x1, x2, y1, y2)
%INTLINE Integer-coordinate line drawing algorithm.
%   [X, Y] = INTLINE(X1, X2, Y1, Y2) computes an
%   approximation to the line segment joining (X1, Y1) and
%   (X2, Y2) with integer coordinates.  X1, X2, Y1, and Y2
%   should be integers.  INTLINE is reversible; that is,
%   INTLINE(X1, X2, Y1, Y2) produces the same results as
%   FLIPUD(INTLINE(X2, X1, Y2, Y1)).

%   Copyright 1993-2003 The MathWorks, Inc.  
%   $Revision: 5.11.6.2 $  $Date: 2003/08/28 18:27:02 $

dx = abs(x2 - x1);
dy = abs(y2 - y1);

% Check for degenerate case.
if ((dx == 0) && (dy == 0))
  x = x1;
  y = y1;
  return;
end

flip = 0;
if (dx >= dy)
  if (x1 > x2)
    % Always "draw" from left to right.
    t = x1; x1 = x2; x2 = t;
    t = y1; y1 = y2; y2 = t;
    flip = 1;
  end
  m = (y2 - y1)/(x2 - x1);
  x = (x1:x2).';
  y = round(y1 + m*(x - x1));
else
  if (y1 > y2)
    % Always "draw" from bottom to top.
    t = x1; x1 = x2; x2 = t;
    t = y1; y1 = y2; y2 = t;
    flip = 1;
  end
  m = (x2 - x1)/(y2 - y1);
  y = (y1:y2).';
  x = round(x1 + m*(y - y1));
end
  
if (flip)
  x = flipud(x);
  y = flipud(y);
end



