How to implement Bresenham’s algorithm as a system of equations?

Bresenham's line drawing algorithm is usually implemented via loops. But in Mathematica we can take advantage of its ability to solve Diophantine equations. From educational viewpoint it is quite interesting to write Bresenham in a form suitable for Solve or other superfunction. Such implementation also potentially can be very terse and sufficiently efficient for most practical applications.

I have implemented Bresenham for the first octant and obtained general solution via reflections (code for version 10):

Clear[bresenhamSolve];
bresenhamSolve[p1_, p2_] /;
GreaterEqual @@ Abs[p2 - p1] && MatchQ[Sign[p2 - p1], {1, 1} | {1, 0}] :=
Block[{ab = First@Solve[{p1, p2}.{a, 1} == b], x, y}, {x, y} /.
Solve[{a x + y + err == b /. ab, -1/2 < err <= 1/2, {x, y} \[Element] Integers,
p1 <= {x, y} <= p2}, {x, y, err}]
];
bresenhamSolve[p1_, p2_] /;
Less @@ Abs[p2 - p1] && MatchQ[Sign[p2 - p1], {1, 1} | {0, 1}] :=
Reverse /@ bresenhamSolve[Reverse[p1], Reverse[p2]];
bresenhamSolve[p1_, p2_] :=
With[{s = 2 UnitStep[p2 - p1] - 1},
Replace[bresenhamSolve[p1 s, p2 s], {x_, y_} :> s {x, y}, {1}]];

Of course this code cannot be called terse.

This implementation is identical to halirutan's implementation:

lines = DeleteCases[Partition[#, 2] & /@ Tuples[Range[-30, 30, 6], {4}], {p_, p_}];
Length[lines]
bresenhamSolve @@@ lines == bresenham @@@ lines
14520

True

Visualization:

p1 = {2, 3}; p2 = {20, 13};
Graphics[{EdgeForm[{Thick, RGBColor[203/255, 5/17, 22/255]}],
FaceForm[RGBColor[131/255, 148/255, 10/17]],
Rectangle /@ (bresenhamSolve[p1, p2] - .5), {RGBColor[0, 43/255, 18/85], Thick,
Line[{p1, p2}]}}, GridLines -> (Range[#1, #2 + 1] & @@@ Transpose[{p1, p2}] - .5),
Frame -> True] Manipulate[Row[{Graphics[{EdgeForm[{Thick, RGBColor[203/255, 5/17, 22/255]}],
FaceForm[RGBColor[131/255, 148/255, 10/17]],
Rectangle /@ (bresenhamSolve @@ Round[pts] - .5), {RGBColor[0, 43/255, 18/85], Thick,
Arrow@pts}}, GridLines -> ({Range[-50, 50], Range[-50, 50]} - .5), Frame -> True,
PlotRange -> {{-20, 20}, {-20, 20}}, ImageSize -> 500],
Column[Round@pts]}], {{pts, {{-11, -13}, {8, 15}}}, Locator}] Checking symmetry:

n = 20; center = {0, 0};
perimeterOfSquare = {x, y} /.
Solve[{x, y} \[Element]
RegionBoundary[Rectangle[{-n, -n} + center, {n, n} + center]], {x, y}, Integers];
ArrayPlot[SparseArray[
Rule @@@ Tally[# - center & /@
Flatten[bresenhamSolve[center, #] & /@ perimeterOfSquare + n + 1, 1]], {2 n + 1,
2 n + 1}], Mesh -> True, PlotRange -> {All, All, {1, 9}}, ClippingStyle -> Red,
PixelConstrained -> True] Timing comparison with halirutan's implementation:

p1 = {2, 3}; p2 = {2001, 1300};
Timing[bresenhamSolve[p1, p2];]
Timing[bresenham[p1, p2];]
{0.171601, Null}

{0.0312002, Null}

My question:

Is it possible to write general implementation without splitting it into special cases and with good performance? Instead of Solve one can use other superfunction(s).

• Is there a reason you want to use the Bresenham algorithm? Wouldn't it be much simpler and more efficient to use floating point arithmetic? I mean, the Bresenham algorithm was useful when floating point operations were expensive and it payed to replace them with several integer instructions. On today's hardware, and in an interpreted language, that is almost certainly not the case anymore. Are you targeting a Raspberry Pi? (+1 for nostalgia, though) – Niki Estner Oct 7 '15 at 8:43
• @nikie I used "Bresenham" more as a brand-name for the method of rendering lines, I'm not interested specifically in restricted to integers operations only algorithm. I'm looking for any efficient analytical solution leading to the same result. My motivation is both educational and finding a way to implement subpixel "Bresenham". – Alexey Popkov Oct 7 '15 at 8:59

If you only want the pixel coordinates of a line, that can be done much simpler:

Clear[pointsOnLine]
pointsOnLine[{p1_, p2_}, nPts_] :=
Array[Round[p1 + # (p2 - p1)] &, nPts, {0., 1.}]
pointsOnLine[l_: {p1_, p2_}] :=
pointsOnLine[l, Round[Max[Abs[p1 - p2]]] + 1]

The result is the same as your algorithm: but it's much faster:

p1 = {2, 3}; p2 = {2001, 1300};
RepeatedTiming[bresenhamSolve[p1, p2];]
RepeatedTiming[bresenham[p1, p2];]
RepeatedTiming[pointsOnLine[{p1, p2}];]

{0.0777, Null}

{0.0120, Null}

{0.000919, Null}

ADD: If you want exactly the same results as the original bresenham function, you have to force Round to round .5 up, instead of to the nearest even integer. Easiest way to do this is to just add a small epsilon to every value:

epsilon = 10^-10;
Clear[pointsOnLine, p1, p2]
pointsOnLine[{p1_, p2_}, nPts_] :=
Array[Round[p1 + # (p2 - p1) - epsilon] &, nPts, {0., 1.}]
pointsOnLine[{p1_, p2_}] :=
pointsOnLine[{p1, p2}, Round[Max[Abs[p1 - p2]]] + 1]
• That is cool but not identical to bresenham: try p1 = {0, 0}; p2 = {20, 10};, p1 = {-1, -2}; p2 = {17, -38};. Also check the code under words "Checking symmetry:" in the question - the pattern is entirely different. – Alexey Popkov Oct 8 '15 at 9:10
• That's because Round rounds .5 towards the nearest even integer. You can use e.g. rnd = Ceiling[# - .5] &;. That gives a symmetric result, but it will still sometimes round different than bresenham, since the calculation is done in machine precision. – Niki Estner Oct 8 '15 at 9:54
• I'm not very familiar with the image processing functionality in Mathematica. Is there a way to do something like ClusteringComponents, but locally, i.e. compensate for the image being darker/brighter in some places? I'm trying to segment a brain into gray matter / white matter. – Szabolcs Oct 8 '15 at 12:08

Original Bresenham

Starting with the code by nikie I was able to implement original Bresenham using exact arithmetic:

Clear[pointsOnLine]
pointsOnLine[{p1_, p2_}] :=
With[{d = p2 - p1},
p1 + IntegerPart[#] + Round[FractionalPart[#]] & /@
Array[# d &, Max[Abs[d]] + 1, {0, 1}]]

Unfortunately I haven't found a way to reproduce original Bresenham with floating point arithmetic.

Proof that this implementation is is identical to halirutan's implementation:

lines = DeleteCases[Partition[#, 2] & /@ Tuples[Range[-30, 30, 6], {4}], {p_, p_}];
Length[lines]
pointsOnLine /@ lines == bresenham @@@ lines
14520

True

Timings:

p1 = {2, 3}; p2 = {2001, 1300};
RepeatedTiming[bresenhamSolve[p1, p2];]
RepeatedTiming[bresenham[p1, p2];]
RepeatedTiming[pointsOnLine[{p1, p2}];]
{0.15, Null}

{0.027, Null}

{0.025, Null}

An alternative to Bresenham

Here is simpler and more efficient implementation which produces a pattern with the same symmetry as Bresenham's using another rounding scheme:

pointsOnLine[{p1_, p2_}] :=
With[{d = p2 - p1}, p1 + # & /@ Round[Array[# d &, Max[Abs[d]] + 1, {0, 1}]]]

The pattern:

n = 20; center = {0, 0};
perimeterOfSquare = {x, y} /.
Solve[{x, y} \[Element]
RegionBoundary[Rectangle[{-n, -n} + center, {n, n} + center]], {x, y}, Integers];
ArrayPlot[SparseArray[
Rule @@@ Tally[# - center & /@
Flatten[pointsOnLine[{center, #}] & /@ perimeterOfSquare + n + 1, 1]], {2 n + 1,
2 n + 1}], Mesh -> True, PlotRange -> {All, All, {1, 9}}, ClippingStyle -> Red,
PixelConstrained -> True] Here is Manipulate which allows to explore the differences between bresenham and different implementation of pointsOnLine:

Manipulate[Row[{Graphics[{{EdgeForm[], FaceForm[GrayLevel[.6]],
Rectangle /@ (pointsOnLine@Round[pts] - .5)}, {EdgeForm[{Thick, Darker@Red}],
FaceForm[], Rectangle /@ (bresenham @@ Round[pts] - .5)}, {Black, Thick,
Arrow@pts}}, GridLines -> ({Range[-20, 20], Range[-20, 20]} - .5), Frame -> True,
PlotRange -> {{-20, 20}, {-20, 20}}, ImageSize -> 500],
Column[Round@pts]}], {{pts, {{-11, -13}, {15, 0}}}, Locator}] 