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.

Answer 1

Let's assume you want to draw a line from P1 = (x1, y1) to P2 = (x2, y2) where all the numbers are floating point pixel coordinates.

1. Calculate the true pixel coordinates of P1 and P2 and paint them: P* = (round(x), round(y)).

2. If abs(x1* - x2*) <= 1 && abs(y1* - y2*) <= 1 then you are finished.

3. Decide whether it is a horizontal (true) or a vertical line (false): abs(x1 - x2) >= abs(y1 - y2).

4. If it is a horizontal line and x1 > x2 or if it is a vertical line and y1 > y2: swap P1 with P2 (and also P1* with P2*).

5. If it is a horizontal line you can get the y-coordinates for all the x-coordinates between x1* and x2* with the following formula:

    y(x) = round(y1 + (x - x1) / (x2 - x1) * (y2 - y1))
   

   If you have a vertical line you can get the x-coordinates for all the y-coordinates between y1* and y2* with this formula:

    x(y) = round(x1 + (y - y1) / (y2 - y1) * (x2 - x1))
   

Here is a demo you can play around with, you can try different points on line 12.

Answer 2

I had a similar problem, with the addition of needing sub-pixel endpoints, I also needed to make sure all pixels which intersect the line are drawn.

I'm not sure that my solution will be helpful to OP, both because its been 4+ years, and because of the sentence "This means while the first and last pixels will remain the same..." For me, that is actually a problem (More on that later). Hopefully this may be helpful to others.

---

I don't know if this can be considered to be Bresenham's algorithm, but it is awful similar. I'll explain it for the (+,+) quadrant. Lets say you wish to draw a line from point (P_x,P_y) to (Q_x,Q_y) over a grid of pixels with width W. Having a grid width W > 1 allows for sub-pixel endpoints.

For a line going in the (+,+) quadrant, the starting point is easy to calculate, just take the floor of (P_x,P_y). As you will see later, this only works if Q_x >= P_x & Q_y >= P_y.

Now you need to find which pixel to go to next. There are 3 possibilities: (x+1,y), (x,y+1), & (x+1,y+1). To make this decision, I use the 2D cross product defined as:

• If this value is negative, vector b is right/clockwise of vector a.
• If this value is positive, vector b is left/anti-clockwise of vector a.
• If this value is zero vector b points in the same direction as vector a.

To make the decision on which pixel is next, compare the cross product between the line P-Q [red in image below] and a line between the point P and the top-right pixel (x+1,y+1) [blue in image below].

The vector between P & the top-right pixel can be calculated as:

So, we will use the value from the 2D cross product:

• If this value is negative, the next pixel will be (x,y+1).
• If this value is positive, the next pixel will be (x+1,y).
• If this value is exactly zero, the next pixel will be (x+1,y+1).

---

That works fine for the starting pixel, but the rest of the pixels will not have a point that lies inside them. Luckily, after the initial point, you don't need a point to be inside the pixel for the blue vector. You can keep extending it like so: The blue vector starts at the starting point of the line, and is updated to the (x+1,y+1) for every pixel. The rule for which pixel to take is the same. As you can see, the red vector is right of the blue vector. So, the next pixel will be the one right of the green pixel.

The value for the cross product needs updated for every pixel, depending on which pixel you took.

Add d_x if the next pixel was (x+1), add d_y if the pixel was (y+1). Add both if the pixel went to (x+1,y+1).

This process is repeated until it reaches the ending pixel, (Q_x / W, Q_y / W).

All combined this leads to the following code:

    int dx = x2 - x2;
    int dy = y2 - y1;
    int local_x = x1 % width;
    int local_y = y1 % width;
    int cross_product = dx*(width-local_y) - dy*(width-local_x);
    int dx_cross = -dy*width;
    int dy_cross = dx*width;

    int x = x1 / width;
    int y = y1 / width;
    int end_x = x2 / width;
    int end_y = y2 / width;
    while (x != end_x || y != end_y) {
        SetPixel(x,y,color);
        int old_cross = cross_product;
        if (old_cross >= 0) {
            x++;
            cross_product += dx_cross;
        }
        if (old_cross <= 0) {
            y++;
            cross_product += dy_cross;
        }
    }

Making it work for all quadrants is a matter of reversing the local coordinates and some absolute values. Heres the code which works for all quadrants:

    int dx = x2 - x1;
    int dy = y2 - y1;
    int dx_x = (dx >= 0) ? 1 : -1;
    int dy_y = (dy >= 0) ? 1 : -1;
    int local_x = x1 % square_width;
    int local_y = y1 % square_width;
    int x_dist = (dx >= 0) ? (square_width - local_x) : (local_x);
    int y_dist = (dy >= 0) ? (square_width - local_y) : (local_y);
    int cross_product = abs(dx) * abs(y_dist) - abs(dy) * abs(x_dist);
    dx_cross = -abs(dy) * square_width;
    dy_cross = abs(dx) * square_width;

    int x = x1 / square_width;
    int y = y1 / square_width;
    int end_x = x2 / square_width;
    int end_y = y2 / square_width;

    while (x != end_x || y != end_y) {
        SetPixel(x,y,color);
        int old_cross = cross_product;
        if (old_cross >= 0) {
            x += dx_x;
            cross_product += dx_cross;
        }
        if (old_cross <= 0) {
            y += dy_y;
            cross_product += dy_cross;
        }
    }

