Here at Wolfram Research and at Wolfram|Alpha we love mathematics and computations. Our favorite topics are algorithms, followed by formulas and equations. For instance, Mathematica can calculate millions of (more precisely, for all practical purposes, infinitely many) integrals, and Wolfram|Alpha knows hundreds of thousands of mathematical formulas (from Euler’s formula and BBP-type formulas for pi to complicated definite integrals containing sin(x)) and plenty of physics formulas (e.g from Poiseuille’s law to the classical mechanics solutions of a point particle in a rectangle to the inverse-distance potential in 4D in hyperspherical coordinates), as well as lesser-known formulas, such as formulas for the shaking frequency of a wet dog, the maximal height of a sandcastle, or the cooking time of a turkey.

Recently we added formulas for a variety of shapes and forms, and the Wolfram|Alpha Blog showed some examples of shapes that were represented through mathematical equations and inequalities. These included fictional character curves:

Shape curves:

And, most popular among our users, person curves:

While these are curves in a mathematical sense, similar to say a lemniscate or a folium of Descartes, they are interesting less for their mathematical properties than for their visual meaning to humans.

After Richard’s blog post was published, a coworker of mine asked me, “How can you make an equation for Stephen Wolfram’s face?” After a moment of reflection about this question, I realized that the really surprising issue is not that there is a formula: a digital image (assume a grayscale image, for simplicity) is a rectangular array of gray values. From such an array, you could build an interpolating function, even a polynomial. But such an explicit function would be very large, hundreds of pages in size, and not useful for any practical application. The real question is how you can make a formula that resembles a person’s face that fits on a single page and is simple in structure. The formula for the curve that depicts Stephen Wolfram’s face, about one page in length, is about the size of a complicated physics formula, such as the gravitational potential of a cube.

In this post, I want to show how to generate such equations. As a “how to calculate…”, the post will not surprisingly contain a fair bit of Mathematica code, but I’ll start with some simple introductory explanations.

Assume you make a line drawing with a pencil on a piece of paper, and assume you draw only lines; no shading and no filling is done. Then the drawing is made from a set of curve segments. The mathematical concept of Fourier series allows us to write down a finite mathematical formula for each of these line segments that is as close as wanted to a drawn curve.

As a simple example, consider the series of functions y_n(x),

which is a sum of sine functions of various frequencies and amplitudes. Here are the first few members of this sequence of functions:

Plotting this sequence of functions suggests that as n increases, y_n(x) approaches a triangular function.

The sine function is an odd function, and as a result all of the sums of terms sin(k x) are also odd functions. If we use the cosine function instead, we obtain even functions. A mixture of sine and cosine terms allows us to approximate more general curve shapes.

Generalizing the above (-1)^{(k – 1)/2)} k^-2 prefactor in front of the sine function to the following even or odd functions,

allows us to model a wider variety of shapes:

It turns out that any smooth curve y(x) can be approximated arbitrarily well over any interval [x₁, x₂] by a Fourier series. And for smooth curves, the coefficients of the sin(k x) and cos(k x) terms approach zero for large k.

Now given a parametrized curve γ(t) = {γ_x(t), γ_y(t)}, we can use such superpositions of sine and cosine functions independently for the horizontal component γ_x(t) and for the vertical component γ_y(t). Using a sum of three sine functions and three cosine functions for each component,

covers a large variety of shapes already, including circles and ellipses. The next demonstration lets us explore the space of possible shapes. The 2D sliders change the corresponding coefficient in front of the cosine function and the coefficient in front of the sine function. (Download this post as a CDF to interact)

If we truncate the Fourier expansion of a curve at, say, n terms, we have 4n free parameters. In the space of all possible curves, most curves will look uninteresting, but some expansion coefficient values will give shapes that are recognizable. However, small changes in the expansion coefficients already quickly change the shapes. The following example allows a modification of the first 4 × 16 Fourier series coefficients of a curve (meaning 16 for the x direction and another 16 for the y direction). Using appropriate values for the Fourier coefficients, we obtain a variety of recognizable shapes.

And if we now take more than one curve, we already have all the ingredients needed to construct a face-like image. The following demonstration uses two eyes, two eye pupils, a nose, and a mouth.

