x^2 + y^2 = r^2
After taking differentials of both sides, we find that
dy = -x/y dx
That is, if we take a step of size dx in the x-direction, we in principle want to take a step of size dy in the y-direction.
Next we start at the top of the circle, so that y=r and x=0. We start increasing x in step sizes of one. We only care about step sizes of one, since our basic unit is now a pixel. The y-coordinate is going to start piling up these dy's, and at some point the integer part of y will increase, and we get a new y-coordinate for the pixel. The idea, then, is to keep adding the dy's together, and once their sum is greater than one, we decrease y (remember that y starts at the top of the circle).
The sneaky way to do this is to treat y as an integer "constant". Then it is very easy to add the dy's together, since they have a common denominator equal to y. So really all we need to do is start adding x-coordinates together, and once their sum is larger than y, we decrease y and hang on to the remaining fractional part of dy. The algorithm then looks like:
y=r
x=0
a=r
loop: x=x+1
a=a-r
if a<=0 then a=a+y:y=y-1
plot (x,y)
if x
Now, Chris McBride pointed something out to me. As you may recall,
the algorithm breaks down for small r. Chris said that if a is initially
set to r/2 instead of r, the algorithm works perfectly. Why is that?
Recall that we add dy to itself until it is greater than one. Wouldn't
it make more sense to add dy to itself until it is greater than 0.5?
That would have the effect of rounding things up. Thus, starting at r/2
is like adding 0.5 to the fractional part of y -- it is the difference
between INT(y) and INT(y+0.5).
Thus, the above line
a=r
should be changed to
a=r/2
for a perfect circle every time. Thus, this corresponds to adding an LSR
to the machine code. Incidentally, this fix appeared in an earlier C=Hacking,
but it was placed in such a crazy place that you probably never saw it.
Ellipses, HO!
Now we can move on to eclipses. Since ellipses are simply a
squashed circle, it seems reasonable that we could modify the above circle
algorithm. So, let's get to it!
Everyone knows the equation of an eclipse:
x^2/a^2 + y^2/b^2 = 1
Upon taking differentials of both sides we have,
2*x*dx/a^2 + 2*y*dy/b^2 = 0
or, equivalently,
dy = -b^2/a^2 * x/y * dx
As you can see, life becomes suddenly becomes more complicated by a factor of
b^2/a^2. Furthermore, with an eclipse we only have reflection symmetries
through the x- and y-axis. In the circle algorithm we could get away with
just drawing an eighth of the circle, but now we have to draw a full quarter
of the eclipse.
We will start drawing the eclipse at x=0, y=b, so that initially x
will increase by one at each step, and y will wait a few steps to increase.
At some point, though, we will want y to increase by one at each step, and
x to wait a few steps before increasing; in the circle algorithm we just quit
once we reached this point, but now we are going to need an equation for dx:
dx = -a^2/b^2 * y/x * dy
In the circle algorithm, we used a single variable to count up and
tell us when it was time to increase y. Perhaps your intuition suggests
that we can do an eclipse with _two_ variables; mine said the same thing,
so that is exactly what we will do.
First, let us assume we have a way of calculating b^2/a^2:
E = b^2/a^2
I will suggest a way to perform this calculation later. Let's write out
the first few terms in the dy summation, starting at x=0, y=b:
dy1 + dy2 + ... = -E * (x0 + x1 + x2 + x3 + ...)/y
= -E * (0 + 1 + 2 + 3 + ...)/b
= - (0 + E + 2E + 3E + ...)/b
So, the basic structure of the algorithm is: add up 0, E, 2E, etc. until
the sum is larger than y. At that point, reset the counter, keeping the
remainder, and decrease y. This is where the two variables come in:
X=X+1
T2=T2+E
T1=T1+T2
IF T1>=Y THEN T1=T1-Y:Y=Y-1
Do you see how it works? T2 simply takes on the values 0, E, 2E, 3E, etc.,
and T1 is the counter. Furthermore, you can see that once T2 is larger
than Y, dy will be larger than one at each step. We need a new algorithm
to continue the calculation, and it turns out to be quite simple.
Look at the expression for dx above. We could calculate a^2/b^2,
but somehow that goes against the spirit of the calculation so far. Let's
instead rewrite dx slightly:
dx = - y/(E*x) * dy
Here we have simply written a^2/b^2 as 1/(b^2/a^2) = 1/E. But E*x is
exactly the variable T2 above, so we can continue the calcuation without
even stopping for breath:
Y=Y-1
T1=T1+Y
IF T1>=T2 THEN T1=T1-T2:X=X+1:T2=T2+E
(remember that T1 keeps track of the fractional part of y). So, we now
have a complete algorithm for drawing an eclipse:
0 REM ELLIPSE ATTEMPT #N SLJ 11/3/95
10 A=150:B=16:E=B*B/(A*A)
20 X=0:Y=B:T1=0:T2=0.5
30 GRAPHIC1,1:SLOW:X0=160:Y0=100:DRAW1,X0+A,Y0:DRAW1,X0,Y0-B
40 X=X+1:T2=T2+E
50 T1=T1+T2
60 IF T1>=Y THEN T1=T1-Y:Y=Y-1
70 DRAW1,X0+X,Y0-Y
80 IF T2=T2 THEN T1=T1-T2:X=X+1:T2=T2+E
120 DRAW1,X0+X,Y0-Y
130 IF Y>0 THEN 90
Lines 40-80 are the top part of the eclipse, and lines 90-130 handle the
bottom part. Note that T2 starts at 0.5, to round off the calculation in the
same spirit as we did in the circle algorithm.
Naturally, this algorithm has a few limitations. In line 30 the start
and end points are plotted, so you can see how close the algorithm really is.
In my experiments it occasionally missed the endpoint by a pixel or two. As
usual, I was a little too lazy to investigate possible ways to get around this.
If you require a perfect eclipse, you need to start the calculation at x=0, y=b
and run it forwards (e.g. lines 40-80 above), and then do another, similar
calcuation, starting at x=a, y=0, and running backwards. That is, for the
second calculation, calculate E2=a^2/b^2, and then run the algorithm just like
lines 40-80, interchanging X and Y.
Now we need to translate this algorithm into assembly. I am going
to make a few assumptions: first, that everything fits in a byte. In
particular, I require that b^2/a < 256. This insures that b^2/a^2 < 256,
and also insures that T2 will not overflow (note that when x=a, T2=E*a,
e.g. T2=b^2/a). What this means is that eclipses can't be too squashed.
Next, we need to deal with the fraction E=b^2/a^2. Any number
like this consists of two parts, an integer part plus a fractional part
(e.g. a number and a decimal). So, let's split E into two parts, EL and EH,
where EL represents the decimal part and EH the integer. Now our addition
consists of adding together the fractional parts, and if there is an overflow,
increasing the integer part. For example, if E=1.62, then EH=1 and EL=0.62.
We add EL to our number, and if it is greater than one, we carry the one to
when we add EH to our number.
The best thing to do is to represent EL as a fractional part of 256.
That is, our EL above should really be 0.62*256. This way, carries and
overflows will be handled automatically (this will become clear in a moment).
Let me give some pseudo-assembly code and we'll push off the
explanation until later:
35 GOTO 200
190 REM ***********************
200 XM=0:YM=B:X=128:Y=0:EH%=INT(E):EL%=INT((E-EH%)*256+0.5)
210 XM=XM+1
220 C=0:A=X:A=A+EL%:IF A>255 THEN A=A-256:C=1
230 X=A:A=Y:A=A+EH%+C:Y=A
235 A=A+T1
240 IF A>=YM THEN A=A-YM:YM=YM-1
250 T1=A:DRAW1, X0+XM, Y0-YM
260 IF Y<=YM THEN 210
265 T2=Y:A=T1
270 YM=YM-1
280 A=A+YM:IF A255 THEN A=A-256:C=1
295 X=A:A=T2:A=A+EH%+C:T2=A:A=T1
300 DRAW1, X0+XM, Y0-YM
310 YM=YM-1:IF YM>=0 THEN 280
XM and YM are the x and y coordinates of the point to be plotted. Note
that in line 200 X starts at 128, and this again is to round up all our
calculations; compare to line 20, where we started T2 at 0.5. In the
above code I store T2 in the X and Y registers for the first part of the
code. Note that in lines 220 and 290 there is some extraneous code to
simulate things that in assembly are taken care of by the 6502. Note
also that the comparison in line 260 has been changed from < to <=. This
makes the branch easier, and I'm not sure how it affects the calculation
(I didn't notice any difference in the few runs I tried it on).
Moving through the code, we increase x, and then add the decimal
part of E to the counter. Then we add the integer part of E to the counter,
along with any carries. If the integer part of the counter is greater than
y, it is time to decrease y and reset the counter.
Moving to the second part of the code, we do a little rearranging
in line 265. Really a better thing to do would be to let A=T1-T2, so that
the compare in line 280 becomes simpler. Anyways, note that the Y register
becomes freed up at this point. From here on, it is pretty much the same
thing as before.
The full assembly code is then:
;Ellipse SLJ 11/3/95 Assumptions:
;0->XM B->YM, x- and y-coordinates
;0->T1
;EL and EH contain remainder and integer parts of E, resp.
LDX #128
LDY #00
CLC
L1 INC XM
TXA
ADC EL
TAX
TYA
ADC EH
TAY
ADC T1
CMP YM
BCC :CONT1
SBC YM
DEC YM
:CONT1 STA T1
JSR PLOT
CPY YM
BCC L1
STY T2
LDA T1
SBC T2
DEC YM
L2 ADC YM
BCC :CONT2
SBC T2
STA T1
INC XM
TXA
ADC EL
TAX
LDA T2
ADC EH
STA T2
LDA T1
:CONT2 JSR PLOT
DEC YM
BPL L2 ;Assuming y<128
Logarithms
Finally, we need a way of calculating b^2/a^2. I suggest using
logarithms for this. I do believe I discussed this concept in an earlier
issue of C=Hacking. Nevertheless, the idea is that if
x = b^2/a^2
then
log(x) = 2*log(b) - 2*log(a)
so that
x = exp(2*(log(b) - log(a))
Thus, three tables need to be created: one for log(x), and one each for
the integer and remainder parts of e^(2*x). Now, to improve accuracy,
the first table might be a table of f(x)=222/log(128) * log(x/2). This
constant is chosen so that f(255) is roughly 255. 222 was chosen because
the inversion (i.e. the e^x part) works best at that value. This pretty
much assumes that x is not zero or one, either. You can of course use
more tables for somewhat better accuracy.
One really nice thing about this is that you don't have to worry
about squaring things, since that part can be taken care of automatically
in logarithm space. On the downside, we are restricted even further by
the types of numbers we can divide (e.g. log(a)-log(b) can't be larger
than 127 or so).
Division then consists of a table lookup, subtraction of another
table lookup, and two more table lookups. Here is a short program to
demonstrate the use of logs in this sort of division, and a very rough
feel for the type of accuracy to expect -- note that it doesn't compare
decimal parts, or convert the decimal parts into fractions of 256, etc.:
1 FAST:PRINT"[CLR]"
10 DIM L(256),EI(256),ER(256):FC=222/LOG(128)
20 FOR I=1 TO 256
25 PRINT "[HOME]"I
30 L(I)= INT(FC*LOG(I/2)+0.5):IF I=1 THEN L(I)=0
40 S=I:IF I>127 THEN S=I-256
50 EX=EXP(2*S/FC):IF EX>256 THEN PRINT"WHOOPS! EX="EX"I="I
60 EI(I)=INT(EX+0.5)
70 ER(I)=EX-EI(I)
80 NEXT I
90 EI(0)=1:ER(0)=0
100 FOR A=2 TO 250
110 FOR B=2 TO 250
120 X=L(B)-L(A)
123 IF X>127 THEN PRINT"OOPS:A="A"B="B"X="X
126 IF X<0 THEN X=X+256
130 A1=EI(X)+ER(X):A2=B*B/(A*A):IF A2>255 THEN B=600
135 BL=INT(A2+0.5)-INT(A1+0.5)
140 PRINT A;B,A1;A2,"ERR="INT(A2+0.5)-INT(A1+0.5)
150 NEXT:NEXT
Conclusion
Sorry, no 3D graphics this time around. Watch for a full-screen, hires
bitmapped solid 3D virtual world sometime in the not too distant future.
Otherwise, may your ellipses never be square :).
Last Updated: 1995-12-4 Rev A