How to use Bresenham's line drawing algorithm with clipping?

Question

When drawing a line with Bresenham line drawing algorithm, where the line may not be within the bounds of the bitmap being written to - it would be useful to clip the results so they fit within the axis aligned bounds of the image being written to.

While its possible to first clip the line to the rectangle, then draw the line. This isn't ideal since it often gives a slightly different slant to the line (assuming int coords are being used).

Since this is such a primitive operation, are there established methods for clipping the line while maintaining the same shape?

In case it helps, here is a reference implementation of the algorithm - it uses int coords, which avoids int/float conversion while drawing the line.

---

I spent some time looking into this:

• Problem is described in detail on virtual-dub's web page.
• Possible solution from Alan Tiedemann
  ... though I'd need to implement code based on the text - to see how well it works.
• Bresenham's Line Generation Algorithm with Built-in Clipping
  (paper from 1995, couldn't find the entire document? - PDF is a single page without the C code it references, looks like its pay-walled).

Manx

Clipping shouldnt change the slant if you use floating point coordinates. Are you looking for a pixel perfect match? – Afterburning

I'd like to be able to keep using int coords if possible to avoid a lot of float/int conversion, since the current code is very efficient with int coords. I assume a paper on this topic was released because simply using float coords isn't as efficient - added a link to reference code in the question. – Manx

I've read the first and second of your reference links and I have to admit to being slightly confused as to what's wrong with the 'cheesy' method. Your bounds are axis aligned; the coordinates of the pixels being set change monotonically; you know exactly when to switch from 'don't put pixel' to 'put pixel' to 'don't put pixel'; not sure where the problem is – Kerley

I don't have a full answer with code, but you should be able to multiply by the gradient as a rational, and use the remainder as the initial delta. IOW, if your line starts x0 pixels to the left of the bounds and has positive gradient dy/dx, then you use x0 * dy (widening e.g. from int to long may be needed). Then divide by dx to get the y value for x0 and use the remainder to initialise delta for that position. Does that help? – Univalence

IMO the trick is to first intersect the line with the bounding box including the accumulated error term, and run the Bresenham between the intersection points. – Backstop

@Backstop "to first intersect the line with the bounding box including the accumulated error term" doesn't this mean "compute the full actual line but only turn on the pixels inside the clipping area"? – Aguirre

@AdrianColomitchi not really, calculating error term is much simpler than computing every point and ignore the ones that are not in the clip region. – Hummel

No, it means: intersecting two lines (the intended line plus one of the edges of the BB). Linear algebra. For the north-edge you know y (exact), and compute x + x_debt. similar for other intersections. – Backstop

Answer 1

Let's reframe the problem so we can see how Bresenham's algorithm really works...

Lets say you are drawing a mostly horizontal line (the method is the same for mostly vertical, but with the axes switched) from (x0,y0) to (x1,y1):

The full line can be described as a function of y in terms of x (all integers):

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

This describes precisely each pixel you will paint when drawing the full line, and when you clip the line consistently, it still describes precisely each pixel you will paint -- you just apply it to a smaller range of x values.

We can evaluate this function using all integer math, by calculating the divisor and remainder separately. For x1 >= x0 and y1 >= y0 (do the normal transformations otherwise):

let dx = (x1-x0);
let dy = (y1-y0);
let remlimit = (dx+1)/2; //round up

rem = (x-x0) * dy % dx;
y = y0 + (x-x0) * dy / dx;
if (rem >= remlimit)
{
    rem-=dx;
    y+=1;
}

Bresenham's algorithm is a just a fast way to update the result of this formula incrementally as you update x.

Here's how we can make use of incremental updates to draw the portion of the very same line from x=xs to x=xe:

let dx = (x1-x0);
let dy = (y1-y0);
let remlimit = (dx+1)/2; //round up

x=xs;
rem = (x-x0) * dy % dx;
y = y0 + (x-x0) * dy / dx;
if (rem >= remlimit)
{
    rem-=dx;
    y+=1;
}
paint(x,y);
while(x < xe)
{
    x+=1;
    rem+=dy;
    if (rem >= remlimit)
    {
        rem-=dx;
        y+=1;
    }
    paint(x,y);
}

If you want do to your remainder comparisons against 0, you can just offset it at the beginning:

let dx = (x1-x0);
let dy = (y1-y0);
let remlimit = (dx+1)/2; //round up

x=xs;
rem = ( (x-x0) * dy % dx ) - remlimit;
y = y0 + (x-x0) * dy / dx;
if (rem >= 0)
{
    rem-=dx;
    y+=1;
}
paint(x,y);
while(x < xe)
{
    x+=1;
    rem+=dy;
    if (rem >= 0)
    {
        rem-=dx;
        y+=1;
    }
    paint(x,y);
}

Swum

While this approach seems correct, I've found an implementation from Kuzmin & Yevgeny P's paper, and its quite a bit more involved then this answer makes out - although its mostly just checking corner cases and accounting for both axis. – Manx

Sorry @ideasman42, it's just not that hard. Well, you have to find xs and xe for your clipping rectangle -- easy for the x bounds but a bit finicky for y bounds due to rounding. if you have trouble with that let me know and I'll show the calculation. – Swum

Right, its not hard to calculate spesific values - however I was attempting to have the clipping and non-clipping versions give exactly the same results for the visible line segment. Turns out there is some difference in the method described by the paper by Kuzmin & Yevgeny P. (unrelated to clipping) I mistook as errors in the algorithm, treating slopes slightly differently - but seems not to be a bug just a very minor difference. I've posted the code as a separate answer. – Manx

@MattTimmermans I've made several attempts at this problem on and off over the past ten years or so. No matter how much energy I have going in, I'm defeated every time. It is definitely hard. It's easier to just not use Bresenham's wretched little algorithm, and instead rasterize in float coords. – Eurythermal

Answer 2

Bresenham's algorithm can be used, taking clipping values into account, based on the paper by Kuzmin & Yevgeny P:

For completeness, here is a working version of the algorithm, a single Python function, though its only using integer arithmetic - so can be easily ported to other languages.

def plot_line_v2v2i(
    p1, p2, callback,
    clip_xmin, clip_ymin,
    clip_xmax, clip_ymax,
):
    x1, y1 = p1
    x2, y2 = p2

    del p1, p2

    # Vertical line
    if x1 == x2:
        if x1 < clip_xmin or x1 > clip_xmax:
            return

        if y1 <= y2:
            if y2 < clip_ymin or y1 > clip_ymax:
                return
            y1 = max(y1, clip_ymin)
            y2 = min(y2, clip_ymax)
            for y in range(y1, y2 + 1):
                callback(x1, y)
        else:
            if y1 < clip_ymin or y2 > clip_ymax:
                return
            y2 = max(y2, clip_ymin)
            y1 = min(y1, clip_ymax)
            for y in range(y1, y2 - 1, -1):
                callback(x1, y)
        return

    # Horizontal line
    if y1 == y2:
        if y1 < clip_ymin or y1 > clip_ymax:
            return

        if x1 <= x2:
            if x2 < clip_xmin or x1 > clip_xmax:
                return
            x1 = max(x1, clip_xmin)
            x2 = min(x2, clip_xmax)
            for x in range(x1, x2 + 1):
                callback(x, y1)
        else:
            if x1 < clip_xmin or x2 > clip_xmax:
                return
            x2 = max(x2, clip_xmin)
            x1 = min(x1, clip_xmax)
            for x in range(x1, x2 - 1, -1):
                callback(x, y1)
        return

    # Now simple cases are handled, perform clipping checks.
    if x1 < x2:
        if x1 > clip_xmax or x2 < clip_xmin:
            return
        sign_x = 1
    else:
        if x2 > clip_xmax or x1 < clip_xmin:
            return
        sign_x = -1

        # Invert sign, invert again right before plotting.
        x1 = -x1
        x2 = -x2
        clip_xmin, clip_xmax = -clip_xmax, -clip_xmin

    if y1 < y2:
        if y1 > clip_ymax or y2 < clip_ymin:
            return
        sign_y = 1
    else:
        if y2 > clip_ymax or y1 < clip_ymin:
            return
        sign_y = -1

        # Invert sign, invert again right before plotting.
        y1 = -y1
        y2 = -y2
        clip_ymin, clip_ymax = -clip_ymax, -clip_ymin

    delta_x = x2 - x1
    delta_y = y2 - y1

    delta_x_step = 2 * delta_x
    delta_y_step = 2 * delta_y

    # Plotting values
    x_pos = x1
    y_pos = y1

    if delta_x >= delta_y:
        error = delta_y_step - delta_x
        set_exit = False

        # Line starts below the clip window.
        if y1 < clip_ymin:
            temp = (2 * (clip_ymin - y1) - 1) * delta_x
            msd = temp // delta_y_step
            x_pos += msd

            # Line misses the clip window entirely.
            if x_pos > clip_xmax:
                return

            # Line starts.
            if x_pos >= clip_xmin:
                rem = temp - msd * delta_y_step

                y_pos = clip_ymin
                error -= rem + delta_x

                if rem > 0:
                    x_pos += 1
                    error += delta_y_step
                set_exit = True

        # Line starts left of the clip window.
        if not set_exit and x1 < clip_xmin:
            temp = delta_y_step * (clip_xmin - x1)
            msd = temp // delta_x_step
            y_pos += msd
            rem = temp % delta_x_step

            # Line misses clip window entirely.
            if y_pos > clip_ymax or (y_pos == clip_ymax and rem >= delta_x):
                return

            x_pos = clip_xmin
            error += rem

            if rem >= delta_x:
                y_pos += 1
                error -= delta_x_step

        x_pos_end = x2

        if y2 > clip_ymax:
            temp = delta_x_step * (clip_ymax - y1) + delta_x
            msd = temp // delta_y_step
            x_pos_end = x1 + msd

            if (temp - msd * delta_y_step) == 0:
                x_pos_end -= 1

        x_pos_end = min(x_pos_end, clip_xmax) + 1
        if sign_y == -1:
            y_pos = -y_pos
        if sign_x == -1:
            x_pos = -x_pos
            x_pos_end = -x_pos_end
        delta_x_step -= delta_y_step

        while x_pos != x_pos_end:
            callback(x_pos, y_pos)

            if error >= 0:
                y_pos += sign_y
                error -= delta_x_step
            else:
                error += delta_y_step

            x_pos += sign_x
    else:
        # Line is steep '/' (delta_x < delta_y).
        # Same as previous block of code with swapped x/y axis.

        error = delta_x_step - delta_y
        set_exit = False

        # Line starts left of the clip window.
        if x1 < clip_xmin:
            temp = (2 * (clip_xmin - x1) - 1) * delta_y
            msd = temp // delta_x_step
            y_pos += msd

            # Line misses the clip window entirely.
            if y_pos > clip_ymax:
                return

            # Line starts.
            if y_pos >= clip_ymin:
                rem = temp - msd * delta_x_step

                x_pos = clip_xmin
                error -= rem + delta_y

                if rem > 0:
                    y_pos += 1
                    error += delta_x_step
                set_exit = True

        # Line starts below the clip window.
        if not set_exit and y1 < clip_ymin:
            temp = delta_x_step * (clip_ymin - y1)
            msd = temp // delta_y_step
            x_pos += msd
            rem = temp % delta_y_step

            # Line misses clip window entirely.
            if x_pos > clip_xmax or (x_pos == clip_xmax and rem >= delta_y):
                return

            y_pos = clip_ymin
            error += rem

            if rem >= delta_y:
                x_pos += 1
                error -= delta_y_step

        y_pos_end = y2

        if x2 > clip_xmax:
            temp = delta_y_step * (clip_xmax - x1) + delta_y
            msd = temp // delta_x_step
            y_pos_end = y1 + msd

            if (temp - msd * delta_x_step) == 0:
                y_pos_end -= 1

        y_pos_end = min(y_pos_end, clip_ymax) + 1
        if sign_x == -1:
            x_pos = -x_pos
        if sign_y == -1:
            y_pos = -y_pos
            y_pos_end = -y_pos_end
        delta_y_step -= delta_x_step

        while y_pos != y_pos_end:
            callback(x_pos, y_pos)

            if error >= 0:
                x_pos += sign_x
                error -= delta_y_step
            else:
                error += delta_x_step

            y_pos += sign_y

Example use:

plot_line_v2v2i(
    # two points
    (10, 2),
    (90, 88),
    # callback
    lambda x, y: print(x, y),
    # xy min
    25, 25,
    # xy max
    75, 75,
)

Notes:

• Clipping min/max values are inclusive
  (so max values should be image_width - 1, image_height - 1)
• Integer divisions // is used everywhere.
• Many languages (C/C++ for example) use floored rounding on division.
  See Fast floor of a signed integer division in C / C++ to avoid having slightly biased results with those languages.

There are some improvements over the code provided in the paper:

• The line will always plot in the direction defined (from p1 to p2).
• There was sometimes a subtle difference in the line gradient, so that rotating of flipping the points, calculating the line, then transforming back - would give slightly different results. The asymmetry was caused by the code swapping the X and Y axis to avoid code duplication.

---

For tests and more example usage, see:

• The Python repository.
• Rust version.

Manx

I've converted it to C, works really nice. One obstacle was understanding that '//' is integer division rather than comment prefix ;) +Kudos ideasman42! – Angulation

@Angulation yw - note that int division in C will round differently based on sign. Will note in answer – Manx

OH! Python brought me down to the floor with that! Thanks! – Angulation

@Angulation can you share the C code? Also, I feel bad, I accidentally down-voted this answer before I realized what was going on. Now, SO won't let me upvote it :(. – Decor

I don't have a C version, the rust version should be fairly simple to port to C though. Also, I made some minor edits. Maybe you can change the vote now :) – Manx