Chapter 4 Recursive Curves

There is a close relationship between trees and recursive algorithms. We saw in the previous chapter how recursion can be used to explore and draw various trees and tree-like shapes, and how the shape of those trees mirrors the way the algorithms operate. In the programs we developed for drawing the Pythagoras tree, for instance, every square of the tree corresponds to a single recursive function call and every branch to an entire cascade of nested function calls. And just as each branch is a smaller version of the whole tree, the algorithm for drawing each branch is a smaller instance of the algorithm for the whole problem.

All recursive algorithms produce such a tree-like hierarchy of problems and sub-problems. Conversely, being able to visualize computational problems as a hierarchy of sub-problems is the key to using recursion effectively. In this chapter we consider a few additional geometric problems whose recursive nature isn’t quite as obvious.

4.1 Hilbert Curves¶

Figure 4.1 shows a so-called Hilbert curve, a curve that starts in the lower-left corner of the square, ends in the lower-right corner, and visits every point on a \(64\times 64\) grid exactly once. It’s clear that the curve is highly symmetrical and consists of a handful of building blocks that are repeated in different orientations and arrangements, but if someone told you that the curve has a recursive structure, you would probably have a hard time working out the construction rules from the final shape.

The Hilbert curve in Fig. 4.1 is just one example from an infinitely large family of curves that are constructed by repeatedly subdividing a simple base shape. This base shape is shown in the left part of Fig. 4.2: it is shaped like an inverted U that visits the corners of a square in clockwise order, starting in the lower-left corner and ending in the lower-right corner. In the first subdivision step, this U-shape is replaced with a more complicated shape that consists of four smaller U-shapes that are connected by three line segments, as shown in the right part of the figure. The two U-shapes at the bottom of the resulting Y-shape turn counterclockwise and the two at the top turn clockwise. Notice that the new curve visits the four corners of a larger square in the same order as before and touches every point on a \(4\times 4\) grid.

The same process can be repeated as many times as desired. In each step, we replace all U-shapes with four appropriately rotated and mirrored Y-shapes. (The subdivision rule for counterclockwise U-shapes can be obtained by reversing the direction of all arrows in Fig. 4.2.) The curve we obtain after \(k\) steps is called a Hilbert curve of order \(k\).

The first few Hilbert curves, scaled down to a common size, are shown in the lower part of Fig. 4.2. All Hilbert curves start in the lower-left corner, end in the lower-right corner, and visit every point on a \(2^{k+1}\times 2^{k+1}\) grid inside the square. It’s easy to see that the curves are self-similar in the sense that every quadrant of a Hilbert curve is itself a Hilbert curve of the previous order. For example, the order-3 curve consists of four order-2 curves, each of which consists of four order-1 curves, each of which consists of four order-0 curves.

Figure 4.2 hints at another remarkable property of Hilbert curves: As the order increases, the curve gets longer and coils up ever more tightly inside the square. Mathematically, as the order tends towards infinity, the curve becomes infinitely long and converges toward the unit square. The Hilbert curve is therefore said to be a space-filling curve, a class of one-dimensional curves that converge toward two-dimensional shapes.

Drawing Hilbert Curves

At the lowest level, every Hilbert curve consists of a sequence of horizontal and vertical line segments of uniform length. We can therefore describe its shape by specifying the curve’s starting point and the list of directions that lead from the start to the end point. For example, the basic U-shape starts in the lower-left corner and can be drawn by going north, east, and south. Similarly, the first-order curve consists of 15 line segments and can be drawn by going

The first three directions on each line describe rotated U-shapes that are connected by additional line segments shown in italics. Notice that the connecting lines by themselves form a basic U-shape north-east-south.

How do we represent these directions in a computer program? One option is to define four integer constants called NORTH, EAST, SOUTH, WEST as shown in Listing 4.1. The constants are defined in clockwise order, which allows us to rotate a direction clockwise by incrementing its value (wrapping around to 0 when going above WEST), and counterclockwise by decrementing it (wrapping around to \(3\) when going below NORTH). Using the modulo operator, we can express both kinds of rotations as follows:

\begin{equation} (\Id {direction} + \Id {steps}) \bmod 4. \end{equation}

Depending on whether steps is positive or negative, this expression rotates direction one or more steps in either clockwise or counterclockwise order. This operation is implemented by the turn() method defined in the same listing; the floorMod() function is used instead of the modulo operator % to ensure that the expression wraps around correctly when steps is negative.

// A class for computing Hilbert curves.
public class HilbertCurve {
    // Predefined directions.
    static final int NORTH = 0;
    static final int EAST = 1;
    static final int SOUTH = 2;
    static final int WEST = 3;

    // Turn the given direction clockwise ('steps' positive)
    // or counterclockwise ('steps' negative).
    static int turn(int direction, int steps) {
        return Math.floorMod(direction + steps, 4);
    }

    // ...
}

