Finding all pixels at least partially within an arbitrarily oriented rectangle

I have a rectangle with real-valued vertices (x0,y0), (x1,y1), (x2,y2), (x3,y3), which may be oriented at any angle in the plane. I'm looking for an efficient way of finding (or iterating over) all pixels (i.e., 1x1 squares) which 开发者_开发百科are at least partially within this rectangle.

It's trivial to do this for rectangles which are oriented orthogonally, and it's also trivial to check whether any particular pixel is within the rectangle. I could check each pixel within the rectangle's bounding box, but in the worst case I would be checking O(n^2) pixels when only O(n) will be within the target rectangle. (This is when the target rectangle is at 45 degrees and is very long and narrow.)

You can compute the range in the x-direction (floor of the minimum x-coordinate, to the ceiling of the maximum x-coordinate). For each x in that range, you can compute the range in the y-direction. You have a few different cases to consider in the generic case, depending on how the rectangle is oriented.

In essence, you have one leftmost point, one rightmost point, one upper point and one lower point. y1 will start at the leftmost, go trough the lower, and end in the rightmost point. y2 will instead go trough the upper point.

To include all touching pixels, we need to look half a pixel in all directions. I chose to use the center of each pixel as the coordinates. This was so that you get a more natural look of the final image.

Here are some F#-code to demonstrate:

let plot_rectangle p0 p1 p2 p3 =
    seq {
        // sort by x-coordinate
        let points = List.sortBy fst [p0; p1; p2; p3]
        let pLeft, pMid1, pMid2, pRight =
            points.[0], points.[1], points.[2], points.[3]

        // sort 2 middle points by y-coordinate
        let points = List.sortBy snd [pMid1; pMid2]
        let pBottom, pTop = points.[0], points.[1]

        // Easier access to the coordinates
        let pLeftX, pLeftY = pLeft
        let pRightX, pRightY = pRight
        let pBottomX, pBottomY = pBottom
        let pTopX, pTopY = pTop
        let pMid1X, pMid1Y = pMid1
        let pMid2X, pMid2Y = pMid2

        // Function: Get the minimum Y for a given X
        let getMinY x0 y0 x1 y1 x =
            let slope = (y1 - y0)/(x1 - x0)
            // Step half a pixel left or right, but not too far
            if slope >= 0.0 then
                let xl = max x0 (x - 0.5)
                y0 + slope * (xl - x0)
                |> round
                |> int
            else
                let xr = min x1 (x + 0.5)
                y0 + slope * (xr - x0)
                |> round
                |> int

        // Function: Get the maximum Y for a given X
        let getMaxY x0 y0 x1 y1 x =
            let slope = (y1 - y0)/(x1 - x0)
            // Step half a pixel left or right, but not too far
            if slope >= 0.0 then
                let xr = min x1 (x + 0.5)
                y0 + slope * (xr - x0)
                |> round
                |> int
            else
                let xl = max x0 (x - 0.5)
                y0 + slope * (xl - x0)
                |> round
                |> int

        let x1 = int (pLeftX + 0.5)
        let x2 = int (pRightX + 0.5)
        for x = x1 to x2 do
            let xf = float x
            if xf < pMid1X then
                // Phase I: Left to Top and Bottom
                // Line from pLeft to pBottom
                let y1 = getMinY pLeftX pLeftY pBottomX pBottomY xf
                // Line from pLeft to pTop
                let y2 = getMaxY pLeftX pLeftY pTopX pTopY xf
                for y = y1 to y2 do
                    yield (x, y)

            elif xf < pMid2X && pMid1Y < pMid2Y then
                // Phase IIa: left/bottom --> top/right
                // Line from pBottom to pRight
                let y1 = getMinY pBottomX pBottomY pRightX pRightY xf
                // Line from pLeft to pTop (still)
                let y2 = getMaxY pLeftX pLeftY pTopX pTopY xf
                for y = y1 to y2 do
                    yield (x, y)

            elif xf < pMid2X && pMid1Y >= pMid2Y then
                // Phase IIb: left/top --> bottom/right
                // Line from pLeft to pBottom (still)
                let y1 = getMinY pLeftX pLeftY pBottomX pBottomY xf
                // Line from pTop to pRight
                let y2 = getMaxY pTopX pTopY pRightX pRightY xf
                for y = y1 to y2 do
                    yield (x, y)

            else
                // Phase III: bottom/top --> right
                // Line from pBottom to pRight
                let y1 = getMinY pBottomX pBottomY pRightX pRightY xf
                // Line from pTop to pRight
                let y2 = getMaxY pTopX pTopY pRightX pRightY xf
                for y = y1 to y2 do
                    yield (x, y)
    }

