Midpoint circle algorithm

Rasterising a circle by the Bresenham algorithm

In computer graphics, the midpoint circle algorithm is an algorithm used to determine the points needed for rasterizing a circle. Bresenham's circle algorithm is derived from the midpoint circle algorithm. The algorithm can be generalized to conic sections.

The algorithm is related to work by Pitteway and Van Aken.

Summary

This algorithm draws all eight octants simultaneously, starting from each cardinal direction (0°, 90°, 180°, 270°) and extends both ways to reach the nearest multiple of 45° (45°, 135°, 225°, 315°). It can determine where to stop because when y = x, it has reached 45°. The reason for using these angles is shown in the above picture: As y increases, it does not skip nor repeat any y value until reaching 45°. So during the while loop, y increments by 1 each iteration, and x decrements by 1 on occasion, never exceeding 1 in one iteration. This changes at 45° because that is the point where the tangent is rise=run. Whereas rise>run before and rise<run after.

The second part of the problem, the determinant, is far trickier. This determines when to decrement x. It usually comes after drawing the pixels in each iteration, because it never goes below the radius on the first pixel. Because in a continuous function, the function for a sphere is the function for a circle with the radius dependent on z (or whatever the third variable is), it stands to reason that the algorithm for a discrete(voxel) sphere would also rely on this Midpoint circle algorithm. But when looking at a sphere, the integer radius of some adjacent circles is the same, but it is not expected to have the same exact circle adjacent to itself in the same hemisphere. Instead, a circle of the same radius needs a different determinant, to allow the curve to come in slightly closer to the center or extend out farther. The circle charts seen relating to Minecraft, like the determinant listed below, only account for one possibility.

Algorithm

The objective of the algorithm is to find a path through the pixel grid using pixels which are as close as possible to solutions of $x^{2}+y^{2}=r^{2}$. At each step, the path is extended by choosing the adjacent pixel which satisfies $x^{2}+y^{2}\leq r^{2}$ but maximizes $x^{2}+y^{2}$. Since the candidate pixels are adjacent, the arithmetic to calculate the latter expression is simplified, requiring only bit shifts and additions.

This algorithm starts with the circle equation. For simplicity, assume the center of the circle is at $(0,0)$. Consider first the first octant only, and draw a curve which starts at point $(r,0)$ and proceeds counterclockwise, reaching the angle of 45.

The fast direction here (the basis vector with the greater increase in value) is the $y$ direction. The algorithm always takes a step in the positive $y$ direction (upwards), and occasionally takes a step in the slow direction (the negative $x$ direction).

From the circle equation is obtained the transformed equation $x^{2}+y^{2}-r^{2}=0$, where $r^{2}$ is computed only once during initialization.

Let the points on the circle be a sequence of coordinates of the vector to the point (in the usual basis). Points are numbered according to the order in which drawn, with $n=1$ assigned to the first point $(r,0)$.

For each point, the following holds:

${\begin{aligned}x_{n}^{2}+y_{n}^{2}=r^{2}\end{aligned}}$

This can be rearranged thus:

${\begin{aligned}x_{n}^{2}=r^{2}-y_{n}^{2}\end{aligned}}$

And likewise for the next point:

${\begin{aligned}x_{n+1}^{2}=r^{2}-y_{n+1}^{2}\end{aligned}}$

Since the next point will always be at least 1 pixel higher than the last, it is true that:

${\begin{aligned}y_{n+1}^{2}&=(y_{n}+1)^{2}\\&=y_{n}^{2}+2y_{n}+1\end{aligned}}$

${\begin{aligned}x_{n+1}^{2}=r^{2}-y_{n}^{2}-2y_{n}-1\end{aligned}}$

So, rework the next-point-equation into a recursive one by substituting $x_{n}^{2}=r^{2}-y_{n}^{2}$:

${\begin{aligned}x_{n+1}^{2}=x_{n}^{2}-2y_{n}-1\end{aligned}}$

Because of the continuity of a circle and because the maxima along both axes is the same, clearly it will not be skipping x points as it advances in the sequence. Usually it stays on the same x coordinate, and sometimes advances by one.

The resulting coordinate is then translated by adding midpoint coordinates. These frequent integer additions do not limit the performance much, as those square (root) computations can be spared in the inner loop in turn. Again, the zero in the transformed circle equation is replaced by the error term.

