4
$\begingroup$

From wikibooks colouring the Julia set, the distance estimate $\delta(z)$ can be calculated by:

$$\begin{aligned} \delta(z) &= \lim_{n \to \infty} \frac{|z_n| \log |z_n|}{\left|\frac{\partial}{\partial z} z_n \right|} &\text{ for } z_n \to \infty \\ \delta(z) &= \lim_{n \to \infty} \frac{|z_{pn} - z_*| |\log |z_{pn} - z_*||}{\left|\frac{\partial}{\partial z} z_{pn} \right|} &\text{ for } z_{pn} \to z_* \text{, super-attracting} \\ \delta(z) &= \lim_{n \to \infty} \frac{|z_{pn} - z_*|}{\left|\frac{\partial}{\partial z} z_{pn} \right|} &\text{ for } z_{pn} \to z_* \text{, attracting} \end{aligned}$$

where the first case is exterior distance estimate, and the latter two are interior distance estimates.

The exterior distance estimate works well, and the interior distance estimate works well for the super-attracting case, but for the attracting case there are lots of small spots in the interior. I presume these occur because the distance estimates are only valid in the limit as $z \to J(f)$, but their existence makes using the interior distance estimate rather impractical - I would like to get rid of them, ideally with an improved interior distance estimate formula.

How to find an interior distance estimate that doesn't suffer from spots?

I am primarily interested in Julia sets for the iterated quadratic polynomial $f_c(z) = z^2 + c$. The following examples and source code are all for these $f_c$.

Example images:

c = i, no interior

quadratic Julia set for i

c = 0, super-attracting

quadratic Julia set for 0

c = -1, super-attracting

quadratic Julia set for -1

c = -0.5, attracting

quadratic Julia set for -0.5

c = -0.01, attracting

quadratic Julia set for -0.01

c = -0.99, attracting

quadratic Julia set for -0.99

Full C99 source code used to render the images:

// gcc -std=c99 -Wall -Wextra -pedantic -O3 -o julia-de julia-de.c -lm

#include <complex.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>

void hsv2rgb(double h, double s, double v, int *rp, int *gp, int *bp) {
  double i, f, p, q, t, r, g, b;
  int ii;
  if (s == 0.0) { r = g = b = v; } else {
    h = 6 * (h - floor(h));
    ii = i = floor(h);
    f = h - i;
    p = v * (1 - s);
    q = v * (1 - (s * f));
    t = v * (1 - (s * (1 - f)));
    switch(ii) {
      case 0: r = v; g = t; b = p; break;
      case 1: r = q; g = v; b = p; break;
      case 2: r = p; g = v; b = t; break;
      case 3: r = p; g = q; b = v; break;
      case 4: r = t; g = p; b = v; break;
      default:r = v; g = p; b = q; break;
    }
  }
  *rp = fmin(fmax(round(r * 255), 0), 255);
  *gp = fmin(fmax(round(g * 255), 0), 255);
  *bp = fmin(fmax(round(b * 255), 0), 255);
}

complex double julia_attractor(complex double c, int maxiters, int *period) {
  double epsilon = nextafter(2, 4) - 2;
  complex double z = c;
  double mzp = 1.0/0.0;
  int p = 0;
  for (int n = 1; n < maxiters; ++n) {
    double mzn = cabs(z);
    if (mzn < mzp) {
      mzp = mzn;
      p = n;
      complex double z0 = z;
      for (int i = 0; i < 64; ++i) {
        complex double f = z0;
        complex double df = 1;
        for (int j = 0; j < p; ++j) {
          df = 2 * f * df;
          f = f * f + c;
        }
        complex double z1 = z0 - (f - z0) / (df - 1);
        if (cabs(z1 - z0) <= epsilon) {
          z0 = z1;
          break;
        }
        if (isinf(creal(z1)) || isinf(cimag(z1)) || isnan(creal(z1)) || isnan(cimag(z1))) {
          break;
        }
        z0 = z1;
      }
      complex double w = z0;
      complex double dw = 1;
      for (int i = 0; i < p; ++i) {
        dw = 2 * w * dw;
        w = w * w + c;
      }
      if (cabs(dw) <= 1) {
        *period = p;
        return z0;
      }
    }
    z = z * z + c;
  }
  *period = 0;
  return 0;
}

