Fractals/fractalzoomer

Fractal zoomer is a fractal program by hrkalona (Kalonakis Christos )

Features

 * Over 150 fractal functions
 * Plane transformations
 * Color options
 * Image filters
 * Julia sets and Julia Map
 * Orbit tracking
 * 3D Heightmap
 * Uses boundary tracing algorithm for faster calculation
 * User formulas
 * User compilable functions
 * Domain Coloring
 * Polar coordinates projection
 * Multiplatform (Java Based)

Run
Windows
 * download exe file
 * Make sure you have JRE 1.8 installed.
 * Run as Administrator (in order to copy some lib files)

Linux
 * download jar file
 * Make sure you openjdk-jre installed.
 * If you want to use the edit/compile code functionality you need to install openjdk-jdk as well.
 * To execute use: java -jar [JarFileName.jar] on a terminal

java -jar ./FractalZoomer.jar

Create a Fractal
In order to create an Image, first we need to translate the pixel coordinates into complex numbers.

For our purposes lets assume that the image has the same size in each dimension.

Then we need to iterate the complex number, and based on the final result we must assign a color to the pixel.

For example this is the generic code for an escaping type fractal.

This is the generic code for a converging type fractal.

Kleinian
Algorithm is based on that described in the paper by Jos Leys The code can be found on Kleinian.java class

Escaping Type Fractals
Fractals that check if they exceeded some predefined boundary, usually use this escape check criterion $$|z_{n}| \ge bailout$$. In most cases the Euclidean norm(2) is used, but you can always use any other norm.

Converging Type Fractals
Fractals that check if they converged into a complex number, usually use this converging criterion $$|z_{n} - z_{n-1}| \le error$$.

Root Finding Methods
The converging criterion is used in root finding methods, in order to determine if a root was found.

Some examples of root finding methods are presented below. All the example images use $$p(z) = z^3 - 1$$ as their function.

Newton Method
description $$z_{n+1}=z_{n}-\frac{f(z_{n})}{f'(z_{n})}$$

Halley Method
$$z_{n+1}=z_{n}-\frac{2f(z_{n})f'(z_{n})}{2f'(z_{n})^2-f(z_{n})f''(z_{n})}$$

Schroder Method
$$z_{n+1}=z_{n}-\frac{f(z_{n})f'(z_{n})}{f'(z_{n})^2-f(z_{n})f''(z_{n})}$$

Householder Method
$$z_{n+1}=z_{n}-\frac{f(z_{n})[2f'(z_{n})^2+f(z_{n})f''(z_{n})]}{2f'(z_{n})^3}$$

Secant Method
$$z_{n+1}=z_{n}-f(z_{n})\frac{z_{n}-z_{n-1}}{f(z_{n})-f(z_{n-1})}$$

Steffensen Method
$$z_{n+1}=z_{n}-\frac{f(z_{n})^2}{f(z_{n}+f(z_{n}))-f(z_{n})}$$

Muller Method
$$q=\frac{z_{n}-z_{n-1}}{z_{n-1}-z_{n-2}}$$

$$A=qf(z_{n})-q(q+1)f(z_{n-1})+q^2f(z_{n-2})$$

$$B=(1+2q)f(z_{n})-(1+q)^2f(z_{n-1})+q^2f(z_{n-2})$$

$$C=(1+q)f(z_{n})$$

$$z_{n+1}=z_{n}-(z_{n}-z_{n-1})\frac{2C}{B\pm\sqrt{B^2-4AC}}$$

In order to select the sign for the denominator you should examine both $$|B+\sqrt{B^2-4AC}|$$ and $$|B-\sqrt{B^2-4AC}|$$ to determine the one with the larger norm.

Parhalley Method
$$z_{n+1}=z_{n}-\frac{2f(z_{n})}{f'(z_{n})\pm\sqrt{f'(z_{n})^2 - 2f(z_{n})f''(z_{n})}}$$

In order to select the sign for the denominator you should examine both $$|f'(z_{n})+\sqrt{f'(z_{n})^2 - 2f(z_{n})f(z_{n})}|$$ and $$|f'(z_{n})-\sqrt{f'(z_{n})^2 - 2f(z_{n})f(z_{n})}|$$ to determine the one with the larger norm.

Laguerre Method
$$z_{n+1}=z_{n}-\frac{f(z_{n})deg}{f'(z_{n})\pm\sqrt{(deg-1)^2f'(z_{n})^2 - deg(deg-1)f(z_{n})f''(z_{n})}}$$

Where deg, is the degree and can be any complex number.

