# mathr / blog / #

## Julia morph orbit in the hairs

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 2N 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.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)
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_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 |
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.