For example, if direction is an integer between 0 and 3, we can construct a clockwise U-shape that starts out in this direction by calling turn() twice to produce a list of three directions:

List.of(direction, turn(direction, 1), turn(direction, 2));

Similarly, a U-shape that bends counterclockwise is produced by

List.of(direction, turn(direction, -1), turn(direction, -2));

Now let’s extend this idea to general Hilbert curves. As indicated in Fig. 4.2, a Hilbert curve of order \(k\) consists of four Hilbert curves of order \(k-1\) that are connected using three line segments. A clockwise Hilbert curve that starts in a given direction can therefore be described recursively as follows:

We can use the same process for counterclockwise curves, except that we have to rotate directions using offsets of the opposite sign (for instance, \(\Id {direction}+2\) becomes \(\Id {direction}-2\)) and that we have to interchange CCW and CW.

Listing 4.2 shows an implementation of this idea. The hilbertCurve() function computes a Hilbert curve of a given order whose orientation is described by a direction and a rotation that is either \(+1\) for a clockwise shape or \(-1\) for a counterclockwise shape. The function outputs the resulting curve as list of directions: Each of the four recursive calls to hilbertCurve() adds the directions that describe the nested Hilbert curve of lower order to result, and the three calls to add() produce the short line segments that connect the nested shapes. The two nested curves in the middle have the same starting direction and orientation as the main curve, whereas the two nested curves at either end are rotated and have the opposite orientation.

public static void hilbertCurve(int order, int direction,
        int rotation, List<Integer> result) {
    if (order >= 0) {
        hilbertCurve(order - 1, turn(direction, rotation),
                -rotation, result);
        result.add(direction);
        hilbertCurve(order - 1, direction, rotation, result);
        result.add(turn(direction, rotation));
        hilbertCurve(order - 1, direction, rotation, result);
        result.add(turn(direction, 2 * rotation));
        hilbertCurve(order - 1, turn(direction, -rotation),
                -rotation, result);
    }
}

For example, the Hilbert curve of order \(2\) can be computed by setting off north in clockwise direction and collecting the resulting steps in an array directions:

var directions = new ArrayList<Integer>();
hilbertCurve(2, NORTH, +1, directions);
System.out.println(directions);

The output of this program is

[1, 0, 3, 0, 0, 1, 2, 1, 0, 1, 2, 2, 3, 2, 1]

which reads east, north, west, north, north, east, south, east, north, east, south, south, west, south, east. As desired, this describes the order 2-curve shown in Fig. 4.2.

Exercises

[4.1]Write a program that prints the coordinates of all points on a Hilbert curve of a given order. For example, the points along the order-1 Hilbert curve might be printed as

(0,0) (1,0) (1,1) (0,1) (0,2) (0,3) (1,3) (1,2) (2,2)
(2,3) (3,3) (3,2) (3,1) (2,1) (2,0) (3,0)

4.2 Bézier Curves¶

In the previous section we studied a few fractal curves that are constructed using recursive subdivision. But fractals aren’t the only kinds of geometric object that consists of smaller versions of themselves. In this section we will meet another example: the Bézier curve. A simple example of an individual Bézier curve is shown in the left part of Fig. 4.3. The curve is specified using four control points: the first and last control points determine where the curve starts and ends, and the two inner points determine how it bends. In the figure, the four control points are labeled \(\vec {p}_0\) through \(\vec {p}_3\). The curve starts at \(\vec {p}_0\) and initially moves toward the first inner control point \(\vec {p}_1\), before it gradually turns toward the second inner control point \(\vec {p}_2\) and finally the end point \(\vec {p}_3\).

Bézier curves are ubiquitous in computer graphics — In fact, there are thousands of these curves right in front of your eyes, regardless of whether you are reading this book on a screen or on paper, because the shapes of all letters and symbols in modern digital typefaces are specified using Bézier curves. For example, the right part of ?? shows how the shape of the letter q can be described by a sequence of line and curve segments that connect the black nodes on the outline. Each of the curved segments is a Bézier curves with four control points, as indicated in the figure.

How do we get from the four control points of a Bézier curve to the actual curve? One way is based on the mathematical idea of a parametric curve, which is a function that maps a number \(t\), the “parameter”, to a point on the curve. In the case of Bézier curves, this function has the following form:

\begin{equation} \label {eq:bezier} \vec {B}(t) = (1-t)^3 \vec {p}_0 + 3t(1-t)^2\vec {p}_1 + 3t^2(1-t)\vec {p}_2 + t^3\vec {p}_3. \end{equation}

The four points \(\vec {p}_0\), \(\vec {p}_1\), \(\vec {p}_2\), and \(\vec {p}_3\) are the control points of the curve, and the parameter \(t\) is a real number that runs from 0 to 1. Smoothly changing \(t\) produces a sequence of points that trace out the Bézier curve, starting at the first control point \(\vec {B}(0)=\vec {p}_0\) and ending at \(\vec {B}(1)=\vec {p}_3\).

