mathr / blog / #

1920x1080p60 MKV (945MB), 4m10s, stereo sound

Julia morphing, an artistic Mandelbrot set zooming technique, gives angled internal addresses that end with a regular structure like:

\[ \cdots p \to_{1/3} 2p+k \to_{1/3} 2(2p+k)+k \cdots \]

Back in 2013, I made an animation of an embedded Julia set orbit in the hairs around several period 7 islands. Embedded Julia sets are (relatively) surface features with simple angled internal addresses, while Julia morphs are much deeper and the addresses are more complicated.

My first attempt simply replaced a prefix of a template Julia morph with 2^N prefixes based on the period patterns in the hairs and while I got images (the angled internal addresses were realizable by actual locations in the Mandelbrot set) they didn't form a coherent animation in any way. I asked a question on mathoverflow about this.

I went back to the graphical Mandelbrot set explorers, and figured out that the suffix needed to change when changing the prefix, because k (as in 2p+k) changes according to which child bulb the hair is attached to. Moreover, there is another case to consider, whether the child bulb's period appears in the address (=> address has more numbers) or not (=> address is the simple 2p+k pattern).

I did the initial graphical analysis in m-perturbator-gtk, heavily using its annotation features such as periodic point marking and external ray tracing to find external angles. Then I used mandelbrot-symbolics in ghci, to convert between angles and addresses. Here's some of the Haskell code:

-- ...

main = do
  [depth] <- map read <$> getArgs
  let prefixes = (".(" ++) . (++ "0111)") . concat <$> replicateM depth [lo, hi]
      julia = 3 * depth + 4
      morph = 4
      lo = "011"
      hi = "100"
      period (Sym.Angled _ _ a) = period a
      period (Sym.Unangled p) = p
  writeFile "a.txt" . unlines . map (\prefix ->
    let Just addressPrefix@(Sym.Angled 1 _ (Sym.Angled 2 _ (Sym.Angled 3 pq _))) = (Sym.angledAddress . Sym.rational) =<< Txt.parse prefix
        addressSuffix (Sym.Angled p t a) = Sym.Angled p t (addressSuffix a)
        addressSuffix (Sym.Unangled p) = Sym.Angled p (1 Sym.% 2) (go morph p0)
          where
            p0 | m < julia = julia + m
               | m > julia = 2 * p
            go 0 p | m < julia = Sym.Unangled p
            go n p | m < julia = Sym.Angled p (1 Sym.% 3) (go (n - 1) (2 * p + m))
            go 0 p | m > julia = Sym.Angled p (1 Sym.% 2) (Sym.Angled (p + 1) (1 Sym.% 2) (Sym.Unangled (p + 2)))
            go n p | m > julia = Sym.Angled p (1 Sym.% 2) (Sym.Angled (p + 1) (1 Sym.% 2) (Sym.Angled (p + 2) (1 Sym.% 3) (go (n - 1) (2 * (p + 2) + m - 2))))
        m = 3 * fromIntegral (Sym.denominator pq)
        address = addressSuffix addressPrefix
        node (Sym.Angled p1 _ (Sym.Unangled p2)) | m < julia = p1 + p2
        node (Sym.Angled p1 _ (Sym.Angled _ _ (Sym.Angled _ _ (Sym.Unangled p2)))) | m > julia = p1 + p2
        node (Sym.Angled _ _ a) = node a
        Just ray = Txt.plain . Sym.binary . <$> Sym.addressAngles address
    in  ray ++ "\t" ++ show (period address) ++ "\t" ++ show (node address)) $ prefixes

That outputs one of the external angles of the angled internal address, of the Julia morph corresponding to each prefix, as well as its period, and the period of a neighbouring node in the Julia morph which I used in some C code to align all the morphs in screen-space. I use the command line tool parallel -k to parallelize tracing these external rays using m-exray-in from mandelbrot-numerics, then found the periodic nucleus of the minibrot islands with m-nucleus, and the approximate zoom level using m-domain-size and some ghc -e scripting (starting ghc is slow, better to run it once to handle multiple inputs). I fed this to some custom C code that aligns the view with the node period:

#include <stdio.h>
#include <mandelbrot-numerics.h>

