- the pixel immediately to the right
- the pixel which is one pixel up and one pixel right

Let's assume we're moving from 00 to 10. When we increase the X coordinate, the error will increase by DY/DX. If the error becomes greater than 0.5, then the pixel above (11) should be selected. 0.5 is the distance between the center of a pixel and its top border. So, our step is:

if (error+dy/dx > 0.5) y = y+1 endifIf we make a right move, the error update is trivial:

error= error + dy/dxOtherwise, the error is calculated from the center of the pixel above. We observe that

E1+DY/DX+E0= 1 E1= 1 - E0 - DY/DXHowever, the sign of the real error is negative, so we take that instead.

E1'= -1 + E0 + DY/DXWhich means that:

error = error + dy/dx - 1In order to get rid of dy/dx terms, we can easily multiply all values with 2dx. The condition and the updates become:

if error + 2dy > dx error = error + 2dy error = error + 2dy - 2dxSumming up, the algorithm becomes:

error= 0 y= y0 for x = x0 .. x1 pixel(x,y) if error + 2dy > dx y= y+ 1 error = error - 2dx end error= error + 2dy endSo, let's put this in the Murphy format for later use:

thin_octant1(x0,y0,dx,dy): error= 0 y= y0 x= x0 threshold = dx - 2dy E_diag= -2dx E_square= 2dy length = dx for p= 1 .. length pixel(x,y) if error > threshold y= y + 1 error = error + E_diag end error = error + E_square x= x + 1 end

The error is again initially zero. The test for diagonal/square move is:

if error - dy/dx < -0.5The error update for square move is:

error= error - dy/dxand for the diagonal move:

E1 + dy/dx - E0 = 1 here, we subtracted E0 because its sign is negative.. E1 = 1 - dy/dx + E0 error = error + 1 - dy/dxSumming up,

if error - dy/dx < -0.5 error= error - dy/dx error = error + 1 - dy/dxSomething strange happens if we multiply all values with -2dx (instead of 2dx in the thin line algorithm).

if error + 2dy > dx error= error + 2dy error= error - 2dx + 2dyThese are exactly the same things as the previous set of equations. So, our algorithm for the left-perpendicular is almost the same, only the coordinate updates are different.

thin_octant1(x0,y0,dx,dy): pleft_octant1(x0,y0,dx,dy): error= 0 error= 0 y= y0 y= y0 x= x0 x= x0 threshold = dx - 2dy threshold = dx - 2dy E_diag= -2dx E_diag= -2dx E_square= 2dy E_square= 2dylength = dx length = ..desired-width..for p= 1 .. length for p= 1 .. length pixel(x,y) pixel(x,y) if error > threshold if error > thresholdy= y + 1 x= x - 1error = error + E_diag error = error + E_diag end end error = error + E_square error = error + E_squarex= x + 1 y= y + 1end end

Here, the green and yellow squares represent pixels painted on alternate runs of the perpendicular loop. As you can see, no line passes thru the red square. This will happen every time the thick_octant1 increases y and thus makes a diagonal move. There are more red squares to the left and top, we could see those if the image was big enough.

Why does this happen? It's because the perpendicular lines start from the initial state (error=0) and don't follow each other. For example, consider the first green and first yellow sequence. They have the same phase at each y coordinate: they are essentially the same sequence since their initial y value is also identical. Whenever one makes a diagonal move, the other does it too. This way, there is no gap left between them.

The same thing doesn't happen with the first yellow and second green sequence. The yellow sequence does a diagonal move at y=1 but the green sequence does a square move since it's the very first move for that sequence. If I could somehow execute the same sequence at the same y coordinate, the problem would be solved.

The solution to this is to start all perpendicular lines not at (x,y) but at (x,y0) and draw only the part that falls to the left side of the blue line. I'm not going to draw it here. Just imagine that you layout copies of the first green line along the x axis. Obviously we wouldn't miss any squares. However, this isn't efficient or elegant. I should figure out how to calculate the perpedicular's error at (x,y) without going from (x,y0) to (x,y).

We could calculate the 'phase' or 'accumulated error' of the perpendicular within the thick_octant1 function using the E_diag, E_square and threshold values. However, we update this error variable only when we make a diagonal move in the thick_octant1 function. When we make a square move, the y coordinate doesn't change and the next perpendicular should start with the same error we have now. When we make a diagonal move, y coordinate has changed and a new error should be calculated to decide whether we will go with the current x coordinate or one less in the perpendicular. So, our thick_octant1 should look like this:

