Wolfram Computation Meets Knowledge

Making Formulas… for Everything—From Pi to the Pink Panther to Sir Isaac Newton

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:

fictional character curves

fictional character curves

Shape curves:

shape curves

shape curves

And, most popular among our users, person curves:

person curve

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.

Finding the formula for the curve that depicts Stephen Wolfram's face

Formula for the curve that depicts Stephen Wolfram's face

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 yn(x),

Example function

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

Finding the sequence of sine functions

Sine functions

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

Plot[sinSums, {x, -Pi, Pi}, PlotRange → All]

Plot of the sequence of sine functions

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,

New functions

allows us to model a wider variety of shapes:

Demonstrating general shape curves

It turns out that any smooth curve y(x) can be approximated arbitrarily well over any interval [x1, x2] 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,

Equations for the horizontal and vertical components

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)

Exploring the space of possible shapes

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.

Example allowing a modification of the first 4 × 16 Fourier series coefficients of a curve

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.

Using more than one curve to construct a face

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.

Allowing the position of a set of points (the blue crosses) to 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.

PythagoreanTheoremTypeset = HoldForm[a^2 + b^2 == c^2]; TraditionalForm[Style[PythagoreanTheoremTypeset, 60]]

a^2 + b^2 = c^2

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

Rasterizing the equation

a^2 + b^2 = c^2

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

EdgeDetect[image] // Show[#, ImageSize → 240] &

Edges of the characters

edgePoints = {#2, -#1} & @@@     Position[ImageData[EdgeDetect[image]], 1, {2}];

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.

Using Monitor to see how the curve construction progresses

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

SeedRandom[22]; hLines = pointListToLines[edgePoints, 6]; Length[hLines]


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.

Graphics[{ColorData["DarkRainbow"][RandomReal[]], Line[#]} & /@    hLines]

a^2 + b^2 = c^2

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.)

Graphics[BSplineCurve[#, SplineDegree → 6, SplineClosed → True] & /@    hLines]

a^2 + b^2 = c^2

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.)

Re-parametrizing curve by arclength

Adding options

fCs = fourierComponents[hLines,     "OpenClose" → Table["Closed", {Length[hLines]}]];

For a continuous function, we expect an average decay rate of 1/k2 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 a2 + b2 = c2 curve above. (The slightly lower than 2 exponent arises from the discretization points in the B-spline curves.)

(Mean[#] ± StandardDeviation[#] ) &[   Coefficient[    Fit[Log[Rest[MapIndexed[{1. #2[[1]], #1} &, #]]], {1, x}, x] & /@      Flatten[Abs[Flatten[#[[2]], 1]] & /@ fCs, 1], x, 1]] //   NumberForm[#, 3] &

-1.74 ± 0.233

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.

Building a log-log-plot of the absolute values of the Fourier series coefficient for the first three curves

Log-log-plot of the absolute values of the Fourier series coefficient for the first three curves

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

Multiplying the Fourier coefficients by cos(k t) and sin(k t) and summing the terms

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

Using the function makeFourierSeriesApproximationManipulate

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

makeFourierSeriesApproximationManipulate[fCs, 50]

Forming 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:

Making the curve coordinates purely imaginary

Plotting the two circles

Plot 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π].

Building a plot showing the real and imaginary parts of the complex-valued parametrization independently

Plot showing the real and imaginary parts of the complex-valued parametrization independently

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.

Using the function sinAmplitudeForm

Separating the individual curve segments

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.

finalCurve = Rationalize[singleParametrization[fCs, t, 12] , 10^-3]; Short[finalCurve, 12] // TraditionalForm

Final formula

ParametricPlot[Evaluate[finalCurve], {t, 0, 12 4 Pi }]

Parametric plot

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.

Searching for Pink Panther images

Images of the Pink Panther

Let’s use the following image:

Importing image of the Pink Panther

Pink Panther

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

Using EdgeDetect

Edges of Pink Panther image

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.)

edgePoints = {#2, -#1} & @@@     Position[ImageData[pinkPantherEdgeImage], 1, {2}];

SeedRandom[2]; hLines = pointListToLines[edgePoints, 16]; Length[hLines]


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.

Graphics[{ColorData["DarkRainbow"][RandomReal[]], Line[#]} &  /@    hLines]

Edges shown in color

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

fCs = fourierComponents[hLines];

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


Representation of the Pink Panther

Building a grid of graphics

Multiple steps of the Pink Panther image

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.

Reducing 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.

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.

Using the function makeSilhouetteFourierResult

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.)

makeSilhouetteFourierResult["spider", 120]


makeSilhouetteFourierResult["spiderman logo spider", 180]

Spiderman logo spider

makeSilhouetteFourierResult["tango", 160]


makeSilhouetteFourierResult["mermaid", 160]


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:

Importing image of Isaac Newton

Isaac Newton

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.

Defining curveAnnotatedNewtonImage

hLines = Reverse[    SortBy[First /@       Cases[curveAnnotatedNewtonImage, _Line, ∞], Length]];  Length[hLines]


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

Separating the image and the lines

Isaac Newton with lines shown

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

fCs = fourierComponents[hLines];

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

Building various orders of the curves

Various orders of the 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.

makeSegmentOrderManipulate[fCs, newtonOrderList =   Last /@ {{1, 12}, {2, 100}, {3, 60}, {4, 24}, {5, 16}, {6, 16}, {7,       12}, {8, 8}, {9, 8}, {10, 10}, {11, 6}, {12, 6},                       {13, 4}, {14, 8}, {15, 8}, {16, 4}}, 120]

Set of equations creating Isaac Newton's face

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.

Compiling the first 50 approximations shown in one graphic

First 50 approximations shown in one graphic

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

newtonCurve =    Rationalize[singleParametrization[fCs, t, newtonOrderList] ,     0.002]; Style[newtonCurve, 6] // TraditionalForm

Final curve for Isaac Newton's portrait

And this is the plotted form of the last formula.

ParametricPlot[Evaluate[N[newtonCurve]], {t, 0, 64 Pi},   PlotPoints → 240]

Isaac Newton curve

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.


Join the discussion

!Please enter your comment (at least 5 characters).

!Please enter your name.

!Please enter a valid email address.


  1. Nice work Michael.

    You could have fun coming up with “A New Formula for Pi”, where the output was just the outline of either the Greek letter or the English equivalent….

  2. The show is quite interesting. Looking at the pink panther original picture, it seems to me that the basis of functions that our brain is using to sketch shapes is not the fourier one (do not know how to clarify this concept, but perhaps you get me anyway…). With MMA capabilities would it be possible to switch to a different basis and look what would happen?
    My compliments.

  3. i wonder what is the average of two faces, so if we write face1+face2 = result face. could this be implemented or it is already available. a more general what is the face of all faces on the globe. ie to what general face all faces will collapse to. i’m sure this is a very amusing subject. and thanks for your great article.

  4. where can this program be found for use?