int main(int argc, char **argv)
{
  if (! (argc > 1)) return 1;
  int Period = atoi(argv[1]);
  // hardcoded precision, should be based on zoom level
  mpfr_t Re, Im, Zoom;
  mpfr_init2(Re, 1000);
  mpfr_init2(Im, 1000);
  mpfr_init2(Zoom, 53);
  // parse KFR input (no error checking)
  getc(stdin); getc(stdin); getc(stdin); getc(stdin);
  mpfr_inp_str(Re, stdin, 10, MPFR_RNDN);
  getc(stdin);
  getc(stdin); getc(stdin); getc(stdin); getc(stdin);
  mpfr_inp_str(Im, stdin, 10, MPFR_RNDN);
  getc(stdin);
  getc(stdin); getc(stdin); getc(stdin); getc(stdin); getc(stdin); getc(stdin);
  mpfr_inp_str(Zoom, stdin, 10, MPFR_RNDN);
  getc(stdin);
  // find nearby node nucleus and use it to align view
  mpfr_d_div(Zoom, 2, Zoom, MPFR_RNDN);
  mpfr_mul_d(Zoom, Zoom, 0.3, MPFR_RNDN);
  mpc_t c1;
  mpc_init2(c1, 1000);
  mpfr_add(mpc_realref(c1), Re, Zoom, MPFR_RNDN);
  mpfr_set(mpc_imagref(c1), Im, MPFR_RNDN);
  m_r_nucleus(c1, c1, Period, 64, 1);
  mpfr_sub(mpc_realref(c1), mpc_realref(c1), Re, MPFR_RNDN);
  mpfr_sub(mpc_imagref(c1), mpc_imagref(c1), Im, MPFR_RNDN);
  mpc_abs(Zoom, c1, MPFR_RNDN);
  mpfr_div_d(Zoom, Zoom, 0.3, MPFR_RNDN);
  mpfr_d_div(Zoom, 2, Zoom, MPFR_RNDN);
  mpfr_printf("Re: %Re\r\nIm: %Re\r\nZoom: %.18e\r\n",
    Re, Im, mpfr_get_d(Zoom, MPFR_RNDN));
  mpc_arg(Zoom, c1, MPFR_RNDN);
  printf("RotateAngle: %.18f\r\nIterations: 15000\r\n",
    mpfr_get_d(Zoom, MPFR_RNDN) / (2 * 3.141592653589793) * 360);
  return 0;
}

Then I rendered the final images with kf.x86_64.exe. All of this was orchestrated by a Bash shell script, a.txt generated by the Haskell with depth 12 (4096 lines) (this took the most time, almost 2 days):

#!/bin/bash
prec=1000

cat a.txt |
while read ray period node
do
  echo "m-exray-in 100 '${ray}' 8 $((2 * 8 * ${period}))"
done |
parallel -k |
paste a.txt - |
tee b.txt

cat b.txt |
while read ray period node re im
do
  echo m-nucleus $prec $re $im $period 64 1
done |
parallel -k |
paste b.txt - |
tee c.txt

cat c.txt |
while read ray period node ere eim nre nim
do
  echo m-domain-size $prec $nre $nim $period
done |
parallel -k |
paste c.txt - | tee d.txt

cat d.txt |
while read ray period node ere eim nre nim ds
do
  echo $ds
done |
ghc -e "interact $ unlines . map (show . (\s -> 2 / (10 * s**1.125)) . read) . lines" |
paste d.txt - |
tee e.txt

cat -n e.txt |
while read n ray period node ere eim nre nim ds zoom
do
  echo -e "Re: $nre\r\nIm: $nim\r\nZoom: $zoom\r" |
  m-align-nodes ${node} > $(printf %04d.kfr $n)
done

for i in ????.kfr
do
  kf.x86_64.exe -s settings.kfs -l "$i" -c palette.kfp -t "$i.tif" 2>&1 |
  grep -v fixme # silence Wine error flood
  convert -colorspace RGB $i.tif -geometry 1920x1080 -colorspace sRGB "${i%kfr}tif" &&
  rm "$i.tif"
done

ffmpeg -r 60 -i "%04d.tif" -pix_fmt yuv420p -profile:v high -level:v 5.1 -crf:v 20 spin.mkv

I asked another question about the magic number 1.125 in this script.

