CS-460/560, Week 3-C
Spring, 1998
R. Eckert



Task--Given pixel coordinates of endpoints P1 (x1,y1) & P2 (x2,y2)
      Determine which pixels need to be turned on.

   Straight as possible
   Constant density (no gaps or bunching)
   Pixel density as independent as possible of orientation
   Must be fast

Line Equations--

First the differential equation--

   dy/dx = m    (m=constant--the slope)

Integrate (indefinite) --> y = m*x + constant
The constant (b) is called the y intercept (value of y when x=0).

   y = m*x + b   (slope-intercept form)

Alternatively we can integrate between the endpoints (definite)-->

   (y2-y1)/(x2-x1) = m    (gives an operational definition of the slope)

Integrate between endpoint (x1,y1) and arbitrary point (x,y) -->

   y - y1 = m*(x-x1)

   y = m*(x-x1) + y1     (point-slope form)

Brute Force Line-Drawing Algorithm (step in x direction, x2>x1)--

   Compute m = (y2-y1)/(x2-x1)
   num-pts = x2-x1+1
   x = x1
   Repeat num-pts times
      y = m*(x-x1) + y1
      x = x+1

Problem--if |y2-y1| > |x2-x1| there will be gaps

In that case should step in y direction. Following assumes y2>y1.

   Compute inv_m = (x2-x1)/(y2-y1)
   num-pts = y2-y1+1
   y = y1
   Repeat num-pts times
      x = inv_m*(y-y1) + x1
      y = y+1

Cases of vertical lines (x2=x1) and horizontal lines (y2=y1) should be
handled separately.

Brute Force Method is TOO SLOW--A floating point multiply and add, as
well as round() operations on each iteration.

INCREMENTAL METHODS--The Digital Differential Analyzer (DDA)

Basic idea--get the new point (x,y) from the previous point plotted.

dy/dx = m   ==>  delta_y/delta_x = m   (for straight lines this is exact)
                 Or delta_y = m*delta_x

where delta_y = ynew - yold, delta_x = xnew-xold

So xnew = xold + delta_x
   ynew = yold + delta_y = yold + m*delta_x

Choose, for example, delta_x = 1 (stepping in x direction)

Then ynew = yold + m

DDA Algorithm (stepping in x direction, x2>x1)--

   Compute m = (y2-y1)/(x2-x1)
   num-pts = x2-x1+1
   x = x1
   y = y1
   Repeat num-pts times
      x = x+1
      y = y+m

As for the Brute force method, if |m|>1, we can step in y--

DDA Algorithm (stepping in y direction, y2>y1):

   Compute inv_m = (x2-x1)/(y2-y1)
   num-pts = y2-y1+1
   x = x1
   y = y1
   Repeat num-pts times
      y = y+1
      x = x+inv_m

The floating point multiply in the loop is gone!! (But there is still
a floating point add and a round()--WE CAN DO BETTER!)

For best performance, we need an algorithm that only uses only integer 
adds or subtracts inside the loop.

Bresenham's Line-drawing Algorithm--
   Used in some form in most graphics packages
   Often implemented in hardware
   Also incremental (Get new pixel from old pixel)
   Uses only integer operations

All lines can be placed in one of four categories:
   A. Steep + slope (m>1)
   B. Gradual + slope (0<m<=1)
   C. Steep - slope (m<-1)
   D. Gradual - slope (0>=m>=-1)

For each case, if we have just plotted the point (x,y), the next point to
be plotted can be one of only two possibilities:

   A. (x,y+1) or (x+1,y+1)
   B. (x+1,y) or (x+1,y+1)
   C. (x,y-1) or (x+1,y-1)
   D. (x+1,y) or (x+1,y-1)

We want to come up with some "Predictor" function that uses only simple
integer math to determine which of the two choices is correct.

In the following we will look at Case-A (Steep positive slope). We also
assume that P1 is to the left of P2 (x1<x2). If that is not true, the
points can be switched.

Since delta_y is greater than delta_x for this case, we are stepping in
the y direction.