thick_octant1(x0,y0,dx,dy): p_error= 0 # the perpendicular error or 'phase' error= 0 y= y0 x= x0 threshold = dx - 2dy E_diag= -2dx E_square= 2dy length = dx for p= 1 .. length pleft_octant1(x,y, dx, dy, p_error) if error > threshold y= y + 1 error = error + E_diag if p_error > threshold p_error= p_error + E_diag end p_error= p_error + E_square end error = error + E_square x= x + 1 endThe new pleft_octant1 is almost the same as the old one. It just doesn't initialize error to 0, it gets the initial value for it from the thick_octant1 function.

pleft_octant1(x0,y0,dx,dy,error): y= y0 x= x0 threshold = dx - 2dy E_diag= -2dx E_square= 2dy length = ..desired-width.. for p= 1 .. length pixel(x,y) if error > threshold x= x - 1 error = error + E_diag end error = error + E_square y= y + 1 endThat's pretty much it, but one more issue remains. Before I go into that, let's look at how this thing runs.

We managed to make all perpendiculars 'in-phase', i.e. all perpendiculars do the same thing (diagonal or square move) at the same y. However, some perpendiculars seem to be missing completely. This situation is what Mr. Murphy calls a 'double square move'.

In our particular example, the very first move is a diagonal one. We draw the first perpendicular (yellow) normally. Then we move on to the next start point using a diagonal move. Since our first move for the first perpendicular was also a diagonal move, we failed to cover the square right above the red square. In order to overcome this, we add another call to pleft_octant1 in the thick_octant1 function. Just where we detect this situation:

thick_octant1(x0,y0,dx,dy): p_error= 0 # the perpendicular error or 'phase' error= 0 y= y0 x= x0 threshold = dx - 2dy E_diag= -2dx E_square= 2dy length = dx for p= 1 .. length pleft_octant1(x,y, dx, dy, p_error) if error > threshold y= y + 1 error = error + E_diag if p_error > thresholdHere, we used the phase of thepleft_octant1(x,y, dx, dy, p_error+E_diag+E_square)p_error= p_error + E_diag end p_error= p_error + E_square end error = error + E_square x= x + 1 end

- The error values will be subtracted instead of added to the current error.
- The error will be compared against -0.5 for diagonal decision.
- x and y coordinates will be updated in the reverse direction.

pleft_octant1(x0,y0,dx,dy,error): pright_octant1(x0,y0,dx,dy,error): y= y0 y= y0 x= x0 x= x0 threshold = dx - 2dy threshold = dx - 2dy E_diag= -2dx E_diag= -2dx E_square= 2dy E_square= 2dy length = ..desired-width.. length = ..desired-width.. for p= 1 .. length for p= 1 .. length pixel(x,y) pixel(x,y) if error > threshold if error < -threshold x= x - 1 x= x + 1 error = error + E_diag error = error - E_diag end end error = error + E_square error = error - E_square y= y + 1 y= y - 1 end endIf I do the multiplication trick using the value -1, I end up with the following definition for pright_octant1:

pright_octant1(x0,y0,dx,dy,error): y= y0 x= x0 threshold = dx - 2dy E_diag= -2dx E_square= 2dy length = ..desired-width..This is almost the same thing as pleft_octant1 except for the directions of the coordinate updates. The underlined statement comes from the multiplication by -1. I feel like I could show this also geometrically if I represented the error and the actual line as vectors.error= -errorfor p= 1 .. length pixel(x,y) if error > threshold x= x + 1 error = error + E_diag end error = error + E_square y= y - 1 end

In order to actually draw the right perpendiculars, we can simply insert the call into the places where the left perpendicular is used:

thick_octant1(x0,y0,dx,dy): p_error= 0 error= 0 y= y0 x= x0 threshold = dx - 2dy E_diag= -2dx E_square= 2dy length = dx for p= 1 .. length pleft_octant1(x,y, dx, dy, p_error)The result is quite nice, but the width errors are even more pronounced.pright_octant1(x,y, dx, dy, p_error)if error > threshold y= y + 1 error = error + E_diag if p_error > threshold pleft_octant1(x,y, dx, dy, p_error+E_diag+E_square)pright_octant1(x,y, dx, dy, p_error+E_diag+E_square)p_error= p_error + E_diag end p_error= p_error + E_square end error = error + E_square x= x + 1 end