In order to select the sign for the denominator you should examine both $$|f'(z_{n})+\sqrt{(deg-1)^2f'(z_{n})^2 - deg(deg-1)f(z_{n})f(z_{n})}|$$ and $$|f'(z_{n})-\sqrt{(deg-1)^2f'(z_{n})^2 - deg(deg-1)f(z_{n})f(z_{n})}|$$ to determine the one with the larger norm.

Bailout Conditions
A bailout condition is defined to be the stop iteration criterion. The second stop iteration criterion is when the number of iterations exceed the maximum number of iterations.

In the current implementation of the software it can only be changed for Escape Type Fractals, but the generalization obviously works in any iterating method.

If you want to set the bailout condition for a Converging Type Fractal, for instance Newton's Method for $$f(z) = z^3-1$$, you need to set the user formula to z - (z^3-1) / 3*z^2 and the algorithm to Escape Type.

At this point you can set a user bailout algorithm with the prefered stopping criterion.

Circle (Euclidean Norm)
$$|z_{n}| \ge bailout$$

Square (Infinity Norm)
$$max(|Re(z)|, |Im(z)|) \ge bailout$$

Rhombus (One Norm)
$$|Re(z)| + |Im(z)| \ge bailout$$

N-Norm
$$(|Re(z)|^n + |Im(z)|^n)^\frac{1}{n} \ge bailout$$

Field Lines
Where zold is the previous value of z.

Rotation
Rotation is technically similar to a plane transformation. Right after we obtain the complex number from the pixel coordinates, we can transform the complex number, using the rotation function.

By using the Euler's formula:

$$\begin{align} e^{i \theta} z &= (\cos \theta + i \sin \theta) (x + i y) \\ &= x \cos \theta + i y \cos \theta + i x \sin \theta - y \sin \theta \\ &= (x \cos \theta - y \sin \theta) + i ( x \sin \theta + y \cos \theta) \\ &= x' + i y' , \end{align}$$

We can also rotate about an arbitrary point. In this case the rotation is defined as:

Polar Projection


Fractal Zoomer's polar projection was implemented by using pfract's implementation. It is an example of Exponential Mapping, a technique that maps the complex plane into radius and angle.

If you want to see part of this large image in Fractal Zoomer, set the real coordinate to -1.786440255563638 If you then start zooming you will see that the zoom resembles https://commons.wikimedia.org/wiki/File:Mandelbrot_set_exponential_mapping_c%3D-1.7864402555636389.png

In order to get the other image, you need to set the rotation to 180 and also set the rotation center to the current center.

In order to render the large polar images (Mercator Maps), you will have to use the Image Expander, and set the number of square images that will be stitched together.

Color

 * color map files = files with map extension, the same format as fractint map files
 * array Color[] describing (color) palette

Palette
A palette is defined as a set of colors. It is usually stored in an array. where N is the number of colors in the palette.

The outcoloring or incoloring algorithms will produce a number, as their result, which will be then translated into a palette index, and therefore into a Color.

In the case that we do not care for continuous colors (Smoothing), the simple translation is to cast the result into an integer, and then use this as our index.

For our purposes, lets assume that the method Iterate is responsible to create its result based on the selected incoloring and outcoloring method.

If Iterate returns the maximum number of iterations as a result, we can choose the designated color, for instance black.

In the case that we want continuous colors (Smoothing), the Iterate method must create a result that will include a fractional part (e.g. 100.73).

The we need to get two consecutive colors, based on the integer part of the result, and then perform interpolation between those colors, using the fractional part.

A simple interpolation would be linear interpolation, but different methods can be used, like cosine interpolation e.t.c.

In Coloring Algorithms
Complex numbers that do not escape the predefined boundary or do not converge are part of the set (inside the set), hence the term in-coloring is used.

In most of the presented algorithms, the last value of the complex number is used, and in some special cases the second and third to last.

The algorithm can be further customized, if the user in-coloring algorithm is used.

Decomposition Like
The constants that were used, like 0.75 or 59π are completely arbitrary.

Squares
The constants that were used, like 0.75 or 59π are completely arbitrary.

Out Coloring Algorithms
Complex numbers that escape the predefined boundary or converge are not part of the set (outside the set), hence the term out-coloring is used.

In most of the presented algorithms, the last value of the complex number is used, and in some special cases the second and third to last.

The algorithm can be further customized, if the user out-coloring algorithm is used.

Escape Time
The escape time algorithm counts the iterations taken for the fractal to either escape a predefined boundary condition (bailout) or to converge into a value, with a small threshold.

