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.)

Answer 3

This is sub-optimal but might give a general idea.

First off treat the special case of the rectangle being aligned horizontally or vertically separately. This is pretty easy to test for and make the rest simpler.

You can represent the rectangle as a set of 4 inequalities a1 x + b1 y >= c1 a1 x + b1 y <= c2 a3 x + b3 y >= c3 a3 x + b3 y <= c4 as the edges of the rectangles are parallel some of the constants are the same. You also have (up to a multiple) a3=b1 and b3=-a1. You can multiply each inequality by a common factor so you are working with integers.

Now consider each scan line with a fixed value of y. For each value of y find the four points where the lines intersect the scan line. That is find the solution with each line above. A little bit of logic will find the minimum and maximum values of x. Plot all pixels between these values.

You condition that you want all partially covered squares makes things a little trickier. You can solve this by considering two adjacent scan lines. You want to plot the points between the minimum x for both lines and the maximum for the both lines. If say a1 x+b1 y>=c is the inequality for the bottom left line in the figure. You want the find the largest x such that a1 x + b1 y < c this will be floor((c-b1 y)/a1) call this minx(y) also find minx(y+1) and the left hand point will be the minimum of these two values.

There is many easy optimisation you can find the y-values of the top and bottom corners reducing the range of y-values to test. You should only need to test two side. For each end point of each line there is one multiplication, one subtraction and one division. The division is the slowest part I think about 4 time slower than other ops. You might be able to remove this with a Bresenham or DDA algorithms others have mentioned.

Answer 4

There is method of Amanatides and Woo to enumerate all intersected cells A Fast Voxel Traversal Algorithm for Ray Tracing.
Here is practical implementation.
As side effect for you - you'll get points of intersection with grid lines - it may be useful if you need areas of partially covered cells (for antialiasing etc).