The initialization of the error term is derived from an offset of ½ pixel at the start. Until the intersection with the perpendicular line, this leads to an accumulated value of $r$ in the error term, so that this value is used for initialization.

The frequent computations of squares in the circle equation, trigonometric expressions and square roots can again be avoided by dissolving everything into single steps and using recursive computation of the quadratic terms from the preceding iterations.

Variant with integer-based arithmetic

Just as with Bresenham's line algorithm, this algorithm can be optimized for integer-based math. Because of symmetry, if an algorithm can be found that only computes the pixels for one octant, the pixels can be reflected to get the whole circle.

We start by defining the radius error as the difference between the exact representation of the circle and the center point of each pixel (or any other arbitrary mathematical point on the pixel, so long as it's consistent across all pixels). For any pixel with a center at $(x_{i},y_{i})$, the radius error is defined as:

$RE(x_{i},y_{i})=\left\vert x_{i}^{2}+y_{i}^{2}-r^{2}\right\vert$

For clarity, this formula for a circle is derived at the origin, but the algorithm can be modified for any location. It is useful to start with the point $(r,0)$ on the positive X-axis. Because the radius will be a whole number of pixels, clearly the radius error will be zero:

$RE(x_{i},y_{i})=\left\vert r^{2}+0^{2}-r^{2}\right\vert =0$

Because it starts in the first CCW positive octant, it will step in the direction with the greatest travel, the Y direction, so it is clear that $y_{i+1}=y_{i}+1$. Also, because it concerns this octant only, the X values have only 2 options: to stay the same as the prior iteration, or decrease by 1. A decision variable can be created that determines if the following is true:

$RE(x_{i}-1,y_{i}+1)<RE(x_{i},y_{i}+1)$

If this inequality holds, then plot $(x_{i}-1,y_{i}+1)$; if not, then plot $(x_{i},y_{i}+1)$. So, how to determine if this inequality holds? Start with a definition of radius error:

${\begin{array}{lcl}RE(x_{i}-1,y_{i}+1)&<&RE(x_{i},y_{i}+1)\\\left\vert (x_{i}-1)^{2}+(y_{i}+1)^{2}-r^{2}\right\vert &<&\left\vert x_{i}^{2}+(y_{i}+1)^{2}-r^{2}\right\vert \\\left\vert (x_{i}^{2}-2x_{i}+1)+(y_{i}^{2}+2y_{i}+1)-r^{2}\right\vert &<&\left\vert x_{i}^{2}+(y_{i}^{2}+2y_{i}+1)-r^{2}\right\vert \\\end{array}}$

The absolute value function does not help, so square both sides, since a square is always positive:

${\begin{array}{lcl}\left[(x_{i}^{2}-2x_{i}+1)+(y_{i}^{2}+2y_{i}+1)-r^{2}\right]^{2}&<&\left[x_{i}^{2}+(y_{i}^{2}+2y_{i}+1)-r^{2}\right]^{2}\\\left[(x_{i}^{2}+y_{i}^{2}-r^{2}+2y_{i}+1)+(1-2x_{i})\right]^{2}&<&\left[x_{i}^{2}+y_{i}^{2}-r^{2}+2y_{i}+1\right]^{2}\\\left(x_{i}^{2}+y_{i}^{2}-r^{2}+2y_{i}+1\right)^{2}+2(1-2x_{i})(x_{i}^{2}+y_{i}^{2}-r^{2}+2y_{i}+1)+(1-2x_{i})^{2}&<&\left[x_{i}^{2}+y_{i}^{2}-r^{2}+2y_{i}+1\right]^{2}\\2(1-2x_{i})(x_{i}^{2}+y_{i}^{2}-r^{2}+2y_{i}+1)+(1-2x_{i})^{2}&<&0\\\end{array}}$

Since x > 0, the term $(1-2x_{i})<0$, so dividing gets:

${\begin{array}{lcl}2\left[(x_{i}^{2}+y_{i}^{2}-r^{2})+(2y_{i}+1)\right]+(1-2x_{i})&>&0\\2\left[RE(x_{i},y_{i})+Y_{Change}\right]+X_{Change}&>&0\\\end{array}}$

Thus, the decision criterion changes from using floating-point operations to simple integer addition, subtraction, and bit shifting (for the multiply by 2 operations). If ⁠$2(RE+Y_{Change})+X_{Change}>0$⁠, then decrement the X value. If ⁠$2(RE+Y_{Change})+X_{Change}\leq 0$⁠, then keep the same X value. Again, by reflecting these points in all the octants, a full circle results.

C example

The above algorithm is implemented in the programming language C, below. Starting from $(r,0)$, it enumerates the points in counterclockwise order for the first octant, while simultaneously mirroring the points to the other octants.

void drawcircle(int x0, int y0, int radius)
{
    int x = radius-1;
    int y = 0;
    int dx = 1;
    int dy = 1;
    int err = dx - (radius << 1);

    while (x >= y)
    {
        putpixel(x0 + x, y0 + y);
        putpixel(x0 + y, y0 + x);
        putpixel(x0 - y, y0 + x);
        putpixel(x0 - x, y0 + y);
        putpixel(x0 - x, y0 - y);
        putpixel(x0 - y, y0 - x);
        putpixel(x0 + y, y0 - x);
        putpixel(x0 + x, y0 - y);

        if (err <= 0)
        {
            y++;
            err += dy;
            dy +=2;
        }
        if (err > 0)
        {
            x--;
            dx += 2;
            err += (-radius << 1) + dx;
        }
    }
}

JavaScript

Implementation that draws a circle in HTML5 canvas (for educational purposes only; there are better ways to draw circles in canvas).

const CHANNELS_PER_PIXEL = 4; //rgba

function drawCircle (x0, y0, radius, canvas) {
	var x = radius-1;
	var y = 0;
	var dx = 1;
	var dy = 1;
	var decisionOver2 = dx - (radius << 1);   // Decision criterion divided by 2 evaluated at x=r, y=0
	var imageWidth = canvas.width;
	var imageHeight = canvas.height;
	var context = canvas.getContext('2d');
	var imageData = context.getImageData(0, 0, imageWidth, imageHeight);
	var pixelData = imageData.data;
	var makePixelIndexer = function (width) {
		return function (i, j) {
			var index = CHANNELS_PER_PIXEL * (j * width + i);
			//index points to the Red channel of pixel 
			//at column i and row j calculated from top left
			return index;
		};
	};
	var pixelIndexer = makePixelIndexer(imageWidth);
	var drawPixel = function (x, y) {
		var idx = pixelIndexer(x,y);
		pixelData[idx] = 255;	//red
		pixelData[idx + 1] = 0;	//green
		pixelData[idx + 2] = 255;//blue
		pixelData[idx + 3] = 255;//alpha
	};

	while (x >= y) {
		drawPixel(x + x0, y + y0);
		drawPixel(y + x0, x + y0);
		drawPixel(-x + x0, y + y0);
		drawPixel(-y + x0, x + y0);
		drawPixel(-x + x0, -y + y0);
		drawPixel(-y + x0, -x + y0);
		drawPixel(x + x0, -y + y0);
		drawPixel(y + x0, -x + y0);
		if (decisionOver2 <= 0) 
		{
		    y++;
			decisionOver2 += dy; // Change in decision criterion for y -> y+1
			dy += 2;
		} 
		if (decisionOver2 > 0)
		{
			x--;
			dx += 2;
			decisionOver2 += (-radius << 1) + dx; // Change for y -> y+1, x -> x-1
		}
	}

	context.putImageData(imageData, 0, 0);
}

Drawing incomplete octants

The implementations above always draw only complete octants or circles. To draw only a certain arc from an angle $\alpha$ to an angle $\beta$, the algorithm needs first to calculate the $x$ and $y$ coordinates of these end points, where it is necessary to resort to trigonometric or square root computations (see Methods of computing square roots). Then the Bresenham algorithm is run over the complete octant or circle and sets the pixels only if they fall into the wanted interval. After finishing this arc, the algorithm can be ended prematurely.

If the angles are given as slopes, then no trigonometry or square roots are necessary: simply check that $y/x$ is between the desired slopes.

Generalizations

It is also possible to use the same concept to rasterize a parabola, ellipse, or any other two-dimensional curve.

External links

• Drawing circles - An article on drawing circles, that derives from a simple scheme to an efficient one
• Midpoint Circle Algorithm in several programming languages