Let’s take a closer look at Eq. (4.2). For a particular value of the parameter \(t\), the point \(\vec {B}(t)\) is a weighted sum of the four control points that uses the following weights:

\begin{equation} \label {eq:bezier-weight} w_0=(1-t)^3,\qquad w_1=3t(1-t)^2,\qquad w_2=3t^2(1-t),\qquad w_3=t^3. \end{equation}

Each of these weights is a polynomial of the parameter \(t\); the curves we are discussing here are also known as cubic Bézier curves because the highest power of \(t\) that occurs in the weights is \(3\). As we move along the curve, the magnitude of the four weights changes, and so does the influence of the four control points. This is the mathematical reason why Bézier curves bend as if they were attracted by the two inner control points. It’s also easy to show that, for any value \(t\in [0,1]\), the four weights satisfy the following two conditions:

\begin{equation} \label {eq:bezierweightsum} 0\le w_0, w_1, w_2, w_3 \le 1,\qquad \text {and}\qquad w_0 + w_1 + w_2 + w_3 = 1. \end{equation}

The geometric interpretation of Eq. (4.4) is that every point on the Bézier curve (and therefore the entire curve) lies inside the region spanned by the four control points. Mathematicians call this region the convex hull of the control points.

We can now write a first program for evaluating and drawing Bézier curves. As shown in Listing 4.3, we represent a Bézier curve as a record that holds its four control points. Arbitrary points on the curve can be computed using pointAt(), which evaluates Eq. (4.2) for a given value of the parameter \(t\) and returns the result as a Point.

// A Bézier curve defined by its four control points.
public record BezierCurve(
        Point p0, Point p1, Point p2, Point p3) {
    // Evaluate the Bézier curve at the given parameter.
    public Point pointAt(double t) {
        double w0 = Math.pow(1 - t, 3);
        double w1 = 3 * t * Math.pow(1 - t, 2);
        double w2 = 3 * Math.pow(t, 2) * (1 - t);
        double w3 = Math.pow(t, 3);
        return new Point(
                w0 * p0.x() + w1 * p1.x() + w2 * p2.x()
                        + w3 * p3.x(),
                w0 * p0.y() + w1 * p1.y() + w2 * p2.y()
                        + w3 * p3.y());
    }

    // ...
}

Drawing a Bézier curve is just as easy: We first construct a sequence of points along the curve and then connect adjacent points with line segments. For example, we can choose a sequence of \(N+1\) parameters \(t_0,t_1,\ldots ,t_{N}\) between 0 and 1, compute the corresponding \(N+1\) points on the curve

\begin{equation*} \vec {B}(t_0), \vec {B}(t_1), \dots , \vec {B}(t_{N}), \end{equation*}

and then draw a line segment between each pair of points \(\vec {B}(t_k)\), \(\vec {B}(t_{k+1})\). How do we choose the parameters \(t_k\)? There are a few obvious constraints. For instance, we want \(t_0=0\) and \(t_N=1\) so that the first line segment starts at \(\vec {p}_0\) and the last one ends at \(\vec {p}_3\). In addition, since we want the sequence of points to advance along the curve, the parameters must increase monotonically, that is, we must have \(t_k < t_{k+1}\).

The simplest way to choose the parameters \(t_k\) is to divide the parameter range \([0,1]\) into \(N\) intervals of the same size, which gives us

\begin{equation*} t_k = \frac {k}{N}. \end{equation*}

Figure 4.4 shows the result of splitting a Bézier curve into ten segments using this approach. As you can see, the line segments approximate the overall shape of the Bézier curve reasonably well; if desired, we can obtain better approximations by increasing the number of line segments \(N\).

The main problem with uniformly spaced parameters is their inflexibility. Even though we can increase the number of intervals \(N\) to produce more points on the Bézier curve, we have no direct control over where those points are placed. For example, it’s not possible to use more points in curved regions than in flat regions. In the following section we will explain a better algorithm for subdividing Bézier curves.

Exercises

[4.2]What kind of curve does the function \(\vec {C}(t) = \lerp (\vec {p}_0, \vec {p}_1, t)\) describe?

4.3 De Casteljau’s Algorithm¶

Bézier curves and Hilbert curves have something important in common: they can both be defined recursively. Similar to the way a Hilbert curve consists of multiple smaller Hilbert curves of lower order, a Bézier curve can be subdivided into multiple smaller Bézier curves with different control points.

One way to perform this subdivision is an elegant algorithm known as de Casteljau’s algorithm. This algorithm constructs points on a Bézier curve geometrically — without evaluating the polynomial in Eq. (4.2) — by repeatedly averaging the four control points \(\vec {p}_0, \dots , \vec {p}_3\). In the first step we take each pair of adjacent control points and compute their average:

\begin{equation*} \vec {q}_0=\frac {\vec {p}_0+\vec {p}_1}{2},\qquad \vec {q}_1=\frac {\vec {p}_1+\vec {p}_2}{2},\qquad \vec {q}_2=\frac {\vec {p}_2+\vec {p}_3}{2}. \end{equation*}

In the second step we average those three points again to obtain two points:

\begin{equation*} \vec {r}_0=\frac {\vec {q}_0+\vec {q}_1}{2},\qquad \vec {r}_1=\frac {\vec {q}_1+\vec {q}_2}{2}. \end{equation*}

And in the last step we average one more time to obtain a single point

\begin{equation*} \vec {s} = \frac {\vec {r}_0 + \vec {r}_1}{2}. \end{equation*}

As illustrated in Fig. 4.5, this final point lies on the Bézier curve; in fact, it is easy to show that we have \(\vec {s}=\vec {B}(1/2)\) (see Exercise 4.8).

The point \(\vec {s}\) splits the original curve into two halves, \(\vec {B}_1(t)\) and \(\vec {B}_2(t)\), both of which are Bézier curves in their own right. As suggested by Fig. 4.5, the control points of the first half are \(\vec {p}_0\), \(\vec {q}_0\), \(\vec {r}_0\), and \(\vec {s}\), and those of the second half are \(\vec {s}\), \(\vec {r}_1\), \(\vec {q}_2\), and \(\vec {p}_3\), so we have

\begin{equation} \label {eq:beziersplit} \begin{aligned} \vec {B}_1(t) &= (1-t)^3 \vec {p}_0 + 3t(1-t)^2 \vec {q}_0 + 3t^2(1-t)\vec {r}_0 + t^3 \vec {s},\\ \vec {B}_2(t) &= (1-t)^3 \vec {s} + 3t(1-t)^2 \vec {r}_1 + 3t^2(1-t)\vec {q}_2 + t^3 \vec {p}_3. \end {aligned} \end{equation}

The split() method in Listing 4.4 uses this equation to split a given BezierCurve into two halves. If we apply de Casteljau’s algorithm repeatedly, we can subdivide a given Bézier curve into 4, 8, 16, and ultimately as many pieces as we like; the process is illustrated in Fig. 4.6.

// Split Bézier curve in the middle, producing two smaller curves.
public BezierCurve[] split() {
    Point q0 = mid(p0, p1), q1 = mid(p1, p2), q2 = mid(p2, p3);
    Point r0 = mid(q0, q1), r1 = mid(q1, q2);
    Point s = mid(r0, r1);
    return new BezierCurve[]{
            new BezierCurve(p0, q0, r0, s),
            new BezierCurve(s, r1, q2, p3)};
}

// Compute the midpoint of a and b.
static Point mid(Point a, Point b) {
    return Point.lerp(a, b, 0.5);
}

Adaptive Subdivision

As we just saw, de Casteljau’s algorithm gives us another way to compute points on a Bézier curve. However, the main advantage of this algorithm is that it doesn’t just give us points on the curve but produces a hierarchy of curve segments that are all Bézier curves in their own right. This makes it possible to apply the divide-and-conquer paradigm to many problems that occur in connection with Bézier curves. In this section, we show how to implement a flexible algorithm for drawing Bézier curves that selectively subdivides curve segments that aren’t “straight enough” until the entire curve can be approximated by a sequence of line segments. This approach is known as adaptive subdivision.

How do we test the straightness of a Bézier curve? A simple criterion is based on the following observation: If the inner control points \(\vec {p}_1\) and \(\vec {p}_2\) lie at \(1/3\) and \(2/3\) of the distance between the two end points \(\vec {p}_0\) and \(\vec {p}_3\), it can be shown that the resulting Bézier curve is mathematically identical to a line segment (see Exercise 4.8). We can therefore consider a Bézier curve almost straight if the first control point \(\vec {p}_1\) is close to \(\lerp (\vec {p}_0, \vec {p}_3, 1/3)\) and the second control point \(\vec {p}_2\) is close to \(\lerp (\vec {p}_0, \vec {p}_3, 2/3)\), which means that their distance is less than a given tolerance \(T\):

\begin{equation} \label {eq:bezier:straight} \bigl |\vec {p}_1 - \lerp (\vec {p}_0, \vec {p}_3, 1/3)\bigr | < T \qquad \text {and}\qquad \bigl |\vec {p}_2 - \lerp (\vec {p}_0, \vec {p}_3, 2/3)\bigr | < T. \end{equation}

By setting \(T\) to a sufficiently small value, say, a value that is smaller than the size of a pixel, we can distinguish between curve segments that are visually indistinguishable from a line segment and segments that need further subdivision.

By combining this straightness test with de Casteljau’s subdivision algorithm, we can write a function subdivide() that recursively subdivides a Bézier curve until all segments satisfy the straightness criterion in Eq. (4.6); see Listing 4.5. The parameter tolerance specifies the desired level of straightness and corresponds to the parameter \(T\) in Eq. (4.6). The function produces a sequence of points on the curve and adds them to the final argument result.

