function pic = polygon(pic, px, py, edge, interior)
% pic = polygon(pic,px,py,edge,interior): draw convex ascii polygon
%     $pic$ is the ascii background
%     $px$, $py$ are vectors of the same length $n$, holding
%     the $n$ vertices of a non-degenerate, convex $n$-gon;
%        coordinates in $px$, $py$ are not necessarily integers
%     if $edge$ is a single character, it is used to draw the edge of the
%        polygon, which the points of distance <.5 from any edge of the polygon
%     if $interior$ is a single character, it is used to fill the interior of
%        the polygon, which does *not* include the edge
%     clipping is performed: no drawing is done outside of $pic$
%
% example:
%
%    clear pic; pic(20,50) = '.'; pic(:,:) = '.'
%    pic = polygon(pic, [2 60 40 10], [5 5 20 15], '+', 'o')
%    pic = polygon(pic, [13 13 47 47], [8 12 12 8], '#', '*')
%    pic = polygon(pic, [13 13 47 47], [8 12 12 8], '', '-')
%    pic = polygon(pic, [1 30 1], [1 20 20], '\', '')

% slower, unoptimized, unvectorized solution
    tic
    pic2 = pic; % don't modify $pic$, which is used in the other solution
    n = length(px); % number of vertices
    for y = 1:size(pic2,1)
        for x = 1:size(pic2,2)
            isin = 1; % == "does $(x,y)$ still look like it is inside so far?"
            isedge = 0; % == "is $(x,y)$ within .5 of any edge so far?"
            for i = 1:n
                j = rem(i,n) + 1; % index of vertex after $i$ with "wraparound"
                k = rem(j,n) + 1; % index of vertex after $j$ with "wraparound"
                A = py(j) - py(i);
                B = - (px(j) - px(i));
                C = - (A * px(i) + B * py(i));
                d = (A * x + B * y + C) / sqrt( A^2 + B^2 );
                isin = isin & ( sign(d) == sign(A * px(k) + B * py(k) + C) );
                inrect = min(px(i),px(j))-.5 < x & x < max(px(i),px(j))+.5 ...
                       & min(py(i),py(j))-.5 < y & y < max(py(i),py(j))+.5;
                isedge = isedge | (abs(d) < .5 & inrect);
            end
            if isedge 
                if length(edge) == 1
                    pic2(y,x) = edge; 
                end
            elseif isin
                if length(interior) == 1
                    pic2(y,x) = interior;
                end
            end
        end % for x
    end % for y
    toc

% faster, optimized, vectorized solution
    tic
    n = length(px); % number of vertices
    % compute positions $(yy,xx)$ within bounding box
        xlo = max(1,floor(min(px))); xhi = min(size(pic,2),ceil(max(px))); 
        ylo = max(1,floor(min(py))); yhi = min(size(pic,1),ceil(max(py)));
        xx = xlo:xhi;
        yy = ylo:yhi;
    
    % duplicate 1st vertex so can connect positions $i$ and $i+1$ for $i=1:n$
        px = [px px(1)];
        py = [py py(1)];
    % compute $A,B,C$ for the equation $Ax+By+C=0$ for each edge
    % such that $A^2 + B^2 = 1$
        A = diff(py);
        B = -diff(px);
        D = sqrt(A .^ 2 + B .^ 2);
        A = A ./ D;
        B = B ./ D;
        C = - ( A .* px(2:end) + B .* py(2:end) );
    % duplicate 2nd vertex so there is a vertex at position $i+2$ for $i=1:n$
        px = [px px(2)];
        py = [py py(2)];
    
    % classifications of points: inside, on boundary, or outside the polygon
        ISIN = 3; ISEDGE = 2; ISOUT = 1;
    % use vectorization tricks from E11 to classify each position:
    % $mini(i,j)$ holds classification for point $(xx(j), yy(i))$
        mini(length(yy), length(xx)) = ISIN; mini(:,:) = ISIN;
        x = repmat(xx, length(yy), 1); % x-position of each position
        y = repmat(yy', 1, length(xx)); % y-position of each position
        for i = 1:n
            d = sign(round(A(i) * x + B(i) * y + C(i)));
            otherside = -sign(A(i) * px(i+2) + B(i) * py(i+2) + C(i));
            mini(d == otherside) = ISOUT;
            mini(d == 0 & mini == ISIN) = ISEDGE;
        end
    % draw polygon onto $pic$: use $orig$ to hold contents within bounding rect
        orig = pic(yy, xx); % copy background
        if length(interior)==1, orig(mini == ISIN) = interior; end
        if length(edge)==1, orig(mini == ISEDGE) = edge; end
        pic(yy, xx) = orig; % write out results
    toc

% verify answers are equal

    if ~isequal(pic, pic2)
        disagree = pic;
        disagree(pic == pic2) = ' ';
        disagree(pic ~= pic2) = 'X';
        pic, pic2, disagree, warning('solutions to $polygon$ disagree');
    end