Example:

Could you use something like a Graham scan?
You could use the set of 5 points (the pixel + the 4 vertexes) and then check to see whether the 4 vertexes define the boundary of the convex hull. This would be at worst O(n log n), which is a marked improvement on your n^2 for large n. Alternatively, a two-dimensional range tree might suffice, though I think this will still be n log n

EDIT: Actually, you could use the angles between the 4 vertexes to create 4 "ranges" where pixels could potentially be located, then just take the intersection of these 4 ranges. That would be a constant time operation, and checking whether a pixel lies within this range is also constant time - just compare the angle it makes with each vertex against the above set of angles.
As another alternative, use the 4 boundary lines (the lines between adjacent vertexes) and just 'walk' between them. Once you hit the line, any further points downwards won't lie within this boundary, etc. That's O(n) on the amount of pixels that lie within the rectangle, and should easily be solved with a trivial breadth-first search

Here's some Python code based on MizardX's answer which does exactly what I was wanting:

#!/usr/bin/python

import math

def minY(x0, y0, x1, y1, x):
  if x0 == x1:
    # vertical line, y0 is lowest
    return int(math.floor(y0))

  m = (y1 - y0)/(x1 - x0)

  if m >= 0.0:
    # lowest point is at left edge of pixel column
    return int(math.floor(y0 + m*(x - x0)))
  else:
    # lowest point is at right edge of pixel column
    return int(math.floor(y0 + m*((x + 1.0) - x0)))

def maxY(x0, y0, x1, y1, x):
  if x0 == x1:
    # vertical line, y1 is highest
    return int(math.ceil(y1))

  m = (y1 - y0)/(x1 - x0)

  if m >= 0.0:
    # highest point is at right edge of pixel column
    return int(math.ceil(y0 + m*((x + 1.0) - x0)))
  else:
    # highest point is at left edge of pixel column
    return int(math.ceil(y0 + m*(x - x0)))


# view_bl, view_tl, view_tr, view_br are the corners of the rectangle
view_bl = (0.16511327500123524, 1.2460844930844697)
view_tl = (1.6091354363329917, 0.6492542948962687)
view_tr = (1.1615128085358943, -0.4337622756706583)
view_br = (-0.2825093527958621, 0.16306792251754265)

pixels = []

# find l,r,t,b,m1,m2
view = [ view_bl, view_tl, view_tr, view_br ]

l, m1, m2, r = sorted(view, key=lambda p: (p[0],p[1]))
b, t = sorted([m1, m2], key=lambda p: (p[1],p[0]))

lx, ly = l
rx, ry = r
bx, by = b
tx, ty = t
m1x, m1y = m1
m2x, m2y = m2

xmin = 0
ymin = 0
xmax = 10
ymax = 10

# outward-rounded integer bounds
# note that we're clamping the area of interest to (xmin,ymin)-(xmax,ymax)
lxi = max(int(math.floor(lx)), xmin)
rxi = min(int(math.ceil(rx)), xmax)
byi = max(int(math.floor(by)), ymin)
tyi = min(int(math.ceil(ty)), ymax)

x1 = lxi 
x2 = rxi 

for x in range(x1, x2):
  xf = float(x)

  if xf < m1x:
    # Phase I: left to top and bottom
    y1 = minY(lx, ly, bx, by, xf)
    y2 = maxY(lx, ly, tx, ty, xf)

  elif xf < m2x:
    if m1y < m2y:
      # Phase IIa: left/bottom --> top/right
      y1 = minY(bx, by, rx, ry, xf)
      y2 = maxY(lx, ly, tx, ty, xf)

    else:
      # Phase IIb: left/top --> bottom/right
      y1 = minY(lx, ly, bx, by, xf)
      y2 = maxY(tx, ty, rx, ry, xf)

  else:
    # Phase III: bottom/top --> right
    y1 = minY(bx, by, rx, ry, xf)
    y2 = maxY(tx, ty, rx, ry, xf)

  y1 = max(y1, byi)
  y2 = min(y2, tyi)

  for y in range(y1, y2):
    pixels.append((x,y))

print pixels

Output:

[(0, 0), (0, 1), (1, 0)]