First of all, there's some mistake in your code, as the sin()
and cos()
functions are radian and not degree based functions. This will make your circle to appear at pseudorandomly selected points as each step draws a point about 57 degrees apart from the previous. This will have no impact as modern graphics work buffered and you'll see the final result, not the workings.
Once this said, nobody draws a circle today by using your exposed algorithm. Have a look at the Bresenham's middle point algorithm, that basically tries to draw a circle by octants, but several times faster.
The idea behind these algorithms is to consider the R^2 = x^2 + y^2
formula and to draw, pixel by pixel, on one of the axis, and considering when the other axis must be followed (you draw by octants, as this way you don't cope with greater than one derivatives and you have only to decide if you move up or not). The routine also takes into account the circle symmetry to calculate only one octant and then draw the eight points at each pass.
As I developed that algorithm from scratch when I was young (without having seen Bresenham's before) probably my reasoning up to the solution will be of help to you.
The first attempt is to take into account that the resolution (the granularity is not angle dependant) for small circles you have to paint less pixels than for big ones, and the one degree approach you have followed has to be redesigned to adapt to finer or coarse resolutions. The idea is to go, pixel by pixel, instead of degree by degree, until you draw the complete thing. We are going to paint only the first octant, and will draw the rest by the symmetry properties of the figure. We part from the (0, -R)
point and will go, pixel by pixel until we get to the (sqrt(2)*R, R - sqrt(2)*R)
point.
The first thing we are going to do is to try to save all operations we have to do. The first place where we can save operations is in calculating squares... We are going to use the R^2 = x^2 + y^2
equation and on it, R
is only used as R^2
all the time, so, lets suppose we want to draw a ten pixels radius circle, we will square things up to 100
(which is the 10 pixels radius squared).
Next, we are going to see one property of squares, that is, they grow odd from one square to the next (0 -> 1(delta is 1) -> 4(delta is 3) -> 9(delta is 5) -> 16(delta is 7) ...
) so if we can arrange to grow by 1 in the x, we can calculate x^2 easily, by just adding two to the odd
variable, and then adding odd
to the last square number, so we'll use two numbers: x
and x2
. We initialize both to 0
, and the first grows by x += 1;
, while the second grows by the relationship x2 += dx2; dx2 += 2;
(we initialize dx2 = 1;
) this makes us to allow x
and x2
growing by only making sums, and no multiplications at all.
If someone thinks we are going to need y2 = 100 - x2
and then whe are forced to calculate y = sqrt(y2)
you are almost right, but the trick here is to be able to manage the y
and y2
sequences backwards the same as the x counterparts. Well, right, y
and y2
can be managed in the reverse direction the same as x
and x2
but this time we have to go backwards, decreasing odd numbers from (what?) to 1
where dy2 -= 2; y2 -= dy2;
and finally reaches 0
. For this to happen, check that the difference between two consecutive squares is precisely the two numbers added, so, for example, the difference between 13^2 = 169
and 14^2 = 196
is 13 + 14 = 27
, and this is the odd number to begin with if we go back from R = 14
to 0
in y
.
The reason of complicating things so much, is that this way we only make additions with integers and no need to make multiplications (well, multiplications are not so expensive, but there was a time that they were) Well, we do multiplications to square initially the radius R
, but we do it only once at the beginning to calculate R^2
.
The idea now is to set the origin at the departure point (0, -R)
and go to the right, pixel by pixel, adding (and modifiying) x
, x2
, and sum
(we substract to sum all the time) until we reach the next square in y
, and then update all the y axis values (we do decrement the y, we have to move a pixel up in that moment) y
, y2
, dy2
, and draw the pixel (or draw it before, as we do in the routine), until... what? (well, the point is until we meet at the 45 degrees point in which the octant is complete, x
and y
coordinates are equal) It is important to stop there, because from that point on, it is possible that one step makes the y coordinate to increment more than one pixel (the derivative is greater than one) and this should complicate the overal algorithm (and we are anyway painting the other symmetrical eight points, so we are drawing the other part of the graphic)
So, suppose we have 100 as radius, and begin with:
x=0, x2= 0, dx2= 1, y=10, y2=100, dy2=19, sum=100 *
x=1, x2= 1, dx2= 3, y= 9, y2= 81, dy2=17, sum= 99
x=2, x2= 4, dx2= 5, y= 9, y2= 81, dy2=17, sum= 96
x=3, x2= 9, dx2= 7, y= 9, y2= 81, dy2=17, sum= 91
x=4, x2=16, dx2= 9, y= 9, y2= 81, dy2=17, sum= 84 *
x=5, x2=25, dx2=11, y= 8, y2= 64, dy2=15, sum= 75 *
x=6, x2=36, dx2=13, y= 7, y2= 49, dy2=13, sum= 64 *
x=7, x2=49, dx2=15, y= 7, y2= 49, dy2=13, sum= 51
The point marked with asterisk are the ones where the sum
value crosses the next y2
value, making the y
value to be decrmented and having to shift the pixel we paint on. The final routine is:
int bh_onlyonepoint(int r, int cx, int cy)
{
int r2 = r*r;
int x = 0, x2 = 0, dx2 = 1;
int y = r, y2 = y*y, dy2 = 2*y - 1;
int sum = r2;
while(x <= y) {
draw(cx + x, cy + y); /* draw the point, see below */
sum -= dx2;
x2 += dx2;
x++;
dx2 += 2;
if (sum <= y2) {
y--; y2 -= dy2; dy2 -= 2;
}
} /* while */
return x; /* number of points drawn */
}
If you want, I have written a simple example to draw a circle in ascii on the screen, given radius as command parameter. It uses ANSI escape sequences to position the cursor in the screen before drawing a single asterisk. The scale is doubled in the X direction to compensate for the character height ("pixels" in ascii are not square) I have included a new drawing function pointer parameter to call back for point drawing and a main
routine to get parameters from the command line:
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
void draw(int x, int y)
{
/* move to position (2*x, y) and plot an asterisk */
printf("\033[%d;%dH*", y, x<<1);
}
int bh(int r, int cx, int cy, void(*draw)(int, int))
{
/* the variables mentioned in the text */
int r2 = r*r;
int x = 0, x2 = 0, dx2 = 1;
int y = r, y2 = y*y, dy2 = 2*y - 1;
int sum = r2;
while(x <= y) {
/* draw the eight points */
draw(cx + x, cy + y);
draw(cx + x, cy - y);
draw(cx - x, cy + y);
draw(cx - x, cy - y);
draw(cx + y, cy + x);
draw(cx + y, cy - x);
draw(cx - y, cy + x);
draw(cx - y, cy - x);
sum -= dx2;
x2 += dx2;
x++;
dx2 += 2;
if (sum <= y2) {
y--; y2 -= dy2; dy2 -= 2;
}
} /* while */
return x; /* number of points drawn */
}
int main(int argc, char **argv)
{
int i;
char *cols = getenv("COLUMNS");
char *lines = getenv("LINES");
int cx, cy;
if (!cols) cols = "80";
if (!lines) lines = "24";
cx = atoi(cols)/4;
cy = atoi(lines)/2;
printf("\033[2J"); /* erase screen */
for (i = 1; i < argc; i++) {
bh(atoi(argv[i]), cx, cy, draw);
}
fflush(stdout);
sleep(10);
puts(""); /* force a new line */
}
And the final result is:
*
* * * * * * * * * *
* * * *
* *
* * * *
* * *
* * * * * * * * * *
* * * *
* * * * * *
* * * *
* * * *
* * * *
* * * *
* * * *
* * * *
* * * *
* * * *
* * * *
* * * *
* * * *
* * * *
* * * *
* * * * * *
* * * *
* * * * * * * * * *
* * *
* * * *
* *
* * * *
* * * * * * * * * *
*
Finally, if you want better results (not those peaks that arise from an exact value of radius that make them only touch the point when x or y are zero) you can pass the routine directly the square value of the radius (that allows to make integer calculations with fractional radii)
filling a circle
If you want to fill the circle, just paint all points between a pair of the calculated points, as in:
lineFromTo(cx - x, cy - y, cx + x, cy - y);
lineFromTo(cx - y, cy + x, cx + y, cy + x);
lineFromTo(cx - y, cy - x, cx + y, cy - x);
lineFromTo(cx - x, cy + y, cx + x, cy + y);
These are all horizontal lines, so perhaps you can get an improvement with something like:
/* X1 X2 Y */
HorizLineX1X2Y(cx - x, cx + x, cy - y);
HorizLineX1X2Y(cx - y, cx + y, cy + x);
HorizLineX1X2Y(cx - y, cx + y, cy - x);
HorizLineX1X2Y(cx - x, cx + x, cy + y);
A new git repository with the final program allowing to fill, draw or trace the run of the algorithm has been created in github