Midpoint circle algorithm 📖 Wikipedia

In computer graphics, the midpoint circle algorithm is an algorithm used to determine the points needed for rasterizing a circle. It is 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 y = x, it has reached 45°. The reason for using these angles is shown in the above picture: As x increases, it neither skips nor repeats any x value until reaching 45°. So during the while loop, x increments by 1 with each iteration, and y 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 y. 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 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.

• One hundred fifty concentric circles drawn with the midpoint circle algorithm.
• 
  On left, all circles are drawn black.
• 
  On right, red, black and blue are used together to demonstrate the concentricity of the circles.

The objective of the algorithm is to approximate a circle or, more formally, the curve ${\displaystyle x^{2}+y^{2}=r^{2}}$ using pixels. For simplicity, we aim to approximate a circle of center ${\displaystyle (0,0)}$ and integer radius ${\displaystyle r\in \mathbb {N} }$. Furthermore, we only draw to the first octant of the plane, for which we draw a curve from the point ${\displaystyle (r,0)}$ and proceed counterclockwise until an angle of 45°; the remaining octants can be trivially filled by rotating and/or mirroring that first segment of the curve around the center. At each step, the path is extended by choosing the adjacent pixel which best satisfies ${\displaystyle x^{2}+y^{2}=r^{2}}$.

The "fast" direction in the first octant (the basis vector with the greater increase in value) is the ${\displaystyle y}$-direction (see differentiation of trigonometric functions). In practice, this means that the algorithm always takes a step in the positive ${\displaystyle y}$ direction (upwards), and occasionally takes a step in the "slow" direction (the negative ${\displaystyle x}$ direction, leftwards). If the point drawn at step ${\displaystyle n}$ is ${\displaystyle (x_{n},y_{n})}$, we thus impose ${\displaystyle y_{n+1}=y_{n}+1}$ and have to decide whether to set ${\displaystyle x_{n+1}}$ to ${\displaystyle x_{n}}$ or ${\displaystyle x_{n}-1}$.

The choice is made by computing the distance between the decision midpoint ${\displaystyle (x_{n}-0.5,y_{n}+1)}$ to the center of the circle, hence the name of the algorithm. With the above assumption of a circle of center ${\displaystyle (0,0)}$, this distance is equal to ${\displaystyle d={\sqrt {(x_{n}-0.5)^{2}+(y_{n}+1)^{2}}}}$. If ${\displaystyle d>r}$ the midpoint mathematically lies outside the circle, and the innermost candidate point at ${\displaystyle (x_{n}-1,y_{n}+1)}$ should be drawn. Inversely if ${\displaystyle d<r}$, the midpoint lies inside the circle and the candidate point ${\displaystyle (x_{n},y_{n}+1)}$ best satisfies the circle equation. The full iteration for the first octant is as follows:

${\displaystyle {\begin{aligned}x_{0}&=r\\y_{0}&=0\\d_{n}&={\sqrt {\left(x_{n}-{\frac {1}{2}}\right)^{2}+(y_{n}+1)^{2}}}\\x_{n+1}&={\begin{cases}x_{n}&{\text{if }}d_{n}<r\\x_{n}-1&{\text{if }}d_{n}>r\end{cases}}\\y_{n+1}&=y_{n}+1\end{aligned}}}$

The iteration ends once the 45° slope line ${\displaystyle y=x}$ is crossed, therefore as soon as the condition ${\displaystyle y>x}$ is met.

Just as with the parent Bresenham's line algorithm, this algorithm can be optimized to use only integer-based math. We start by defining the radius error ${\displaystyle P_{n}}$ as the difference between ${\displaystyle d_{n}^{2}}$ and the squared radius ${\displaystyle r^{2}}$ to avoid having to compute an expensive and non-integer square-root. This also slightly changes the condition for determining ${\displaystyle x_{n+1}}$:

${\displaystyle {\begin{aligned}P_{n}&=d_{n}^{2}-r^{2}=\left(x_{n}-{\frac {1}{2}}\right)^{2}+(y_{n}+1)^{2}-r^{2}\\x_{n+1}&={\begin{cases}x_{n}&{\text{if }}P_{n}<0\\x_{n}-1&{\text{if }}P_{n}>0\end{cases}}\end{aligned}}}$

To eliminate the ${\displaystyle 1/2}$ term, the radius error can be computed recursively:

${\displaystyle {\begin{aligned}P_{n+1}&=\left(x_{n+1}-{\frac {1}{2}}\right)^{2}+(y_{n+1}+1)^{2}-r^{2}\\&=\left(x_{n+1}-{\frac {1}{2}}\right)^{2}+(y_{n}+1)^{2}+2(y_{n}+1)+1-r^{2}\\&=\left(x_{n+1}-{\frac {1}{2}}\right)^{2}+P_{n}-\left(x_{n}-{\frac {1}{2}}\right)^{2}+2(y_{n}+1)+1\\&=P_{n}+\left(x_{n+1}^{2}-x_{n}^{2}\right)-(x_{n+1}-x_{n})+2(y_{n}+1)+1\end{aligned}}}$

Together with the condition on ${\displaystyle x_{n+1}}$, computing ${\displaystyle P_{n+1}}$ becomes

${\displaystyle P_{n+1}=P_{n}+{\begin{cases}2(y_{n}+1)+1&=2y_{n+1}+1&{\text{if }}P_{n}<0\\2(y_{n}+1)-2(x_{n}-1)+1&=2y_{n+1}-2x_{n+1}+1&{\text{if }}P_{n}>0\end{cases}}}$

An initial value ${\displaystyle P_{0}}$ for ${\displaystyle (x_{0},y_{0})=(r,0)}$ is found by approximation:

${\displaystyle {\begin{aligned}P_{0}&=\left(r-{\frac {1}{2}}\right)^{2}+(0+1)^{2}-r^{2}\\&=r^{2}-r+{\frac {1}{4}}+1-r^{2}\\&=1.25-r\\&\approx 1-r\end{aligned}}}$

This optimized algorithm can be implemented in Python as follows:

import numpy as np										# only for the image array

r = 67													# radius of the circle

img = np.zeros([2 * r + 1]*2, dtype=int)				# image of size (2r + 1) x (2r + 1) to fit the circle
x, y, p = r, 0, 1 - r									# initial values x0 = r, y0 = 0, p0 = 1 - r
while x >= y:											# while the point (x, y) is in the first octant
	for j, k in [(1, 1), (1, -1), (-1, 1), (-1, -1)]:	# draw to all quadrants
		img[j * x + r, k * y + r] = 1
		img[k * y + r, j * x + r] = 1
	x, y = x - (1 if p > 0 else 0), y + 1				# update x and y according to the radius error p
	p += 1 - 2 * x + 2 * y if p > 0 else 1 + 2 * y		# update the radius error p with updated x and y

An improved midpoint circle algorithm^[4] only requires 5 arithmetic operations per step (for 8 pixels) and is thus best suitable for low-performance systems. The counted operations in the main loop are:

1. The comparison x >= y (is counted as a subtraction: x - y >= 0)
2. y=y+1 [y++]
3. t1 + y
4. t1 - x

   The comparison t2 >= 0 is not counted as no real arithmetic takes place. In two's complement representation of the variables, only the sign bit has to be compared.

5. x=x-1 [x--]

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 ${\displaystyle \alpha }$ to an angle ${\displaystyle \beta }$, the algorithm needs first to calculate the ${\displaystyle x}$ and ${\displaystyle 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 ${\displaystyle y/x}$ is between the desired slopes.

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

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