// Recursively subdivide curve until all segments are straight.
void subdivide(double tolerance, List<Point> result) {
    if (isStraight(tolerance)) {
        result.add(p3);
    } else {
        BezierCurve[] segments = split();
        segments[0].subdivide(tolerance, result);
        segments[1].subdivide(tolerance, result);
    }
}

boolean isStraight(double tolerance) {
    double d1 = Point.distance(p1, Point.lerp(p0, p3, 1.0 / 3.0));
    double d2 = Point.distance(p2, Point.lerp(p0, p3, 2.0 / 3.0));
    return d1 < tolerance && d2 < tolerance;
}

The function first uses isStraight() to test whether the current curve must be subdivided further. If the current curve segment is sufficiently straight, subdivide() stops the subdivision process and appends the end point of the current segment p3 to the list of points. Otherwise, it splits the curve in half using de Casteljau’s algorithm and recursively calls itself to further subdivide the two halves. Notice that subdivide() only adds one of the two end points (p3) of a segment to result. The reason is that the end point of the each segment is also the starting point of the next, and we don’t want to add every interior point twice. The starting point of the original curve must be added to result before starting the subdivision process, as shown in Listing 4.6.

public List<Point> draw(double tolerance) {
    List<Point> result = new ArrayList<>();
    result.add(p0);
    subdivide(tolerance, result);
    return result;
}

The images in Fig. 4.7 illustrate how the straightness tolerance influences the final subdivision of a given Bézier curve. For \(T=1\), we obtain four line segments which corresponds to two applications of Casteljau’s algorithm, but as we decrease the value of \(T\), additional points on the curve are generated to produce a finer subdivision. Notice that decreasing the tolerance always produces a refinement of the previous subdivision: More points are added in certain parts of the curve as necessary, but existing points are never removed. It’s also easy to see that the subdivision is indeed adaptive, so the point density is highest in regions where the curve is strongly bent.

There is one important question we haven’t addressed yet: With this straightness test, can we be sure that the subdivision process always terminates, or are there any Bézier curves that can be subdivided endlessly, without ever being considered straight enough? Fortunately, Bézier curves are sufficiently well-behaved that this cannot happen. Every time we split a Bézier curve into two smaller segments, the control points of the new curves move closer together (and closer to the curve itself) because they are computed by averaging the old control points; this effect can be easily observed in Fig. 4.6. For every fixed tolerance \(T\), the conditions in Eq. (4.6) must therefore be satisfied after a finite number of iterations.

Exercises

⋆ [4.3]How do we prove that the two Bézier curves in Eq. (4.5) together form the original curve \(\vec {B}(t)\)? The key is to realize that the first half of the original curve can be written as \(\vec {B}(t/2)\) and the second half as \(\vec {B}(t/2 + 1/2)\). We then have to show that

\begin{equation*} \vec {B}(t/2) = \vec {B}_1(t)\qquad \text {and}\qquad \vec {B}(t/2+1/2) = \vec {B}_2(t). \end{equation*}

Complete the proof by expanding both sides of these equations using Eq. (4.2) for \(\vec {B}(t)\) and Eq. (4.5) for \(\vec {B}_1(t)\) and \(\vec {B}_2(t)\).

[4.4]Modify de Casteljau’s algorithm to compute the approximate length of a Bézier curve.

4.4 Problems¶

The Koch Snowflake

[4.5]The Koch snowflake in Fig. 4.8 is a fractal curve that resembles a snowflake. The curve can be defined recursively:

• • The order-0 snowflake is an equilateral triangle.

• • The order-\(n\) snowflake is obtained from the order-\(n-1\) snowflake by dividing each edge into three parts, and replacing the middle part with the tip of a small triangle.

Write a program that computes the Koch snowflake of arbitrary order \(n\).

Closest Points on Lines and Bézier Curves

[4.6]A common problem in graphical applications is to determine the distance between a given point \(\vec {m}\) (for example the mouse cursor) and certain geometric objects. In this exercise, we consider the distance to lines, line segments, and Bézier curves.

Let’s start with lines. A line that goes through two points \(\vec {a}=(a_x,a_y)\) and \(\vec {b}=(b_x,b_y)\) can be described as a parametric curve of the following form:

\begin{equation} \label {eq:line} \vec {x}(t) = \vec {a} + t (\vec {b}-\vec {a}). \end{equation}

The point \(\vec {x}_c\) on this line that is closest to \(\vec {m}\) is the one for which the segment connecting the two points is perpendicular to the line itself:

Equivalently, the vector \(\vec {m}-\vec {x}_c\) must be perpendicular to the direction of the line \(\vec {b}-\vec {a}\), and since two vectors are perpendicular if their dot product is zero, we have the following condition:

\begin{equation} \label {eq:lineperp} (\vec {m} - \vec {x}_c)\cdot (\vec {b}-\vec {a}) = 0. \end{equation}

• a) Plug Eq. (4.7) into Eq. (4.8) and solve for \(t\) to find a formula for the closest point.

For line segments we can use the same approach with a small modification. The formula for a line segment from \(\vec {a}\) to \(\vec {b}\) is identical to that of a line, but now the parameter is restricted to the interval between 0 and 1:

\begin{equation} \label {eq:lineseg} \vec {x}(t) = \vec {a} + t (\vec {b}-\vec {a}),\qquad 0\le t\le 1. \end{equation}

To find the closest point on the line segment to a point \(\vec {m}\), we first compute the parameter \(t\) of the closest point along the line through \(\vec {a}\) and \(\vec {b}\), as above. If the resulting parameter satisfies \(0\le t\le 1\), the corresponding point \(\vec {x}_c\) lies on the line segment and is the desired closest point. Otherwise, if \(t<0\) or \(t>1\), the closest point is one of the end points \(\vec {a}\) or \(\vec {b}\).

• b) Write a class LineSegment that represents a line segment between two Points. Define the following two methods:

  Point pointAt(double t);
  Point closestPoint(Point m);
  

  The first method pointAt() implements Eq. (4.9) and returns the point corresponding to a given parameter t. The second method closestPoint() computes the closest point to m on the line segment.

Finding the closest point on a Bézier curve is more difficult. While it is possible to derive a mathematical criterion similar to Eq. (4.8), the resulting equation is a quintic (order 5) polynomial in \(t\) that cannot be solved analytically. But if we are content with an approximate result, we can find the closest point to \(\vec {m}\) by subdividing the curve into line segments, finding the closest point on each, and returning the overall minimum.

• c) Devise a recursive algorithm for computing the point on a Bézier curve that is closest to a given point \(\vec {m}\) by extending de Casteljau’s algorithm.

Bézier Curves Continued

[4.7]As we saw in the text, cubic Bézier curves can be constructed by repeatedly averaging four given control points using de Casteljau’s algorithm. If we start with a different number of control points, we obtain Bézier curves of a different order. In particular, if we start three control points \(\vec {p}_0\), \(\vec {p}_1\), and \(\vec {p}_2\), de Casteljau’s algorithm produces a quadratic Bézier curve. An example of a quadratic Bézier curve is shown in Fig. 4.9.

• a) Show that the construction produces a point that lies on the curve

  \(\seteqnumber{0}{4.}{9}\)

  \begin{equation} \label {eq:qbezier} \vec {B}(t)=(1-t)^2\vec {p}_0 + 2t(1-t)\vec {p}_1+t^2\vec {p}_2 \end{equation}

• b) Devise a criterion for determining whether a quadratic Bézier curve is approximately straight.

• c) Develop a recursive algorithm for adaptively subdividing quadratic Bézier curves.

⋆ [4.8]In the text we stated and used a few properties of Bézier curves without proof, most importantly the relationship between the parametric form of Bézier curves in Eq. (4.2) and de Casteljau’s algorithm.

• a) Show that de Casteljau’s method shown in Fig. 4.5 yields the point \(\vec {B}(1/2)\) on the Bézier curve.

• b) Use a similar derivation to show that any point \(\vec {B}(t)\) can be computed by starting with \(\vec {p}_0, \dots , \vec {p}_3\) and repeatedly combining adjacent points \(\vec {a},\vec {b}\) using \(\lerp (\vec {a}, \vec {b}, t)\).

• c) Show that a Bézier curve whose two control points are placed \(1/3\) and \(2/3\) of the way between its two end points is a line segment.

• d) Given an affine transformation \(A(\vec {x})=\vec {A}\vec {x}+\vec {a}\) and a Bézier curve \(\vec {B}(t)\) specified by four control points \(\vec {p}_0,\dots ,\vec {p}_3\), show that the transformed Bézier curve \(\vec {B}'(t)=A(\vec {B}(t))\) is equal to the Bézier curve that corresponds to the transformed control points \(\vec {p}_k'=A(\vec {p}_k)\). (What this means is that Bézier curves can be rotated, scaled, and sheared simply by transforming their control points.) Hint: Use Eq. (4.4).

Penrose Tilings

[4.9]An arrangement of shapes that covers a surface seamlessly and without overlap is called a tiling or a tessellation. The simplest tilings — commonly used in bathrooms and board games — use regular polygons such as triangles, squares, or hexagons to cover the plane. A famous example of a more complicated tiling is the Penrose tiling shown in Fig. 4.10.

