Fractalforums

I was reading quite a number of articles about Julia sets, but kept wondering, how one can be sure, if the image computed is actually correct (besides software glitches)? (Or am I just reading the wrong texts?)

Or more bluntly spoken: If I compute a Julia set to be a square, how can I be sure that after a trillion iterations, it won't turn into a triangle? I am referring to what is seen on the screen: a finite set of discrete, rational points, that is.

I was computing some sets and increasing the iteration count sometimes showed that the set actually was dust or even empty, although the initial image looked quite nice.

In case a point escapes, one knows its fate for sure. If it gets stuck in a cycle, definite knowledge as well. But the rest? Is there a specific Julia set (or more precisely, an image thereof at a given resolution), where all points are absolutely certainly determined in their fate? (Aside from the circle at c=0, or the line at c=-2).


Linkback: https://fractalforums.org/fractal-mathematics-and-new-theories/28/julia-sets-true-shape-and-escape-time/2725/

http://webdoc.sub.gwdg.de/ebook/serien/e/IMPA_A/721.pdf paper
https://webserver2.tecgraf.puc-rio.br/~lhf/ftp/doc/oral/oktobermat.pdf slides

Images of Julia sets that you can trust
LUIZ HENRIQUE DE FIGUEIREDO and DIEGO NEHABIMPA, JORGE STOLF, JOAO BATISTA OLIVEIRA
Abstract: We present an algorithm for computing images of quadratic Julia sets thatcan be trusted in the sense that they contain numerical guarantees against sampling artifacts and rounding errors in floating-point arithmetic. We use cell mapping and color propagation in graphs to avoid function iteration and rounding errors. As a result, our algorithm avoids point sampling and can robustly classify entire rectangles in the complex plane as being on either side of the Julia set. The union of the regions that cannot be so classified is guaranteed to contain the Julia set. Our algorithm computes a refinable quadtree decomposition of the complex plane adapted to the Julia set which can be used for rendering and for approximating geometric properties such as the area of the filled Julia set and the fractal dimension of the Julia set.
Keywords: Fractals, Julia sets, adaptive refinement, cellular models, cell mapping, computer-assisted proofs

Creating images of fractals is an numerical computation. It's precission is finite.

https://mathoverflow.net/questions/200161/rounding-errors-in-images-of-julia-sets
If you want to be sure then use math:
https://mathoverflow.net/questions/248167/why-are-the-julia-sets-so-simple-quadratic-family
 
See also
https://en.wikibooks.org/wiki/Fractals/Mathematics/Numerical#Escaping_test
https://en.wikipedia.org/wiki/Shadowing_lemma

and maybe :
https://gitlab.com/adammajewski/repeated-square-root-
https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/cremer


HTH

Thanks for the links! I started with the Figueiredo article. It was quite interesting, although unfortunately they did not perform a comparison between their outcome and one of the "old" techniques. And I was very surprised to see that before 2013 (as they stated) no one had done any research upon that topic. If I understood it right, they peform interval arithmetics and do only need a fixed amount of precision because they enlarge their bounding box for f(A) onto the discrete rational screen grid.

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?

I used the Figuereido algorithm to construct an image of the Mandelbrot set. It's nice to see what can be done with this idea, however my implementation is computationally very expensive. a small image that's not even that sharp, took several hours.

Thanks to pauldelbrot's suggestion, I could omit testing the points in the period-2-bulb and a rectangle lying truely inside the main cardioid (it proved to cumbersome for me to reverse the exact equation for that structure, so I settled with a rectangle).

The process is quite straight forward though: Walk the plane, compute individual Julia sets for the quadratic case (hence any Mandelbrot image pixel corresponds to one single point in the complex plane, not a square like for the Julia sets) - only the white parts and a small region around the origin were actually computed - and check whether the origin can be colored black or white.

To speed things up a bit, the size of the Mandelbrot image was fixed, the size of the Julia sets incremented with each iteration. It startet with a tiny 64 pixels, and for every pixel remaining gray it was doubled in the next round. It's nice to see new black arising and the shape of the rectangle turn slowly into the cardioid.

The image below was computed in the white, black and gray parts, the other colors above the x-axis were from mirroring the lower pixels.

But for my final goal (a list of complex numbers with rational coordinates of fixed digit length, that are definitely within or outside of the Mandebort set), I learnt it's better
to work with distance estimation (just starting there).

