Julia walk around the Mandelbrot lake

on 21 Oct 2010 by Mukund (@muks)

I wanted to make a video of a Julia walk around the Mandelbrot lake. It is now finished and available to watch and download.

This is a video of various Julia sets that are produced by making the value of k walk the boundary of the Mandelbrot set. These are some of the most interesting places which create good Julias. By linearly walking the boundary, we can get a smoothly transitioning Julia. But this is not possible as the boundary is infinite. Hence this Julia animation is slightly discontinuous, as even small changes in k in the region of the boundary produce vastly different Julias.

The program that generates the frames of this video is a hack. To find interesting points, it shoots a ray from (0,0) outwards at an angle to intersect with the boundary of the Mandelbrot set. There are better approaches, but this is the one used here. Once the boundary is found, it creates a frame for that point. It repeats this counter-clockwise for all 2\pi angles to produce this 8 minute 24 fps animation.

As the frames are not continuous, to enjoy this video, I suggest the following:

  1. Watch the changing boundaries of the Julias rather than concentrating at the centre of the video. The boundaries display the most interesting changes.
  2. As the video takes 8 minutes to travel 360 degrees, if you are bored by the lack of action, run it at a very high fps (600 fps) to speed it up: mplayer -fixed-vo -fs -fps 600 -loop 0 tifa.webm
  3. Running the animation at 6 fps (instead of the default 24 fps) at 1:1 zoom lets you watch details in the video: mplayer -fps 6 tifa.webm

The dark grayscale colormap is used because this began as a tribute to Mandelbrot. Enjoy.

Sample frames from the video

Julia frame x00000000 Julia frame x00000480 Julia frame x00000960 Julia frame x00001440 Julia frame x00001920 Julia frame x00002400 Julia frame x00002885 Julia frame x00003360 Julia frame x00003840 Julia frame x00004320 Julia frame x00004800 Julia frame x00005280 Julia frame x00005760

How this video was made

The code is a modified version of the Julia program presented a few days ago. It is a quick hack, and no effort was made in optimizing it or making it numerically stable. We are moving and I wanted to spend as little time as possible on it.

The process (single thread) runs for about 21 hours on an Intel i7 920. Memory consumption is anywhere between 15–35 MB. It is CPU intensive code. The computer on which this program was written doesn't like rendering fractals though. :) After about 20 minutes into the run, the CPU probably heats up too much and the computer suddenly powers off. Instead, the entire video was produced on a server in a datacenter.

/* This code was placed in the public-domain by Mukund Sivaraman on
   October 16, 2010. Do whatever you want with it. */

#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <png.h>

/* Boundaries of the render */
#define X1 -1.894
#define X2  1.894
#define Y1 -1.0
#define Y2  1.0

/* Dimensions of the output PPM in pixels */
#define WIDTH  (1920 * 4)
#define HEIGHT (1080 * 4)

/* 8 minutes, 24 frames per second */
#define DURATION (60 * 8)
#define FPS      24
#define FRAMES   (DURATION * FPS)

static int map[256];

static void
make_map (void)
{
  int i;

  for (i = 0; i < 256; i++) {
    double a = M_PI * 2 * i / 256;
    map[i] = (int) fabs (sin (a) * 255);
  }
}

static void
find_interesting_point (double  angle,
                        double *x,
                        double *y)
{
  double d, a, ym;
  double cx;
  int n;

  d = (X2 - X1) / WIDTH;

  if ((angle >=0.0) && (angle < (M_PI * 0.5))) {
    a = angle;
    ym = 1.0;
  } else if ((angle >= (M_PI * 0.5)) && (angle < M_PI)) {
    a = M_PI - angle;
    d = (-d);
    ym = 1.0;
  } else if ((angle >= M_PI) && (angle < (M_PI * 1.5))) {
    a = angle - M_PI;
    d = (-d);
    ym = -1.0;
  } else {
    a = (2 * M_PI) - angle;
    ym = -1.0;
  }

  /* printf ("a = %f\n", a); */
  /* printf ("d = %f\n", d); */
  /* printf ("ym = %f\n", ym); */

  for (cx = 0.0, n = 0; ; cx += d) {
    double cy;
    double zx, zy;
    int k;

    cy = ym * cx * tan (angle);

    zx = cx;
    zy = cy;

    for (k = 0; k < 256; k++) {
      double tx;

      if (((zx * zx) + (zy * zy)) >= 4.0)
        break;

      tx = cx + ((zx * zx) - (zy * zy));
      zy = cy + (2 * zx * zy);
      zx = tx;
    }

    /* printf ("cx = %f, cy = %f, k = %d\n", */
    /*         cx, cy, k); */

    if (k < 256) {
      if (n > 1) {
        *x = cx;
        *y = cy;
        break;
      }
      n++;
    }
  }
}

