This question is about the derivation / explanation of the Distance Estimation Method (DEM) for rendering Julia and Mandelbrot fractals.
I have not succeeded in finding an explanation online so I thought I'd have a go myself - so I am seeking feedback and corrections.
Standard (Divergence Testing) Algorithm
The standard method for rendering Julia and Mandelbrot fractals is to iterate $z_{n+1} = z_n^2 +c$ for complex $z$ and $c$ and see if the orbit of $z_n$ diverges. If it does, the point being tested lies outside the set.
The method is simple and popular, but doesn't reveal fine structured detail. This is because the points being tested correspond to the centre (or corner) of a pixel, and this grid is too coarse for the fine "thin" structure of the fractal sets.
Distance Estimation (by Gradient)
The Distance Estimation Method determines how far a pixel test point is from the boundary of the Mandelbrot or Julia set. Because of this, the fine structures of the fractal sets only need to be under or near a pixel to be detected.
The method relies on the assumption that the gradient of the orbit of $z_n$ is higher closer to the fractal boundary. The further away we are from that boundary, the smaller that gradient. Note, the gradient is not the same as the magnitude of the values of $z_n$, so further away the values grow large under iteration but the rate of change in the magnitude of $z_n$ between spatially nearby points is small.
Question 1: On what basis can we assume the gradient is higher closer to the boundary? Observation clearly confirms this but a mathematical reason would be useful.
Derivation
The following is my attempt to derive the algorithm calculation steps for DEM.
The iteration mapping is
$$z_{n+1} = z_n^2 + c$$
For the Julia set, we have fixed $c$ and are testing each $z_0$. Therefore the gradient of $z_n$ is with respect to $z$.
Question 2: Is it wrt to $z$ or $z_0$?
We can differentiate the above expression using the chain rule, noting that $c$ is constant so its differential is zero.
$$ \frac{d}{dz} z_{n+1} = 2 \times z_n \times \frac{d}{dz} z_{n} + 0$$
We can write this more simply as
$$z_{n+1}' = 2 z_n z_n'$$
Whilst this may look unhelpful, this form actually works well with the iterative algorithm. That is, at each iteration:
- we calculate $z_{n+1}$ from the previous $z_n$ using $z_{n+1} = z_n^2 + c$
- we can also calculate $z_{n+1}'$ from the previous $z_n'$ using $z_{n+1}' = 2 z_n z_n'$
Question 3: What should be the value of $z_0'$ - and why?
Many set $z_0'$ to 1 without explanation. I can see it is multiplicatively neutral but there must be a better reason. A value of zero would clearly not work as it would set all $z_n'$ to zero.
For the Mandelbrot set, we have fixed $z_0 = 0$ and are testing different $c$. Therefore the gradient of $z_n$ is with respect to $c$.
The differentiation is slightly different now.
$$ \frac{d}{dc} z_{n+1} = 2 \times z_n \times \frac{d}{dc} z_{n} + 1$$
We can write this more simply as
$$z_{n+1}' = 2 z_n z_n' + 1$$
Question 4: Again, what should the initial value $z_0'$ be set to - and why? Some suggest it should be zero but without explanation. There is some discussion online that some implementations have it erroneously set to 0, although I don't think in practice it make a significant difference to rendering detail.
Question 5: What is a good rule or heuristic to determine how many iterations should be run with this method?
With the standard "divergence test" we can stop when the escape condition is reached or when a tolerable (time) maximum number of iteration is reached. Should the heuristic test the second derivative and see if it falls below a small number? All I have seen elsewhere is that the usual $|z_n|>2$ test is loosened to $|z_n|>4$ to, quote, "allow the gradient to develop" which is far from scientific in approach.
Q1: an informal argument: the gradient at escape is inversely proportional to the width of bands of equal iteration count for a fixed escape radius; as iteration count goes to $\infty$ as you approach the boundary, and there is a finite distance to the boundary, the width of the bands must go to $0$, and so the gradient at escape goes to $\infty$.
Q2: $z_0$
Q3: $z_0' = \frac{dz_0}{dz_0} = 1$ for Julia set, because $\frac{dx}{dx} = 1$ for all $x$
Q4: $z_0' = \frac{dz_0}{dc} = 0$ for Mandelbrot set, because $z_0$ does not vary with $c$. (Your paragraph before this question has $\frac{d}{dz}$ that should be $\frac{d}{dc}$)
Q5: the formula is only valid in the limit as $n\to\infty$, but after escape from a large radius then $z^2+c$ is roughly the same as $z^2$ (assuming $c$ is small), which means that the estimated distance will stay roughly the same from that point on: increasing the escape radius reduces the error but doesn't eliminate it completely.