Computing circles without Sin() or Cos()

May 22, 2019
651
312
63
So I was thinking about how I'd draw a circle, but without using the SIN and COS functions in the software. Instead, I'm using the Pythagorean theorem.

So what do we know about circles? A circle is the graph of all points that are equidistant from the center point. In other words, all of the points on a circle are the same distance from the center.

In graphics, any time we talk "distance", we might invoke Pythagoras. Anyone who has taken high school geometry probably remembers the formula:
a²+b²=c²
Re-written for a 2D graphics plane, we can say
x²+y²=d²
where d is the distance from the starting point, and x and y are the horizontal and vertical distance between the start and end points.

Normally, we solve this for d, getting the distance to a point, but what if we wanted to solve this for the x component? So I transformed the formula into...
d² = x² + y²
d² - y² = x² + y² - y²
d² - y² = x²
x² = d² - y²
x = √(d²-y²)
finally, we convert this to program code:
x = sqrt(d*d - y*y)

What this actually gives us is an arc...., like this:
1573067716732.png


So how to turn that into a circle? We just plot the same point by negating the x offset, then repeat with the y offset negated:

x = sqrt(d*d - y*y)
x = -sqrt(d*d - y*y)
x = sqrt(d*d - -y*-y)
x = -sqrt(d*d - -y*-y)

and some BASIC code:
1573067970048.png

This creates a filled circle, by drawing from the X axis out to the edge. If we want a filled circle, we should instead draw a line from the previous point to the next one. This prevents the broken arc we saw in the first image:

1573068386324.png


Here's the program listing:
Code:
10 SCREEN $80
20 RECT 0,0,319,199,12
30 X=160:Y=100:R=75:C=2
100 FOR YO=0 TO R
110 XO = SQR(R*R - YO*YO)
130 LINE X-XO,Y+YO,X+XO,Y+YO,C
140 LINE X+XO,Y-YO,X-XO,Y-YO,C
170 LX=XO:LY=YO
190 NEXT
200 C=5
210 LX=R:LY=Y
220 FOR YO=0 TO R
230 XO = SQR(R*R - YO*YO)
240 LINE X+LX,Y+YO,X+XO,Y+YO,C
250 LINE X-LX,Y+YO,X-XO,Y+YO,C
260 LINE X+LX,Y-YO,X+XO,Y-YO,C
270 LINE X-LX,Y-YO,X-XO,Y-YO,C
280 LX=XO:LY=YO
290 NEXT
300 C=1
310 PSET X,Y,C
And, for fun, here is an animated GIF showing the whole sequence:

circle.gif
 
Last edited:
  • Like
Reactions: mobluse
May 22, 2019
651
312
63
So I wanted to look closer at the broken arc in my first figure. Let's take another look at that:



What's going on here is that once we pass the 45° mark, the X offset is changing by more than 1 unit per row. This is entirely expected, as 45° is the point at which the X and Y offsets match.

In my first example, I saved the last point and drew a line from that to the next one, to guarantee contiguous lines. There is another method, one that uses more lines of code, but is a little faster.

So in the first example, we draw 4 pixels for each calculation by negating the X and Y offsets, in turn:

x = sqrt(d*d - y*y)
x = -sqrt(d*d - y*y)
x = sqrt(d*d - -y*-y)
x = -sqrt(d*d - -y*-y)

Of course, there's no rule saying we can't also solve this equation for Y. So what happens if we solve for Y, instead of X?

y = sqrt(d*d - x*x)


We get the same broken circle, but turned 90°
1573074311132.png


So can we take advantage of this?

As it turns out, we can... we can draw half of the arc by solving for X, then draw half of the arc solving for Y. But the trick is to figure out where the 45° mark actually is.... it's not at 50% down. It's actually at 71% down, or more precisely, 0.707.

Here's a graphic illustration of the 45° point:
1573075070179.png

What this means is we need to draw 71% of the way down solving for X, then draw 71% of the way across solving for Y.

