Bezier Curve Turning Points

2023-08-08

In graphics applications, it is mainly Cubic Bezier curves that are used. Cubic Bezier curves are defined by the following equation:

$$\mathbf{B}(t)=(1-t{)}^3{\mathbf{P}}_0+3(1-t{)}^{2}t{\mathbf{Q}}_0+3(1-t)t^2{\mathbf{Q}}_1+t^3{\mathbf{P}}_1$$

def bezier(p0, q0, q1, p1):
    return lambda t : (1-t)**3 * p0 + 3*(1-t)**2 * t * q0 \
                      + 3*(1-t) * t**2 * q1 + t**3 * p1

These are two dimensional curves that go from a point ${\mathbf{P}}_0$ to some other point ${\mathbf{P}}_1$. The points ${\mathbf{Q}}_0$ and ${\mathbf{Q}}_1$ are the control points, and they manipulate the shape of the curve.

When trying to code vector shapes such as Bezier curves, computing their bounding rectangle is important, as this allows the graphics computations to efficient understand where to begin evaluating the curve and its coloring from. If it didn't have a boundary rectangle, then each individual pixel on the screen would have to be checked to see it's part of the curve or not. Furthermore, these bounding rectangle computations are also used to determine the selection rectangle for when the curve is selected, or if its bounds are within the bounds of the selection attempted by the cursor.

To compute the bounding rectangle, we find the turning points of the curve. These are the points where the derivative of the curve is zero. The derivative of the cubic bezier is:

$$\begin{array}{rl}{\mathbf{B}}^{\prime }(t)& =-3(1-t{)}^2{\mathbf{P}}_0+3(-2(1-t)t+(1-t{)}^2){\mathbf{Q}}_0+3(-t^2+2t(1-t)){\mathbf{Q}}_1+3t^2{\mathbf{P}}_1\\ & =-3(1-t{)}^2{\mathbf{P}}_0+3(-2t+t^2+t^2-2t+1){\mathbf{Q}}_0+3(-t^2+2t-2t^2){\mathbf{Q}}_1+3t^2{\mathbf{P}}_1\\ & =(-3t^2+6t-3){\mathbf{P}}_0+(6t^2-12t+3){\mathbf{Q}}_0+(-9t^2+6t){\mathbf{Q}}_1+3t^2{\mathbf{P}}_1\\ & =(-3{\mathbf{P}}_0+6{\mathbf{Q}}_0-9{\mathbf{Q}}_1+3{\mathbf{P}}_1)t^2+(6{\mathbf{P}}_0-12{\mathbf{Q}}_0+6{\mathbf{Q}}_1)t+(-3{\mathbf{P}}_0+3{\mathbf{Q}}_0)\\ & =3[(-{\mathbf{P}}_0+2{\mathbf{Q}}_0-3{\mathbf{Q}}_1+{\mathbf{P}}_1)t^2+2({\mathbf{P}}_0-2{\mathbf{Q}}_0+{\mathbf{Q}}_1)t+(-{\mathbf{P}}_0+{\mathbf{Q}}_0)]\\ & =3[\mathbf{A}{t}^2+2\mathbf{B}t+\mathbf{C}]\end{array}$$

where $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$ are the vector coefficients simplfied from their respected brackets:

$$\begin{array}{rl}\mathbf{A}& =-{\mathbf{P}}_0+2{\mathbf{Q}}_0-3{\mathbf{Q}}_1+{\mathbf{P}}_1\\ \mathbf{B}& ={\mathbf{P}}_0-2{\mathbf{Q}}_0+{\mathbf{Q}}_1\\ \mathbf{C}& =-{\mathbf{P}}_0+{\mathbf{Q}}_0\end{array}$$

From this simplification, we have a quadratic formula which will be computed through element-wise vector operations. That is, products of two vectors will be defined as $(a_x,a_y)\cdot (b_x,b_y)=(a_x{b}_x,a_y{b}_y)$, as opposed to either the dot product or cross product. This is because this definition is simply more convenient for our derivation. We similarly define an element-wise division of vectors, for the duration of this derivation: $(a_x,a_y)/(b_x,b_y)=(a_x/b_x,a_y/b_y)$.

By the quadratic formula, we want that ${\mathbf{B}}^{\prime }(t)=0$, so we have:

$$\begin{array}{rl}0& =3[\mathbf{A}{t}^2+2\mathbf{B}t+\mathbf{C}]\\ 0& =\mathbf{A}{t}^2+2\mathbf{B}t+\mathbf{C}\\ 0& =t^2+2\frac{\mathbf{B}}{\mathbf{A}}t+\frac{\mathbf{C}}{\mathbf{A}}=t^2+2{\mathbf{B}}^{\ast }t+{\mathbf{C}}^{\ast }\end{array}$$

where ${\mathbf{B}}^{\ast }=\mathbf{B}/\mathbf{A}$ and ${\mathbf{C}}^{\ast }=\mathbf{C}/\mathbf{A}$. We'll complete the square to simplify:

$$\begin{array}{rl}0& =(t+{\mathbf{B}}^{\ast }{)}^2-({\mathbf{B}}^{\ast }{)}^2+{\mathbf{C}}^{\ast }\\ (t+{\mathbf{B}}^{\ast }{)}^2& =({\mathbf{B}}^{\ast }{)}^2-{\mathbf{C}}^{\ast }\\ t& ={\mathbf{B}}^{\ast }\pm \sqrt{({\mathbf{B}}^{\ast }{)}^2-{\mathbf{C}}^{\ast }}=\frac{\mathbf{B}\pm \sqrt{{\mathbf{B}}^2-\mathbf{AC}}}{\mathbf{A}}\end{array}$$

Let's clarify what it is we've done this for. This above derivation will give us 4 separate equations for $t$, two for the horizontal turning points, and two for the vertical turning points. We then compare these turning points with ${\mathbf{P}}_0$ and ${\mathbf{P}}_1$ to find the minimum and maximum values for the curve. The turning points, can be paired up for the $x$ and $y$ coordinates to give us vector equations.

Since the $t$ values will essentially give us the turning points of the Bezier curve, one pair for each coordinate, let us rewrite it in terms of a vector $\mathbf{T}$:

$$\mathbf{T}=\frac{\mathbf{B}\pm \sqrt{{\mathbf{B}}^2-\mathbf{AC}}}{\mathbf{A}}$$

Not much of a change, besides the fact we're now looking at system of two equations which we've condensed into vector form. Again, keep in mind all arithmetic operations here are just element-wise vector operations on numbers: the subtractions, additions, square roots, products, and divisions.

The code from here is straightforward. We apply the above values and compute $\mathbf{T}$ directly in code. Further simplification doesn't help us much. Though I've left the below derivation for those who are interested.

$$\begin{array}{rl}{\mathbf{B}}^2& =({\mathbf{P}}_0-2{\mathbf{Q}}_0+{\mathbf{Q}}_1{)}^2=({\mathbf{P}}_0-2{\mathbf{Q}}_0+{\mathbf{Q}}_1)({\mathbf{P}}_0-2{\mathbf{Q}}_0+{\mathbf{Q}}_1)\\ & =({\mathbf{P}}_0^2-2{\mathbf{P}}_0{\mathbf{Q}}_0+{\mathbf{P}}_0{\mathbf{Q}}_1)-2({\mathbf{P}}_0{\mathbf{Q}}_0-2{\mathbf{Q}}_0^2+{\mathbf{Q}}_0{\mathbf{Q}}_1)+({\mathbf{Q}}_1{\mathbf{P}}_0-2{\mathbf{Q}}_0{\mathbf{Q}}_1+{\mathbf{Q}}_1^2)\\ & =({\mathbf{P}}_0^2+4{\mathbf{Q}}_0^2+{\mathbf{Q}}_1^2)+2(-2{\mathbf{P}}_0{\mathbf{Q}}_0+{\mathbf{P}}_0{\mathbf{Q}}_1-2{\mathbf{Q}}_0{\mathbf{Q}}_1)\\ \mathbf{A}\mathbf{C}& =(-{\mathbf{P}}_0+2{\mathbf{Q}}_0-3{\mathbf{Q}}_1+{\mathbf{P}}_1)(-{\mathbf{P}}_0+{\mathbf{Q}}_0)\\ & ={\mathbf{P}}_0^2-2{\mathbf{P}}_0{\mathbf{Q}}_0+3{\mathbf{P}}_0{\mathbf{Q}}_1-{\mathbf{P}}_0{\mathbf{P}}_1-{\mathbf{P}}_0{\mathbf{Q}}_0+2{\mathbf{Q}}_0^2-3{\mathbf{Q}}_0{\mathbf{Q}}_1+{\mathbf{Q}}_0{\mathbf{P}}_1\\ & =({\mathbf{P}}_0^2+2{\mathbf{Q}}_0^2)+(-3{\mathbf{P}}_0{\mathbf{Q}}_0+3{\mathbf{P}}_0{\mathbf{Q}}_1-{\mathbf{P}}_0{\mathbf{P}}_1-3{\mathbf{Q}}_0{\mathbf{Q}}_1+{\mathbf{Q}}_0{\mathbf{P}}_1)\\ & =({\mathbf{P}}_0^2+2{\mathbf{Q}}_0^2)+3(-{\mathbf{P}}_0{\mathbf{Q}}_0+{\mathbf{P}}_0{\mathbf{Q}}_1-{\mathbf{Q}}_0{\mathbf{Q}}_1)-({\mathbf{P}}_0{\mathbf{P}}_1+{\mathbf{Q}}_0{\mathbf{P}}_1)\\ {\mathbf{B}}^2-\mathbf{A}\mathbf{C}& =(2{\mathbf{Q}}_0^2+{\mathbf{Q}}_1^2)+(-7{\mathbf{P}}_0{\mathbf{Q}}_0+5{\mathbf{P}}_0{\mathbf{Q}}_1-7{\mathbf{Q}}_0{\mathbf{Q}}_1)+({\mathbf{P}}_0{\mathbf{P}}_1+{\mathbf{Q}}_0{\mathbf{P}}_1)\end{array}$$

