Midpoint circle algorithm and decision parameter

The midpoint circle algorithm draws a circle using only integer additions. Only one octant (from 90° to 45°) is computed and the rest is obtained by symmetry. The algorithm decides between two candidate pixels at each step by testing the midpoint between them against the circle equation.

Derivation of the decision parameter

Take a circle centred at the origin with radius r:

f(x, y) = x² + y² − r²
f(x, y) < 0  : point is inside the circle
f(x, y) = 0  : point is on the circle
f(x, y) > 0  : point is outside the circle

Start at (0, r) and move in the +x direction (clockwise) in the octant where the slope has magnitude less than 1 (x increases by 1 each step). If the pixel (x_k, y_k) has been plotted, the next pixel is either (x_k+1, y_k) or (x_k+1, y_k−1). The midpoint of these two candidates is (x_k+1, y_k−½). The decision parameter is the circle function at this midpoint:

p_k = f(x_k+1, y_k−½)
    = (x_k+1)² + (y_k−½)² − r²

• If p_k < 0, the midpoint is inside the circle, so the pixel at y_k is closer to the circle: choose (x_k+1, y_k).
• If p_k ≥ 0, the midpoint is outside or on the circle, so the pixel at y_k−1 is closer: choose (x_k+1, y_k−1).

Incremental form. The next parameter is

p_k+1 = f(x_k+1 + 1, y_k+1 − ½)
      = [(x_k+1)+1]² + (y_k+1 − ½)² − r²

Subtracting p_k:

p_k+1 = p_k + 2(x_k+1) + (y²_k+1 − y²_k) − (y_k+1 − y_k) + 1

• Case p_k < 0 (y_k+1 = y_k): p_k+1 = p_k + 2x_k + 3
• Case p_k ≥ 0 (y_k+1 = y_k − 1): p_k+1 = p_k + 2x_k + 3 − 2y_k + 2 = p_k + 2(x_k − y_k) + 5

Initial value. At the start point (0, r) the first midpoint is (1, r−½):

p_0 = f(1, r−½) = 1 + (r−½)² − r² = 5/4 − r

For integer radius, 5/4 is rounded to 1 (all other terms are integers and p is only compared with 0), so p_0 = 1 − r.

Midpoint circle algorithm (centre (xc, yc), radius r)

1. Input r and centre (xc, yc).
2. x = 0, y = r.
3. p = 1 − r.
4. Plot the 8 symmetric points of (x, y) shifted by (xc, yc):
   (xc±x, yc±y) and (xc±y, yc±x).
5. While x < y:
      x = x + 1
      if p < 0:
          p = p + 2x + 1          // 2x_k + 3 in terms of old x
      else:
          y = y − 1
          p = p + 2(x − y) + 1    // 2(x_k − y_k) + 5 in terms of old values
      plot the 8 symmetric points of (x, y)
6. Stop.

(Here the updates use the new x, y: with x_new = x_k+1, 2x_k+3 = 2x_new+1 and 2(x_k−y_k)+5 = 2(x_new − y_new)+1.)

Example, r = 10

p_0 = 1 − 10 = −9. The octant points are (0,10), (1,10), (2,10), (3,10), (4,9), (5,9), (6,8), (7,7).