How to use Bresenham's line drawing algorithm with sub pixel bias?

Question

Bresenham's line drawing algorithm is well known and quite simple to implement.

While there are more advanced ways to draw anti-ailesed lines, Im interested in writing a function which draws a single pixel width non anti-aliased line, based on floating point coordinates.

This means while the first and last pixels will remain the same, the pixels drawn between them will have a bias based on the sub-pixel position of both end-points.

In principle this shouldn't be all that complicated, since I assume its possible to use the sub-pixel offsets to calculate an initial error value to use when plotting the line, and all other parts of the algorithm remain the same.

No sub pixel offset:

X###
    ###X

Assuming the right hand point has a sub-pixel position close to the top, the line could look like this:

With sub pixel offset for example:

X######
       X

Is there a tried & true method of drawing a line that takes sub-pixel coordinates into account?

Note:

This seems like a common operation, I've seen OpenGL drivers take this into account for example - using GL_LINE, though from a quick search I didn't find any answers online - maybe used wrong search terms?
At a glance this question looks like it might be a duplicate of:
Precise subpixel line drawing algorithm (rasterization algorithm)
However that is asking about drawing a wide line, this is asking about offsetting a single pixel line.
If there isn't some standard method, I'll try write this up to post as an answer.

I suppose might apply a fixed-point scaling factor and evalute the initial error term from the fractional part of the starting point. I can't see the benefit of Bresenham over easier algorithms though, say directly evaluating the function or a fixed-point DDA? Optimizing for the setup seems unlikely to be a win unless your average lines are tiny and I'd be hard-pressed to avoid taking a division in computing the error seed in any event. — doynax
While fixed-point DDA is probably a good solution, I'd like to be satisfied that applying bias to Bresenham's method as described in the question - is for some reason impractical. I can't really make too many assumptions about how many and what length the lines might be, but I would like to keep roughly the same performance, apart from a small speed hit from the initial setup. — ideasman42
I guess you'd seed with e = frac(x1) * (y2 - y1) + frac(y1) * (x2 - x1), albeit with the coordinates in fixed-point to avoid the rounding errors. My point is that the DDA innerloop (basically an addition and a shift) is usually faster than Bresenham, with the drawback being a small initialization hit which is typically irrelevant unless your lines are very short or division exceedingly expensive. — doynax

Dave Dave · Accepted Answer · 2017-01-25T23:58:13

Having just encountered the same challenge, I can confirm that this is possible as you expected.

First, return to the simplest form of the algorithm: (ignore the fractions; they'll disappear later)

x = x0
y = y0
dx = x1 - x0
dy = y1 - y0
error = -0.5
while x < x1:
    if error > 0:
        y += 1
        error -= 1
    paint(x, y)
    x += 1
    error += dy/dx

This means that for integer coordinates, we start half a pixel above the pixel boundary (error = -0.5), and for each pixel we advance in x, we increase the ideal y coordinate (and therefore the current error) by dy/dx.

First let's see what happens if we stop forcing x0, y0, x1 and y1 to be integers: (this will also assume that instead of using pixel centres, the coordinates are relative to the bottom-left of each pixel¹, since once you support sub-pixel positions you can simply add half the pixel width to the x and y to return to pixel-centred logic)

x = x0
y = y0
dx = x1 - x0
dy = y1 - y0
error = (0.5 - (x0 % 1)) * dy/dx + (y0 % 1) - 1
while x < x1:
    if error > 0:
        y += 1
        error -= 1
    paint(x, y)
    x += 1
    error += dy/dx

The only change was the initial error calculation. The new value comes from simple trig to calculate the y coordinate when x is at the pixel centre. It's worth noting that you can use the same idea to clip the line's start position to be within some bound, which is another challenge you'll likely face when you want to start optimising things.

Now we just need to convert this into integer-only arithmetic. We'll need some fixed multiplier for the fractional inputs (scale), and the divisions can be handled by multiplying them out, just as the standard algorithm does.

# assumes x0, y0, x1 and y1 are pre-multiplied by scale
x = x0
y = y0
dx = x1 - x0
dy = y1 - y0
error = (scale - 2 * (x0 % scale)) * dy + 2 * (y0 % scale) * dx - 2 * dx * scale
while x < x1:
    if error > 0:
        y += scale
        error -= 2 * dx * scale
    paint(x / scale, y / scale)
    x += scale
    error += 2 * dy * scale

Note that x, y, dx and dy keep the same scaling factor as the input variables (scale), whereas error has a more complex scaling factor: 2 * dx * scale. This allows it to absorb the division and fraction in its original formulation, but means we need to apply the same scale everywhere we use it.

Obviously there's a lot of room to optimise here, but that's the basic algorithm. If we assume scale is a power-of-two (2^n), we can start to make things a little more efficient:

dx = x1 - x0
dy = y1 - y0
mask = (1 << n) - 1
error = (2 * (y0 & mask) - (2 << n)) * dx - (2 * (x0 & mask) - (1 << n)) * dy
x = x0 >> n
y = y0 >> n
while x < (x1 >> n):
    if error > 0:
        y += 1
        error -= 2 * dx << n
    paint(x, y)
    x += 1
    error += 2 * dy << n

As with the original, this only works in the (x >= y, x > 0, y >= 0) octant. The usual rules apply for extending it to all cases, but note that there are a few extra gotchyas due to the coordinates no-longer being centred in the pixel (i.e. reflections become more complex).

You'll also need to watch out for integer overflows: error has twice the precision of the input variables, and a range of up to twice the length of the line. Plan your inputs, precision, and variable types accordingly!

^{1: Coordinates are relative to the corner which is closest to 0,0. For an OpenGL-style coordinate system that's the bottom left, but it could be the top-left depending on your particular scenario.}

How to use Bresenham's line drawing algorithm with sub pixel bias?

3 Answers

When exactly is a pixel drawn?