Also, the right side looks wider. This is because I first draw the left side, then the right side in both cases starting from the same point. So, this middle point looks like it belongs to the right side of the thick line. This also causes the middle point for each perpendicular pair to be painted twice. I will fix this in the real implementation and skip the first pixel for the left or right side.

The angle m is the angle of our base line. So

D= sqrt(dx^2+dy^2) tan(m)= dy/dx = U/1 cos(m)= K/1 = T/(U+1) = dx/D K= dx/D T= dx/D * (dy/dx+1) = dx/D * (dy+dx)/dx = (dy+dx)/DWhen we scale these equations with 2D, we get two things:

- We get rid of the floating point numbers.
- The scale matches the scale of the error variables.

Here, we're interested in the length of the red segment within the pixel. e is the error of the (lower) blue base line. The green lines are parallel to the x axis. Therefore, they form triangles with the base line angle m.

tan(m)= p/0.5 = dy/dx p= 0.5dy/dx r= 0.5 - e p+r= 0.5 - e + 0.5dy/dx cos(m)= U/(p+r) = dx/D U= dx/D*(0.5 -e + 0.5dy/dx) scaled with 2D 2dx * 0.5 - 2dx*e + 2dx*0.5dy/dx = dx - 2dx*e + dyWe already know what 2dx*e is: it's the error variable of the base line. For the thickness accumulators in perpendiculars, the initial values and updates are:

LEFT RIGHT INIT w= dx+dy-error_A w= dx+dy+error_A SQUR w= w + 2dx w= w+2dx DIAG w= w + 2dx + 2dy w= w+2dx+2dyInterestingly, Murphy doesn't include the term dx+dy for the initial thickness value. It seems to be getting it wrong for this reason. When I draw a 45 degree line with width 10, my algorithm draws one with width 9.9 (7 pixels across). Murphy's algorithm draws the line with width 11.3 (8 pixels across). Here is the final result of the algorithm.

Below is the final version for octant 1. Main improvements over previous versions are:

- Left and right perpendiculars are collected into one function.
- Proper thickness calculation is done.

perp_octant1(x0,y0,dx,dy,einit,width,winit): threshold = dx - 2*dy E_diag= -2*dx E_square= 2*dy wthr= 2*width*sqrt(dx*dx+dy*dy) x,y= x0,y0 error= einit tk= dx+dy-winit while tk<=wthr: pixel(x,y) if error > threshold: x= x - 1 error = error + E_diag tk= tk + 2*dy end error = error + E_square y= y + 1 tk= tk + 2*dx end x,y= x0,y0 error= -einit tk= dx+dy+winit while tk<=wthr: pixel(x,y) if error > threshold: x= x + 1 error = error + E_diag tk= tk + 2*dy end error = error + E_square y= y - 1 tk= tk + 2*dx end thick_octant1(x0,y0,x1,y1,width): dx= x1 - x0 dy= y1 - y0 p_error= 0 error= 0 y= y0 x= x0 threshold = dx - 2*dy E_diag= -2*dx E_square= 2*dy length = dx+1 for p= 1 .. length: perp_octant1(x,y, dx, dy, p_error,width,error) if error > threshold: y= y + 1 error = error + E_diag if p_error > threshold: perp_octant1(x,y, dx, dy, p_error+E_diag+E_square,width,error) p_error= p_error + E_diag end p_error= p_error + E_square end error = error + E_square x= x + 1 endWithout realizing why I did it, I put the second call to perpendicular

error= 0 y= y0 for x = x0 .. x1 pixel(x,y) if error + 2dy > dx y= y+ 1 error = error - 2dx end error= error + 2dy endIf we always use non-negative numbers for dx and dy, obviously the error computation doesn't depend on the direction of the line.

Given non-negative dx and dy, the error code computes a particular sequence of square and diagonal moves. You can use this sequence to move in any of the four directions. For an x-based line our code would look like:

error= 0 x,y= x0,y0 for p= 1 .. dx+1 pixel(x,y) if error + 2dy > dx y= y+ ystep error = error - 2dx end error= error + 2dy x= x + xstep endBy choosing xstep and ystep from {-1,1}, we can make this code paint any x-based line. The same argument applies to all functions we've written so far. Just change the increments and you cover the other 3 octants.

For handling y-based lines where abs(dy)>abs(dx), we need a different function because now we're looping over the y coordinates and occasionally updating x coordinates. This could still be stuffed into the same function using even more variables and adding some occasionally-nop code but we will have to sacrifice code readability.

I realized that the line p0->p1 isn't necessarily the same line as p1->p0. The error is defined to be zero for p0. When the code ends, it isn't necessarily zero at p1. The non-zero situation happens whenever dy doesn't divide dx (for octant1). However, it would be zero for p1 if we had switched the points around.

To paraphrase it: the error calculation computes a particular sequence of square and diagonal moves from p0 to p1. This sequence doesn't have to be symmetric. Therefore, if you applied the sequence from p1 to p0, you'd end up with a different set of pixels painted.

With a minor change to pright_octant1 and its correspondent in perp_octant1, the algorithm can be made to traverse each pixel only once. The only pixels drawn twice are the ones on the base line. This could be eliminated easily. This would allow the algorithm to be used for situations where the objective is not to directly paint but execute some operation on the pixels, such as negating the background, doing an intersection etc.

Since we have a measure of the current thickness in the perpendicular function, we could draw anti-aliased lines using this algorithm. We could vary the brightness as we move away from the center of the thick line.

Let's say that we're doing the very first example, basic Bresenham
line in the first octant. When going forward, the **error≥threshold**
condition selects the upper pixel when error=threshold. When going in
reverse from the upper pixel, if error=threshold, we must stay in the
upper y coordinate. We should switch to the lower y coordinate when error
exceeds threshold. This is necessary because otherwise the left and
right sides of a perpendicular get out of phase and artifacts occur.

Another failed assumption was that adapting the code for other octants would be trivial. Not really so. The xstep and ystep values for the perpendiculars have to be chosen carefully such that what I mention as 'left' in this article really follows the normal we've discussed. i.e. the normal coming in from bottom-right to top-left for a line in the first octant.

Actually, trying to go left in the first perpendicular doesn't really work. What we have really discussed in the introduction to perpendicular part is the normal starting from the y-side of the origin and continuing towards to y-side of the target (for an x-based line). This isn't always left. However, I made arrangements in the code such that left and right of a vector are still distinguishable and can be given different thicknesses.

Here is an implementation of the algorithm in this article. There is only one entry function, which takes two functions as arguments. These functions should return the width for the left and right sides. The arguments to these functions are:

- A user pointer.
- The position along the line.
- The total length of the line.

Here are the figures I used,
in case I need to modify them. Below you can see a test run, showing
the different cases I mentioned in the "Variants" section.

The code could be improved in several ways. First of all, there is
no simple fixed-width line drawing function. A wrapper should be
implemented for that.

Secondly, the outer and inner loops aren't in a state suitable for function inlining. The two calls to 'perpendicular' could be reduced to one by playing with the outer loop control a little bit. This would increase the likelyhood of inlining the inner loop since it's a static function called only from one function.

The code I provide here is part of my drawing package I'm working on. You can't compile it out of the box, but adapting it is simple. Simply modify the pbuffer_t* argument to your image type and replace the six calls to 'pixel' accordingly.

When doing this, it's very important to initialize K properly. Let's make an example, the basic Bresenham line stuff:

error= 0 .. if error+dy/dx > 0.5: y = y+1 end .. error= error + dy/dx .. error = error + dy/dx - 1Setting EB= 2dx*error and multiplying everything above with 2dx:

EB= 0 .. if EB+2dy > dx y = y+1 end .. EB= EB + 2dy .. EB = EB + 2dy - 2dxThen I replace EB with 'error' again.

- Program to draw a line using Bresenham's Line Algorithm (BLA)
- The Bresenham Line-Drawing Algorithm
- Bresenham-based supercover line algorithm - paints all the points the line touches.
- Bresenham - Several C functions implementing the Bresenham algorithm for different objects suchs as lines, circles etc.
- Algorithms for drawing thick lines and curves on raster devices
- Murphy's Modified Bresenham Line Drawing Algorithm - The page by Alan Murphy