Assume that the last pixel plotted is at (r,s). The next pixel then has
to be at either (r,s+1) (the "left" pixel) or (r+1,s+1) (the "right pixel).

The actual mathematical point would lie on the line at (rp, s+1).

We define two horizontal distances--dl, the distance between the
mathematical point and the "left" pixel, and dr, the distance between
the mathematical point and the "right" pixel.

   dl = rp-r
   dr = r+1-rp

If dl<dr, the "left" pixel is closer to the actual point than the "right"
pixel, so we should plot it. Instead of looking at that inequality, we
look at dl-dr. If dl-dr<0, we need to choose the "left" pixel. Our
criterion for choosing it is:

   dl-dr = rp-r - (r+1-rp) < 0,   or:    

   2*rp - 2*r -1 < 0

The equation for the line is y=m*x+b, or y = (del_y/del_x)*x + b, where
we have defined del_y=y2-y1 and del_X=x2-x1. (Recall that (x1,y1) and
(x2,y2) are the endpoints, which are given.)

Since the mathematical point lies on the line, it must satisfy the equation:

   s+1 = (del_y/del_x)*rp + b,   so:

   rp = (s+1-b)*del_x/del_y,     so:

   dl-dr = (2*s*del_x + 2*del_x - 2*b*del_x)/del_y - 2*r - 1

If this is negative, we need to choose the "left" pixel.

But we have gained nothing--this "predictor" is too complex.

Let's multiply through by del_y (always positive for Case-A lines), and call
the result the "predictor", P

   P = del_y*(dl-dr) = 2*del_x*(s+1-b) - 2*r*del_y - del_y 

At least the divide is gone--but it's still too complex.

Bresenham's contribution is at this point. He reasoned that, just as we're
trying to get the new point from the old point, perhaps we can get a new
value of the predictor from its old value in a simple way. In other words,
we need a recurrence relation for the predictor. In the following, sn means
the new value of s; so the old value; similarly rn and ro for r, and Pn
and Po for P. If we can find the recurrence relation, we will have:

   del_P = Pn - Po, or in other words:

   Pn = Po + del_P, which will give us the new predictor from the old one.

There are two cases:  the left case in which rn=ro and sn=so+1
                      the right case in which rn=ro+1 and sn=so+1

For both cases:

   Po = 2*del_x*(so+1-b) - 2*ro*del_y - del_y

Left case:

   Pn = 2*del_x*((so+1)+1-b) - 2*ro*del_y - del_y

   Subtracting Po from Pn gives us del_P, and a lot of stuff drops out!
   The result is: 

   del_P = 2*del_x

Right case:

   Pn = 2*del_x*((so+1)+1-b) - 2*((ro+1)*del_y - del_y
   Again subtracting gives:

   del_P = 2*(del_x - del_y)

Both these results are very simple--and ALL INTEGERS!

So we look at the current value of the predictor

   If P < 0    // left case
      P = P + 2*del_x
      x = x
      y = y + 1

   If P>0      // right case
      P = P + 2*(del_x-del_y)
      x = x + 1
      y = y + 1

But for this to work, we need to have an initial value P0 of the predictor.
Substituting the left-hand endpoint (x1,y1) into the predictor definition:

   P0 = 2*del_x*(y1+1-b) -2*x1*del_y - del_y

And using the fact that y1 = (del_y/del_x)*x1 + b, the result is:

   P0 = 2*del_x - del_y

So our Bresenham algorithm for plotting Case-A lines (|y2-y1| > |x2-x1| and
assuming that x1<x2) is:

   del_x = x2-x1
   del_y = y2-y1
   P = 2*del_x - del_y
   cleft = 2*delx
   cright = 2*del_x - 2*del_y
   num_pts = |del_y| + 1
   x = x1
   y = y1
   Repeat num_pts times
      y = y + 1
      If (P < 0)
         P = P + cleft
         P = P + cright
         x = x + 1

This can be generalized to handle Case-C (steep - slope) lines by computing

   sdy = sign(del_y) = 1 if y2>y1
                     = -1 if not

Then, in the algorithm's definition of P and cright, replace del_y with


In addition, the statement y = y + 1 would be replaced with:

   y = y + sdy

With these minor changes the algorithm handles both Case-A and Case-C lines.


See the Hearn & Baker Text Book, Section 3-1 (pages 88-95) for more information
on the Bresenham Line-drawing algorithm, specifically Case-B lines.