graphics - Bresenham's line algorithm

Bresenham's line algorithm is producing discretized line for given two points for purpose of plotting for example.

Like that:

I have to stress that I'm interested in positions, not a plot.

Wikipedia link I've provided includes an algorithm of course. I've just rewritten it thoughtlessly, I don't have time now or special need to work on neat implementation.

But if someone want to improve it (compile e.g.), got it already or know something more, I think this thread may be usefull for future visitors.

Very nice implementation can be found on rossetacode :P, according to that this algorithm should be built in so maybe someone knows how to get it.

Anyway, here's that code:

bresenham[{x1_, y1_}, {x2_, y2_}] := 
 Module[{dx, xi, dy, yi, ai, bi, x, y, d},
     If[x1 < x2, {xi, dx} = {1, x2 - x1};, {xi, dx} = {-1, x1 - x2};];

     If[y1 < y2, {yi, dy} = {1, y2 - y1};, {yi, dy} = {-1, y1 - y2};];
     x = x1; y = y1;
     Sow[{x, y}];
     If[dx > dy,
      (ai = 2 (dy - dx); bi = 2 dy; d = bi - dx;
       While[
        If[d >= 0,
         {x, y, d} += {xi, yi, ai},
         {x, d} += {xi, bi}];
        Sow[{x, y}];

        x != x2])
      ,
      (ai = 2 (dx - dy); bi = 2 dx; d = bi - dy;
       While[
        If[d >= 0,
         {x, y, d} += {xi, yi, ai},
         {y, d} += {yi, bi}];
        Sow[{x, y}];
        y != y2])
      ]] // Reap // Last // First

Answer

Original Bresenham

I guess I can come of with a somewhat shorter implementation without using Reap and Sow. If someone is interested, it follows almost exactly the pseudo-code here

bresenham[p0_, p1_] := Module[{dx, dy, sx, sy, err, newp},
  {dx, dy} = Abs[p1 - p0];
  {sx, sy} = Sign[p1 - p0];
  err = dx - dy;
  newp[{x_, y_}] := 
   With[{e2 = 2 err}, {If[e2 > -dy, err -= dy; x + sx, x], 
     If[e2 < dx, err += dx; y + sy, y]}];

  NestWhileList[newp, p0, # =!= p1 &, 1]
]

To test this I use the setup given by the comment of Kuba under this answer:

p1 = {17, 1}; p2 = {7, 25}; 
Graphics[{EdgeForm[{Thick, RGBColor[203/255, 5/17, 22/255]}], 
  FaceForm[RGBColor[131/255, 148/255, 10/17]], 
  Rectangle /@ (bresenham[p1, p2] - .5), {RGBColor[0, 43/255, 18/85], 
   Thick, Line[{p1, p2}]}}, 
 GridLines -> {Range[150], Range[150]} - .5]

Exercise implementation

What follows was only a fun project I did with my wife. Actually, this is not the original Bresenham algorithm. The task for this weekend-fun was to re-invent the algorithm (the iterative steps and the required correction steps) on the blackboard.

For simplicity this algorithm only makes correction steps in one direction (meaning the points stay always on one half-plane of the line) and therefore, the final points are not as close to the original line as with the real Bresenham algorithm.

Anyway, this is my Mathematica implementation of what my wife had to do in Python:

bresenham[p1 : {x1_, y1_}, p2 : {x2_, y2_}] := 
 Module[{dx, dy, dir, corr, test, side},
  {dx, dy} = p2 - p1;
  dir = If[Abs[dx] > Abs[dy], {Sign[dx], 0}, {0, Sign[dy]}];

  test[{x_, y_}] := dy*x - dx*y + dx*y1 - dy*x1;
  side = Sign[test[p1 + dir]];
  corr = side*{-1, 1}*Reverse[dir];
  NestWhileList[
   Block[{new = # + dir}, If[Sign[test[new]] == side, new += corr];
     new] &, p1, #1 =!= p2 &, 1, 500]]

Here a small dynamic test whether the calculated points do indeed look like a line:

DynamicModule[{p = {{0, 0}, {50, 40}}},
 LocatorPane[Dynamic[p],

  Dynamic@
   Graphics[{Line[bresenham @@ Round[p]], Red, PointSize[Large], 
     Dynamic[Point[p]]}, PlotRange -> {{-200, 200}, {-200, 200}}, 
    ImageSize -> 400, Frame -> True, FrameTicks -> False, 
    GridLines -> True],
  Appearance -> None
  ]
 ]