In the case that continuous colors (Smoothing) is not enabled we only need to return the number of iterations, which will be integer. In the case that continuous colors (Smoothing) is enabled we need to produce a result that will include the number of iterations plus a fractional part.

Calculating Continuous Iteration Count for Escaping Type Fractals

$$A = n + \frac{log(bailout^2) - log(|z_{n-1}|^2)}{log(|z_{n}|^2) - log(|z_{n-1}|^2)}$$ Keep in mind that all the norms are squared just so we can avoid calculating the square root, the norm is defined as $$\sqrt{(re^2 + im^2)}$$
 * Method 1

$$B = n + 1 - \frac{log(\frac{log(|z_{n}|^2)}{log(bailout^2)})}{log(\frac{log(|z_{n}|^2)}{log(|z_{n-1}|^2)})}$$
 * Method 2

Increasing the bailout value will produce smoother images.

Calculating Continuous Iteration Count for Converging Type Fractals

$$error = 1e-9$$
 * Method 1

$$C = n + \frac{log(error) - log(|z_{n-1} - z_{n-2}|^2)}{log(|z_{n} - z_{n-1}|^2) - log(|z_{n-1} - z_{n-2}|^2)}$$

$$error = 1e-9$$
 * Method 2

$$D = n + \frac{log(\frac{log(error)}{log(|z_{n} - z_{n-1}|^2)})}{log(\frac{log(|z_{n} - z_{n-1}|^2)}{log(|z_{n-1} - z_{n-2}|^2)})}$$

The escapeTime function above should be adjusted to include the required arguments. As you can see it must always return a floating point number as a result.

For reference we will use the names A, B, C, and D in some of the following outcoloring methods.

In special cases like the Magnet function, that uses both Escaping and Converging methods, we must use the corresponding smoothing method.

If the final complex value escaped, use the escaping smoothing method. If the final complex value converged, use the converging smoothing method.

If you know the exact root or that a point can converge into you can change the smoothing methods to incorporate this root for better results.

For instance in Magnet functions the root it at $$1 + 0i$$.

$$error = 1e-9$$
 * Method 1

$$n + \frac{log(error) - log(|z_{n-1} - root|^2)}{log(|z_{n} - root|^2) - log(|z_{n-1} - root|^2)}$$

$$error = 1e-9$$
 * Method 2

$$n + \frac{log(\frac{log(error)}{log(|z_{n} - root|^2)})}{log(\frac{log(|z_{n} - root|^2)}{log(|z_{n-1} - root|^2)})}$$

Binary Decomposition
Binary decomposition is based on the final complex number, we can add an offset to the iterations.

Continuous Iteration Count can also be used on this algorithm.

Binary Decomposition 2
Based on the final complex number, we can add an offset to the iterations.

Continuous Iteration Count can also be used on this algorithm.

Biomorph
Based on the final complex number, we can add an offset to the iterations.

Continuous Iteration Count can also be used on this algorithm.

Color Decomposition
The constants that were used, like 0.75 or 59π are completely arbitrary.

This algorithm is slightly modified for converging type fractals, just so different roots can have different colors. Obviously the method that is used cannot map a complex number into a real number, but the results were usable.

The constants that were used, like 0.75 or 59π are completely arbitrary.

Escape Time + Color Decomposition
The constants that were used, like 0.75 or 59π are completely arbitrary.

This algorithm is slightly modified for converging type fractals, just so different roots can have different colors. Obviously the method that is used cannot map a complex number into a real number, but the results were usable.

Continuous Iteration Count can also be used on this algorithm. The constants that were used, like 0.75 or 59π are completely arbitrary.

Escape Time + Gaussian Integer
Where gaussian integer can be obtained with this function:

Escape Time + Grid
For smoother grid lines the following version can be used.

Entropy mandelbrot coloring
" rendering the M-set by computing the local entropy. For each point I compute the entropy of the normalized iteration count on a 2n+1×2n+1 grid around this point, n=2 in the examples below. As expected it behaves somewhat like a distance estimation, points close to the boundary are more disordered. But in addition it shows interesting looking "solar flares" which come from the normalization of the iteration count, i.e, they capture information about the escape."

perturbation method
Perturbation method

Threads
Since the coloring of each pixel is completely independent from each other, we can speed up the process by diving the image into segments and assiging different threads to handle each segment.

There are many different ways to segment the image, for example horizontally, vertically, in a grid e.t.c.

Lets assume the image is split into a grid. (1x1, 2x2, 3x3, ...)

Each thread will have a different gridI and gridJ assigned to it. gridSize will be 1 or 2 or 3...

