In computer graphics, the midpoint circle algorithm is an algorithm used to determine the points needed for rasterizing a circle. It's a generalization of Bresenham's line algorithm. The algorithm can be further generalized to conic sections.[1] [2] [3]
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 =, it has reached 45°. The reason for using these angles is shown in the above picture: As increases, it neither skips nor repeats any value until reaching 45°. So during the while loop, increments by 1 with each iteration, and decrements by 1 on occasion, never exceeding 1 in one iteration. This changes at 45° because that is the point where the tangent is =. Whereas before and after.
The second part of the problem, the determinant, is far trickier. This determines when to decrement . 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 (or whatever the third variable is), it stands to reason that the algorithm for a discrete (voxel) sphere would also rely on the 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 objective of the algorithm is to approximate a circle, more formally put, to approximate the curve
x2+y2=r2
x2+y2\leqr2
x2+y2
This algorithm starts with the circle equation. For simplicity, assume the center of the circle is at
(0,0)
(r,0)
The fast direction here (the basis vector with the greater increase in value) is the
y
y
x
From the circle equation is obtained the transformed equation
x2+y2-r2=0
r2
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
(r,0)
For each point, the following holds:
\begin{align}
2 | |
x | |
n |
+
2 | |
y | |
n |
=r2\end{align}
\begin{align}
2 | |
x | |
n |
=r2-
2 | |
y | |
n |
\end{align}
\begin{align}
2 | |
x | |
n+1 |
=r2-
2 | |
y | |
n+1 |
\end{align}
Since for the first octant the next point will always be at least 1 pixel higher than the last (but also at most 1 pixel higher to maintain continuity), it is true that:
\begin{align}
2 | |
y | |
n+1 |
&=(yn+1)2\ &=
2 | |
y | |
n |
+2yn+1\end{align}
\begin{align}
2 | |
x | |
n+1 |
=r2-
2 | |
y | |
n |
-2yn-1\end{align}
2 | |
x | |
n |
=r2-
2 | |
y | |
n |
\begin{align}
2 | |
x | |
n+1 |
=
2 | |
x | |
n |
-2yn-1\end{align}
Because of the continuity of a circle and because the maxima along both axes are the same, clearly it will not be skipping points as it advances in the sequence. Usually it stays on the same coordinate, and sometimes advances by one to the left.
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
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.
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
(xi,yi)
RE(xi,yi)=\left\vert
2 | |
x | |
i |
+
2 | |
y | |
i |
-r2\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)
RE(xi,yi)=\left\vert
2 | |
x | |
i |
+02-r2\right\vert=0
Because it starts in the first counter-clockwise positive octant, it will step in the direction with the greatest travel, the Y direction, so it is clear that
yi+1=yi+1
RE(xi-1,yi+1)<RE(xi,yi+1)
If this inequality holds, then plot
(xi-1,yi+1)
(xi,yi+1)
\begin{align} RE(xi-1,yi+1)&<RE(xi,yi+1)\\ \left\vert
2 | |
(x | |
i-1) |
+
2 | |
(y | |
i+1) |
-r2\right\vert&<\left\vert
2 | |
x | |
i |
+
2 | |
(y | |
i+1) |
-r2\right\vert\\ \left\vert
2 | |
(x | |
i |
-2xi+1)+
2 | |
(y | |
i |
+2yi+1)-r2\right\vert&<\left\vert
2 | |
x | |
i |
+
2 | |
(y | |
i |
+2yi+1)-r2\right\vert\\ \end{align}
The absolute value function does not help, so square both sides, since a square is always positive:
\begin{align} \left[
2 | |
(x | |
i |
-2xi+1)+
2 | |
(y | |
i |
+2yi+1)-r2\right]2&<\left[
2 | |
x | |
i |
+
2 | |
(y | |
i |
+2yi+1)-r2\right]2\\ \left[
2 | |
(x | |
i |
+
2 | |
y | |
i |
-r2+2yi+1)+(1-2xi)\right]2&<\left[
2 | |
x | |
i |
+
2 | |
y | |
i |
-r2+2yi+1\right]2\\ \left(
2 | |
x | |
i |
+
2 | |
y | |
i |
-r2+2yi+1\right)2+2(1-2xi)
2 | |
(x | |
i |
+
2 | |
y | |
i |
-r2+2yi+1)+(1-2
2 | |
x | |
i) |
&<\left[
2 | |
x | |
i |
+
2 | |
y | |
i |
-r2+2yi+1\right]2\\ 2(1-2xi)
2 | |
(x | |
i |
+
2 | |
y | |
i |
-r2+2yi+1)+(1-2
2 | |
x | |
i) |
&<0\\ \end{align}
Since > 0, the term
(1-2xi)<0
\begin{align} 2\left[
2 | |
(x | |
i |
+
2 | |
y | |
i |
-r2)+(2yi+1)\right]+(1-2xi)&>0\\ 2\left[RE(xi,yi)+YChange\right]+XChange&>0\\ \end{align}
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, then decrement the value. If, then keep the same value. Again, by reflecting these points in all the octants, a full circle results.
We may reduce computation by only calculating the delta between the values of this decision formula from its value at the previous step. We start by assigning as which is the initial value of the formula at, then as above at each step if we update it as (and decrement), otherwise thence increment as usual.
The algorithm has already been explained to a large extent, but there are further optimizations.
The new presented method[4] gets along with only 5 arithmetic operations per step (for 8 pixels) and is thus best suitable for low-performate systems. In the "if" operation, only the sign is checked (positive? Yes or No) and there is a variable assignment, which is also not considered an arithmetic operation.The initialization in the first line (shifting by 4 bits to the right) is only due to beauty and not really necessary.
So we get countable operations within main-loop:
The comparison t2 >= 0 is not counted as no real arithmetic takes place. In two's complement representation of the vars only the sign bit has to be checked.
Operations: 5
t1 = r / 16 x = r y = 0 Repeat Until x < y Pixel (x, y) and all symmetric pixels are colored (8 times) y = y + 1 t1 = t1 + y t2 = t1 - x If t2 >= 0 t1 = t2 x = x - 1
The implementations above always draw only complete octants or circles. To draw only a certain arc from an angle
\alpha
\beta
x
y
If the angles are given as slopes, then no trigonometry or square roots are necessary: simply check that
y/x
It is also possible to use the same concept to rasterize a parabola, ellipse, or any other two-dimensional curve.[5]