function Rasterizer( vertices )
% vertices is a 2-by-3 matrix; 
% each column is the coordinates of a vertex

minX = floor(min(vertices(1,:)));
maxX = ceil(max(vertices(1,:)));
minY = floor(min(vertices(2,:)));
maxY = ceil(max(vertices(2,:)));

% change it to true to make it a normal rasterizer
normalMode = false;

if ~normalMode
    segments = [13,14,10,11,12,13,10,11,12,10,11,12,11;...
                22,22,23,23,23,23,24,24,24,25,25,25,26];
end

hold on

axis equal on
axis([minX maxX minY maxY])
set(gca, 'xtick', minX:maxX)
set(gca, 'ytick', minY:maxY)
grid on

plot([vertices(1,:), vertices(1,1)], [vertices(2,:), vertices(2,1)], '-k', 'LineWidth', 3);

S = Area(vertices(:,1), vertices(:,2), vertices(:,3));

for cx=minX:maxX
    for cy=minY:maxY
        pt = [cx, cy];
        inTriangle = Area(pt, vertices(:,1), vertices(:,2)) * S >= 0 ...
            && Area(pt, vertices(:,2), vertices(:,3)) * S >= 0 ...
            && Area(pt, vertices(:,3), vertices(:,1)) * S >= 0;
        if normalMode
            if inTriangle
                scatter(cx, cy, 'or', 'filled');
            else
                scatter(cx, cy, 'ob', 'filled');
            end
        else
            inSegments = sum(segments(1,:)==cx & segments(2,:)==cy) > 0;
            if inTriangle && inSegments
                scatter(cx, cy, 'or', 'filled');
            elseif inTriangle && ~inSegments
                scatter(cx, cy, 'xr', 'LineWidth', 3);
            elseif ~inTriangle && inSegments
                scatter(cx, cy, 'xb', 'LineWidth', 3);
            else
                scatter(cx, cy, 'ob', 'filled');
            end
        end
    end
end

if ~normalMode
    depths = vertices(1,:)*0.2+vertices(2,:)*0.3+10;
    title(sprintf(['Depth:\n\tz(%.3f,%.3f)=%.3f\n\tz(%.3f,%.3f)' ... 
        '=%.3f\n\tz(%.3f,%.3f)=%.3f'], ...
        vertices(1,1), vertices(2,1), depths(1), ...
        vertices(1,2), vertices(2,2), depths(2), ...
        vertices(1,3), vertices(2,3), depths(3)));
end

hold off

function S=Area(v0, v1, v2)
S = 0.5 * (v0(1)*v1(2)+v1(1)*v2(2)+v2(1)*v0(2) ...
    -v0(2)*v1(1)-v1(2)*v2(1)-v2(2)*v0(1));