While the rays were tracing (using the CPU) I ran some of the commands on the first ray to get the first KFR file, and I rendered a zoom video (using the GPU) to the minibrot at its center, controlling the zoom speed and depth in zoomasm so that the zoom video slowed down and paused at the exact zoom depth of the spin video (I tried to use maths to work out the numbers to enter into zoomasm, but couldn't figure it out so I eventually used binary search by hand and got it close enough).

I also made a soundtrack, using the pairs of rays of each spin video frame, angles expressed in binary, converted to wavetables and stretched to the length of the frame using FFT/zeropad/IFFT. I processed the result of this in audacity, to extend it to the length of the whole video (making sure to keep the spin section correctly aligned in time):

#include <complex.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sndfile.h>
#include <fftw3.h>

int main(int argc, char **argv)
{
  const int SR = 192000;
  const int FPS = 60;
  const int length = SR / FPS;
  float audio[length][2];
  double _Complex *ifft_in = fftw_malloc(length * sizeof(*ifft_in));
  double _Complex *ifft_out = fftw_malloc(length * sizeof(*ifft_out));
  fftw_plan backward = fftw_plan_dft_1d(length, ifft_in, ifft_out, FFTW_BACKWARD, FFTW_PATIENT | FFTW_DESTROY_INPUT);
  int old_period = 0;
  double _Complex *fft_in = 0;
  double _Complex *fft_out = 0;
  fftw_plan forward = 0;
  SF_INFO onfo = { 0, SR, 2, SF_FORMAT_WAV | SF_FORMAT_FLOAT, 0, 0 };
  SNDFILE *ofile = sf_open("audio.wav", SFM_WRITE, &onfo);
  while (true)
  {
    int period = 0;
    char *ray[2] = { 0, 0 };
    if (3 != scanf("%d\t%ms\t%ms\n", &period, &ray[0], &ray[1]))
    {
      break;
    }
    if (period != old_period)
    {
      if (fft_in)
      {
        fftw_free(fft_in);
      }
      if (fft_out)
      {
        fftw_free(fft_out);
      }
      if (forward)
      {
        fftw_destroy_plan(forward);
      }
      fft_in  = fftw_malloc(period * sizeof(*fft_in));
      fft_out = fftw_malloc(period * sizeof(*fft_out));
      forward = fftw_plan_dft_1d(period, fft_in, fft_out, FFTW_FORWARD, FFTW_MEASURE | FFTW_DESTROY_INPUT);
    }
    for (int c = 0; c < 2; ++c)
    {
      int v = 0;
      for (int p = 0; p < period; ++p)
      {
        fft_in[p] = v;
        v += ray[c][p] == '0' ? -1 : 1;
      }
      double s = v / (double) period;
      double _Complex dc = 0;
      for (int p = 0; p < period; ++p)
      {
        dc += (fft_in[p] -= p * s);
      }
      dc /= period;
      for (int p = 0; p < period; ++p)
      {
        fft_in[p] -= dc;
      }
      fftw_execute(forward);
      memset(ifft_in, 0, length * sizeof(*ifft_in));
      for (int p = 0; p < period/2; ++p)
      {
        ifft_in[p] = fft_out[p];
        ifft_in[length - 1 - p] = fft_out[period - 1 - p];
      }
      fftw_execute(backward);
      double rms = 0;
      for (int p = 0; p < length; ++p)
      {
        audio[p][c] = creal(ifft_out[p]);
        rms += audio[p][c] * audio[p][c];
      }
      rms = 0.25 / sqrt(rms / length);
      for (int p = 0; p < length; ++p)
      {
        audio[p][c] *= rms;
      }
      free(ray[c]);
    }
    sf_writef_float(ofile, &audio[0][0], length);
    old_period = period;
  }
  sf_close(ofile);
  return 0;
}

Then I split the zoom video at this point and splicing in the spin video:

ffmpeg -i zoom.mkv -t  136 -codec:v copy zoom-1.mkv
ffmpeg -i zoom.mkv -ss 136 -codec:v copy zoom-2.mkv
cat > zoom.txt <<EOF
file 'zoom-1.mkv'
file 'spin.mkv'
file 'zoom-2.mkv'
EOF
ffmpeg -f concat -i zoom.txt -i soundtrack.m4a \
  -codec:v copy -codec:a copy "Julia morph orbit in the hairs.mkv"

Then two-pass encoding to a lower bitrate version for web streaming:

ffmpeg -i "Julia morph orbit in the hairs.mkv" \
  -codec:a copy -profile:v high -level:v 5.1 -g 120 \
  -x264opts keyint=120:no-scenecut -b:v 6M -movflags +faststart \
  -pass 1 -y "Julia morph orbit in the hairs.mp4"
ffmpeg -i "Julia morph orbit in the hairs.mkv" \
  -codec:a copy -profile:v high -level:v 5.1 -g 120 \
  -x264opts keyint=120:no-scenecut -b:v 6M -movflags +faststart \
  -pass 2 -y "Julia morph orbit in the hairs.mp4"

And that's the video linked from the image at the top of the page.

Meanwhile, a test of motion blur:

1920x360p10 MP4 (17MB), 13s, no sound

I was considering whether to add motion blur, and if so, how much. 1024 frames in total were rendered, 128 frames in the output video at 10fps. The left video takes groups of 8 frames, outputting the first (0% shutter speed). The middle video takes groups of 8 frames, averaging the first 4 of them (50% shutter speed). The right video takes groups of 8 frames, averaging all of them (100% shutter speed).

Also, a test of soundtrack algorithm:

1920x1080p15 MP4 (37MB), 1m08s, stereo sound

More work in progress, this time on soundtrack. 2^12 pairs of external angles of various periods are sonified as wavetables, all stretched to the same length using FFT, zero-pad, IFFT. Meanwhile the video as 2^10 frames at 15fps, the coordinates of the final set of 2^12 frames at 60fps are still calculating and are yet to be rendered.