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.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)) $ prefixesThat 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.mkvI 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.