And here is a quick demonstration of the reverse: we allow the position of a set of points (the blue crosses) that form a line to be changed and plot the Fourier approximations of this line.

Side note: Fourier series are not the only way to encode curves. We could use wavelet bases or splines, or encode the curves piecewise through circle segments. Or, with enough patience, using the universality of the Riemann zeta function, we could search for any shape in the critical strip. (Yes, any possible [sufficiently smooth] image, such as Jesus on a toast, appears somewhere in the image of the Riemann zeta function ζ(s) in the strip 0 ≤ Re(s) ≤ 1, but we don’t have a constructive way to search for it.)

To demonstrate how to find simple, Fourier series-based formulas that approximate given shapes, we will start with an example: a shape with sharp, well-defined boundaries—a short formula. More concretely, we will use a well-known formula: the Pythagorean theorem.

Rasterizing the equation gives the starting image that we will use.

It’s easy to get a list of all points on the edges of the characters using the function EdgeDetect.

Now that we have the points that form the edges, we want to join them into straight-line (or curved) segments. The following function pointListToLines carries out this operation. We start with a randomly chosen point and find all nearby points (using the function Nearest to be fast). We continue this process as long as we find points that are not too far away. We also try to continue in a “straight” manner by slightly penalizing sharp turns. To see how the curve construction progresses, we use Monitor.

For the Pythagorean theorem, we obtain 13 individual curves from the edge points.

Joining the points and coloring each segment differently shows that we obtained the expected curves: the outer boundaries of the letters, the inner boundaries of the letters a and b, the three squares, and the plus and equal signs.

Now for each curve segment we want to find a Fourier series (of the x and y component) that approximates the segment. The typical textbook definition of the Fourier coefficients of a function f(x) are integrals of the function multiplied by cos(k x) and sin(k x). But at this point we have sets of points, not functions. To turn them into functions that we can integrate, we make a B-spline curve of each curve segment. The parametrization variable of the B-spline curve will be the integration variable. (Using B-splines instead of piecewise linear interpolations between the points will have the additional advantage of making jagged curves smoother.)

We could find the integrals needed to obtain the Fourier coefficients by numerical integration. A faster way is to use the fast Fourier transform (FFT) to get the Fourier coefficients.

To get more uniform curves, we perform one more step: re-parametrize the spline interpolated curve of the given curve segments by arclength. The function fourierComponents implements the B-spline curve making, the re-parametrization by arclength, and the FTT calculation to obtain the Fourier coefficient. We also take into account if a curve segment is open or closed to avoid Gibbs phenomena-related oscillations. (The above demonstration of approximating the pentagram nicely shows the Gibbs phenomenon in case the “Closed” checkbox is unchecked.)

For a continuous function, we expect an average decay rate of 1/k² for the kth Fourier series coefficient. This is the case for the just-calculated Fourier series coefficient. This means that on average the 10th Fourier coefficient is only 1% in magnitude compared with the first one. This decay allows us to truncate the Fourier series at a not too high order, as we do not want to obtain formulas that are too large. This expression gives the exponent in the decay rate of the Fourier components for the a² + b² = c² curve above. (The slightly lower than 2 exponent arises from the discretization points in the B-spline curves.)

Here is a log-log-plot of the absolute values of the Fourier series coefficient for the first three curves. In addition to the general trend of an approximately quadratic decay of the Fourier coefficients, we see that the magnitude of nearby coefficients often fluctuates by more than an order of magnitude.

Multiplying the Fourier coefficients by cos(k t) and sin(k t) and summing the terms gives us the desired parametrizations of the curves.

The function makeFourierSeriesApproximationManipulate visualizes the resulting curve approximations as a function of the series truncation order.

For the Pythagorean theorem, starting with a dozen ellipses, we quickly form the characters of the inequality with increasing Fourier series order.

We want a single formula for the whole equation, even if the formula is made from disjoint curve segments. To achieve this, we use the 2π periodicity of the Fourier series of each segment to plot the segments for the parameter ranges [0, 2π], [4π, 6π], [8π, 10π], …, and in the interleaving intervals (2π, 4π), (6π, 8π), …, we make the curve coordinates purely imaginary. As a result, the curve cannot be drawn there, and we obtain disjoint curve segments. Here this construction is demonstrated for the case of two circles:

The next plot shows the real and imaginary parts of the complex-valued parametrization independently. The red line shows the purely imaginary values from the parameter interval [2π, 4π].

As we want the final formula for the curves to look as short and as simple as possible, we change sums of the form a cos(k t) + b sin(k t) to A sin(k t + φ) using the function sinAmplitudeForm and round the floating-point Fourier series coefficients to nearby rationals. Instead of Piecewise, we use UnitStep in the final formula to separate the individual curve segments. The real segments we list in explicit form, and all segments that should not be drawn are encoded through the θ(sgn(sin(t/2)^(1/2))) term.

Now we have everything together to write down the final parametrization {x(t), y(t)} of the typeset form of the Pythagorean theorem as a mathematical formula.

After having discussed the principal construction idea for the parametrizations, let’s look at a more fun example, say the Pink Panther. Looking at the image search results of the Bing search engine, we quickly find an image that seems desirable for a “closed form” parametrization.

Let’s use the following image:

We apply the function EdgeDetect to find all edges on the panther’s face.

Connecting the edges to curve segments yields about 20 segments. (Changing the optional second and third argument of pointListToLines, we obtain fewer or more segments.)

Here is each segment shown with a different color. We see that some closed curves arise from two joined curve segments; we could separate them by changing the second argument of pointListToLines. But for the goal of sketching a line drawing, the joined curve will work just fine.

Proceeding now as above, it is straightforward to approximate the curve segments by trigonometric series.

Plotting the series shows that with 20 terms per segment, we obtain a good representation of the Pink Panther.

As some of the segments of the panther’s face are more intricate than others, we define a function makeSegmentOrderManipulate that allows the number of terms of the Fourier series for each segment to be varied. This lets us further reduce the size of the resulting parametrizations.

We use initial settings for the number of Fourier coefficients that yield a clearly recognizable drawing of the Pink Panther.

For simple cases, we can now roll up all of the above function into a single function. The next function makeSilhouetteFourierResult takes a string as an argument. The function then 1) performs a search on Bing’s image site for this string; 2) selects an image that seems appropriate from the algorithmic point of view; 3) calculates the Fourier series; and 4) returns as the result plots of the Fourier series and an interactive version that lets us change the order of the Fourier series. For simplicity, we restrict the program to only single curves. (In the web search, we use the word “silhouette” to mostly retrieve images that are formed by just a single curve.) As the function relies on the result of an external search engine, there is no guarantee that the function will always return the wanted result.

Here are three examples showing the function at work. We build the Fourier series for a generic spider, Spiderman’s spider, a couple dancing the tango, and a mermaid. (Evaluating these functions might give different results, as the search engine results might change over time.)

So far, the initial line segments were computed from images. But we can also start with hand-drawn curves. Assume we want a formula for Isaac Newton. As I am not good at drawing faces, I cheated a bit and used the curve drawing tool to draw characteristic facial and hair curves on top of an image of Sir Isaac. (For algorithmic approaches on how to extract feature lines from faces, see the recent paper by Zhao and Zhu.) Here is the image that we will use:

Fortunately, small random wiggles in the hand-drawn curve will not matter, as they will be removed by omitting higher-order terms in the Fourier series.

To better see the hand-drawn curves, we separate the image and the lines.

This time, we have 16 segments. We build their Fourier series.

And here are again various approximation orders of the resulting curves.

We use different series orders for the various segments. For the hair, we use relatively high orders, and for the eyes relatively low orders. This will make sure that the resulting equations for the face will not be larger than needed.

Here are the first 50 approximations shown in one graphic with decreasing opacity of the curves, and with each set of corresponding curve segments shown in a different color.

This gives the following final curve for Sir Isaac’s portrait.

And this is the plotted form of the last formula.

This ends today’s discussion about how to make curves that resemble people’s faces, fictional characters, animals, or other shapes. Next time, we will discuss the endless graphics capabilities that arise from these formulas and how you can use these types of equations in a large variety of images.

To interact with the examples above, first, download the Wolfram CDF Player. Then, download this post as a Computable Document Format (CDF) file.