From the true shape image for some Julia sets I constructed circumferencing polygons that separate the exterior from the rest. My first goal was to devise a polygon and use a point-in-polygon-test to judge pixels as either "exterior" or "not definitely known". In a second step the true interior is to be described as well - however that will be quite difficult because it is usually split in several parts.

Finding a polygon proved rather difficult, and I'm not sure my approach is mathematically sound for all shapes - if someone knows of a validated algorithm for that type of problem, help is greatly appreciated!
(I have to inflate the gray boundary into the white, since otherwise sometimes pixels are falsely classified because polygon edges cross the actual gray shape.)

However I perform a computational check on the Figuereido output image. For every gray or black pixel, the above test-in-polygon must reply with "not definitely known". I might extend that test to the 8 neighbours, so that they must not be cklassifed as "exterior" either, otherwise the polygon is rejected.

If, in the future, describing interior and exterior with several polygons should work (more or less tightly attached to the actual boundary), one could greatly reduce the space needed to store a very large image of a true shape Julia set.

As for the point-in-polygon test, I am currently using the brilliant even-odd-algorithm in the implmenetation by Glauner et. al https://github.com/pglauner/point_in_polygon.

Here are the first and preliminary results:
The left image shows the Julia set (a downscaled version), the right one the output of the check (white stands for exterior, gray for "not definitely known".)

The zipped file contains the original 4096 pixel version and a list of the coordinates of the polygon in p/4096 form, where p is a  whole numbers (the 4096 is omitted).

The complex plane ranges from -2 .. +2 in both axis.

Cool
can you put image with polygons marked ?

If I 'm not wrong I call such polygons: trap or target set
https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/def_cqp#Target_set


Here is the image of Julia set with interior trap
https://commons.wikimedia.org/wiki/File:Parabolic_Julia_set_for_internal_angle_1_over_7_with_flower_trap_and_critical_orbit.png

https://commons.wikimedia.org/wiki/File:Target_set_for_internal_ray_0.ogv

See also :
https://commons.wikimedia.org/wiki/File:Discrete_dynamic_in_case_of_complex_quadratic_polynomial.gif


Here is "my" polygon test
https://en.wikibooks.org/wiki/Fractals/Computer_graphic_techniques/2D/algorithms#Testing_if_point_is_inside_polygon

HTH

Thanks for the links!

If I 'm not wrong I call such polygons: trap or target set

But aren't traps artificially introduced geometric shapes upon the image? The (yet to be constructed) black polygons of the Julia sets I'm interested in, encompass regions of definite fate (here: bounded forever), so no iterational outcome.

Here is "my" polygon test

Quite sophisticated - a test, that only needs one line of code basically. In my implementation there are a lot of preliminary verifications (point of interest is a vertex of the polygon, ray is colinear etc.)

Here's the image with the blue-lined polygon marked. Sometimes it gets very close to the gray boundary (it just passes the computational verification), in other parts there are loops and the polygon crosses itself. I'm thinking about manually shaping the polygon a bit (some images have only a couple of pixels that did not pass the test, but otherwise the polygon is fine).

how do you prove that your polygon is correct?  distance estimation is prone to rounding errors...

Since the Figuereido algorithm works on squares in the plane and represents them as leaves (pixels), the fate of every complex point in a white or black pixel is known. So currently I am trying to construct the polygon in a way that there is at least one white pixel between any polygon line and the beginning of the gray area, so that I have "true" distance between the polygon and the unknown area. Since the black is even farther away, this assures that no non-white pixel is falsely identified as "outside".

However currently my polygon sometimes touches the gray area. I tried to avoid that be expanding the gray area before creating the polygon, but sometimes I need a dozen or so expansion rounds which defeats the refinement process, so I have to come up with a better way to create a polygon in the first place.

Current way is: Enlarge the gray area (every gray pixel gets its 8 eighbours set to gray, no matter what color they have). this is done 1-n times. Afterwards I identify "border pixels" as those white pixels that have at least one white and one gray neighbour. Those are marked. I start with the bottommost marked pixel (randomly chose one if there are more) and go to its nearest (by added pixel coordinate difference) marked neighbour, the rightmost one if there are more than one (and if still more than one, I choose one randomly). Visited pixels are unmarked. Then I compress the polygon by removing unnecessary points that lie colinear with the predecessor and successor in my list.