Penrose tilings consist of two different shapes, a thick rhomb whose inner angles are \(72\) and \(108\) degrees and a thin rhomb whose inner angles are \(36\) and \(144\) degrees (Fig. 4.10). To obtain a valid Penrose tiling, these two kinds of rhombs must be combined in such a way that the arrow symbols on all shared edges point in the same direction and have the same number of arrow heads. It’s these placement rules that give Penrose tilings their many fascinating properties. For instance, it can be shown that it is possible to cover the entire plane with Penrose tiles and that the resulting pattern is always aperiodic — it extends indefinitely but never repeats itself.¹

• a) One recurring pattern in Fig. 4.10 is a five-pointed star that consists of five thick rhombs. In such a star, do the arrows on the shared edges point inwards or outwards? Is it possible to build a similar star out of ten thin rhombs?

An important property of Penrose tilings is that they can be subdivided recursively: Given any legal arrangement of Penrose tiles, it is possible to split and reassemble the rhombs in such a way that the result is again a valid arrangement of smaller tiles. This process is known as deflation since it decreases the size of the rhombs.

The deflation rules are illustrated in Fig. 4.11. First, we split each rhomb into two triangles that are mirror images of each other and mark the common edge with a new arrow symbol that allows us to later reassemble the triangles into rhombs. We then subdivide each triangle into multiple smaller triangles as shown in the bottom part of the figure. This subdivision splits existing edges according to the golden ratio \(1:\bigl (1+\sqrt {5}\bigr )/2\). Notice that each new triangle is a scaled-down version of one of the original triangles.

You can convince yourself that this subdivision scheme is consistent with the rules for Penrose tilings. For example, every outer edge with a single arrow “\(\rangle \)” is transformed into two edges, one with three reversed arrows “\(\langle \langle \langle \)” and one with two reversed arrows “\(\langle \langle \)”. We can therefore take any valid arrangement of rhombs, perform any number of subdivision steps, and obtain an arrangement of triangles that can be recombined into a valid Penrose tiling. (The only exception are triangles that lie on the boundary, which may not have a counterpart to combine with.)

• b) Sketch the result of deflating a thick rhomb twice; in other words, extend the left part of Fig. 4.11 with a second subdivision step. How many complete rhombs do you get?

The deflation process can be repeated multiple times, each time producing a more detailed tiling. The tiling shown in Fig. 4.10, for example, was obtained by starting with a single thick rhomb (the outline of the tiling) and deflating it seven times.

• c) Write a program for deflating thick and thin rhombs a given number of times. Hint: It is easier to work with triangles instead of rhombs.

Plotting Parametric Curves

[4.10]A parametric curve is a function that maps an interval of real numbers \([a,b]\) to points in the plane; the parametric form of Bézier curves in Eq. (4.2) is an example of a parametric curve. Figure 4.12 shows two examples of more complicated parametric curves. The first is the so-called “mystery curve” which owes its name to its unexpected five-fold symmetry.² Its shape is described by the following formula:

\begin{equation} \label {eq:mystery} \vec {C}(t) = \begin{pmatrix} \cos t\\ \sin t \end {pmatrix} + \frac {1}{2} \begin{pmatrix} \cos 6t\\ \sin 6t \end {pmatrix} - \frac {1}{3} \begin{pmatrix} \sin (-14t) \\ \cos (-14t) \end {pmatrix},\qquad 0\le t\le 2\pi . \end{equation}

As the parameter \(t\) goes through the range \(0\le t\le 2\pi \), the points \(\vec {C}(t)\) continuously trace out the left curve in the figure. The other curve is the “Butterfly curve,” which is described by the following formula:

\begin{equation} \label {eq:butterflycurve} \vec {C}(t) = \begin{pmatrix} \sin t\\ \cos t \end {pmatrix}\cdot \left [ e^{\cos t} - 2\cos 4t - \sin ^5 (t/12) \right ],\qquad 0\le t\le 12\pi . \end{equation}

(The notation \(\sin ^5x\) is a shorthand for \((\sin x)^5\).)

• a) Design an interface for representing and evaluating parametric curves and use it to implement the mystery curve.

To draw a given parametric curve \(\vec {C}(t)\), we compute a sequence of points \(\vec {p}_1,\dots ,\vec {p}_N\) on the curve and connect them using line segments. Each point \(\vec {p}_k\) is obtained by evaluating the parametric curve at a certain parameter \(t_k\) so that \(\vec {p}_k = C(t_k)\).

The main challenge is to place the points \(\vec {p}_k\) in such a way that connecting them using line segments produces a reasonably accurate image of the curve. As in our discussion of Bézier curves, the simplest approach is to divide the parameter range \([a,b]\) into \(N\) pieces of the same size, which gives us \(t_k = (b-a)/N\) or

\begin{equation*} \vec {p}_k = C\left (k\frac {b-a}{N}\right ). \end{equation*}

The only parameter in this case is the number \(N\) of line segments.

• b) Experiment with drawing the mystery and butterfly curves. How large do you have to choose \(N\) in each case to obtain a reasonably smooth rendition of the curve?