double julia_exterior_de(complex double c, complex double z, int maxiters, double escape_radius) {
  complex double dz = 1;
  for (int n = 0; n < maxiters; ++n) {
    if (cabs(z) >= escape_radius) {
      return cabs(z) * log(cabs(z)) / cabs(dz);
    }
    dz = 2 * z * dz;
    z = z * z + c;
  }
  return 0;
}

double julia_interior_de(complex double c, complex double z, int maxiters, double escape_radius, double pixel_size, complex double z0, int period, bool superattracting, int *fatou) {
  complex double dz = 1;
  for (int n = 0; n < maxiters; ++n) {
    if (cabs(z) >= escape_radius) {
      *fatou = -1;
      return cabs(z) * log(cabs(z)) / cabs(dz);
    }
    if (cabs(z - z0) <= pixel_size) {
      *fatou = n % period;
      if (superattracting) {
        return cabs(z - z0) * fabs(log(cabs(z - z0))) / cabs(dz);
      } else {
        return cabs(z - z0) / cabs(dz);
      }
    }
    dz = 2 * z * dz;
    z = z * z + c;
  }
  *fatou = -2;
  return 0;
}

int main(int argc, char **argv) {
  int size = 512;
  double radius = 2;
  double escape_radius = 1 << 10;
  int maxiters = 1 << 13;
  if (! (argc > 2)) { return 1; }
  complex double c = atof(argv[1]) + I * atof(argv[2]);

  int period = 0;
  bool superattracting = false;
  complex double z0 = julia_attractor(c, maxiters, &period);
  if (period > 0) {
    double epsilon = nextafter(1, 2) - 1;
    complex double z = z0;
    complex double dz = 1;
    for (int i = 0; i < period; ++i) {
      dz = 2 * z * dz;
      z = z * z + c;
    }
    superattracting = cabs(dz) <= epsilon;
  }

  double pixel_size = 2 * radius / size;
  printf("P6\n%d %d\n255\n", size, size);
  for (int j = 0; j < size; ++j) {
    for (int i = 0; i < size; ++i) {
      double x = 2 * radius * ((i + 0.5) / size - 0.5);
      double y = 2 * radius * (0.5 - (j + 0.5) / size);
      complex double z = x + I * y;
      double de = 0;
      int fatou = -1;
      if (period > 0) {
        de = julia_interior_de(c, z, maxiters, escape_radius, pixel_size, z0, period, superattracting, &fatou);
      } else {
        de = julia_exterior_de(c, z, maxiters, escape_radius);
      }
      int r, g, b;
      hsv2rgb(fatou / (double) period, 0.25 * (0 <= fatou), tanh(de / pixel_size), &r, &g, &b);
      putchar(r);
      putchar(g);
      putchar(b);
    }
  }
  return 0;
}

EDIT

Using a complex-valued distance estimate (just remove all the $|.|$ from each formula) shows something interesting.

magnitude

magnitude of complex de

magnitude with phase

magnitude and phase of complex de

notice that the complex de points away from the Julia set, but towards the spots. however, the gradient of the magnitude of distance estimate points away from the Julia set and away from the spots as well:

magnitude of de with phase of gradient

magnitude of de with phase of gradient

computing the phase difference gives some kind of metric of distance estimate validity (defining invalid as near a spot), red is valid and cyan is invalid:

phase difference between de and gradient

thresholding this value isolates the spots:

thresholded phase difference

the red spot regions could have their distance estimate replaced by interpolating from their boundary, but here they are simply replaced with white:

thresholded phase difference white

This works well for the far-interior, but near the Julia set the visible effects of the spots extend further than their invalidity region - between the spot and the Julia set the distance estimate is reduced. The effect seems to be stronger for more attracting Julia sets, giving bad glitches even though the spots themselves are gone:

thresholded phase difference white c = -0.01

$\endgroup$
  • $\begingroup$ shadertoy.com/view/4dVSRw $\endgroup$ – Adam Aug 22 '18 at 18:21
  • $\begingroup$ hse.ru/data/2014/06/24/1310210120/moscow-2014.pdf $\endgroup$ – Adam Dec 9 '18 at 15:06
  • 1
    $\begingroup$ IMHO these sposts are periodic points and it's preimages. Algorithm works well as an fractional iteration count for the interior. $\endgroup$ – Adam Dec 11 '18 at 17:00
  • $\begingroup$ @Adam so, if it was a smooth iteration count I should be able to differentiate it across image space to get an approximate distance estimate, but I tried and it worked ok in some cases but not in all of them $\endgroup$ – Claude Dec 12 '18 at 20:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.