However there is a problem! This code will not stop in some cases. To understand why, you need to really look into exactly what conditions count as the intersection between a line and a pixel.

When exactly is a pixel drawn?

I said I need to make that all pixels which intersect a line need to be drawn. But there's some ambiguity in the edge cases.

Here is a list of all possible intersections in which a pixel will be drawn for a line where Q_x >= P_x & Q_y >= P_y:

• A - If a line intersects the pixel completely, the pixel will be drawn.
• B - If a vertical line intersects the pixel completely, the pixel will be drawn.
• C - If a horizontal line intersects the pixel completely, the pixel will be drawn.
• D - If a vertical line perfectly touches the left of the pixel, the pixel will be drawn.
• E - If a horizontal line perfectly touches the bottom of the pixel, the pixel will be drawn.
• F - If a line endpoint starts inside of a pixel going (+,+), the pixel will be drawn.
• G - If a line endpoint starts exactly on the left side of a pixel going (+,+), the pixel will be drawn.
• H - If a line endpoint starts exactly on the bottom side of a pixel going (+,+), the pixel will be drawn.
• I - If a line endpoint starts exactly on the bottom left corner of a pixel going (+,+), the pixel will be drawn.

---

And here are some pixels which do NOT intersect the line:

• A' - If a line obviously doesn't intersect a pixel, the pixel will NOT be drawn.

• B' - If a vertical line obviously doesn't intersect a pixel, the pixel will NOT be drawn.

• C' - If a horizontal line obviously doesn't intersect a pixel, the pixel will NOT be drawn.

• D' - If a vertical line exactly touches the right side of a pixel, the pixel will NOT be drawn.

• E' - If a horizontal line exactly touches the top side of a pixel, the pixel will NOT be drawn.

• F' - If a line endpoint starts exactly on the top right corner of a pixel going in the (+,+) direction, the pixel will NOT be drawn.

• G' - If a line endpoint starts exactly on the top side of a pixel going in the (+,+) direction, the pixel will NOT be drawn.

• H' - If a line endpoint starts exactly on the right side of a pixel going in the (+,+) direction, the pixel will NOT be drawn.

• I' - If a line exactly touches a corner of the pixel, the pixel will NOT be drawn. This applies to all corners.

---

Those rules apply as you would expect (just flip the image) for the other quadrants. The problem I need to highlight is when an endpoint lies exactly on the edge of a pixel. Take a look at this case:

This is like image G' above, except the y-axis is flipped because the Q_y < P_y. There are 4x4 red dots because W is 4, making the pixel dimensions 4x4. Each of the 4 dots are the ONLY endpoints a line can touch. The line drawn goes from (1.25, 1.0) to (somewhere).

This shows why it's incorrect (at least how I defined pixel-line intersections) to say the pixel endpoints can be calculated as the floor of the line endpoints. The floored pixel coordinate for that endpoint seems to be (1,1), but it is clear that the line never really intersects that pixel. It just touches it, so I don't want to draw it.

Instead of flooring the line endpoints, you need to floor the minimal endpoints, and ceil the maximal endpoints minus 1 across both x & y dimensions.

So finally here is the complete code which does this flooring/ceiling:

    int dx = x2 - x1;
    int dy = y2 - y1;
    int dx_x = (dx >= 0) ? 1 : -1;
    int dy_y = (dy >= 0) ? 1 : -1;
    int local_x = x1 % square_width;
    int local_y = y1 % square_width;
    int x_dist = (dx >= 0) ? (square_width - local_x) : (local_x);
    int y_dist = (dy >= 0) ? (square_width - local_y) : (local_y);
    int cross_product = abs(dx) * abs(y_dist) - abs(dy) * abs(x_dist);
    dx_cross = -abs(dy) * square_width;
    dy_cross = abs(dx) * square_width;

    int x = x1 / square_width;
    int y = y1 / square_width;
    int end_x = x2 / square_width;
    int end_y = y2 / square_width;

    // Perform ceiling/flooring of the pixel endpoints
    if (dy < 0)
    {
        if ((y1 % square_width) == 0)
        {
            y--;
            cross_product += dy_cross;
        }
    }
    else if (dy > 0)
    {
        if ((y2 % square_width) == 0)
            end_y--;
    }

    if (dx < 0)
    {
        if ((x1 % square_width) == 0)
        {
            x--;
            cross_product += dx_cross;
        }
    }
    else if (dx > 0)
    {
        if ((x2 % square_width) == 0)
            end_x--;
    }

    while (x != end_x || y != end_y) {
        SetPixel(x,y,color);
        int old_cross = cross_product;
        if (old_cross >= 0) {
            x += dx_x;
            cross_product += dx_cross;
        }
        if (old_cross <= 0) {
            y += dy_y;
            cross_product += dy_cross;
        }
    }

This code itself hasn't been tested, but it comes slightly modified from my GitHub project where it has been tested.

Answer 3

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