[2.0.x] Fix Bresenham rounding errors

[ejtagle]

Fix Bresenham algorithm (rounding errors), add theory of operation and deduction of the employed algorithm.

On the Bresenham algorithm as implemented by Marlin:

(Taken from https://www.cs.helsinki.fi/group/goa/mallinnus/lines/bresenh.html)

The basic Bresenham algorithm:

Consider drawing a line on a raster grid where we restrict the allowable slopes of the line to the range 0 <= m <= 1

If we further restrict the line-drawing routine so that it always increments x as it plots, it becomes clear that, having plotted a point at (x,y), the routine has a severely limited range of options as to where it may put the next point on the line:

-It may plot the point (x+1,y), or:
-It may plot the point (x+1,y+1).

So, working in the first positive octant of the plane, line drawing becomes a matter of deciding between two possibilities at each step. We can draw a diagram of the situation which the plotting program finds itself in having plotted (x,y).

y+1 +--------------*
    |             /
    |            /
    |           /
    |          /
    |    y+e+m*--------+-
    |        /| ^    |
    |       / | |m   |
    |      /  | |    |
    |     /   | v    |
    | y+e*----|-----   |m+ε
    |   /|    | ^    |
    |  / |    | |ε   |
    | /  |    | |    |
    |/   |    | v    v
  y *----+----+----------+--
         x   x+1

In plotting (x,y) the line drawing routine will, in general, be making a compromise between what it would like to draw and what the resolution of the stepper motors actually allows it to draw. Usually the plotted point (x,y) will be in error, the actual, mathematical point on the line will not be addressable on the pixel grid. So we associate an error, ε , with each y ordinate, the real value of y should be y+ε . This error will range from -0.5 to just under +0.5.

In moving from x to x+1 we increase the value of the true (mathematical) y-ordinate by an amount equal to the slope of the line, m. We will choose to plot (x+1,y) if the difference between this new value and y is less than 0.5

y + ε + m < y + 0.5

Otherwise we will plot (x+1,y+1). It should be clear that by so doing we minimise the total error between the mathematical line segment and what actually gets drawn on the display.

The error resulting from this new point can now be written back into ε, this will allow us to repeat the whole process for the next point along the line, at x+2.

The new value of error can adopt one of two possible values, depending on what new point is plotted. If (x+1,y) is chosen, the new value of error is given by:

ε[new] = (y + ε + m) - y

Otherwise, it is:

ε[new] = (y + ε + m) - (y + 1)

This gives an algorithm for a DDA which avoids rounding operations, instead using the error variable ε to control plotting:

ε = 0, y = y[1]
for x = x1 to x2 do
  Plot point at (x,y)
  if (ε + m < 0.5)
    ε = ε + m
  else
    y = y + 1, ε = ε + m - 1
  endif
endfor

This still employs floating point values. Consider, however, what happens if we multiply across both sides of the plotting test by Δx and then by 2:

        ε + m < 0.5
    ε + Δy/Δx < 0.5
2.ε.Δx + 2.Δy < Δx

All quantities in this inequality are now integral.

Substitute ε' for ε.Δx . The test becomes:

2.(ε' + Δy) < Δx

This gives an integer-only test for deciding which point to plot.

The update rules for the error on each step may also be cast into ε' form:
Consider the floating-point versions of the update rules:

ε = ε + m
ε = ε + m - 1

Multiplying through by Δx yields:

ε.Δx = ε.Δx + Δy
ε.Δx = ε.Δx + Δy - Δx

Which is in ε' form:

ε' = ε' + Δy
ε' = ε' + Δy - Δx

Using this new "error" value, ε' with the new test and update equations gives Bresenham's integer-only line drawing algorithm:

ε' = 0, y = y[1]
for x = x1 to x2 do
  Plot point at (x,y)
  if (2.(ε' + Δy) < Δx)
    ε' = ε' + Δy
  else
    y = y + 1, ε' = ε' + Δy - Δx
  endif
endfor

It is an Integer-only algorithm - hence efficient (fast). And the Multiplication by 2 can be implemented by left-shift. 0 <= m <= 1

A possible implementation in C:

void bresenham(Screen &s,
          unsigned x1, unsigned y1,
          unsigned x2, unsigned y2,
          unsigned char colour )
{
  int dx  = x2 - x1,
      dy  = y2 - y1,
      y   = y1,
      eps = 0;

  for ( int x = x1; x <= x2; x++ )  {
    s.Plot(x,y,colour);
    eps += dy;
    if ( (eps << 1) >= dx )  {
      y++;  eps -= dx;
    }
  }
}

We can optimize it a bit, by noting that

(eps << 1) >= dx

Is equivalent to:

 eps >= (dx / 2)

Even if it would seem that right shifting dx would lose the least significant bit, we can handle it by proper rounding: eps is an integer variable, so only has integer values. dx is also an integer value, so after dividing by 2, either the value was an integer, or we missed by 0.5 less than true value.

1] If, after the division, the value was exactly an integer, then the following holds true:

(eps << 1) >= dx is equal to eps >= (dx>>1)       [ for dx even ]

2] If, after the division, the result was truncation of the last bit of dx, then

(eps << 1) >= dx is equal to eps >= (dx>>1)+0.5   [ for dx odd ]

But, as eps is an integer, we can rewrite the expression as:

(eps << 1) >= dx is equal to eps >   (dx>>1)       [ for dx odd ]
(eps << 1) >= dx is equal to eps >=  (dx>>1)+1     [ for dx odd ]

Now, the algorithm can be rewritten as:

void bresenham2(Screen &s,
          unsigned x1, unsigned y1,
          unsigned x2, unsigned y2,
          unsigned char colour )
{
  int dx  = x2 - x1,
      dy  = y2 - y1,
      y   = y1,
      eps = 0,
      limit = (dx>>1) + (dx&1);

  for ( int x = x1; x <= x2; x++ )  {
    s.Plot(x,y,colour);
    eps += dy;
    if ( eps >= limit )  {
      y++;  eps -= dx;
    }
  }
}

The last thing we can do to improve the algorithm, is instead of initializing eps to 0, we can initialize it to -((dx>>1) + (dx&1)) , so the algorithm becomes:

void bresenham3(Screen &s,
          unsigned x1, unsigned y1,
          unsigned x2, unsigned y2,
          unsigned char colour )
{
  int dx  = x2 - x1,
      dy  = y2 - y1,
      y   = y1,
      eps = -((dx>>1) + (dx&1));

  for ( int x = x1; x <= x2; x++ )  {
    s.Plot(x,y,colour);
    eps += dy;
    if ( eps >= 0 )  {
      y++;  eps -= dx;
    }
  }
}

And this is exactly the algorithm implemented by Marlin

[ejtagle]

@thinkyhead : Be careful, as the quoted article is missing the last piece of the puzzle, the one that explains the rounding needed to get it right. And that was the problem of the previous implementation!!

[ejtagle]

As they say, the devil is in the details. Grbl implementation is not exactly this one, but is also properly handling rounding

[thinkyhead]

I appreciate the effort to document the code, but quoting an entire article is serious overkill in the middle of our cpp files. We need to be able to work with this code on a daily basis, and concision that saves us the effort of scrolling is appreciated.

If we very much want to document the operational theory that goes into making Marlin work and we're concerned that the pages we're linking to might go away, then we should add a new section to marlinfw.org for that purpose and link to that instead.

[ejtagle]

@thinkyhead : That would be great, or at least, add a new folder to the solution (called docs/) and place all those comments there.
In fact, it has already happened: The paper that was used to write LA is not available anymore. Eventually, using webarchive,org, i was able to grab a copy...

[thinkyhead]

Was the previous code deduced to produce an error of one step sometimes?

I'd be curious to see an example where initializing to one less (with the rounding) and using >= 0 in place of > 0 produces a more correct number of steps than the previous code.

I could experiment with that, as I have a handy C file already set up for testing Bresnham variants.

[thinkyhead]

Hmm, I see that the simple Bresenham class I wrote to play around with also initializes the same way and does the same <= test. So, that's a good sign.

count[0] = -((inDividend + 1) >> 1);
. . .
const bool over = (count >= 0);

[thinkyhead]

Well, so far in my testing, rounding down and using > gives exactly the same pattern as rounding up and using >=. At least, so far I haven't hit on any magic values that result in a different pattern using this:

count[0] = -(inDividend >> 1);
. . .
const bool over = (count > 0);

Can you come up with any ratio where a step would be missed using this over the former?

[ejtagle]

It depends on how long the line is, and it has to be "nearly exact" multiple of the lenght. Then, at the end a step will be lost (it is just a rounding error)

[thinkyhead]

#include <iostream>
#include <string>
#include <math.h>
using namespace std;

class Bresenham {
  public:
    long divisor, dividend, count, ticker, invocations;

  public:

    Bresenham() { Bresenham(1, 1, 0); }

    Bresenham(const long &inDividend, const long &inDivisor, const long &inTicker=0) {
      cout << "Bresenham(" << inDividend << ", " << inDivisor << ")" << endl;
      invocations = 0;
      divisor = inDivisor;
      dividend = inDividend;
      count = -((inDividend + 1) >> 1);
      ticker = inTicker;
    }

    bool tick() {
      invocations++;
      count += dividend;
      const bool over = (count >= 0);
      if (over) {
        ticker++;
        count -= divisor;
      }
      report();
      return over;
    }

    long operator++() {
      tick();
      return ticker;
    }

    void report() {
      static long lastTicker = 0;
      rj(invocations, 3);
      cout << ": " << "count: ";
      rj(count[0], 3);
      cout << ", ticker: ";
      rj(ticker, 3);
      if (ticker != lastTicker) {
        cout << " *";
        lastTicker = ticker;
      }
      cout << endl;
    }
};

Bresenham bres(13, 14);

int main(int argc, char const *argv[]) {
  do { bres.tick(); } while (bres.invocations < bres.divisor);
  return 0;
}

[ejtagle]

I do a agree, it is a rare situation 👍 -- But math does not lie here 🥇

[thinkyhead]

It depends on how long the line is, and it has to be "nearly exact" multiple of the [length]. Then, at the end a step will be lost (it is just a rounding error)

I think we're already guaranteed to get the right number of ticks. Rounding closer to 0 (on the negative side of the number line) at the start, like we currently do, simply means we can use > rather than shifting our test one to the left with >=.

[ejtagle]

Ok, I´ll send a test program... 👍

[thinkyhead]

Surely you can see all you've done is shifted the count values one to the left,
necessitating >= rather than >. The shape of the line is no different.

[thinkyhead]

I realized that with positive counts there may be a difference. Here's an example where it differs.

Bresenham(8, 17)
  No rounding                     Rounding
  1: count: -13, ticker:   1 *	  count: -13, ticker:   1 *
  2: count:  -5, ticker:   1	  count:  -5, ticker:   1
  3: count: -14, ticker:   2 *	  count: -14, ticker:   2 *
  4: count:  -6, ticker:   2	  count:  -6, ticker:   2
  5: count: -15, ticker:   3 *	  count: -15, ticker:   3 *
  6: count:  -7, ticker:   3	  count:  -7, ticker:   3
  7: count: -16, ticker:   4 *	  count: -16, ticker:   4 *
  8: count:  -8, ticker:   4	  count:  -8, ticker:   4
  9: count:   0, ticker:   4	  count: -17, ticker:   5 *
 10: count:  -9, ticker:   5 *	  count:  -9, ticker:   5
 11: count:  -1, ticker:   5	  count:  -1, ticker:   5
 12: count: -10, ticker:   6 *	  count: -10, ticker:   6 *
 13: count:  -2, ticker:   6	  count:  -2, ticker:   6
 14: count: -11, ticker:   7 *	  count: -11, ticker:   7 *
 15: count:  -3, ticker:   7	  count:  -3, ticker:   7
 16: count: -12, ticker:   8 *	  count: -12, ticker:   8 *
 17: count:  -4, ticker:   8	  count:  -4, ticker:   8

A tick at iteration 9 instead of iteration 10.

[ejtagle]

The problem , as was stated in the report, is that > is not the same as >= due to the error rounding. At least, that causes 1isr pulse jitter (depending if dx is even or odd) ...

[ejtagle]

(and there is also something you don´t know: It is much cheaper to test for >= than to test for >. Reason is that the first one only needs to read the sign bit, and the other one need to read the whole value... ;)

[thinkyhead]

It is much cheaper to test for >= than to test for >

Ah, this is true!

[thinkyhead]

So, as long as we're counting cycles, which of these is faster?
-((dividend + 1) >> 1) or -((dividend >> 1) + (dividend & 1))

[ejtagle]

Is nearly the same. Maybe the first one, as there is one less operation 👍 -- Yes, the first one. One addition, then one shift.
The other one requires one addition, one shift and one AND operation

[thinkyhead]

Yep, the first one does less, compiles 14 bytes smaller, and uses 4 fewer registers. Excelsior!