As in our discussion of Bézier curves, there are two main problems with rigid subdivision schemes like this: They provide no control over the placement of points along the curve, and choosing a suitable number of line segments requires a lot of trial and error. A better solution is to perform the subdivision adaptively, by repeatedly splitting curve segments into smallers parts until they are straight enough to be replaced by line segments. Splitting a parametric curve is easily done by splitting its parameter range: If the current parameter range is \(a\le t<b\), we split it into two halves \(a \le t<(a+b)/2\) and \((a+b)/2\le t < b\).

Unfortunately, there is no reliable way to measure the straightness of a general parametric curve since we have no a priori knowledge about the shape of the curve. The best we can do is devise a straightness test that works well for smooth curves and fails gracefully for curves that aren’t well-behaved.

• c) Given a parametric curve \(\vec {C}(t)\) and a parameter range \([a,b]\), suggest a simple heuristic that decides whether the curve segment between \(\vec {C}(a)\) and \(\vec {C}(b)\) is relatively straight.

• d) Write a program that plots a parametric curve for a given parameter range. The algorithm should adaptively subdivide the parameter range until all curve segments are almost straight. What parameters does your algorithm have? How do you ensure that it always terminates?

• e) Use your algorithm to draw the mystery curve and the Butterfly curve. Compare the resulting curves to those you obtained in part (b).

Here are the formulas of two additional parametric curves that you can use for testing:

\begin{align*} \vec {C}(t) &= \begin{pmatrix} A\cos (at + \delta )\\ B\sin (bt) \end {pmatrix} && \text {Lissajous curve}\\ \vec {C}(t) &= \begin{pmatrix} A\cos (at) \cos (t)\\ A\cos (at) \sin (t) \end {pmatrix} && \text {Rhodonea curve} \end{align*} As usual, \(t\) is the curve parameter, and \(A, a, B, b, \delta \) are additional parameters that determine the shape of the curve.

4.5 Chapter Notes¶

Hilbert curves were first described by the German mathematician David Hilbert in 1891, after a similar family of space-filling curves had been discovered the year before by Giuseppe Peano [50, 75]. In general, space-filling curves twist and turn in such a way that they can be made to converge towards two-dimensional shapes. In the case of the Hilbert curve, this limiting shape is the unit square. There are many other kinds of space-filling curves — including the z-order curve, the Peano curve, the Gosper curve, and the dragon curve — and it’s also possible to construct space-filling curves that converge towards cubes and hypercubes in higher dimensions [2].

In addition to being intriguing and pleasant to look at, space-filling curves also have unexpected practical applications [8]. What makes space-filling curves useful in many applications is their ability to translate between a one-dimensional quantity (the distance from the start of the curve) and two-dimensional coordinates (the location of the corresponding point on the curve) in a locality-preserving manner: Due to the way the curves bend and twist, nearby points in the plane tend to be also close to each other along the curve and vice versa [43]. For example, this can be used to rearrange the memory layout of multidimensional arrays so that entries that are close to each other in the array are also close to each other in memory; this usually isn’t the case when storing the array row by row or column by column.

The usefulness of Bézier curves for modeling smooth shapes was recognized early on and helped kick-start the development of computer-aided design during the late 1950s and 1960s. The foundations were laid independently by two French researchers working at competing car companies: Paul de Faget de Casteljau, working at Citroën, and Pierre Étienne Bézier, working at Renault [34]. Even though the curves themselves are now named after Bézier, it was actually de Casteljau who employed this particular type of curve first, and who came up with the crucial idea of shaping them using control points that lie next to the curve, not on it.

It is possible to define Bézier curves of any order. Curves of order \(n\) are specified using \(n+1\) control points, and by averaging adjacent control points \(n\) times we obtain a general version of de Casteljau’s algorithm. Even though higher-order curves are more flexible, cubic (\(n=3\)) and quadratic (\(n=2\)) Bézier curves prevail in most applications since they are easier to work with and not as prone to undesirable oscillations. Cubic Bézier curves are the most common type of curve used in computer graphics and are the basis of many graphics standards such as OpenType fonts and the PDF and SVG file formats. Quadratic Bézier curves, which we discuss in Exercise 4.7, are the basis of TrueType fonts.

In 3D graphics, similar techniques are used to model and render smooth surfaces: These so-called subdivision surfaces are obtained by recursively subdividing an initial mesh of control points into small polygons. There are many different kinds of subdivision surfaces: Bézier surfaces, NURBS, Catmull-Clark surfaces, Loop surfaces, and many more. For more on this subject, see the books by Farin [35] and Akenine-Möller et al. [1].

The Penrose tilings discussed in Exercise 4.9 are named after the British mathematician Roger Penrose who discovered the first variants in the 1970s. The first two chapters of Martin Gardner’s book Penrose Tiles to Trapdoor Ciphers [41] are a good introduction to Penrose tilings and discuss many of their fundamental properties. For a comprehensive account of tilings in general, see the book by Grünbaum and Shephard [45].