Below the code for polar projection with the Threaded optimization is presented:

Greedy Drawing Algorithms
These types of algorithms try to reduce the calculated pixels, by skipping all the pixels of some area. This area must have a boundary with the same color.

These algorithms can significantly speedup the drawing process, but they can introduce errors

Boundary Tracing
This algorithm follows a boundary of a single color, and then uses a flood fill algorithm in order to paint all the pixel contained in that boundary.

Divide and Conquer
By iteratively subdividing the image into four parts (using a cross), this algorithm examines rectangular boundaries in order to skip an entire area.

This algorithm is also referenced as Mariani/Silver Algorithm.

Periodicity Checking
To prevent having to do huge numbers of iterations for points in the set, one can perform periodicity checking. Check whether a point reached in iterating a pixel has been reached before. If so, the pixel cannot diverge and must be in the set.

If periodicity checks succeeds, it is safe to assume that we will reach the maximum number of iterations, so we can stop iterating.

Super Sampling
Supersampling is a spatial anti-aliasing method, i.e. a method used to remove aliasing (jagged and pixelated edges, colloquially known as "jaggies") from images.

Basically we create a larger image, that includes more details, and then we can downscale it to a smaller size. The downscaling is performed by averaging the colors. In that way areas that contain alot of noise get smoothed.

The method presented below, does not store the data of the larger image, in order to save up space.

For each pixel of the original image, it also calculates some nearby points by using a smaller step.



This method also includes the Threaded optimization.

Below the code for polar projection, using supersampling, with the Threaded optimization is presented:

If you want to minimize the use of trigonometric and exponential functions when supersampling, you can use some of the identities of the functions.

For the Exponential Functions:

We will use the following identities:

$$e^{\alpha + \beta} = e^{\alpha}e^{\beta}$$

$$e^{\alpha\beta} = (e^{\alpha})^{\beta}$$

$$e^{-\alpha} = \frac{1}{e^{\alpha}}$$

Therefore:

$$e^{x*dx + a} = e^{x*d_{x}}e^a$$

$$e^{x*dx + 2a} = e^{x*d_{x}}(e^a)^2$$

$$e^{x*dx - a} = e^{x*d_{x}}\frac{1}{e^a}$$

$$e^{x*dx - 2a} = e^{x*d_{x}}\frac{1}{(e^a)^2}$$

You must precalculate $$e^a$$ once, and replace it to all expressions by using the identities.

For the Trigonometric Functions:

We will use the following identities:

$$\sin(\alpha \pm \beta) = \sin \alpha \cos \beta \pm \cos \alpha \sin \beta $$

$$\cos(\alpha \pm \beta) = \cos \alpha \cos \beta \mp \sin \alpha \sin \beta$$

$$\sin(2\alpha) = 2\sin(\alpha)\cos(\alpha)$$

$$\cos(2\alpha) = 2\cos(\alpha)^2 - 1$$

Therefore:

$$\cos(y*d_{y} + a) = \cos(y*d_{y}) \cos(a) - \sin(y*d_{y}) \sin(a)$$

$$\cos(y*d_{y} - a) = \cos(y*d_{y}) \cos(a) + \sin(y*d_{y}) \sin(a)$$

$$\cos(y*d_{y} + 2a) = \cos(y*d_{y}) \cos(2a) - \sin(y*d_{y}) \sin(2a)$$

$$\cos(y*d_{y} - 2a) = \cos(y*d_{y}) \cos(2a) + \sin(y*d_{y}) \sin(2a)$$

$$\sin(y*d_{y} + a) = \sin(y*d_{y}) \cos(a) + \cos(y*d_{y}) \sin(a)$$

$$\sin(y*d_{y} - a) = \sin(y*d_{y}) \cos(a) - \cos(y*d_{y}) \sin(a)$$

$$\sin(y*d_{y} + 2a) = \sin(y*d_{y}) \cos(2a) + \cos(y*d_{y}) \sin(2a)$$

$$\sin(y*d_{y} - 2a) = \sin(y*d_{y}) \cos(2a) - \cos(y*d_{y}) \sin(2a)$$

You must precalculate $$\cos(a)$$, $$\sin(a)$$ once, and replace them to all expressions by using the identities.

Files
File types
 * jar file = program ( executable file)
 * *.fzs file Fractal Zoomer Settings ( binary file)
 * map extension is used for the color maps
 * images

Images

 * Category: Fractals_created_with_Fractal_Zoomer in commons

Source Code
Repositories:
 * sourceforge repository
 * github repository
 * An Expression Parser/Evaluator for Complex Numbers