The paper "Images of Julia sets that you can trust" from 2013 was the most interesting one I've read in recent times. And after some valuable help from claude (thanks!) I was finally able to understand the details.

Here's the first example for c=-1+0i in the standard quadratic case, in a resolution of 16384x 16384 pixels (zipped 4bit BMP) and a 4096 pixel JPG-version. The 16k image took only 24 minutes to compute (the white/gray border was finished after some minutes, the rest was black detection), the 4k took 30 sec.

Is anybody else currently implementing that algorithm and aiming for speed?

However my Quadtree data stucture was so slow, I switched to a pixel array and an unordered list since - currently - I'm not computing images where space is the biggest issue, so I optimized for speed.

I start with leaves that are as small as a pixel, so I omitted the refinement process that was (with my Quadtree) immensely time-consuming (probably because of repeated cache misses and page faults).

The first thing I noticed, was, one gets very quickly a nice border if only checking for white by the rule: a gray pixel A can be coloured white if its bounding box for f(A) lies completely outside the -2..2 x -2..2 square or only hits white pixels (so no path following).

If that repeated process did not produce new white pixels, I switch to black detection. The authors state in their paper that black occurs quite late in the refinement process, but once the first is found, it grows pretty quickly. So the fewer gray pixels to test, the better. Since black can only occur in connected sets, it must be around the plane origin. So I started to check gray pixels in a small square around the origin. I'm not computing the whole cell graph (way too slow), I'm checking one gray pixel at a time and use a stack to store the yet-to-be-visited points, remove points, check their bounding box for f(A) until the stack is empty. If no white pixel was hit, then ALL the pixels ever in the stack can be coloured black at once. Second rule for new black is the same as for white: if one pixel only points (via its bounding box) to black pixels, itself is black. If these two rules do not produce new black pixels, I widen the initial sought-in black square, checking only the newly covered regions.

There are a few things I plan to do next:

- Let the program sit for a while (on the harddrive and in my head) and check whether it's actually correct (I did a plausibility check with the image produced, that was passed: white and black pixels do not touch.

- Introduce integer arithmetic. I am using C++ double now. The plane covers -2 .. 2 (so 4 in width, 2^2) and the image size is also a power of 2. So from my understanding of floating point representation, these numbers (and their ratio as a scale for converting plane coordinates to pixels) can accurately be represented.

- A better data representation since there still are (I assume) a lot of page faults, as the routine searches for a gray pixel somehwere in the image, then checks the bounding box, which is (probably) somehwere completely different, so a new segment has to be read from memory.

- Sharpen the bounding box. I make them a bit larger than probably necessary, especially on their right side (I add 1 to the pixel coordinate after calculating it from the plane coordinates).

- Other methods to mathematically exactly detect black pixels, especially at the start.

I've read somewhere that around the boundary of Julia sets the orbit is very sensitive to the actual starting point.

Does that mean in reverse, that for points deep in a black (interior) region, the bounding boxes f(A) and f(B) for two adjacent pixels A and B lie in close vicinity as well?