That results in two loops, which look something like this:
Code:
R1=R * 0.707
FOR YO=0 TO R1
  XO = SQR(R*R - YO*YO)
  PSET X+XO,Y+YO,C
NEXT
FOR XO=0 TO R1
  YO = SQR(R*R - XO*XO)
  PSET X+XO,Y+YO,C
NEXT
That, indeed draws a full 90° arc.
1573075401515.png

But the two loop approach is a little ugly... so let's combine the two loops into one block of code

Code:
R1=R * 0.707
FOR YO=0 TO R1
  XO = SQR(R*R - YO*YO)
  PSET X+XO,Y+YO,C
  PSET X+=YO,Y+XO,C
NEXT
And when we then draw all four quadrants using this method, we get the program below

Code:
10 SCREEN $80
20 RECT 0,0,319,199,15
30 X=160:Y=100
40 R=75:R1=R*0.71
300 C=0
305 T1=TI
310 FOR YO=0 TO R1
320 XO = SQR(R*R - YO*YO)
330 PSET X+XO,Y+YO,C
340 PSET X-XO,Y+YO,C
350 PSET X-XO,Y-YO,C
360 PSET X+XO,Y-YO,C
370 PSET X+YO,Y+XO,C
380 PSET X-YO,Y+XO,C
390 PSET X-YO,Y-XO,C
400 PSET X+YO,Y-XO,C
410 NEXT
420 T2=TI
430 PRINTT2-T1
500 C=1
510 PSET X,Y,C
 
  • Like
Reactions: JayLewis

Schol-R-LEA

New Member
Sep 20, 2019
20
15
3
Perhaps you aren't aware (I am guessing that you are and are just using the opening paragraph as a framing device) that generally speaking, no one has been using trig functions in drawing complete 2D circles since the early 1970s (partial arcs are another matter, and even then they are only used for the endpoints).

The algorithm you are describing is similar to the Midpoint Circle Algorithm (explained somewhat less abstractly on Geeks4Geeks and given in several languages on Rosetta Code), except that the MPC manages to avoid the need for the square roots as well, and the most commonly used version even works with integer-only math.

