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
c = 0, super-attracting
c = -1, super-attracting
c = -0.5, attracting
c = -0.01, attracting
c = -0.99, attracting
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 with phase
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
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:
thresholding this value isolates the spots:
the red spot regions could have their distance estimate replaced by interpolating from their boundary, but here they are simply replaced with 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:
