mathr / blog / #

Interior distance estimates can be calculated for the Mandelbrot set using the formula:

\[\frac{1-\left|\frac{\partial}{\partial{z}}f_c^p(z_0)\right|^2}{\left|\frac{\partial}{\partial{c}}\frac{\partial}{\partial{z}}f_c^p(z_0) + \frac{\frac{\partial}{\partial{z}}\frac{\partial}{\partial{z}}f_c^p(z_0) \frac{\partial}{\partial{c}}f_c^p(z_0)} {1-\frac{\partial}{\partial{z}}f_c^p(z_0)}\right|}\]

where \(f_c^p(z_0) = z_0\). Obtaining \(z_0\) by iteration of \(f_c\) is impractical, requiring a huge number of iterations, moreover that leaves \(p\) still to be determined. Assuming \(p\) is known, it's possible to find \(z_0\) more directly by using Newton's method to solve the equation:

\[z_0^{m+1} = z_0^{m} - \frac{f_c^p(z_0^m) - z_0^m}{\frac{\partial}{\partial{z}}f_c^p(z_0^m) - 1}\]

which can be implemented in C99 like this:

#include <complex.h>

int attractor
  ( complex double *z_out
  , complex double *dz_out
  , complex double z_in
  , complex double c
  , int period
  ) {
  double epsilon = 1e-10;
  complex double zz = z_in;
  for (int j = 0; j < 64; ++j) {
    complex double z = zz;
    complex double dz = 1;
    for (int i = 0; i < period; ++i) {
      dz = 2 * z * dz;
      z = z * z + c;
    }
    complex double zz1 = zz - (z  - zz) / (dz - 1);
    if (cabs(zz1 - zz) < epsilon) {
      *z_out = z;
      *dz_out = dz;
      return 1;
    }
    zz = zz1;
  }
  return 0;
}

The derivative is returned because a simple test for interior-hood is \(\left|\frac{\partial}{\partial{z}}f_c^p(z_0)\right| \le 1\) (this follows directly from the interior distance formula, which gives a non-negative distance if the point is interior). Now the problem is finding the period \(p\). Recalling the atom domain representation function, a point \(c\) is in a domain of period \(p\) when \(\left|f_c^p(0)\right| < \left|f_c^q(0)\right| \text{ for all } 1 \le q < p\). These domains completely enclose the hyperbolic components of the same period, so good guesses for the period of an enclosing hyperbolic component of \(c\) are the partials, namely the values of \(p\) for which \(\left|f_c^p(0)\right|\) reach a new minimum. Here's what the first few atom domains look like:

Pseudo-code for interior distance estimate rendering now looks something like:

dc := 0
z := 0
m := infinity
p := 0
for (n := 1; n <= maxiters; ++n) {
  dc := 2 * z * dc + 1
  z := z^2 + c
  if (|z| < m) {
    m := |z|
    p := n
    if (attractor(&z0, &dz0, z, c, p)) {
      if (|dz0| <= 1) {
        // point is interior with period p and known z0
        // compute interior distance estimate
        break
      }
    }
  }
  if (|z| > R) {
    // point is exterior
    // compute exterior distance estimate from z and dc
    break
  }
}

Computing the interior distance estimate with known \(p\) and \(z_0\) is quite simple:

double interior_distance(complex double z0, complex double c, int period) {
  complex double z = z0;
  complex double dz = 1;
  complex double dzdz = 0;
  complex double dc = 0;
  complex double dcdz = 0;
  for (int p = 0; p < period; ++p) {
    dcdz = 2 * (z * dcdz + dz * dc);
    dc = 2 * z * dc + 1;
    dzdz = 2 * (dz * dz + z * dzdz);
    dz = 2 * z * dz;
    z = z * z + c;
  }
  return (1 - cabs(dz) * cabs(dz)) / cabs(dcdz + dzdz * dc / (1 - dz));
}

So far we have two algorithms for rendering the Mandelbrot set, the first plain just computes exterior distance estimates with no interior checking, and the second unbiased checks for interior-hood every time a partial is reached. One would think the extra interior checks would slow down the rendering, but unbiased is actually faster than plain for the default view because lots of pixels are interior with low periods, reducing the total number of iterations required (with plain all interior points are iterated up to the maximum iteration limit). In the benchmark graph here, plain is red and unbiased is green.

However, zooming in to a region with no interior pixels visible shows a different result entirely: the interior checks are all useless, and do indeed slow down the rendering dramatically.

The solution is a modified rendering algorithm biased, which uses local connectedness properties to postpone interior checking in exterior regions. We keep track of the outcome of the previous pixel (whether it was interior or exterior) and use that to adjust the calculations - if the previous pixel was interior, perform interior checking as in unbiased, but if the previous pixel was exterior, instead of performing interior checking for each period as it is encountered, store the inputs to the interior checking and carry on. Postponing the expensive interior checks until the maximum iteration count has been reached means that for exterior points that escape before the maximum iteration count has been reached, we don't need to perform the interior checking at all. In the benchmark graphs, biased is blue, and you can see that it improves performance to near that of plain for the exterior-dominated embedded julia set view.

Similar speed improvements result for most views, and this can be explained by the mostly-green images - they are plots of the algorithm performance - points are green when the bias from the previous pixel was accurate to determine the outcome from the current pixel, other colours indicate that the wrong algorithm was chosen (such as iterating all the way to maxiters only to find later that the point was interior, or performing expensive interior checks only to find later that the point was exterior).

You can download the full source code for the program used to render the images in this post. The code is slightly awkward because it includes runtime checks for which algorithm is being used - but now we have shown that the biased method is best we could strip out the checks and use the biased method always.

However, there is a flaw - for deep images (close to the resolution of double) where exterior distance estimates have no problems, the interior distance technique can fail in spectacularly ugly fashion:

It's usually possible to work around this by throwing extra precision at the problem (whether with double-double or MPFR-style arbitrary precision software floating point). Alternatively, with perturbation technique based renderers, given a reference at the nucleus of the dominating island, it's possible to compute interior distances using perturbed double-precision orbits. But this post is long enough, more on that another time perhaps.