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

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);
}

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.