static void
make_julia (double kx,
            double ky,
            long   frame)
{
  char *filename;
  FILE *out;
  int i, j, k;

  double x1 = X1;
  double x2 = X2;
  double y1 = Y1;
  double y2 = Y2;

  double cx, cy, zx, zy, tx;
  double dx, dy;

  png_structp png_ptr;
  png_infop info_ptr;
  png_bytep *rowp;

  printf ("Making at (%f, %f) for frame %ld\n", kx, ky, frame);

  asprintf (&filename, "x%08ld.png", frame);
  out = fopen (filename, "wb");
  if (!out) {
    fprintf (stderr, "Unable to open output file: %s\n", filename);
    abort ();
  }

  png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING,
                                     NULL, NULL, NULL);
  if (!png_ptr) {
    fprintf (stderr, "Unable to create PNG write struct\n");
    abort ();
  }

  info_ptr = png_create_info_struct (png_ptr);
  if (!info_ptr) {
    fprintf (stderr, "Unable to create PNG info struct\n");
    png_destroy_write_struct(&png_ptr, NULL);
    abort ();
  }

  if (setjmp (png_jmpbuf (png_ptr))) {
    fprintf (stderr, "Error writing PNG file\n");
    png_destroy_write_struct (&png_ptr, &info_ptr);
    abort ();
  }

  rowp = png_malloc (png_ptr, HEIGHT * sizeof (png_bytep));
  for (i = 0; i < HEIGHT; i++)
    rowp[i] = png_malloc (png_ptr, WIDTH);

  png_set_rows (png_ptr, info_ptr, rowp);

  png_init_io (png_ptr, out);
  png_set_IHDR (png_ptr, info_ptr, WIDTH, HEIGHT,
                8, PNG_COLOR_TYPE_GRAY,
                PNG_INTERLACE_NONE,
                PNG_COMPRESSION_TYPE_DEFAULT,
                PNG_FILTER_TYPE_DEFAULT);

  dx = x2 - x1;
  dy = y2 - y1;

  for (j = 0; j < HEIGHT; j++)
    {
      cy = y1 + ((j / (HEIGHT * 1.0)) * dy);

      for (i = 0; i < WIDTH; i++)
        {
          cx = x1 + ((i / (WIDTH * 1.0)) * dx);

          zx = cx;
          zy = cy;

          for (k = 0; k < 255; k++)
            {
              if (((zx * zx) + (zy * zy)) >= 4.0)
                break;

              tx = kx + ((zx * zx) - (zy * zy));
              zy = ky + (2 * zx * zy);
              zx = tx;
            }

          rowp[j][i] = map[k];
        }
    }

  png_write_png (png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL);

  for (i = 0; i < HEIGHT; i++)
    png_free (png_ptr, rowp[i]);
  png_free (png_ptr, rowp);

  png_destroy_write_struct (&png_ptr, &info_ptr);

  fclose (out);
  free (filename);
}

int
main (int argc, char *argv[])
{
  double angle, step;
  long frame;
  double x, y;

  setbuf (stdout, NULL);

  make_map ();

  angle = M_PI;
  step = (2 * M_PI) / FRAMES;

  for (frame = 0; frame < FRAMES; frame++) {
    printf ("Finding interesting point for angle %f\n",
            (360 / (2 * M_PI)) * angle);
    find_interesting_point (angle, &x, &y);

    make_julia (x, y, frame);

    angle += step;
    if (angle >= (2 * M_PI))
      angle = 0.0;
  }

  return 0;
}

Download tifabare.c, tifasmooth.c, tifacolor.c.

The program creates 11520 grayscale PNG files totalling 53 GB in disk usage. These files are 7680×4320 px each in size. They are meant to run for 8 minutes at 24 FPS.

# compile the program
gcc -Wall -msse4 -O2 -funroll-loops -ffast-math -o tifabare tifabare.c -lm -lpng

# run it, producing PNG files in the working dir, and log messages in tifabare-log.txt
time ./tifabare > tifabare-log.txt

# produce 1080p sized frames
mkdir 1080p
for f in x*.png; do echo $f; convert -geometry 1920x1080 $f 1080p/$f; done

# produce webm video file
cd 1080p
ffmpeg -i x%08d.png -r 24 tifa.webm

BTW, for what it's worth, all Banu services including DNS, web, email, lists, XMPP, rsync, git, etc. are now available on native non-tunneled IPv6.