Search code examples
c#algorithmrasterizingrange-tree

Algorithm for finding grid cells contained in arbitrary rotated rectangle (rasterizing)

I'm looking for an algorithm that can compute all grid cells occupied by an arbitrary rectangle in 2d space in a defined area. The rectangle is defined by it's 4 corner coordinates. In the picture below I marked two of them as red, the coordinates of the grid cells are in black. I'm looking for a generic case that also covers unaligned grids (center of grid != center of underlying coodinate system), but converting between Cell coordinates <=> Coordinate system is trivial and already implemented. All cells are squares.

The calculation is done thousands in very short amounts of time and needs to be as fast as possible.

enter image description here

What I'm doing right now:

Right now I'm calculating the edges of the rectangle (AB, BC, CD, DA) and sample them at sampleRate = cellWidth / 4 intervals. To find the middle cells I construct new lines from AB + sampleRate * x to CD + sampleRate * x and sample those lines again at sampleRate. I put all cells that I find at each sampled point into a HashSet to remove duplicates. But I feel this is increadibly inefficient.

I also thought about grabbing all cells in my relevant area into a buffer and generate a range tree or index tree. Then I could queue the bounds of the rectangle and get all contained cells at O(log n+m), but I'm not sure how to implement that and I can't even find any C# range tree implementations.

My knowledge in computer graphics is very rusty, but shouldn't there be an easy geometrical approach that works without sampling?

Solution

  • The points you care about are marked with a dot in the image below. Each point represents either the minimum X value or the maximum X value for a given Y value. For each Y value, the X value is easily calculated from the equation for a line:

    x = x0 + (y - y0)(x1 - x0) / (y1 - y0)

    Note that an axis-aligned rectangle (where (y1 - y0) is 0) needs to be handled as a special (but trivial) case.

    Also note that once you've computed the first X value along an edge of the rectangle, the other X values are equally spaced. The spacing is the inverse of the slope of the line. So the division
    M = (x1 - x0) / (y1 - y0)
    only needs to be done once, and then repeatedly added to the X value. For example, after calculating the X coordinate for point A in the image, the X coordinate for point B is Bx = Ax + M. And the X coordinate for point C is Cx = Bx + M.

    Once you have the min and max X values for each Y value, it's a simple for loop to get all the cells that have that Y value.

    enter image description here