The Bézier circle

How close does a Bézier curve come to a circle? Fonts, SVG paths and most drawing programs draw a circle as four cubic Bézier curves, one per quarter. The two inner control points of each quarter sit on the tangents at its ends, a distance \(k r\) out. A cubic polynomial cannot follow a circle exactly. The usual choice of handle,

\[k = \tfrac{4}{3}\left(\sqrt{2} - 1\right) \approx 0.5523,\]

puts the middle of each quarter on the circle, and everywhere else the curve runs outside it, by at most 0.027 % of the radius. CAD systems store circles as rational B-splines instead, and a rational quadratic quarter is exactly the circle.

The example draws the four cubic quarters with px.BSpline, measures how far they run from the circle, finds the handle with the smallest largest deviation and checks the rational quarter.

This page is part of a collection of geometry and design examples using planex2d.

Open this example in the playground, the interactive code editor in the browser, to change the code and run it.

import math

import numpy as np
from scipy.optimize import minimize_scalar
import planex as px

Four cubic quarters

With the radius \(r = 10\), the quarter from \((r, 0)\) to \((0, r)\) has the control points \((r, 0)\), \((r, k r)\), \((k r, r)\) and \((0, r)\). A cubic Bézier curve is a B-spline of degree 3 with the knots 0, 0, 0, 0, 1, 1, 1, 1. The other three quarters are the first one turned by right angles, and together they close the circle. planex2d computes the area inside them from the exact curves (Figure 1):

r = 10
k = 4 / 3 * (math.sqrt(2) - 1)


def quarter(k):
    """the control points of the cubic quarter from (r, 0) to (0, r) with handles k r long"""
    return np.array([(r, 0), (r, k * r), (k * r, r), (0, r)])


def turned(points, j):
    """the points turned j right angles about the center"""
    c, s = math.cos(j * math.pi / 2), math.sin(j * math.pi / 2)
    return [(c * x - s * y, s * x + c * y) for x, y in points]


def circle(k):
    """the face inside four cubic quarters"""
    return px.Face([px.BSpline(turned(quarter(k), j), [0, 0, 0, 0, 1, 1, 1, 1], 3) for j in range(4)])


bezier = circle(k)
print("k", k, "| area", bezier.mass, "| π r²", math.pi * r**2, "| relative", bezier.mass / (math.pi * r**2) - 1)
# px.Plot(bezier)

k 0.5522847498307935 | area 314.2472332656507 | π r² 314.1592653589793 | relative 0.00028001054360404787

Figure 1: The Bézier circle, its control points and handles (dashed)

How far from the circle

A point of the curve is the sum of its control points weighted by the Bernstein polynomials \((1 - t)^3\), \(3t(1 - t)^2\), \(3t^2(1 - t)\) and \(t^3\), so its distance from the center is known at every \(t\). At the middle, \(t = 1/2\), the curve is on the circle. Everywhere else it is outside, most at \(t \approx 0.21\) and \(0.79\) (Figure 2):

t = np.linspace(0, 1, 20001)
bernstein = np.array([(1 - t)**3, 3 * t * (1 - t)**2, 3 * t**2 * (1 - t), t**3]).T


def deviation(k):
    """the distance from the center along the quarter, relative to r, less 1"""
    return np.hypot(*(bernstein @ quarter(k)).T) / r - 1


d = deviation(k)
print("middle", d[10000], "| largest", d.max(), "at t", t[d.argmax()], "| smallest", d.min())

middle 0.0 | largest 0.00027253000136240324 at t 0.7887000000000001 | smallest 0.0

A better handle

The largest deviation is smallest when the curve runs as far inside the circle as outside. SciPy’s minimize_scalar finds that handle. It brings the largest deviation down by a quarter and the area error by more than two thirds (Figure 2):

best = minimize_scalar(lambda k: abs(deviation(k)).max(), bounds=(0.54, 0.56), method="bounded",
                       options={"xatol": 1e-12})
d_best = deviation(best.x)
print("k", best.x, "| largest", d_best.max(), "| smallest", d_best.min())
print("area relative", circle(best.x).mass / (math.pi * r**2) - 1)
# plt.plot(t, 100 * d, t, 100 * d_best); plt.show()

k 0.5519150250750289 | largest 0.00019607657118991817 | smallest -0.00019607616148020757
area relative 7.553164089002529e-05

Figure 2: How far a cubic quarter runs from the circle

The exact quarter

A rational B-spline divides by weights. The quadratic quarter with the control points \((r, 0)\), \((r, r)\), \((0, r)\) and the weights \(1\), \(\sqrt{2}/2\), \(1\) is exactly the circle. planex2d measures its length and the area inside four of them from the exact curves:

w = math.sqrt(2) / 2
exact = [px.BSpline(turned([(r, 0), (r, r), (0, r)], j), [0, 0, 0, 1, 1, 1], 2, weights=[1, w, 1]) for j in range(4)]
print("quarter length", exact[0].mass, "| π r / 2", math.pi * r / 2)
print("area", px.Face(exact).mass, "| π r²", math.pi * r**2)

quarter length 15.707963267948966 | π r / 2 15.707963267948966
area 314.15926535897927 | π r² 314.1592653589793

See Splines for px.BSpline and rational splines.