How can I draw a Homer with epicycloids? [duplicate]

Question

You may have seen this video: https://www.youtube.com/watch?v=QVuU2YCwHjw.

How does this work, and how to do it with Mathematica?

Answer 1

The main challenge to make something general is to create a list of points that, when joined, is as close as possible to the original image. The first example (elephant) focuses on the mathematical idea. In the second example (Bart Simpson), the code is wrapped in a function. The last example gives some ideas on how to apply the method to a photograph instead of a simple curve.

Explanations on the elephant example

The beginning is pretty similar to this answer: we import an image, extract the points:

img = Import["https://i.sstatic.net/wtJoA.png"];
img = Binarize[img~ColorConvert~"Grayscale"~ImageResize~500~Blur~3];
pts = DeleteDuplicates@
   Cases[Normal@
      ListContourPlot[Reverse@ImageData[img], 
       Contours -> {0.5}], _Line, -1][[1, 1]];
center = Mean@MinMax[pts] & /@ Transpose@pts;
pts = # - center & /@ pts[[;; ;; 20]];
elephantPlot = ListPlot[pts, AspectRatio -> Automatic]

This one is a good example because the points (in the right order!) form a nicely closed curve. Then, the idea is to transform each point as a complex number (z) and take the Discrete Fourier Transform (cn) up to a prescribed order m. This provides a parametric curve of the form:

$$z(t)=\sum_{j=-m}^m c_j e^{-2\pi ij/n}$$

than can be computed and plotted (here, for $m=5$) with:

SetAttributes[toPt, Listable]
toPt[z_] := ComplexExpand[{Re@z, Im@z}] // Chop;
cf = Compile[{{z, _Complex, 1}}, 
    Module[{n = Length@z}, 
1/n*Table[Sum[z[[k]]*Exp[-I*i*k*2 Pi/n], {k, 1, n}], {i, -m, m}]]];
z = pts[[All, 1]] + I*pts[[All, 2]];
m = 5;
cn = cf[z];
{f[t_], g[t_]} = 
Sum[cn[[j]]*Exp[I*(j - m - 1)*t], {j, 1, 2 m + 1}] // toPt;
ParametricPlot[{f[t], g[t]}, {t, 0, 2 Pi}, AspectRatio -> Automatic]

The above sum can be understood as the sum of circles or radii $|c_i|$ with a phase $\arg(c_i)$ and a algebraic angular velocity $i$.

Here is an animation of the result. Note that the sum is commutative to the representation is not unique---e.g. one could plot from the largest circle to the smaller one. Here the order is given by increasing angular velocity.

r = Abs /@ cn;
theta = Arg /@ cn;
index = {m + 1}~Join~
   Riffle[Range[m + 2, 2 m + 1], Reverse[Range[1, m]]];
p[t_] = Accumulate@Table[cn[[j]]*Exp[I*(j - m - 1)*t], {j, index}] // toPt;
circles[t_] = 
  Table[Circle[p[t][[i-1]], r[[index[[i]]]]], {i, 2, 2 m + 1}];

anims = ParallelTable[
   ParametricPlot[{f[s], g[s]}, {s, 0, t}, AspectRatio -> Automatic,
     Epilog -> {circles[t], Line[p[t]], Point[p[t]]}, 
     PlotRange -> {{-400, 400}, {-400, 200}}, 
     ImageSize -> 500],
   {t, Subdivide[0.1, 4 Pi, 100]}];
ListAnimate@anims

Of course this is to illustrate the principle. With m=50, that is 101 circles, one get a much finer result:

Wraping things up with Bart

We'll use this image:

img0 = Import[
  "http://drawinghowtodraw.com/stepbystepdrawinglessons/wp-content/\
uploads/2009/12/finished-bart-simpson.png"]

Then, we extract the contour and using FindShortestTour, we order the points such that joining them gives a sensible image:

center = Mean@MinMax[pts] & /@ Transpose@pts;
pts = # - center & /@ pts;
shortest = (Last@FindShortestTour@pts)[[;; ;; 5]];
pts = pts[[shortest]];
ListLinePlot[pts, AspectRatio -> Automatic]

The eyes or not very good but there is not much to be done: we need to work on a single closed curve. Here pts has more than 1500 points: that can be much improved! But I'll leave that for the gurus...

Then, we can wrap the previous commands in a function that returns the list of points and circles:

compute[pts_, m_] := Block[{z, cn, r, theta, index, tab, p, circles},
  z = pts[[All, 1]] + I*pts[[All, 2]];
  cn = cf[z, m];
  r = Abs /@ cn;
  theta = Arg /@ cn;
  index = {m + 1}~Join~
    Riffle[Range[m + 2, 2 m + 1], Reverse[Range[1, m]]];
  tab = Table[toPt[cn[[j]]*Exp[I*(j - m - 1)*t]], {j, index}];
  p[t_] = Accumulate@tab;
  circles[t_] = Circle @@@ Transpose[{p[t], r[[index]]}][[;; 2 m + 1]];
  {p[t], circles[t]}]

With 300 modes (601 circles):

{p[t_], circles[t_]} = compute[pts, 300];
anims = ParallelTable[
   ParametricPlot[Evaluate[p[s][[-1]]], {s, 0, t}, 
    AspectRatio -> Automatic, Epilog -> {
      circles[t], Line[p[t]], Point[p[t]]}, 
    PlotRange -> {{-200, 100}, {-200, 200}}, ImageSize -> 500, 
    Axes -> False], {t, Subdivide[0.1, 4 Pi, 100]}];

The curve can be plotted with ParametricPlot[Evaluate@p[t][[-1]], {t, 0, 2 Pi}].

Third example: a photography

I chose a picture that is very difficult for this exercise because you can't see an obvious closed curve that would depict the content. But that's why it's an interesting case (and additionally because I love this photography of Bill Evans).

img0 = Import["https://jazzdagama.com/wp-content/uploads/2016/09/Bill-Evans.jpg"]

After a bit of fiddling, I got this:

img = Binarize[
   img0~ColorConvert~"Grayscale"~ImageResize~190~Blur~0, .15];
img = DeleteSmallComponents@img;
pts = DeleteDuplicates@
   Flatten[Cases[
      Normal@ListContourPlot[Reverse@ImageData[img], 
        Contours -> {0.5}], _Line, -1][[All, 1]], 1];
center = Mean@MinMax[pts] & /@ Transpose@pts;
pts = # - center & /@ pts;
plot = ListPlot[pts, AspectRatio -> Automatic, PlotRange -> Full]

Maybe it looks like an abstraction... This might help:

Then we re-order the points to make a closed curve:

shortest = Last@FindShortestTour@pts;
pts = pts[[shortest]];
ListLinePlot[pts, options]

And, finally:

{p[t_], circles[t_]} = compute[pts, 300];
anims = ParallelTable[
   ParametricPlot[Evaluate[p[s][[-1]]], {s, 0, t}, 
    AspectRatio -> Automatic, Epilog -> {
      circles[t], Line[p[t]], Point[p[t]]}, 
    PlotRange -> {{-100, 150}, {-130, 50}}, ImageSize -> 250, 
    Axes -> False], {t, Subdivide[0.1, 4 Pi, 100]}];

To be improved

• automatic detection of the PlotRange
• extraction of points: deletion of points that are too close to each other, making something more robust.
• speed up the computations of the animations by not plotting circles that are smaller than one or two pixels.
• use Fourier instead of the in-house compiled function.