I can try to get a version of MPC (and the more general Bresenham's ellipse algorithm) written RSN, if you like.
 
May 22, 2019
651
312
63
At this point, I’m just experimenting with trig and geometry. As it turns out, a 32 point ellipse with SIN() and COS() is still faster to plot than even using Bresenham.

I was really mostly just interested in the trigonometry, and solving for a and b grew out of that. I’ve actually written several implementations of the circle algorithm, but those are moot, since the intent was to demonstrate some algebra, rather than MPC or Bresenham. The thing I don’t like about those algorithms is they kind of divorce the math from the underlying geometry, which is fine when you’re optimizing programs for speed, but not when you’re explaining math and basic computing science.

I guess this is really an answer to a question asked elsewhere... why do I need to know algebra (and other maths) to write computer software? This is one application. This is a super simple use of algebra, but it’s an important one.
 
Last edited:

JoshuaScholar

New Member
Nov 6, 2019
24
10
3
I vaguely remember there was a trick that would draw almost circles, probably without any multiplies. They were SLIGHTLY elliptical diagonally.
 

gertk

New Member
Oct 7, 2019
6
9
3
Try searching for Bresenham circles, they can be done with integer math

 
May 22, 2019
651
312
63
Try searching for Bresenham circles, they can be done with integer math

Yes, I tried Bresenham, too, but while it's roughly 2x as fast as using the Pythagorean method, it's not very useful to demonstrate algebra, trigonometry, or geometry, which is more what I was going for.

Regardless, the SIN/COS method is still faster than any other method I've tried, depending on how many vertices I plot. (I usually go with 32.)
 

Schol-R-LEA

New Member
Sep 20, 2019
20
15
3
Regardless, the SIN/COS method is still faster than any other method I've tried, depending on how many vertices I plot. (I usually go with 32.)
I'm really quite surprised by that, especially for a system without dedicated hardware FP. My understanding has always been that Bresenham - and the somewhat more complex Wu - is significantly faster; it is often cited that Bresenham is usually included in most GPU's firmware (or even hardware) because it is still useful for 2D acceleration in windowing environments even today.

I seem to recall (though I don't know the source offhand - maybe somewhere in the DDJ back issue collections? I read a bunch of those in college during the early 1990s) it being stated that, for an assembly implementation on the 6502 (or maybe it was Z80, not sure), the integer version of Bresenham would complete before a single sine calculation similarly implemented would. While that claim seems dubious, the point that trig functions done purely in software are very expensive still makes your claim surprising.

Would you mind sharing your implementations and benchmark results, please? I would interested in see how it works out.

One thing I will note is that for an implementation in interpreted BASIC, it might make some more sense - in the case using the trig functions, most of the heavy lifting is done by the ROM routines, whereas the Bresenham implementation would spend most of its time in the interpreter itself processing each line of the core loop (even with the lines of code pre-tokenized by the editor, the interpretation time for each line is significant).
 
May 22, 2019
651
312
63
I'm really quite surprised by that, especially for a system without dedicated hardware FP.

Would you mind sharing your implementations and benchmark results, please? I would interested in see how it works out.
I believe MS BASIC uses a lookup table for the SIN/COS stuff. I did something similar in a flight simulator I wrote a few years ago: I used a lookup table for each 1° increment, then interpolated the value for in-between numbers.

My Bresenham version runs in 28 ticks (0.46 seconds), and my SIN/COS version takes 19 ticks to plot a circle using 32 line segments. At 48 segments, it breaks even at 28 ticks, but anything more than 32 seems to be wasted time (there just aren't enough pixels at that resolution to see more curvature with > 32 points.)

Here's the SIN/COS circle, with some benchmark code wrapped around it:
Code:
5 CP=8
10 X=160:Y=100
20 R=75
30 EA={PI}*2
40 RD=EA/CP
50 SCREEN $80
90 TI$="000000"
100 YO=R:XO=0:FOR A=0 TO EA STEP RD
110 LX=XO:XO=SIN(A)*R
120 LY=YO:YO=COS(A)*R
130 LINE X+LX,Y+LY,X+XO,Y+YO
140 NEXT
150 PRINT CP,TI:INPUT A$
160 CP=CP+8:IF CP<=360 GOTO 10
1574362457799.png

and here's my Bresenham version:
Code:
0 TI$="000000"
10 SCREEN $80
20 X=160:Y=100
100 R=75
110 XO=0:YO=R
120 DP=3-2*R
130 GOTO 200
140 XO=XO+1
150 IF DP<0 THEN 190
160 YO=YO-1
170 DP=DP+4*(XO-YO)+10
180 GOTO 200
190 DP=DP+4*XO+6
200 PSET X+XO,Y+YO
210 PSET X-XO,Y+YO
220 PSET X+XO,Y-YO
230 PSET X-XO,Y-YO
240 PSET X+YO,Y+XO
250 PSET X-YO,Y+XO
260 PSET X+YO,Y-XO
270 PSET X-YO,Y-XO
280 IF YO>=XO THEN 140
999 PRINT TI
1574362523847.png


Finally, here's a version that will superimpose the two circles, so you can compare the curve. I also made the circle as large as possible, but which would still fit on the screen. This should tell us how well the algorithms scale. Note that we only gained one tick with the SIN/COS but gained 10 ticks with the Bresenham. This is more the fault of BASIC than the Bresenham algorithm.

The BASIC portion of the SIN/COS version actually takes exactly the same amount of time, no matter how big the circle is, since this version only computes 32 points. So the only slowdown is in the amount of time it takes the runtime to draw the line segments. This only holds true to a point, of course, since on a higher resolution display, you'd want to compute more vertices.

On the other hand, the Bresenham algorithm has to compute and plot every point, so it scales linearly with the radius. When the radius is 50, it takes 19 ticks. When the radius is 100, it takes 38 ticks.



1574365921141.png

Note that the curve is a little more pronounced at the 45° mark in the SIN/COS plot. I suspect the differences here are due to rounding issues and not either of the algorithms themselves. Also note that the Bresenham algorithm scales

Code:
10 SCREEN $80
20 REM PYTHAGEREAN CIRCLE
30 X=160:Y=100:R=99:C=7
40 FOR YO=0 TO R
50 XO = SQR(R*R - YO*YO)
60 REM PSET X+XO,Y+YO,C
70 LINE X-XO,Y+YO,X+XO,Y+YO,C
80 LINE X+XO,Y-YO,X-XO,Y-YO,C
90 LX=XO:LY=YO
100 NEXT
1000 REM SIN/COS CIRCLE
1010 CP=32
1020 EA={PI}*2
1030 RD=EA/CP
1040 TI$="000000"
1050 YO=R:XO=0
1060 FOR A=0 TO EA STEP RD
1070 LX=XO:XO=SIN(A)*R
1080 LY=YO:YO=COS(A)*R
1090 LINE X+LX,Y+LY,X+XO,Y+YO,5
1100 NEXT
2000 PRINT TI
2010 REM DRAW TICKS (TO SHOW VERTICES OF SIN/COS CIRCLE)
2020 YO=R:XO=0
2030 FOR A=0 TO EA STEP RD
2040 XO=SIN(A)*R
2050 YO=COS(A)*R
2060 LINE X+XO,Y+YO,X+XO*0.9,Y+YO*0.9,15
2070 NEXT
3000 REM BRESENHAM CIRCLE
3010 TI$="000000"
3020 XO=0:YO=R
3030 DP=3-2*R
3040 GOTO 3110
3050 XO=XO+1
3060 IF DP<0 THEN 3100
3070 YO=YO-1
3080 DP=DP+4*(XO-YO)+10
3090 GOTO 3110
3100 DP=DP+4*XO+6
3110 PSET X+XO,Y+YO
3120 PSET X-XO,Y+YO
3130 PSET X+XO,Y-YO
3140 PSET X-XO,Y-YO
3150 PSET X+YO,Y+XO
3160 PSET X-YO,Y+XO
3170 PSET X+YO,Y-XO
3180 PSET X-YO,Y-XO
3190 IF YO>=XO THEN 3050
3200 PRINT TI
3210

Finally, just for fun, here's another approximation algorithm that runs slightly faster than Bresenham. It uses pure addition, but uses floating point numbers, rather than integers.

Code:
10 SCREEN $80
20 X=160:Y=100:R=99:R1=R/2
30 TI$="000000"
40 YO=R:DP=0:XO=0:Y1=0:Y2=0.007074:Y3=0.0002
50 FOR XO=0 TO 70
60 YO=YO-Y1:Y1=Y1+Y2:Y2=Y2+Y3
70 PSET X+XO,Y+YO
80 PSET X-XO,Y+YO
90 PSET X+XO,Y-YO
100 PSET X-XO,Y-YO
110 PSET X+YO,Y+XO
120 PSET X-YO,Y+XO
130 PSET X+YO,Y-XO
140 PSET X-YO,Y-XO
150 NEXT
160 PRINT TI
There would probably need to be some math on the various offsets (Y2 and Y3) to adjust for different circle sizes, but this is faster than Bresenham only because of quirks in BASIC, not because the algorithm itself is superior.
1574369784450.png
 

Attachments

777ismyname

New Member
Nov 26, 2019
1
1
3
Thank you for the code examples, TomXP411. I was working on something yesterday in the integer only version of Fast BASIC on the Atari 8 bit machines and this will help me a good bit.
 
  • Like
Reactions: TomXP411