Due to the square roots and the above not simplifying into a nice square, the computations will only get uglier if we try to simplify further. So we'll just apply the values for $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$ directly in code.

def turning_points_t(p0, q0, q1, p1):
    a = - p0 + 2 * q0 - 3 * q1 + p1
    b = p0 - 2 * q0 + q1
    c = - p0 + q0
    d = sqrt(b * b - a * c)
    return ((b - d) / a, (b + d) / a);

Again, as a reminder, the above code gets applied twice: once for the x-coordinate, and once for the y-coordinate. It however, only provides the $t$-values for the turning points. These turning points need to then be compared with the start and end points: ${\mathbf{P}}_0$ and ${\mathbf{P}}_1$ respectively.

Furthermore, these values must be between 0 and 1. If they are not, then they are not part of the drawn Bezier curve, and are therefore ignored. Thus, we need to compare the following 4 points, and get the minimum and maximum in both x- and y-coordinates:

$$\{{\mathbf{P}}_0,\mathbf{B}(\mathbf{T}[0]),\mathbf{B}(\mathbf{T}[1]),{\mathbf{P}}_1\}$$

We can however, ignore if they are below zero or above one. We will work around this as shown:

def bounds(p0, r, s, p1):
    # r and s are bezier(t), where t are the turning points
    lo = min(p0, r, s, p1)
    hi = max(p0, r, s, p1)

    return (lo, hi)

And yet again, the above code applies once to the x-coordinate, and once to the y-coordinate.

Putting it all together, we get the following:

def bezier_bounds(p0, q0, q1, q1):
    bezier = lambda t: (1 - t) ** 3 * p0 + 3 * (1 - t) ** 2 * t * q0 \
                        + 3 * (1 - t) * t ** 2 * q1 + t ** 3 * p1

    def turning_points_t(p_0, q_0, q_1, p_1):
        a = - p_0 + 2 * q_0 - 3 * q_1 + p_1
        b = p_0 - 2 * q_0 + q_1
        c = - p_0 + q_0
        d = sqrt(b * b - a * c)
        return ((b - d) / a, (b + d) / a);

    t_x = turning_points_t(p0.x, q0.x, q1.x, p1.x)
    t_y = turning_points_t(p0.y, q0.y, q1.y, p1.y)

    r_x, s_x = bezier(t_x[0]).x, bezier(t_x[1]).x
    r_y, s_y = bezier(t_y[0]).y, bezier(t_y[1]).y

    x_min = max(min(p0.x, r_x, s_x, p1.x), min(p0.x, p1.x))
    x_max = min(max(p0.x, r_x, s_x, p1.x), max(p0.x, p1.x))
    y_min = max(min(p0.y, r_y, s_y, p1.y), min(p0.y, p1.y))
    y_max = min(max(p0.y, r_y, s_y, p1.y), max(p0.y, p1.y))

    return (x_min, x_max, y_min, y_max)

Our workaround for when $t$ is less than 0 or greater than 1 is to simply use the min/max values of the start and end points (${\mathbf{P}}_0$ and ${\mathbf{P}}_1$). This is because the start and end points are the limit points of the curve. Either the remainder of the curve will be between them, or exceed their bounds. If it exceeds their bounds, then bound will be at one of the turning points. Otherwise, the bound will be at the start or end point. This is the purpose of the second min/max computation at the end, in the above code.