The objective of the algorithm is to approximate a
circle, more formally put, to approximate the curve x^2 + y^2 = r^2 using pixels; in
layman's terms every pixel should be approximately the same distance from the center, as is the definition of a circle. 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. But a simplification can be done in order to understand the bitshift. Keep in mind that a left bitshift of a
binary number is the same as multiplying with 2. Ergo, a left bitshift of the radius only produces the
diameter which is defined as radius times two. This algorithm starts with the
circle equation. For simplicity, assume the center of the circle is at (0,0). To start with, consider 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 (see
Differentiation of trigonometric functions). 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{align} x_n^2 + y_n^2 = r^2 \end{align} This can be rearranged thus: :\begin{align} x_n^2 = r^2 - y_n^2 \end{align} And likewise for the next point: :\begin{align} x_{n+1}^2 = r^2 - y_{n+1}^2 \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} y_{n+1}^2 &= (y_n + 1)^2 \\ &= y_n^2 + 2y_n + 1 \end{align} :\begin{align} x_{n+1}^2 = r^2 - y_n^2 - 2y_n - 1 \end{align} So, rework the next-point-equation into a recursive one by substituting x_n^2 = r^2 - y_n^2: :\begin{align} x_{n+1}^2 = x_n^2 - 2y_n - 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 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 is 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 x_i^2 + 0^2 - r^2 \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 y_{i+1} = y_i + 1. Also, because it concerns this octant only, the 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) 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{align} RE(x_i-1, y_i+1) & The
absolute value function does not help, so square both sides, since a square is always positive: : \begin{align} \left [ (x_i^2 - 2 x_i + 1) + (y_i^2 + 2 y_i + 1) - r^2 \right ]^2 & Since > 0, the term (1 - 2 x_i) , so dividing gets: : \begin{align} 2 \left [ (x_i^2 + y_i^2 - r^2) + (2 y_i + 1) \right ] + (1 - 2 x_i) & > 0 \\ 2 \left [ RE(x_i,y_i) + Y_\text{Change} \right ] + X_\text{Change} & > 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 {{tmath|2(RE+Y_\text{Change})+X_\text{Change} > 0}}, then decrement the value. If {{tmath|2(RE+Y_\text{Change})+X_\text{Change} \le 0}}, 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.
Jesko's method The algorithm has already been explained to a large extent, but there are further optimizations. The new presented method gets along with only 5 arithmetic operations per step (for 8 pixels) and is thus best suitable for low-performance 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 x >= y (is counted as a subtraction: x - y >= 0) • y=y+1 [y++] • t1 + y • t1 - x • : 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. • x=x-1 [x--] Operations: 5
t1 = r / 16
x = r
y = 0
Repeat Until x = 0
t1 =
t2 x =
x - 1 == Drawing incomplete octants ==