I'm however, not very satisfied with the "nearest neighbour" choice part. But I figured since I only want one polygon that works, I can check it afterwards - if the "one white pixel as true distance" part is correct.

I wanted to see whether one could extend the algorithm to Lyapunov-varied Julia sets. So I implemented the normal quadratic version f(z)=z^2+c(sequence) where the sequence is AB and c(sequence) alternates between two different c values cA and cB.

For a square Q in the complex plane, the bounding box is a rectangle around f(Q). Here a two step approach was taken to perform both the cA and the cB iteration at once, so the bounding box for Q was: rectangle around f( rectangle around f(Q,cA), cB).

This can easily be extended to arbitrary sequences, just performing one whole sequence turn around. However, the bounding box will get larger in width by each sequence kength increase, so I suspect that it becomes more difficult to judge a pixel as white or black.

The first image below contains a very nice looking Julia set with lots of separated black, however since it is disconnected I currently cannot use the polygon feature.

The second image probably does not have the most interesting shape, but it contains the smallest gray border I've seen so far. The bottom most tip of the structure has only one  gray pixel that separates black from white.

But I noticed that the algorithm is either very fast (white can always be judged very fast, but black only with some parameters) or it takes hours to come up with the answer "no black found". Maybe the time-consuming parameters are parabolic ones? I think I'll write them down next time I encounter some.

(attachments were missing from previous post, sorry)

I was finally able to figure out a better way of constructing the polygon. I check every 3x3 white grid now and send rays in every four directions until they hit a gray pixel, step back and mark the border.

Now the final polyline fulfills the requirements: at least one white pixel distance everywhere to the gray area, the complex plane squares pass the point-in-polygon-test and the polygon only consists of vertical and horizontal parts (the pip algorithm works directly on diagonal lines, the lines drawn on a screen however are series of horizontal and vertical segments with diagonal jumps, now they both "see" the same).

There's still the need for some manual correction, mainly on the last coordinates in the polygon when it is short of closing up. It often goes in circles, but that can be fixed easily.

The example shows a very fine-structured shape with lots of turns and changes in directions, a good test for the polygon routine (quadratic Julia set with sequence AB,
cA=0.65625+0.953125i, cB=0.578125-0.984375). The red rectangle shows a copy of what the bottommost tip looked like before the manual correction.

The polygon coordinates file is in succession order of X,Y coordinates as rationals with common (and omitted in the file) denominator 4096.

To check an arbitrary complex number a+bi I use a conservative approach: I write a and b as a rational with the denominator of the image, rounding their enumerator down and then check this number an its 8 "whole-number neighbours". If all lie outside the polygon, the initial number is known to be there too, in all other cases the result is "unknown".

I did a comparison of the outcomes of the true shape algorithm (TSA) and the standard escape time (ETA). Parameters were: Range -2..+2 in both axis, 4096x4096 pixels, bailout 10^6, max it 2500. (Images are paired (blue rectangle), upper one ETA, lower one
TSA). All are connected Julia sets (as judged by the ETA).

When there is a lot of black to be found by TSA, shapes of both routines were pretty similar and have a small gray border. But as soon as some spirals or swirls appear in the ETA image, the TSA gets stuck, takes hours and comes up with no black (therefore the red marked images could only be computed in 1/4 the size).

So I was wondering, is a spiral somehow indicative of a parabolic or otherwise slowly computable Julia set?

And interestingly sometimes although the ETA images are pretty different, what the TSA produces is rather similar (upper pair row: rightmost 2 pictures, and bottom pair row). I think I have to go to a higher resolution there. It seems that some parameters follow the same refinement path.

I compared the outcomes whether there were discrepancies: How often does the definite outcome of the TSA (black/white) differ from that of the ETA?

Well, short answer: Never at the computed resolution. I think only at the Julia set itself discrcpancies might occur due to small changes in initial parameters leading to big changes in the orbit. But with the current resolution, the gray area is so thick that the boundary itself is way within it.

I would like to know what happens if the TSA discovers the spirals of the 3rd image in the top row. But I have no idea at what resolution that will occur. I guess that's the first thing I take a look at after I finished the new data structure/algorithm that only looks for white but in a higher resolution (I hope to see a gray spiral then).