Rotated rectangle rasterisation algorithm

Question

In a nutshell: I want to do a non-approximate version of Bresenham's line algorithm, but for a rectangle rather than a line, and whose points aren't necessarily aligned to the grid.

Given a square grid, and a rectangle comprising four non-grid-aligned points, I want to find a list of all grid squares that are covered, partially or completely, by the rectangle.

Bresenham's line algorithm is approximate – not all partially covered squares are identified. I'm looking for a "perfect" algorithm, that has no false positives or negatives.

Answer 1

It's an old question, but I have solved this issue (C++)

https://github.com/feelinfine/tracer

Maybe it will be usefull for someone

(sorry for my poor english)

Single line tracing

template <typename PointType>
std::set<V2i> trace_line(const PointType& _start_point, const PointType& _end_point, size_t _cell_size)
{   
    auto point_to_grid_fnc = [_cell_size](const auto& _point)
    {
        return V2i(std::floor((double)_point.x / _cell_size), std::floor((double)_point.y / _cell_size));
    };

    V2i start_cell = point_to_grid_fnc(_start_point);
    V2i last_cell = point_to_grid_fnc(_end_point);

    PointType direction = _end_point - _start_point;

    //Moving direction (cells)
    int step_x = (direction.x >= 0) ? 1 : -1;
    int step_y = (direction.y >= 0) ? 1 : -1;

    //Normalize vector
    double hypot = std::hypot(direction.x, direction.y);
    V2d norm_direction(direction.x / hypot, direction.y / hypot);

    //Distance to the nearest square side
    double near_x = (step_x >= 0) ? (start_cell.x + 1)*_cell_size - _start_point.x : _start_point.x - (start_cell.x*_cell_size);    
    double near_y = (step_y >= 0) ? (start_cell.y + 1)*_cell_size - _start_point.y : _start_point.y - (start_cell.y*_cell_size);

    //How far along the ray we must move to cross the first vertical (ray_step_to_vside) / or horizontal (ray_step_to_hside) grid line
    double ray_step_to_vside = (norm_direction.x != 0) ? near_x / norm_direction.x : std::numeric_limits<double>::max();
    double ray_step_to_hside = (norm_direction.y != 0) ? near_y / norm_direction.y : std::numeric_limits<double>::max();

    //How far along the ray we must move for horizontal (dx)/ or vertical (dy) component of such movement to equal the cell size
    double dx = (norm_direction.x != 0) ? _cell_size / norm_direction.x : std::numeric_limits<double>::max();
    double dy = (norm_direction.y != 0) ? _cell_size / norm_direction.y : std::numeric_limits<double>::max();

    //Tracing loop
    std::set<V2i> cells;
    cells.insert(start_cell);

    V2i current_cell = start_cell;

    size_t grid_bound_x = std::abs(last_cell.x - start_cell.x);
    size_t grid_bound_y = std::abs(last_cell.y - start_cell.y);

    size_t counter = 0;

    while (counter != (grid_bound_x + grid_bound_y))
    {
        if (std::abs(ray_step_to_vside) < std::abs(ray_step_to_hside))
        {
            ray_step_to_vside = ray_step_to_vside + dx; //to the next vertical grid line
            current_cell.x = current_cell.x + step_x;
        }
        else
        {
            ray_step_to_hside = ray_step_to_hside + dy;//to the next horizontal grid line
            current_cell.y = current_cell.y + step_y;
        }

        ++counter;

        cells.insert(current_cell);
    };

    return cells;
}

Get all cells

