The first thing to make note of is that the algorithm listed in the question is the integer algorithm. So no floating point math and the plotting increments by +1 or -1 in the x and y directions.
The function takes for input the line starting position (x0,y0) and ending position (x1,y1). The starting position variables(x0,y0) are reused when drawing the line.
Using inputs the algorithm computes the x and y increment directions. These values get held in variables sx and sy and are always either +1 or -1. They can be negative or positive depending on the line. For example we should get consistent results whether or not (x0,y0) and (x1,y1) are switched.
sx = x0 < x1 ? 1 : -1; // step x left or right
sy = y0 < y1 ? 1 : -1; // step y left or right
To actually draw the line we need to determine whether to step in the x direction, the y direction or both x and y (ie diagonally). (that is why the algorithm uses two separate if statements instead of an if then else.) To step in the x direction it is simply:
x0 = x0 + sx; // move x0 in the x direction setup in sx
y0 = y0 + sy; // move y0 in the y direction setup in sy
Our end condition happens when x0==x1 or y0==y1 or (x0,y0) == (x1,y1). This is handled in the algorithm by the "if(..) then break;" statements.
if x0 == x1 && y0 == y1 break; // break the while true loop
if x0 == x1 break; // only x meets the break condition
if y0 == y1 break // only y meets the break condition
Finally we need to be able to decide when to step in the x direction or the y direction or both. That is where the variables e2, dx, dy, and error come in.
Notice that dx is always positive ie dx = abs(x1-x0) and dy is always negative ie dy = -abs(y1-y0) ... and error starts out as the sum of a positive and negative number ie error = dx+dy. The algorithm adds dx(which is postive) to error when the line is plotted in the x direction, and it adds dy(which is negative) to error when it is plotted in the y direction.
dx = abs(x1-x0); // distance x needs to move in increments of 1
dy = -abs(y1-y0); // negative distance y needs to move in increments of 1
error = dx + dy; // Start state for difference between dx and dy
e2 = error * 2; // look ahead
If for example the line is diagonal such as the line (0,0):(5,5). The value of error will be zero, and so will e2. Making both if statements true for each step. And since dx is added and dy is subtracted both if statements will remain true until x0,y0==x1,y1!
dx and dy can be thought of as the numerator and denominator of the slope of a line. When dx is 3 units bigger then dy, it will increment 3 times for every 1 dy.
So how do we figure out the next step? Simple look ahead one step which is e2=2*error When either x0 or y0 is greater then e2 then they have gone too far, ie they are in an "error state".
One way to look at this is: As the algorithm plots in the x direction error will increase in the positive direction eventually making e2 greater then dx...which will stop it from plotting in the x direction ie if(e2<=dx) will become false. Also as the algorithm plots in the y direction error will increment in the negative direction eventually making e2 less then dy, ie if(e2>=dy) will become false.
e2 isn't being pulled out of the air, it is the calculation of the next step.
So that's a pretty long description...the way I originally came to understand it was to code it up in C and then step through various lines in the debugger looking at the values until that light bulb moment happened for me.