template <typename Container>
std::set<V2i> pick_cells(Container&& _points, size_t _cell_size)
{
    if (_points.size() < 2 || _cell_size <= 0)
        return std::set<V2i>();

    Container points = std::forward<Container>(_points);

    auto add_to_set = [](auto& _set, const auto& _to_append)
    {
        _set.insert(std::cbegin(_to_append), std::cend(_to_append));
    };

    //Outline
    std::set<V2i> cells;

    /*
    for (auto it = std::begin(_points); it != std::prev(std::end(_points)); ++it)
        add_to_set(cells, trace_line(*it, *std::next(it), _cell_size));
    add_to_set(cells, trace_line(_points.back(), _points.front(), _cell_size));
    */

    //Maybe this code works faster
    std::vector<std::future<std::set<V2i> > > results;

    using PointType = decltype(points.cbegin())::value_type;

    for (auto it = points.cbegin(); it != std::prev(points.cend()); ++it)           
        results.push_back(std::async(trace_line<PointType>, *it, *std::next(it), _cell_size));

    results.push_back(std::async(trace_line<PointType>, points.back(), points.front(), _cell_size));    

    for (auto& it : results)
        add_to_set(cells, it.get());

    //Inner
    std::set<V2i> to_add;

    int last_x = cells.begin()->x;
    int counter = cells.begin()->y;

    for (auto& it : cells)
    {
        if (last_x != it.x)
        {
            counter = it.y;
            last_x = it.x;
        }

        if (it.y > counter) 
        {
            for (int i = counter; i < it.y; ++i)
                to_add.insert(V2i(it.x, i));
        }

        ++counter;
    }

    add_to_set(cells, to_add);

    return cells;
}

Types

template <typename _T>
struct V2
{
    _T x, y;

    V2(_T _x = 0, _T _y = 0) : x(_x), y(_y)
    {
    };

    V2 operator-(const V2& _rhs) const
    {
        return V2(x - _rhs.x, y - _rhs.y);
    }

    bool operator==(const V2& _rhs) const
    {
        return (x == _rhs.x) && (y == _rhs.y);
    }

    //for std::set sorting
    bool operator<(const V2& _rhs) const
    {
        return (x == _rhs.x) ? (y < _rhs.y) : (x < _rhs.x);
    }
};

using V2d = V2<double>;
using V2i = V2<int>;

Usage

std::vector<V2d> points = { {200, 200}, {400, 400}, {500,100} };
size_t cell_size = 30;
auto cells = pick_cells(points, cell_size);
for (auto& it : cells)
    ...                 //do something with cells

Answer 2

You can use a scanline approach. The rectangle is a closed convex polygon, so it is sufficient to store the leftmost and rightmost pixel for each horizontal scanline. (And the top and bottom scanlines, too.)

The Bresenham algorithm tries to draw a thin, visually pleasing line without adjacent cells in the smaller dimension. We need an algorithm that visits each cell that the edges of the polygon pass through. The basic idea is to find the starting cell (x, y) for each edge and then to adjust x whenever the edge intersects a vertical border and to adjust y when it intersects a horizontal border.

We can represent the intersections by means of a normalised coordinate s that travels along the edge and that is 0.0 at the first node n1 and 1.0 at the second node n2.

    var x = Math.floor(n1.x / cellsize);
    var y = Math.floor(n1.y / cellsize);
    var s = 0;

The vertical insersections can the be represented as equidistant steps of with dsx from an initial sx.

    var dx = n2.x - n1.x;

    var sx = 10;            // default value > 1.0

    // first intersection
    if (dx < 0) sx = (cellsize * x - n1.x) / dx;
    if (dx > 0) sx = (cellsize * (x + 1) - n1.x) / dx;

    var dsx = (dx != 0) ? grid / Math.abs(dx) : 0;

Likewise for the horizontal intersecions. A default value greater than 1.0 catches the cases of horizontal and vertical lines. Add the first point to the scanline data:

    add(scan, x, y);

Then we can visit the next adjacent cell by looking at the next intersection with the smallest s.

    while (sx <= 1 || sy <= 1) {
        if (sx < sy) {
            sx += dsx;
            if (dx > 0) x++; else x--;
        } else {
            sy += dsy;
            if (dy > 0) y++; else y--;
        }
        add(scan, x, y);
    }

Do this for all four edges and with the same scanline data. Then fill all cells:

    for (var y in scan) {
        var x = scan[y].min;
        var xend = scan[y].max + 1;
        while (x < xend) {
            // do something with cell (x, y)
            x++;
        }
    }

(I have only skimmed the links MBo provided. It seems that the approach presented in that paper is essentially the same as mine. If so, please excuse the redundant answer, but after working this out I thought I could as well post it.)